function M = euclideansparsefactory(A)
% Returns a manifold struct to optimize over real matrices with given sparsity pattern.
%
% function M = euclideansparsefactory(A)
%
% Returns M, a structure describing the Euclidean space of real matrices
% with a fixed sparsity pattern. This linear manifold is equipped with
% the standard Frobenius distance and associated trace inner product,
% and is usable as a Riemannian manifold for Manopt.
%
% The matrices are represented as sparse matrices. Their sparsity pattern 
% is fixed. The tangent vectors are represented in the same way as points 
% since this is a Euclidean space. Point and vectors in the embedding space,
% that is, in the space of (possibly full) matrices of the same size as A,
% are represented as matrices of the same size as A, full or sparse, real.
%
% The current code relies on Matlab's built-in representation of sparse
% matrices, which has the inconvenient effect that we cannot control the
% sparsity structure: if entries of points or tangent vectors which are
% allowed to be nonzero (by the sparsity structure) happen to be zero,
% then Matlab internally restructures the sparse matrix, which may be
% costly, and which may increase computation time when using that matrix
% in combination with other sparse matrices. There is also no built-in way
% to let Matlab know that two matrices have the same sparsity structure.
% For this reason, in a future update, it will be good to try to represent
% points and tangent vectors internally as vectors of nonzeros, with truly
% fixed sparsity pattern. In the meantime, this factory is provided for
% convenience and prototyping, bearing in mind it is likely not efficient.
%
% See also: euclideanfactory euclideancomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Mar. 28, 2019.
% Change log: 
%    May 3, 2019 (NB): adapted many functions to take advantage of sparsity a bit more.
    
    dimensions_vec = size(A);
    assert(length(dimensions_vec) == 2, 'A should be a matrix (or a vector).');
    [I, J] = find(A);
    nvals = length(I);
    S = sparse(I, J, ones(nvals, 1), dimensions_vec(1), dimensions_vec(2), nvals);
      
    M.size = @() dimensions_vec;
    
    M.name = @() sprintf('Euclidean space R^(%dx%d) with fixed sparsity pattern containg %d non-zero entries', ...
                                        dimensions_vec(1), dimensions_vec(2), nvals);
    
    M.dim = @() nvals;
    
    M.inner = @(x, d1, d2) d1(:).'*d2(:); % nonzeros(d1).'*nonzeros(d2); might not work since d1, d2 might have extra zeros
    
    M.norm = @(x, d) norm(d, 'fro');
    
    M.dist = @(x, y) norm(x-y, 'fro');
    
    M.typicaldist = @() sqrt(prod(dimensions_vec));
    
    M.proj = @(x, d) S.*d; % could replace with: d(ind) where ind = find(S); which is faster?
    
    M.egrad2rgrad = @(x, g) S.*g;
    
    M.ehess2rhess = @(x, eg, eh, d) S.*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(nonzeros(x))];
    
    M.rand = @() sprandn(S);
    
    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = sprandn(S);
        u = u / norm(u, 'fro');
    end
    
    M.lincomb = @matrixlincomb;
    
    M.zerovec = @(x) spalloc(dimensions_vec(1), dimensions_vec(2), nvals);
    
    M.transp = @(x1, x2, d) d;
    M.isotransp = M.transp; % the transport is isometric
    
    M.pairmean = @(x1, x2) .5*(x1+x2);
    
    M.vec = @(x, u_mat) nonzeros(u_mat);
    M.mat = @(x, u_vec) sparse(I, J, u_vec, m, n, nvals);
    M.vecmatareisometries = @() true;

end
