"""
FRAMERATE FROM FIRST PRINCIPLES — PERCOLATION THRESHOLD + CASCADE
================================================================================
Date: 2026-04-03
Author: Jonathan Shelton (theory), Claude (computation)
Status: EXPLORATORY CALCULATION

The hypothesis:
  1. Each dimension fills from a floor frequency to a ceiling (f_max)
  2. The floor of D+1 = the ceiling of D (cascade)
  3. The overflow is PULSED, not continuous
  4. The pulse fires at the percolation threshold (~29%)
  5. C_potential IS the bandwidth constraint (not separate from frequency)
  6. Framerate = Nyquist sampling rate for the frequency range

Question: does this produce the Fibonacci framerate ratios?

Target ratios (from Fibonacci budget):
  c₂/c₁ = 0.625/0.250 = 2.500
  c₃/c₂ = 1.000/0.625 = 1.600
  c₄/c₃ = 1.625/1.000 = 1.625
================================================================================
"""

import numpy as np

# Known constants
phi_c = 0.2895  # 3D continuum percolation threshold (Lorenz & Ziff 2001)
phi_c_sc = 0.3116  # Simple cubic site percolation
r_eq = 0.5  # Decoherence equilibrium ceiling

# Target Fibonacci framerates (normalized to c₃ = 1)
fib_targets = {
    1: 0.250,
    2: 0.625,
    3: 1.000,
    4: 1.625,
    5: 2.625,
    6: 4.250,
}

fib_ratios = {
    '2/1': 2.500,
    '3/2': 1.600,
    '4/3': 1.625,
    '5/4': 1.615,
    '6/5': 1.619,
}

print("=" * 70)
print("FRAMERATE DERIVATION — PERCOLATION + CASCADE")
print("=" * 70)
print(f"\nPercolation threshold (continuum): φ_c = {phi_c}")
print(f"Percolation threshold (SC site):   φ_c = {phi_c_sc}")
print(f"Decoherence ceiling:               r   = {r_eq}")
print(f"\nTarget Fibonacci ratios: {fib_ratios}")

# ============================================================================
# DENSITY OF STATES IN D DIMENSIONS
# ============================================================================
# Mode count from f_floor to f_max in D-dimensional space:
#   M_D = ∫_{f_floor}^{f_max} g_D(f) df
#
# Where g_D(f) ∝ f^(D-1) (standard result)
#
# So: M_D = [f_max^D - f_floor^D] / D
#
# This is the number of independent modes — the information capacity.
#
# FRAMERATE = NYQUIST = 2 × BANDWIDTH
# At r = 0.5: S/N = 1, Shannon capacity = bandwidth
# So: c_D ∝ 2 × M_D
# ============================================================================

def mode_count(D, f_floor, f_max):
    """Mode count in D dimensions from f_floor to f_max"""
    return (f_max**D - f_floor**D) / D


# ============================================================================
# MODEL 1: SELF-SIMILAR CASCADE
# ============================================================================
# Each dimension has the same expansion ratio α = f_max / f_floor
# The cascade compounds: f_floor(D+1) = f_max(D) = α × f_floor(D)
#
# The question: what α produces Fibonacci ratios?
# ============================================================================

print("\n" + "=" * 70)
print("MODEL 1: SELF-SIMILAR CASCADE (same α for all dimensions)")
print("=" * 70)

def model1_framerates(alpha, n_dims=6):
    """Compute framerate ratios for a self-similar cascade"""
    f_floor = np.zeros(n_dims + 1)
    f_max = np.zeros(n_dims + 1)
    modes = np.zeros(n_dims + 1)

    f_floor[1] = 1.0  # 1D floor = base frequency (normalized)

    for D in range(1, n_dims + 1):
        f_max[D] = alpha * f_floor[D]
        modes[D] = mode_count(D, f_floor[D], f_max[D])
        if D < n_dims:
            f_floor[D + 1] = f_max[D]

    # Framerates proportional to mode count
    c = 2 * modes

    # Normalize to c₃ = 1
    if c[3] > 0:
        c = c / c[3]

    return c, modes, f_floor, f_max

# Scan α values
print(f"\n{'α':>8} | {'c₁':>8} | {'c₂':>8} | {'c₃':>8} | {'c₄':>8} | {'c₂/c₁':>8} | {'c₃/c₂':>8} | {'c₄/c₃':>8}")
print("-" * 80)

best_alpha = None
best_error = float('inf')

for alpha in np.arange(1.5, 5.0, 0.05):
    c, _, _, _ = model1_framerates(alpha)
    if c[1] > 0 and c[2] > 0:
        r21 = c[2] / c[1]
        r32 = c[3] / c[2]
        r43 = c[4] / c[3]
        error = abs(r21 - 2.5) + abs(r32 - 1.6) + abs(r43 - 1.625)
        if error < best_error:
            best_error = error
            best_alpha = alpha
        if abs(alpha - round(alpha, 1)) < 0.001:
            print(f"{alpha:8.3f} | {c[1]:8.4f} | {c[2]:8.4f} | {c[3]:8.4f} | {c[4]:8.4f} | {r21:8.4f} | {r32:8.4f} | {r43:8.4f}")

print(f"\nBest α = {best_alpha:.4f} (total error = {best_error:.4f})")
c, modes, f_fl, f_mx = model1_framerates(best_alpha)
r21 = c[2] / c[1]
r32 = c[3] / c[2]
r43 = c[4] / c[3]
print(f"  c₁={c[1]:.4f}, c₂={c[2]:.4f}, c₃={c[3]:.4f}, c₄={c[4]:.4f}")
print(f"  c₂/c₁={r21:.4f} (target 2.500)")
print(f"  c₃/c₂={r32:.4f} (target 1.600)")
print(f"  c₄/c₃={r43:.4f} (target 1.625)")


# ============================================================================
# MODEL 2: PERCOLATION-MODIFIED CASCADE
# ============================================================================
# The overflow pulse carries φ_c fraction of the boundary energy.
# The floor of the next dimension is NOT f_max — it's the frequency
# corresponding to the percolation-capture energy.
#
# f_floor(D+1) = f_max(D) × φ_c^(1/D)  ← energy scales as f^D,
#   so the frequency corresponding to φ_c fraction of f_max energy
#   is f_max × φ_c^(1/D)
#
# Wait — if the overflow carries the FULL frequency but only φ_c
# of the energy, the floor frequency is still f_max(D). The φ_c
# affects the AMPLITUDE of the founding pulse, not its frequency.
#
# Let me try: the overflow fires at f_max(D), but the next dimension's
# effective starting energy = φ_c × E_D. This determines how quickly
# the next dimension can fill.
#
# ALTERNATIVE: φ_c constrains the expansion ratio α.
# At r = 0.5, the bubble fills φ_c of the boundary before it fires.
# The expansion ratio α = 1 / (1 - φ_c)?
# Or: α is set by how much energy the dimension can hold before
# the boundary percolates.
# ============================================================================

print("\n" + "=" * 70)
print("MODEL 2: PERCOLATION THRESHOLD SETS EXPANSION RATIO")
print("=" * 70)

# If the bubble before overflow represents φ_c of the TOTAL capacity,
# then the usable capacity is (1 - φ_c) and the bubble is φ_c.
# The expansion ratio could be: α = 1 / (1 - φ_c)

alpha_perc = 1.0 / (1.0 - phi_c)
print(f"\nα = 1/(1 - φ_c) = 1/(1 - {phi_c}) = {alpha_perc:.4f}")

c2, _, _, _ = model1_framerates(alpha_perc)
if c2[1] > 0 and c2[2] > 0:
    print(f"  c₂/c₁ = {c2[2]/c2[1]:.4f} (target 2.500)")
    print(f"  c₃/c₂ = {c2[3]/c2[2]:.4f} (target 1.600)")
    print(f"  c₄/c₃ = {c2[4]/c2[3]:.4f} (target 1.625)")

# What about α = 1 + φ_c?
alpha_perc2 = 1.0 + phi_c
print(f"\nα = 1 + φ_c = {alpha_perc2:.4f}")
c2b, _, _, _ = model1_framerates(alpha_perc2)
if c2b[1] > 0 and c2b[2] > 0:
    print(f"  c₂/c₁ = {c2b[2]/c2b[1]:.4f} (target 2.500)")
    print(f"  c₃/c₂ = {c2b[3]/c2b[2]:.4f} (target 1.600)")
    print(f"  c₄/c₃ = {c2b[4]/c2b[3]:.4f} (target 1.625)")


# ============================================================================
# MODEL 3: DIMENSION-DEPENDENT EXPANSION FROM MODE GEOMETRY
# ============================================================================
# The expansion ratio α_D is NOT the same for all dimensions.
# In D dimensions, the mode count scales as f^D, so the same
# "percolation fraction" of modes corresponds to different frequency
# ranges in different dimensions.
#
# If φ_c fraction of the total mode space at f_max must be in the
# "bubble" before overflow:
#   modes_bubble / modes_total = φ_c
#   [f_max^D - f_bubble^D] / [f_max^D - f_floor^D] = φ_c
#
# Where f_bubble is the frequency at which the bubble starts
# (the equilibrium ceiling). Setting f_bubble/f_max = β:
#   [1 - β^D] / [1 - (f_floor/f_max)^D] = φ_c
#
# The usable range is f_floor to f_bubble.
# The bubble range is f_bubble to f_max.
# The framerate is determined by the USABLE range, not the bubble.
# ============================================================================

print("\n" + "=" * 70)
print("MODEL 3: PERCOLATION PARTITIONS THE MODE SPACE")
print("=" * 70)
print("""
  The mode space from f_floor to f_max is partitioned:
    (1 - φ_c) of modes = USABLE (contribute to framerate)
    φ_c of modes = BUBBLE (accumulates before overflow)

  The framerate is determined by the usable modes only.
  The bubble's energy becomes the founding pulse of D+1.
""")

def model3_framerates(phi_c_val, alpha, n_dims=6):
    """
    Framerate from usable modes (excluding percolation bubble).

    Total mode space: f_floor to f_max (expansion α).
    Usable fraction: (1 - φ_c) of total modes.
    Bubble fraction: φ_c of total modes.

    The bubble sits at the TOP of the frequency range.
    f_bubble: where the bubble starts.
    [f_max^D - f_bubble^D] / [f_max^D - f_floor^D] = φ_c
    → f_bubble^D = f_max^D - φ_c × (f_max^D - f_floor^D)
    → f_bubble = [f_max^D - φ_c × (f_max^D - f_floor^D)]^(1/D)
    → f_bubble = [(1-φ_c) × f_max^D + φ_c × f_floor^D]^(1/D)
    """
    f_floor = np.zeros(n_dims + 1)
    f_max = np.zeros(n_dims + 1)
    f_bubble = np.zeros(n_dims + 1)
    modes_usable = np.zeros(n_dims + 1)
    modes_bubble = np.zeros(n_dims + 1)

    f_floor[1] = 1.0

    for D in range(1, n_dims + 1):
        f_max[D] = alpha * f_floor[D]

        # Find f_bubble: where bubble starts
        fm_D = f_max[D]**D
        ff_D = f_floor[D]**D
        fb_D = (1 - phi_c_val) * fm_D + phi_c_val * ff_D
        f_bubble[D] = fb_D ** (1.0/D)

        # Usable modes (floor to bubble)
        modes_usable[D] = mode_count(D, f_floor[D], f_bubble[D])

        # Bubble modes (bubble to max)
        modes_bubble[D] = mode_count(D, f_bubble[D], f_max[D])

        if D < n_dims:
            f_floor[D + 1] = f_max[D]

    # Framerate from USABLE modes
    c = 2 * modes_usable

    # Normalize
    if c[3] > 0:
        c = c / c[3]

    return c, modes_usable, modes_bubble, f_bubble

# Scan α for Model 3
print(f"Scanning α with φ_c = {phi_c}:")
print(f"\n{'α':>8} | {'c₂/c₁':>8} | {'c₃/c₂':>8} | {'c₄/c₃':>8} | {'error':>8}")
print("-" * 55)

best_alpha3 = None
best_error3 = float('inf')

for alpha in np.arange(1.2, 8.0, 0.01):
    c3, mu, mb, fb = model3_framerates(phi_c, alpha)
    if c3[1] > 0 and c3[2] > 0:
        r21 = c3[2] / c3[1]
        r32 = c3[3] / c3[2]
        r43 = c3[4] / c3[3]
        error = abs(r21 - 2.5) + abs(r32 - 1.6) + abs(r43 - 1.625)
        if error < best_error3:
            best_error3 = error
            best_alpha3 = alpha

        # Print at intervals
        if abs(alpha * 100 - round(alpha * 100)) < 0.5 and alpha % 0.5 < 0.01:
            print(f"{alpha:8.3f} | {r21:8.4f} | {r32:8.4f} | {r43:8.4f} | {error:8.4f}")

print(f"\nBest α = {best_alpha3:.4f} (total error = {best_error3:.6f})")
c3, mu3, mb3, fb3 = model3_framerates(phi_c, best_alpha3)
r21 = c3[2] / c3[1]
r32 = c3[3] / c3[2]
r43 = c3[4] / c3[3]
print(f"  c₂/c₁ = {r21:.6f} (target 2.500)")
print(f"  c₃/c₂ = {r32:.6f} (target 1.600)")
print(f"  c₄/c₃ = {r43:.6f} (target 1.625)")

# Check if best α has a clean relationship to φ_c
print(f"\n  α / φ_c = {best_alpha3/phi_c:.4f}")
print(f"  α × φ_c = {best_alpha3*phi_c:.4f}")
print(f"  α - 1 = {best_alpha3 - 1:.4f}")
print(f"  1/(1-φ_c) = {1/(1-phi_c):.4f}")
print(f"  φ_c/(1-φ_c) = {phi_c/(1-phi_c):.4f}")
print(f"  1/φ_c = {1/phi_c:.4f}")
print(f"  (1+φ_c)/(1-φ_c) = {(1+phi_c)/(1-phi_c):.4f}")
print(f"  π × φ_c = {np.pi * phi_c:.4f}")
print(f"  e × φ_c = {np.e * phi_c:.4f}")
print(f"  φ (golden ratio) = {(1+np.sqrt(5))/2:.4f}")
print(f"  α / φ = {best_alpha3 / ((1+np.sqrt(5))/2):.4f}")


# ============================================================================
# MODEL 4: CAN WE FIND φ_c SUCH THAT α HAS A CLEAN RELATIONSHIP?
# ============================================================================
print("\n" + "=" * 70)
print("MODEL 4: SCANNING φ_c × α RELATIONSHIPS")
print("=" * 70)
print("Looking for φ_c values where best α has a clean form...\n")

phi_golden = (1 + np.sqrt(5)) / 2

for pc in np.arange(0.25, 0.35, 0.005):
    best_a = None
    best_e = float('inf')
    for alpha in np.arange(1.2, 8.0, 0.001):
        c4, _, _, _ = model3_framerates(pc, alpha)
        if c4[1] > 0 and c4[2] > 0:
            r21 = c4[2] / c4[1]
            r32 = c4[3] / c4[2]
            r43 = c4[4] / c4[3]
            error = abs(r21 - 2.5) + abs(r32 - 1.6) + abs(r43 - 1.625)
            if error < best_e:
                best_e = error
                best_a = alpha

    if best_a:
        print(f"  φ_c={pc:.3f}: best α={best_a:.4f}, "
              f"α/φ={best_a/phi_golden:.4f}, "
              f"α-1={best_a-1:.4f}, "
              f"1/(1-φ_c)={1/(1-pc):.4f}, "
              f"err={best_e:.6f}")


# ============================================================================
# MODEL 5: WHAT IF α ITSELF IS phi (golden ratio)?
# ============================================================================
print("\n" + "=" * 70)
print("MODEL 5: α = φ (golden ratio) — WHAT φ_c PRODUCES FIBONACCI?")
print("=" * 70)

alpha_phi = (1 + np.sqrt(5)) / 2
print(f"\nUsing α = φ = {alpha_phi:.6f}")
print(f"Scanning φ_c...\n")

for pc in np.arange(0.00, 0.50, 0.01):
    c5, _, _, _ = model3_framerates(pc, alpha_phi)
    if c5[1] > 0 and c5[2] > 0:
        r21 = c5[2] / c5[1]
        r32 = c5[3] / c5[2]
        r43 = c5[4] / c5[3]
        error = abs(r21 - 2.5) + abs(r32 - 1.6) + abs(r43 - 1.625)
        if error < 0.1:
            print(f"  φ_c={pc:.3f}: c₂/c₁={r21:.4f}, c₃/c₂={r32:.4f}, "
                  f"c₄/c₃={r43:.4f}, err={error:.6f}")


# ============================================================================
# MODEL 6: PURE INFORMATION — NO GEOMETRIC CONSTANTS
# ============================================================================
# Strip out π and L. Use pure mode count ratios.
# In D dimensions from f_floor to f_max:
#   M_D = (f_max^D - f_floor^D) / D
# With f_max = α × f_floor and cascade:
#   M_D = f_floor^D × (α^D - 1) / D
# And f_floor(D) = α^(D-1) × f₀ (from cascade)
# So: M_D = α^(D(D-1)) × f₀^D × (α^D - 1) / D
#
# Framerate ∝ M_D. Ratio:
#   c_{D+1}/c_D = [α^((D+1)D) × f₀^(D+1) × (α^(D+1)-1)/(D+1)]
#                / [α^(D(D-1)) × f₀^D × (α^D - 1) / D]
#               = α^(D(D+1) - D(D-1)) × f₀ × D/(D+1) × (α^(D+1)-1)/(α^D-1)
#               = α^(2D) × f₀ × D/(D+1) × (α^(D+1)-1)/(α^D-1)
#
# This still has f₀. But if we normalize to c₃, f₀ cancels.
# ============================================================================

print("\n" + "=" * 70)
print("MODEL 6: PURE MODE COUNT CASCADE (no geometric constants)")
print("=" * 70)

def model6_ratios(alpha, phi_c_val=0.0, n_dims=6):
    """
    Pure mode-count cascade with optional percolation partition.
    No geometric prefactors (π, L). Just D-dimensional mode scaling.
    """
    f_floor = np.zeros(n_dims + 1)
    f_max = np.zeros(n_dims + 1)
    modes = np.zeros(n_dims + 1)

    f_floor[1] = 1.0

    for D in range(1, n_dims + 1):
        f_max[D] = alpha * f_floor[D]

        if phi_c_val > 0:
            # Usable modes only (exclude bubble)
            fm_D = f_max[D]**D
            ff_D = f_floor[D]**D
            fb_D = (1 - phi_c_val) * fm_D + phi_c_val * ff_D
            f_bubble = fb_D ** (1.0 / D)
            modes[D] = mode_count(D, f_floor[D], f_bubble)
        else:
            modes[D] = mode_count(D, f_floor[D], f_max[D])

        if D < n_dims:
            f_floor[D + 1] = f_max[D]

    c = 2 * modes
    if c[3] > 0:
        c = c / c[3]
    return c

# Fine scan around the percolation threshold
print(f"\nFine scan: α and φ_c that give Fibonacci ratios\n")
print(f"{'α':>8} | {'φ_c':>6} | {'c₁':>8} | {'c₂':>8} | {'c₃':>8} | {'c₂/c₁':>8} | {'c₃/c₂':>8} | {'c₄/c₃':>8} | {'err':>8}")
print("-" * 90)

overall_best = (None, None, float('inf'))

for pc in np.arange(0.0, 0.50, 0.005):
    for alpha in np.arange(1.5, 6.0, 0.001):
        c6 = model6_ratios(alpha, pc)
        if c6[1] > 0 and c6[2] > 0:
            r21 = c6[2] / c6[1]
            r32 = c6[3] / c6[2]
            r43 = c6[4] / c6[3]
            error = abs(r21 - 2.5) + abs(r32 - 1.6) + abs(r43 - 1.625)
            if error < overall_best[2]:
                overall_best = (alpha, pc, error)

a_best, pc_best, err_best = overall_best
print(f"\nOVERALL BEST: α = {a_best:.4f}, φ_c = {pc_best:.4f}, error = {err_best:.6f}")
c_best = model6_ratios(a_best, pc_best)
print(f"  c₁={c_best[1]:.6f}, c₂={c_best[2]:.6f}, c₃={c_best[3]:.6f}, c₄={c_best[4]:.6f}")
print(f"  c₂/c₁ = {c_best[2]/c_best[1]:.6f} (target 2.500)")
print(f"  c₃/c₂ = {c_best[3]/c_best[2]:.6f} (target 1.600)")
print(f"  c₄/c₃ = {c_best[4]/c_best[3]:.6f} (target 1.625)")

# Check relationships
print(f"\n  Notable relationships at best fit:")
print(f"  α = {a_best:.6f}")
print(f"  φ_c = {pc_best:.6f}")
print(f"  α × φ_c = {a_best * pc_best:.6f}")
print(f"  α - 1 = {a_best - 1:.6f}")
print(f"  α / φ_golden = {a_best / phi_golden:.6f}")
print(f"  (α - 1) / φ_c = {(a_best - 1) / pc_best:.6f}" if pc_best > 0 else "")
print(f"  α² = {a_best**2:.6f}")
print(f"  √α = {np.sqrt(a_best):.6f}")
print(f"  α / π = {a_best / np.pi:.6f}")
print(f"  α / e = {a_best / np.e:.6f}")


# ============================================================================
# IMPORTANT CHECK: Without percolation (φ_c = 0), can ANY α match?
# ============================================================================
print("\n" + "=" * 70)
print("CONTROL: WITHOUT PERCOLATION (φ_c = 0)")
print("=" * 70)

best_a0 = None
best_e0 = float('inf')
for alpha in np.arange(1.5, 8.0, 0.001):
    c0 = model6_ratios(alpha, 0.0)
    if c0[1] > 0 and c0[2] > 0:
        r21 = c0[2] / c0[1]
        r32 = c0[3] / c0[2]
        r43 = c0[4] / c0[3]
        error = abs(r21 - 2.5) + abs(r32 - 1.6) + abs(r43 - 1.625)
        if error < best_e0:
            best_e0 = error
            best_a0 = alpha

print(f"\nBest without percolation: α = {best_a0:.4f}, error = {best_e0:.6f}")
c0 = model6_ratios(best_a0, 0.0)
print(f"  c₂/c₁ = {c0[2]/c0[1]:.6f} (target 2.500)")
print(f"  c₃/c₂ = {c0[3]/c0[2]:.6f} (target 1.600)")
print(f"  c₄/c₃ = {c0[4]/c0[3]:.6f} (target 1.625)")

print(f"\nBest WITH percolation:    error = {err_best:.6f}")
print(f"Best WITHOUT percolation: error = {best_e0:.6f}")
if err_best < best_e0:
    print(f"  → Percolation IMPROVES the fit by {best_e0 - err_best:.6f}")
else:
    print(f"  → Percolation does NOT improve the fit")


# ============================================================================
# ANALYTIC CHECK: What do the ratios HAVE to look like?
# ============================================================================
print("\n" + "=" * 70)
print("ANALYTIC: RATIO STRUCTURE OF THE CASCADE")
print("=" * 70)
print("""
For a cascade with expansion α and no percolation:
  M_D = (α^D - 1) / D  (modes, with f_floor normalized to 1)

  Since f_floor(D) = α^(D-1) (cascade), and modes ∝ f_floor^D:

  c_D ∝ α^(D(D-1)) × (α^D - 1) / D

Ratio c_{D+1}/c_D = α^(2D) × D/(D+1) × (α^(D+1)-1)/(α^D-1)
""")

print("With α = best value, checking analytic formula:")
alpha = best_a0

for D in range(1, 5):
    ratio = (alpha**(2*D)) * (D/(D+1)) * (alpha**(D+1) - 1) / (alpha**D - 1)
    print(f"  c_{D+1}/c_{D} = {ratio:.6f}")

print("\nThese ratios CANNOT all equal phi simultaneously because")
print("the α^(2D) term grows with D while D/(D+1) shrinks.")
print("The ratio is inherently DIMENSION-DEPENDENT.")
print("\nFibonacci ratios (2.5, 1.6, 1.625, 1.615...) are also")
print("dimension-dependent — they converge to φ from alternating sides.")
print("The question: does the mode-count cascade's D-dependence")
print("match the Fibonacci sequence's D-dependence?")