% Quantum State Tomography & Quantum Process Tomography of an LOQC CNOT gate
% by Peter Rohde
% Centre for Quantum Computer Technology
% University of Queensland
% 23rd September, 2004

% clear away all the homeless gutter trash
clear all;

% single qubit input states
H = [1 0 ; 0 0];
V = [0 0 ; 0 1];
D = 1/2 * [1 1 ; 1 1];
A = 1/2 * [1 -1 ; -1 1];
L = 1/2 * [1 i ; -i 1];
R = 1/2 * [1 -i ; i 1];

% single qubit operators 
I = [1 0 ; 0 1];  % Pauli-I
X = [0 1 ; 1 0];  % Pauli-X
Y = [0 -i ; i 0]; % Pauli-Y
Z = [1 0 ; 0 -1]; % Pauli-Z

% CNOT operator
CNOT = [
    1 0 0 0;
    0 1 0 0;
    0 0 0 1;
    0 0 1 0
];

% two qubit tomographic measurement expectation values
%   rows correspond to input settings
%   columns correspond to tomographic settings
%   order: HH,HV,HD,HR,VH,VV,VD,VR,DH,DV,DD,DR,RH,RV,RD,RR
expVals = [0.111111E0,0.E-308,0.555556E-1,0.555556E-1,0.E-308,0.E-308, ...\n  \
0.E-308,0.E-308,0.555556E-1,0.E-308,0.277778E-1,0.277778E-1, ...\n  \
0.555556E-1,0.E-308,0.277778E-1,0.277778E-1;0.E-308,0.111111E0, ...\n  \
0.555556E-1,0.555556E-1,0.E-308,0.E-308,0.E-308,0.E-308,0.E-308, ...\n  \
0.555556E-1,0.277778E-1,0.277778E-1,0.E-308,0.555556E-1, ...\n  \
0.277778E-1,0.277778E-1;0.555556E-1,0.555556E-1,0.111111E0, ...\n  \
0.555556E-1,0.E-308,0.E-308,0.E-308,0.E-308,0.277778E-1, ...\n  \
0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1,0.277778E-1, ...\n  \
0.555556E-1,0.277778E-1;0.555556E-1,0.555556E-1,0.555556E-1, ...\n  \
0.111111E0,0.E-308,0.E-308,0.E-308,0.E-308,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.555556E-1;0.E-308,0.E-308,0.E-308,0.E-308,0.E-308, ...\n  \
0.111111E0,0.555556E-1,0.555556E-1,0.E-308,0.555556E-1, ...\n  \
0.277778E-1,0.277778E-1,0.E-308,0.555556E-1,0.277778E-1, ...\n  \
0.277778E-1;0.E-308,0.E-308,0.E-308,0.E-308,0.111111E0,0.E-308, ...\n  \
0.555556E-1,0.555556E-1,0.555556E-1,0.E-308,0.277778E-1, ...\n  \
0.277778E-1,0.555556E-1,0.E-308,0.277778E-1,0.277778E-1;0.E-308, ...\n  \
0.E-308,0.E-308,0.E-308,0.555556E-1,0.555556E-1,0.111111E0, ...\n  \
0.555556E-1,0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1;0.E-308,0.E-308, ...\n  \
0.E-308,0.E-308,0.555556E-1,0.555556E-1,0.555556E-1,0.E-308, ...\n  \
0.277778E-1,0.277778E-1,0.277778E-1,0.E-308,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.E-308;0.555556E-1,0.E-308,0.277778E-1, ...\n  \
0.277778E-1,0.E-308,0.555556E-1,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.E-308;0.E-308,0.555556E-1,0.277778E-1, ...\n  \
0.277778E-1,0.555556E-1,0.E-308,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1;0.277778E-1,0.277778E-1, ...\n  \
0.555556E-1,0.277778E-1,0.277778E-1,0.277778E-1,0.555556E-1, ...\n  \
0.277778E-1,0.555556E-1,0.555556E-1,0.111111E0,0.555556E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1;0.277778E-1, ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.E-308,0.277778E-1,0.277778E-1,0.555556E-1, ...\n  \
0.277778E-1,0.555556E-1,0.E-308,0.277778E-1,0.277778E-1; ...\n  \
0.555556E-1,0.E-308,0.277778E-1,0.277778E-1,0.E-308,0.555556E-1, ...\n  \
0.277778E-1,0.277778E-1,0.277778E-1,0.277778E-1,0.277778E-1, ...\n  \
0.555556E-1,0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1; ...\n  \
0.E-308,0.555556E-1,0.277778E-1,0.277778E-1,0.555556E-1,0.E-308, ...\n  \
0.277778E-1,0.277778E-1,0.277778E-1,0.277778E-1,0.277778E-1, ...\n  \
0.E-308,0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1; ...\n  \
0.277778E-1,0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1, ...\n  \
0.277778E-1,0.555556E-1,0.277778E-1,0.277778E-1,0.277778E-1, ...\n  \
0.555556E-1,0.277778E-1,0.555556E-1,0.555556E-1,0.111111E0, ...\n  \
0.555556E-1;0.277778E-1,0.277778E-1,0.277778E-1,0.555556E-1, ...\n  \
0.277778E-1,0.277778E-1,0.277778E-1,0.E-308,0.E-308,0.555556E-1, ...\n  \
0.277778E-1,0.277778E-1,0.277778E-1,0.277778E-1,0.555556E-1, ...\n  \
0.277778E-1];

% ------------------------
% QUANTUM STATE TOMOGRAPHY
% ------------------------

% two qubit input states
rho(1,:)  = reshape(transpose(kron(H,H)),1,16); % HH
rho(2,:)  = reshape(transpose(kron(H,V)),1,16); % HV
rho(3,:)  = reshape(transpose(kron(H,D)),1,16); % HD
rho(4,:)  = reshape(transpose(kron(H,R)),1,16); % HR

rho(5,:)  = reshape(transpose(kron(V,H)),1,16); % VH
rho(6,:)  = reshape(transpose(kron(V,V)),1,16); % VV
rho(7,:)  = reshape(transpose(kron(V,D)),1,16); % VD
rho(8,:)  = reshape(transpose(kron(V,R)),1,16); % VR

rho(9,:)  = reshape(transpose(kron(D,H)),1,16); % DH
rho(10,:) = reshape(transpose(kron(D,V)),1,16); % DV
rho(11,:) = reshape(transpose(kron(D,D)),1,16); % DD
rho(12,:) = reshape(transpose(kron(D,R)),1,16); % DR

rho(13,:) = reshape(transpose(kron(R,H)),1,16); % RH
rho(14,:) = reshape(transpose(kron(R,V)),1,16); % RV
rho(15,:) = reshape(transpose(kron(R,D)),1,16); % RD
rho(16,:) = reshape(transpose(kron(R,R)),1,16); % RR

% single qubit observables associated with tomographic measurments
mu0 = H;
mu1 = V;
mu2 = D;
mu3 = R;

% two qubit observables
observable(1,:)  = reshape(transpose(kron(mu0,mu0)),1,16);
observable(2,:)  = reshape(transpose(kron(mu0,mu1)),1,16);
observable(3,:)  = reshape(transpose(kron(mu0,mu2)),1,16);
observable(4,:)  = reshape(transpose(kron(mu0,mu3)),1,16);

observable(5,:)  = reshape(transpose(kron(mu1,mu0)),1,16);
observable(6,:)  = reshape(transpose(kron(mu1,mu1)),1,16);
observable(7,:)  = reshape(transpose(kron(mu1,mu2)),1,16);
observable(8,:)  = reshape(transpose(kron(mu1,mu3)),1,16);

observable(9,:)  = reshape(transpose(kron(mu2,mu0)),1,16);
observable(10,:) = reshape(transpose(kron(mu2,mu1)),1,16);
observable(11,:) = reshape(transpose(kron(mu2,mu2)),1,16); 
observable(12,:) = reshape(transpose(kron(mu2,mu3)),1,16);

observable(13,:) = reshape(transpose(kron(mu3,mu0)),1,16);
observable(14,:) = reshape(transpose(kron(mu3,mu1)),1,16);
observable(15,:) = reshape(transpose(kron(mu3,mu2)),1,16);
observable(16,:) = reshape(transpose(kron(mu3,mu3)),1,16);

% construct ideal output density matrices
for m = 1:16
    rho_M = transpose(reshape(rho(m,:),4,4));
    eps_rho_M = CNOT * rho_M * CNOT;
    ideal_eps_rho(m,:) = reshape(transpose(eps_rho_M),1,16);
end

% transfer matrices (relate measurements to coefficients in density matrices)
%   for input state: p = alpha*|H><H| + beta*|V><V| + gamma*|D><D| + delta*|R><R|
%   tr(p|H><H|) = alpha                     + 1/2 * gamma + 1/2 * delta
%   tr(p|V><V|) = beta                      + 1/2 * gamma + 1/2 * delta
%   tr(p|D><D|) = 1/2 * alpha + 1/2 * beta  + gamma       + 1/2 * delta
%   tr(p|R><R|) = 1/2 * alpha + 1/2 * beta  + 1/2 * gamma + delta
transferOne = [
%   <H> <V> <D> <R>
    1   0   1/2 1/2; % |H><H|
    0   1   1/2 1/2; % |V><V|
    1/2 1/2 1   1/2; % |D><D|
    1/2 1/2 1/2 1;   % |R><R|
];

% transfer matrix for two qubits is given by the square of the single qubit
% matrix where we employ the Kronecker product
transferTwo = kron(transferOne,transferOne);

% coefficients of the observables
obsComps = (inv(transferTwo) * expVals')';

% reconsturct density matrices from basis of observables and their coefficients
actual_eps_rho = obsComps * observable / 4; 
 
% normalise each of the density matrices
for m = 1:16
    actual_eps_rho(m,:) = actual_eps_rho(m,:) / trace(reshape(actual_eps_rho(m,:),4,4));
end    

% --------------------------
% QUANTUM PROCESS TOMOGRAPHY
% --------------------------

% basis of two qubit operators
%pre_mult = CNOT;     % use this pre-multiplier to view the chi-matrix in the CNOT-basis
pre_mult = kron(I,I); % use this pre-multiplier to view the chi-matrix in the Pauli-basis

operator(1,:)  = reshape(transpose(pre_mult * kron(I,I)),1,16); % II
operator(2,:)  = reshape(transpose(pre_mult * kron(I,X)),1,16); % IX
operator(3,:)  = reshape(transpose(pre_mult * kron(I,Y)),1,16); % IY
operator(4,:)  = reshape(transpose(pre_mult * kron(I,Z)),1,16); % IZ

operator(5,:)  = reshape(transpose(pre_mult * kron(X,I)),1,16); % XI
operator(6,:)  = reshape(transpose(pre_mult * kron(X,X)),1,16); % XX
operator(7,:)  = reshape(transpose(pre_mult * kron(X,Y)),1,16); % XY
operator(8,:)  = reshape(transpose(pre_mult * kron(X,Z)),1,16); % XZ

operator(9,:)  = reshape(transpose(pre_mult * kron(Y,I)),1,16); % YI
operator(10,:) = reshape(transpose(pre_mult * kron(Y,X)),1,16); % YX
operator(11,:) = reshape(transpose(pre_mult * kron(Y,Y)),1,16); % YY
operator(12,:) = reshape(transpose(pre_mult * kron(Y,Z)),1,16); % YZ

operator(13,:) = reshape(transpose(pre_mult * kron(Z,I)),1,16); % ZI
operator(14,:) = reshape(transpose(pre_mult * kron(Z,X)),1,16); % ZX
operator(15,:) = reshape(transpose(pre_mult * kron(Z,Y)),1,16); % ZY
operator(16,:) = reshape(transpose(pre_mult * kron(Z,Z)),1,16); % ZZ

% verify normalisation condition for the basis operators
opNorm = zeros(4,4);
for m = 1:16
    opNorm = opNorm + transpose(reshape(operator(m,:),4,4)) * transpose(reshape(operator(m,:),4,4))';    
end

% determine lambda matrix
ideal_lambda  = transpose(inv(transpose(rho)) * transpose(ideal_eps_rho));
actual_lambda = transpose(inv(transpose(rho)) * transpose(actual_eps_rho));

% determine beta, kappa and chi matrices
ideal_chi  = zeros(17,17);
actual_chi = zeros(17,17);
for m = 1:16
    for n = 1:16
        ERE = [];
        for j = 1:16
            Em = transpose(reshape(operator(m,:),4,4));
            En = transpose(reshape(operator(n,:),4,4));
            rhoj = transpose(reshape(rho(j,:),4,4));
            ERE(j,:) = reshape(transpose(Em * rhoj * En'),1,16);
        end
        % beta matrix
        beta(:,:,m,n)   = inv(transpose(rho)) * transpose(ERE);
        
        % kappa matrix
        kappa(:,:,m,n)  = inv(beta(:,:,m,n));
        
        % ideal chi-matrix
        ideal_chi(m,n)  = sum(sum(kappa(:,:,m,n) .* ideal_lambda));
        
        % actual chi-matrix
        actual_chi(m,n) = sum(sum(kappa(:,:,m,n) .* actual_lambda));
    end
end

% normalise chi matrices (i.e. must have trace 1)
ideal_chi  = ideal_chi  / trace(ideal_chi);
actual_chi = actual_chi / trace(actual_chi);

% verify normalisation of process matrices
idealNorm  = 0;
actualNorm = 0;
for m = 1:16
    for n = 1:16
        idealNorm  = idealNorm  + ideal_chi(m,n)  * transpose(reshape(operator(m,:),4,4)) * conj(reshape(operator(n,:),4,4));
        actualNorm = actualNorm + actual_chi(m,n) * transpose(reshape(operator(m,:),4,4)) * conj(reshape(operator(n,:),4,4));    
    end
end

% process fidelity
Fp = trace(real(ideal_chi) * real(actual_chi));

% plot real part of ideal chi-matrix
subplot(2,2,1);
pcolor(real(ideal_chi)); axis square; colorbar;
title('Re(\chi_{ideal})');
set(gca,'XTick',[1:16],'YTick',[1:16]);
set(gca,'XTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});
set(gca,'YTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});

% plot imaginary ideal chi-matrix
subplot(2,2,2);
pcolor(imag(ideal_chi)); axis square; colorbar;
title('Im(\chi_{ideal})');
set(gca,'XTick',[1:16],'YTick',[1:16]);
set(gca,'XTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});
set(gca,'YTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});

% plot real part of actual chi-matrix
subplot(2,2,3);
pcolor(real(actual_chi)); axis square; colorbar;
title('Re(\chi_{actual})');
set(gca,'XTick',[1:16],'YTick',[1:16]);
set(gca,'XTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});
set(gca,'YTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});

% plot imaginary part of actual chi-matrix
subplot(2,2,4);
pcolor(imag(actual_chi)); axis square; colorbar;
title('Im(\chi_{actual})');
set(gca,'XTick',[1:16],'YTick',[1:16]);
set(gca,'XTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});
set(gca,'YTickLabel',{'II';'IX';'IY';'IZ';'XI';'XX';'XY';'XZ';'YI';'YX';'YY';'YZ';'ZI';'ZX';'ZY';'ZZ'});

colormap gray;