Listing: Listing

#!/usr/bin/env python3
"""Hawk02-C finite-proxy diagnostics. New, unexecuted artifact."""
import datetime
import hashlib
import json
import math
from pathlib import Path
import sys
import numpy as np

D = 2
C = 1.0
OMEGA = 3.0
TIME = 1.0
COLLECTED = 8
REFERENCE = 24
TOL = 2e-10

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

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

def chain(n):
    H = OMEGA * np.eye(n + 1)
    for j in range(n):
        H[j, j + 1] = H[j + 1, j] = C
    values, vectors = np.linalg.eigh(H)
    initial = np.zeros(n + 1)
    initial[0] = 1.0
    psi = vectors @ (
        np.exp(-1j * TIME * values) * (adj(vectors) @ initial))
    return H, psi

def tail_bound(n):
    # Bounds 2 * sum_{m=n+1}^\infty z^m/m! by its first
    # term and a geometric tail. Defaults meet the ratio condition.
    z = 2 * C * TIME
    m = n + 1
    ratio = z / (m + 1)
    if ratio >= 1:
        return 1.0
    first = math.exp(m * math.log(z) - math.lgamma(m + 1))
    return min(1.0, 2 * first / (1 - ratio))

def full_isometry(psi):
    n = len(psi) - 1
    rad = 1 + n * D
    W = np.zeros(((D + 1) * rad, D), complex)
    for a in range(D):
        W[(a + 1) * rad, a] = psi[0]
        for j in range(1, n + 1):
            W[1 + (j - 1) * D + a, a] = psi[j]
    return W

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

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

def main():
    if len(sys.argv) != 2:
        raise SystemExit("usage: script.py NEW_OUTPUT_DIRECTORY")
    if not (1 <= COLLECTED < REFERENCE):
        raise ValueError("require 1 <= COLLECTED < REFERENCE")
    root = Path(sys.argv[1])
    root.mkdir(parents=True, exist_ok=False)
    report = {
        "status": "RUNNING",
        "model": "FINITE_CHAIN_PROXY_NOT_CONTINUUM",
        "time_utc": datetime.datetime.now(
            datetime.timezone.utc).isoformat(),
        "python": sys.version,
        "numpy": np.__version__,
        "script_sha256": hashlib.sha256(
            Path(__file__).read_bytes()).hexdigest(),
        "parameters": {
            "d": D, "c": C, "omega0": OMEGA, "T": TIME,
            "collected_packets": COLLECTED,
            "reference_packets": REFERENCE, "tolerance": TOL
        },
        "checks": {},
        "diamond_optimization_executed": False,
        "proof_machine_verified": False,
        "latex_compilation_verified": False
    }

    def check(name, residual, bound=TOL):
        value = float(residual)
        limit = float(bound)
        passed = math.isfinite(value) and value <= limit
        report["checks"][name] = {
            "residual": value, "bound": limit, "passed": passed
        }
        (root / "manifest.json").write_text(
            json.dumps(report, indent=2, allow_nan=False))
        if not passed:
            raise AssertionError(name)

    try:
        Hr, pr = chain(REFERENCE)
        Hn, pn = chain(COLLECTED)
        W = full_isometry(pr)
        padded = np.zeros(REFERENCE + 1, complex)
        padded[:COLLECTED + 1] = pn
        Wn = full_isometry(padded)
        check("reference_isometry", opnorm(adj(W) @ W - np.eye(D)))
        check("finite_predictive_comparison",
              opnorm(W - Wn), tail_bound(COLLECTED) + TOL)

        q = float(abs(pr[0])**2)
        ell = float(np.sum(abs(pr[COLLECTED + 1:])**2))
        success = float(np.sum(abs(pr[1:COLLECTED + 1])**2))
        retained = q + ell
        check("finite_proxy_weight_sum",
              abs(q + ell + success - 1.0))

        # Complete collection split: residual B first, accessible A second.
        A_dim = 1 + COLLECTED * D
        miss_dim = 1 + (REFERENCE - COLLECTED) * D
        B_dim = (D + 1) * miss_dim
        L = np.zeros((A_dim, D), complex)
        R = np.zeros((B_dim, D), complex)
        avac = np.zeros((A_dim, 1), complex)
        bvac = np.zeros((B_dim, 1), complex)
        avac[0, 0] = bvac[0, 0] = 1.0
        for a in range(D):
            R[(a + 1) * miss_dim, a] = pr[0]
            for j in range(1, COLLECTED + 1):
                L[1 + (j - 1) * D + a, a] = pr[j]
            for j in range(COLLECTED + 1, REFERENCE + 1):
                k = j - COLLECTED - 1
                R[1 + k * D + a, a] = pr[j]

        Wsplit = np.kron(R, avac) + np.kron(bvac, L)
        check("complete_split_isometry",
              opnorm(adj(Wsplit) @ Wsplit - np.eye(D)))
        tensor = Wsplit.reshape(B_dim, A_dim, D)
        Ka = [tensor[b, :, :] for b in range(B_dim)]
        Kb = [tensor[:, a, :] for a in range(A_dim)]
        J = L / math.sqrt(success)
        PJ = J @ adj(J)
        sigma = np.eye(D) / D

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

        for a, b, X in units(D):
            target_a = (L @ X @ adj(L)
                        + retained * np.trace(X) * avac @ adj(avac))
            target_b = (R @ X @ adj(R)
                        + success * np.trace(X) * bvac @ adj(bvac))
            actual_a = channel(Ka, X)
            actual_b = channel(Kb, X)
            check(f"finite_proxy_accessible_{a}_{b}",
                  opnorm(actual_a - target_a))
            check(f"finite_proxy_joint_residual_{a}_{b}",
                  opnorm(actual_b - target_b))
            check(f"finite_proxy_decoded_{a}_{b}",
                  opnorm(decode(actual_a)
                         - (success * X
                            + retained * np.trace(X) * sigma)))

        # Complete packet-normalized proxy comparison.
        projected = pr.copy()
        projected[COLLECTED + 1:] = 0
        projected[1:COLLECTED + 1] *= math.sqrt((1 - q) / success)
        What = full_isometry(projected)
        h = q + math.sqrt((1 - q) * success)
        check("finite_proxy_scalar_overlap",
              opnorm(adj(W) @ What - h * np.eye(D)))
        check("finite_proxy_packet_norm_bound",
              opnorm(W - What), math.sqrt(2 * ell) + TOL)

        z = 2 * C * TIME
        eta = z**(COLLECTED + 1) / (
            2 * math.factorial(COLLECTED + 1))
        report["analytical_continuum_allowances"] = {
            "reference_complete_bound": tail_bound(REFERENCE),
            "predictive_complete_bound": tail_bound(COLLECTED),
            "packet_missed_weight_upper_bound": eta**2,
            "survival_upper_bound": min(
                1.0, 4 / (math.pi**2 * C**2 * TIME**2))
        }
        report["finite_proxy_values"] = {
            "q": q, "ell": ell, "success": success,
            "recovery_and_leakage_formula": (1 - D**-2) * retained,
            "packet_complete_formula": math.sqrt(
                max(0.0, 1 - h**2))
        }
        np.savez(root / "matrices.npz",
                 H_reference=Hr, H_finite=Hn,
                 psi_reference=pr, psi_finite=pn,
                 W_reference=W, W_finite=Wn,
                 collected_embedding=L, joint_residual_embedding=R)
        report["status"] = "FINITE_PROXY_DIAGNOSTICS_COMPLETED"
    except BaseException as exc:
        report["status"] = "FAILED"
        report["failure"] = repr(exc)
        raise
    finally:
        (root / "manifest.json").write_text(
            json.dumps(report, indent=2, allow_nan=False))
    print(json.dumps(report, indent=2, allow_nan=False))

if __name__ == "__main__":
    main()