Review index · Download raw source
reference/analytic/count_requirements.py
1"""Exact Poisson planning requirements; no measured receiver is asserted.23Only Python's standard library and numpy are required. No normal approximation4is used in the threshold/95%-power or one-sided count interval calculations.5"""6from pathlib import Path7import math, json, csv8from functools import lru_cache9import numpy as np1011HERE=Path(__file__).resolve().parent12ROOT=HERE.parents[2]13EPH=6.62607015e-34*299792458/337.1e-914ALPHA_SINGLE=2.866515718791933e-715TESTS_PER_DAY=50000*86400*316ALPHA_OPERATIONAL=.1/TESTS_PER_DAY1718def poisson_sf(n,mu):19 """P(N>=n). Positive recurrences in the smaller tail avoid cancellation."""20 n=int(n)21 if n<=0:return 1.22 if mu==0:return 0.23 if mu<0:raise ValueError('negative Poisson mean')24 if n-1<mu:25 k=n-1;term=math.exp(-mu+k*math.log(mu)-math.lgamma(k+1));s=term26 while k>0:27 term*=k/mu;k-=1;s+=term28 if term<1e-15*s:break29 return max(0.,min(1.,1-s))30 k=n;term=math.exp(-mu+k*math.log(mu)-math.lgamma(k+1));s=term31 if term==0:return 0.32 for _ in range(1000000):33 k+=1;term*=mu/k;s+=term34 if term<1e-15*s:break35 else:raise RuntimeError('Poisson recurrence did not converge')36 return max(0.,min(1.,s))3738def poisson_cdf(n,mu):39 if n<0:return 0.40 if mu==0:return 1.41 n=int(n)42 if n<mu:43 k=n;term=math.exp(-mu+k*math.log(mu)-math.lgamma(k+1));s=term44 while k>0:45 term*=k/mu;k-=1;s+=term46 if term<1e-15*s:break47 return max(0.,min(1.,s))48 return 1-poisson_sf(n+1,mu)4950@lru_cache(None)51def critical_count(bg,alpha):52 low=0;high=max(1,int(bg+12*math.sqrt(bg+1)+100))53 while poisson_sf(high,bg)>alpha:high*=254 while high-low>1:55 mid=(high+low)//256 if poisson_sf(mid,bg)<=alpha:high=mid57 else:low=mid58 return high5960@lru_cache(None)61def count_at_power(bg,alpha,power=.95):62 n=critical_count(bg,alpha);lo=0.;hi=max(10.,n-bg+8*math.sqrt(n+1)+100)63 while poisson_sf(n,bg+hi)<power:hi*=264 for _ in range(55):65 mid=(lo+hi)/266 if poisson_sf(n,bg+mid)>=power:hi=mid67 else:lo=mid68 return n,hi6970@lru_cache(None)71def poisson_upper(n,alpha):72 lo=0.;hi=max(10.,n+12*math.sqrt(n+1)+100)73 while poisson_cdf(n,hi)>alpha:hi*=274 for _ in range(55):75 mid=(lo+hi)/276 if poisson_cdf(n,mid)<alpha:hi=mid77 else:lo=mid78 return hi7980@lru_cache(None)81def poisson_lower(n,alpha):82 if n==0:return 0.83 lo=0.;hi=float(n)84 for _ in range(55):85 mid=(lo+hi)/286 if poisson_sf(n,mid)>alpha:hi=mid87 else:lo=mid88 return lo8990def conservative_ratio_upper(nuv,nred,muv,mred,on_s,off_s,kappa_lower,alpha=.05):91 """Four marginal one-sided intervals; union bound gives >=1-alpha coverage.9293 Coverage is unconditional for a prespecified window. A UV-trigger-selected94 catalog needs a separate selection-aware coverage argument.95 kappa_lower is a deterministic validated lower bound, not a nominal estimate96 or a confidence interval silently used without an additional error budget.97 """98 a=alpha/4;scale=on_s/off_s99 sr=max(0.,poisson_upper(nred,a)-scale*poisson_lower(mred,a))100 su=max(0.,poisson_lower(nuv,a)-scale*poisson_upper(muv,a))101 return (None if su==0 else sr/su/kappa_lower),sr,su102103def write(name,rows):104 with (HERE/name).open('w',newline='',encoding='utf-8')as f:105 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)106107def main():108 thresholds=[]109 for diameter in [1.,5.,25.]:110 area=math.pi*(diameter*.001)**2/4111 for fov in [60.,90.]:112 omega=math.pi*math.sin(math.radians(fov/2))**2113 quiet=.2*area*omega*4*3e-9/EPH114 for aurora_kr in [0.,1.]:115 aurora=.2*area*omega*aurora_kr*1000*1e10/(4*math.pi)116 for dark in [50.,1e4,1e6]:117 rate=quiet+aurora+dark118 for window in [30e-6,3e-3]:119 bg=rate*window120 for name,alpha in [('single_5sigma_tail',ALPHA_SINGLE),('0p1_false_per_observing_day',ALPHA_OPERATIONAL)]:121 n,s=count_at_power(bg,alpha)122 pfa=poisson_sf(n,bg)123 assert pfa<=alpha*(1+1e-9) and poisson_sf(n-1,bg)>alpha124 assert abs(poisson_sf(n,bg+s)-.95)<2e-10125 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))126 write('exact_detection_requirements.csv',thresholds)127 ratio=[]128 for kappa in [.1,.6223529411764706,1.,2.]:129 for confidence in [.95,.99]:130 for n in [10,30,100,300,1000,10000]:131 a=1-confidence;upper=math.expm1(-math.log(a)/n)/kappa132 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))133 write('zero_red_exact_ratio_bounds.csv',ratio)134 count_targets=[]135 for kappa in [.6223529411764706,1.]:136 for r in [.01,.1,1.]:137 n=math.floor(-math.log(.05)/math.log1p(kappa*r))+1138 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))139 write('null_channel_count_requirements.csv',count_targets)140 onoff=[]141 for kappa in [.6223529411764706,1.]:142 for s in [10,30,100,300,1000,10000]:143 for bg in [0.,1.,10.,100.,1000.]:144 nuv=int(round(s+bg));nred=int(round(bg));muv=mred=int(round(3*bg))145 ru,sru,sul=conservative_ratio_upper(nuv,nred,muv,mred,.003,.009,kappa)146 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))147 write('background_aware_ratio_bounds.csv',onoff)148 # A directly implementable guard for an estimated, stationary baseline.149 # This trades power for a proof by two one-sided probability bounds.150 guards=[]151 for row in thresholds:152 if row['diameter_mm']==5 and row['full_fov_deg']==60 and row['dark_cps']==50 and row['criterion']=='0p1_false_per_observing_day':153 off=.009;on=row['window_s'];m=int(round(row['background_rate_cps']*off));alpha=row['per_test_alpha']154 bg_upper=on/off*poisson_upper(m,alpha/2);crit=critical_count(bg_upper,alpha/2)155 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))156 write('estimated_baseline_trigger_guard.csv',guards)157 # Validate Poisson recurrence against a separately evaluated 50-digit series.158 from decimal import Decimal,localcontext159 check=[]160 with localcontext()as context:161 context.prec=50162 for mu in [.001,.15,1.,5.,20.,100.]:163 for n in [1,2,5,10,25,100,150]:164 m=Decimal(str(mu));term=(-m).exp();cdf=Decimal(0)165 for k in range(n):166 if k:term*=m/Decimal(k)167 cdf+=term168 sf=float(1-cdf);actual=poisson_sf(n,mu)169 assert abs(sf-actual)<=max(2e-13,abs(sf)*1e-10),(mu,n,sf,actual)170 check.append(dict(mean=mu,n=n,decimal_sf=sf,recurrence_sf=actual))171 assert abs(poisson_upper(0,.05)+math.log(.05))<1e-12172 assert abs(poisson_upper(1,.05)-4.743864518390577)<1e-10173 # Coverage audit for the conservative ratio bound. This is an implementation174 # check; the union-bound argument supplies coverage under the stated model.175 rng=np.random.default_rng(2026102201);coverage=[];coverage_raw=[]176 for su_true in [30.,300.]:177 for bg in [0.,10.,100.]:178 for r_true in [.01,.1,1.]:179 kappa=.6223529411764706;rejected=0;finite=0;ntrial=1024180 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)])181 for trial,(nu,nr,mu,mr) in enumerate(draws):182 ru,_,_=conservative_ratio_upper(int(nu),int(nr),int(mu),int(mr),.003,.009,kappa)183 finite+=ru is not None;rejected+=ru is not None and ru<r_true184 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))185 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))186 write('ratio_coverage_audit.csv',coverage)187 write('ratio_coverage_trials.csv',coverage_raw)188 import gzip189 with gzip.open(HERE/'ratio_coverage_trials.csv.gz','wb')as f:f.write((HERE/'ratio_coverage_trials.csv').read_bytes())190 shields=[]191 for thickness in [1.,2.,5.,10.]:192 inner=.1-2*thickness*.001;mass=2700*(.1**3-inner**3)193 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))194 write('shield_mass_geometry.csv',shields)195 calibrations=[]196 for rho in [0.,.5,1.]:197 variance=.05**2+.05**2-2*rho*.05*.05+.03**2+.03**2+.02**2+.02**2198 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))199 write('calibration_allocation.csv',calibrations)200 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.'])201 (HERE/'count_requirements.json').write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf-8')202 print(json.dumps(result['counts'],indent=2));print(json.dumps(result['threshold_examples'],indent=2))203204if __name__=='__main__':main()205