function varargout = chanconv(E,from,varargin)

% CHANCONV   Convert between representations of a quantum channel
% requires: TrX, schmidt, purification
% author: Toby Cubitt
% license: GPL2
%
%    S = chanconv(E,TO) converts the quantum channel E to the TO
%    representation, determining the input representation automatically.
%
%    S = chanconv(E,FROM,TO) converts the quantum channel E in the FROM
%    representation to the TO representation.
%
%    S = chnconv(E,TO,DIM) converts the quantum channel E to the TO
%    representation, determining the input representation
%    automatically. The input and output dimensions DIMA and DIMB of the
%    channel are specified by the vector DIM = [DIMA, DIMB].
%
%    S = chnconv(E,FROM,TO,DIM) converts the quantum channel E in the
%    FROM representation to the TO representation. The input and output
%    dimensions DIMA and DIMB of the channel are specified by the vector
%    DIM = [DIMA, DIMB].
%
%    [S,dim] = chancon(E, ...) returns the input, output and environment
%    dimensions of the channel in a vector DIM = [DIMA,DIMB,DIME].
%
%    FROM and TO must be one of:
%      'kraus' - a cell array containing a set of Kraus operators
%      'choi'  - a Choi-Jamiolkowski state representative density matrix
%      'linop' - a matrix representation of the linear operator
%      'isom'  - an isometry matrix
%      'state' - a purification of the Choi-Jamiolkowski state
%      'func'  - a function on density matrices which implements a channel
%      'pauli' - a Pauli representation of a qubit channel
%                (FIXME: not implemented)
%
%    The Kraus operators are the operators A_k acting from the *left* in
%    the Kraus decomposition (i.e. the action of the channel is given by
%    sum_k A_k rho A_k').
%
%    The entries of a Choi-Jamiolkowski state representative matrix
%    should be ordered so that the input system comes first, the output
%    second. I.e. it should be the state resulting from applying the
%    channel to the *second* half of an entangled state. Non-standard
%    Choi-Jamiolkowski states are allowed as input representations; there
%    is no requirement that the reduced state on the first subsystem be
%    the identity (i.e. it can result from applying the channel to the
%    first half of any full-Schmidt-rank entangled state, which need not
%    necessarily be maximally entangled).
%
%    Output Choi-Jamiolkowski state representatives follow the same
%    ordering convention, but the reduced state on the first subsystem is
%    always the identity (i.e. they are always the state resulting from
%    applying the channel to the second half of a canonical maximally
%    entangled state).
%
%    The entries of a state vector representation of a channel should be
%    ordered so that the subsystems occur in the following order: input
%    system, output system, environment. Ouput state vectors follow the
%    same ordering convention.
%
%    It is possible for the same matrix to simultaneously be a valid
%    Choi-Jamiolkowski state representative, a valid matrix
%    representation of the linear operator, and even a valid isometry.
%    Therefore the input representation can not always be inferred
%    unambiguously, and a reasonable guess must be made. The input
%    representation is assumed to be a Choi-Jamiolkowski state
%    representative if E is an (unnormalized) density matrix (i.e. is
%    Hermitian and positive, to within a tolerance of 10E-10). Otherwise
%    it is assumed to be a matrix representation of the linear
%    operator. However, to be absolutely sure the correct input
%    representation is used, you should specify it explicitly.
%
%    See also chan2R, randChan and applychan.


%% Copyright (C) 2004-2010, 2012 Toby Cubitt
%%
%% 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.



%----- initialization -----

% sort out function arguments
if length(varargin) == 2
  to = varargin{1};
  dim = varargin{2};
elseif length(varargin) == 1
  if  ischar(varargin{1});
    to = varargin{1};
  else
    dim = varargin{1};
    to = from;
    from = '';
  end
else
  to = from;
  from = '';
end


% if input representation is not explicitly specified, infer it from channel
if strcmp(from,'')
  if iscell(E)
    from = 'kraus';
  elseif isa(E,'function_handle')
    from = 'func';
  elseif isnumeric(E)
    if size(E,2) == 1
      from = 'state';
    elseif size(E,1) ~= size(E,2)
      from = 'isom';
    elseif norm(E - E') < 10E-10 && all(real(eig(E)) >= 0)
      warning(['Can not unambiguously infer input representation; guessing Choi-Jamiolkowski state representative']);
      from = 'choi';
    else
      warning(['Can not unambiguously infer input representation; guessing linear operator matrix representation']);
      from = 'linop';
    end
  else
    error('E is not a valid representation of a channel');
  end
end

% check FROM and TO representation arguments
if ~(strcmp(to,'choi') || strcmp(to,'linop') ...
     || strcmp(to,'kraus') || strcmp(to,'isom') ...
     || strcmp(to,'state') || strcmp(to,'func')) ...
   || ~(strcmp(from,'choi') || strcmp(from,'linop') ...
     || strcmp(from,'kraus') || strcmp(from,'isom') ...
     || strcmp(from,'state') || strcmp(from,'func'))
  error(['FROM and TO arguments must be one of: ',...
         'choi, kraus, linop, isom, state or func']);
end


% determine dimensions
if strcmp(from,'choi')
  if exist('dim')
    dimA = dim(1);
    dimB = dim(2);
  else
    dimA = sqrt(size(E,1));  % assume dimA = dimB
    dimB = dimA;
  end
elseif strcmp(from,'linop')
  if exist('dim')
    dimA = dim(1);
    dimB = dim(2);
  else
    dimA = sqrt(size(E,2));
    dimB = sqrt(size(E,1));
  end
elseif strcmp(from,'kraus')
  dimA = size(E{1},2);
  dimB = size(E{1},1);
  dimE = length(E);
elseif strcmp(from,'isom')
  dimA = size(E,2);
  if exist('dim')
    dimB = dim(2);
  else
    dimB = dimA;  % assume dimB = dimA
  end
  dimE = size(E,1)/dimB;
elseif strcmp(from,'state')
  if exist('dim');
    dimA = dim(1);
    dimB = dim(2);
    dimE = length(E)/(dimA*dimB);
  else
    error('State input representation requires explicit DIM argument');
  end
elseif strcmp(from,'func')
  if exist('dim');
    dimA = dim(1);
    if length(dim) > 1
      dimB = dim(2);
    end
  else
    error(['Function input representation requires explicit DIM argument']);
  end
end

% check consistency of dimensions
if mod(dimA,1) ~= 0 || (exist('dimB') && mod(dimB,1)) ~= 0 || ...
      (exist('dimE') && mod(dimE,1) ~= 0)
  error(['Dimensions could not be inferred automatically, ', ...
         'or E is not a valid representation of a channel']);
end
if exist('dim') && (dimA ~= dim(1) || (exist('dimB') && dimB ~= dim(2)))
  error(['DIM argument does not match channel dimensions']);
end




%----- algorithms -----

% set ouput to input, in case we're converting from and to the same
% representation
S = E;


%--- kraus ---
if strcmp(from,'kraus')
  if strcmp(to,'func')
    % convert from kraus to func by building a function that calls 'applychan'
    S = @(rho) applychan(S,rho,'kraus',[dimA,dimB]);

  elseif strcmp(to,'isom') || strcmp(to,'state')
    % convert to isometry by tensoring with projectors on computational
    % basis states of the environment; if converting to state, convert to
    % isometry and leave rest till later
    S = zeros(dimB*dimE,dimA);
    for j = 1:dimE
      S([j:dimE:dimB*dimE],:) = E{j};
    end
    from = 'isom';
    E = S;

  elseif strcmp(to,'choi') || strcmp(to,'kraus')
    % convert from kraus to choi by applying the channel to second half
    % of a maximally entangled state; if converting from kraus to kraus,
    % convert to choi and convert back later, to produce a minimal kraus
    % representation
    w = zeros(dimA^2,1);
    w([1:dimA+1:dimA^2]) = 1;
    w = w*w';
    S = zeros(dimA*dimB);
    for j = 1:dimE
      S = S + kron(eye(dimA),E{j})*w*kron(eye(dimA),E{j}');
    end
    from = 'choi';
    E = S;

  elseif strcmp(to,'linop')
    % convert kraus to linop by summing all tensor products of kraus
    % operators with their complex conjugates; if converting to anything
    % else, convert to linop and leave rest till later
    S = zeros(dimB^2,dimA^2);
    for j = 1:length(E)
      S = S + kron(conj(E{j}),E{j});
    end
    from = 'linop';
    E = S;
  end
end



%--- func ---
if strcmp(from,'func') && ~strcmp(to,'func')
  % construct linear operator by applying function to complete operator
  % basis, and leave rest till later
  clear('S');
  for j = 1:dimA
    for k = 1:dimA
      e = zeros(dimA);
      e(k,j) = 1;
      e = E(e);
      S(:,dimA*(j-1)+k) = reshape(e,[prod(size(e)),1]);
    end
  end

  % calculate / check consistency of output dimension
  if ~exist('dimB')
    dimB = sqrt(size(S,1));
    if mod(dimB,1) ~= 0
      error(['function E does not implement a channel']);
    end
  elseif dimB ~= sqrt(size(S,1))
    error(['DIM argument does not match channel output dimension']);
  end

  from = 'linop';
  E = S;
end



%--- linop ---
if strcmp(from,'linop') && ~strcmp(to,'linop')
  if strcmp(to,'func')
    % convert from linop to func by building a function that calls 'applychan'
    S = @(rho) applychan(S,rho,'linop',[dimA,dimB]);

  else
    % if converting to anything else, convert to choi by "gamma" involution and
    % leave rest till later (note that because of the order in which matlab
    % takes elements when reshaping, the inner reshape gives the multi-indices
    % in reverse order, i.e. E_{j,i,l,k} *not* E_{i,j,k,l}, which needs to be
    % converted to E_{jl,ik} for a choi matrix whose first subsystem is the
    % input, but because of the outer reshape we have to permute to
    % E_{l,j,k,i})
    S = reshape(permute(reshape(E,[dimB,dimB,dimA,dimA]),[1,3,2,4]),...
                dimA*dimB,dimA*dimB);
    from = 'choi';
    E = S;
  end
end


%--- choi ---
if strcmp(from,'choi')
  % convert choi to standard Choi-Jamiolkowski form, i.e. ensure reduced
  % state on first system is identity
  s = kron(TrX(E,2,[dimA,dimB])^(-1/2),eye(dimB));
  S = s*E*s;
  S = dimA*S/trace(S);
  E = S;

  if strcmp(to,'linop')
    % convert from choi to linop by "gamma" involution (E_{j,i,l,k} from inner
    % reshape becomes E_{i,k,j,l}, which the outer reshape converts to
    % E_{ki,lj})
    S = reshape(permute(reshape(E,[dimB,dimA,dimB,dimA]),[1,3,2,4]),...
                dimB^2,dimA^2);

  elseif strcmp(to,'func')
    % convert from choi to func by building a function that calls 'applychan'
    S = @(rho) applychan(S,rho,'choi',[dimA,dimB]);

  elseif ~strcmp(to,'choi')
    % convert from choi to state by taking purification; if converting to
    % anything else, convert to state and leave rest till later
    [S,dimE] = purification(E);
    from = 'state';
    E = S;
  end
end



%--- state ---
% FIXME: should convert state to "standard" state a la choi to choi
if strcmp(from,'state') && ~strcmp(to,'state')
  if strcmp(to,'func')
    % convert isometry to func by building function that calls 'applychan'
    S = @(rho) applychan(S,rho,'state',[dimA,dimB]);

  else
    % convert state to isometry by finding Schmidt decomposition for A|BE
    % partition, and "reversing" the ket on input system into a bra; if
    % converting to anything else, convert to isometry and leave rest till
    % later
    [U,D,V] = schmidt(E,[dimA,dimB*dimE]);
    S = V*diag(D)*U.';
    from = 'isom';
    E = S;
  end
end



%--- isom ---
if strcmp(from,'isom')
  if strcmp(to,'kraus')
    % convert isometry to kraus by left-multiplying by computational basis
    % state environment bras
    S = cell(dimE,1);
    for j = 1:dimE
      S{j} = E([j:dimE:dimB*dimE],:);
    end

  elseif strcmp(to,'choi') || strcmp(to,'linop')
    % convert isometry to choi by applying isometry to second half of a
    % maximally entangled state and tracing over environment
    w = zeros(dimA^2,1);
    w([1:dimA+1:dimA^2]) = 1;
    w = w*w';
    S = TrX(kron(eye(dimA),E)*w*kron(eye(dimA),E'),2,[dimA*dimB,dimE]);

    % convert choi to linop by "gamma" involution
    if strcmp(to,'linop')
      S = reshape(permute(reshape(S,[dimB,dimA,dimB,dimA]),[1,3,2,4]),...
                  [dimB^2,dimA^2]);
    end

  elseif strcmp(to,'state')
    % convert isometry to state by "reversing" the bra on the input
    % system into a ket
    S = reshape(E,dimA*dimB*dimE,1);

  elseif strcmp(to,'func')
    % convert isometry to func by building function that calls 'applychan'
    S = @(rho) applychan(S,rho,'isom',[dimA,dimB]);
  end
end


% output results
varargout{1} = S;
if nargout > 1;
   varargout{2} = [dimA,dimB,dimE];
end