r'''Bound residual norms of next deflated system.

    :param arnoldifyer: an instance of
      :py:class:`~krypy.deflation.Arnoldifyer`.
    :param Wt: coefficients :math:`\tilde{W}\in\mathbb{C}^{n+d,k}` of the
      considered deflation vectors :math:`W` for the basis :math:`[V,U]`
      where ``V=last_solver.V`` and ``U=last_P.U``, i.e.,
      :math:`W=[V,U]\tilde{W}` and
      :math:`\mathcal{W}=\operatorname{colspan}(W)`. Must fulfill
      :math:`\tilde{W}^*\tilde{W}=I_k`.
    :param g_norm: norm :math:`\|g\|` of difference :math:`g=c-b` of
      right hand sides. Has to fulfill :math:`\|g\|<\|b\|`.
    :param G_norm: norm :math:`\|G\|` of difference
      :math:`G=B-A` of operators.
    :param GW_norm: Norm :math:`\|G|_{\mathcal{W}}\|` of difference
      :math:`G=B-A` of operators restricted to :math:`\mathcal{W}`.
    :param WGW_norm: Norm :math:`\|\langle W,GW\rangle\|_2`.
    :param pseudo_type: One of

      * ``'auto'``: determines if :math:`\hat{H}` is non-normal, normal or
        Hermitian and uses the corresponding mode (see other options below).
      * ``'nonnormal'``: the pseudospectrum of the Hessenberg matrix
        :math:`\hat{H}` is used (involves one computation of a pseudospectrum)
      * ``'normal'``: the pseudospectrum of :math:`\hat{H}` is computed
        efficiently by the union of circles around the eigenvalues.
      * ``'hermitian'``: the pseudospectrum of :math:`\hat{H}` is computed
        efficiently by the union of intervals around the eigenvalues.
      * ``'contain'``: the pseudospectrum of the extended Hessenberg matrix
        :math:`\begin{bmatrix}\hat{H}\\S_i\end{bmatrix}` is used
        (pseudospectrum has to be re computed for each iteration).
      * ``'omit'``: do not compute the pseudospectrum at all and just use the
        residual bounds from the approximate Krylov subspace.
    :param pseudo_kwargs: (optional) arguments that are passed to the method
      that computes the pseudospectrum.
    :param terminate_factor: (optional) terminate the computation if the ratio
      of two subsequent residual norms is larger than the provided factor.
      Defaults to 1.

def bound_pseudo(arnoldifyer, Wt,
                 g_norm=0.,
                 G_norm=0.,
                 GW_norm=0.,
                 WGW_norm=0.,
                 tol=1e-6,
                 pseudo_type='auto',
                 pseudo_kwargs=None,
                 delta_n=20,
                 terminate_factor=1.
                 ):
    r'''Bound residual norms of next deflated system.

    :param arnoldifyer: an instance of
      :py:class:`~krypy.deflation.Arnoldifyer`.
    :param Wt: coefficients :math:`\tilde{W}\in\mathbb{C}^{n+d,k}` of the
      considered deflation vectors :math:`W` for the basis :math:`[V,U]`
      where ``V=last_solver.V`` and ``U=last_P.U``, i.e.,
      :math:`W=[V,U]\tilde{W}` and
      :math:`\mathcal{W}=\operatorname{colspan}(W)`. Must fulfill
      :math:`\tilde{W}^*\tilde{W}=I_k`.
    :param g_norm: norm :math:`\|g\|` of difference :math:`g=c-b` of
      right hand sides. Has to fulfill :math:`\|g\|<\|b\|`.
    :param G_norm: norm :math:`\|G\|` of difference
      :math:`G=B-A` of operators.
    :param GW_norm: Norm :math:`\|G|_{\mathcal{W}}\|` of difference
      :math:`G=B-A` of operators restricted to :math:`\mathcal{W}`.
    :param WGW_norm: Norm :math:`\|\langle W,GW\rangle\|_2`.
    :param pseudo_type: One of

      * ``'auto'``: determines if :math:`\hat{H}` is non-normal, normal or
        Hermitian and uses the corresponding mode (see other options below).
      * ``'nonnormal'``: the pseudospectrum of the Hessenberg matrix
        :math:`\hat{H}` is used (involves one computation of a pseudospectrum)
      * ``'normal'``: the pseudospectrum of :math:`\hat{H}` is computed
        efficiently by the union of circles around the eigenvalues.
      * ``'hermitian'``: the pseudospectrum of :math:`\hat{H}` is computed
        efficiently by the union of intervals around the eigenvalues.
      * ``'contain'``: the pseudospectrum of the extended Hessenberg matrix
        :math:`\begin{bmatrix}\hat{H}\\S_i\end{bmatrix}` is used
        (pseudospectrum has to be re computed for each iteration).
      * ``'omit'``: do not compute the pseudospectrum at all and just use the
        residual bounds from the approximate Krylov subspace.
    :param pseudo_kwargs: (optional) arguments that are passed to the method
      that computes the pseudospectrum.
    :param terminate_factor: (optional) terminate the computation if the ratio
      of two subsequent residual norms is larger than the provided factor.
      Defaults to 1.
    '''
    if pseudo_kwargs is None:
        pseudo_kwargs = {}

    # Arnoldify!
    Hh, Rh, q_norm, vdiff_norm, PWAW_norm = arnoldifyer.get(Wt)

    # get original linear system
    ls_orig = arnoldifyer._deflated_solver.linear_system

    k = Wt.shape[1]
    if k > 0:
        # smallest singular value of W^*AW
        WAW = Wt.T.conj().dot(arnoldifyer.J.dot(arnoldifyer.L.dot(Wt)))
        sigma_min = numpy.min(scipy.linalg.svdvals(WAW))

        if sigma_min <= WGW_norm:
            raise utils.AssumptionError(
                'sigma_min(W^*AW) > ||W^*GW|| not satisfied.')
        eta = GW_norm/(sigma_min - WGW_norm)
    else:
        eta = 0.
    b_norm = ls_orig.MMlb_norm
    beta = PWAW_norm*(eta*(b_norm + g_norm) + g_norm) + vdiff_norm

    # check assumption on g_norm and b_norm
    if g_norm >= b_norm:
        raise utils.AssumptionError(
            '||g_norm|| < ||b_norm|| not satisfied')

    # compute residual norms of Hh*z=e_1*b_norm
    ls_small = linsys.LinearSystem(Hh,
                                   numpy.eye(Hh.shape[0], 1) * q_norm,
                                   normal=ls_orig.normal,
                                   self_adjoint=ls_orig.self_adjoint,
                                   positive_definite=ls_orig.positive_definite
                                   )

    Solver = type(arnoldifyer._deflated_solver)
    if issubclass(Solver, linsys.Minres) or issubclass(Solver, linsys.Gmres):
        aresnorms = utils.get_residual_norms(Hh,
                                             self_adjoint=ls_orig.self_adjoint)
    else:
        # TODO: compute residuals more efficiently for CG
        try:
            solver = Solver(ls_small, tol=tol, maxiter=Hh.shape[0])
        except utils.ConvergenceError as e:
            # use all residuals that have been computed
            # (useful for short recurrences)
            solver = e.solver
        aresnorms = numpy.array(solver.resnorms)
    # absolute residual norm
    aresnorms = aresnorms * q_norm

    if pseudo_type == 'omit':
        return aresnorms / (b_norm - g_norm)

    # spectrum of Hh
    evals, evecs = scipy.linalg.eig(Hh)
    if ls_small.self_adjoint:
        evals = numpy.real(evals)

    # norm of Hh
    Hh_norm = numpy.linalg.norm(Hh, 2)

    def _auto():
        '''determine pseudo automatically'''
        # is Hh Hermitian?
        if numpy.linalg.norm(Hh-Hh.T.conj(), 2) < 1e-14*Hh_norm:
            return 'hermitian'

        # is Hh normal?
        if numpy.linalg.cond(evecs, 2) < 1+1e-14:
            return 'normal'

        return 'nonnormal'

    if pseudo_type == 'auto':
        pseudo_type = _auto()

    # for delta >= min(|\lambda|), the pseudospectrum will contain zero and
    # the thus polymax > 1. nevertheless, the bound may provide useful
    # information in early iterations with large values of delta.
    # Therefore, the maximal perturbation is chosen as the maximal
    # eigenvalue of Hh
    delta_max = 1e2*numpy.max(numpy.abs(evals))

    # minimal delta is defined via Rh
    # HACK until numpy.linal.svd (and thus numpy.linalg.norm) is fixed
    from scipy.linalg import svd
    _, Rhsvd, _ = svd(Rh[:, :1])
    delta_min = PWAW_norm*(eta*(Hh_norm + G_norm) + G_norm) + numpy.max(Rhsvd)
    if delta_min == 0:
        delta_min = 1e-16

    import pseudopy
    if not ls_small.normal:
        # construct pseudospectrum for the expected range
        pseudo = pseudopy.NonnormalAuto(Hh, delta_min*0.99, delta_max*1.01,
                                        **pseudo_kwargs)
    elif not ls_small.self_adjoint:
        pseudo = pseudopy.NormalEvals(evals)
    else:
        pseudo = None

    bounds = [aresnorms[0]]
    for i in range(1, len(aresnorms)):
        # compute roots of polynomial
        if issubclass(Solver, linsys.Cg):
            roots = scipy.linalg.eigvalsh(Hh[:i, :i])
        else:
            # TODO: more stable way of computing the roots of the MINRES
            #       poly with exploitation of symmetry?
            HhQ, HhR = scipy.linalg.qr(Hh[:i+1, :i], mode='economic')
            roots_inv = scipy.linalg.eigvals(HhQ[:i, :].T.conj(), HhR)
            roots = 1./roots_inv[numpy.abs(roots_inv) > 1e-14]

        if ls_small.self_adjoint:
            roots = numpy.real(roots)

        # compute polynomial
        p = utils.NormalizedRootsPolynomial(roots)
        if ls_small.self_adjoint:
            p_minmax_candidates = p.minmax_candidates()

        # absolute residual
        aresnorm = aresnorms[i]

        # perturbation
        # HACK until numpy.linal.svd (and thus numpy.linalg.norm) is fixed
        from scipy.linalg import svd
        _, Rhsvd, _ = svd(Rh[:, :i])
        Rhnrm = numpy.max(Rhsvd)
        epsilon = PWAW_norm*(eta*(Hh_norm + G_norm) + G_norm) \
            + Rhnrm
            #+ numpy.linalg.norm(Rh[:, :i], 2)
        if epsilon == 0:
            epsilon = 1e-16

        if pseudo_type == 'contain':
            raise NotImplementedError('contain not yet implemented')

        # exit if epsilon >= delta_max
        if epsilon >= delta_max:
            break

        delta_log_range = numpy.linspace(numpy.log10(1.01*epsilon),
                                         numpy.log10(delta_max),
                                         delta_n+2
                                         )[0:-1]

        def compute_pseudo(delta_log):
            delta = 10**delta_log
            if ls_small.self_adjoint:
                # pseudospectrum are intervals
                pseudo_intervals = utils.Intervals(
                    [utils.Interval(ev-delta, ev+delta) for ev in evals])

                # add roots of first derivative of p
                candidates = []
                for candidate in p_minmax_candidates:
                    if pseudo_intervals.contains(candidate):
                        candidates.append(candidate)
                all_candidates = numpy.r_[pseudo_intervals.get_endpoints(),
                                          numpy.array(candidates)]

                # evaluate polynomial
                polymax = numpy.max(numpy.abs(p(all_candidates)))
                pseudolen = 2 * delta

            else:
                # get pseudospectrum paths
                pseudo_path = pseudo.contour_paths(delta)

                # length of boundary
                pseudolen = pseudo_path.length()

                # evaluate polynomial on points of path
                if pseudolen > 0:
                    polymax = numpy.max(numpy.abs(p(pseudo_path.vertices())))
                else:
                    polymax = numpy.Inf

            # compute THE bound
            return pseudolen/(2*numpy.pi*delta) \
                * (epsilon/(delta-epsilon)*(q_norm + beta) + beta) \
                * polymax

        # minimization
        from scipy.optimize import minimize_scalar
        opt_res = minimize_scalar(compute_pseudo,
                                  bounds=(delta_log_range[0],
                                          delta_log_range[-1]),
                                  method='bounded',
                                  options={'maxiter': delta_n}
                                  )
        # the delta with minimal value is min_delta = 10**opt_res.x
        min_val = opt_res.fun

        # minimal bound value
        boundval = aresnorm + min_val

        # if not increasing: append to bounds
        if i > 1 and boundval/bounds[-1] > terminate_factor:
            break
        else:
            bounds.append(numpy.min([boundval, bounds[-1]]))
    return numpy.array(bounds) / (b_norm - g_norm)