clear all
%profile off
%profile on
%% April 5/14 : version 2 of code started
%%     Solving for Choi matrix and then attempt to reduce the rank
%% April 7/14: Further debugging done. 
%%             Still need to use matlab lu command to get proper Lnew
%%             with proper row rank.
%% april23/14  using only WITH TP constraint! AND unital constraint
%% april25/14  Further polishing done. Removed redundant options since
%%               both TP and unital constraints are always included.
%%% may 11/14   some further polishing to output
%%  jun21/14    Removing global vrble statements/polishing
%%              Also removing the save/load at start
%%              Being careful to distinguish between new r value and 
%%              finding max rank P when restarting DRalg


%% comment out next two lines to repeat the same 'random' problem
rng('shuffle')
%rng('default')
saverng=rng;
rng(saverng)
fprintf('\n\nStarting a *****new problem***** \ngenerating and solving TPCP problem; using **max rank** adaptation \n')
fprintf('We assume that both the TP and unital constraint are included so order of A_i m = order of B_i n \n')
m=input('USER: enter desired order of A_i and B_i , m=n value [m=n=15]');
if isempty(m),
	 m=15;
end
n=m;
k=input('enter desired number of constraints, k value [k=10]');
if isempty(k),
	 k=10;
end
kinit=k;  % save k to avoid increasing by 1 each time for unital constr.
Pmax=input('enter 1,0 whether maxrank P is desired or not [=1]');
if isempty(Pmax),
	 Pmax=1;
end
fprintf('now enter the desired rank values r - one value or a vector \n')
rrange=input('enter desired original ranks of feasible TPCP Ptr, r value [rrange=19]');
if isempty(rrange),
	 rrange=19;
end
fprintf('input values [m n k Pmax] are  %i %i %i %i %1 \n',m,n,k,Pmax);
fprintf('r values are:   ')
for rr=rrange,
	fprintf(' %i ',rr)
end
fprintf('\n')

totaltimes=[];
ranklengths=[];
minfinaliters=[];
maxfinaliters=[];
meanfinaliters=[];
maxcos=[];
maxranks=[];
for r=rrange,
%%%%% the save, clear, load is needed due to the 'global variables conflict'
%save('datamnkr.mat','maxcos','maxranks','minfinaliters','maxfinaliters','meanfinaliters','Pmax','m','n','k','r','rrange','saverng','totaltimes','ranklengths')
%clear all
%load('datamnkr.mat','maxcos','maxranks','minfinaliters','maxfinaliters','meanfinaliters','Pmax','m','n','k','r','rrange','saverng','totaltimes','ranklengths')
rng(saverng)   % use seed/state from last end before save/clear

mn=m*n;
toler=1e-14;   % tolerance for residual in alt proj
iterlimit=3500;
fprintf('current tolerance for desired-accuracy; and iterlimit are  %g , %i \n',...
          toler,iterlimit)
%%%%test without knowledge of original rank 
P0=[];
%i=sqrt(-1);   % better to use    1i
  % generate Fi's using unitary Fi and prob distrib. p to get Ptr
	fprintf('generating cell F with unitaries for generating initial feasible Ptr \n')
	F=cell(r,1);
	p=rand(r,1)+1;
	p=p/sum(p);             % prob distrib.
	for rr=1:r,
		[temp,~]=qr(randn(m,n)+1i*randn(m,n));
		F{rr}=temp;
	end
	Ptr=FtoP(F,p,toler);
rankPs{1}=rank(Ptr);
% add I for last cell for unital constraint
k=kinit;    % initialize to original k if in for loop
k=k+1;
A=cell(k,1);   
B=cell(k,1);
for ii=1:k,
    if ii==k,
      A{k}=eye(n)/n;
      B{k}=eye(m)/m;
      fprintf('last constraint k= %i is the unital constraint \n',k);
    else
	[q,~]=qr(randn(n)+1i*randn(n));
	d=10*rand(n,1);
	%d=[1;zeros(n-1,1)];
	d=d/sum(d);
	A{ii}=q*diag(d)*q';
	A{ii}=(A{ii}+A{ii}')/2;
	%[q,~]=qr(randn(m)+1i*randn(m));
	%d=10*rand(m,1);
	%d=d/sum(d);
	%B{ii}=q*diag(d)*q';
	%B{ii}=(B{ii}+B{ii}')/2;
	B{ii}=TP(m,n,Ptr,A{ii});  % use image to try and get feasible problem
%testAB
     end
end
starttotal=cputime;
[Pcurrs,ranks,totaliters,maxcosvalue]=...
               maxrankTPCP(n,m,r,A,B,Pmax,toler,iterlimit,P0);
endtotaltime=cputime-starttotal;
fprintf('total time from maxrankTPCP steps is  %g   \n',endtotaltime)



fprintf('feasible Ps are in cell Pcurrs \n')
fprintf('%i rank values are in variable ranks  ',length(ranks))
for rr=1:length(ranks),
    fprintf(' %i ',ranks(rr))
end
fprintf('\n')
fprintf('%i iteration values are in variable totaliters ',length(totaliters))
for rr=1:length(totaliters),
    fprintf(' %i ',totaliters(rr))
end
fprintf('\n')


maxcos=[maxcos maxcosvalue]; % save max of the iterations
maxfinaliters=[maxfinaliters max(totaliters)]; % save max of the iterations
minfinaliters=[minfinaliters min(totaliters)]; % save min of the iterations
meanfinaliters=[meanfinaliters mean(totaliters)]; % save mean of the iterations
totaltimes=[totaltimes endtotaltime];  % need?
ranklengths=[ranklengths length(ranks)];
maxranks=[maxranks max(ranks)];

end   % for r=rrange
fprintf('\n table created for max rank for each row r value \n')
%input.tableData=[rrange' finaliters' ranklengths'];
%disp([num2str(input.tableData)])    % changed in latexoutputfromrun.m
latexoutputfromrun   % change data to latex table
fileID=fopen('outputtable.tex','w');
for rr=1:size(latex,1)
	fprintf(fileID,'%s \n',latex{rr});
end
fclose('all');
!pdflatex outputtable > temp
!pdflatex outputtable > temp
disp('files: outputtable.tex/outputtable.pdf for latex/pdf table output')

%profile report
