function gradfd = getGradientFD(problem, x, storedb, key)
% Computes an approx. of the gradient w/ finite differences of the cost.
%
% function gradfd = getGradientFD(problem, x)
% function gradfd = getGradientFD(problem, x, storedb)
% function gradfd = getGradientFD(problem, x, storedb, key)
%
% Returns a finite difference approximation of the gradient at x for
% the cost function described in the problem structure. The finite
% difference is based on M.dim()+1 computations of the cost.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% If the cost cannot be computed, an exception is thrown.
%
% See also: approxgradientFD

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

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % This gradient approximation is based on the cost:
    % check availability.
    if ~canGetCost(problem)
        up = MException('manopt:getGradientFD:nocost', ...
            'getGradientFD requires the cost to be computable.');
        throw(up);
    end
    
    
    % Default parameters. See approxgradientFD for explicit user access to
    % these parameters.
    stepsize = 2^-23;
    subspacedim = [];
    
    
    % 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);
        % Evaluate the cost there
        fxk = getCost(problem, xk, storedb);
        % Finite difference
        df(k) = (fxk - fx)/stepsize;
    end
    
    % Build the gradient approximation.
    gradfd = lincomb(problem.M, x, B, df);
    
end
