#!/usr/bin/env python3
"""
ACCELERATION RAMP COMPUTATION — FRAMERATE GRADIENT AT EACH ELEMENT
================================================================================
Date: 2026-04-04
Author: Jonathan Shelton (theory), Claude (computation)
Status: OBSERVATIONAL — compute and observe

For each element, compute:
  1. The slope of the energy quadratic at d_eff (global framerate pressure)
  2. The local ramp position within the period's sub-range
  3. The fraction past the 29% bulge (how much of the orbit gets caught)
  4. The estimated "sweep fraction" — what portion of electron orbits
     pass through the acceleration zone

Then observe: do the elements needing corrections cluster at specific
sweep fractions? Does the sweep fraction correlate with known SO values?
================================================================================
"""

import numpy as np
import json

# Constants
amu_to_kg = 1.6605e-27
c_light = 2.998e8
h_planck = 6.626e-34
phi_c = 0.2895

def compton_freq(mass_amu):
    m_kg = mass_amu * amu_to_kg
    return m_kg * c_light**2 / h_planck

# Load sweep data
with open('/root/rhetroiluso/project_prometheus/time_ledger_theory/alchemical_geometry/results/unaudited/BLIND_SWEEP_PREDICTIONS_2026-04-03.json') as f:
    elements = json.load(f)

# Known SO values from NIST (the "answer key" — for COMPARISON only, not input)
NIST_SO = {
    21: 20, 22: 30, 23: 35, 24: 45, 25: 50, 26: 52, 27: 65, 28: 85,
    29: 83, 30: 90,
    39: 60, 40: 75, 41: 90, 42: 120, 43: 210, 44: 240, 45: 275,
    46: 200, 47: 180, 48: 250,
    72: 600, 73: 700, 74: 800, 75: 860, 76: 940, 77: 1020,
    78: 800, 79: 700, 80: 1300,
}

# Elements needing correction
need_correction = {43, 44, 45, 62, 75, 76, 77, 80, 84}

# =============================================================================
# COMPUTE PERIOD SUB-RANGES
# =============================================================================
period_data = {}
for el in elements:
    Z = el['Z']
    if Z > 118 or el['block'] != 'd':
        continue
    per = el['period']
    mass = el['mass_amu']
    log_nu = np.log10(compton_freq(mass))

    if per not in period_data:
        period_data[per] = {'elements': [], 'min_log': log_nu, 'max_log': log_nu}
    else:
        period_data[per]['min_log'] = min(period_data[per]['min_log'], log_nu)
        period_data[per]['max_log'] = max(period_data[per]['max_log'], log_nu)

    period_data[per]['elements'].append({
        'Z': Z, 'sym': el['symbol'], 'log_nu': log_nu,
        'd_pos': el.get('d_pos', 0), 'mass': mass,
        'arch': el['archetype'], 'period': per,
    })

# Sort elements within each period
for per in period_data:
    period_data[per]['elements'].sort(key=lambda x: x['Z'])
    period_data[per]['range'] = period_data[per]['max_log'] - period_data[per]['min_log']

# =============================================================================
# THE GLOBAL ACCELERATION CURVE
# =============================================================================
print("=" * 90)
print("GLOBAL ACCELERATION CURVE (from energy quadratic)")
print("=" * 90)
print()
print("  slope(d_eff) = 2 × 0.1964 × d_eff + 8.0932")
print("  curvature    = 2 × 0.1964 = 0.3928 (constant)")
print("  10^0.3928    = 2.471 ≈ first Fibonacci ratio")
print()

# Compute slope at key d_eff values
print(f"  {'d_eff':>6}  {'slope':>8}  {'Δslope from 3.5':>16}  {'Context':>20}")
print("  " + "-" * 60)
for d in [3.0, 3.2, 3.4, 3.50, 3.52, 3.54, 3.56, 3.58, 3.60, 4.0]:
    slope = 2 * 0.1964 * d + 8.0932
    delta = slope - (2 * 0.1964 * 3.5 + 8.0932)
    context = ""
    if abs(d - 3.0) < 0.01: context = "3D→4D boundary"
    elif abs(d - 3.5) < 0.01: context = "mid-interior"
    elif abs(d - 4.0) < 0.01: context = "4D→5D boundary"
    elif abs(d - 3.56) < 0.01: context = "Period 6 d-block"
    print(f"  {d:6.2f}  {slope:8.4f}  {delta:+16.4f}  {context:>20}")

# =============================================================================
# THE LOCAL RAMP — WITHIN EACH PERIOD
# =============================================================================
print()
print("=" * 90)
print("LOCAL ACCELERATION RAMP — WITHIN-PERIOD POSITION AND SWEEP FRACTION")
print("=" * 90)
print()
print("The 'ramp fraction' is how far through the period's sub-range.")
print("The 'sweep fraction' is how far PAST the 29% bulge (0 = at bulge, >0 = past).")
print("Elements past the bulge have electrons being caught by the ramp.")
print()

all_ramp_data = []

for per in sorted(period_data.keys()):
    pd = period_data[per]
    sub_range = pd['range']
    bulge_log = pd['min_log'] + phi_c * sub_range

    print(f"PERIOD {per} d-block:")
    print(f"  Sub-range: {pd['min_log']:.6f} to {pd['max_log']:.6f} ({sub_range:.6f} decades)")
    print(f"  Bulge at:  {bulge_log:.6f} (29% of sub-range)")
    print()

    print(f"  {'Z':>3} {'Sym':>4} {'d_pos':>5} {'ramp_frac':>10} {'past_bulge':>10} "
          f"{'sweep_frac':>11} {'NIST_SO':>8} {'need_corr':>9} {'Arch':>6}")
    print(f"  " + "-" * 80)

    for el in pd['elements']:
        Z = el['Z']
        log_nu = el['log_nu']

        # Ramp fraction (0 to 1 within period)
        ramp_frac = (log_nu - pd['min_log']) / sub_range if sub_range > 0 else 0

        # Past the bulge?
        past = ramp_frac > phi_c

        # Sweep fraction: how far past the bulge (0 at bulge, normalized)
        if past:
            sweep_frac = (ramp_frac - phi_c) / (1.0 - phi_c)
        else:
            sweep_frac = 0.0

        nist = NIST_SO.get(Z, '')
        corr = '←CORR' if Z in need_correction else ''

        print(f"  {Z:3d} {el['sym']:>4s} {el['d_pos']:5d} {ramp_frac:10.4f} "
              f"{'YES' if past else 'no':>10s} {sweep_frac:11.4f} "
              f"{str(nist):>8s} {corr:>9s} {el['arch']:>6s}")

        all_ramp_data.append({
            'Z': Z, 'sym': el['sym'], 'period': per, 'd_pos': el['d_pos'],
            'ramp_frac': ramp_frac, 'past_bulge': past,
            'sweep_frac': sweep_frac,
            'nist_so': NIST_SO.get(Z, None),
            'needs_correction': Z in need_correction,
        })

    print()

# =============================================================================
# CORRELATION: SWEEP FRACTION vs NIST SO
# =============================================================================
print("=" * 90)
print("CORRELATION CHECK: SWEEP FRACTION vs NIST SO (for comparison, not input)")
print("=" * 90)
print()
print("Does the sweep fraction (computed from Z alone) track with")
print("measured SO coupling (from NIST)?")
print()

print(f"{'Z':>3} {'Sym':>4} {'Per':>3} {'d_pos':>5} {'sweep_frac':>11} {'NIST_SO':>8} {'corr':>5}")
print("-" * 50)

sweep_vs_so = []
for rd in all_ramp_data:
    if rd['nist_so'] is not None:
        print(f"{rd['Z']:3d} {rd['sym']:>4s} {rd['period']:3d} {rd['d_pos']:5d} "
              f"{rd['sweep_frac']:11.4f} {rd['nist_so']:8d} "
              f"{'YES' if rd['needs_correction'] else '':>5s}")
        sweep_vs_so.append((rd['sweep_frac'], rd['nist_so'], rd['needs_correction']))

# Compute correlation
if sweep_vs_so:
    sweeps = np.array([x[0] for x in sweep_vs_so])
    sos = np.array([x[1] for x in sweep_vs_so])

    # Pearson correlation
    if np.std(sweeps) > 0 and np.std(sos) > 0:
        corr = np.corrcoef(sweeps, sos)[0, 1]
        print(f"\nPearson correlation (sweep_frac vs NIST_SO): r = {corr:.4f}")
    else:
        print("\nCannot compute correlation (zero variance)")

    # Check: does sweep_frac SEPARATE corrected from uncorrected?
    corr_sweeps = [s for s, so, c in sweep_vs_so if c]
    uncorr_sweeps = [s for s, so, c in sweep_vs_so if not c]

    if corr_sweeps and uncorr_sweeps:
        print(f"\nCorrected elements:   sweep_frac = {min(corr_sweeps):.4f} to {max(corr_sweeps):.4f}")
        print(f"Uncorrected elements: sweep_frac = {min(uncorr_sweeps):.4f} to {max(uncorr_sweeps):.4f}")
        gap = min(corr_sweeps) - max(uncorr_sweeps)
        if gap > 0:
            print(f"CLEAN SEPARATION: gap = {gap:.4f}")
        else:
            print(f"OVERLAP: {-gap:.4f}")
            # But check within each period
            print(f"\n  Within-period check:")
            for per in sorted(period_data.keys()):
                c_sw = [rd['sweep_frac'] for rd in all_ramp_data
                        if rd['period'] == per and rd['needs_correction']]
                u_sw = [rd['sweep_frac'] for rd in all_ramp_data
                        if rd['period'] == per and not rd['needs_correction']
                        and rd['nist_so'] is not None]
                if c_sw and u_sw:
                    gap_p = min(c_sw) - max(u_sw)
                    print(f"  Period {per}: corrected [{min(c_sw):.4f}-{max(c_sw):.4f}] "
                          f"uncorrected [{min(u_sw):.4f}-{max(u_sw):.4f}] "
                          f"gap={gap_p:+.4f} {'SEPARATED' if gap_p > 0 else 'OVERLAP'}")
                elif c_sw:
                    print(f"  Period {per}: corrected [{min(c_sw):.4f}-{max(c_sw):.4f}], "
                          f"no uncorrected with NIST data")
                elif u_sw:
                    print(f"  Period {per}: no corrected, "
                          f"uncorrected [{min(u_sw):.4f}-{max(u_sw):.4f}]")


# =============================================================================
# THE KEY PREDICTION: What sweep_frac threshold separates corrections?
# =============================================================================
print()
print("=" * 90)
print("THRESHOLD SEARCH: What sweep_frac value best separates corrected from uncorrected?")
print("=" * 90)
print()

# Test thresholds from 0 to 1
best_threshold = None
best_accuracy = 0

for threshold in np.arange(0.0, 1.0, 0.01):
    correct = 0
    total = 0
    for rd in all_ramp_data:
        if rd['period'] == 4:
            continue  # Period 4 is ground state, skip
        predicted_corr = rd['sweep_frac'] > threshold
        actual_corr = rd['needs_correction']
        if predicted_corr == actual_corr:
            correct += 1
        total += 1

    acc = correct / total if total > 0 else 0
    if acc > best_accuracy:
        best_accuracy = acc
        best_threshold = threshold

print(f"Best threshold (Periods 5+6+7): sweep_frac > {best_threshold:.2f}")
print(f"Accuracy at this threshold: {best_accuracy:.1%}")
print()

# Show what this threshold predicts
print(f"At threshold = {best_threshold:.2f}:")
for per in sorted(period_data.keys()):
    if per == 4:
        continue
    print(f"\n  Period {per}:")
    for rd in all_ramp_data:
        if rd['period'] != per:
            continue
        predicted = rd['sweep_frac'] > best_threshold
        actual = rd['needs_correction']
        match = "✓" if predicted == actual else "✗"
        pred_label = "CORR" if predicted else "stable"
        act_label = "CORR" if actual else "stable"
        print(f"    {rd['sym']:>4s} d{rd['d_pos']}: sweep={rd['sweep_frac']:.4f} "
              f"→ pred={pred_label:>6s} actual={act_label:>6s} {match}")