% Skript zum Vortrag "Quantenfehlerkorrektur", Julian Lorenz (lorenz@in.tum.de)
% 
% Hauptseminar "Quantum Computing" im WS 02/03 am Lehrstuhl Prof. Mayr,
% Institut für Informatik der TUM
%
% Dieses Skript demonstriert den 9-Qubit-Code von P. Shor 
%
% Benötigt wird das Softwarepaket von Prof. Bornemann, Zentrum für Mathematik der TUM,
% das zusammen mit anderen Materialien zu dessen Vorlesung "Quantum Computing" 
% unter http://www-m3.mathematik.tu-muenchen.de/m3/teaching/QCOMP0102 verfügbar ist
 
clear all;
SetBasicGates;
NOTNOT = Tensor(X,X);

% Initialisiere zu schützenden Zustand
psi = full(1/sqrt(2)*(Ket ('0') + Ket ('1')))

% Konstruiere Schaltkreis für Kodierung
EncodeCircuit = {
  1    1    H    1    1    w    w    w    w
  w    w    w   CNOT  w    w    w    w    w
  w    w    w    w   CNOT  w    w    w    w
 CNOT  w    H    w    w    1    1    w    w
  w    w    w    w    w   CNOT  w    w    w
  w    w    w    w    w    w   CNOT  w    w
  w   CNOT  H    w    w    w    w    1    1 
  w    w    w    w    w    w    w   CNOT  w
  w    w    w    w    w    w    w    w   CNOT};
  
% Bilde hieraus einen Operator auf 9 Qubit  
U = Circuit2Operator (EncodeCircuit);

% Kodiere den zu schützenden Zustand unter Zuhilfenahme von 8 Ancilla Qubits
psi_encoded = U*Tensor(psi,Ket('00000000'));

% Erzeuge zufälligen unitären Fehler auf einem Qubit (=2x2 unitäre Matrix)
[Error,R] = qr(rand(2,2)+i*rand(2,2));

% Würfle Qubit aus, auf dem der Fehler auftritt
nr = mod(round(rand(1)*9),9)+1;

% ... nimm einen bestimmten Fehler
Error = X;
nr = 7;

%Error
%nr

% Wende den Fehler an
psi_encoded = ApplySingleQubitGate (Error, psi_encoded, nr);

% Führe Fehlerdiagnose durch
% Statt den ganzen Schaltkreis in einen Operator umzuwandeln, was in sehr großen
% Raumdimensionen resultiert, führen wir die Diagnose in Schritten durch
% und messen die Syndrome in Zweier-Packs

psi_encoded = Tensor (Ket('00'), psi_encoded);
psi_encoded = Place (Ctrl(NOT), [3 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [4 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [4 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [5 2], 11) * psi_encoded;
[psi_encoded,synd12] = Measure(psi_encoded,0,2,9);

psi_encoded = Tensor (Ket('00'), psi_encoded);
psi_encoded = Place (Ctrl(NOT), [6 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOTNOT), [7 1 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [8 2], 11) * psi_encoded;
[psi_encoded,synd34] = Measure(psi_encoded,0,2,9);

psi_encoded = Tensor (Ket('00'), psi_encoded);
psi_encoded = Place (Ctrl(NOT), [9 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOTNOT), [10 1 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [11 2], 11) * psi_encoded;
[psi_encoded,synd56] = Measure(psi_encoded,0,2,9);

psi_encoded = Hadamard (psi_encoded);

psi_encoded = Tensor (Ket('00'), psi_encoded);
psi_encoded = Place (Ctrl(NOT), [3 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [4 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [5 1], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOTNOT), [6 1 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOTNOT), [7 1 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOTNOT), [8 1 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [9 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [10 2], 11) * psi_encoded;
psi_encoded = Place (Ctrl(NOT), [11 2], 11) * psi_encoded;
[psi_encoded,synd78] = Measure(psi_encoded,0,2,9);

psi_encoded = Hadamard (psi_encoded);

% Fehlersyndrom nun gemessen
syndrom = [sscanf(dec2bin(synd12,2),'%1d'); sscanf(dec2bin(synd34,2),'%1d'); 
 sscanf(dec2bin(synd56,2),'%1d'); sscanf(dec2bin(synd78,2),'%1d')]'

% Korrigiere nun entsprechend

switch synd12
case 1
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 3);
case 2
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 1);
case 3
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 2);
end

switch synd34
case 1
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 6);
case 2
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 4);
case 3
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 5);
end

switch synd56
case 1
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 9);
case 2
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 7);
case 3
   psi_encoded = ApplySingleQubitGate (X, psi_encoded, 8);
end

switch synd78
case 1
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 7);
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 8);
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 9);
case 2
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 1);
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 2);
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 3);
case 3
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 4);
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 5);
   psi_encoded = ApplySingleQubitGate (Z, psi_encoded, 6);
end

% Dekodiere Zustand (Beachte: U^(-1) = U') und "ignoriere" die 8 Ancilla Qubits
% (werden nicht mehr benötigt) durch Messung
psi_decoded = Measure(U'*psi_encoded,1,8,0)

% Beachte: Übrig bleibt unter Umständen eine unerhebliche globale Phase
% Die Dichtematrizen stimmen aber in jedem Fall überein
Equal(psi_decoded*psi_decoded',psi*psi')