function gradfun = approxgradientFD(problem, options)
% Gradient approx. fnctn handle based on finite differences of the cost.
%
% function gradfun = approxgradientFD(problem)
% function gradfun = approxgradientFD(problem, options)
%
% Input:
%
% A Manopt problem structure (already containing the manifold and enough
% information to compute the cost) and an options structure (optional),
% containing one option:
%    options.stepsize (positive double; default: 2^-23).
%    options.subspacedim (positive integer; default: [], for M.dim()).
%
% If the cost cannot be computed on 'problem', a warning is issued.
%
% Output:
% 
% Returns a function handle, encapsulating a generic finite difference
% approximation of the gradient of the problem cost. The finite difference
% is based on M.dim()+1 computations of the cost.
% 
% The returned gradfun has this calling pattern:
% 
%   function gradfd = gradfun(x)
%   function gradfd = gradfun(x, storedb)
%   function gradfd = gradfun(x, storedb, key)
% 
% x is a point on the manifold problem.M, storedb is a StoreDB object,
% and key is the StoreDB key to point x.
%
% Usage:
%
% Typically, the user will set problem.M and other fields to define the
% cost (typically, problem.cost). Then, to use this generic purpose
% gradient approximation:
%
%   problem.approxgrad = approxgradientFD(problem, options);
%
% See also: steepestdescent conjugategradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Nov. 1, 2016.
% Contributors: 
% Change log: 

    % This gradient approximation is based on the cost:
    % check availability.
    if ~canGetCost(problem)
        warning('manopt:approxgradFD:nocost', ...
                'approxgradFD requires the cost to be computable.');
    end

    % Set local defaults here, and merge with user options, if any.
    localdefaults.stepsize = 2^-23;
    localdefaults.subspacedim = [];
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);
    
    % % Finite-difference parameters
    % How far do we look?
    stepsize = options.stepsize;
    % Approximate the projection of the gradient on a random subspace of
    % what dimension? If [], uses full tangent space.
    subspacedim = options.subspacedim;
                   
    % Build and return the function handle here. This extra construct via
    % funhandle makes it possible to make storedb and key optional.
    gradfun = @funhandle;
    function gradfd = funhandle(x, storedb, key)
        % Allow omission of the key, and even of storedb.
        if ~exist('key', 'var')
            if ~exist('storedb', 'var')
                storedb = StoreDB();
            end
            key = storedb.getNewKey();
        end
        gradfd = gradientFD(stepsize, subspacedim, problem, x, storedb, key);
    end
    
end


function gradfd = gradientFD(stepsize, subspacedim, problem, x, storedb, key)
% This function does the actual work.
%
% Original code: Nov. 1, 2016 (NB).
    
    % Evaluate the cost at the root point
    fx = getCost(problem, x, storedb, key);

    % Pick an orthonormal basis for the tangent space at x, or a subspace
    % thereof. The default is a full subspace. If a strict subspace is
    % picked, the returned vector approximates the orthogonal projection of
    % the gradient to that subspace.
    B = tangentorthobasis(problem.M, x, subspacedim);
    
    % Use finite differences to approximate the directional derivative
    % along each direction in the basis B.
    df = zeros(size(B));
    for k = 1 : numel(B)
        % Move in the B{k} direction
        xk = problem.M.retr(x, B{k}, stepsize);
        keyk = storedb.getNewKey();
        % Evaluate the cost there
        fxk = getCost(problem, xk, storedb, keyk);
        % Don't keep this point in cache
        storedb.remove(keyk);
        % Finite difference
        df(k) = (fxk - fx)/stepsize;
    end
    
    % Build the gradient approximation.
    gradfd = lincomb(problem.M, x, B, df);
    
end
