Review index · Download raw source
reference/statistics/oracle_reference/calculate_ideal_poisson_rise.py
1from pathlib import Path2import numpy as np3import math,csv,json,datetime,hashlib4P=Path(__file__).resolve().parent5P.mkdir(parents=True,exist_ok=True)6T=0.0037RISES=[30e-6,80e-6,200e-6]8N_EQ=[100,300,1000,3000,10000]9BACKGROUND=[50.0,1e4,1e6,1e8]10BINS=[20e-6,100e-6]11PHASES=[0.0,.25,.5,.75]12START=datetime.datetime.now(datetime.timezone.utc).isoformat()13def integral(t,d):14 # Exact antiderivative of a linear ramp beginning at0, then plateau,15 # with the entire optical pulse ending at T (not T after the ramp).16 z=np.clip(np.asarray(t,dtype=float),0,T)17 return np.where(z<d,z*z/(2*d),z-d/2)18def experiment(n_eq,b,dt,phase,rise):19 d=rise/.820 amp=n_eq/T21 nbin=int(math.ceil(T/dt+phase-1e-12))22 edges=(np.arange(nbin+1,dtype=float)-phase)*dt23 source=amp*np.diff(integral(edges,d))24 means=b*dt+source25 assert np.min(source)>-1e-1026 assert abs(source.sum()-amp*(T-d/2))<1e-827 return means,source,edges28def affinity_bounds(mu,nu):29 h=float(.5*np.sum((np.sqrt(mu)-np.sqrt(nu))**2))30 bc=math.exp(-h) if h<746 else 0.31 sqrt1=math.sqrt(-math.expm1(-2*h))32 log_lower=-2*h-math.log(2)-math.log1p(sqrt1)33 log_upper=-h-math.log(2)34 lower=math.exp(log_lower) if log_lower>-746 else 0.35 upper=math.exp(log_upper) if log_upper>-746 else 0.36 return {'minus_log_affinity':h,'affinity':bc,'bayes_error_lower':lower,37 'bayes_error_upper':upper,'log_bayes_error_lower':log_lower,38 'log_bayes_error_upper':log_upper,'numeric_underflow':log_lower<-745 or log_upper<-745}39early=[]40for n in N_EQ:41 for rise in RISES:42 d=rise/.8;A=n/T43 early.append({'N_plateau_equivalent_3ms_PE':n,'rise_10_90_us':rise*1e6,44 'ramp_0_100_us':d*1e6,'peak_source_cps':A,45 'actual_total_source_PE':A*(T-d/2),46 'expected_source_0_30us_PE':float(A*integral(30e-6,d)),47 'expected_source_0_200us_PE':float(A*integral(200e-6,d)),48 'probability_zero_source_PE_0_200us':math.exp(-float(A*integral(200e-6,d)))})49pairs=[];detect=[];mc_inputs=[]50for n in N_EQ:51 for b in BACKGROUND:52 for dt in BINS:53 for phase in PHASES:54 slow,_,edges=experiment(n,b,dt,phase,200e-6)55 for fast_rise in RISES[:2]:56 fast,_,_=experiment(n,b,dt,phase,fast_rise)57 stats=affinity_bounds(fast,slow)58 row={'N_plateau_equivalent_3ms_PE':n,'background_cps':b,'bin_us':dt*1e6,59 'onset_phase_fraction':phase,'fast_rise_10_90_us':fast_rise*1e6,60 'slow_rise_10_90_us':200.0,'observation_left_us':edges[0]*1e6,61 'observation_right_us':edges[-1]*1e6,'hardware_bins':len(fast),**stats}62 pairs.append(row)63 if n in [100,1000,10000] and b in [50.,1e6,1e8] and phase in [0.,.5]:64 mc_inputs.append((dict(row),fast,slow))65 for rise in RISES:66 signal,_,edges=experiment(n,b,dt,phase,rise)67 pure_bg=np.full(len(signal),b*dt)68 detect.append({'N_plateau_equivalent_3ms_PE':n,'background_cps':b,69 'bin_us':dt*1e6,'onset_phase_fraction':phase,'rise_10_90_us':rise*1e6,70 'known_signal_versus_background_only':True,**affinity_bounds(signal,pure_bg)})71# Exact Poisson controls, without Gaussian count approximations.72controls=[]73for mu,nu in [(0.,.3),(.1,.2),(.1,3.),(1.,10.),(10.,30.)]:74 def pmf(k,m):75 return (1. if k==0 else 0.) if m==0 else math.exp(-m+k*math.log(m)-math.lgamma(k+1))76 p=np.array([pmf(k,mu) for k in range(300)])77 q=np.array([pmf(k,nu) for k in range(300)])78 bc=float(np.sqrt(p*q).sum()); pe=float(.5*np.minimum(p,q).sum())79 analytic=affinity_bounds(np.array([mu]),np.array([nu]))80 assert abs(bc-analytic['affinity'])<1e-1381 assert analytic['bayes_error_lower']-1e-13<=pe<=analytic['bayes_error_upper']+1e-1382 assert abs(p.sum()-1)<1e-13 and abs(q.sum()-1)<1e-1383 controls.append({'mu':mu,'nu':nu,'summed_affinity':bc,'exact_truncated_bayes_error':pe,**analytic})84# Coarse bins are sums of five aligned fine bins for the four chosen phases.85# Extra bins outside pulse have identical background under both hypotheses.86coarse_checks=087lookup={(r['N_plateau_equivalent_3ms_PE'],r['background_cps'],r['bin_us'],88 r['onset_phase_fraction'],r['fast_rise_10_90_us']):r for r in pairs}89for r in pairs:90 if r['bin_us']==100.:91 f=lookup[(r['N_plateau_equivalent_3ms_PE'],r['background_cps'],20.,92 r['onset_phase_fraction'],r['fast_rise_10_90_us'])]93 assert r['minus_log_affinity']<=f['minus_log_affinity']+1e-1094 coarse_checks+=195# Limited Monte Carlo with exact NumPy Poisson generators and exact log likelihood.96# These estimates are not called exact Bayes errors.97rng=np.random.default_rng(20261022004838)98N_PER_HYPOTHESIS=2000099M=len(mc_inputs)100# Two-sided Hoeffding bound for the balanced average of40000 bounded error scores.101# Union bound gives99% coverage across all M displayed estimates.102mc_eps=math.sqrt(math.log(2*M/.01)/(4*N_PER_HYPOTHESIS))103mc=[]104for row,mu,nu in mc_inputs:105 active=np.abs(mu-nu)>1e-12*np.maximum(1.,np.maximum(mu,nu))106 a=mu[active];c=nu[active]107 logratio=np.log(a/c);constant=float((a-c).sum())108 n1=rng.poisson(a,size=(N_PER_HYPOTHESIS,len(a)))109 n2=rng.poisson(c,size=(N_PER_HYPOTHESIS,len(c)))110 ll1=n1@logratio-constant;ll2=n2@logratio-constant111 err1=(ll1<0).astype(float)+.5*(ll1==0)112 err2=(ll2>0).astype(float)+.5*(ll2==0)113 estimate=float(.5*(err1.mean()+err2.mean()))114 # Deterministic bounds need not enclose a finite random estimate, but115 # its stated conservative simultaneous interval must intersect them.116 assert estimate+mc_eps>=row['bayes_error_lower']-1e-12117 assert estimate-mc_eps<=row['bayes_error_upper']+1e-12118 mc.append({**row,'MC_samples_each_hypothesis':N_PER_HYPOTHESIS,119 'MC_error_estimate':estimate,'MC_errors_score_fast':float(err1.sum()),120 'MC_errors_score_slow':float(err2.sum()),'MC_simultaneous_99pct_Hoeffding_halfwidth':mc_eps,121 'MC_interval_lower':max(0.,estimate-mc_eps),'MC_interval_upper':min(.5,estimate+mc_eps),122 'seed':20261022004838})123def save_csv(name,rows):124 with (P/name).open('w',newline='',encoding='utf-8') as f:125 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)126save_csv('source_count_expectations.csv',early)127save_csv('rise_pair_affinity_bounds.csv',pairs)128save_csv('pulse_detection_affinity_bounds.csv',detect)129save_csv('limited_poisson_MC.csv',mc)130save_csv('single_bin_exact_controls.csv',controls)131summary={'status':'EXECUTED_IDEAL_POISSON_KNOWN_HYPOTHESES_TOY_NOT_FLIGHT_PERFORMANCE',132'actual_research_start_utc':START,'calculation_finished_utc':datetime.datetime.now(datetime.timezone.utc).isoformat(),133'parameters':{'pulse_end_s':T,'rise_10_90_s':RISES,'N_plateau_equivalent_PE':N_EQ,134'background_cps':BACKGROUND,'hardware_bin_s':BINS,'onset_phases_fraction':PHASES,135'fixed_scale':'Peak amplitude A=N_plateau_equivalent/3ms is fixed across competing rise hypotheses. Actual integrated source PE differs because ramp areas differ.'},136'counts':{'rise_comparison_rows':len(pairs),'detection_rows':len(detect),137'source_count_rows':len(early),'Monte_Carlo_rows':len(mc),'Monte_Carlo_samples_each_hypothesis_per_row':N_PER_HYPOTHESIS,138'single_bin_exact_distribution_controls':len(controls),'coarse_bin_information_monotonicity_checks':coarse_checks},139'MC_familywise_bound':{'confidence_at_least':.99,'halfwidth':mc_eps,'method':'Hoeffding for40000 bounded independent error scores per scenario, followed by a union bound over all72 scenarios. No Gaussian count or confidence approximation.'},140'proofs':[141'For one Poisson bin, sum_k sqrt(P_mu(k)P_nu(k))=exp[-(mu+nu)/2+sqrt(mu*nu)]=exp[-(sqrt(mu)-sqrt(nu))^2/2]. Independence multiplies affinities, so log affinity sums over bins.',142'Equal-prior optimal error Pe=(1-TV)/2. Since min(p,q)<=sqrt(p*q), 1-TV<=BC; hence Pe<=BC/2.',143'Cauchy-Schwarz applied to |p-q|=|sqrt(p)-sqrt(q)|*(sqrt(p)+sqrt(q)) gives TV<=sqrt(1-BC^2); hence Pe>=(1-sqrt(1-BC^2))/2.',144'Bounds concern fixed simple hypotheses with known parameters. They are not Gaussian approximations, significance thresholds, or measured classification errors. Very small displayed numeric bounds can underflow; natural-log bounds are retained.'145],146'waveform_reference_plane':'The templates specify expected PE rate at the hypothetical linear receiver input, after any assumed atmospheric transport. They are not asserted intrinsic NBE emission histories or cloud-depth labels.',147'high_background_registration_limit':'The1e8cps background is only an ideal Poisson-registration scenario; it is not a realizable17ns nonparalyzable counter operating point. Count-loss response is studied separately.',148'limitations':[149'Known pulse start, rise templates, peak amplitude, background and pulse end. No nuisance fitting or training/generalization error.',150'Poisson PE arrivals, independent hardware-bin counts and ideal linear photon registration. No deadtime, saturation, gain fluctuations, excess noise, afterpulsing, analog filtering or threshold crossings.',151'End-to-end photon detection efficiency is absorbed into PE amplitude. Aperture/source brightness distribution and cloud depth are not specified.',152'Background is known and constant. Actual source and total counts are separate; early source-only counts are not directly labeled in the detected mixture.',153'The exponential optical decay after plateau is not modeled; pulse stops at3ms identically under each hypothesis. This arbitrary shared cutoff adds no artificial inter-hypothesis timing signal.',154'Binning, unobserved timing/nuisance parameters, or additional hypothesis-independent stochastic processing cannot add information relative to an oracle retaining the known waveform and photon times. Correlated physical variation could encode other information, so no universal real-population misclassification claim is made.',155'Binary equal-prior fast-versus-slow error does not establish a three-class classifier, cloud-top/depth inference, event-rate completeness or actual instrument performance.',156'All paired means include full hardware bins straddling onset/termination. Expected source counts in physical0–30/200us windows are analytic integrals, not a claim that a coarse sampled detector can recover partial-bin photon times.'157],158'all_checks_passed':True}159(P/'calculation_summary.json').write_text(json.dumps(summary,indent=2),encoding='utf-8')160examples=[r for r in pairs if r['N_plateau_equivalent_3ms_PE'] in [100,1000] and r['background_cps'] in [50.,1e6] and r['bin_us']==20. and r['onset_phase_fraction']==0.]161print(json.dumps({'counts':summary['counts'],'checks_pass':True,'examples':examples,'early_counts_for_N100':[r for r in early if r['N_plateau_equivalent_3ms_PE']==100]},indent=2))162