Estimate evals based on eval inclusion theorem + heuristic.

        :returns: Intervals object with inclusion intervals for eigenvalues

def _estimate_eval_intervals(ritz, indices, indices_remaining,
                                 eps_min=0,
                                 eps_max=0,
                                 eps_res=None):
        '''Estimate evals based on eval inclusion theorem + heuristic.

        :returns: Intervals object with inclusion intervals for eigenvalues
        '''
        if len(indices) == 0:
            return utils.Intervals(
                [utils.Interval(mu-resnorm, mu+resnorm)
                 for mu, resnorm in zip(ritz.values, ritz.resnorms)])
        if len(ritz.values) == len(indices):
            raise utils.AssumptionError(
                'selection of all Ritz pairs does not allow estimation.')
        if eps_res is None:
            eps_res = numpy.max(numpy.abs([eps_min, eps_max]))

        # compute quantities for bound
        delta_sel = numpy.linalg.norm(ritz.resnorms[indices], 2)
        delta_non_sel = numpy.linalg.norm(ritz.resnorms[indices_remaining], 2)
        delta = utils.gap(ritz.values[indices],
                          ritz.values[indices_remaining])
        mu_ints = utils.Intervals(
            [utils.Interval(mu+eps_min, mu+eps_max)
             for mu in ritz.values[indices]])
        mu_min = mu_ints.min_abs()

        # check gap assumption
        #if delta_sel + delta_non_sel + eps_max - eps_min >= delta:
        if delta_sel + eps_max - eps_min >= delta:
            raise utils.AssumptionError(
                'delta_sel + delta_non_sel + eps_max - eps_min >= delta'
                + '({0} >= {1}'.format(
                    delta_sel + delta_non_sel + eps_max - eps_min,
                    delta)
                )

        # check assumption on mu_min
        if mu_min == 0:
            raise utils.AssumptionError('mu_min == 0 not allowed')

        # compute eta
        #eta = (delta_sel+eps_res)**2 * (
        #    1/(delta-delta_non_sel-eps_max+eps_min)
        #    + 1/mu_min
        #    )
        #left = - delta_non_sel + eps_min - eta
        #right = delta_non_sel + eps_max + eta
        eta = (delta_sel+eps_res)**2 * (
            1/(delta-eps_max+eps_min)
            + 1/mu_min
            )
        left = eps_min - eta
        right = eps_max + eta

        # return bound
        return utils.Intervals(
            [utils.Interval(mu+left, mu+right)
             for mu in ritz.values[indices_remaining]])