from pathlib import Path import numpy as np import math,csv,json,datetime,hashlib P=Path(__file__).resolve().parent P.mkdir(parents=True,exist_ok=True) T=0.003 RISES=[30e-6,80e-6,200e-6] N_EQ=[100,300,1000,3000,10000] BACKGROUND=[50.0,1e4,1e6,1e8] BINS=[20e-6,100e-6] PHASES=[0.0,.25,.5,.75] START=datetime.datetime.now(datetime.timezone.utc).isoformat() def integral(t,d): # Exact antiderivative of a linear ramp beginning at0, then plateau, # with the entire optical pulse ending at T (not T after the ramp). z=np.clip(np.asarray(t,dtype=float),0,T) return np.where(z-1e-10 assert abs(source.sum()-amp*(T-d/2))<1e-8 return means,source,edges def affinity_bounds(mu,nu): h=float(.5*np.sum((np.sqrt(mu)-np.sqrt(nu))**2)) bc=math.exp(-h) if h<746 else 0. sqrt1=math.sqrt(-math.expm1(-2*h)) log_lower=-2*h-math.log(2)-math.log1p(sqrt1) log_upper=-h-math.log(2) lower=math.exp(log_lower) if log_lower>-746 else 0. upper=math.exp(log_upper) if log_upper>-746 else 0. return {'minus_log_affinity':h,'affinity':bc,'bayes_error_lower':lower, 'bayes_error_upper':upper,'log_bayes_error_lower':log_lower, 'log_bayes_error_upper':log_upper,'numeric_underflow':log_lower<-745 or log_upper<-745} early=[] for n in N_EQ: for rise in RISES: d=rise/.8;A=n/T early.append({'N_plateau_equivalent_3ms_PE':n,'rise_10_90_us':rise*1e6, 'ramp_0_100_us':d*1e6,'peak_source_cps':A, 'actual_total_source_PE':A*(T-d/2), 'expected_source_0_30us_PE':float(A*integral(30e-6,d)), 'expected_source_0_200us_PE':float(A*integral(200e-6,d)), 'probability_zero_source_PE_0_200us':math.exp(-float(A*integral(200e-6,d)))}) pairs=[];detect=[];mc_inputs=[] for n in N_EQ: for b in BACKGROUND: for dt in BINS: for phase in PHASES: slow,_,edges=experiment(n,b,dt,phase,200e-6) for fast_rise in RISES[:2]: fast,_,_=experiment(n,b,dt,phase,fast_rise) stats=affinity_bounds(fast,slow) row={'N_plateau_equivalent_3ms_PE':n,'background_cps':b,'bin_us':dt*1e6, 'onset_phase_fraction':phase,'fast_rise_10_90_us':fast_rise*1e6, 'slow_rise_10_90_us':200.0,'observation_left_us':edges[0]*1e6, 'observation_right_us':edges[-1]*1e6,'hardware_bins':len(fast),**stats} pairs.append(row) if n in [100,1000,10000] and b in [50.,1e6,1e8] and phase in [0.,.5]: mc_inputs.append((dict(row),fast,slow)) for rise in RISES: signal,_,edges=experiment(n,b,dt,phase,rise) pure_bg=np.full(len(signal),b*dt) detect.append({'N_plateau_equivalent_3ms_PE':n,'background_cps':b, 'bin_us':dt*1e6,'onset_phase_fraction':phase,'rise_10_90_us':rise*1e6, 'known_signal_versus_background_only':True,**affinity_bounds(signal,pure_bg)}) # Exact Poisson controls, without Gaussian count approximations. controls=[] for mu,nu in [(0.,.3),(.1,.2),(.1,3.),(1.,10.),(10.,30.)]: def pmf(k,m): return (1. if k==0 else 0.) if m==0 else math.exp(-m+k*math.log(m)-math.lgamma(k+1)) p=np.array([pmf(k,mu) for k in range(300)]) q=np.array([pmf(k,nu) for k in range(300)]) bc=float(np.sqrt(p*q).sum()); pe=float(.5*np.minimum(p,q).sum()) analytic=affinity_bounds(np.array([mu]),np.array([nu])) assert abs(bc-analytic['affinity'])<1e-13 assert analytic['bayes_error_lower']-1e-13<=pe<=analytic['bayes_error_upper']+1e-13 assert abs(p.sum()-1)<1e-13 and abs(q.sum()-1)<1e-13 controls.append({'mu':mu,'nu':nu,'summed_affinity':bc,'exact_truncated_bayes_error':pe,**analytic}) # Coarse bins are sums of five aligned fine bins for the four chosen phases. # Extra bins outside pulse have identical background under both hypotheses. coarse_checks=0 lookup={(r['N_plateau_equivalent_3ms_PE'],r['background_cps'],r['bin_us'], r['onset_phase_fraction'],r['fast_rise_10_90_us']):r for r in pairs} for r in pairs: if r['bin_us']==100.: f=lookup[(r['N_plateau_equivalent_3ms_PE'],r['background_cps'],20., r['onset_phase_fraction'],r['fast_rise_10_90_us'])] assert r['minus_log_affinity']<=f['minus_log_affinity']+1e-10 coarse_checks+=1 # Limited Monte Carlo with exact NumPy Poisson generators and exact log likelihood. # These estimates are not called exact Bayes errors. rng=np.random.default_rng(20261022004838) N_PER_HYPOTHESIS=20000 M=len(mc_inputs) # Two-sided Hoeffding bound for the balanced average of40000 bounded error scores. # Union bound gives99% coverage across all M displayed estimates. mc_eps=math.sqrt(math.log(2*M/.01)/(4*N_PER_HYPOTHESIS)) mc=[] for row,mu,nu in mc_inputs: active=np.abs(mu-nu)>1e-12*np.maximum(1.,np.maximum(mu,nu)) a=mu[active];c=nu[active] logratio=np.log(a/c);constant=float((a-c).sum()) n1=rng.poisson(a,size=(N_PER_HYPOTHESIS,len(a))) n2=rng.poisson(c,size=(N_PER_HYPOTHESIS,len(c))) ll1=n1@logratio-constant;ll2=n2@logratio-constant err1=(ll1<0).astype(float)+.5*(ll1==0) err2=(ll2>0).astype(float)+.5*(ll2==0) estimate=float(.5*(err1.mean()+err2.mean())) # Deterministic bounds need not enclose a finite random estimate, but # its stated conservative simultaneous interval must intersect them. assert estimate+mc_eps>=row['bayes_error_lower']-1e-12 assert estimate-mc_eps<=row['bayes_error_upper']+1e-12 mc.append({**row,'MC_samples_each_hypothesis':N_PER_HYPOTHESIS, 'MC_error_estimate':estimate,'MC_errors_score_fast':float(err1.sum()), 'MC_errors_score_slow':float(err2.sum()),'MC_simultaneous_99pct_Hoeffding_halfwidth':mc_eps, 'MC_interval_lower':max(0.,estimate-mc_eps),'MC_interval_upper':min(.5,estimate+mc_eps), 'seed':20261022004838}) def save_csv(name,rows): with (P/name).open('w',newline='',encoding='utf-8') as f: w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows) save_csv('source_count_expectations.csv',early) save_csv('rise_pair_affinity_bounds.csv',pairs) save_csv('pulse_detection_affinity_bounds.csv',detect) save_csv('limited_poisson_MC.csv',mc) save_csv('single_bin_exact_controls.csv',controls) summary={'status':'EXECUTED_IDEAL_POISSON_KNOWN_HYPOTHESES_TOY_NOT_FLIGHT_PERFORMANCE', 'actual_research_start_utc':START,'calculation_finished_utc':datetime.datetime.now(datetime.timezone.utc).isoformat(), 'parameters':{'pulse_end_s':T,'rise_10_90_s':RISES,'N_plateau_equivalent_PE':N_EQ, 'background_cps':BACKGROUND,'hardware_bin_s':BINS,'onset_phases_fraction':PHASES, 'fixed_scale':'Peak amplitude A=N_plateau_equivalent/3ms is fixed across competing rise hypotheses. Actual integrated source PE differs because ramp areas differ.'}, 'counts':{'rise_comparison_rows':len(pairs),'detection_rows':len(detect), 'source_count_rows':len(early),'Monte_Carlo_rows':len(mc),'Monte_Carlo_samples_each_hypothesis_per_row':N_PER_HYPOTHESIS, 'single_bin_exact_distribution_controls':len(controls),'coarse_bin_information_monotonicity_checks':coarse_checks}, '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.'}, 'proofs':[ '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.', 'Equal-prior optimal error Pe=(1-TV)/2. Since min(p,q)<=sqrt(p*q), 1-TV<=BC; hence Pe<=BC/2.', '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.', '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.' ], '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.', '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.', 'limitations':[ 'Known pulse start, rise templates, peak amplitude, background and pulse end. No nuisance fitting or training/generalization error.', '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.', 'End-to-end photon detection efficiency is absorbed into PE amplitude. Aperture/source brightness distribution and cloud depth are not specified.', 'Background is known and constant. Actual source and total counts are separate; early source-only counts are not directly labeled in the detected mixture.', '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.', '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.', '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.', '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.' ], 'all_checks_passed':True} (P/'calculation_summary.json').write_text(json.dumps(summary,indent=2),encoding='utf-8') examples=[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.] print(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))