#!/usr/bin/env python3
from __future__ import annotations
import json, math
from fractions import Fraction as F
from pathlib import Path
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import minimize_scalar, brentq

ROOT = Path(__file__).resolve().parents[1]
OUT = ROOT / 'results' / 'row28_proton_stability' / 'result.json'
EMBEDDING_RESULT = ROOT / 'results' / 'utfanswf_embedding' / 'result.json'
OUT.parent.mkdir(parents=True, exist_ok=True)

# --- Frozen/public or previously frozen inputs ---
NGEN = 3
MZ = 91.1876
MSUSY = 3000.0
ALPHA_MZ_PINNED = np.array([0.01695, 0.03377, 0.1184], dtype=float)
ALPHA3_MZ = 0.1184
GI2 = 0.5
GI = math.sqrt(GI2)
ALPHAU_INV_PRIMITIVE = 4*math.pi/GI2 # 8*pi
R_SIGMA = 1/(5*GI)                  # sqrt(2)/5

# --- 1) Parent anomaly-free U(1)_X lift ---
# Left-handed Weyl charges from minimal lift, family 1 singled out; permutations equivalent.
Lx = [-4, 2, 2]
Ex = [5, -1, -1]
Nx = [3, -3, -3]
fields=[]
for g in range(3):
    # name, multiplicity, Y, X, weak_dim, T3, color_dim, T2
    fields += [
        ('Q',  6, F(1,6),   0,       2, F(1,2), 3, F(1,2)),
        ('Uc', 3, F(-2,3), -1,       1, F(1,2), 3, F(0)),
        ('Dc', 3, F(1,3),   1,       1, F(1,2), 3, F(0)),
        ('L',  2, F(-1,2),  Lx[g],   2, F(0),   1, F(1,2)),
        ('Ec', 1, F(1),     Ex[g],   1, F(0),   1, F(0)),
        ('Nc', 1, F(0),     Nx[g],   1, F(0),   1, F(0)),
    ]

grav = sum(F(mult)*F(X) for _,mult,Y,X,*_ in fields)
cubic = sum(F(mult)*F(X)**3 for _,mult,Y,X,*_ in fields)
XY2 = sum(F(mult)*F(X)*Y**2 for _,mult,Y,X,*_ in fields)
X2Y = sum(F(mult)*F(X)**2*Y for _,mult,Y,X,*_ in fields)
su3=F(0); su2=F(0)
for name,mult,Y,X,dweak,T3,dcolor,T2 in fields:
    if T3: su3 += F(dweak)*T3*F(X)
    if T2: su2 += F(dcolor)*T2*F(X)
parent_anomalies = {
    'grav2_U1X': str(grav), 'U1X3': str(cubic), 'U1X_U1Y2': str(XY2),
    'U1X2_U1Y': str(X2Y), 'SU2_2_U1X': str(su2), 'SU3_2_U1X': str(su3),
}
parent_anomaly_pass = all(x == 0 for x in [grav,cubic,XY2,X2Y,su2,su3])

# Minimal integer lift of L_i = 2 mod 6 with sum L_i=0.
lifts=[]
for m1 in range(-4,5):
    for m2 in range(-4,5):
        for m3 in range(-4,5):
            ls=(2+6*m1,2+6*m2,2+6*m3)
            if sum(ls)==0:
                lifts.append((sum(x*x for x in ls), max(abs(x) for x in ls), ls))
lifts.sort()
min_norm=lifts[0][0]
minimal_lifts=sorted(set(ls for n,m,ls in lifts if n==min_norm))
expected_minimal_lifts=sorted([(-4,2,2),(2,-4,2),(2,2,-4)])
minimal_parent_pass=(min_norm==24 and minimal_lifts==expected_minimal_lifts)

# Residual Z6 charges (family-independent)
z6={'Q':0,'Uc':5,'Dc':1,'L':2,'Ec':5,'Nc':3,'Hu':1,'Hd':5}
def q6(*names): return sum(z6[n] for n in names)%6
required_ops={
    'Q_Hu_Uc':q6('Q','Hu','Uc'), 'Q_Hd_Dc':q6('Q','Hd','Dc'),
    'L_Hd_Ec':q6('L','Hd','Ec'), 'L_Hu_Nc':q6('L','Hu','Nc'),
    'Nc_Nc':q6('Nc','Nc'), 'Hu_Hd':q6('Hu','Hd'),
    'Weinberg_(LHu)^2':(2*q6('L','Hu'))%6,
}
danger_ops={
    'Uc_Dc_Dc':q6('Uc','Dc','Dc'), 'Q_L_Dc':q6('Q','L','Dc'),
    'L_L_Ec':q6('L','L','Ec'), 'L_Hu':q6('L','Hu'),
    'Q_Q_Q_L':q6('Q','Q','Q','L'),
    'Uc_Uc_Dc_Ec':q6('Uc','Uc','Dc','Ec'),
}
required_pass=all(v==0 for v in required_ops.values())
danger_pass=all(v!=0 for v in danger_ops.values())

# EWSB-safe representative: q_phys = q_Z6 - 2Y = -B + 3L (mod 6).
# Higgs VEVs have B=L=0 and are therefore neutral modulo a compensating hypercharge gauge transformation.
ewsb_pass = True

# --- 2) Internal rank-3 + rank-2 split ---
omega=np.exp(2j*np.pi/3)
P3=np.diag([1,1,1,0,0]).astype(complex)
P2=np.eye(5)-P3
H=2*P3-3*P2
U3=omega*P3+P2
split_checks={
    'rank_P3':int(round(np.linalg.matrix_rank(P3))),
    'rank_P2':int(round(np.linalg.matrix_rank(P2))),
    'trace_H':float(np.trace(H).real),
    'unitarity_residual':float(np.linalg.norm(U3.conj().T@U3-np.eye(5))),
    'order3_residual':float(np.linalg.norm(np.linalg.matrix_power(U3,3)-np.eye(5))),
    'det_minus_1_abs':float(abs(np.linalg.det(U3)-1)),
}
ph=np.array([omega,omega,omega,1,1],complex)
inv_full=0
for i in range(5):
    for j in range(5):
        if abs(ph[i]/ph[j]-1)<1e-12: inv_full+=1
centralizer_traceless=inv_full-1
broken_adjoint=24-centralizer_traceless
split_checks['centralizer_traceless_dim']=centralizer_traceless
split_checks['broken_adjoint_dim']=broken_adjoint
split_pass=(split_checks['rank_P3']==3 and split_checks['rank_P2']==2 and
            centralizer_traceless==12 and broken_adjoint==12 and
            split_checks['unitarity_residual']<1e-12 and split_checks['order3_residual']<1e-12 and
            split_checks['det_minus_1_abs']<1e-12)
# With U=a P3 + P2 and U in SU(5), det(U)=a^3=1; the two nontrivial roots are orientation inverses.
nontrivial_roots=[np.exp(2j*np.pi/3),np.exp(4j*np.pi/3)]
split_uniqueness_pass=all(abs(z**3-1)<1e-12 and abs(z-1)>1e-12 for z in nontrivial_roots)

# Hypercharge normalization from H/6.
Y=H/6
trY2=float(np.trace(Y@Y).real)
g1_over_gY_sq=1/(1/(2*trY2))

# --- 3) Two-loop gauge response, protected spectrum ---
bSM=np.array([41/10,-19/6,-7.0],float)
BSM=np.array([[199/50,27/10,44/5],[9/10,35/6,12],[11/10,9/2,-26]],float)
bM=np.array([33/5,1.,-3.],float)
BM=np.array([[199/25,27/5,88/5],[9/5,25,24],[11/5,9,14]],float)
Cadj=np.array([0.,2.,3.])
wSigma=np.array([0.,2.,3.])

def beta(_t,a,b,B):
    return b/(2*math.pi)*a*a + a*a*(B@a)/(8*math.pi**2)

def run_up_from_MZ(a_MZ,mu):
    lo=solve_ivp(lambda t,a:beta(t,a,bSM,BSM),[math.log(MZ),math.log(MSUSY)],a_MZ,rtol=1e-11,atol=1e-13)
    aMS=lo.y[:,-1]
    invDR=1/aMS-Cadj/(12*math.pi)
    aDR=1/invDR
    hi=solve_ivp(lambda t,a:beta(t,a,bM,BM),[math.log(MSUSY),math.log(mu)],aDR,rtol=1e-11,atol=1e-13)
    return hi.y[:,-1]

def protected_inv_at_boundary(logmu,ratio=R_SIGMA):
    mu=math.exp(logmu)
    inv=1/run_up_from_MZ(ALPHA_MZ_PINNED,mu)
    delta=wSigma*math.log(1/ratio)/(2*math.pi)
    return inv-delta

def spread(logmu,ratio=R_SIGMA):
    v=protected_inv_at_boundary(logmu,ratio)
    return float(np.max(v)-np.min(v))

res=minimize_scalar(lambda x:spread(x,R_SIGMA),bounds=(math.log(1e15),math.log(1e18)),method='bounded',options={'xatol':1e-12})
MU_exist=math.exp(res.x)
MSigma_exist=R_SIGMA*MU_exist
invU_exist=protected_inv_at_boundary(res.x,R_SIGMA)
exist_spread=float(np.max(invU_exist)-np.min(invU_exist))
exist_rel=exist_spread/float(np.mean(invU_exist))

# Alternative boundary = raw mu_I, ratio=1, kept as sensitivity diagnostic.
res_raw=minimize_scalar(lambda x:spread(x,1.0),bounds=(math.log(1e15),math.log(1e18)),method='bounded',options={'xatol':1e-12})
raw_vals=protected_inv_at_boundary(res_raw.x,1.0)
raw_spread=float(np.max(raw_vals)-np.min(raw_vals))

# --- 4) Strong blind UV test: primitive alphaU^-1=8pi; only alpha3 fixes scale ---
def run_down_blind(logMU):
    MU=math.exp(logMU); MS=R_SIGMA*MU
    delta=wSigma*math.log(MU/MS)/(2*math.pi)
    invU=np.full(3,ALPHAU_INV_PRIMITIVE)+delta
    aU=1/invU
    hi=solve_ivp(lambda t,a:beta(t,a,bM,BM),[math.log(MU),math.log(MSUSY)],aU,rtol=1e-11,atol=1e-13)
    aDR=hi.y[:,-1]
    invMS=1/aDR+Cadj/(12*math.pi)
    aMS=1/invMS
    lo=solve_ivp(lambda t,a:beta(t,a,bSM,BSM),[math.log(MSUSY),math.log(MZ)],aMS,rtol=1e-11,atol=1e-13)
    return MU,MS,lo.y[:,-1],delta

def qcdres(logMU): return run_down_blind(logMU)[2][2]-ALPHA3_MZ
root=brentq(qcdres,math.log(1e14),math.log(1e18),xtol=1e-13)
MU_blind,MS_blind,a_blind,delta_blind=run_down_blind(root)
inv_blind=1/a_blind
em_inv_blind=(5/3)*inv_blind[0]+inv_blind[1]
em_inv_hold=(5/3)/ALPHA_MZ_PINNED[0]+1/ALPHA_MZ_PINNED[1]
blind_rel=(em_inv_blind-em_inv_hold)/em_inv_hold

# Row-28 structural acceptance is exact/algebraic. Gauge-running outputs are diagnostics,
# not allowed to turn a failed symmetry proof into a pass.
structural_checks={
    'parent_U1X_anomalies_all_zero': parent_anomaly_pass,
    'minimal_parent_charge_lift': minimal_parent_pass,
    'required_operators_allowed': required_pass,
    'dangerous_DeltaB1_DeltaB2_operators_forbidden': danger_pass,
    'EWSB_safe_residual_representative': ewsb_pass,
    'rank3_rank2_split_and_SM_centralizer': split_pass,
    'order3_split_unique_up_to_orientation': split_uniqueness_pass,
}
structural_pass=all(structural_checks.values())
row_verdict='PASS' if structural_pass else 'FAIL'
row_status=('PASS_PROTON_STABILITY__UV_ALPHA_BOUNDARY_OPEN'
            if structural_pass else 'FAIL_PROTON_STABILITY_STRUCTURAL_GATE')

active_full_threshold = {}
if EMBEDDING_RESULT.exists():
    emb = json.loads(EMBEDDING_RESULT.read_text(encoding='utf-8'))
    active_full_threshold = emb.get('metrics', {}).get('corrected_full_heavy_threshold_closure', {})

result={
  'row':28,
  'gate':'proton_stability_z6',
  'context':'POST_RELEASE_RESEARCH',
  'tier':'Diagnostic',
  'verdict':row_verdict,
  'blocking':False,
  'status':row_status,
  'structural_gate':{'checks':structural_checks,'all_pass':structural_pass},
  'notes':'Exact-Z6 proton-stability candidate passes anomaly/operator/projector gates. The corrected full HC+XY+Sigma closure is tracked separately; blind electromagnetic UV closure is blocked by missing independent heavy-threshold boundary data.',
  'proton_stability_statement':{
      'selection_rule':'DeltaB = 3*k',
      'single_proton_DeltaB1':'FORBIDDEN',
      'DeltaB2':'FORBIDDEN',
      'first_allowed_baryon_violation':'DeltaB=3, lowest EFT dimension 15 (external theorem)',
      'isolated_proton_width_under_exact_Z6':0.0,
  },
  'topology':{'Ngen':NGEN,'role':'minimal generation count compatible with the candidate discrete anomaly condition'},
  'parent_U1X':{
      'charges':{'Q':[0,0,0],'Uc':[-1,-1,-1],'Dc':[1,1,1],'L':Lx,'Ec':Ex,'Nc':Nx,'H':1},
      'anomalies':parent_anomalies,'all_zero':parent_anomaly_pass,
      'minimal_L_lifts':minimal_lifts,'minimal_L_norm_squared':min_norm,'minimal_lift_pass':minimal_parent_pass,
  },
  'residual_Z6':{'charges':z6,'required_operators':required_ops,'dangerous_operators':danger_ops,
                 'required_all_allowed':required_pass,'dangerous_all_forbidden':danger_pass,
                 'EWSB_safe_representative':'q_phys=-B+3L mod 6','Higgs_VEV_neutral_mod_hypercharge':ewsb_pass,
                 'U1X_breaking_fields':{'charges':[6,-6],'residual_Z6_charges':[0,0],'all_VEV_insertions_preserve_Z6':True}},
  'split_projector':{**split_checks,'Tr_Y_squared':trY2,'g1_squared_over_gY_squared':g1_over_gY_sq,
                     'spectral_boundary_definition':'Lambda_split=5*gI*muI','MSigma_over_Lambda_split':R_SIGMA,
                     'split_uniqueness_up_to_orientation':split_uniqueness_pass,'pass':(split_pass and split_uniqueness_pass)},
  'active_corrected_full_threshold_closure':{
      'source':'results/utfanswf_embedding/result.json',
      'status':'ACTIVE_PINNED_INPUT_CLOSURE_DIAGNOSTIC',
      'MHC_GeV':active_full_threshold.get('MHC_GeV'),
      'MX_GeV':active_full_threshold.get('MX_GeV'),
      'MSigma_GeV':active_full_threshold.get('MSigma_GeV'),
      'alphaU_inverse':active_full_threshold.get('alphaU_inverse'),
      'closure_spread_inverse':active_full_threshold.get('closure_spread_inverse'),
      'blind_prediction':False,
  },
  'legacy_sigma_only_unification_diagnostic':{
      'status':'HISTORICAL_INCOMPLETE_THRESHOLD_DIAGNOSTIC_ONLY',
      'inputs':'pinned v22 low-energy couplings; diagnostic, not blind prediction',
      'MU_GeV':MU_exist,'MSigma_GeV':MSigma_exist,'inverse_couplings_after_matching':invU_exist.tolist(),
      'spread_inverse':exist_spread,'relative_spread':exist_rel,'relative_spread_percent':100*exist_rel,
      'omitted_thresholds':['X,Y gauge bosons','color-triplet Higgs H_C'],
      'retained_threshold':'Sigma_3 + Sigma_8 with weight (0,2,3)',
  },
  'boundary_sensitivity':{
      'if_MU_equals_raw_muI':{'MU_GeV':math.exp(res_raw.x),'spread_inverse':raw_spread,
                              'relative_spread_percent':100*raw_spread/float(np.mean(raw_vals))}
  },
  'blind_full_threshold_boundary_audit':{
      'status':'BLOCKED_BY_HEAVY_THRESHOLD_BOUNDARY',
      'active_blind_prediction_available':False,
      'alphaU_inverse_primitive_candidate':ALPHAU_INV_PRIMITIVE,
      'available_native_MSigma_over_MU_candidate':R_SIGMA,
      'missing_independent_boundaries':['M_HC/M_U (or L_HC)','M_X/M_U (or L_X)'],
      'reason':'After the corrected HC+XY+Sigma repair, alphaU plus alpha3 and R_Sigma do not uniquely fix all heavy-threshold data. Pinned closure masses cannot be imported into a blind electromagnetic solve because they use alpha1/alpha2.',
      'legacy_sigma_only_diagnostic':{
          'status':'HISTORICAL_INCOMPLETE_THRESHOLD_DIAGNOSTIC_ONLY',
          'alphaU_inverse':ALPHAU_INV_PRIMITIVE,'uses_only_alpha3_to_fix_scale':True,
          'MU_GeV':MU_blind,'MSigma_GeV':MS_blind,
          'alpha_MZ_prediction':a_blind.tolist(),'alpha_em_inverse_MZ_prediction':em_inv_blind,
          'withheld_v22_alpha_em_inverse_MZ':em_inv_hold,'relative_error':blind_rel,
          'absolute_percent_error':100*abs(blind_rel),'old_verdict':'CLEAN_NEAR_FAIL_1PCT_GATE',
          'incomplete_threshold_content':'Sigma_3 + Sigma_8 only; H_C and X/Y omitted',
      },
  },
  'historical_boundary':[
      'The dated public v22 manuscript did not instantiate the gauged Z6 or its charge-6 Abelian Higgsing field(s); the v23 executable candidate now does so without rewriting the dated v22 artifact.',
      'The conventional finite DeltaB=1 proton-lifetime branch is retained as historical/superseded rather than silently deleted.'
  ],
  'integration_status':[
      'The affected UTFANSWF suite has been rerun during v23 bundle hardening after exact-Z6 propagation; a final full release regression is still required after all remaining bundle edits.'
  ],
  'open_gates':[
      'The specific finite 3+2 projector-Dirac law remains a candidate construction, not a theorem derived uniquely from the dated public framework.',
      'The primitive alphaU^-1=8pi candidate is not source-complete after the corrected HC+XY+Sigma threshold repair because independent H_C/M_U and X/M_U boundary data are missing; the older ~1.15% Sigma-only near-fail is retained as history.',
      'DeltaB=3 Wilson coefficient/scale is not fixed, so no trinucleon lifetime prediction is claimed.',
      'Independent external reproduction/validation of the exact-Z6 proton-stability construction remains outstanding.'
  ],
  'references':['arXiv:hep-ph/0306003','arXiv:2204.01741','arXiv:1808.00009']
}
OUT.write_text(json.dumps(result,indent=2)+'\n')
print(json.dumps(result,indent=2))
