%% function [cvx_optval,X,W,Y] = RobCPTNI(J,dim)
% computes the CPTNI robustness R of the Hermiticity-preserving map J
% with input/output dimensions optionally provided in dim=[dimA,dimB]

function [cvx_optval,X,W,Y] = RobCPTNI(J,dim)
    
J = (J+J')/2; % make sure that J is Hermitian, avoid numerical issues

if nargin<2
    d = max(size(J));
    dim = [1,1]*sqrt(d);
else
    d = dim(1)*dim(2);
end
 
cvx_begin sdp quiet

   variable X(d,d) hermitian semidefinite
   variable lam nonnegative
   dual variable Z
   dual variable Y
   
   PartialTrace(X,2,dim) <= lam*eye(dim(1));
   Y: PartialTrace(J+X,2,dim) <= (1+lam)*eye(dim(1));
   Z: J + X >=0
   
   minimize lam
   
cvx_end

W = kron(Y,eye(dim(2)))-Z; % notation from the paper
   
end