QCQI — Chapter 9: Quantum Error Correction & Noise (Expanded)

# Backend selector
backend = "cirq"  # change: "pennylane", "braket", "pyquil", "qiskit"
print("Selected backend:", backend)
# Optional installs:
# !pip install cirq pennylane amazon-braket-sdk pyquil qiskit qiskit-aer numpy

import numpy as np, random, math
I = np.eye(2, dtype=complex)
X = np.array([[0,1],[1,0]], complex)
Y = np.array([[0,-1j],[1j,0]], complex)
Z = np.array([[1,0],[0,-1]], complex)
def apply_kraus(rho, Ks):
    out = np.zeros_like(rho, dtype=complex)
    for K in Ks: out += K @ rho @ K.conj().T
    return out
def depolarizing(p):
    return [np.sqrt(1-p)*I, np.sqrt(p/3)*X, np.sqrt(p/3)*Y, np.sqrt(p/3)*Z]
def phase_damping(lmb):
    K0 = np.sqrt(1-lmb) * I
    K1 = np.sqrt(lmb) * np.diag([1,0])
    K2 = np.sqrt(lmb) * np.diag([0,1])
    return [K0, K1, K2]
def amplitude_damping(g):
    E0 = np.array([[1,0],[0,np.sqrt(1-g)]], complex)
    E1 = np.array([[0,np.sqrt(g)],[0,0]], complex)
    return [E0,E1]

import numpy as np, random, math
def majority_decode(bits): return 1 if bits.count(1)>=2 else 0
def sample_logical_error_rate(p, trials=3000):
    err=0
    for _ in range(trials):
        phys=[0,0,0]
        flips=[random.random()<p for _ in range(3)]
        for i,f in enumerate(flips):
            if f: phys[i]^=1
        dec = majority_decode(phys)
        if dec!=0: err+=1
    return err/trials
for p in [0.01,0.05,0.1,0.2]:
    print(p, "-> logical error ≈", round(sample_logical_error_rate(p,1500),4), "  (theory ~ 3p^2 - 2p^3)")

if backend == "cirq":
    import cirq, numpy as np, random
    q = cirq.LineQubit.range(5)  # 0..2 data, 3..4 ancilla
    c = cirq.Circuit()
    c.append([cirq.CNOT(q[0],q[1]), cirq.CNOT(q[0],q[2])])
    # inject X errors with probability p
    p=0.12
    for j in range(3):
        if random.random()<p: c.append(cirq.X(q[j]))
    # Z1Z2 on ancilla 3; Z2Z3 on ancilla 4
    c.append([cirq.CNOT(q[0],q[3]), cirq.CNOT(q[1],q[3]), cirq.measure(q[3], key='s12')])
    c.append([cirq.CNOT(q[1],q[4]), cirq.CNOT(q[2],q[4]), cirq.measure(q[4], key='s23')])
    res = cirq.Simulator().run(c, repetitions=1)
    s12 = int(res.measurements['s12'][0][0]); s23=int(res.measurements['s23'][0][0])
    print("Cirq syndrome (s12,s23)=", (s12,s23), " -> map (01->X1, 11->X2, 10->X3)")
elif backend == "pennylane":
    import pennylane as qml, numpy as np, random
    dev = qml.device("default.qubit", wires=5, shots=1)
    @qml.qnode(dev)
    def circ(p=0.12):
        qml.CNOT(wires=[0,1]); qml.CNOT(wires=[0,2])
        for j in [0,1,2]:
            if np.random.rand()<p: qml.PauliX(j)
        qml.CNOT(wires=[0,3]); qml.CNOT(wires=[1,3])
        qml.CNOT(wires=[1,4]); qml.CNOT(wires=[2,4])
        return qml.sample(qml.PauliZ(3)), qml.sample(qml.PauliZ(4))
    s12,s23 = circ(); s12=int((1-s12)/2); s23=int((1-s23)/2)
    print("PennyLane syndrome:", (s12,s23))
elif backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np, random
    p=0.12
    c = Circuit().cnot(0,1).cnot(0,2)
    for j in [0,1,2]:
        if np.random.rand()<p: c = c.x(j)
    c = c.cnot(0,3).cnot(1,3).cnot(1,4).cnot(2,4).measure(3).measure(4)
    res = LocalSimulator().run(c, shots=50).result().measurement_counts
    print("Braket syndrome counts:", res)
elif backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import CNOT, X, MEASURE
    from pyquil.api import get_qc
    import numpy as np
    p=0.12
    prog = Program(); ro = prog.declare('ro','BIT',2)
    prog += CNOT(0,1); prog += CNOT(0,2)
    for j in [0,1,2]:
        if np.random.rand()<p: prog += X(j)
    prog += CNOT(0,3); prog += CNOT(1,3); prog += MEASURE(3,ro[0])
    prog += CNOT(1,4); prog += CNOT(2,4); prog += MEASURE(4,ro[1])
    counts = get_qc('5q-qvm').run(prog.wrap_in_numshots_loop(50))
    bits = counts.get_register_map()['ro']
    print("PyQuil first 5 samples:", bits[:5])
elif backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np
    p=0.12
    qc = QuantumCircuit(5,2)
    qc.cx(0,1); qc.cx(0,2)
    for j in [0,1,2]:
        if np.random.rand()<p: qc.x(j)
    qc.cx(0,3); qc.cx(1,3); qc.measure(3,0)
    qc.cx(1,4); qc.cx(2,4); qc.measure(4,1)
    counts = AerSimulator().run(qc, shots=256).result().get_counts()
    print("Qiskit syndrome counts:", counts)

import numpy as np
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], complex)
HXH = H @ np.array([[0,1],[1,0]], complex) @ H
print("H X H =", np.round(HXH,3), " (≈ Z)")

# Binary symplectic form for quick syndrome calc on Pauli errors
import numpy as np
# For repetition-3 with checks Z1Z2, Z2Z3, represent checks as rows of [X|Z] (6 cols)
# Z1Z2 -> x-part 000, z-part 110 ; Z2Z3 -> x 000, z 011
H = np.array([[0,0,0, 1,1,0],
              [0,0,0, 0,1,1]], dtype=int)
# Error e as [x|z] length 6, e.g., X on qubit 2 -> x=010, z=000
def syndrome(e):
    e = np.array(e, dtype=int).reshape(6)
    x = e[:3]; z = e[3:]
    # commutation s = H_X z + H_Z x (mod 2)
    s = (H[:,:3] @ z + H[:,3:] @ x) % 2
    return tuple(int(v) for v in s)
# Test: X errors on each qubit
for q in range(3):
    e = [0,0,0,0,0,0]; e[q]=1
    print("X on",q+1,"-> syndrome", syndrome(e))