function [y, lambda, info] = hessianextreme(problem, x, side, y0, options, storedb, key)
% Compute an extreme eigenvector / eigenvalue of the Hessian of a problem.
%
% [u, lambda, info] = hessianextreme(problem, x)
% [u, lambda, info] = hessianextreme(problem, x, side)
% [u, lambda, info] = hessianextreme(problem, x, side, u0)
% [u, lambda, info] = hessianextreme(problem, x, side, u0, options)
% [u, lambda, info] = hessianextreme(problem, x, side, u0, options, storedb)
% [u, lambda, info] = hessianextreme(problem, x, side, u0, options, storedb, key)
% 
% (For side, u0 and options, pass [] to omit any.)
%
% Given a Manopt problem structure and a point x on the manifold problem.M,
% this function computes a tangent vector u at x of unit norm such that the
% Hessian quadratic form is minimized or maximized:
%
%    minimize or maximize <u, Hess f(x)[u]> such that <u, u> = 1,
%
% where <.,.> is the Riemannian metric on the tangent space at x. Choose
% between minimizing and maximizing by setting side = 'min' or 'max', with
% 'min' being the default. The value attained is returned as lambda, and
% is the minimal or maximal eigenvalue of the Hessian (actually, the last
% value attained when the solver stopped). This is a real number since the
% Hessian is a symmetric operator.
%
% If u0 is specified, it should be a unit-norm tangent vector at x. It is
% then used as initial guess to solve the above problem. Pass [] to omit.
%
% The options structure, if provided, will be passed along to manoptsolve.
% As such, you may choose which solver to use to solve the above
% optimization problem by setting options.solver. See manoptsolve's help.
% The other options will be passed along to the chosen solver too.
% Pass [] to omit.
%
% Often times, it is only necessary to compute a vector u such that the
% quadratic form is negative, if that is at all possible. To do so, set the
% following stopping criterion: options.tolcost = -1e-10; (for example)
% and side = 'min'. The solver will return as soon as the quadratic cost
% defined above drops below the set value (or sooner if another stopping
% criterion triggers first.)
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% info is the info struct-array returned by the solver.
%
% See also: hessianspectrum manoptsolve tangentspherefactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Aug. 13, 2014.
% Contributors: 
% Change log: 
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   May 7, 2015 (NB):
%       Default solver options: verbosity = 0 and defaults to trustregions.
%
%   Nov 27, 2015 (NB):
%       The function now also returns the info struct-array.

    
    % By default, minimize
    if ~exist('side', 'var') || isempty(side)
        side = 'min';
    end
    
    % If no initial guess was specified, prepare the empty one.
    if ~exist('y0', 'var')
        y0 = [];
    end

    % Merge default solver options with potential user-specified options.
    % Set local defaults here
    localdefaults.verbosity = 0;
    localdefaults.solver = @trustregions;
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end
    
    % Convert the side into a sign.
    % Since Manopt minimizes, 'min' asks for no sign change.
    switch lower(side)
        case 'min'
            sign = +1;
        case 'max'
            sign = -1;
        otherwise
            error('The side should be either ''min'' or ''max''.');
    end

    % We define a manifold that is actually the unit sphere on the tangent
    % space to problem.M at x. A generalization would be to consider
    % Stiefel or Grassmann on the tangent space, but this would require
    % manipulating collections of tangent vectors, which in full generality
    % may be more complex (from a programming point of view).
    % Points are represented as tangent vectors of unit norm.
    % Tangent vectors are represented as tangent vectors orthogonal to the
    % root point, with respect to the Riemannian metric on the tangent
    % space.
    
    % M is the original manifold. x is a point on M.
    M = problem.M;
    
    % N is the manifold we build. y will be a point on N, thus also a
    % tangent vector to M at x. This is a typical Riemannian submanifold of
    % a Euclidean space, hence it is easy to describe in terms of the tools
    % available for M.
    N = tangentspherefactory(M, x);
    
    % It is usually a good idea to force a gradient computation to make
    % sure precomputable things are precomputed.
    if canGetGradient(problem)
        [unused1, unused2] = getCostGrad(problem, x, storedb, key); %#ok
    end
    
    % This is the star operator of this party.
    hessian = @(y) getHessian(problem, x, y, storedb, key);
    
    % Start a Manopt problem structure for the quadratic optimization
    % problem on the sphere N.
    new_problem.M = N;
    
    % Define the cost function, its gradient and its Hessian.

    new_problem.cost = @cost;
    function [f, store] = cost(y, store)
        store = prepare(y, store);
        f = sign*store.f;
    end

    new_problem.grad = @grad;
    function [g, store] = grad(y, store)
        store = prepare(y, store);
        g = N.lincomb(y, sign*2, store.Hy, sign*(-2)*store.f, y);
    end

    new_problem.hess = @hess;
    function [h, store] = hess(y, ydot, store)
        store = prepare(y, store);
        Hydot = hessian(ydot);
        h = N.lincomb(y, sign*2, Hydot, sign*(-2)*store.f, ydot);
        h = N.proj(y, h);
    end

    % This helper makes sure we do not duplicate Hessian computations.
    function store = prepare(y, store)
        if ~isfield(store, 'ready')
            Hy = hessian(y);
            store.f = M.inner(x, y, Hy);
            store.Hy = Hy;
            store.ready = true;
        end
    end
    
    % Call a Manopt solver to solve the quadratic optimization problem on
    % the abstract sphere N.
    [y, lambda, info] = manoptsolve(new_problem, y0, options);
    lambda = sign*lambda;

end
