function checkretraction(M, x, v)
% Check the order of agreement of a retraction with an exponential.
% 
% function checkretraction(M)
% function checkretraction(M, x)
% function checkretraction(M, x, v)
%
% checkretraction performs a numerical test to check the order of agreement
% between the retraction and the exponential map in a given Manopt
% manifold structure M. The test is performed at the point x if it is
% provided (otherwise, the point is picked at random) and along the tangent
% vector v at x if one is provided (otherwise, a tangent vector at x is
% picked at random.)
%
% See also: checkmanifold checkdiff checkgradient checkhessian

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

    if ~isfield(M, 'exp')
        error(['This manifold has no exponential (M.exp): ' ...
               'no reference to compare the retraction.']);
    end
    if ~isfield(M, 'dist')
        error(['This manifold has no distance (M.dist): ' ...
               'this is required to run this check.']);
    end

    if ~exist('x', 'var') || isempty(x)
        x = M.rand();
        v = M.randvec(x);
    end
    
    if ~exist('v', 'var') || isempty(v)
        v = M.randvec(x);
    end
    
    % Compare the retraction and the exponential over steps of varying
    % length, on a wide log-scale.
    tt = logspace(-12, 0, 251);
    ee = zeros(size(tt));
    for k = 1 : numel(tt)
        t = tt(k);
        ee(k) = M.dist(M.exp(x, v, t), M.retr(x, v, t));
    end
    
    % Plot the difference between the exponential and the retration over
    % that span of steps, in log-log scale.
    loglog(tt, ee);
    
    % We hope to see a slope of 3, to confirm a second-order retraction. If
    % the slope is only 2, we have a first-order retration. If the slope is
    % less than 2, this is not a retraction.
    % Slope 3
    line('xdata', [1e-12 1e0], 'ydata', [1e-30 1e6], ...
         'color', 'k', 'LineStyle', '--', ...
         'YLimInclude', 'off', 'XLimInclude', 'off');
    % Slope 2
    line('xdata', [1e-14 1e0], 'ydata', [1e-20 1e8], ...
         'color', 'k', 'LineStyle', ':', ...
         'YLimInclude', 'off', 'XLimInclude', 'off');
     

    % Figure out the slope of the error in log-log, by identifying a piece
    % of the error curve which is mostly linear.
    window_len = 10;
    [range, poly] = identify_linear_piece(log10(tt), log10(ee), window_len);
    hold all;
    loglog(tt(range), 10.^polyval(poly, log10(tt(range))), 'LineWidth', 3);
    hold off;
    
    xlabel('Step size multiplier t');
    ylabel('Distance between Exp(x, v, t) and Retr(x, v, t)');
    title(sprintf('Retraction check.\nA slope of 2 is required, 3 is desired.'));
    
    fprintf('Check agreement between M.exp and M.retr. Please check the\n');
    fprintf('factory file of M to ensure M.exp is a proper exponential.\n');
    fprintf('The slope must be at least 2 to have a proper retraction.\n');
    fprintf('For the retraction to be second order, the slope should be 3.\n');
    fprintf('It appears the slope is: %g.\n', poly(1));

end
