function [Pcurrs,ranks,totaliters,maxcos]=maxrankTPCP(n,m,r,A,B,Pmax,toler,iterlimit,P0)
% [Pcurrs,ranks,totaliters,maxcos]=maxrankTPCP(n,m,r,A,B,Pmax,toler,iterlimit,P0)
%    Main file for finding the TPCP matrix with both TP and unital
%        constfraints included; and the largest rank is desired.
%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
% Pmax   Indicator; 1 for finding maxrank feasible solution
%                   0 for just one iteration for finding a feasible soln
%  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
%%%%%%
%OUTPUT: Pcurrs      cell, Hermitian size m*n for compl pos. T optimal Ps
%                          (the last P has the largest rank found)
%        ranks       vector with corresponding ranks for the optimal Ps
%      totaliters    vector of iterations for all DRalg steps
%      maxcos        maximum of cosine values from DRalg ... uses last max



if ~exist('P0','var'),
	P0=[];
end
mn=m*n;  % order of P - so max rank possible
k=length(A);
rankdiff=1;
ranks(1)=0;
rankstep=1;
Ic=eye(mn);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%use global variables for NOW!!!! fix this later!!!!???
% the global variables our evaluated/described in matrixrep.m
% they are used between matrixrep and DRalg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
totaliters=[];
while rankdiff>0 && ranks(end)<mn, % keep iterating while the rank increases
startcpu=cputime;
[P,noiter,maxcos]=DRalg(n,m,A,B,toler,iterlimit,P0);
totaliters=[totaliters noiter];
fprintf('total cpusecs for this step of alt proj algor/Doug.-Ra. is  %g \n',...
                  cputime-startcpu)
Popt{rankstep}=P;    % latest optimal P
Pcurr=zeros(mn);
for ii=1:length(Popt),
	Pcurr=Pcurr+Popt{ii}/length(Popt);
end
Pcurrs{rankstep}=Pcurr;
ranks(rankstep)=rank(Pcurr);
if rankstep>1,
      rankdiff=ranks(rankstep)-ranks(rankstep-1);
else
      rankdiff=ranks(rankstep);
end
if Pmax==0,   % ensure only one iteration since rankmax not desired at input
	rankdiff=0;
end
direction=Ic-trace(Pcurr)*Pcurr;   % direction to find new starting point
alpha=2*norm(direction);
P0=Pcurr+alpha*direction;   % new starting point
[~,Pp]=chol(P0);
while Pp==0,
     alpha=2*alpha;
     P0=Pcurr+alpha*norm(direction)*direction;   % new starting point
     [~,Pp]=chol(P0);
end
rankstep=rankstep+1;
fprintf('\n - for [m n k r] =  %i %i %i %i; current ranks of Ps =:  ',m,n,k,r)
for rr=1:length(ranks),
    fprintf(' %i ',ranks(rr))
end
end   % while
fprintf('\n the final ranks for optimal Ps are  ')
for rr=1:length(ranks),
    fprintf(' %i ',ranks(rr))
end
fprintf('\n')
fprintf('From maxrankTPCP: final rank is %i and order of P  is %i \n\n',...
            ranks(end),mn);


%%%%%%%%%%%%%%%check on equations for last full size facially reduced P
P=Pcurrs{end};
for kk=1:k,
	if norm(B{kk}-TP(m,n,P,A{kk}))>(max(m,n)*toler),
		if kk==k,
		   fprintf('error for last 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 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('error check for the mean of P=Pcurrs{%i} for orig. constrs.  completed in maxrankTPCP \n',length(Pcurrs));
fprintf('%i Ps are in cell Pcurrs \n',length(Pcurrs));
if Pmax==0,
	fprintf('ONLY 1 iteration; maxrank was ***NOT*** requested \n')
end
