%% QUANTUMSEDUMI     Quantum SDP front-end for SeDuMi
%  requires: SeDuMi (http://sedumi.ie.lehigh.edu/)
%  author: Nathaniel Johnston
%  license: GPL2
%
%  [PRIM,DUAL,Y] = QUANTUMSEDUMI(A,B,PHIA,PHIB) returns the primal and
%  dual optimal values, as well as the dual optimal operator Y of the
%  "quantum" semidefinite program associated with the primal operator A,
%  the dual operator B, and the Hermicity-preserving map defined by the
%  left and right Choi-Kraus operators given by the arrays PHIA and PHIB,
%  respectively.
%
%  An optional fifth argument, PARS, may be supplied to QUANTUMSEDUMI,
%  which performs the same functions as the PARS argument for SeDuMi
%  itself. See the SeDuMi documentation to learn more.
%
%  See [1] and its appendix for a description of this "quantum"
%  semidefinite program form and how its conversion to the standard
%  semidefinite program form takes place.
%
%  You must have SeDuMi installed for this script to function. Use SeDuMi
%  1.1 if at all possible, since there is a bug in SeDuMi 1.21 and 1.3
%  when dealing with complex matrices that makes it fail sometimes.
%
%  [1] N. Johnston and D. W. Kribs. A family of norms with applications in
%  quantum information theory II. Quantum Inf. Comput., 11:104-123, 2011.
%
%  http://www.nathanieljohnston.com/quantumsedumi/

%% Copyright (c) 2009-2011 Nathaniel Johnston
% 
%  This program is free software; you can redistribute it and/or
%  modify it under the terms of the GNU General Public License
%  as published by the Free Software Foundation; either version 2
%  of the License, or (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
%  MA 02110-1301, USA.

%% Code starts here
function [prim,dual,Y] = QuantumSeDuMi(A,B,PhiA,PhiB,pars)

% check user input for errors
if (nargin < 4)
    error('QuantumSeDuMi requires at least four input arguments: A,B,PhiA,PhiB');
elseif (nargin < 5)
    pars.fid = 1;   % use noisy SeDuMi output by default
end

sA = size(A);
sB = size(B);
sPA = size(PhiA);
sPB = size(PhiB);

% the following error checks are provided for ease of use, but should be
% removed if speed is extremely important

% BEGIN REMOVABLE ERROR CHECKS
if (sA(1) ~= sA(2))
    error('A must be a square matrix');
elseif (sB(1) ~= sB(2))
    error('B must be a square matrix');
elseif (sPA(1) ~= sB(1) || sPA(2) ~= sA(1))
    error('PhiA size mismatch');
elseif (sPB(1) ~= sA(1) || sPB(2) ~= sB(1))
    error('PhiB size mismatch');
elseif (sPA(3) ~= sPB(3))
    error('PhiA and PhiB must have the same number of Kraus operators');
end
% END OF REMOVABLE ERROR CHECKS

len_arr = [size(A,1),size(B,1)];
max_len = sum(len_arr);
num_kraus = size(PhiA,3);

% pad the problem appropriately with zeros and add an appropriate negative
% identity Kraus operator to make this a legitimate SDP
newB = zeros(max_len);
newB(1:len_arr(2),1:len_arr(2)) = B;
newPhiA = zeros(max_len,len_arr(1),num_kraus+1);
newPhiA(1:len_arr(2),:,1:num_kraus) = PhiA;
newPhiA((len_arr(2)+1):max_len,:,num_kraus+1) = -eye(len_arr(1));
newPhiB = zeros(len_arr(1),max_len,num_kraus+1);
newPhiB(:,1:len_arr(2),1:num_kraus) = PhiB;
newPhiB(:,(len_arr(2)+1):max_len,num_kraus+1) = eye(len_arr(1));

B = newB;
PhiA = newPhiA;
PhiB = newPhiB;

% vectorize the input matrices A and B
tot_len = max_len*max_len;
b = reshape(A,prod(sA),1);
c = reshape(B,tot_len,1);

% vectorize the Kraus operators
F = zeros(max_len,max_len,len_arr(1)*len_arr(1));
for i = 1:max_len
    for j = 1:max_len
        for k = 0:(len_arr(1)*len_arr(1)-1)
            F(i,j,k+1) = permute(PhiB(mod(k,len_arr(1))+1,i,:),[1,3,2])*permute(PhiA(j,floor(k/len_arr(1))+1,:),[3,1,2]);
        end
    end
end

% set up equality constraints
At = zeros(len_arr(1)*len_arr(1),tot_len);
for i = 1:(len_arr(1)*len_arr(1))
    At(i,:) = reshape(F(:,:,i),tot_len,1)';
end

% remove potential linear dependancies in the constraints
[As,r] = GaussElim([At b]);     % faster than MATLAB's rref function
At = As(1:r,1:(size(As,2)-1));
b = As(1:r,size(As,2));

% we want to work in the cone of PSD complex matrices
K.s = max_len;

% If the script halts with the following error message...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ??? Error using ==> minus
% Matrix dimensions must agree.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ...then you are experiencing a bug in SeDuMi 1.21 or 1.3 that affects
% problems with complex-valued data. For a potential workaround, see
% http://sedumi.ie.lehigh.edu/index.php?option=com_kunena&Itemid=78&func=view&catid=9&id=4623
% Try to restate your problem using all real variables if possible. The
% following four lines are the cause of the problem.
if(~isreal(A) || ~isreal(B) || ~isreal(PhiA) || ~isreal(PhiB))
    K.scomplex = 1;
    K.ycomplex = 1:length(b);  % uncomment once sedumi problem is figured out
end

% solve the SDP
if pars.fid == 1
    disp('QuantumSeDuMi front-end by Nathaniel Johnston, 2009-2011.');
end
[x,y] = sedumi(At,b,c,K,pars);
dual = real(x'*c);
prim = real(y'*b);
Y = reshape(x,size(B));
Y = Y(1:sB(1),1:sB(2));


% Gaussian elimination helper function since we need to remove linear
% dependancies before passing the constraints to SeDuMi. This function is
% much faster than MATLAB's built-in rref function for large systems. The
% output argument A is the reduced row echelon form of the input argument
% A. The output argument R is the rank of A.
function [A,r] = GaussElim(A)
[m,n] = size(A);
tol = max([m,n])*eps(size(A))*norm(A,'inf');
i = 1;
j = 1;
while (i <= m) && (j <= n)
   [p,k] = max(abs(A(i:m,j)));
   k = k+i-1;
   if (p <= tol)
      A(i:m,j) = zeros(m-i+1,1);
      j = j + 1;
   else
      A([i k],j:n) = A([k i],j:n);
      Ai = A(i,j:n)/A(i,j);    
      A(:,j:n) = A(:,j:n) - A(:,j)*Ai;
      A(i,j:n) = Ai;
      i = i + 1;
      j = j + 1;
   end
end
r = i - 1;