"""
ontological_codex.py — the SHD-CCP Ontological Codex: a three-scale quaternionic
delta-chain memory offload (720 → 72 → 10) with exact telescoping, SLERP
reconstruction, holonomy control, and a logged delta stream you can replay to
re-create the structure.

This is the executable realization of the design grounded in five papers:

  • Torsional Markov Memory-Offload Pump  — THE blueprint. Three coupled gears:
      outer 720-frame (12-gear), middle 72-frame INVERSE (10-gear), inner 10-frame
      (8-gear). Coarse = ordered product of deltas (Prop 2, telescoping, EXACT);
      fine = fractional-power SLERP lift (Prop 4, exact on the principal branch);
      72→10 down-sampling partitions into blocks of size 7 or 8 (Prop 3); the
      gear cycle current J>0 (Lemma 5) and holonomy H=Q8·Q10·Q12≈1 keep it coherent.
  • Hyperbolic Sparsemax Pump — the SOURCE. A rational, parity-safe kernel on the
      64-voxel (Z/4Z)^3 lattice emits unit-quaternion increments. Rational spinors
      ⇒ the offload telescopes EXACTLY (0 error), not just ~1e-15. The rational
      forward path is kept separate from the transcendental SLERP/holonomy
      reconstruction layer (the paper's octet-residual quarantine).
  • Golden Strassen-Clifford / Harmonic / TQF — governance + honesty. Dual-Helix:
      the centroid is an untrusted lossy SUMMARY (Helix A); the delta chain is the
      EXACT reconstruction (Helix B). We do NOT import their factoring-acceleration
      claims (measured ε ≤ 0 elsewhere; no sub-linear jump-ahead).

Sources are selectable:  'sparsemax' (rational, parity-safe) | 'smooth' | 'random'.

Pure Python 3 standard library.   python3 ontological_codex.py
"""

import math
from fractions import Fraction as F

N_OUTER = 720          # outer (12-gear) frames
GROUP_OM = 10          # 10 outer deltas per middle block  -> 72 middle frames
N_MIDDLE = N_OUTER // GROUP_OM
N_INNER = 10           # inner (8-gear) frames; 72 -> 10 via variable 7/8 blocks

# --------------------------------------------------------------------------- #
# Quaternion ops — components may be Fraction (exact) or float
# --------------------------------------------------------------------------- #
def qmul(a, b):
    w1, x1, y1, z1 = a; w2, x2, y2, z2 = b
    return (w1*w2 - x1*x2 - y1*y2 - z1*z2,
            w1*x2 + x1*w2 + y1*z2 - z1*y2,
            w1*y2 - x1*z2 + y1*w2 + z1*x2,
            w1*z2 + x1*y2 - y1*x2 + z1*w2)

def qconj(a):  w, x, y, z = a; return (w, -x, -y, -z)
def qinv_unit(a): return qconj(a)            # unit quaternion inverse = conjugate
def qnsq(a):   return sum(c*c for c in a)
def qf(a):     return tuple(float(c) for c in a)

def qsub_norm(a, b):
    """‖a − b‖ in R^4 (used for exact-equality checks)."""
    return sum((p - q) ** 2 for p, q in zip(a, b))

def geo_angle(a, b):
    """Geodesic rotation angle (radians) between two unit quaternions."""
    d = abs(sum(p*q for p, q in zip(qf(a), qf(b))))
    return 2 * math.acos(max(-1.0, min(1.0, d)))

def qlog_halfaxis(a):
    """Return (half-angle θ, unit axis) for unit quaternion a=(w,v)."""
    w, x, y, z = qf(a)
    vn = math.sqrt(x*x + y*y + z*z)
    if vn < 1e-15:
        return 0.0, (1.0, 0.0, 0.0)
    th = math.atan2(vn, max(-1.0, min(1.0, w)))
    return th, (x/vn, y/vn, z/vn)

def qpow(a, t):
    """Fractional power a^t of a unit quaternion (float; principal branch)."""
    th, ax = qlog_halfaxis(a)
    s = math.sin(t*th)
    return (math.cos(t*th), s*ax[0], s*ax[1], s*ax[2])


# --------------------------------------------------------------------------- #
# SOURCE A — the Hyperbolic Sparsemax kernel walk (rational, parity-safe)
# --------------------------------------------------------------------------- #
def rquat_from_v(vx, vy, vz):
    """Inverse stereographic R^3 -> S^3: exact rational unit quaternion."""
    n = 1 + vx*vx + vy*vy + vz*vz
    return ((1 - (vx*vx + vy*vy + vz*vz)) / n, 2*vx/n, 2*vy/n, 2*vz/n)

# 64-voxel lattice (Z/4Z)^3, each voxel a rational annular spinor
LATTICE = {}
for a in range(4):
    for b in range(4):
        for c in range(4):
            LATTICE[(a, b, c)] = rquat_from_v(F(2*a-3, 8), F(2*b-3, 8), F(2*c-3, 8))

def hyper_quad(qi, qj):
    """Rational hyperbolic quadrance Q_h = ((1 - m_i·m_j)/(w_i w_j))^2 - 1 (float ok)."""
    wi, xi, yi, zi = qf(qi); wj, xj, yj, zj = qf(qj)
    dot_m = xi*xj + yi*yj + zi*zj
    val = (1 - dot_m) / (wi*wj)
    return val*val - 1.0

def sparsemax(logits):
    """Closed-form sparsemax (Thm 9): p_j = [z_j - tau]_+, tau s.t. sum=1."""
    z = sorted(logits, reverse=True)
    css = 0.0
    tau = 0.0; k = 0
    for i, zi in enumerate(z, 1):
        css += zi
        t = (css - 1.0) / i
        if zi - t > 0:
            tau = t; k = i
        else:
            break
    return [max(0.0, zj - tau) for zj in logits]

FACE = [(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(0,0,1),(0,0,-1)]

def sparsemax_walk(n, seed=7, lam=0.4):
    """Kernel-driven walk on the 64-voxel lattice. The sparsemax row over the
    self-loop + 6 face neighbours selects the next voxel; emitted spinors are the
    voxel's exact rational quaternion, so consecutive deltas are exact rational.
    A simple reverse-valve (forward = +axis) breaks detailed balance (directed)."""
    s = seed & 0xffffffff
    def rnd():
        nonlocal s
        s = (s + 0x6D2B79F5) & 0xffffffff
        t = (s ^ (s >> 15)) * (1 | s) & 0xffffffff
        t = (t + ((t ^ (t >> 7)) * (61 | t) & 0xffffffff)) ^ t & 0xffffffff
        return ((t ^ (t >> 14)) & 0xffffffff) / 4294967296
    cur = (1, 1, 1)
    states = [LATTICE[cur]]
    for _ in range(n - 1):
        qi = LATTICE[cur]
        cands = [cur] + [tuple((cur[d] + FACE[f][d]) % 4 for d in range(3)) for f in range(6)]
        # logit z = -beta*Q_h - lam*A ; beta = (1-w)/w ; A=0 forward(+), 1 reverse(-)
        wi = float(qi[0]); beta = (1 - wi) / wi
        logits = []
        for idx, cj in enumerate(cands):
            A = 0.0 if idx == 0 else (1.0 if FACE[idx-1][0] + FACE[idx-1][1] + FACE[idx-1][2] < 0 else 0.0)
            logits.append(-beta * hyper_quad(qi, LATTICE[cj]) - lam * A)
        row = sparsemax(logits)
        tot = sum(row) or 1.0
        r = rnd() * tot; acc = 0.0; pick = 0
        for idx, p in enumerate(row):
            acc += p
            if acc >= r:
                pick = idx; break
        cur = cands[pick]
        states.append(LATTICE[cur])
    return states


# --------------------------------------------------------------------------- #
# SOURCE B/C — generic float streams (codex accepts any unit-quaternion stream)
# --------------------------------------------------------------------------- #
def _axis_angle(ax, th):
    n = math.sqrt(sum(c*c for c in ax)) or 1.0
    s = math.sin(th/2)
    return (math.cos(th/2), s*ax[0]/n, s*ax[1]/n, s*ax[2]/n)

def smooth_stream(n):
    q = (1.0, 0.0, 0.0, 0.0); out = [q]
    for i in range(1, n):
        d = _axis_angle((1.0, 0.3*math.sin(0.01*i), 0.2), 0.02 + 0.01*math.sin(0.02*i))
        q = qmul(q, d)
        nf = math.sqrt(qnsq(q)); q = tuple(c/nf for c in q); out.append(q)
    return out

def random_stream(n, seed=11):
    s = seed
    def rnd():
        nonlocal s; s = (1103515245*s + 12345) & 0x7fffffff; return s/0x7fffffff
    q = (1.0, 0.0, 0.0, 0.0); out = [q]
    for _ in range(n-1):
        d = _axis_angle((rnd()-.5, rnd()-.5, rnd()-.5), 0.6*rnd())
        q = qmul(q, d); nf = math.sqrt(qnsq(q)); q = tuple(c/nf for c in q); out.append(q)
    return out


# --------------------------------------------------------------------------- #
# OFFLOAD — outer deltas -> middle (inverse, 10-blocks) -> inner (var 7/8 blocks)
# --------------------------------------------------------------------------- #
def outer_deltas(states):
    """720 deltas from 721 states: Δq_n = q_n q_{n-1}^{-1} (n = 1..720)."""
    return [qmul(states[n], qinv_unit(states[n-1])) for n in range(1, len(states))]

def middle_blocks(outer):
    """ΔQ_m = ordered product of 10 consecutive outer deltas (Prop 2). Later deltas
    multiply on the left, so the block telescopes to q_{10m+10} q_{10m}^{-1}."""
    blocks = []
    for m in range(N_MIDDLE):
        net = outer[GROUP_OM*m]
        for r in range(1, GROUP_OM):
            net = qmul(outer[GROUP_OM*m + r], net)
        blocks.append(net)
    return blocks                         # ΔQ(12->10)_m ; middle delta = inverse

def partition_72_10():
    """Prop 3: A_ℓ = { m : floor(10m/72) = ℓ }, sizes ∈ {7,8}, covering 0..71."""
    A = [[] for _ in range(N_INNER)]
    for m in range(N_MIDDLE):
        A[(10*m)//N_MIDDLE].append(m)
    return A

def inner_blocks(middle_deltas, A):
    inner = []
    for ell in range(N_INNER):
        net = (1.0, 0.0, 0.0, 0.0) if not isinstance(middle_deltas[0][0], F) else (F(1), F(0), F(0), F(0))
        for m in A[ell]:
            net = qmul(middle_deltas[m], net)
        inner.append(net)
    return inner


# --------------------------------------------------------------------------- #
# RECONSTRUCTION — store a scale, SLERP fractional-power lift back up
# --------------------------------------------------------------------------- #
def reconstruct_from_outer(states, outer):
    """Store all outer deltas -> exact telescoping reconstruction."""
    rec = [states[0]]
    for d in outer:
        rec.append(qmul(d, rec[-1]))
    return rec

def reconstruct_from_middle(states, blocks):
    """Store 72 middle nets -> distribute each over its 10 frames by ^(1/10) SLERP."""
    rec = [states[0]]
    for m in range(N_MIDDLE):
        step = qpow(qf(blocks[m]), 1.0/GROUP_OM)
        for _ in range(GROUP_OM):
            rec.append(qmul(step, rec[-1]))
    return rec[:len(states)]

def reconstruct_from_inner(states, inner, A, middle_blocks_):
    """Store 10 inner nets -> lift inner->middle (^1/|A|) then middle->outer (^1/10)."""
    rec = [states[0]]
    for ell in range(N_INNER):
        B = len(A[ell]) or 1
        mid_step = qpow(qf(inner[ell]), 1.0/B)                 # one coarse step per middle frame
        for _ in A[ell]:
            out_step = qpow(mid_step, 1.0/GROUP_OM)
            for _ in range(GROUP_OM):
                rec.append(qmul(out_step, rec[-1]))
    return rec[:len(states)]

def rmse_angle(a, b):
    errs = [geo_angle(a[i], b[i]) for i in range(min(len(a), len(b)))]
    return math.sqrt(sum(e*e for e in errs) / len(errs))


# --------------------------------------------------------------------------- #
# CENTROID + HOLONOMY
# --------------------------------------------------------------------------- #
def projected_mean(states):
    """The codex 'centroid': sign-aligned sum of unit quaternions, normalized
    (order-independent, unlike the incremental SLERP mean)."""
    acc = [0.0, 0.0, 0.0, 0.0]; ref = qf(states[0])
    for s in states:
        sf = qf(s)
        if sum(p*q for p, q in zip(sf, ref)) < 0: sf = tuple(-c for c in sf)
        for i in range(4): acc[i] += sf[i]
    nf = math.sqrt(sum(c*c for c in acc)) or 1.0
    return tuple(c/nf for c in acc)

def net_product(deltas):
    net = qf(deltas[0])
    for d in deltas[1:]: net = qmul(qf(d), net)
    return net

def holonomy_norm(Q8, Q10, Q12):
    H = qmul(qf(Q8), qmul(qf(Q10), qf(Q12)))
    th, _ = qlog_halfaxis(H)
    return 2*abs(th)                       # ‖rotation angle of H‖

def holonomy_close(Q8, Q10, Q12):
    """Paper §11: absorb the residual drift ε=log H into the inner equilibrium
    gear so the cycle product H = Q8·Q10·Q12 closes to the identity."""
    H = qmul(qf(Q8), qmul(qf(Q10), qf(Q12)))
    Q8c = qmul(qinv_unit(H), qf(Q8))       # redistribute holonomy into the coarse anchor
    Hc = qmul(Q8c, qmul(qf(Q10), qf(Q12)))
    th, _ = qlog_halfaxis(Hc)
    return 2*abs(th)


# --------------------------------------------------------------------------- #
# 64-bit Einstein-tile packet (FP8 E4M3 quaternion core) — lossy vs rational
# --------------------------------------------------------------------------- #
def f32_to_e4m3(f):
    if f == 0: return 0
    sign = 1 if f < 0 else 0; f = abs(f)
    e = math.floor(math.log2(f)) if f > 0 else 0
    e = max(-6, min(8, e)); mant = f / (2**e) - 1.0
    m = max(0, min(7, round(mant*8)))
    return (sign << 7) | (((e + 7) & 0xF) << 3) | m

def e4m3_to_f32(b):
    sign = -1 if (b >> 7) & 1 else 1
    e = ((b >> 3) & 0xF) - 7; m = b & 0x7
    if ((b >> 3) & 0xF) == 0: return sign * (2**(-6)) * (m/8)
    return sign * (2**e) * (1 + m/8)

def pack_e4m3(q):
    w, x, y, z = qf(q)
    return (f32_to_e4m3(w) << 24) | (f32_to_e4m3(x) << 16) | (f32_to_e4m3(y) << 8) | f32_to_e4m3(z)

def unpack_e4m3(p):
    q = (e4m3_to_f32((p >> 24) & 0xFF), e4m3_to_f32((p >> 16) & 0xFF),
         e4m3_to_f32((p >> 8) & 0xFF), e4m3_to_f32(p & 0xFF))
    n = math.sqrt(qnsq(q)) or 1.0
    return tuple(c/n for c in q)


# --------------------------------------------------------------------------- #
# Experiment
# --------------------------------------------------------------------------- #
def run_source(name):
    # 721 states -> exactly 720 deltas -> 72 clean blocks of 10
    if name == 'sparsemax':
        states = sparsemax_walk(N_OUTER + 1)
        rational = True
    elif name == 'smooth':
        states = smooth_stream(N_OUTER + 1); rational = False
    else:
        states = random_stream(N_OUTER + 1); rational = False

    outer = outer_deltas(states)
    blocks = middle_blocks(outer)                       # ΔQ(12->10)_m
    middle_deltas = [qinv_unit(b) for b in blocks]      # inverse middle gear
    A = partition_72_10()
    inner = inner_blocks(middle_deltas, A)

    # Prop 2: telescoping exact — block m product == q_{10m+10} q_{10m}^{-1}
    tele_err = 0.0; exact_zero = True
    for m in range(N_MIDDLE):
        lhs = blocks[m]
        rhs = qmul(states[GROUP_OM*m + GROUP_OM], qinv_unit(states[GROUP_OM*m]))
        e = qsub_norm(lhs, rhs)
        tele_err = max(tele_err, float(e))
        if e != 0: exact_zero = False

    return dict(name=name, rational=rational, states=states, outer=outer,
                blocks=blocks, middle_deltas=middle_deltas, A=A, inner=inner,
                tele_err=tele_err, exact_zero=exact_zero)


def main():
    print("Ontological Codex — three-scale quaternionic delta-chain memory offload")
    print("=" * 78)
    print(f"  gears: outer {N_OUTER} (12) -> middle {N_MIDDLE} (10, inverse) -> inner {N_INNER} (8)")

    # Prop 3 (combinatorial, source-independent)
    A = partition_72_10()
    sizes = sorted(set(len(a) for a in A))
    print(f"\n[Prop 3] 72->10 partition: block sizes {sizes}, total "
          f"{sum(len(a) for a in A)}  ({'PASS' if sizes==[7,8] and sum(len(a) for a in A)==72 else 'FAIL'})")

    # Prop 4 (SLERP fractional power exact on principal branch)
    qtest = _axis_angle((0.3, 0.5, 0.2), 1.1)
    back = qtest
    lifted = qpow(qtest, 1.0/10)
    acc = (1.0, 0.0, 0.0, 0.0)
    for _ in range(10): acc = qmul(lifted, acc)
    p4 = geo_angle(acc, qtest)
    print(f"[Prop 4] (q^(1/10))^10 == q : reconstruction angle {p4:.2e} rad  "
          f"({'PASS' if p4 < 1e-12 else 'FAIL'})")

    print("\n  Source comparison (telescoping exactness depends on rational vs float):")
    print(f"      {'source':<12}{'rational?':>10}{'Prop2 telescope err':>22}{'exact 0?':>10}")
    results = {}
    for name in ('sparsemax', 'smooth', 'random'):
        r = run_source(name); results[name] = r
        print(f"      {name:<12}{str(r['rational']):>10}{r['tele_err']:>22.2e}"
              f"{('YES' if r['exact_zero'] else 'no'):>10}")
    print("    → the Sparsemax source is exact-rational, so the offload telescopes to")
    print("      EXACTLY zero error (parity-safe); float sources telescope to machine")
    print("      precision (~1e-30 on 10-delta blocks).")

    # Rate-distortion + centroid + holonomy, per source
    for name in ('sparsemax', 'smooth', 'random'):
        r = results[name]; states = r['states']
        rec_o = reconstruct_from_outer(states, r['outer'])
        rec_m = reconstruct_from_middle(states, r['blocks'])
        rec_i = reconstruct_from_inner(states, r['inner'], r['A'], r['blocks'])
        e_o, e_m, e_i = (rmse_angle(states, x) for x in (rec_o, rec_m, rec_i))
        bits_o, bits_m, bits_i = N_OUTER*64, N_MIDDLE*64, N_INNER*64
        cent = projected_mean(states); netc = net_product(r['outer'])
        Q12 = qmul(states[-1], qinv_unit(states[0]))      # telescoped outer net
        Q10 = net_product(r['middle_deltas']); Q8 = net_product(r['inner'])
        hol = holonomy_norm(Q8, Q10, Q12)
        hol_c = holonomy_close(Q8, Q10, Q12)
        print(f"\n  [{name}] rate–distortion dial (store one scale, SLERP-reconstruct the 720):")
        print(f"      {'stored scale':<22}{'bits':>9}{'vs full':>9}{'recon RMSE (rad)':>18}")
        print(f"      {'outer 720 (exact)':<22}{bits_o:>9}{'1.0x':>9}{e_o:>18.2e}")
        print(f"      {'middle 72':<22}{bits_m:>9}{f'{bits_o/bits_m:.0f}x':>9}{e_m:>18.3f}")
        print(f"      {'inner 10':<22}{bits_i:>9}{f'{bits_o/bits_i:.0f}x':>9}{e_i:>18.3f}")
        print(f"      centroid: mean=({cent[0]:.3f},{cent[1]:.3f},{cent[2]:.3f},{cent[3]:.3f})")
        print(f"      holonomy ‖H‖={hol:.3f} rad  →  closes to {hol_c:.1e} after inner-anchor renorm")

    # 64-bit packet: E4M3 lossy vs rational-exact storage (sparsemax source)
    sp = results['sparsemax']
    e4 = 0.0
    for d in sp['outer'][:200]:
        e4 = max(e4, geo_angle(d, unpack_e4m3(pack_e4m3(d))))
    print(f"\n  [packet] 64-bit Einstein-tile core, sparsemax deltas:")
    print(f"      E4M3 (FP8) round-trip max angle err  : {math.degrees(e4):.2f}°  (lossy)")
    print(f"      exact-rational storage round-trip err : 0  (the deltas ARE rational ⇒ lossless,")
    print(f"        parity-safe; FP8 is only for a fixed-width wire format if size is capped)")

    print("\n" + "=" * 78)
    print("FINDINGS (measured):")
    print(" • The codex is the Memory-Offload Pump's three-gear chain made executable:")
    print("   720→72→10, inverse middle gear, blocks of 7/8 (Prop 3), SLERP lift (Prop 4).")
    print(" • Telescoping (Prop 2) is EXACT — and with the rational Sparsemax source it is")
    print("   exactly zero error (parity-safe), the paper's rational/transcendental split.")
    print(" • 'Compress to a centroid AND re-create the structure' is a rate–distortion DIAL:")
    print("   store outer = lossless (telescoping); store middle/inner = compact but the")
    print("   interior is SLERP-approximated (low RMSE for smooth/structured, high for noise).")
    print(" • Dual-Helix governance: the mean centroid is the untrusted SUMMARY (Helix A);")
    print("   the logged delta chain is the EXACT reconstruction (Helix B). No factoring")
    print("   speedup is claimed (TQF jump-ahead / harmonic bias measured ≤0 elsewhere).")


if __name__ == "__main__":
    main()