function Y=krausdual(X,Mcell)
% Y=krausdual(X,Mcell)
% 
% Y = Sum_1^k M_i^* X M_i
%
% X,Y are observables, Mi are square n by n
[k,m]=size(Mcell); % m is supposed to be 1
n=length(cell2mat(Mcell(1)));

Y=zeros(n,n);
for i=1:k,
    Mi=cell2mat(Mcell(i));
    Y=Y+Mi'*X*Mi;
end
end