QCQI — Chapter 2: Qubits, States & the Bloch Sphere

# 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
theta = np.pi * 0.37   # polar angle
phi   = np.pi * 0.73   # azimuth
shots = 2000
print("theta, phi =", theta, phi, "| shots =", shots)

import numpy as np
X = np.array([[0,1],[1,0]], dtype=complex)
Y = np.array([[0,-1j],[1j,0]], dtype=complex)
Z = np.array([[1,0],[0,-1]], dtype=complex)
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], dtype=complex)
S = np.array([[1,0],[0,1j]], dtype=complex)
Sd = np.array([[1,0],[0,-1j]], dtype=complex)
def Ry(a): return np.array([[np.cos(a/2), -np.sin(a/2)], [np.sin(a/2), np.cos(a/2)]], dtype=complex)
def Rz(a): return np.array([[np.exp(-1j*a/2), 0], [0, np.exp(1j*a/2)]], dtype=complex)
def bloch_from_angles(theta, phi):
    return np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
print("Bloch vector (theory):", np.round(bloch_from_angles(theta, phi), 4))

# Prepare |psi> = Ry(theta) Rz(phi) |0>  (global phase ignored)
ket0 = np.array([[1],[0]], dtype=complex)
psi = Ry(theta) @ Rz(phi) @ ket0
rho = psi @ psi.conj().T
ex = float(np.real(np.trace(rho @ X)))
ey = float(np.real(np.trace(rho @ Y)))
ez = float(np.real(np.trace(rho @ Z)))
print("<X>,<Y>,<Z> (analytic) =", np.round([ex,ey,ez], 4))

if backend == "cirq":
    import cirq, numpy as np, collections
    q = cirq.LineQubit(0)
    def prep():
        return [cirq.ry(theta)(q), cirq.rz(phi)(q)]
    def measure_along(axis):
        circuit = cirq.Circuit()
        circuit.append(prep())
        if axis == 'x':
            circuit.append(cirq.H(q))
        elif axis == 'y':
            circuit.append(cirq.MatrixGate(np.array([[1,0],[0,-1j]])).on(q))  # S^\dagger
            circuit.append(cirq.H(q))
        circuit.append(cirq.measure(q, key='m'))
        sim = cirq.Simulator()
        res = sim.run(circuit, repetitions=shots)
        h = res.histogram(key='m')
        N0 = h.get(0,0); N1 = h.get(1,0)
        return (N0 - N1) / shots
    print({k: round(measure_along(k),4) for k in ['x','y','z']})

if backend == "pennylane":
    import pennylane as qml, numpy as np
    dev = qml.device("default.qubit", wires=1, shots=shots)
    @qml.qnode(dev)
    def expvals():
        qml.RY(theta, wires=0); qml.RZ(phi, wires=0)
        return qml.expval(qml.PauliX(0)), qml.expval(qml.PauliY(0)), qml.expval(qml.PauliZ(0))
    print(np.round(expvals(), 4))

if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np
    def counts(axis):
        circ = Circuit().ry(0, float(theta)).rz(0, float(phi))
        if axis == 'x':
            circ.h(0)
        elif axis == 'y':
            circ.phaseshift(0, -np.pi/2).h(0)   # S^\dagger then H
        circ = circ.measure(0)
        dev = LocalSimulator()
        res = dev.run(circ, shots=shots).result()
        c = res.measurement_counts
        N0 = c.get('0', 0); N1 = c.get('1', 0)
        return (N0 - N1) / shots
    print({k: round(counts(k), 4) for k in ['x','y','z']})

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import RY, RZ, H as H_gate, PHASE, MEASURE
    from pyquil.api import get_qc
    import numpy as np
    def est(axis):
        p = Program()
        ro = p.declare('ro', 'BIT', 1)
        p += RY(theta, 0); p += RZ(phi, 0)
        if axis == 'x':
            p += H_gate(0)
        elif axis == 'y':
            p += PHASE(-np.pi/2, 0); p += H_gate(0)  # S^\dagger then H
        p += MEASURE(0, ro[0])
        qc = get_qc('1q-qvm')
        res = qc.run(p.wrap_in_numshots_loop(shots))
        bits = res.get_register_map()['ro']
        N1 = int((bits == 1).sum()); N0 = shots - N1
        return (N0 - N1) / shots
    print({k: round(est(k), 4) for k in ['x','y','z']})

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np
    sim = AerSimulator()
    def est(axis):
        qc = QuantumCircuit(1,1)
        qc.ry(float(theta), 0); qc.rz(float(phi), 0)
        if axis == 'x':
            qc.h(0)
        elif axis == 'y':
            qc.sdg(0); qc.h(0)
        qc.measure(0,0)
        res = sim.run(qc, shots=shots).result().get_counts()
        N0 = res.get('0',0); N1 = res.get('1',0)
        return (N0 - N1) / shots
    print({k: round(est(k), 4) for k in ['x','y','z']})

# Show that adding Rz(gamma) before *immediate Z-measure* doesn't change counts,
# but it changes interference if followed by H.
gamma = np.pi/3
# Analytic check:
psi0 = np.array([[1],[0]], dtype=complex)
A = Rz(gamma) @ psi0
pZ_direct = [abs(A[0,0])**2, abs(A[1,0])**2]
# Now add H after Rz:
B = H @ (Rz(gamma) @ (H @ psi0))  # H,Rz,H around |0> ~ phase->relative
pZ_afterH = np.abs(B.flatten())**2
print("Z after Rz only:", np.round(pZ_direct,4), " | Z after H+Rz+H:", np.round(pZ_afterH,4))