function costs = surfprofile(problem, x, d1, d2, t1, t2)
% Plot the cost function as a surface over a 2-dimensional subspace.
%
% function surfprofile(problem, x, d1, d2, t1, t2)
% function costs = surfprofile(problem, x, d1, d2, t1, t2)
%
% Evaluates the cost function at points
%
%   gamma(t1, t2) = exponential_x(t1*d1 + t2*d2)
% 
% where the exponential map at x is specified by problem.M.exp (retr is
% used instead if needed). d1 and d2 are two tangent vectors to problem.M
% at the point x. The values assigned to t1 and t2 are as specified in the
% two input vectors t1 and t2.
% 
% If the function is called with an output, the plot is not drawn and the
% values of the cost are returned in a matrix of size
% length(t1)*length(t2). To plot a surf, call surf(t1, t2, costs.') (notice
% the transpose).
%
% If x is omitted, a point is generated at random. If d1 is omitted, a
% random tangent vector at x is generated. If d2 is omitted, a random
% tangent vector at x is generated, orthogonally to d1. If t1, t2 are
% omitted, they are generated with linspace's in [-1, 1].

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sep. 1, 2014.
% Contributors: 
% Change log: 
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   Nov. 12, 2016 (NB):
%       Most inputs are now optional.

    % Verify that the problem description is sufficient.
    if ~canGetCost(problem)
        error('It seems no cost was provided.');  
    end
    

    if ~exist('x', 'var') || isempty(x)
        x = problem.M.rand();
        if (exist('d1', 'var') && ~isempty(d1)) || ...
           (exist('d2', 'var') && ~isempty(d2))
            error('If x is omitted, d1, d2 should not be specified.');
        end
    end
    if ~exist('d1', 'var') || isempty(d1)
        d1 = problem.M.randvec(x);
    end
    if ~exist('d2', 'var') || isempty(d2)
        d2 = problem.M.randvec(x);
        % Make it orthogonal to d1
        coeff = problem.M.inner(x, d1, d2) / problem.M.inner(x, d1, d1);
        d2 = problem.M.lincomb(x, 1, d2, -coeff, d1);
    end
    if ~exist('t1', 'var') || isempty(t1)
        t1 = linspace(-1, 1, 51);
    end
    if ~exist('t2', 'var') || isempty(t2)
        t2 = linspace(-1, 1, 51);
    end
    
    
    if isfield(problem.M, 'exp')
        expo = problem.M.exp;
        str = 'Exp';
    else
        expo = problem.M.retr;
        str = 'Retr';
    end
    
    storedb = StoreDB();
    linesearch_fun = @(ta, tb) getCost(problem, ...
                         expo(x, problem.M.lincomb(x, ta, d1, tb, d2)), ...
                         storedb);
    
    costs = zeros(length(t1), length(t2));
    for i = 1 : length(t1)
        for j = 1 : length(t2)
            costs(i, j) = linesearch_fun(t1(i), t2(j));
        end
    end
    
    if nargout == 0
        surf(t1, t2, costs.');
        xlabel('t1');
        ylabel('t2');
        zlabel(['f(' str '_x(t1*d1+t2*d2))']);
    end
    
end
