"""
CONVERGENCE TEST: QUADRATIC ENERGY FORMULA vs FIBONACCI COMPRESSION
================================================================================
Date: 2026-04-03
Author: Jonathan Shelton (theory), Claude (computation)
Status: CRITICAL VALIDATION

Two independent derivation paths:

PATH 1 (Paper 1, empirical):
  log₁₀(E/eV) = 0.1964 d² + 8.0932 d - 20.0373
  Fit from 3 measured boundary energies:
    d=2: 0.86 meV   (He superfluid transition)
    d=3: 1.022 MeV  (pair production threshold)
    d=4: 3.0 PeV    (cosmic ray proton knee)

PATH 2 (this session, derived):
  c_D = Budget(D) × v_write
  Budget(D) = Fibonacci sum at dimension D
  v_write = c/8 (constant frame-writing speed)
  Cascade: f_floor(D+1) = f_max(D)
  Overflow at percolation threshold φ_c ≈ 0.29

QUESTION: Do they produce the same numbers?
================================================================================
"""

import numpy as np

phi = (1 + np.sqrt(5)) / 2
c_light = 2.998e8  # m/s
h_planck = 6.626e-34  # J·s
eV = 1.602e-19  # J per eV

print("=" * 70)
print("PATH 1: QUADRATIC EQUATION (empirical, Paper 1)")
print("=" * 70)

def E_quadratic(d):
    """Paper 1 energy boundary formula"""
    log_E = 0.1964 * d**2 + 8.0932 * d - 20.0373
    return 10**log_E  # in eV

print("\nDimensional boundary energies:")
for d in range(1, 8):
    E = E_quadratic(d)
    log_E = np.log10(E)
    units = "eV"
    E_display = E
    if E < 1e-3:
        E_display = E * 1e6
        units = "μeV"
    elif E < 1:
        E_display = E * 1e3
        units = "meV"
    elif E > 1e15:
        E_display = E / 1e15
        units = "PeV"
    elif E > 1e12:
        E_display = E / 1e12
        units = "TeV"
    elif E > 1e9:
        E_display = E / 1e9
        units = "GeV"
    elif E > 1e6:
        E_display = E / 1e6
        units = "MeV"
    elif E > 1e3:
        E_display = E / 1e3
        units = "keV"

    status = ""
    if d == 2: status = " ← He superfluid (measured)"
    elif d == 3: status = " ← pair production (measured)"
    elif d == 4: status = " ← cosmic ray knee (measured)"
    else: status = " ← PREDICTION"

    print(f"  d={d}: log₁₀(E) = {log_E:8.4f}  →  {E_display:.3f} {units}{status}")


print("\n" + "=" * 70)
print("PATH 2: FIBONACCI COMPRESSION (derived this session)")
print("=" * 70)

# Fibonacci budgets
fib = [0, 2, 5, 8, 13, 21, 34, 55]  # Budget at each dimension
# These are F(D+1) + F(D+2) where F is standard Fibonacci

v_write_fraction = 1.0 / 8.0  # v_write = c/8

print("\nFramerates from Fibonacci budget:")
for d in range(1, 8):
    c_d = fib[d] * v_write_fraction
    print(f"  d={d}: Budget={fib[d]:3d}, c_d = {fib[d]}/8 × c = {c_d:.4f}c")


# ============================================================================
# THE CONNECTION: How does framerate relate to boundary energy?
# ============================================================================
print("\n" + "=" * 70)
print("THE CONNECTION: FRAMERATE → BOUNDARY ENERGY")
print("=" * 70)

print("""
The boundary energy E(D) is the total energy a dimension holds
at overflow. This is NOT the framerate — it's the CAPACITY.

Framerate = throughput (frames per second)
Energy = total accumulated content at overflow

The relationship:
  E(D) = (total modes at equilibrium) × (average energy per mode)

Total modes M(D) = integral of density of states from f_floor to f_max.
Average energy per mode at equilibrium: ε = h × f_mean (quantum)
  or ε = kT (thermal), where T is the temperature at r = 0.5.

The CASCADE means:
  f_floor(D) = f_max(D-1)  (each ceiling becomes next floor)
  f_max(D) > f_floor(D)    (new modes fill above the floor)

The COMPRESSION means:
  c_D = Budget(D) × v_write  (framerate = compressed throughput)
  M_unique(D) = Budget(D) / 2  (from Nyquist: c = 2B)
  M_raw(D) = M_unique(D) × S(D)  (raw modes = unique × symmetry)

The ENERGY is determined by raw modes × frequency:
  E(D) ∝ M_raw(D) × h × f_mean(D)
""")

# ============================================================================
# APPROACH 1: Direct frequency cascade
# ============================================================================
print("=" * 70)
print("APPROACH 1: Energy from frequency cascade")
print("=" * 70)

print("""
If we know the floor frequency at each dimension, the energy is:
  E(D) ~ h × f_floor(D)  (the characteristic energy scale)

The floor frequency cascades:
  f_floor(1) = f_base ≈ 8 Hz
  f_floor(D+1) = f_max(D)

The expansion at each dimension:
  f_max(D) = f_floor(D) × expansion(D)

The expansion(D) is determined by how many modes fit
between floor and ceiling. This depends on:
  - Dimensional mode density (∝ f^(D-1))
  - Compression at that level (Budget(D))
  - Information capacity per frame (constant)

Key insight: if info per frame = Budget(D) / 2 (from Nyquist),
and M_raw = info × S(D), and M_raw fills the frequency space,
then the expansion encodes the mode count.
""")

# Work backwards from measured data to extract the cascade
print("Working backwards from measured boundary energies:\n")

E_boundaries = {
    2: 0.86e-3,    # eV (He superfluid)
    3: 1.022e6,    # eV (pair production)
    4: 3.0e15,     # eV (cosmic ray knee)
}

# Convert to frequency: E = hf → f = E/h
for d in [2, 3, 4]:
    E_eV = E_boundaries[d]
    E_J = E_eV * eV
    f = E_J / h_planck
    print(f"  d={d}: E = {E_eV:.3e} eV → f = {f:.3e} Hz")

# The frequency RATIOS between boundaries
f_boundaries = {}
for d in [2, 3, 4]:
    f_boundaries[d] = E_boundaries[d] * eV / h_planck

print(f"\n  f(3)/f(2) = {f_boundaries[3]/f_boundaries[2]:.3e}")
print(f"  f(4)/f(3) = {f_boundaries[4]/f_boundaries[3]:.3e}")
print(f"  f(4)/f(2) = {f_boundaries[4]/f_boundaries[2]:.3e}")

# Log ratios
print(f"\n  log₁₀(f(3)/f(2)) = {np.log10(f_boundaries[3]/f_boundaries[2]):.4f}")
print(f"  log₁₀(f(4)/f(3)) = {np.log10(f_boundaries[4]/f_boundaries[3]):.4f}")

# The quadratic equation predicts these log ratios as:
# log₁₀(E(d+1)/E(d)) = 0.1964(2d+1) + 8.0932
# For d=2→3: 0.1964×5 + 8.0932 = 0.982 + 8.0932 = 9.0752
# For d=3→4: 0.1964×7 + 8.0932 = 1.3748 + 8.0932 = 9.4680
print(f"\n  Quadratic prediction for log₁₀(E(3)/E(2)): {0.1964*5 + 8.0932:.4f}")
print(f"  Quadratic prediction for log₁₀(E(4)/E(3)): {0.1964*7 + 8.0932:.4f}")
print(f"  Actual log₁₀(E(3)/E(2)): {np.log10(E_boundaries[3]/E_boundaries[2]):.4f}")
print(f"  Actual log₁₀(E(4)/E(3)): {np.log10(E_boundaries[4]/E_boundaries[3]):.4f}")


# ============================================================================
# APPROACH 2: Can the Fibonacci budget reproduce the quadratic?
# ============================================================================
print("\n" + "=" * 70)
print("APPROACH 2: Does log(Budget) follow a quadratic in D?")
print("=" * 70)

print("\nFibonacci budgets and their logs:")
for d in range(1, 8):
    log_b = np.log10(fib[d])
    print(f"  d={d}: Budget={fib[d]:5d}, log₁₀(Budget) = {log_b:.6f}")

# Fit quadratic to log(Budget) vs d
d_vals = np.array([1, 2, 3, 4, 5, 6, 7])
log_budgets = np.array([np.log10(fib[d]) for d in d_vals])

# Quadratic fit: log₁₀(Budget) = a×d² + b×d + c
coeffs = np.polyfit(d_vals, log_budgets, 2)
print(f"\nQuadratic fit to log₁₀(Budget(D)):")
print(f"  log₁₀(Budget) = {coeffs[0]:.6f} d² + {coeffs[1]:.6f} d + {coeffs[2]:.6f}")
print(f"\nCompare to Paper 1:")
print(f"  Paper 1:  0.1964 d² + 8.0932 d + const")
print(f"  Budget:   {coeffs[0]:.6f} d² + {coeffs[1]:.6f} d + {coeffs[2]:.6f}")
print(f"\n  Quadratic coeff ratio: {coeffs[0]/0.1964:.4f}")
print(f"  Linear coeff: Budget has {coeffs[1]:.4f} vs Paper 1 has 8.0932")
print(f"  (Linear coeff encodes absolute energy scale, not ratio)")

# The budget alone is too small to match. The quadratic term in the
# energy formula comes from the CASCADE, not just the budget.
# Let's compute the CUMULATIVE cascade.

print("\n" + "=" * 70)
print("APPROACH 3: CUMULATIVE CASCADE ENERGY")
print("=" * 70)

print("""
The boundary energy isn't just Budget(D). It's the CASCADE:
each dimension's overflow frequency becomes the next floor.

The total energy at overflow(D) involves:
  1. The floor frequency (from cascade): compounds multiplicatively
  2. The mode count at that level: determined by expansion + dimension
  3. The compression efficiency: Budget(D)

If the expansion at each dimension scales with the budget:
  f_max(D) / f_floor(D) = g(Budget(D))

Then the floor cascades as:
  f_floor(D+1) = f_floor(D) × g(Budget(D))

And log(f_floor) accumulates the sum of log(g(Budget)):
  log(f_floor(D)) = Σ_{k=1}^{D-1} log(g(Budget(k)))

If g(Budget) ∝ Budget, then:
  log(f_floor(D)) = Σ log(Budget(k)) for k=1..D-1
""")

# Compute cumulative log(Budget) as proxy for log(f_floor)
cum_log_budget = np.zeros(8)
for d in range(2, 8):
    cum_log_budget[d] = cum_log_budget[d-1] + np.log10(fib[d-1])

print("Cumulative log₁₀(Budget) as cascade proxy:")
for d in range(1, 8):
    print(f"  d={d}: cum_log = {cum_log_budget[d]:.6f}")

# The total log-energy should be:
# log(E(D)) ∝ log(f_floor(D)) + log(modes(D))
#            ∝ cum_log_budget(D) + D × log(Budget(D))
#
# But this is getting complicated. Let's try the direct product.

# ENERGY MODEL: E(D) ∝ f_floor(D) × Budget(D)
# where f_floor cascades through the Fibonacci budgets

print("\n  Energy model: E(D) ∝ f_floor(D) × Budget(D)")
print(f"  Starting from f_floor(1) = 8 Hz (the base pulse)\n")

f_floor_cascade = np.zeros(8)
f_max_cascade = np.zeros(8)
E_model = np.zeros(8)

f_floor_cascade[1] = 8.0  # Hz, the base pulse

for d in range(1, 8):
    # Each dimension expands by its budget
    # (budget = compressed modes, expansion = how far frequency stretches)
    f_max_cascade[d] = f_floor_cascade[d] * fib[d]
    E_model[d] = h_planck * f_max_cascade[d] / eV  # energy in eV

    if d < 7:
        f_floor_cascade[d+1] = f_max_cascade[d]

    print(f"  d={d}: f_floor={f_floor_cascade[d]:.3e} Hz, "
          f"×Budget({fib[d]})={f_max_cascade[d]:.3e} Hz, "
          f"E={E_model[d]:.3e} eV")

print(f"\nComparison: Model vs Quadratic vs Measured")
print(f"{'d':>3} | {'Model E (eV)':>14} | {'Quad E (eV)':>14} | {'Measured E':>14} | {'log ratio':>10}")
print("-" * 72)
for d in range(1, 7):
    E_q = E_quadratic(d)
    E_m = E_model[d]
    ratio = np.log10(E_m / E_q) if E_q > 0 and E_m > 0 else 0
    measured = ""
    if d == 2: measured = f"{0.86e-3:.3e}"
    elif d == 3: measured = f"{1.022e6:.3e}"
    elif d == 4: measured = f"{3.0e15:.3e}"
    print(f"{d:3d} | {E_m:14.3e} | {E_q:14.3e} | {measured:>14s} | {ratio:10.4f}")


# ============================================================================
# APPROACH 4: What if expansion = Budget^D (not just Budget)?
# ============================================================================
print("\n" + "=" * 70)
print("APPROACH 4: EXPANSION = Budget^(something)")
print("=" * 70)

# Try different power laws for expansion
for power_name, get_expansion in [
    ("Budget",     lambda d: fib[d]),
    ("Budget²",    lambda d: fib[d]**2),
    ("Budget^D",   lambda d: fib[d]**d),
    ("Budget^(D-1)", lambda d: fib[d]**(d-1)),
    ("2^Budget",   lambda d: 2**fib[d]),
]:
    ff = np.zeros(8)
    fm = np.zeros(8)
    Em = np.zeros(8)
    ff[1] = 8.0

    for d in range(1, 7):
        expansion = get_expansion(d)
        fm[d] = ff[d] * expansion
        Em[d] = h_planck * fm[d] / eV
        if d < 7: ff[d+1] = fm[d]

    # Compute log-energy and fit quadratic
    d_arr = np.array([2, 3, 4])
    logE = np.array([np.log10(Em[d]) if Em[d] > 0 else 0 for d in d_arr])

    # Compare to measured at d=2,3,4
    logE_measured = np.array([np.log10(0.86e-3), np.log10(1.022e6), np.log10(3.0e15)])
    error = np.sum(np.abs(logE - logE_measured))

    # Also fit quadratic to the model
    if all(Em[d] > 0 for d in range(1, 7)):
        d_fit = np.array(range(1, 7))
        logE_fit = np.array([np.log10(Em[d]) for d in d_fit])
        qcoeffs = np.polyfit(d_fit, logE_fit, 2)
        print(f"\n  {power_name:15s}: quad coeff = {qcoeffs[0]:.4f} "
              f"(target 0.1964), lin = {qcoeffs[1]:.4f} "
              f"(target 8.0932), error = {error:.2f}")
        print(f"    d=2: {Em[2]:.3e} eV (target 8.6e-4)")
        print(f"    d=3: {Em[3]:.3e} eV (target 1.022e+6)")
        print(f"    d=4: {Em[4]:.3e} eV (target 3.0e+15)")


# ============================================================================
# APPROACH 5: Self-consistent cascade with v_write
# ============================================================================
print("\n" + "=" * 70)
print("APPROACH 5: SELF-CONSISTENT CASCADE FROM v_write = c/8")
print("=" * 70)

print("""
The fundamental quantities:
  v_write = c/8 (frame writing speed, constant)
  Budget(D) = Fibonacci sum (compressed modes per frame)
  c_D = Budget(D) × v_write (framerate)

Energy at dimensional boundary:
  E(D) = h × f_overflow(D)

Where f_overflow is the highest frequency at which the dimension
overflows. This is related to the framerate and the number of
frames that have been written:

  f_overflow(D) = c_D × k_D / λ_D

But more directly: the boundary energy IS the energy equivalent
of the framerate's bandwidth.

  E(D) = h × c_D / λ_min(D)

Where λ_min(D) is the shortest wavelength at dimension D.
λ_min = c_D / f_max(D).

So: E(D) = h × f_max(D).

And f_max cascades: f_max(D) = f_floor(D) × expansion(D).
""")

# The key question: what is the expansion at each dimension?
# From the constant-info-per-frame model:
# M_unique = Budget(D) / 2 (Nyquist) = constant
# M_raw = M_unique × S(D) = Budget(D) × S(D) / 2
# M_raw = (f_max^D - f_floor^D) / D (from density of states)
# So: (expansion^D - 1) × f_floor^D / D = Budget(D) × S(D) / 2

S = {1: 2, 2: 12, 3: 48, 4: 1152, 5: 55296, 6: 15925248}

print("\nSolving for expansion at each dimension:")
print("(expansion^D - 1) × f_floor^D / D = Budget(D) × S(D) / 2\n")

ff = [0, 8.0, 0, 0, 0, 0, 0]  # Hz
expansion = [0, 0, 0, 0, 0, 0, 0]
fm = [0, 0, 0, 0, 0, 0, 0]
E_self = [0, 0, 0, 0, 0, 0, 0]

for d in range(1, 7):
    if d > 1:
        ff[d] = fm[d-1]

    # Solve: (α^D - 1) × ff^D / D = Budget × S / 2
    # α^D = 1 + D × Budget × S / (2 × ff^D)
    rhs = d * fib[d] * S[d] / (2 * ff[d]**d) if ff[d] > 0 else 0
    alpha_d = (1 + rhs) ** (1.0/d) if rhs > 0 else 1

    expansion[d] = alpha_d
    fm[d] = alpha_d * ff[d]
    E_self[d] = h_planck * fm[d] / eV

    print(f"  d={d}: f_floor={ff[d]:.6e} Hz, α={alpha_d:.6f}, "
          f"f_max={fm[d]:.6e} Hz, E={E_self[d]:.3e} eV")

print(f"\n{'d':>3} | {'Self-consistent':>16} | {'Quadratic':>16} | {'Measured':>16} | {'log₁₀ diff':>12}")
print("-" * 75)
for d in range(1, 7):
    E_q = E_quadratic(d)
    measured = ""
    if d == 2: measured = "8.6e-4"
    elif d == 3: measured = "1.022e+6"
    elif d == 4: measured = "3.0e+15"
    diff = np.log10(E_self[d] / E_q) if E_self[d] > 0 and E_q > 0 else 0
    print(f"{d:3d} | {E_self[d]:16.3e} | {E_q:16.3e} | {measured:>16s} | {diff:12.4f}")

# Fit quadratic to the self-consistent model
d_arr = np.array(range(1, 7))
logE_self = np.array([np.log10(E_self[d]) if E_self[d] > 0 else 0 for d in d_arr])
qc = np.polyfit(d_arr, logE_self, 2)
print(f"\nQuadratic fit to self-consistent model:")
print(f"  log₁₀(E) = {qc[0]:.4f} d² + {qc[1]:.4f} d + {qc[2]:.4f}")
print(f"\nPaper 1 quadratic:")
print(f"  log₁₀(E) = 0.1964 d² + 8.0932 d + (-20.0373)")
print(f"\nCoefficient comparison:")
print(f"  Quadratic: model {qc[0]:.4f} vs paper {0.1964:.4f} (ratio: {qc[0]/0.1964:.4f})")
print(f"  Linear:    model {qc[1]:.4f} vs paper {8.0932:.4f} (ratio: {qc[1]/8.0932:.4f})")
print(f"  Constant:  model {qc[2]:.4f} vs paper {-20.0373:.4f}")