%----------------------------------
%  MATLAB code for channel coding SDP
%  Corresponding to paper: arXiv:1506.08810
%  This code computes an upper bound on the optimal
%  success probability for entanglement assisted 
%  coding over a channel. 
%  This is via the first level SDP relaxation 
%  defined in arXiv:1506.08810.
%----------------------------------

k=2; % number of messages
X=4; % number of input symbols
Y=6; % number of output symbols

% The channel 
% p(x,y) is the probability that the output is y given input x
p = zeros(X,Y);

%
% The channel from the Prevedel et al. paper
%
p = [1/3 0 1/3 0 1/3 0; 1/3 0 0 1/3 0 1/3; 0 1/3 1/3 0 0 1/3; 0 1/3 0 1/3 1/3 0];


%
% Precomputing the constraints and objective
%

% F defines the objective function
F=zeros(1+k*X+k*Y,1+k*X+k*Y);
for i=1:k
    for x=1:X
        for y=1:Y
            F(1+to_num(i,x,k), 1+k*X + to_num(i,y,k)) = p(x,y)*1/k;
            F(1+k*X+ to_num(i,y,k), 1+to_num(i,x,k)) = p(x,y)*1/k;            
        end
    end
end



% The input normalization constraint
input_normalization = zeros(k, 1+k*X + k*Y);
for i=1:k
    for x=1:X
       input_normalization(i,1+to_num(i,x,k)) = 1;
    end
end

% The output normalization constraint
output_normalization = zeros(Y, 1 + k*X + k*Y);
for y=1:Y
    for i=1:k
       output_normalization(y, 1+k*X+to_num(i,y,k)) = 1;
    end
end   


%cvx_solver sedumi

cvx_begin sdp
    variable Omega(1+k*X+k*Y,1+k*X+k*Y) symmetric 
    maximize 1/2*trace(F*Omega)
    subject to
        Omega(1,1) == 1;
        Omega >= 0;
        
        % Entrywise nonnegativity for Omega
        for w=1:(1+k*X+k*Y)
            for wp=1:(1+k*X+k*Y)
                Omega(w,wp) >= 0;
            end
        end
        
        input_normalization*Omega == ones(k,1)*Omega(1,:);
        output_normalization*Omega == ones(Y,1)*Omega(1,:);           
cvx_end

fprintf('The maximum success probability for this channel is at most %d\n', cvx_optval);
