function M = symfixedrankYYcomplexfactory(n, k)
% Manifold of n x n complex Hermitian pos. semidefinite matrices of rank k.
%
% function M = symfixedrankYYcomplexfactory(n, k)
%
% Manifold of n-by-n complex Hermitian positive semidefinite matrices of 
% fixed rank k. This follows the quotient geometry described 
% in Sarod Yatawatta's 2013 paper:
% "Radio interferometric calibration using a Riemannian manifold", ICASSP.
%
% Paper link: http://dx.doi.org/10.1109/ICASSP.2013.6638382.
%
% A point X on the manifold M is parameterized as YY^*, where 
% Y is a complex matrix of size nxk. For any point Y on the manifold M, 
% given any kxk complex unitary matrix U, we say Y*U  is equivalent to Y, 
% i.e., YY^* does not change. Therefore, M is the set of equivalence 
% classes and is a Riemannian quotient manifold C^{nk}/SU(k). 
% The metric is the usual real-trace inner product, that is, 
% it is the usual metric for the complex plane identified with R^2.
%
% Notice that this manifold is not complete: if optimization leads Y to be
% rank-deficient, the geometry will break down. Hence, this geometry should
% only be used if it is expected that the points of interest will have rank
% exactly k. Reduce k if that is not the case.
%
% The geometry is based on the following papers (and references therein).
% Please cite the Manopt paper as well as the research papers:
%
% @INPROCEEDINGS{Yatawatta2013A,
%  author={Yatawatta, S.},
%  booktitle={Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on},
%  title={Radio interferometric calibration using a {R}iemannian manifold},
%  year={2013},
%  month={May},
%  pages={3866--3870},
%  doi={10.1109/ICASSP.2013.6638382},
%  ISSN={1520-6149},
% }
%
% @article{Yatawatta2013B,
%  author = {Yatawatta, S.}, 
%  title = {On the interpolation of calibration solutions obtained in radio interferometry},
%  volume = {428}, 
%  number = {1}, 
%  pages = {828--833}, 
%  year = {2013}, 
%  doi = {10.1093/mnras/sts069}, 
%  journal = {Monthly Notices of the Royal Astronomical Society} 
% }
%
% See also: symfixedrankYYfactory sympositivedefinitefactory


% This file is part of Manopt: www.manopt.org.
% Original author: Sarod Yatawatta, June 29, 2015.
% Contributors: Bamdev Mishra.
% Change log:
%
%   June 28, 2016 (NB):
%       Metric scaled down by factor 2 to match the metric used in
%       euclideancomplexfactory.
%
%   Apr. 18, 2018 (NB):
%       Removed lyap depedency.
%
%   Sep.  6, 2018 (NB):
%       Removed M.exp() as it was not implemented.
    
    M.name = @() sprintf('YY'' quotient manifold of Hermitian %dx%d complex matrices of rank %d.', n, n, k);
    
    M.dim = @() 2*k*n - k*k; % SY: dim of ambient space (2*k*n) - dim of kxk unitary matrix  (k^2).
    
    % Euclidean metric on the total space.
    % BM: equivalent to real(trace(eta'*zeta)), but more efficient.
    M.inner = @(Y, eta, zeta) real(eta(:)'*zeta(:));
    
    M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));
    
    % Find unitary U to minimize ||Y - Z*U||,
    % i.e., the Procrustes problem, with svd(Y'*Z).
    M.dist = @(Y, Z) distance;
    function distval = distance(Y, Z)
        [u, ignore, v] = svd(Z'*Y); %#ok<ASGLU>
        E = Y - Z*u*v'; % SY: checked.
        distval = real(E(:)'*E(:));
    end
    
    M.typicaldist = @() 10*k; % BM: To do.
    
    M.proj = @projection;
    function etaproj = projection(Y, eta)
        % Projection onto the horizontal space
        xx = Y'*Y;
        rr = Y'*eta - eta'*Y;
        Omega = lyapunov_symmetric(xx, rr);
        etaproj = eta - Y*Omega;
    end
    
    M.tangent = M.proj;
    M.tangent2ambient = @(Y, eta) eta;
    
    M.retr = @retraction;
    function Ynew = retraction(Y, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Ynew = Y + t*eta;
    end
    
    
    M.egrad2rgrad = @(Y, eta) eta;
    M.ehess2rhess = @(Y, egrad, ehess, U) M.proj(Y, ehess);
    
    
    % Notice that the hash of two equivalent points will be different...
    M.hash = @(Y) ['z' hashmd5([real(Y(:)); imag(Y(:))])];
    
    M.rand = @random;
    function Y = random()
        Y = randn(n, k) + 1i*randn(n,k);
    end
    
    M.randvec = @randomvec;
    function eta = randomvec(Y)
        eta = randn(n, k) + 1i*randn(n,k);
        eta = projection(Y, eta);
        nrm = M.norm(Y, eta);
        eta = eta / nrm;
    end
    
    M.lincomb = @matrixlincomb;
    
    M.zerovec = @(Y) zeros(n, k);
    
    M.transp = @(Y1, Y2, d) projection(Y2, d);
    
    M.vec = @(Y, u_mat) [real(u_mat(:)); imag(u_mat(:))];
    M.mat = @(Y, u_vec) reshape(u_vec(1 : n*k), [n, k]) + 1i*reshape(u_vec(n*k + 1: end), [n, k]);
    M.vecmatareisometries = @() true; 
    
end
