%% The example uses a specific choice for hat2phi and phi2hat
% namely,
% phi2hat and hat2phi solve
%
% $\rho = \phi^{1/2} \hat\phi \phi^{1/2}$
%
% in the corresponding direction, e.g., given \rho and \phi
% routine phi2hat computes \phihat, etc.

format compact
clc

disp('******* the example uses n=2, k=3, i.e., a Kraus map with three 2 by 2 coefficient matrices *********')
disp('******* modify below the choice of n,m **************************************************************')
n=2;  % size of matrices
k=3; % number of summands in the Kraus map

% Selecting randomly coefficients for the Kraus map
Mcell=krausdata(n,k);

disp('randomly selected a Kraus map with coeffients [E1; E2; ...]')
E=cell2mat(Mcell)

% Normalize to a doubly stochastic Kraus map
Mcell=krausnormalize(Mcell);

disp('we normalize to a doubly stochastic Kraus map; the coeffients are [F1; F2; ...]')
F=cell2mat(Mcell)
% You can test that this is indeed normalized, by typing:
% kraus(eye(n),Mcell)
% krausdual(eye(n),Mcell)
% both give the identity (or at least very close to it)

% Selecting $\rho_0$, $\rho_1$ randomly
rho0=sort(rand(n,1),'ascend'); rho0=rho0/sum(rho0);
rho1=sort(rand(n,1),'ascend'); rho1=rho1/sum(rho1);

disp('randomly selected'), rhosub0=diag(rho0), clear rho_0
disp('and'), rhosub1=diag(rho1),

%%
% Starting the iteration
% 
% here
% x=$\hat\phi_0$
% y=$\hat\phi_1$
% X=$\phi_0$
% Y=$\phi_1$
x=eye(n);

for i=1:100,
    y=kraus(x,Mcell);
    Y=hat2phi(y,rho1);
    X=krausdual(Y,Mcell);
    xnew=phi2hat(X,rho0);


    x=xnew;
end

phi0hat=x;
phi1hat=y;
phi0=X;
phi1=Y;

%%
% validating (if desired uncomment below)
%
% kraus(phi0hat,Mcell)
% phi1hat
% 
% hat2phi(phi1hat,rho1)
% phi1
% 
% krausdual(phi1,Mcell)
% phi0
% 
% phi2hat(phi0,rho0)
% phi0hat

Mcell_out=mtransform(Mcell,sqrtm(phi0),sqrtm(phi1));
disp('The fixed point of the iteration gives a new Kraus map G=[G1; G2,...]')
G=cell2mat(Mcell_out)

% compare $\rho_1$ with $\cal F(\rho_0)$
disp('Verifying that G(\rho_0)=\rho_1')
rhosub1
G_of_rhosub0=kraus(diag(rho0),Mcell_out)