function M = multinomialsymmetricfactory(n)
% Manifold of n-by-n symmetric stochastic matrices with positive entries.
%
% function M = multinomialsymmetricfactory(n) 
%
% M is a Manopt manifold structure to optimize over the set of n-by-n
% symmetric matrices with (strictly) positive entries and such that the
% entries of each column and each row sum to one.
%
% Points on the manifold and tangent vectors are represented naturally as
% symmetric matrices of size n. The Riemannian metric imposed on the
% manifold is the Fisher metric, that is, if X is a point on the manifold
% and U, V are two tangent vectors:
%
%     M.inner(X, U, V) = <U, V>_X = sum(sum(U.*V./X)).
%
% The  retraction here provided is only first order. Consequently, the
% slope test in the checkhessian tool is only valid at points X where the
% gradient is zero. Furthermore, if some entries of X are very close to
% zero, this may cause numerical difficulties that can also lead to a
% failed slope test. More generally, it is important the the solution of
% the optimization problem should have positive entries, sufficiently far
% away from zero to avoid numerical issues.
%
% Link to the paper: https://arxiv.org/abs/1802.02628.
%
% Please cite the Manopt paper as well as the research paper:
% @Techreport{Douik2018Manifold,
%   Title   = {Manifold Optimization Over the Set of Doubly Stochastic 
%              Matrices: {A} Second-Order Geometry},
%   Author  = {Douik, A. and Hassibi, B.},
%   Journal = {Arxiv preprint ArXiv:1802.02628},
%   Year    = {2018}
% }
%
% See also: multinomialdoublystochasticfactory multinomialfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Ahmed Douik, March 06, 2018.
% Contributors:
% Change log:
%
%    Sep. 6, 2018 (NB):
%        Removed M.exp() as it was not implemented.

    % Helpers
    e = ones(n, 1);
    symm = @(X) .5*(X+X');
    maxDSiters = 100 + 2*n;

    M.name = @() sprintf('%dx%d symmetric doubly-stochastic matrices with positive entries', n, n);

    M.dim = @() n*(n-1)/2;

    % Fisher metric
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        ip = sum((eta(:).*zeta(:))./X(:));
    end

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(X, Y) error('multinomialsymmetricfactory.dist not implemented yet.');

    % The manifold is not compact as a result of the choice of the metric,
    % thus any choice here is arbitrary. This is notably used to pick
    % default values of initial and maximal trust-region radius in the
    % trustregions solver.
    M.typicaldist = @() n;
   
    % Pick a random point on the manifold
    M.rand = @random;
    function X = random()
        Z = symm(abs(randn(n, n)));     % Random point in the ambient space
        X = symm(doubly_stochastic(Z, maxDSiters)); % Projection on the manifold
    end

    % Pick a random vector in the tangent space at X, of norm 1
    M.randvec = @randomvec;
    function eta = randomvec(X)
        % A random vector in the ambient space
        Z = symm(randn(n, n)) ; 
        % Projection to the tangent space
        alpha = sum((eye(n) + X)\(Z),2) ;
        eta = Z - (alpha*e' + e*alpha').*X ;
        % Normalizing the vector
        nrm = M.norm(X, eta);
        eta = eta / nrm;
    end

    % Orthogonal projection of the vector eta in the ambient space to the
    % tangent space.
    M.proj = @projection; 
    function etaproj = projection(X, eta)
        alpha = sum((eye(n) + X)\(eta), 2);
        etaproj = eta - (alpha*e' + e*alpha').*X;
    end

    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;

    % Conversion of Euclidean to Riemannian gradient
    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)
        mu = sum((X.*egrad), 2);
        alpha = (eye(n) + X)\mu;
        rgrad = (egrad - alpha*e' - e*alpha').*X; 
    end

    % First-order retraction
    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y = X.*exp(t*(eta./X));
        Y = symm(doubly_stochastic(Y, maxDSiters));
        Y = max(Y, eps);
    end

    % Conversion of Euclidean to Riemannian Hessian
    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, eta)

        % Computing the directional derivative of the Riemannian
        % gradient
        gamma = egrad.*X ;
        gammadot = ehess.*X + egrad.*eta;
        alpha = sum((eye(n) + X)\(gamma), 2);
        m = (eye(n)+X)\eta;
        alphadot = sum((eye(n) + X)\(gammadot - m*gamma), 2);
        S = (alpha*e' + e*alpha');
        deltadot = gammadot - (alphadot*e' + e*alphadot').*X - S.*eta;

        % Projecting gamma
        delta = gamma - S.*X; 

        % Computing and projecting nabla
        nabla = deltadot - 0.5*(delta.*eta)./X;
        w = sum((eye(n) + X)\(nabla), 2);
        rhess = nabla - (w*e' + e*w').*X; 
        
    end

    
    % Miscellaneous manifold functions    
    M.hash = @(X) ['z' hashmd5(X(:))];
    M.lincomb = @matrixlincomb;
    M.zerovec = @(X) zeros(n, n);
    M.transp = @(X1, X2, d) projection(X2, d);
    M.vec = @(X, U) U(:);
    M.mat = @(X, u) reshape(u, n, n);
    M.vecmatareisometries = @() false;
        
end
