"""Independent, lightweight count-statistics review. No flight/source-depth validation."""
from pathlib import Path
import sys, math, json, csv, ast, hashlib, datetime, argparse, shutil, subprocess
HERE=Path(__file__).resolve().parent
PACKAGE=HERE.parent
if (PACKAGE/"runtime/deps").is_dir():
    sys.path.insert(0,str(PACKAGE/"runtime/deps"))
sys.path.insert(0,str(HERE))
from common import ROOT,H,C,PHOTON_ENERGY_337,PHOTON_ENERGY_777
import numpy as np
import scipy
from scipy.stats import poisson,binom,chi2,beta
from scipy.special import xlogy,gammaln
from scipy.optimize import brentq,minimize
T_PULSE=.003
RISES_US=(30,80,200)
ALPHA=2.866515718791933e-7
MARGIN=math.log(19.)

def _integer(n):
    if isinstance(n,(bool,np.bool_)) or int(n)!=n or n<0:raise ValueError("Nonnegative integer count required")
    return int(n)
def _finite(x,nonnegative=False):
    x=float(x)
    if not math.isfinite(x) or (nonnegative and x<0):raise ValueError("Finite nonnegative value required")
    return x
def _alpha(a):
    a=_finite(a)
    if not 0<a<1:raise ValueError("0<alpha<1 required")
    return a

def log_poisson_mass(k,mu):
    """Stable log-PMF; Stirling for k>=50. Not formal interval arithmetic."""
    k=_integer(k);mu=_finite(mu,True)
    if mu==0:return 0. if k==0 else -math.inf
    if k==0:return -mu
    if k<50:return -mu+k*math.log(mu)-math.lgamma(k+1)
    delta=(k-mu)/mu
    return (k-mu)-k*math.log1p(delta)-.5*math.log(2*math.pi*k)-1/(12*k)+1/(360*k**3)-1/(1260*k**5)

def direct_poisson_tail(k,mu,tol=2e-14):
    """P(N>=k), positive recurrence plus decreasing geometric remainder.
    For k<=mu sum the lower tail and subtract; k1 uses stable expm1.
    Includes no Gaussian approximation. PMF roundoff is separate from tail truncation.
    """
    k=_integer(k);mu=_finite(mu,True)
    if k==0:return 1.
    if mu==0:return 0.
    if k==1:return -math.expm1(-mu)
    if k-1<mu:
        j=k-1;term=math.exp(log_poisson_mass(j,mu));terms=[term];s=term
        while j>0:
            ratio=j/mu
            if ratio<1 and term*ratio/(1-ratio)<=tol*s:break
            term*=ratio;j-=1;terms.append(term);s+=term
        return max(0.,min(1.,1-math.fsum(terms)))
    j=k;term=math.exp(log_poisson_mass(j,mu));terms=[term];s=term
    while True:
        ratio=mu/(j+1)
        if term==0 or term*ratio/(1-ratio)<=tol*s:break
        j+=1;term*=mu/j;terms.append(term);s+=term
        if len(terms)>200000:raise RuntimeError("Use a qualified large-parameter method")
    return max(0.,min(1.,math.fsum(terms)))

def critical_known(b,alpha=ALPHA,method="scipy"):
    """Smallest integer k>=1 satisfying specified single-test tail."""
    b=_finite(b,True);alpha=_alpha(alpha)
    tail=(lambda k:float(poisson.sf(k-1,b))) if method=="scipy" else (lambda k:direct_poisson_tail(k,b))
    lo=0;hi=max(1,math.ceil(b+12*math.sqrt(b+1)+60))
    while tail(hi)>alpha:hi*=2
    while hi-lo>1:
        m=(hi+lo)//2
        if tail(m)<=alpha:hi=m
        else:lo=m
    return hi

def required_signal_known(b,alpha=ALPHA,power=.95):
    """Expected detected source count s; known b in same ON window."""
    b=_finite(b,True)
    if not 0<power<1:raise ValueError("0<power<1 required")
    k=critical_known(b,alpha)
    hi=max(b+10,k+10*math.sqrt(k+1))
    while poisson.sf(k-1,hi)<power:hi*=2
    total=b if poisson.sf(k-1,b)>=power else brentq(lambda m:float(poisson.sf(k-1,m))-power,b,hi,xtol=1e-11,rtol=1e-14)
    return {"critical":k,"source_mean_PE":total-b,"background_mean_PE":b,
            "power":float(poisson.sf(k-1,total)),"tail":float(poisson.sf(k-1,b))}

def onoff_pvalue(n,m,r):
    """N~Pois(s+b), M~Pois(r*b); r=T_off/T_on."""
    n=_integer(n);m=_integer(m);r=_finite(r)
    if r<=0:raise ValueError("r=T_off/T_on must be positive")
    return float(binom.sf(n-1,n+m,1/(1+r)))

def direct_onoff_tail(n,m,r):
    n=_integer(n);m=_integer(m);r=_finite(r)
    if r<=0:raise ValueError("Positive exposure ratio required")
    t=n+m;p=1/(1+r)
    if n==0:return 1.
    terms=[math.exp(math.lgamma(t+1)-math.lgamma(j+1)-math.lgamma(t-j+1)
           +j*math.log(p)+(t-j)*math.log1p(-p)) for j in range(n,t+1)]
    return math.fsum(terms)

def source_upper_known(n,b,alpha=.05):
    """At-least nominal unconditional coverage for expected source PE."""
    return max(0.,float(chi2.isf(_alpha(alpha),2*(_integer(n)+1))/2)-_finite(b,True))

def source_upper_off(n,m,r,alpha=.05,split=.5):
    """Conservative union-bound UL, fixed alpha split; r=T_off/T_on."""
    n=_integer(n);m=_integer(m);r=_finite(r);alpha=_alpha(alpha);split=_alpha(split)
    if r<=0:raise ValueError("Positive exposure ratio required")
    U=float(chi2.isf(alpha*split,2*(n+1))/2)
    L=0. if m==0 else float(chi2.ppf(alpha*(1-split),2*m)/2)
    return max(0.,U-L/r)

def mc_intervals(k,n,alpha=.05):
    """Individual Wilson and exact Clopper-Pearson binomial intervals."""
    k=_integer(k);n=_integer(n);alpha=_alpha(alpha)
    if k>n:raise ValueError("k<=n required")
    if n==0:return {"fraction":None,"Wilson":[None,None],"Clopper_Pearson":[None,None]}
    from scipy.stats import norm
    z=float(norm.isf(alpha/2));p=k/n;den=1+z*z/n
    mid=(p+z*z/(2*n))/den;half=z*math.sqrt(p*(1-p)/n+z*z/(4*n*n))/den
    cp=[0. if k==0 else float(beta.ppf(alpha/2,k,n-k+1)),
        1. if k==n else float(beta.isf(alpha/2,k+1,n-k))]
    return {"fraction":p,"Wilson":[max(0.,mid-half),min(1.,mid+half)],"Clopper_Pearson":cp}

def primitive(t,rise_us):
    z=np.clip(np.asarray(t,float),0,T_PULSE);d=rise_us*1e-6/.8
    return np.where(z<d,z*z/(2*d),z-d/2)

def q_integrals(left,right,onset,rise_us):
    return primitive(np.asarray(right)-onset,rise_us)-primitive(np.asarray(left)-onset,rise_us)

def poisson_loglik(n,mu):
    n=np.asarray(n,float);mu=np.asarray(mu,float)
    if np.any(mu<0):return -math.inf
    return float(np.sum(xlogy(n,mu)-mu-gammaln(n+1)))

def nuisance_profile(n,d,q):
    """Independent full-bin concave mixture solution, not imported R4 code."""
    n=np.asarray(n,float);d=np.asarray(d,float);q=np.asarray(q,float)
    if np.any(n<0) or np.any(d<=0) or np.any(q<0):raise ValueError("Invalid counts/windows")
    N=float(n.sum());D=float(d.sum());Q=float(q.sum())
    if Q<=0:raise ValueError("Nonzero source support required")
    if N==0:return {"theta":0.,"A_cps":0.,"b_cps":0.,"logL":0.,"mu_hat":np.zeros(len(n))}
    h=q/Q;w=d/D;a=h/w-1
    def derivative(x):
        use=n>0
        return float(np.sum(n[use]*a[use]/(1+x*a[use])))
    if derivative(0)<=0:theta=0.
    elif np.all(n[q==0]==0) and derivative(1-1e-14)>=0:theta=1.
    else:theta=brentq(derivative,0,1-1e-14,xtol=1e-14,rtol=1e-14)
    A=theta*N/Q;b=(1-theta)*N/D;mu=A*q+b*d
    return {"theta":theta,"A_cps":A,"b_cps":b,"logL":poisson_loglik(n,mu),"mu_hat":mu}

def worked_example(seed=202610221110):
    """Actual generated trial with explicit n,d,q and an independent 2D fit."""
    left=np.r_[np.arange(-100,400,20)*1e-6,.001,-.0105]
    right=np.r_[np.arange(-80,420,20)*1e-6,.0025,-.0005]
    d=right-left;truth={"rise_us":80,"onset_s":35e-6,"A_cps":1000/.003,"b_cps":1e4}
    q=q_integrals(left,right,truth["onset_s"],truth["rise_us"])
    q[-1]=0.
    rng=np.random.default_rng(seed);mu=truth["A_cps"]*q+truth["b_cps"]*d
    n=rng.poisson(mu);fit=nuisance_profile(n,d,q)
    # Direct original lambda=A*q+b*d optimized in two independent rate coordinates.
    scales=np.array([1e5,1e4])
    def obj(x):return -poisson_loglik(n,(x[0]*scales[0])*q+(x[1]*scales[1])*d)
    def jac(x):
        m=(x[0]*scales[0])*q+(x[1]*scales[1])*d
        ratio=np.divide(n,m,out=np.zeros_like(m),where=m>0)-1
        return -np.array([np.sum(ratio*q)*scales[0],np.sum(ratio*d)*scales[1]])
    opt=minimize(obj,[2.,1.],jac=jac,method="L-BFGS-B",bounds=[(1e-12,None)]*2,
                 options={"ftol":1e-15,"gtol":1e-9,"maxiter":300})
    direct_rates=opt.x*scales
    assert abs(fit["logL"]+opt.fun)<2e-7
    assert np.max(abs(direct_rates/np.array([fit["A_cps"],fit["b_cps"]])-1))<2e-5
    rows=[]
    for i in range(len(n)):
        rows.append({"bin":i,"label":"prewindow" if i==len(n)-1 else "late_plateau" if i==len(n)-2 else "early",
                     "left_s":float(left[i]),"right_s":float(right[i]),"d_s":float(d[i]),"q_s":float(q[i]),
                     "n_PE":int(n[i]),"truth_mean_PE":float(mu[i]),"fitted_mean_PE":float(fit["mu_hat"][i])})
    candidates=[]
    for rise in RISES_US:
      for onset in np.arange(-50,100.1,5)*1e-6:
        qq=q_integrals(left,right,onset,rise);qq[-1]=0
        f=nuisance_profile(n,d,qq)
        candidates.append({"rise_us":rise,"onset_us":float(onset*1e6),"A_cps":f["A_cps"],"b_cps":f["b_cps"],"logL":f["logL"]})
    candidates.sort(key=lambda x:x["logL"],reverse=True)
    result={"status":"EXECUTED_ONE_SYNTHETIC_TRIAL_NOT_EMPIRICAL_SOURCE_DATA","seed":seed,"truth":truth,
            "observed_total_PE":int(n.sum()),"window_total_s":float(d.sum()),"source_integral_q_s":float(q.sum()),
            "profile":{"theta":fit["theta"],"A_cps":fit["A_cps"],"b_cps":fit["b_cps"],"logL":fit["logL"]},
            "direct_SciPy_two_parameter_fit":{"success":bool(opt.success),"message":str(opt.message),"A_cps":float(direct_rates[0]),"b_cps":float(direct_rates[1]),"logL":float(-opt.fun)},
            "absolute_loglikelihood_difference":abs(fit["logL"]+opt.fun),
            "fitted_total_minus_observed_PE":float(fit["mu_hat"].sum()-n.sum()),
            "bins":rows,"top_grid_candidates":candidates[:9],
            "qualification":"The true80us candidate is used for the explicit2D/profile identity check. The separate3rise/onset search illustrates unknown morphology and onset; fitted labels are not source depth."}
    return result

def photon_planes(source_PE,efficiency=.2,diameter_mm=5):
    """Aperture fluence link; expected primary PE before multiplication."""
    eta=_finite(efficiency);s=_finite(source_PE,True);d=_finite(diameter_mm)
    if not 0<eta<=1 or d<=0:raise ValueError("Positive aperture/efficiency required")
    A=math.pi*(d*1e-3)**2/4
    return {"expected_source_primary_PE":s,"assumed_efficiency":eta,"area_m2":A,
            "expected_aperture_photons":s/eta,"aperture_photon_fluence_m2":s/(eta*A),
            "aperture_energy_fluence_J_m2_at337":s/(eta*A)*PHOTON_ENERGY_337,
            "E337_J":PHOTON_ENERGY_337,"E777_J":PHOTON_ENERGY_777,
            "qualification":"Aperture photons after propagation; no emitted-source yield or second cloud attenuation; no PMT/APD gain factor."}

def count_checks():
    rows=[]
    for mu in [0.,1e-12,.01,.1,1.,10.,100.,1000.,1e5,1e6]:
      ns=sorted({0,1,max(0,int(mu)),int(mu+5*math.sqrt(mu+1))})
      for n in ns:
        a=direct_poisson_tail(n,mu);b=float(poisson.sf(n-1,mu))
        assert abs(a-b)<1e-8
        rows.append({"kind":"Poisson_tail","mean":mu,"n_inclusive":n,"direct_sum":a,"SciPy":b,
                     "abs_difference":abs(a-b),"relative_difference":abs(a/b-1) if b else None})
    thresholds=[]
    for b in [0.,.01,1.,10.,100.,1000.,1e5,1e6]:
      for a in [ALPHA,.1/(50000*86400*3)]:
        t=required_signal_known(b,a);k2=critical_known(b,a,method="sum")
        assert k2==t["critical"]
        direct_power=direct_poisson_tail(k2,b+t["source_mean_PE"])
        assert abs(direct_power-.95)<1e-8
        thresholds.append({"background_PE":b,"alpha_per_test":a,**t,
                           "critical_direct_sum":k2,"power_direct_sum":direct_power})
    onoff=[]
    for n,m,r in [(0,0,10),(1,0,10),(10,0,10),(5,2,3),(30,40,10),(100,200,10)]:
        a=onoff_pvalue(n,m,r);b=direct_onoff_tail(n,m,r);assert abs(a-b)<1e-12
        onoff.append({"n_on":n,"m_off":m,"r_off_over_on":r,"SciPy":a,"finite_sum":b})
    assert abs(source_upper_known(0,0)+math.log(.05))<1e-13
    assert abs(source_upper_off(0,0,10)+math.log(.025))<1e-13
    upper=[]
    for n,b in [(0,0),(0,2),(1,0),(3,2),(10,5),(100,80)]:
        u=source_upper_known(n,b)
        total=float(chi2.isf(.05,2*(n+1))/2)
        cdf_direct=1-direct_poisson_tail(n+1,total)
        assert abs(cdf_direct-.05)<1e-11
        upper.append({"n_on":n,"known_b_PE":b,"source95_upper_PE":u,"total_mean_inversion_CDF":cdf_direct})
    # Coverage is an independent finite enumeration, not a Monte Carlo claim.
    coverage=[]
    for src,b,r in [(0,0,10),(1,0,10),(3,2,10),(10,3,5)]:
        n=np.arange(100);m=np.arange(200)
        pn=poisson.pmf(n,src+b);pm=poisson.pmf(m,r*b)
        UK=np.maximum(0.,chi2.isf(.05,2*(n+1))/2-b)
        up_total=chi2.isf(.025,2*(n+1))/2
        lower=np.zeros(len(m));lower[1:]=chi2.ppf(.025,2*m[1:])/2
        UF=np.maximum(0.,up_total[:,None]-lower[None,:]/r)
        pk=float(pn[UK>=src-1e-12].sum())
        pf=float(np.sum((UF>=src-1e-12)*pn[:,None]*pm[None,:]))
        omitted=max(0.,1-float(pn.sum()*pm.sum()))
        assert pk>=.95-1e-10 and pf>=.95-1e-10
        coverage.append({"source_PE":src,"b_on_PE":b,"r":r,"known_b_enumerated_coverage":pk,
                         "finite_off_enumerated_coverage":pf,"omitted_probability_bound":omitted})
    pins=[]
    for mu,k,pinned in [(1000000.,1005000,2.934034048031641e-7),
                        (1000000.,1007000,1.3595737503819116e-12),
                        (1000000.,1005004,2.8739183629417515e-7),
                        (1000000.,1005005,2.859076013349398e-7)]:
        val=direct_poisson_tail(k,mu);rel=abs(val/pinned-1)
        assert rel<1e-8
        pins.append({"mean":mu,"n_inclusive":k,"pinned_60digit_Decimal_tail":pinned,
                     "direct_recurrence":val,"relative_difference":rel,
                     "SciPy_tail":float(poisson.sf(k-1,mu)),
                     "provenance":"R4 qa/browser_large_mean_tail_control.json;60digit Decimal/Stirling/geometric-remainder calculation, not empirical data."})
    return {"tail_checks":rows,"pinned_high_mean_controls":pins,"threshold_power_checks":thresholds,"onoff_checks":onoff,
            "upper_inversions":upper,"coverage_finite_enumerations":coverage,
            "zero_count_controls":{"known_zero_background95":source_upper_known(0,0),"finite_baseline_zero_zero95":source_upper_off(0,0,10)},
            "qualification":"Exact probability constructions, floating-point implementations. At mean1e6 the SciPy incomplete-gamma tail differs from direct positive summation by several parts per million near5sigma; integer checks here agree. This is not formal interval arithmetic."}

def confusion_checks():
    report=json.loads((ROOT/"data/statistics_results.json").read_text(encoding="utf-8-sig"))
    pr=report["practical_morphology"];rows=[]
    metric_count=0;max_ci_difference=0.
    for case in pr["scenarios"]:
        trial=ROOT/"data/trials"/f"case_{case['case_id']:03d}.npz"
        if not trial.is_file():raise FileNotFoundError("Missing synthetic trial summaries: "+str(trial))
        with np.load(trial,allow_pickle=False) as t:
            scores=t["scores"];truth=RISES_US.index(case["truth_receiver_rise_10_90_us"])
            best=scores.argmax(axis=1);group=(scores[:,2]>scores[:,:2].max(axis=1)).astype(int)
            detected=np.asarray(t["detection_pvalue"])<=ALPHA
            pcheck=binom.sf(t["gate_count"]-1,t["gate_count"]+t["pre_count"],.0034/(.0034+.010))
            assert np.allclose(pcheck,t["detection_pvalue"],rtol=1e-8,atol=1e-14)
            sort=np.sort(scores,axis=1);margin3=sort[:,-1]-sort[:,-2]
            margin2=np.abs(scores[:,2]-scores[:,:2].max(axis=1))
            a3=detected&(margin3>=MARGIN);a2=detected&(margin2>=MARGIN)
            pred3=np.where(~detected,4,np.where(~a3,3,best))
            pred2=np.where(~detected,3,np.where(~a2,2,group))
            c3=np.bincount(pred3,minlength=5).tolist();c2=np.bincount(pred2,minlength=4).tolist()
            assert c3==case["three_class"]["confusion_counts"] and c2==case["fast_slow"]["confusion_counts"]
            assert len(best)==case["trials"]==1000
            metrics=[(case["detection"],int(detected.sum()),len(best)),
                     (case["three_class"]["accepted"],int(a3.sum()),len(best)),
                     (case["three_class"]["wrong_among_accepted"],int((a3&(best!=truth)).sum()),int(a3.sum())),
                     (case["fast_slow"]["accepted"],int(a2.sum()),len(best)),
                     (case["fast_slow"]["wrong_among_accepted"],int((a2&(group!=(truth==2))).sum()),int(a2.sum())),
                     (case["three_class"]["rejected_among_detected"],int((detected&~a3).sum()),int(detected.sum())),
                     (case["fast_slow"]["rejected_among_detected"],int((detected&~a2).sum()),int(detected.sum())),
                     (case["unthresholded_practical_three_class_error"],int((best!=truth).sum()),len(best))]
            rebuilt=[]
            for old,k,n in metrics:
                assert old["count"]==k and old["denominator"]==n
                ii=mc_intervals(k,n);metric_count+=1
                if n:
                    for field,newkey in [("individual_95pct_Wilson_interval","Wilson"),
                                         ("individual_95pct_Clopper_Pearson_interval","Clopper_Pearson")]:
                        if field in old:
                            diff=float(np.max(abs(np.array(old[field])-ii[newkey])))
                            max_ci_difference=max(max_ci_difference,diff);assert diff<2e-12
                rebuilt.append({"k":k,"n":n,**ii})
            rows.append({"case_id":case["case_id"],"N_plateau_equivalent_PE":case["N_plateau_equivalent_3ms_PE"],
                         "background_cps":case["background_cps"],"sample_khz":case["sample_rate_khz"],
                         "truth_rise_us":case["truth_receiver_rise_10_90_us"],"confusion3":c3,"confusion_fastslow":c2,
                         "metrics_order":["detected","accepted3","wrong_given_accepted3","accepted_fastslow","wrong_given_accepted_fastslow","rejected3_given_detected","rejected_fastslow_given_detected","unthresholded3_error"],
                         "rebuilt_metrics":rebuilt,"trial_sha256":hashlib.sha256(trial.read_bytes()).hexdigest()})
    assert len(rows)==81
    return {"status":"ALL_SAVED_TRIAL_SUMMARIES_RECOMPUTED","scenarios":81,"signal_trials":81000,
            "metrics_recomputed":metric_count,"max_interval_absolute_difference":max_ci_difference,"rows":rows,
            "qualification":"Rebuilds decisions and denominators from saved synthetic fitted scores, not raw experimental traces. Does not independently refit every pulse because full per-bin count arrays were not retained in these summaries. No depth-label or flight validation."}

def original_source_map():
    requested={"oracle_reference/calculate_ideal_poisson_rise.py":["experiment","affinity_bounds"],"practical_morphology.py":["primitive","fit_one","draw_case","analyze_case"],
               "count_inference/count_inference.py":["count_threshold_known_background","exact_onoff_pvalue",
                  "poisson_total_mean_upper","poisson_total_mean_lower","red_source_upper_known_background",
                  "red_source_upper_finite_baseline"]}
    result=[]
    base=ROOT/"reference/statistics"
    for rel,names in requested.items():
        p=base/rel
        if not p.is_file():result.append({"file":str(p.relative_to(ROOT)),"present":False});continue
        tree=ast.parse(p.read_text(encoding="utf-8-sig"))
        line={n.name:n.lineno for n in ast.walk(tree) if isinstance(n,(ast.FunctionDef,ast.AsyncFunctionDef))}
        result.append({"file":str(p.relative_to(ROOT)),"present":True,"sha256":hashlib.sha256(p.read_bytes()).hexdigest(),
                       "functions":{n:line.get(n) for n in names}})
    return result

def javascript_crosscheck(node_path=None):
    p=ROOT/"reference/statistics/browser_math.js"
    if not p.is_file():return {"status":"NOT_RUN","reason":"Reference browser_math.js not present"}
    lines=p.read_text(encoding="utf-8-sig").splitlines()
    mapping={name:next((i+1 for i,line in enumerate(lines) if "function "+name+"(" in line),None)
             for name in ["sf","critical","threshold","upper"]}
    node=node_path or shutil.which("node")
    receipt={"file":str(p.relative_to(ROOT)),"sha256":hashlib.sha256(p.read_bytes()).hexdigest(),"function_lines":mapping}
    if not node:return {**receipt,"status":"NOT_RUN","reason":"Node optional and unavailable; Python checks do not execute JavaScript"}
    cases=[]
    for b in [0.,.01,1.,10.,1000.,1e6]:
      for alpha in [ALPHA,.1/(50000*86400*3)]:cases.append([b,alpha])
    red_cases=[[0,0],[0,2],[1,0],[3,2],[10,5],[100,80]]
    tiny_cases=[]
    for alpha in [ALPHA,.1/(50000*86400*3)]:
        for factor in [1-1e-8,1+1e-8]:
            tiny_cases.append([-math.log1p(-alpha)*factor,alpha])
    script="const M=require(process.argv[1]);const q=JSON.parse(process.argv[2]);process.stdout.write(JSON.stringify({threshold:q.threshold.map(a=>M.threshold(...a)),upper:q.upper.map(a=>M.upper(...a)),tiny:q.tiny.map(a=>({sf1:M.sf(1,a[0]),critical:M.critical(...a)}))}));"
    payload={"threshold":cases,"upper":red_cases,"tiny":tiny_cases}
    run=subprocess.run([node,"-e",script,str(p),json.dumps(payload)],text=True,capture_output=True,check=True)
    outputs=json.loads(run.stdout);checks=[]
    for (b,a),js in zip(cases,outputs["threshold"]):
        ref=required_signal_known(b,a)
        assert js["critical"]==ref["critical"] and abs(js["signal"]-ref["source_mean_PE"])<5e-5
        checks.append({"b_PE":b,"alpha":a,"JS":js,"Python":ref,"signal_abs_difference":abs(js["signal"]-ref["source_mean_PE"])})
    red=[]
    for (n,b),js in zip(red_cases,outputs["upper"]):
        ref=source_upper_known(n,b);assert abs(js-ref)<1e-8
        red.append({"n_on":n,"b_on_PE":b,"JS_upper95":js,"Python_upper95":ref,"abs_difference":abs(js-ref)})
    tiny=[]
    for (b,a),js in zip(tiny_cases,outputs["tiny"]):
        ref=-math.expm1(-b);critical=critical_known(b,a,method="sum")
        assert js["critical"]==critical and abs(js["sf1"]/ref-1)<1e-12
        tiny.append({"b_PE":b,"alpha":a,"JS":js,"expected_sf1":ref,"expected_critical":critical})
    version=subprocess.run([node,"--version"],text=True,capture_output=True,check=True).stdout.strip()
    return {**receipt,"status":"EXECUTED_OPTIONAL_NODE","node_version":version,"checks":checks,
            "red_upper_checks":red,"tiny_background_boundary_checks":tiny}

def profile_boundary_checks():
    rows=[]
    for n,theta in [([0,0],0.),([10,0],0.),([0,10],1.),([10,20],1/3)]:
        fit=nuisance_profile(n,[1.,1.],[0.,1.])
        assert abs(fit["theta"]-theta)<1e-12
        assert abs(float(fit["mu_hat"].sum())-sum(n))<1e-12
        rows.append({"counts":n,"expected_theta":theta,"fitted_theta":fit["theta"],
                     "A_cps":fit["A_cps"],"b_cps":fit["b_cps"],"logL":fit["logL"]})
    return rows

def affinity_bound_values(mu,nu):
    mu=np.asarray(mu,float);nu=np.asarray(nu,float)
    if mu.shape!=nu.shape or np.any(mu<0) or np.any(nu<0):raise ValueError("Matched nonnegative bin means required")
    h=float(.5*np.sum((np.sqrt(mu)-np.sqrt(nu))**2))
    bc=math.exp(-h)
    root=math.sqrt(max(0.,1-bc*bc))
    return {"affinity":bc,"Bayes_error_lower":bc*bc/(2*(1+root)),"Bayes_error_upper":bc/2}

def oracle_affinity_controls():
    rows=[]
    for mu,nu in [([.2,1.],[1.3,2.5]),([0.,.3],[.2,1.]),([1.,.1],[1.,3.])]:
        pa=[np.array([math.exp(log_poisson_mass(k,m)) for k in range(60)]) for m in mu]
        qa=[np.array([math.exp(log_poisson_mass(k,m)) for k in range(60)]) for m in nu]
        p=np.outer(pa[0],pa[1]);q=np.outer(qa[0],qa[1])
        bc=float(np.sqrt(p*q).sum());pe=float(.5*np.minimum(p,q).sum())
        formula=affinity_bound_values(mu,nu)
        assert abs(p.sum()-1)<1e-13 and abs(q.sum()-1)<1e-13
        assert abs(bc-formula["affinity"])<1e-13
        assert formula["Bayes_error_lower"]-1e-13<=pe<=formula["Bayes_error_upper"]+1e-13
        rows.append({"mu_PE_per_bin":mu,"nu_PE_per_bin":nu,"explicit_joint_PMF_affinity":bc,
                     "explicit_equal_prior_Bayes_error":pe,**formula,"finite_grid_mass_P":float(p.sum()),"finite_grid_mass_Q":float(q.sum())})
    return {"status":"THREE_NEW_TWO_BIN_CONTROLS_EXECUTED","rows":rows,
            "retained320_grid":"NOT_INDEPENDENTLY_REBUILT_IN_THIS_LIGHT_REVIEW",
            "assumptions":"Known simple hypotheses, independent Poisson bins, equal priors and0-1loss. Unknown nuisance fits or rejected/accepted-only errors do not inherit these bounds."}

def _numeric_assignments(path):
    """Read simple numeric constants with AST; never execute retained scripts."""
    env={}
    def calc(n):
        if isinstance(n,ast.Constant) and isinstance(n.value,(int,float)):return n.value
        if isinstance(n,ast.Name):return env[n.id]
        if isinstance(n,ast.BinOp):
            x=calc(n.left);y=calc(n.right)
            if isinstance(n.op,ast.Mult):return x*y
            if isinstance(n.op,ast.Div):return x/y
            if isinstance(n.op,ast.Add):return x+y
            if isinstance(n.op,ast.Sub):return x-y
        raise ValueError("Not a simple numeric expression")
    for node in ast.parse(path.read_text(encoding="utf-8-sig")).body:
        if isinstance(node,ast.Assign) and len(node.targets)==1 and isinstance(node.targets[0],ast.Name):
            try:env[node.targets[0].id]=calc(node.value)
            except (ValueError,KeyError):pass
    return env

def downstream_threshold_bindings():
    data=json.loads((ROOT/"data/analytic_results.json").read_text(encoding="utf-8-sig"))
    ref=ROOT/"reference/analytic/count_requirements.py"
    fovref=ROOT/"reference/analytic/fov_population.py"
    constants=_numeric_assignments(ref)
    operational=constants["ALPHA_OPERATIONAL"]
    assert abs(operational-.1/(50000*86400*3))<1e-25
    source=fovref.read_text(encoding="utf-8-sig")
    assert "m.count_at_power((rate+dark)*dt,m.ALPHA_OPERATIONAL,.95)" in source
    dark=float(data["fov_population"]["report"]["inputs"]["dark_cps"])
    cache={};rows=[]
    for study in ["fov_population","auroral_counts"]:
      for index,row in enumerate(data[study]["rows"]):
        if study=="fov_population":
            b=(row["background_photoelectron_cps"]+dark)*row["window_s"]
            alpha=operational;storedS=row["required_source_count_for_95pct_detection"]
            alpha_source="reference/analytic/fov_population.py:41 -> count_requirements.py ALPHA_OPERATIONAL"
        else:
            b=row["known_background_mean_count"];alpha=row["per_test_alpha"]
            storedS=row["signal_mean_for_95pct_power"];alpha_source="row.per_test_alpha"
        key=(b,alpha)
        if key not in cache:cache[key]=required_signal_known(b,alpha)
        check=cache[key]
        delta=abs(storedS-check["source_mean_PE"])
        assert row["critical_observed_count"]==check["critical"]
        assert delta<=5e-5+1e-9*abs(storedS)
        achieved=float(poisson.sf(check["critical"]-1,b+storedS))
        assert abs(achieved-.95)<1e-8
        rows.append({"study":study,"row_index":index,"background_mean_PE":b,
                     "alpha_per_test":alpha,"alpha_source":alpha_source,
                     "stored_critical":row["critical_observed_count"],"recomputed_critical":check["critical"],
                     "stored_S95_PE":storedS,"recomputed_S95_PE":check["source_mean_PE"],
                     "S95_abs_difference_PE":delta,"power_at_stored_S95":achieved})
    assert sum(x["study"]=="fov_population" for x in rows)==216
    assert sum(x["study"]=="auroral_counts" for x in rows)==192
    return {"status":"ALL408_STATISTICAL_COLUMN_BINDINGS_CHECKED",
            "rows_checked":len(rows),"unique_background_alpha_computations":len(cache),
            "FOV_rows":216,"auroral_rows":192,
            "maximum_S95_absolute_difference_PE":max(x["S95_abs_difference_PE"] for x in rows),
            "maximum_power_absolute_difference":max(abs(x["power_at_stored_S95"]-.95) for x in rows),
            "constant_provenance":{"file":str(ref.relative_to(ROOT)),"sha256":hashlib.sha256(ref.read_bytes()).hexdigest(),
                                   "ALPHA_OPERATIONAL":operational,"dark_cps_from_FOV_report":dark},
            "qualification":"Checks only known-background Poisson critical/S95columns and their binding to actualB/declared false-alarm alpha. Population-slope alpha is a different variable. Does not validate FOVsource geography, brightness law, background scene or mission rate.",
            "rows":rows}



def equation_function_map():
    mapping={
      "S01":["photon_planes"],"S02":["photon_planes"],
      "S03":["onoff_pvalue"],"S04":["onoff_pvalue"],"S05":["critical_known"],
      "S06":["required_signal_known"],"S07":["onoff_pvalue","direct_onoff_tail"],
      "S08":["count_checks"],"S09":["source_upper_known"],"S10":["source_upper_known"],
      "S11":["source_upper_off"],"S12":["source_upper_off"],
      "S13":["q_integrals","poisson_loglik"],"S14":["nuisance_profile"],
      "S15":["nuisance_profile"],"S16":["nuisance_profile","worked_example"],
      "S17":["confusion_checks"],"S18":["mc_intervals"],"S19":["mc_intervals"],"S20":["affinity_bound_values"],"S21":["oracle_affinity_controls"]}
    tree=ast.parse(Path(__file__).read_text(encoding="utf-8-sig"))
    lines={n.name:n.lineno for n in ast.walk(tree) if isinstance(n,ast.FunctionDef)}
    return {eq:[{"file":"modules/statistics_review.py","function":f,"line":lines[f]} for f in funs]
            for eq,funs in mapping.items()}

def run_checks(write=False,node_path=None):
    checks={"schema":"INDEPENDENT_STATISTICS_REVIEW_V1","executed_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),
            "versions":{"python":sys.version.split()[0],"numpy":np.__version__,"scipy":scipy.__version__},
            "photon_planes":photon_planes(100),
            "counts":count_checks(),"oracle_affinity":oracle_affinity_controls(),"downstream_threshold_bindings":downstream_threshold_bindings(),"profile_boundaries":profile_boundary_checks(),"worked_example":worked_example(),"confusions":confusion_checks(),
            "source_map":original_source_map(),"equation_function_map":equation_function_map(),"javascript":javascript_crosscheck(node_path),
            "theory_scope":"Own derivations are in methods/statistics.md. Independent numerical comparisons check implementation, not physical assumptions or source classes.",
            "mission_rate":"NOT_COMPUTABLE_FROM_THESE_INPUTS","passed":True}
    if write:
        out=ROOT/"checks";out.mkdir(exist_ok=True)
        (out/"statistics_results.json").write_text(json.dumps(checks,indent=2,allow_nan=False),encoding="utf-8")
        rows=checks["worked_example"]["bins"]
        with (out/"statistics_worked_bins.csv").open("w",newline="",encoding="utf-8") as f:
            w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
        (out/"statistics_summary.json").write_text(json.dumps({
            "passed":True,"worked_logL_abs_difference":checks["worked_example"]["absolute_loglikelihood_difference"],
            "scenarios_rebuilt":checks["confusions"]["scenarios"],"metrics_rebuilt":checks["confusions"]["metrics_recomputed"],
            "max_interval_abs_difference":checks["confusions"]["max_interval_absolute_difference"],
            "zero_count_controls":checks["counts"]["zero_count_controls"],"javascript_status":checks["javascript"]["status"]},indent=2),encoding="utf-8")
    return checks

if __name__=="__main__":
    arg=argparse.ArgumentParser();arg.add_argument("--write",action="store_true");arg.add_argument("--node")
    a=arg.parse_args();result=run_checks(a.write,a.node)
    print(json.dumps({"passed":result["passed"],"worked_logL_difference":result["worked_example"]["absolute_loglikelihood_difference"],
                     "scenarios":result["confusions"]["scenarios"],"zero_counts":result["counts"]["zero_count_controls"],
                     "JavaScript":result["javascript"]["status"]},indent=2))