function [P,L,R]=altproj(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres)
%function [P,L,R]=altproj(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres)
%INPUT: 
%  n,m  input matrix sizes of A,B
%  A,B  cells of size k containing input of Hermitian matrices
%        for constraints   T(Ai)=Bi  i=1...k
%        We assume the last of these constraint is  unital T(In)=Im
%        We are including the trace preserving constraint
%  toler desired stopping tolerance   default .01
% iterlimit   iteration limit    default 30
% P0          initial estimate for compl pos Choi matrix P, 
%                    use P0=[] if initial estimate not available
% rankP     desired rank for final P - means projection on faces of rank
%                                          rankP are used. (default mn)
% unital    is 1 (default) if unital constraint is included and 0 otherwise
%           If 1, then it is assumed that the last constraint in A is the
%           unital constraint A{k}=eye(n)/n;B{k}=eye(m)/m;
% tracepres is 1 (default) if trace preserving constraint 
%                    should be included and 0 otherwise.
%               If 0, then it is added in matrix representation in
%               matrixrep.m
%%%%%%
%OUTPUT: P      Hermitian size m*n for compl pos. T
%        L,R  the matrix representation and RHS after removing
%                redundant constraints but before changing column order
%                This is needed when trying to reduce rank
%%%%%%%%%%%%%%%%%%%
%%



fromstart=cputime;   % time complete algorithm
%%%%%%%%%%%%%Initialize
mn=m*n;
tmnm=mn*(mn-1)/2;
svindu=logical(triu(ones(mn),1));
sqrt2=sqrt(2);
k=length(A);
if ~exist('unital','var')
	unital=1;
end
if ~exist('tracepres','var')
	tracepres=1;
end
if ~exist('toler','var')
	toler=.01;
end
if ~exist('iterlimit','var')
        iterlimit=30;
end
if ~exist('rankP','var'),
	rankP=mn;
	fprintf('setting rankP to mn =  %i  \n',mn);
end
if ~exist('P0','var')   ||   size(P0,1)==0,
        P0=mn*blkdiag(eye(rankP),zeros(mn-rankP));  % initial starting point
	fprintf('using [mn*I_r 0] as initial starting point for P \n')
end
if rankP<mn,
fprintf('solving using alt. proj. with fixed rank <=   %i \n',...
     rankP)
else
fprintf('solving using alt. proj. with NO fixed rank  \n')
end
P=P0;
clear P0;
%find the matrix representation of the constraints and RHS
startLtime=cputime;
L=matrixrep(m,n,A,tracepres);  % find matrix repr. of comp. pos.  T(A)
if tracepres,
       R=zeros(k*m^2+n^2,1);
else
       R=zeros(k*m^2,1);
end
for ii=1:k,
	R((ii-1)*m^2+1:ii*m^2)=sHvec(B{ii});
end
if tracepres,
    temp=m^2*(1:k);
    temp=setdiff(1:size(L,1),temp);
    %% rank is too expensive
    %fprintf('rank of L before and after removing k redund rows  %i  %i \n',...
    %       rank(full(L)),rank(full(L(temp,:))));
           %rank(full(L)),rank(full(L(1:end-n,:))));
    R(k*m^2+1:end)=[zeros(n^2-n,1);ones(n,1)];
    R=R(temp);  % removes redundant last constr for k blocks of constrs
%   if rank(full(L))>rank(full(L(temp,:))),
%	    fprintf('ranks L are different ... enter   return   to continue \n')
%	    keyboard
%    end
    L=L(temp,:);  % removes redund. last constraint of each of k blocks
    fprintf('solving problem with the trace preserving constraint \n')
else
    fprintf('solving problem withOUT the trace preserving constraint \n')
end
%%%???? change rows to get better block diagonal structure????? here???
%%% remove condition number check as too expensive
%fprintf(' condition number of new matrix repres. L is %g \n',cond(full(L)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('time for generating L and RHS is  %g   \n',cputime-startLtime);

residual=R-(L*sHvec(P));
noiter=0;
norms=[];
cosines=[];
normPdiffold=norm((P));    % using P_-1 = 0
%%%%%%%%%%%%%%%%%%%test permutation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%permute rows first%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
indsrows=[];
if tracepres,
   for ii=1:m^2-1,
       indsrows=[indsrows ((1:m^2-1:(k*(m^2-1)))+(ii-1))];
   end
   indsrows=[indsrows length(indsrows)+1:size(L,1)];
else
   for ii=1:m^2,
       indsrows=[indsrows ((1:m^2:(k*(m^2)))+(ii-1))];
   end
end
L=L(indsrows,:);
R=R(indsrows,:);
%figure(1)
%spy(L)
%title('matrix repres. of  L after row permutations/before colamd')
tic
pcols=colamd(L);    % permutation of cols for sparsity
Lp=L(:,pcols);
%figure(2)
%spy(Lp)
%xlabel(['matrix repr. L after row permut. and colamd, m,n,k is  ' num2str([n m n k])])
[~,pcolsinv]=sort(pcols);   % inverse perm for cols
%%%%%%%%%%%%%%%%%testing  permute rows first%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% temporary fix!!!!! for roundoff ... get/find structure!!!
fprintf('temporarily!!!! permute cols using colamd for sparsity!!!! \n')
%pinvL=(pinv(full(L(:,pcols))));  % what if L is sparse ????? and it is!!!!
%pinvL=Lp'*inv(Lp*Lp');
pinvL=(inv(Lp*Lp')*Lp)';   % taking ' at end appears to be faster
fprintf('time for colamd and pinv(L) is  %g   \n',toc);
indstemp=find((abs(pinvL)<1e-10)& (abs(pinvL)>0));
fprintf('temporarily!!! removing %i small elements in pinvL - roundoff \n',...
           length(indstemp));
pinvL(indstemp)=0;
%figure(3)
%spy(pinvL)
%xlabel(['gen inv after removing roundoff, m,n,k is  ' num2str([n m n k])]);
%%Lt=L*L';
%%figure(4)
%%spy(Lt)
%%title('L L'' gives ')
if ~issparse(pinvL),
     fprintf('changing pinvL to sparse(pinvL)')
     pinvL=sparse(pinvL);
end
%%%% temporary fix!!!!! for roundoff ... get/find structure!!! of L!!!
fprintf('nnz of L, pinvL and density pinvL are  %i %i %g \n',...
       nnz(L),nnz(pinvL),nnz(pinvL)/numel(L));
%%%possible 'large' sparse case for pinvL???
%LL=sparse(L); [uu,ss,vv]=svds(LL,size(LL,1)); pinvsvd=vv*(ss\uu');
feasdetect=0;
fprintf('\n iter norms-ratio     cos-vectors   norm-residual eig-proj \n')
startwhiletime=cputime;   % start iterations
while noiter < iterlimit &&...
        max(norm(residual),normPdiffold)>toler && feasdetect<4,
	noiter=noiter+1;
	Pvec=pinvL*residual; % Pl-P;
	Pvec=Pvec(pcolsinv); % using inverse cols perm. for perm in pinvL
	% construct Pdiffold=sHMat(Pvec) here
	Sympart=zeros(mn);
	Sympart(svindu)=Pvec(1:tmnm)/sqrt2;
	Sympart=Sympart+Sympart';
	Skewpart=zeros(mn);
	Skewpart(svindu)=Pvec(tmnm+1:2*tmnm)/sqrt2;
	Skewpart=Skewpart-Skewpart';
	Pdiffold=Sympart+1i*Skewpart+diag(Pvec(2*tmnm+1:end)); % 1i is sqrt(-1)
	%Pdiffoldtemp2=...
	%     usHMat(Pvec(1:tmnm),mn) +...
	%           i*uskewHMat(Pvec(tmnm+1:2*tmnm),mn) +diag(Pvec(2*tmnm+1:end));
	%Pdiffoldtemp=sHMat(pinvL*residual); % Pl-P
	Pl=P+Pdiffold;  % projection onto linear manifold
	%%%check Pl is on linear manifold
	normPdiffnew=norm(Pl-P,'fro');
	norms=[norms normPdiffnew/normPdiffold];
	%fprintf('norm of residual after proj onto lin manifold  %g \n',...
	%                norm(full(R-(L*sHvec(Pl)))));
	%fprintf('ratio of norms after proj onto lin manifold  %g \n',...
	%                normPdiffnew/normPdiffold);
	[u,d]=eig(Pl);
	d=real(d);  % to avoid roundoff in eig calculation
	[dd,indsall]=sort(diag(d),'descend');   
        indsneg=find(dd<=0);
        % project using largest pos. rankP eigs
	indspos=indsall(1:min(mn-length(indsneg),rankP)); 
        P=u(:,indspos)*d(indspos,indspos)*u(:,indspos)';
	           % projection onto psd
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if size(P,1)==0,           % proj. of nsd is 0
		P=zeros(mn);
	end
	P=(P+P')/2;    
	P=(n/trace(P))*P;  % normalize for trace
	Pdiffnew=P-Pl;
	cosvalue=real(trace((Pdiffold/norm(Pdiffold,'fro'))'*...
	                  (Pdiffnew/norm(Pdiffnew,'fro'))));
	cosines=...
	   [cosines cosvalue];
	RP=real(P);
	IP=imag(P);
	vP=[sqrt2*RP(svindu);sqrt2*IP(svindu);diag(P)];
	if norm(vP-sHvec(P))>1e-12,
		'errror in sHvec P'
		keyboard
	end
        residual=R-(L*vP); % update residual
        %residual=R-(L*sHvec(P)); % update residual
	normratio=normPdiffnew/normPdiffold;
	%%%?????need to decide theoretically how to detect infeas
	%%%/??? what accuracy to use below
	if abs(1-normratio)+abs(-1-cosvalue)<1e-6, % not converging
	           feasdetect=feasdetect+1;
	end
	if 100*round(noiter/100)==noiter,
	  fprintf('%i  %10.4g  %14.4g  %14.4g  %8.4g \n',...
	  noiter,normPdiffnew/normPdiffold,cosvalue,norm(residual),norm(P-Pl));
	end
	normPdiffold=normPdiffnew;
	%%% check for separating hyperplane after ?? iterations???
	%%%%?????
end
fprintf('%i  %10.4g  %14.4g  %14.4g  %8.4g \n',...
         noiter,normratio,cosvalue,norm(residual),norm(P-Pl));
fprintf('time for iterations/while loop is  %g   \n',cputime-startwhiletime);
if feasdetect>=4,
	fprintf('algorithm shows - infeasibile???\n')
	fprintf('a separating hyperplane is ????\n')
end
%%%%%%%%%%%remove comments to plot ratios of norms and cosines
%figure(1)
%plot(norms)
%title('ratios of norms of consecutive iterates')
%figure(2)
%plot(cosines)
%title('cosines of angles between consecutive iterates')
fprintf('\n final norms ratio is  %g ; final cosine is  %g  \n',...
                  norms(end),cosines(end));
%%%%%%%%%%%%%%%%%verify/check constraints
for kk=1:k,
	if norm(B{kk}-TP(m,n,P,A{kk}))>(max(m,n)*toler),
		if kk==k && unital,
		   fprintf('error for unital k= %i for ||TP(A_k)-B_k||,  is  %g \n',...
		   kk,norm(B{kk}-TP(m,n,P,A{kk})));
		else
		   fprintf('error for k= %i for ||TP(A_k)-B_k||,  is  %g \n',...
		   kk,norm(B{kk}-TP(m,n,P,A{kk})));
	        end
	end
end
if tracepres && norm(eye(n)-partialtrace(m,n,P))>((mn)*toler),
          fprintf('error for partial trace ||parttr(P)-I_n||,  is  %g \n',...
		   norm(partialtrace(m,n,P)-eye(n)));
end
fprintf('total time for for algorithm from altproj.m file is  %g \n',...
cputime-fromstart)
