QCQI — Chapter 10: Entanglement, Entropy & Information Tasks

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

import numpy as np
from numpy.linalg import eigvalsh
def dag(M): return M.conj().T
def ket(i, N): 
    v = np.zeros((N,1), complex); v[i,0]=1; return v
def dm(psi): 
    psi = psi.reshape(-1,1); return psi @ psi.conj().T
def partial_trace(rho, dims, keep):
    # dims: list of subsystem dims; keep: tuple of indices to keep
    import numpy as np
    N = np.prod(dims)
    rho = rho.reshape(*(dims+dims))
    axes = tuple(i for i in range(len(dims)) if i not in keep)
    for ax in sorted(axes, reverse=True):
        rho = np.trace(rho, axis1=ax, axis2=ax+len(dims))
    dkeep = int(np.prod([dims[i] for i in keep]))
    return rho.reshape(dkeep, dkeep)
def vonNeumann_entropy(rho, tol=1e-12):
    vals = np.clip(eigvalsh((rho+dag(rho))/2), 0, 1)
    vals = vals[vals>tol]
    return float(-np.sum(vals*np.log2(vals)))

# Bell |Phi+>
Phi = (1/np.sqrt(2))*np.array([1,0,0,1], complex)
rhoPhi = dm(Phi)
# Product |00>
rho00 = dm(np.array([1,0,0,0], complex))
# Werner: rho = p |Phi+><Phi+| + (1-p) I/4
def werner(p):
    return p*rhoPhi + (1-p)*np.eye(4)/4

def I_mutual(rhoAB):
    rhoA = partial_trace(rhoAB, [2,2], keep=(0,))
    rhoB = partial_trace(rhoAB, [2,2], keep=(1,))
    return vonNeumann_entropy(rhoA)+vonNeumann_entropy(rhoB)-vonNeumann_entropy(rhoAB)
for name, rho in [("Bell", rhoPhi), ("Product", rho00), ("Werner p=0.8", werner(0.8))]:
    SA = vonNeumann_entropy(partial_trace(rho,[2,2],keep=(0,)))
    SB = vonNeumann_entropy(partial_trace(rho,[2,2],keep=(1,)))
    SAB= vonNeumann_entropy(rho)
    print(name, "S(A),S(B),S(AB), I(A:B) ->", round(SA,3), round(SB,3), round(SAB,3), round(I_mutual(rho),3))

sy = np.array([[0,-1j],[1j,0]], complex)
Y = np.kron(sy, sy)
def concurrence(rho):
    R = rho @ Y @ rho.conj() @ Y
    vals = np.sort(np.sqrt(np.real(np.linalg.eigvals(R))))[::-1]
    C = max(0.0, float(vals[0]-vals[1]-vals[2]-vals[3]))
    return C
def EoF_from_C(C):
    x = (1+np.sqrt(1-C**2))/2
    if x <= 0 or x >= 1:
        return 0.0
    h = -x*np.log2(x) - (1-x)*np.log2(1-x)
    return float(h)
for p in [0.0,0.5,0.8,1.0]:
    rho = werner(p)
    C = concurrence(rho)
    print(f"Werner p={p:.1f}: concurrence {C:.3f}, EoF {EoF_from_C(C):.3f}")

def partial_transpose(rho, dims=(2,2), sys=1):
    # transpose on subsystem 'sys' (0 or 1) for 2x2
    rho = rho.reshape(2,2,2,2)
    if sys==0:
        rho = rho.transpose(1,0,2,3)
    else:
        rho = rho.transpose(0,1,3,2)
    return rho.reshape(4,4)
def negativity(rho):
    evals = np.linalg.eigvals(partial_transpose(rho, sys=1))
    return float((np.sum(np.abs(evals)) - 1)/2)
for p in [0.3,0.6,1.0]:
    print(f"Werner p={p:.1f}: negativity {negativity(werner(p)):.3f}")

backend_selected = None
if backend == "cirq":
    import cirq, numpy as np
    q = cirq.LineQubit.range(3)  # 0:Q, 1:A, 2:B
    c = cirq.Circuit()
    # prepare random |ψ> on Q
    c.append(cirq.ry(0.7).on(q[0])); c.append(cirq.rz(1.2).on(q[0]))
    # share Bell on A-B
    c.append([cirq.H(q[1]), cirq.CNOT(q[1], q[2])])
    # Bell-measure Q-A
    c.append(cirq.CNOT(q[0], q[1])); c.append(cirq.H(q[0]))
    c.append(cirq.measure(q[0], key='m0')); c.append(cirq.measure(q[1], key='m1'))
    # classically-controlled correction (simulate by decomposing into branches)
    # For demo, we just append unconditional corrections to show structure
    c.append([cirq.CNOT(q[1], q[2]), cirq.CZ(q[0], q[2])])  # stand-in "frame" idea
    sim = cirq.Simulator()
    res = sim.run(c, repetitions=1)
    print("Cirq teleportation circuit built (Pauli frame placeholder).")
    backend_selected = "cirq"
if backend == "pennylane":
    import pennylane as qml, numpy as np
    dev = qml.device("default.qubit", wires=3, shots=None)
    @qml.qnode(dev)
    def telep():
        qml.RY(0.7, wires=0); qml.RZ(1.2, wires=0)
        qml.Hadamard(1); qml.CNOT(wires=[1,2])
        qml.CNOT(wires=[0,1]); qml.Hadamard(0)
        # emulate perfect classical correction by directly applying controlled Paulis with projectors
        # (omitted for brevity in a no-control-flow setting)
        return qml.state()
    psi = telep(); print("PennyLane teleportation state vector generated.")
    backend_selected = "pennylane"
if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    c = Circuit().ry(0,0.7).rz(0,1.2).h(1).cnot(1,2).cnot(0,1).h(0)
    c = c.measure(0).measure(1)
    dev = LocalSimulator(); _ = dev.run(c, shots=5).result()
    print("Braket teleportation skeleton executed (measurements).")
    backend_selected = "braket"
if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import RY, RZ, H, CNOT, MEASURE
    from pyquil.api import get_qc
    p = Program(); ro = p.declare('ro','BIT',2)
    p += RY(0.7,0); p += RZ(1.2,0); p += H(1); p += CNOT(1,2); p += CNOT(0,1); p += H(0)
    p += MEASURE(0, ro[0]); p += MEASURE(1, ro[1])
    qc = get_qc('3q-qvm'); _ = qc.run(p.wrap_in_numshots_loop(5))
    print("PyQuil teleportation skeleton ran (measurements).")
    backend_selected = "pyquil"
if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    qc = QuantumCircuit(3,2)
    qc.ry(0.7,0); qc.rz(1.2,0)
    qc.h(1); qc.cx(1,2); qc.cx(0,1); qc.h(0)
    qc.measure([0,1],[0,1])
    sim = AerSimulator(); _ = sim.run(qc, shots=8).result()
    print("Qiskit teleportation skeleton executed (measurements).")
    backend_selected = "qiskit"
if backend_selected is None:
    print("Select a valid backend to run the teleportation demo.")

# Similar to teleportation, construct a simple dense-coding skeleton per framework.
# For brevity, we provide one illustrative backend (Cirq) and leave others as exercise hooks.
if backend == "cirq":
    import cirq
    q = cirq.LineQubit.range(2)  # 0:Alice, 1:Bob
    c = cirq.Circuit()
    c.append([cirq.H(q[0]), cirq.CNOT(q[0], q[1])])  # share Bell
    # Encode two bits 'b1 b0' -> apply Z^b1 X^b0 on Alice
    b1,b0 = 1,0
    if b1: c.append(cirq.Z(q[0]))
    if b0: c.append(cirq.X(q[0]))
    # Send qubit 0; Bob performs Bell meas.
    c.append(cirq.CNOT(q[0], q[1])); c.append(cirq.H(q[0]))
    c.append(cirq.measure(q[0], key='b0')); c.append(cirq.measure(q[1], key='b1'))
    sim = cirq.Simulator(); res = sim.run(c, repetitions=1)
    print("Dense coding recovered bits:", int(res.measurements['b1'][0][0]), int(res.measurements['b0'][0][0]))

# Bell basis projectors
Phi = (1/np.sqrt(2))*np.array([1,0,0,1], complex); Psi = (1/np.sqrt(2))*np.array([0,1,1,0], complex)
Phim = (1/np.sqrt(2))*np.array([1,0,0,-1], complex); Psim = (1/np.sqrt(2))*np.array([0,1,-1,0], complex)
P = [dm(Phi), dm(Psi), dm(Phim), dm(Psim)]
def twirl_to_bell_diagonal(rho):
    # Simple average over local Paulis (not exact in code but ok as a pedagogical placeholder)
    return sum(Pi * np.real(np.trace(Pi @ rho)) for Pi in P)
def bbpssw_round(werner_p):
    rho = werner(werner_p)
    rho = twirl_to_bell_diagonal(rho)
    # idealized map of fidelities (textbook formula for Werner -> Werner')
    F = np.real(np.trace(P[0] @ rho))
    Fp = (F**2 + ((1-F)**2)/9) / (F**2 + 2*F*(1-F)/3 + 5*((1-F)**2)/9)
    return F, Fp
for p in [0.6, 0.7, 0.8, 0.9]:
    F,Fp = bbpssw_round(p)
    print(f"Werner p={p:.1f}: Bell fidelity ~{F:.3f} -> after round ~{Fp:.3f}")