function [Mcellds,Mds,x,y,X,Y]=krausnormalize(Mcell)
% [Mcellds,M,x,y,X,Y]=doublyKraus(Mcell)
%
% Mcellds : data for the Schroedinger-normalized
%                doubly stochastic Kraus map.
%
% Additional variables of possible interest:
%
% x=phihat0, y=phihat1, X=phi0, Y=phi1.
%
[k,~]=size(Mcell);
n=length(cell2mat(Mcell(1)));

% Next we go through phihat0 - phihat1 - phi1 - phi0
% with rho's identities to normalize Mcell
x=eye(n,n);  
tolerence=.001;
test=1;

while test>tolerence,
y=kraus(x,Mcell); Y=inv(y); X=krausdual(Y,Mcell);
xnext=inv(X);
test=hilbertmetric(x,xnext);
x=xnext;
end

% Normalize Mcell
M=cell2mat(Mcell);
invyhalf=inv(sqrtm(y));
Mds=kron(eye(k),invyhalf)*M*sqrtm(x);
Mcellds=mat2cell(Mds,n*ones(k,1),n);
end