clear all
profile off
profile on
%% TPCP and unital (m=n)
%% k=1 (feasible when A majorizes B)
%% randomly generate A_1 and B_1
%$ A_2=I_n, B_2=I_m
%% rank of final P is at most mn
m=3;n=3;
k=1; 
mn=m*n;
r=mn;
toler=1e-12;   % tolerance for residual in alt proj
iterlimit=10000;
rankP=8;     % project onto rank r+1 ...  use mn otherwise ????
unital=1; tracepres=1;
fprintf('generating/solving unital TPCP problem with  m,n,k,r given by %i %i %i %i  \n',m,n,k,r)
i=sqrt(-1);   % better to use    1i
%% comment out next two lines to repeat the same 'random' problem
%rng('shuffle')
%saverng=rng;
%rng(saverng)

% add I for last cell for unital constraint
k=k+1;
A=cell(k,1);   
B=cell(k,1);
for ii=1:k,
    if ii==k && unital,
      A{k}=eye(n)/n;
      B{k}=eye(m)/m;
      fprintf('last constraint k= %i is the unital constraint \n',k);
    else
%     q=randn(n)+1i*randn(n);
%         %q=q(:,1); %if we want rank A=1
%         %tq=(q')*q;
%     q=q*(q');
%     A{ii}=q/trace(q);
%         %qb=randn(n)+1i*randn(n);
%         %qb=qb*(qb');
%         %B{ii}=qb/trace(qb);
%     [q,~]=qr(randn(m)+1i*randn(m));
%     d=[10*rand(m,1)];
%     d=d/sum(d);
%     B{ii}=q*diag(d)*q';
%         
    %EXAMPLE 1 where A majorizes B but the algorithm doesnt seem to work when rankP is set to certain value
        A{ii}=[0.3210 + 0.0000i, 0.1743 - 0.0657i,  -0.3102 + 0.0027i; 0.1743 + 0.0657i,   0.2081 + 0.0000i,  -0.0934 - 0.1455i;-0.3102 - 0.0027i,  -0.0934 + 0.1455i,   0.4709 + 0.0000i];
        A{ii}=A{ii}/trace(A{ii});
        B{ii}=[0.1304 + 0.0000i,  -0.0548 - 0.0460i,  -0.0564 - 0.1229i;-0.0548 + 0.0460i,   0.5758 + 0.0000i,  0.0976 - 0.1704i;  -0.0564 + 0.1229i,   0.0976 + 0.1704i,   0.2937 + 0.0000i];
        B{ii}=B{ii}/trace(B{ii});    
    %EXAMPLE 2 where A majorizes B but the algorithm does seem to work
        %A{ii}=[0.3249, -0.0838+0.1909i,-0.0329-0.0679i;
        %-0.0838-0.1909i, 0.3556, -0.2663+0.0604i;
        %-0.0329+0.0679i, -0.2663-0.0604i, 0.3194];
        %A{ii}=A{ii}/trace(A{ii});
        %B{ii}=[0.2405, 0.0364-0.1729i, -0.0642+0.0679i;
        %0.0364+0.1729i,0.5314, -0.1436+0.1128i;
        %-0.0642-0.0679i, -0.1436-0.1128i, 0.2281];
        %B{ii}=B{ii}/trace(B{ii});    
     end
end
P0=mn*blkdiag(eye(rankP),zeros(mn-rankP));

tic
[P,L,R]=altproj(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres);
fprintf('total time for L generation and alt proj algor is  %g   \n',toc)
%profile report

%CHECK IF A MAJORIZES B
v=[sort(eig(A{1}),'descend'),sort(eig(B{1}),'descend')];
v=[v,v(:,1)-v(:,2)]
m=1;
for i=1:n-1
    if (sum(v(1:i,1))>=sum(v(1:i,2)))
    m=m*1;
    else 
    m=m*0;
    end
end
if m==1
    fprintf('A majorizes B, this problem is feasible \n')
else
    fprintf('A does not majorize B, this problem is infeasible \n')
end
