#!/usr/bin/env python3
"""
extension_benchmarks.py
=======================
Deterministic extension benchmarks accompanying the revised technical edition of
"First-Principles Validation and Extension of the S-I-C-T Category-Theoretic
Reconstruction" (Roth Complexity Lab working paper, 21 October 2026).

Executable research source code. Run locally with:

    python -m pip install numpy
    python extension_benchmarks.py

Requires Python 3.10+ and NumPy. The deterministic assertions verify:
  * the reversible cyclic-shift target and the irreversible lossy simulator,
  * the six-step logical-erasure budget (Theorem 3.13 / Table 7),
  * the finite-factor permutation theorem (Theorem 6.1),
  * the near-sharp reachable-orbit Euler error ratio (Section 6.3).

The finite-state results (entropy, erasure, support) are exact and
architecture-independent. Floating-point ratios can show last-digit noise
across runtimes; the assertion tolerance reflects that.
"""
import json
import math

RESULTS = {}


# ---------------------------------------------------------------------------
# 1. Reversible target vs lossy simulator on {0,1}^6  (Section 6.2, Table 7)
# ---------------------------------------------------------------------------
N = 6
SIZE = 1 << N  # 64 states


def shift(state):
    """Reversible cyclic shift F(z)_i = z_{i+1 mod N} (a permutation)."""
    return tuple(state[(i + 1) % N] for i in range(N))


def lossy_shift(state):
    """Apply the reversible shift, then clear the last register bit (non-injective)."""
    s = shift(state)
    return s[:-1] + (0,)


def shannon_entropy_bits(dist):
    """H(mu) in bits for a probability distribution (dict state -> prob)."""
    h = 0.0
    for p in dist.values():
        if p > 0:
            h -= p * math.log2(p)
    return h


# All 64 six-bit states, uniform initial ensemble.
states = [tuple(int(b) for b in format(i, f"0{N}b")) for i in range(SIZE)]
mu = {s: 1.0 / SIZE for s in states}

# Verify the target is a permutation (bijection).
target_images = [shift(s) for s in states]
assert len(set(target_images)) == SIZE, "target shift must be a permutation"

# Verify the simulator is non-injective (merges states).
sim_images = [lossy_shift(s) for s in states]
assert len(set(sim_images)) < SIZE, "lossy simulator must merge states"

# Evolve the ensemble under the lossy simulator and record the entropy budget.
entropy_series = [shannon_entropy_bits(mu)]
erasure_series = []
support_series = [sum(1 for p in mu.values() if p > 0)]

for _step in range(1, 11):
    new_mu = {}
    for s, p in mu.items():
        img = lossy_shift(s)
        new_mu[img] = new_mu.get(img, 0.0) + p
    erasure_series.append(entropy_series[-1] - shannon_entropy_bits(new_mu))
    mu = new_mu
    entropy_series.append(shannon_entropy_bits(mu))
    support_series.append(sum(1 for p in mu.values() if p > 0))

# Table 7 expected values.
expected_entropy = [6, 5, 4, 3, 2, 1, 0, 0, 0, 0, 0]
expected_erasure = [1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
expected_support = [64, 32, 16, 8, 4, 2, 1, 1, 1, 1, 1]

for i, (e, ee) in enumerate(zip(entropy_series, expected_entropy)):
    assert abs(e - ee) < 1e-12, f"entropy step {i}: {e} != {ee}"
for i, (e, ee) in enumerate(zip(erasure_series, expected_erasure)):
    assert abs(e - ee) < 1e-12, f"erasure step {i}: {e} != {ee}"
for i, (s, es) in enumerate(zip(support_series, expected_support)):
    assert s == es, f"support step {i}: {s} != {es}"

total_erased = sum(erasure_series)
assert abs(total_erased - 6.0) < 1e-12, "total logical erasure must be 6 bits"

RESULTS["corrected_six_bit_benchmark"] = {
    "N": N,
    "description": "Reversible cyclic-shift target with a lossy (bit-clearing) simulator on {0,1}^6",
    "entropy_bits": entropy_series,
    "erasure_bits": erasure_series,
    "support_size": support_series,
    "total_erased_bits": total_erased,
}


# ---------------------------------------------------------------------------
# 2. Finite-factor permutation theorem (Theorem 6.1)
# ---------------------------------------------------------------------------
# If F: A->A is surjective, h: A->Z is surjective onto a finite Z, and
# h o F = G o h, then G must be a permutation of Z.
def check_finite_factor_theorem():
    # A = {0,...,7}, F(x) = (x+1) mod 8 (surjective / bijective).
    A = list(range(8))
    def F(x):
        return (x + 1) % 8
    # h: A -> Z = {0,1,2,3} surjective, h(x) = x mod 4.
    Z = list(range(4))
    def h(x):
        return x % 4
    # Induced G: for each z, pick a with h(a)=z, then G(z) = h(F(a)).
    G = {}
    for z in Z:
        a = next(a for a in A if h(a) == z)
        G[z] = h(F(a))
    # Verify the factor relation h o F = G o h on all of A.
    for a in A:
        assert h(F(a)) == G[h(a)], "factor relation h o F = G o h must hold"
    # G must be a permutation of the finite set Z.
    images = [G[z] for z in Z]
    assert sorted(images) == sorted(Z), "G must be a permutation of Z"
    return True


assert check_finite_factor_theorem()
RESULTS["finite_factor_theorem"] = {
    "verified": True,
    "statement": "A finite exact surjective factor of a surjective dynamics is a permutation (Theorem 6.1).",
}


# ---------------------------------------------------------------------------
# 3. Near-sharp reachable-orbit Euler error ratio (Section 6.3)
# ---------------------------------------------------------------------------
# Harmonic oscillator, dt = 0.01, identified with a complex number.
# f = exp(-i dt) (exact rotation), g = 1 - i dt (explicit Euler), z0 = 1.
dt = 0.01
f = complex(math.cos(dt), -math.sin(dt))
g = complex(1.0, -dt)
n = 200000

actual_error = abs(g ** n - f ** n)
# Path-weighted reachable bound B_n = |g - f| * sum_{k=0}^{n-1} |g|^k
# = |g - f| * (|g|^n - 1) / (|g| - 1)   for |g| != 1.
abs_g = abs(g)
geom_sum = (abs_g ** n - 1.0) / (abs_g - 1.0)
bound = abs(g - f) * geom_sum
ratio = actual_error / bound if bound else float("nan")

RESULTS["reachable_orbit_euler_ratio"] = {
    "dt": dt,
    "n": n,
    "actual_error": actual_error,
    "reachable_bound": bound,
    "ratio": ratio,
}

# The paper reports ratio = 0.9999778799 at n = 200000 (near-sharp, not equality).
assert 0.9999 < ratio < 1.0001, f"ratio {ratio} not near-sharp"


# ---------------------------------------------------------------------------
# Emit results
# ---------------------------------------------------------------------------
print(json.dumps(RESULTS, indent=2))
with open("extension_results.json", "w") as fh:
    json.dump(RESULTS, fh, indent=2)

print("\nAll extension assertions passed.")