%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.
%% Jul8/14 adding option for a constraint on the rank
%%         the rank constraint is implemented by changing the PSD proj.
%% Dec28/14 rerunning tests for the revision


m=input('enter desired m=n value [m=n=8]');
if isempty(m),
	 m=8;
end
n=m;
k=input('enter desired k value [k=5]');
if isempty(k),
	 k=5;
end
r=input('enter desired r value [r=9]');
if isempty(r),
	 r=9;
end
fprintf('To resolve the current problem with A_i,B_i')
fprintf('enter (rankP *pos. value*) or (*0*) if desired bound on \n')
rankP=input('rank of final P is known/or not [rankP=0]');
if isempty(rankP),
	 rankP=m*n;   % used in altproj for PSD projection
end
fprintf('input values [m n k r] are  %i %i %i %i \n',m,n,k,r);
pause(1)
if rankP<m*n,   % reuse existing problem
fprintf('reusing the existing problem with rank constraint\n')
else
Fmats=1;    %  generate test P0 using sum_ p_iF_i F_i^*   formula for m=n
if Fmats==1&& m~=n,
	fprintf('WARNING!!: generating P0 using sum_pFF^* in file rankrunTPCP.m, BUT m=%i not equal to n=%i  \n',m,n);
	fprintf('resetting m=n to min(m,n) \n')
	if m<n,
		n=m;
	else
		m=n;
	end
end
mn=m*n;
toler=1e-14;   % tolerance for residual in alt proj
fprintf('decreased tolerance for reducing rank is  %g \n',toler)
iterlimit=3500;
unital=1;   
tracepres=1;
%%%%test without knowledge of original rank ?????
P0=[];
fprintf('generating/solving 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)
if Fmats,   % generate Fi's using unitary Fi and prob distrib. p to get Ptr
	fprintf('generating cell F with unitaries for generating P \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);
else
if m==n,
Ptr=zeros(mn);
Ptrr=zeros(mn);
for rr=1:r,
[q,~]=qr(randn(m)+1i*randn(m));
for ii=1:n,
	for jj=ii:n,
		if ii==jj,
			tempii=q(:,ii)*q(:,ii)';
			Ptr((ii-1)*m+1:ii*m,(ii-1)*m+1:ii*m)=tempii;
			%tempii=tempii-diag(diag(tempii));
		else
			tempii=q(:,ii)*q(:,jj)';
			%temp=temp-diag(diag(temp));
			Ptr((ii-1)*m+1:ii*m,(jj-1)*m+1:jj*m)=tempii;
			Ptr((jj-1)*m+1:jj*m,(ii-1)*m+1:ii*m)=tempii';
		end
	end
end
Ptrr=Ptrr+(1/r)*Ptr;
end
Ptr=(Ptrr+Ptrr')/2;
else
[q,~]=qr(randn(m*n)+1i*randn(m*n));
d=10*rand(m*n,1);
Ptr=q*diag(d)*q';
Ptr=(Ptr+Ptr')/2;
end
end
rankPs{1}=rank(Ptr);
% add I for last cell for unital constraint
if unital,
	k=k+1;
end
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,~]=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
     end
end

end   % end for checking if we keep the existing problem


starttotal=cputime;
[Poptfull,rankPs,Fstep]=...
   rankredTPCP(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres,rankPs{1});
Poptfull{1}=Ptr;
%[P,L,R]=altproj(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres);
endtotaltime=cputime-starttotal;
fprintf('total time from rankredTPCP steps is  %g   \n',endtotaltime)



fprintf('initial P and %i optimal Ps are in cell Poptfull \n',(Fstep-1))
fprintf('%i rank values are in cell rankPs \n',(Fstep))



%profile report
