"""Exact Poisson planning requirements; no measured receiver is asserted.

Only Python's standard library and numpy are required. No normal approximation
is used in the threshold/95%-power or one-sided count interval calculations.
"""
from pathlib import Path
import math, json, csv
from functools import lru_cache
import numpy as np

HERE=Path(__file__).resolve().parent
ROOT=HERE.parents[2]
EPH=6.62607015e-34*299792458/337.1e-9
ALPHA_SINGLE=2.866515718791933e-7
TESTS_PER_DAY=50000*86400*3
ALPHA_OPERATIONAL=.1/TESTS_PER_DAY

def poisson_sf(n,mu):
    """P(N>=n). Positive recurrences in the smaller tail avoid cancellation."""
    n=int(n)
    if n<=0:return 1.
    if mu==0:return 0.
    if mu<0:raise ValueError('negative Poisson mean')
    if n-1<mu:
        k=n-1;term=math.exp(-mu+k*math.log(mu)-math.lgamma(k+1));s=term
        while k>0:
            term*=k/mu;k-=1;s+=term
            if term<1e-15*s:break
        return max(0.,min(1.,1-s))
    k=n;term=math.exp(-mu+k*math.log(mu)-math.lgamma(k+1));s=term
    if term==0:return 0.
    for _ in range(1000000):
        k+=1;term*=mu/k;s+=term
        if term<1e-15*s:break
    else:raise RuntimeError('Poisson recurrence did not converge')
    return max(0.,min(1.,s))

def poisson_cdf(n,mu):
    if n<0:return 0.
    if mu==0:return 1.
    n=int(n)
    if n<mu:
        k=n;term=math.exp(-mu+k*math.log(mu)-math.lgamma(k+1));s=term
        while k>0:
            term*=k/mu;k-=1;s+=term
            if term<1e-15*s:break
        return max(0.,min(1.,s))
    return 1-poisson_sf(n+1,mu)

@lru_cache(None)
def critical_count(bg,alpha):
    low=0;high=max(1,int(bg+12*math.sqrt(bg+1)+100))
    while poisson_sf(high,bg)>alpha:high*=2
    while high-low>1:
        mid=(high+low)//2
        if poisson_sf(mid,bg)<=alpha:high=mid
        else:low=mid
    return high

@lru_cache(None)
def count_at_power(bg,alpha,power=.95):
    n=critical_count(bg,alpha);lo=0.;hi=max(10.,n-bg+8*math.sqrt(n+1)+100)
    while poisson_sf(n,bg+hi)<power:hi*=2
    for _ in range(55):
        mid=(lo+hi)/2
        if poisson_sf(n,bg+mid)>=power:hi=mid
        else:lo=mid
    return n,hi

@lru_cache(None)
def poisson_upper(n,alpha):
    lo=0.;hi=max(10.,n+12*math.sqrt(n+1)+100)
    while poisson_cdf(n,hi)>alpha:hi*=2
    for _ in range(55):
        mid=(lo+hi)/2
        if poisson_cdf(n,mid)<alpha:hi=mid
        else:lo=mid
    return hi

@lru_cache(None)
def poisson_lower(n,alpha):
    if n==0:return 0.
    lo=0.;hi=float(n)
    for _ in range(55):
        mid=(lo+hi)/2
        if poisson_sf(n,mid)>alpha:hi=mid
        else:lo=mid
    return lo

def conservative_ratio_upper(nuv,nred,muv,mred,on_s,off_s,kappa_lower,alpha=.05):
    """Four marginal one-sided intervals; union bound gives >=1-alpha coverage.

    Coverage is unconditional for a prespecified window. A UV-trigger-selected
    catalog needs a separate selection-aware coverage argument.
    kappa_lower is a deterministic validated lower bound, not a nominal estimate
    or a confidence interval silently used without an additional error budget.
    """
    a=alpha/4;scale=on_s/off_s
    sr=max(0.,poisson_upper(nred,a)-scale*poisson_lower(mred,a))
    su=max(0.,poisson_lower(nuv,a)-scale*poisson_upper(muv,a))
    return (None if su==0 else sr/su/kappa_lower),sr,su

def write(name,rows):
    with (HERE/name).open('w',newline='',encoding='utf-8')as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)

def main():
    thresholds=[]
    for diameter in [1.,5.,25.]:
      area=math.pi*(diameter*.001)**2/4
      for fov in [60.,90.]:
       omega=math.pi*math.sin(math.radians(fov/2))**2
       quiet=.2*area*omega*4*3e-9/EPH
       for aurora_kr in [0.,1.]:
        aurora=.2*area*omega*aurora_kr*1000*1e10/(4*math.pi)
        for dark in [50.,1e4,1e6]:
         rate=quiet+aurora+dark
         for window in [30e-6,3e-3]:
          bg=rate*window
          for name,alpha in [('single_5sigma_tail',ALPHA_SINGLE),('0p1_false_per_observing_day',ALPHA_OPERATIONAL)]:
           n,s=count_at_power(bg,alpha)
           pfa=poisson_sf(n,bg)
           assert pfa<=alpha*(1+1e-9) and poisson_sf(n-1,bg)>alpha
           assert abs(poisson_sf(n,bg+s)-.95)<2e-10
           thresholds.append(dict(diameter_mm=diameter,full_fov_deg=fov,efficiency=.2,bandwidth_nm=4.,quiet_proxy_W_m2_sr_nm=3e-9,aurora_kR=aurora_kr,dark_cps=dark,window_s=window,background_rate_cps=rate,known_background_mean=bg,criterion=name,per_test_alpha=alpha,critical_observed_count=n,actual_poisson_false_alarm=pfa,expected_signal_for_95pct_detection=s,expected_aperture_signal_photons=s/.2,required_aperture_fluence_ph_m2=s/(.2*area),required_aperture_energy_J_m2=s*EPH/(.2*area),assumed_deadtime_s=20e-9,background_count_loss_nonparalyzable=rate*20e-9/(1+rate*20e-9),background_known_exactly=True))
    write('exact_detection_requirements.csv',thresholds)
    ratio=[]
    for kappa in [.1,.6223529411764706,1.,2.]:
     for confidence in [.95,.99]:
      for n in [10,30,100,300,1000,10000]:
       a=1-confidence;upper=math.expm1(-math.log(a)/n)/kappa
       ratio.append(dict(uv_observed_count=n,red_observed_count=0,kappa=kappa,confidence=confidence,exact_conditional_aperture_red_uv_upper=upper,background_free=True))
    write('zero_red_exact_ratio_bounds.csv',ratio)
    count_targets=[]
    for kappa in [.6223529411764706,1.]:
     for r in [.01,.1,1.]:
      n=math.floor(-math.log(.05)/math.log1p(kappa*r))+1
      count_targets.append(dict(kappa=kappa,thermal_aperture_ratio_tested=r,confidence=.95,minimum_uv_count_with_zero_red_for_upper_strictly_below_ratio=n,background_free=True))
    write('null_channel_count_requirements.csv',count_targets)
    onoff=[]
    for kappa in [.6223529411764706,1.]:
     for s in [10,30,100,300,1000,10000]:
      for bg in [0.,1.,10.,100.,1000.]:
       nuv=int(round(s+bg));nred=int(round(bg));muv=mred=int(round(3*bg))
       ru,sru,sul=conservative_ratio_upper(nuv,nred,muv,mred,.003,.009,kappa)
       onoff.append(dict(representative_uv_signal_count=s,background_per_channel_on=bg,kappa_validated_lower=kappa,confidence_lower_bound=.95,nuv=nuv,nred=nred,muv=muv,mred=mred,on_s=.003,off_s=.009,red_signal_upper=sru,uv_signal_lower=sul,aperture_red_uv_upper=ru,upper_is_finite=ru is not None,representative_counts_not_detection_probability=True))
    write('background_aware_ratio_bounds.csv',onoff)
    # A directly implementable guard for an estimated, stationary baseline.
    # This trades power for a proof by two one-sided probability bounds.
    guards=[]
    for row in thresholds:
     if row['diameter_mm']==5 and row['full_fov_deg']==60 and row['dark_cps']==50 and row['criterion']=='0p1_false_per_observing_day':
      off=.009;on=row['window_s'];m=int(round(row['background_rate_cps']*off));alpha=row['per_test_alpha']
      bg_upper=on/off*poisson_upper(m,alpha/2);crit=critical_count(bg_upper,alpha/2)
      guards.append(dict(aurora_kR=row['aurora_kR'],on_s=on,off_s=off,representative_off_count=m,per_test_error_allowance=alpha,background_upper_for_this_off_count=bg_upper,critical_on_count=crit,known_background_critical_count=row['critical_observed_count'],single_test_false_alarm_bound=alpha,power_not_integrated_over_random_off_count=True))
    write('estimated_baseline_trigger_guard.csv',guards)
    # Validate Poisson recurrence against a separately evaluated 50-digit series.
    from decimal import Decimal,localcontext
    check=[]
    with localcontext()as context:
     context.prec=50
     for mu in [.001,.15,1.,5.,20.,100.]:
      for n in [1,2,5,10,25,100,150]:
       m=Decimal(str(mu));term=(-m).exp();cdf=Decimal(0)
       for k in range(n):
        if k:term*=m/Decimal(k)
        cdf+=term
       sf=float(1-cdf);actual=poisson_sf(n,mu)
       assert abs(sf-actual)<=max(2e-13,abs(sf)*1e-10),(mu,n,sf,actual)
       check.append(dict(mean=mu,n=n,decimal_sf=sf,recurrence_sf=actual))
    assert abs(poisson_upper(0,.05)+math.log(.05))<1e-12
    assert abs(poisson_upper(1,.05)-4.743864518390577)<1e-10
    # Coverage audit for the conservative ratio bound. This is an implementation
    # check; the union-bound argument supplies coverage under the stated model.
    rng=np.random.default_rng(2026102201);coverage=[];coverage_raw=[]
    for su_true in [30.,300.]:
     for bg in [0.,10.,100.]:
      for r_true in [.01,.1,1.]:
       kappa=.6223529411764706;rejected=0;finite=0;ntrial=1024
       draws=np.column_stack([rng.poisson(su_true+bg,ntrial),rng.poisson(kappa*r_true*su_true+bg,ntrial),rng.poisson(3*bg,ntrial),rng.poisson(3*bg,ntrial)])
       for trial,(nu,nr,mu,mr) in enumerate(draws):
        ru,_,_=conservative_ratio_upper(int(nu),int(nr),int(mu),int(mr),.003,.009,kappa)
        finite+=ru is not None;rejected+=ru is not None and ru<r_true
        coverage_raw.append(dict(true_uv_signal=su_true,true_aperture_red_uv_ratio=r_true,background_mean_per_channel=bg,trial=trial,nuv=int(nu),nred=int(nr),muv=int(mu),mred=int(mr),kappa_lower=kappa,aperture_ratio_upper=ru,finite_upper=ru is not None,false_exclusion=ru is not None and ru<r_true))
       coverage.append(dict(true_uv_signal=su_true,true_aperture_red_uv_ratio=r_true,background_mean_per_channel=bg,trials=ntrial,false_exclusion_count=rejected,finite_upper_count=finite,false_exclusion_fraction=rejected/ntrial,nominal_familywise_error_bound=.05))
    write('ratio_coverage_audit.csv',coverage)
    write('ratio_coverage_trials.csv',coverage_raw)
    import gzip
    with gzip.open(HERE/'ratio_coverage_trials.csv.gz','wb')as f:f.write((HERE/'ratio_coverage_trials.csv').read_bytes())
    shields=[]
    for thickness in [1.,2.,5.,10.]:
     inner=.1-2*thickness*.001;mass=2700*(.1**3-inner**3)
     shields.append(dict(outer_cube_mm=100,aluminium_density_kg_m3=2700,wall_thickness_mm=thickness,closed_six_wall_mass_kg=mass,internal_side_mm=inner*1000,remaining_cavity_litre=inner**3*1000,areal_density_g_cm2=.27*thickness,no_aperture_or_mounts=True,no_radiation_transport_prediction=True))
    write('shield_mass_geometry.csv',shields)
    calibrations=[]
    for rho in [0.,.5,1.]:
     variance=.05**2+.05**2-2*rho*.05*.05+.03**2+.03**2+.02**2+.02**2
     calibrations.append(dict(single_channel_reference_standard_fraction=.05,reference_correlation=rho,relative_field_term=.03,relative_spectral_term=.03,relative_drift_term=.02,relative_temperature_term=.02,combined_ratio_calibration_standard_fraction=math.sqrt(variance),illustrative_total_ratio_target=.1,minimum_equal_background_free_signal_count_per_channel=math.ceil(2/(.1**2-variance)),allocation_is_assumed_not_measured=True))
    write('calibration_allocation.csv',calibrations)
    result=dict(analysis_date='2026-10-22',test_definition=dict(tested_window_start_rate_s=50000,search_window_families=3,observing_seconds_per_day=86400,test_starts_per_day=TESTS_PER_DAY,expected_false_triggers_per_observing_day_allocation=.1,per_test_alpha=ALPHA_OPERATIONAL,union_bound_not_independence_assumption=True),counts=dict(threshold_scenarios=len(thresholds),ratio_examples=len(ratio),onoff_examples=len(onoff),coverage_trials=sum(r['trials']for r in coverage)),threshold_examples=[r for r in thresholds if r['diameter_mm']==5 and r['full_fov_deg']==60 and r['dark_cps']==50 and r['criterion']=='0p1_false_per_observing_day'],estimated_baseline_guard=guards,null_ratio_count_targets=count_targets,onoff_examples=[r for r in onoff if r['kappa_validated_lower']==1 and r['representative_uv_signal_count']in [100,1000]and r['background_per_channel_on']in [0,100,1000]],shield_mass=shields,calibration_allocation=calibrations,validation=dict(decimal_crosschecks=len(check),poisson_zero_count_95pct_upper=poisson_upper(0,.05),all_threshold_tail_brackets_passed=True,all_power_targets_passed=True,ratio_coverage=coverage),limitations=['Known-background count thresholds are optimistic when the background must be estimated. The separate conservative ratio calculation includes independent off-source background counts.','The false-trigger allocation is chosen for illustration; particle bursts, background evolution, detector dead time and adaptive selection are excluded. The expected-number union bound does not require independent overlapping tests.','Sky radiance, efficiency, dark rates and 1 kR auroral filling are planning inputs, not measured mission values.','Ratio intervals require independent Poisson bins, stable on/off backgrounds, common event acceptance and a validated deterministic lower bound on relative photon response. Correlated avalanches/analogue gains require another likelihood.','No physical minimum thermal source red/UV ratio is supplied. An upper limit cannot by itself identify an event as nonthermal.','Shield geometry is only a mass/volume calculation, not a dose, displacement-damage or single-event model.'])
    (HERE/'count_requirements.json').write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf-8')
    print(json.dumps(result['counts'],indent=2));print(json.dumps(result['threshold_examples'],indent=2))

if __name__=='__main__':main()