function M = shapefitfactory(VJt)
% Linear manifold structure for optimization over the ShapeFit search space
%
% function M = shapefitfactory(VJt)
%
% Input: VJt is a matrix of size dxn, such that VJt * ones(n, 1) = 0.
%
% Returns M, a structure describing the Euclidean space of d-by-n matrices
% equipped with the standard Frobenius distance and associated trace inner
% product, as a manifold for Manopt. Matrices on M, denoted by T, have size
% dxn and obey T*ones(n, 1) = 0 (centered columns) and <VJt, T> = 1, where
% <A, B> = Trace(A' * B).
%
% See this paper: http://arxiv.org/abs/1506.01437
% ShapeFit: Exact location recovery from corrupted pairwise directions, 2015
% Paul Hand, Choongbum Lee, Vladislav Voroninski
%
% See also: shapefit_smoothed

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 18, 2015.
% Contributors: 
% Change log: 
%
%   Jan. 25, 2017 (NB):
%       M.tangent = M.proj now, instead of being identity. This is notably
%       necessary so that checkgradient will pick up on gradients that do
%       not lie in the appropriate tangent space.
    
    [d, n] = size(VJt);

    M.name = @() sprintf('ShapeFit space of size %d x %d', d, n);
    
    M.dim = @() d*n - d - 1;
    
    M.inner = @(x, d1, d2) d1(:).'*d2(:);
    
    M.norm = @(x, d) norm(d, 'fro');
    
    M.dist = @(x, y) norm(x-y, 'fro');
    
    M.typicaldist = @() sqrt(d*n);
    
    M.proj = @(T, U) projection(U);
    VJt_normed = VJt / norm(VJt, 'fro');
    function PU = projection(U)
        % Center the columns
        PU = bsxfun(@minus, U, mean(U, 2));
        % Remove component along VJt
        % Note: these two actions can be executed separately, without
        % interference, owing to VJt having centered columns itself.
        PU = PU - (VJt_normed(:)'*U(:))*VJt_normed;
    end
    
    M.egrad2rgrad = M.proj;
    
    M.ehess2rhess = @(x, eg, eh, d) projection(eh);
    
    M.tangent = M.proj;
    
    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end
    
    M.retr = M.exp;
    
    M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];
    
    M.randvec = @(x) randvec();
    function u = randvec()
        u = projection(randn(d, n));
        u = u / norm(u, 'fro');
    end
    
    % We exploit the fact that VJt_normed belongs to the manifold
    M.rand = @() VJt_normed + randn(1) * randvec();
    
    M.lincomb = @matrixlincomb;
    
    M.zerovec = @(x) zeros(d, n);
    
    M.transp = @(x1, x2, d) d;
    
    M.pairmean = @(x1, x2) .5*(x1+x2);
    
    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [d, n]);
    M.vecmatareisometries = @() true;

end
