QCQI — Chapter 3: Multiple Qubits, Entanglement & Measurement

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

import numpy as np
def bell_phi_plus():
    psi = (1/np.sqrt(2))*np.array([[1,0,0,1]], dtype=complex).T
    return psi / np.linalg.norm(psi)
def rho_from_psi(psi): return psi @ psi.conj().T
def ptrace_rhoA(rho):
    # basis ordering: 00,01,10,11
    R00 = rho[0:2, 0:2]
    R11 = rho[2:4, 2:4]
    return R00 + R11
def von_neumann_entropy(rho, base=2):
    w = np.linalg.eigvalsh(rho)
    w = w[w > 1e-14]
    return float(-(w*np.log(w)/np.log(base)).sum())
def concurrence_pure(a,b,c,d):
    return float(np.clip(2*abs(a*d - b*c), 0, 1))

psi = bell_phi_plus()
rho = rho_from_psi(psi)
rhoA = ptrace_rhoA(rho)
print("rho_A =", np.round(rhoA,3))
print("S(rho_A) =", round(von_neumann_entropy(rhoA), 4))

# Measure along axes in XY-plane with angles alpha (Alice) and beta (Bob) relative to X
def rot_to_z(theta):
    # rotate measurement axis in XY-plane to Z then measure Z (Rz then Ry)
    # For XY-plane axis at angle theta from X, measurement operator is cos(theta) X + sin(theta) Y
    # To map to Z measurement, apply Rz(-theta) then Rx(-pi/2); implement via Ry equivalents for frameworks as needed.
    return theta

if backend == "cirq":
    import cirq, numpy as np
    q0, q1 = cirq.LineQubit.range(2)
    def bell_circuit():
        c = cirq.Circuit()
        c.append([cirq.H(q0), cirq.CNOT(q0, q1)])
        return c
    def measure_axis(circ, a_angle, b_angle, shots=2000):
        c = circ.copy()
        # rotate XY-axis to Z: Rz(-theta) then Rx(-pi/2) maps X to Z; implement with Rx and Rz
        c.append([cirq.rz(-a_angle)(q0), cirq.rx(-np.pi/2)(q0)])
        c.append([cirq.rz(-b_angle)(q1), cirq.rx(-np.pi/2)(q1)])
        c.append([cirq.measure(q0, key='a'), cirq.measure(q1, key='b')])
        sim = cirq.Simulator()
        res = sim.run(c, repetitions=shots)
        A = 1 - 2*res.measurements['a'][:,0]  # map {0,1}->{+1,-1}
        B = 1 - 2*res.measurements['b'][:,0]
        return float(np.mean(A*B))
    circ = bell_circuit()
    # Near-optimal CHSH angles (XY-plane): 0, pi/2 for Alice; pi/4, -pi/4 for Bob
    A0, A1 = 0.0, np.pi/2
    B0, B1 = np.pi/4, -np.pi/4
    E00 = measure_axis(circ, A0, B0)
    E01 = measure_axis(circ, A0, B1)
    E10 = measure_axis(circ, A1, B0)
    E11 = measure_axis(circ, A1, B1)
    S = E00 + E01 + E10 - E11
    print("E00,E01,E10,E11 =", [round(x,3) for x in [E00,E01,E10,E11]])
    print("CHSH S ~", round(S,3))

if backend == "pennylane":
    import pennylane as qml, numpy as np
    dev = qml.device("default.qubit", wires=2, shots=None)
    def measure_corr(alpha, beta):
        @qml.qnode(dev)
        def circuit():
            qml.H(0); qml.CNOT(wires=[0,1])
            # rotate measurement axes in XY plane: observable cos t X + sin t Y
            A = np.cos(alpha)*qml.PauliX(0) + np.sin(alpha)*qml.PauliY(0)
            B = np.cos(beta)*qml.PauliX(1) + np.sin(beta)*qml.PauliY(1)
            return qml.expval(A @ B)
        return circuit()
    A0, A1 = 0.0, np.pi/2
    B0, B1 = np.pi/4, -np.pi/4
    E00 = measure_corr(A0,B0); E01 = measure_corr(A0,B1)
    E10 = measure_corr(A1,B0); E11 = measure_corr(A1,B1)
    S = E00 + E01 + E10 - E11
    print("E00,E01,E10,E11 =", [round(float(x),3) for x in [E00,E01,E10,E11]])
    print("CHSH S ~", round(float(S),3))

if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np
    def bell():
        return Circuit().h(0).cnot(0,1)
    def est(alpha, beta, shots=2000):
        # rotate XY axis to Z: Rz(-theta) then Rx(-pi/2)
        circ = bell().rz(0, float(-alpha)).rx(0, float(-np.pi/2)).rz(1, float(-beta)).rx(1, float(-np.pi/2))
        circ = circ.measure(0).measure(1)
        dev = LocalSimulator()
        res = dev.run(circ, shots=shots).result().measurement_counts
        N = sum(res.values())
        def bit_to_pm1(x): return 1 if x=='0' else -1
        E = 0.0
        for bitstr, count in res.items():
            a = bit_to_pm1(bitstr[-1]); b = bit_to_pm1(bitstr[-2])  # qubit0 -> last char
            E += a*b*count
        return E/N
    A0, A1 = 0.0, np.pi/2
    B0, B1 = np.pi/4, -np.pi/4
    E00 = est(A0,B0); E01 = est(A0,B1); E10 = est(A1,B0); E11 = est(A1,B1)
    S = E00 + E01 + E10 - E11
    print("E00,E01,E10,E11 =", [round(x,3) for x in [E00,E01,E10,E11]])
    print("CHSH S ~", round(S,3))

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, CNOT, RZ, RX, MEASURE
    from pyquil.pyqvm import PyQVM
    import numpy as np
    def est(alpha, beta, shots=2000):
        p = Program()
        ro = p.declare('ro','BIT',2)
        p += H(0); p += CNOT(0,1)
        p += RZ(-alpha, 0); p += RX(-np.pi/2, 0)
        p += RZ(-beta, 1);  p += RX(-np.pi/2, 1)
        p += MEASURE(0, ro[0]); p += MEASURE(1, ro[1])
        qvm = PyQVM(n_qubits = 2)
        qvm.execute(p.wrap_in_numshots_loop(shots))
        res = qvm.read_memory(region_name = 'ro')
        A = 1 - 2*res[:,0]; B = 1 - 2*res[:,1]
        return float(np.mean(A*B))
    A0, A1 = 0.0, np.pi/2
    B0, B1 = np.pi/4, -np.pi/4
    vals = [est(A0,B0), est(A0,B1), est(A1,B0), est(A1,B1)]
    print("E00,E01,E10,E11 =", [round(v,3) for v in vals])
    print("CHSH S ~", round(vals[0]+vals[1]+vals[2]-vals[3],3))

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np
    sim = AerSimulator()
    def est(alpha, beta, shots=2000):
        qc = QuantumCircuit(2,2)
        qc.h(0); qc.cx(0,1)
        qc.rz(-float(alpha), 0); qc.rx(-np.pi/2, 0)
        qc.rz(-float(beta), 1);  qc.rx(-np.pi/2, 1)
        qc.measure([0,1],[0,1])
        res = sim.run(qc, shots=shots).result().get_counts()
        N = sum(res.values())
        E = 0.0
        for bitstr, count in res.items():
            a = 1 if bitstr[-1]=='0' else -1
            b = 1 if bitstr[-2]=='0' else -1
            E += a*b*count
        return E/N
    A0, A1 = 0.0, np.pi/2
    B0, B1 = np.pi/4, -np.pi/4
    E00 = est(A0,B0); E01 = est(A0,B1); E10 = est(A1,B0); E11 = est(A1,B1)
    S = E00 + E01 + E10 - E11
    print("E00,E01,E10,E11 =", [round(x,3) for x in [E00,E01,E10,E11]])
    print("CHSH S ~", round(S,3))