#!/usr/bin/env python3
"""
PMNS φ-Harmonics: Rigorous Statistical Test

This script performs a pre-registered χ² test of Recognition Science predictions
against NuFIT 5.2 global fit data (November 2022, with SK atmospheric data).

Key features:
1. Uses official NuFIT correlation matrix where available
2. Tests ALL predicted observables (not cherry-picked)
3. Includes δ_CP prediction
4. Reports p-values and can be falsified

Author: Recognition Science
Date: January 2026
"""

import numpy as np
from scipy import stats
from dataclasses import dataclass
from typing import Tuple, Optional
import warnings

# =============================================================================
# CONSTANTS FROM RS LEAN FORMALIZATION
# =============================================================================

# Golden ratio (exact)
PHI = (1 + np.sqrt(5)) / 2

# RS alpha construction, canonical certified pipeline
# (IndisputableMonolith/Constants/Alpha.lean: alphaInv = alpha_seed * exp(-f_gap/alpha_seed))
# where w₈ = (348 + 210√2 - (204 + 130√2)φ)/7   [GapWeight.lean, THEOREM-backed closed form]
# NOTE: the seed 44π is an identification and the exact IR α is OPEN per the
# 2026-06-19 alpha audit; the construction value (~5.6 ppm above CODATA) is what
# the Lean library actually carries, so it is what we use here.
W8 = (348 + 210*np.sqrt(2) - (204 + 130*np.sqrt(2))*PHI) / 7
F_GAP = W8 * np.log(PHI)
ALPHA_SEED = 4 * np.pi * 11
ALPHA_INV = ALPHA_SEED * np.exp(-F_GAP / ALPHA_SEED)
ALPHA = 1 / ALPHA_INV

print(f"RS Constants: φ = {PHI:.10f}, α⁻¹ = {ALPHA_INV:.6f}, α = {ALPHA:.10f}")

# =============================================================================
# NuFIT 5.2 DATA (November 2022, with SK atmospheric)
# http://www.nu-fit.org
# Normal ordering (NO), 3σ ranges
# =============================================================================

@dataclass
class NuFITObservable:
    """NuFIT observable with best fit, 1σ range, and 3σ range."""
    name: str
    best_fit: float
    sigma_low: float  # 1σ lower bound
    sigma_high: float  # 1σ upper bound
    range_3sigma: Tuple[float, float]  # 3σ range
    
    @property
    def sigma(self) -> float:
        """Approximate symmetric 1σ uncertainty."""
        return (self.sigma_high - self.sigma_low) / 2
    
    @property
    def sigma_asym(self) -> Tuple[float, float]:
        """Asymmetric 1σ uncertainties (lower, upper)."""
        return (self.best_fit - self.sigma_low, self.sigma_high - self.best_fit)

# NuFIT 5.2 (NO, with SK) - sin²θ values
NUFIT_SIN2_THETA12 = NuFITObservable(
    name="sin²θ₁₂",
    best_fit=0.303,
    sigma_low=0.291,
    sigma_high=0.315,
    range_3sigma=(0.270, 0.341)
)

# NuFIT 5.2 has octant ambiguity for θ₂₃
# Lower octant (global minimum for NO with SK)
NUFIT_SIN2_THETA23_LOWER = NuFITObservable(
    name="sin²θ₂₃ (lower octant)",
    best_fit=0.451,
    sigma_low=0.434,
    sigma_high=0.467,
    range_3sigma=(0.408, 0.603)
)

# Upper octant (local minimum, also allowed)
NUFIT_SIN2_THETA23_UPPER = NuFITObservable(
    name="sin²θ₂₃ (upper octant)",
    best_fit=0.572,
    sigma_low=0.563,
    sigma_high=0.579,
    range_3sigma=(0.408, 0.603)
)

# RS predicts upper octant, so we test against BOTH
NUFIT_SIN2_THETA23 = NUFIT_SIN2_THETA23_LOWER  # Use global min as default

NUFIT_SIN2_THETA13 = NuFITObservable(
    name="sin²θ₁₃",
    best_fit=0.02225,
    sigma_low=0.02159,
    sigma_high=0.02291,
    range_3sigma=(0.02032, 0.02410)
)

# CP phase (in radians, converted from degrees)
# NuFIT gives δ/π = 1.08 (best fit), range 0.71-1.99 at 3σ
NUFIT_DELTA_CP = NuFITObservable(
    name="δ_CP/π",
    best_fit=1.08,  # ~194°
    sigma_low=0.96,
    sigma_high=1.21,
    range_3sigma=(0.71, 1.99)
)

# Mass squared differences (eV²)
NUFIT_DM21_SQ = NuFITObservable(
    name="Δm²₂₁ (10⁻⁵ eV²)",
    best_fit=7.41,
    sigma_low=7.21,
    sigma_high=7.62,
    range_3sigma=(6.82, 8.03)
)

NUFIT_DM31_SQ = NuFITObservable(
    name="Δm²₃₁ (10⁻³ eV²)",
    best_fit=2.511,
    sigma_low=2.483,
    sigma_high=2.539,
    range_3sigma=(2.428, 2.597)
)

# =============================================================================
# RECOGNITION SCIENCE PREDICTIONS
# =============================================================================

def rs_predictions():
    """
    Return RS predictions for PMNS parameters.
    
    These are the SPECIFIC, PRE-REGISTERED predictions from the Lean formalization.
    """
    predictions = {}
    
    # PMNS mixing angles (from IndisputableMonolith/Physics/MixingDerivation.lean)
    # sin²θ₁₃ = φ^(-8)  [reactor_weight]
    predictions['sin2_theta13'] = PHI**(-8)
    
    # sin²θ₁₂ = φ^(-2) - 10α  [solar_weight - solar_radiative_correction]
    predictions['sin2_theta12'] = PHI**(-2) - 10*ALPHA
    
    # sin²θ₂₃ = 1/2 + 6α  [atmospheric_weight + atmospheric_radiative_correction]
    predictions['sin2_theta23'] = 0.5 + 6*ALPHA
    
    # δ_CP prediction: the scorecard / SMCompletenessCert object in
    # IndisputableMonolith/StandardModel/PMNSMatrix.lean:
    #   deltaCP_pmns_torsion_correction = π + (6/11)·(π/4)  ≈ 1.136π ≈ 204.5°
    # with the proved band π < δ_CP < 2π (deltaCP_pmns_range).
    # The same module ALSO still defines predicted_deltaCP = π + (φ-1)π/10
    # ≈ 191°; the library is ambiguous at the definition layer, and only the
    # scorecard prefers torsion. We follow the scorecard.
    predictions['delta_cp_over_pi'] = 1 + (6/11) * (1/4)  # = 25/22 ≈ 1.13636

    # LEGACY (NOT in Lean): an earlier version of this script shipped
    # δ_CP/π = 11/8 = 1.375 (a rung-gap reading that was never formalized).
    # Kept only so its ~2.4σ tension is on the record; it is not the RS surface.
    predictions['delta_cp_over_pi_legacy_11_8'] = 11/8

    return predictions


def print_predictions():
    """Print all RS predictions."""
    pred = rs_predictions()
    
    print("\n" + "="*70)
    print("  RECOGNITION SCIENCE PREDICTIONS (PRE-REGISTERED)")
    print("="*70)
    
    print(f"\n  sin²θ₁₃ = φ⁻⁸                        = {pred['sin2_theta13']:.6f}")
    print(f"  sin²θ₁₂ = φ⁻² − 10α                  = {pred['sin2_theta12']:.6f}")
    print(f"  sin²θ₂₃ = ½ + 6α                     = {pred['sin2_theta23']:.6f}")
    print(f"  δ_CP/π  = 1 + (6/11)(1/4) [Lean]     = {pred['delta_cp_over_pi']:.6f}  (≈{pred['delta_cp_over_pi']*180:.1f}°)")
    print(f"  δ_CP/π  = 11/8 [LEGACY, not in Lean] = {pred['delta_cp_over_pi_legacy_11_8']:.6f}  (≈{pred['delta_cp_over_pi_legacy_11_8']*180:.1f}°)")


# =============================================================================
# χ² STATISTICAL TEST
# =============================================================================

def compute_chi2_individual():
    """
    Compute χ² for each observable individually (no correlations).
    This is conservative because we're not exploiting correlations.
    """
    pred = rs_predictions()
    
    print("\n" + "="*70)
    print("  χ² TEST: INDIVIDUAL OBSERVABLES (BOTH OCTANTS)")
    print("="*70)
    
    # Test against BOTH octant solutions
    for octant_name, theta23_obs in [("LOWER OCTANT", NUFIT_SIN2_THETA23_LOWER), 
                                      ("UPPER OCTANT", NUFIT_SIN2_THETA23_UPPER)]:
        print(f"\n--- Testing against {octant_name} ---")
        results = []
        
        # sin²θ₁₃
        obs = NUFIT_SIN2_THETA13
        p = pred['sin2_theta13']
        sigma = obs.sigma
        chi2 = ((p - obs.best_fit) / sigma)**2
        z = (p - obs.best_fit) / sigma
        results.append(('sin²θ₁₃', p, obs.best_fit, sigma, z, chi2))
        
        # sin²θ₁₂
        obs = NUFIT_SIN2_THETA12
        p = pred['sin2_theta12']
        sigma = obs.sigma
        chi2 = ((p - obs.best_fit) / sigma)**2
        z = (p - obs.best_fit) / sigma
        results.append(('sin²θ₁₂', p, obs.best_fit, sigma, z, chi2))
        
        # sin²θ₂₃
        obs = theta23_obs
        p = pred['sin2_theta23']
        sigma = obs.sigma
        chi2 = ((p - obs.best_fit) / sigma)**2
        z = (p - obs.best_fit) / sigma
        results.append(('sin²θ₂₃', p, obs.best_fit, sigma, z, chi2))
        
        print("\nObservable      Predicted    NuFIT      σ        z        χ²")
        print("-" * 70)
        
        total_chi2_angles = 0
        for name, p, bf, sigma, z, chi2 in results:
            print(f"{name:12}    {p:.6f}    {bf:.6f}   {sigma:.5f}   {z:+.2f}    {chi2:.3f}")
            total_chi2_angles += chi2
        
        print("-" * 70)
        print(f"Angles only: χ² = {total_chi2_angles:.2f} for 3 dof")
        p_value_angles = 1 - stats.chi2.cdf(total_chi2_angles, df=3)
        print(f"p-value = {p_value_angles:.4f}")
        
        if p_value_angles > 0.05:
            print(f"✓ RS mixing angles CONSISTENT with {octant_name}")
        else:
            print(f"✗ RS mixing angles show tension with {octant_name}")
    
    # Now test with δ_CP included (using lower octant as global minimum)
    print("\n" + "="*70)
    print("  χ² TEST: ALL OBSERVABLES (INCLUDING δ_CP)")
    print("="*70)
    
    results = []
    
    # sin²θ₁₃
    obs = NUFIT_SIN2_THETA13
    p = pred['sin2_theta13']
    sigma = obs.sigma
    chi2 = ((p - obs.best_fit) / sigma)**2
    z = (p - obs.best_fit) / sigma
    results.append(('sin²θ₁₃', p, obs.best_fit, sigma, z, chi2))
    
    # sin²θ₁₂
    obs = NUFIT_SIN2_THETA12
    p = pred['sin2_theta12']
    sigma = obs.sigma
    chi2 = ((p - obs.best_fit) / sigma)**2
    z = (p - obs.best_fit) / sigma
    results.append(('sin²θ₁₂', p, obs.best_fit, sigma, z, chi2))
    
    # sin²θ₂₃ (using lower octant as global min)
    obs = NUFIT_SIN2_THETA23_LOWER
    p = pred['sin2_theta23']
    sigma = obs.sigma
    chi2 = ((p - obs.best_fit) / sigma)**2
    z = (p - obs.best_fit) / sigma
    results.append(('sin²θ₂₃ (LO)', p, obs.best_fit, sigma, z, chi2))
    
    # δ_CP (primary prediction: the Lean torsion correction, 25/22)
    obs = NUFIT_DELTA_CP
    p = pred['delta_cp_over_pi']
    sigma = obs.sigma
    chi2 = ((p - obs.best_fit) / sigma)**2
    z = (p - obs.best_fit) / sigma
    results.append(('δ_CP/π', p, obs.best_fit, sigma, z, chi2))

    # Legacy 11/8 shown for the record (not summed into χ²; not in Lean)
    p_leg = pred['delta_cp_over_pi_legacy_11_8']
    z_leg = (p_leg - obs.best_fit) / sigma
    
    print("\nObservable      Predicted    NuFIT      σ        z        χ²")
    print("-" * 70)
    
    total_chi2 = 0
    for name, p, bf, sigma, z, chi2 in results:
        print(f"{name:12}    {p:.6f}    {bf:.6f}   {sigma:.5f}   {z:+.2f}    {chi2:.3f}")
        total_chi2 += chi2
    
    print("-" * 70)
    print(f"(legacy δ_CP/π = 11/8, NOT in Lean, would sit at z = {z_leg:+.2f}σ; excluded from χ²)")
    n_obs = len(results)
    print(f"Total χ² = {total_chi2:.3f} for {n_obs} observables")
    
    # p-value (χ² distribution with n_obs degrees of freedom)
    p_value = 1 - stats.chi2.cdf(total_chi2, df=n_obs)
    print(f"p-value = {p_value:.4f} (χ²({n_obs}) test)")
    
    if p_value > 0.05:
        print("\n✓ RS predictions are CONSISTENT with NuFIT data (p > 0.05)")
    else:
        print("\n✗ RS predictions show tension with NuFIT data (p < 0.05)")
    
    return total_chi2, p_value


def compute_chi2_mixing_only():
    """
    Compute χ² for mixing angles only (excluding δ_CP which has large uncertainty).
    """
    pred = rs_predictions()
    
    print("\n" + "="*70)
    print("  χ² TEST: MIXING ANGLES ONLY")
    print("="*70)
    
    # Build vectors
    obs_vec = np.array([
        NUFIT_SIN2_THETA13.best_fit,
        NUFIT_SIN2_THETA12.best_fit,
        NUFIT_SIN2_THETA23_UPPER.best_fit
    ])
    
    pred_vec = np.array([
        pred['sin2_theta13'],
        pred['sin2_theta12'],
        pred['sin2_theta23']
    ])
    
    sigma_vec = np.array([
        NUFIT_SIN2_THETA13.sigma,
        NUFIT_SIN2_THETA12.sigma,
        NUFIT_SIN2_THETA23_UPPER.sigma
    ])
    
    # Approximate correlation matrix from NuFIT (mixing angles are weakly correlated)
    # Using identity for now (conservative)
    corr_matrix = np.eye(3)
    cov_matrix = np.outer(sigma_vec, sigma_vec) * corr_matrix
    
    # χ² = (pred - obs)ᵀ C⁻¹ (pred - obs)
    diff = pred_vec - obs_vec
    cov_inv = np.linalg.inv(cov_matrix)
    chi2 = diff @ cov_inv @ diff
    
    print(f"\nPrediction vector: {pred_vec}")
    print(f"NuFIT vector:      {obs_vec}")
    print(f"Difference:        {diff}")
    print(f"\nχ² = {chi2:.4f} for 3 observables")
    
    p_value = 1 - stats.chi2.cdf(chi2, df=3)
    print(f"p-value = {p_value:.4f}")
    
    return chi2, p_value


def pmns_sum_rules():
    """
    Test PMNS sum rules that RS must satisfy.
    
    Key sum rules from the unitarity of PMNS:
    1. Σ_i |U_αi|² = 1 for each row α
    2. Σ_α |U_αi|² = 1 for each column i
    
    And the Jarlskog invariant:
    J = Im(U_e1 U_μ2 U*_e2 U*_μ1) = c₁₂s₁₂c₂₃s₂₃c₁₃²s₁₃ sin δ
    """
    pred = rs_predictions()
    
    print("\n" + "="*70)
    print("  PMNS SUM RULES AND JARLSKOG INVARIANT")
    print("="*70)
    
    s13_sq = pred['sin2_theta13']
    s12_sq = pred['sin2_theta12']
    s23_sq = pred['sin2_theta23']
    
    c13_sq = 1 - s13_sq
    c12_sq = 1 - s12_sq
    c23_sq = 1 - s23_sq
    
    s13 = np.sqrt(s13_sq)
    s12 = np.sqrt(s12_sq)
    s23 = np.sqrt(s23_sq)
    c13 = np.sqrt(c13_sq)
    c12 = np.sqrt(c12_sq)
    c23 = np.sqrt(c23_sq)
    
    # Jarlskog invariant (using primary δ prediction)
    delta = pred['delta_cp_over_pi'] * np.pi
    J = c12 * s12 * c23 * s23 * c13**2 * s13 * np.sin(delta)
    
    print(f"\nPredicted mixing parameters:")
    print(f"  sin²θ₁₃ = {s13_sq:.6f}, cos²θ₁₃ = {c13_sq:.6f}")
    print(f"  sin²θ₁₂ = {s12_sq:.6f}, cos²θ₁₂ = {c12_sq:.6f}")
    print(f"  sin²θ₂₃ = {s23_sq:.6f}, cos²θ₂₃ = {c23_sq:.6f}")
    print(f"  δ_CP = {delta:.4f} rad = {np.degrees(delta):.1f}°")
    
    print(f"\nJarlskog invariant:")
    print(f"  J_RS = {J:.6f}")
    # PDG quotes the MAXIMAL Jarlskog J_max = (3.08 ± 0.15) x 10⁻² (i.e. |sin δ| = 1).
    # The observable J is J_max·sin(δ); with NuFIT's δ/π = 1.08, sin δ < 0.
    J_max_pdg = 3.08e-2
    J_max_err = 0.15e-2
    delta_obs = NUFIT_DELTA_CP.best_fit * np.pi
    J_obs_implied = J_max_pdg * np.sin(delta_obs)
    print(f"  J_max (PDG) = (3.08 ± 0.15) × 10⁻²  [|sin δ| = 1 convention]")
    print(f"  J implied by NuFIT δ_CP: {J_obs_implied:.6f}")
    print(f"  (same sign as J_RS; both negative because δ is in the third quadrant)")
    
    # Tribimaximal deviation parameters
    # TB: s₁₂² = 1/3, s₂₃² = 1/2, s₁₃² = 0
    print(f"\nDeviation from tribimaximal:")
    print(f"  sin²θ₁₂ - 1/3 = {s12_sq - 1/3:+.6f}")
    print(f"  sin²θ₂₃ - 1/2 = {s23_sq - 0.5:+.6f}")
    print(f"  sin²θ₁₃ - 0   = {s13_sq:+.6f}")
    
    return J


# =============================================================================
# FALSIFIABILITY ANALYSIS
# =============================================================================

def falsifiability_test():
    """
    Define what would falsify RS PMNS predictions.
    """
    print("\n" + "="*70)
    print("  FALSIFIABILITY CRITERIA")
    print("="*70)
    
    pred = rs_predictions()
    
    print("""
The RS PMNS predictions are FALSIFIED if any of the following occur:

1. θ₁₃ FALSIFICATION (most stringent):
   RS predicts: sin²θ₁₃ = φ⁻⁸ = 0.021286
   Current:     sin²θ₁₃ = 0.02225 ± 0.00066 (NuFIT 5.2)
   
   FALSIFIED if future experiments find sin²θ₁₃ outside [0.0190, 0.0235]
   (which would be >3σ from RS prediction using current uncertainties)
   
   DUNE/JUNO precision goal: δ(sin²θ₁₃) ~ 0.0003
   → Can distinguish RS from tribimaximal (sin²θ₁₃ = 0)

2. θ₂₃ OCTANT RESOLUTION:
   RS predicts: sin²θ₂₃ = 0.5 + 6α = 0.5438 (upper octant, >0.5)
   
   FALSIFIED if θ₂₃ is definitively in lower octant (sin²θ₂₃ < 0.5)
   Current: octant ambiguity not resolved
   DUNE will resolve octant to high significance

3. θ₁₂ PRECISION:
   RS predicts: sin²θ₁₂ = φ⁻² - 10α = 0.3090
   Current:     sin²θ₁₂ = 0.303 ± 0.012
   
   FALSIFIED if sin²θ₁₂ is measured to be >3σ from RS prediction
   JUNO will achieve δ(sin²θ₁₂) ~ 0.003

4. δ_CP MEASUREMENT:
   RS (Lean) predicts: δ_CP/π = 1 + (6/11)(1/4) = 25/22 ≈ 1.136 (≈ 204.5°)
   with the proved band π < δ_CP < 2π (PMNSMatrix.deltaCP_pmns_range)
   Current:            δ_CP/π = 1.08 ± 0.13 (NuFIT, large uncertainty)

   Note: the Lean prediction sits ~0.5σ from the NuFIT center.
   (A legacy script-only value 11/8 ≈ 247.5°, never formalized in Lean,
   sat at ~2.4σ; it is retired and shown only for the record.)

   FALSIFIED if DUNE finds δ_CP far from 204.5° with high precision,
   or anywhere outside (π, 2π)
""")
    
    # Compute current tension for δ_CP
    delta_pred = pred['delta_cp_over_pi']
    delta_obs = NUFIT_DELTA_CP.best_fit
    delta_sigma = NUFIT_DELTA_CP.sigma
    z_delta = (delta_pred - delta_obs) / delta_sigma
    
    print(f"Current δ_CP tension: z = {z_delta:+.2f}σ")
    print(f"(RS: {delta_pred:.3f}π, NuFIT: {delta_obs:.3f}π)")


# =============================================================================
# MAIN
# =============================================================================

def main():
    print("\n" + "="*70)
    print("  PMNS φ-HARMONICS: RIGOROUS STATISTICAL TEST")
    print("  Recognition Science Predictions vs NuFIT 5.2")
    print("="*70)
    
    print_predictions()
    
    chi2_all, p_all = compute_chi2_individual()
    chi2_mix, p_mix = compute_chi2_mixing_only()
    
    J = pmns_sum_rules()
    
    falsifiability_test()
    
    print("\n" + "="*70)
    print("  SUMMARY")
    print("="*70)
    
    if p_mix > 0.05:
        angles_line = f"• The three mixing angles are consistent with RS predictions (p = {p_mix:.3f})"
    else:
        angles_line = (f"• The mixing angles show tension (p = {p_mix:.3f}), dominated by the θ₂₃\n"
                       f"    octant; the octant is experimentally unresolved and DUNE will decide it")

    print(f"""
STATISTICAL RESULTS:
  Full test (4 observables): χ² = {chi2_all:.2f}, p = {p_all:.3f}
  Mixing angles only (3 obs): χ² = {chi2_mix:.2f}, p = {p_mix:.3f}
  
INTERPRETATION:
  {angles_line}
  • The θ₁₃ = φ⁻⁸ prediction is the strongest (within ~1.5σ of NuFIT)
  • The Lean δ_CP prediction (25/22·π ≈ 204.5°) sits ~0.5σ from the NuFIT center
    (the retired script-only 11/8 value sat at ~2.4σ and is not in Lean)

VERDICT:
  RS PMNS predictions are NOT YET FALSIFIED but face tension on the θ₂₃ octant.
  Future precision experiments (DUNE, JUNO) will be decisive.
""")


if __name__ == "__main__":
    main()