Listing: Listing

#!/usr/bin/env python3
"""Hawk02-B: optional diagnostic artifact; unexecuted in manuscript."""
import datetime
import hashlib
import itertools
import json
from pathlib import Path
import sys
import numpy as np

D, N = 2, 3
OMEGA = 1.0
THETA = np.full(N, np.pi / 4)
PERTURBATION = 1e-4
TOL = 2e-10

def adj(A):
    return A.conj().T

def opnorm(A):
    return float(np.linalg.norm(A, 2))

def td(A, B):
    X = (A - B + adj(A - B)) / 2
    return float(np.sum(np.abs(np.linalg.eigvalsh(X))) / 2)

def exp_hermitian(H):
    w, V = np.linalg.eigh(H)
    return (V * np.exp(-1j * w)) @ adj(V)

def channel(Ks, X):
    out = np.zeros((Ks[0].shape[0], Ks[0].shape[0]), complex)
    for A in Ks:
        out += A @ X @ adj(A)
    return out

def matrix_units(d):
    for a in range(d):
        for b in range(d):
            X = np.zeros((d, d), complex)
            X[a, b] = 1
            yield a, b, X

def referenced(chan, d):
    # Normalized Choi state: channel applied to half of Phi_d.
    sample = chan(np.eye(d) / d)
    out = np.zeros((sample.shape[0] * d,) * 2, complex)
    for a, b, X in matrix_units(d):
        out += np.kron(chan(X), X) / d
    return out

def bell(d):
    v = np.eye(d).reshape(-1) / np.sqrt(d)
    return np.outer(v, v.conj())

def check(name, condition, manifest):
    passed = bool(condition)
    manifest["checks"][name] = passed
    if not passed:
        raise AssertionError(name)

def diagnostics(root, manifest):
    d, n = D, N
    r, m, z = 1 + n * d, d + 1, d * (n + 1)
    Id = np.eye(d)
    sigma = Id / d

    # One-excitation sector ordering: memory, then radiation bins.
    S = []
    U = np.eye(z, dtype=complex)
    Ua = np.eye(z, dtype=complex)
    detuning = np.zeros((z, z), complex)
    detuning[:d, :d] = np.diag(np.linspace(-1, 1, d))

    for j, theta in enumerate(THETA):
        A = np.zeros((z, z), complex)
        for a in range(d):
            b = d + j * d + a
            A[a, b] = A[b, a] = 1
        S.append(A)
        P = A @ A
        step = np.eye(z) + (np.cos(theta) - 1) * P
        step = step - 1j * np.sin(theta) * A
        U = step @ U
        actual = exp_hermitian(theta * A +
                               PERTURBATION * detuning)
        Ua = actual @ Ua

    E = np.zeros((z, d), complex)
    E[:d] = Id
    T = np.zeros((m * r, z), complex)
    for a in range(d):
        T[(a + 1) * r, a] = 1
    for j in range(n):
        for a in range(d):
            T[1 + j * d + a, d + j * d + a] = 1
    W, Wa = T @ U @ E, T @ Ua @ E

    survive = 1.0
    alpha = []
    for theta in THETA:
        alpha.append(-1j * np.sin(theta) * survive)
        survive *= np.cos(theta)
    alpha = np.asarray(alpha)
    q, p = float(survive**2), float(1 - survive**2)
    cd = 1 - 1 / d**2

    JC = np.zeros((m, d), complex)
    JC[1:] = Id
    vac = np.zeros((r, 1), complex)
    vac[0, 0] = 1
    ground = np.zeros((m, 1), complex)
    ground[0, 0] = 1
    J = np.zeros((r, d), complex)
    for j in range(n):
        J[1 + j*d:1 + (j+1)*d] = alpha[j] * Id / np.sqrt(p)
    Wan = np.sqrt(q) * np.kron(JC, vac)
    Wan += np.sqrt(p) * np.kron(ground, J)

    check("isometry", opnorm(adj(W) @ W - Id) < TOL, manifest)
    check("Hamiltonian_to_isometry", opnorm(W - Wan) < TOL,
          manifest)
    check("amplitude_normalization",
          abs(np.vdot(alpha, alpha).real - p) < TOL, manifest)

    # Accessible and residual Kraus operators.
    tensor = W.reshape(m, r, d)
    Kr = [tensor[a, :, :] for a in range(m)]
    Km = [tensor[:, b, :] for b in range(r)]

    PJ = J @ adj(J)
    def decode(Y):
        return adj(J) @ Y @ J + np.trace((np.eye(r)-PJ) @ Y)*sigma

    # Explicit decoder Kraus matrices.
    vals, vecs = np.linalg.eigh(np.eye(r) - PJ)
    Kd = [adj(J)]
    for k in np.flatnonzero(vals > 0.5):
        for a in range(d):
            B = np.zeros((d, r), complex)
            B[a, :] = vecs[:, k].conj() / np.sqrt(d)
            Kd.append(B)

    check("decoder_Kraus_complete",
          opnorm(sum(adj(A) @ A for A in Kd)-np.eye(r)) < TOL,
          manifest)

    for a, b, X in matrix_units(d):
        rad = channel(Kr, X)
        mem = channel(Km, X)
        target_r = p * J @ X @ adj(J)
        target_r += q * np.trace(X) * (vac @ adj(vac))
        target_m = q * JC @ X @ adj(JC)
        target_m += p * np.trace(X) * (ground @ adj(ground))
        check(f"radiation_unit_{a}_{b}",
              opnorm(rad-target_r) < TOL, manifest)
        check(f"residual_unit_{a}_{b}",
              opnorm(mem-target_m) < TOL, manifest)
        check(f"decoded_unit_{a}_{b}",
              opnorm(decode(rad)-(p*X+q*np.trace(X)*sigma)) < TOL,
              manifest)

    tau = p * ground @ adj(ground) + q * JC @ sigma @ adj(JC)
    rho_dec = referenced(lambda X: decode(channel(Kr, X)), d)
    rho_mem = referenced(lambda X: channel(Km, X), d)
    check("recovery_Bell_witness",
          abs(td(rho_dec, bell(d))-q*cd) < TOL, manifest)
    check("leakage_Bell_witness",
          abs(td(rho_mem, np.kron(tau, sigma))-q*cd) < TOL,
          manifest)

    Hm = np.diag([0] + [OMEGA]*d)
    Hr = np.diag([0] + [OMEGA]*(n*d))
    Htotal = np.kron(Hm, np.eye(r)) + np.kron(np.eye(m), Hr)
    check("prepared_energy_sector",
          opnorm(adj(W) @ Htotal @ W - OMEGA*Id) < TOL,
          manifest)
    check("perturbed_isometry",
          opnorm(adj(Wa) @ Wa - Id) < TOL, manifest)
    check("Duhamel_operator_allowance",
          opnorm(Wa-W) <= n*abs(PERTURBATION)+TOL, manifest)

    # Time pinching: preserve the internal labels in each bin.
    blocks = [slice(0, 1)]
    blocks += [slice(1+j*d, 1+(j+1)*d) for j in range(n)]
    def time_pinch(Y):
        out = np.zeros_like(Y)
        for block in blocks:
            out[block, block] = Y[block, block]
        return out

    def time_decode(Y):
        out = Y[0, 0] * sigma
        for j in range(n):
            block = slice(1+j*d, 1+(j+1)*d)
            out = out + Y[block, block]
        return out

    coherent_weight = float(np.sum(np.abs(alpha)**4) / p)
    for a, b, X in matrix_units(d):
        Y = time_pinch(channel(Kr, X))
        check(f"time_preserves_code_{a}_{b}",
              opnorm(time_decode(Y) -
                     (p*X+q*np.trace(X)*sigma)) < TOL, manifest)
        check(f"fixed_mode_after_time_pinch_{a}_{b}",
              opnorm(decode(Y) -
                     (coherent_weight*X +
                      (1-coherent_weight)*np.trace(X)*sigma))
              < TOL, manifest)

    # Internal-label dephasing.
    label_chan = lambda X: decode(channel(Kr, np.diag(np.diag(X))))
    label_error = 1 - p/d - q/d**2
    check("label_dephasing_Bell_witness",
          abs(td(referenced(label_chan, d), bell(d))-label_error)
          < TOL, manifest)

    # Keep only the first emitted cell: CP loss map on compressed R.
    keep = np.zeros((r, r), complex)
    keep[:1+d, :1+d] = np.eye(1+d)
    Kloss = [keep]
    for b in range(1+d, r):
        B = np.zeros((r, r), complex)
        B[0, b] = 1
        Kloss.append(B)
    Jfirst = np.zeros((r, d), complex)
    Jfirst[1:1+d] = Id
    def decode_first(Y):
        Pfirst = Jfirst @ adj(Jfirst)
        return (adj(Jfirst) @ Y @ Jfirst +
                np.trace((np.eye(r)-Pfirst) @ Y)*sigma)
    collected = lambda X: decode_first(
        channel(Kloss, channel(Kr, X)))
    pfirst = float(abs(alpha[0])**2)
    check("incomplete_collection_Bell_witness",
          abs(td(referenced(collected, d), bell(d)) -
              (1-pfirst)*cd) < TOL, manifest)

    # Product decoder, K=2: witness for the analytic many-cell optimum.
    phase = np.diag(np.exp(2j*np.pi*np.arange(d)/d))
    shift = np.roll(np.eye(d), 1, axis=0)
    single = [np.sqrt(p)*Id]
    for a in range(d):
        for b in range(d):
            single.append(np.sqrt(q)/d *
                          np.linalg.matrix_power(shift, a) @
                          np.linalg.matrix_power(phase, b))
    product = [np.kron(A, B) for A, B in
               itertools.product(single, repeat=2)]
    target_many = 1 - (1-cd*q)**2
    check("two_cell_Bell_witness",
          abs(td(referenced(lambda X: channel(product, X), d*d),
                 bell(d*d))-target_many) < TOL, manifest)

    arrays = {
        "S_sector": np.stack(S), "W_ideal": W, "W_actual": Wa,
        "J_radiation": J, "J_memory": JC,
        "radiation_Kraus": np.stack(Kr),
        "memory_Kraus": np.stack(Km),
        "decoder_Kraus": np.stack(Kd), "tau_memory": tau,
        "alpha": alpha, "theta": THETA,
        "H_memory": Hm, "H_radiation_compressed": Hr
    }
    np.savez(root/"matrices.npz", **arrays)
    manifest["analytic_targets"] = {
        "q": q, "single_recovery": q*cd,
        "single_leakage": q*cd, "two_cell_recovery": target_many,
        "detuning_Duhamel_allowance": n*abs(PERTURBATION)
    }

def main():
    if len(sys.argv) != 2:
        raise SystemExit("usage: script.py NEW_OUTPUT_DIRECTORY")
    root = Path(sys.argv[1])
    root.mkdir(parents=True, exist_ok=False)
    manifest = {
        "status": "RUNNING",
        "time_utc": datetime.datetime.now(
            datetime.timezone.utc).isoformat(),
        "python": sys.version,
        "numpy": np.__version__,
        "parameters": {
            "d": D, "N": N, "omega": OMEGA,
            "angles": THETA.tolist(), "event_duration": 1,
            "perturbation": PERTURBATION, "tolerance": TOL
        },
        "checks": {},
        "diamond_optimization_executed": False,
        "proof_machine_verified": False,
        "latex_compilation_verified": False
    }
    source = Path(__file__).read_bytes()
    manifest["script_sha256"] = hashlib.sha256(source).hexdigest()
    try:
        diagnostics(root, manifest)
        manifest["status"] = "OPTIONAL_DIAGNOSTICS_COMPLETED"
    except BaseException as exc:
        manifest["status"] = "FAILED"
        manifest["failure"] = repr(exc)
        raise
    finally:
        (root/"manifest.json").write_text(
            json.dumps(manifest, indent=2, allow_nan=False))
    print(json.dumps(manifest, indent=2, allow_nan=False))

if __name__ == "__main__":
    main()