docs: 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]])