from pathlib import Path import os for _k in ("OMP_NUM_THREADS","OPENBLAS_NUM_THREADS","MKL_NUM_THREADS"): os.environ[_k]="1" import sys, math, json, csv, time, datetime, hashlib, gzip, argparse import numpy as np P=Path(__file__).resolve().parent OLD=Path("D:/Photometer_research_20261022/research") sys.path.insert(0,str(OLD/"cloud_transport/private_readers")) from scipy.stats import binom from scipy.optimize import minimize_scalar T=0.003 T_PRE=0.010 RISES_US=[30,80,200] COUNTS=[100,1000,10000] BACKGROUNDS=[50.,1e4,1e6] SAMPLE_KHZ=[10,50,100] TRIALS=1000 BASE_SEED=20261022081001 ALPHA_PER_TEST=2.866515718791933e-7 MARGIN=math.log(19.) ONSET_GRID=np.arange(-150.,150.0001,5.)*1e-6 START_UTC=datetime.datetime.now(datetime.timezone.utc).isoformat() def primitive(t,d): z=np.clip(np.asarray(t,dtype=float),0,T) return np.where(z.001-1e-12)&(centers<.0025-1e-12))[0] assert len(early)==round(.0008/dt) and len(late)==round(.0015/dt) durations=np.r_[np.full(len(early),dt),len(late)*dt,T_PRE] models=[] for rise in RISES_US: q_all=source_integrals(edges,ONSET_GRID,rise) for k,onset in enumerate(ONSET_GRID): q=np.r_[q_all[k,early],q_all[k,late].sum(),0.] f=q/durations full=f>=1-1e-11 zero=f<=1e-11 active=~(full|zero) Q=float(q.sum()); D=float(durations.sum()) K=D/Q a=f[active]*K-1. models.append(dict(rise_us=rise,onset_s=float(onset),q=q,f=f,Q=Q,D=D, full=full,zero=zero,active=active,a=a,K=K)) return dict(dt=dt,edges=edges,early=early,late=late,durations=durations,models=models) def fit_one(counts,model,max_iter=40,tol=2e-10): """Exact joint Poisson profile over nonnegative peak rate and background. Normalized source/background weights sum to1. At the MLE their total expectation equals observed N, leaving a concave 1-D mixture fraction theta. """ n=counts.sum(axis=1).astype(float) n0=counts[:,model["zero"]].sum(axis=1).astype(float) np_=counts[:,model["full"]].sum(axis=1).astype(float) nr=counts[:,model["active"]].astype(float) a=model["a"]; c=model["K"]-1. g0=-n0+np_*c+(nr*a).sum(axis=1) theta=np.zeros(len(n)) signal=g0>1e-10 if np.any(signal): ix=np.flatnonzero(signal) N0=n0[ix]; NP=np_[ix]; NR=nr[ix] # Boundary theta=1 is possible only with no counts where q=0. g1=np.full(len(ix),-np.inf) support=N0==0 if np.any(support): g1[support]=NP[support]*c/(1+c)+(NR[support]*a/(1+a)).sum(axis=1) upper=g1>=0 theta[ix[upper]]=1. inner=~upper if np.any(inner): jx=ix[inner]; N0=N0[inner]; NP=NP[inner]; NR=NR[inner] lo=np.zeros(len(jx)); hi=np.ones(len(jx)) # One Newton step from theta=0 is a stable nonnegative initializer. curvature=N0+NP*c*c+(NR*a*a).sum(axis=1) x=np.clip(g0[jx]/np.maximum(curvature,1e-300),1e-12,.99) done=np.zeros(len(jx),dtype=bool) for iteration in range(max_iter): den=1+x[:,None]*a g=-N0/(1-x)+NP*c/(1+x*c)+(NR*a/den).sum(axis=1) h=-N0/(1-x)**2-NP*c*c/(1+x*c)**2-(NR*a*a/(den*den)).sum(axis=1) lo=np.where(g>0,x,lo); hi=np.where(g<=0,x,hi) proposal=x-g/np.minimum(h,-1e-300) bad=(~np.isfinite(proposal))|(proposal<=lo)|(proposal>=hi) proposal=np.where(bad,(lo+hi)/2,proposal) step=np.abs(proposal-x) converged=(np.abs(g)/(n[jx]+1)0 v0[nz]=n0[nz]*np.log1p(-theta[nz]) llr=v0+np_*np.log1p(theta*c)+(nr*np.log1p(theta[:,None]*a)).sum(axis=1) peak=theta*n/model["Q"] b=(1-theta)*n/model["D"] return llr,theta,peak,b def profile(counts,design): ntrial=len(counts) profiles=np.full((ntrial,3),-np.inf) onsets=np.full((ntrial,3),np.nan) peaks=np.full((ntrial,3),np.nan) bs=np.full((ntrial,3),np.nan) for model in design["models"]: k=RISES_US.index(model["rise_us"]) llr,theta,peak,b=fit_one(counts,model) update=llr>profiles[:,k] profiles[update,k]=llr[update] onsets[update,k]=np.where(theta[update]>0,model["onset_s"],np.nan) peaks[update,k]=peak[update];bs[update,k]=b[update] return profiles,onsets,peaks,bs def wilson(k,n,z=1.959963984540054): if n==0: return [None,None] p=k/n; den=1+z*z/n middle=(p+z*z/(2*n))/den half=z*math.sqrt(p*(1-p)/n+z*z/(4*n*n))/den return [max(0.,middle-half),min(1.,middle+half)] def diagnostic_metric(k,n): return {"count":int(k),"denominator":int(n),"fraction":float(k/n) if n else None, "individual_95pct_Wilson_interval":wilson(k,n)} def draw_case(design,n_eq,b,rise_us,seed,ntrial=TRIALS,signal=True): rng=np.random.default_rng(seed) onset=rng.uniform(-100e-6,100e-6,size=ntrial) q=source_integrals(design["edges"],onset,rise_us) source=(n_eq/T)*q if signal else np.zeros_like(q) means=b*design["dt"]+source full=rng.poisson(means) pre=rng.poisson(b*T_PRE,size=ntrial) fitcounts=np.column_stack([full[:,design["early"]],full[:,design["late"]].sum(axis=1),pre]) # Exact conditional on/off test of equal unknown constant background rates. total=full.sum(axis=1) p=binom.sf(total-1,total+pre,.0034/(.0034+T_PRE)) detected=p<=ALPHA_PER_TEST return dict(onset=onset,full=full,pre=pre,fitcounts=fitcounts,total=total, pvalue=p,detected=detected,source=source) def analyze_case(design,n_eq,b,rise_us,seed,ntrial=TRIALS): draw=draw_case(design,n_eq,b,rise_us,seed,ntrial) scores,onsets,peaks,bs=profile(draw["fitcounts"],design) n=len(scores); ar=np.arange(n) best=scores.argmax(axis=1) sorted_scores=np.sort(scores,axis=1) margin3=sorted_scores[:,-1]-sorted_scores[:,-2] accept3=draw["detected"]&(margin3>=MARGIN) fastscore=scores[:,:2].max(axis=1) slow=scores[:,2] bestgroup=(slow>fastscore).astype(int) #0 fast morphology,1 slow morphology margin2=np.abs(slow-fastscore) accept2=draw["detected"]&(margin2>=MARGIN) truth=RISES_US.index(rise_us) truthgroup=0 if truth<2 else 1 pred3=np.where(~draw["detected"],4,np.where(~accept3,3,best)) pred2=np.where(~draw["detected"],3,np.where(~accept2,2,bestgroup)) actual=n_eq*(1-rise_us*1e-6/.8/(2*T)) early={f"expected_source_first_{us}us_PE":float((n_eq/T)*primitive(us*1e-6,rise_us*1e-6/.8)) for us in [30,80,200]} row={"N_plateau_equivalent_3ms_PE":n_eq,"background_cps":b,"sample_rate_khz":round(1/design["dt"]/1000), "truth_receiver_rise_10_90_us":rise_us,"trials":ntrial,"seed":seed, "actual_total_source_PE":actual,"peak_source_cps":n_eq/T,**early, "expected_background_first_200us_PE":b*.0002, "expected_source_in_gate_PE_mean":float(draw["source"].sum(axis=1).mean()), "expected_background_in_gate_PE":b*.0034, "expected_background_prewindow_PE":b*T_PRE, "peak_total_rate_times_17ns":(b+n_eq/T)*17e-9, "count_loss_over_1pct_stationary_mean_flag":(b+n_eq/T)*17e-9>.01/ .99, "detection":diagnostic_metric(draw["detected"].sum(),ntrial), "three_class":{"columns":["rise30","rise80","rise200","rejected_detected","undetected"], "confusion_counts":np.bincount(pred3,minlength=5).tolist(), "accepted":diagnostic_metric(accept3.sum(),ntrial), "wrong_among_accepted":diagnostic_metric((accept3&(best!=truth)).sum(),accept3.sum()), "rejected_among_detected":diagnostic_metric((draw["detected"]&~accept3).sum(),draw["detected"].sum())}, "fast_slow":{"columns":["fast30or80","slow200","rejected_detected","undetected"], "confusion_counts":np.bincount(pred2,minlength=4).tolist(), "accepted":diagnostic_metric(accept2.sum(),ntrial), "wrong_among_accepted":diagnostic_metric((accept2&(bestgroup!=truthgroup)).sum(),accept2.sum()), "rejected_among_detected":diagnostic_metric((draw["detected"]&~accept2).sum(),draw["detected"].sum())}, "oracle_same_data_three_class":{"error":diagnostic_metric((best!=truth).sum(),ntrial), "qualification":"Name retained only for compatibility; this is unthresholded practical profile choice, NOT an oracle."}, "fitted_diagnostics":{"median_peak_rate_cps":float(np.median(peaks[ar,best])), "median_background_cps":float(np.median(bs[ar,best])), "median_onset_error_us":float(np.nanmedian((onsets[ar,best]-draw["onset"])*1e6)), "median_abs_onset_error_us":float(np.nanmedian(abs(onsets[ar,best]-draw["onset"])*1e6)), "onset_grid_boundary_fraction":float(np.mean(abs(onsets[ar,best])>=149.99e-6)), "median_three_class_score_margin":float(np.median(margin3)), "median_fast_slow_score_margin":float(np.median(margin2))}} del row["oracle_same_data_three_class"] row["unthresholded_practical_three_class_error"]=diagnostic_metric((best!=truth).sum(),ntrial) trials={"onset_s":draw["onset"],"pre_count":draw["pre"],"gate_count":draw["total"], "detection_pvalue":draw["pvalue"],"detected":draw["detected"], "scores":scores,"best_rise_index":best,"margin3":margin3,"accept3":accept3, "best_group":bestgroup,"margin2":margin2,"accept2":accept2, "fitted_onset_s":onsets[ar,best],"fitted_peak_cps":peaks[ar,best], "fitted_background_cps":bs[ar,best]} return row,trials def check_math(): checks=[] for khz in SAMPLE_KHZ: design=setup(khz) for rise in RISES_US: onsets=np.array([-.0001,-.000037,0,.000099]) q=source_integrals(design["edges"],onsets,rise) expected=T-rise*1e-6/.8/2 assert np.max(abs(q.sum(axis=1)-expected))<1e-14 checks.append({"test":"full_gate_exact_source_integral","sample_rate_khz":khz,"passed":True}) design=setup(50) rng=np.random.default_rng(5821031) test_counts=np.vstack([np.zeros(len(design["durations"])), rng.poisson(np.r_[np.full(40,3.),200.,500.],size=(8,42))]) # Independent scalar bounded optimization uses the full normalized mixture. errors=[] for model in design["models"][::13]: ll,theta,_,_=fit_one(test_counts,model) h=model["q"]/model["Q"]; w=design["durations"]/model["D"] for j,n in enumerate(test_counts): a=h/w-1 def objective(x): valid=n>0 return -float((n[valid]*np.log1p(x*a[valid])).sum()) opt=minimize_scalar(objective,bounds=(0,1-1e-13),method="bounded", options={"xatol":1e-12}) best=min(0.,opt.fun,objective(1-1e-13)) errors.append(abs(ll[j]+best)) assert max(errors)<2e-6 checks.append({"test":"profile_MLE_vs_independent_scalar_optimizer","comparisons":len(errors), "max_absolute_loglikelihood_difference":max(errors),"passed":True}) # Same fitted parameters imply expected count sum equals actual total. model=design["models"][60] ll,theta,A,b=fit_one(test_counts,model) lam=A[:,None]*model["q"]+b[:,None]*design["durations"] assert np.max(abs(lam.sum(axis=1)-test_counts.sum(axis=1)))<1e-8 assert np.all(ll>=-1e-8) checks.append({"test":"joint_MLE_count_scale_identity_and_LLR_nonnegative","passed":True}) # Time shifts wholly inside fit windows leave late plateau flat. assert all(np.allclose(m["f"][-2],1) and m["f"][-1]==0 for m in design["models"]) checks.append({"test":"late_plateau_and_prewindow_source_exclusion","passed":True}) return checks def csv_rows(rows): for r in rows: yield {k:r[k] for k in ["N_plateau_equivalent_3ms_PE","background_cps","sample_rate_khz", "truth_receiver_rise_10_90_us","trials","seed","actual_total_source_PE", "expected_source_first_30us_PE","expected_source_first_80us_PE", "expected_source_first_200us_PE","expected_background_first_200us_PE", "peak_total_rate_times_17ns"]} | { "detected":r["detection"]["count"],"detection_efficiency":r["detection"]["fraction"], "accepted3":r["three_class"]["accepted"]["count"], "wrong3":r["three_class"]["wrong_among_accepted"]["count"], "error_among_accepted3":r["three_class"]["wrong_among_accepted"]["fraction"], "accepted_fast_slow":r["fast_slow"]["accepted"]["count"], "wrong_fast_slow":r["fast_slow"]["wrong_among_accepted"]["count"], "error_among_accepted_fast_slow":r["fast_slow"]["wrong_among_accepted"]["fraction"], "rejected_detected_fast_slow":r["fast_slow"]["confusion_counts"][2], "undetected":r["fast_slow"]["confusion_counts"][3]} def save_csv(path,rows): rows=list(rows) with path.open("w",newline="",encoding="utf-8") as f: w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows) def run(mode): P.mkdir(parents=True,exist_ok=True) start=time.perf_counter() checks=check_math() print(json.dumps({"QA":checks}),flush=True) if mode=="qa": (P/"math_QA.json").write_text(json.dumps({"checks":checks,"passed":True},indent=2)) return if mode=="benchmark": d=setup(100) r,tr=analyze_case(d,1000,1e6,80,BASE_SEED,1000) print(json.dumps({"runtime_s":time.perf_counter()-start,"example":r}),flush=True) return rows=[];case=0 (P/"synthetic_trial_summaries").mkdir(exist_ok=True) for khz in SAMPLE_KHZ: d=setup(khz) for n_eq in COUNTS: for b in BACKGROUNDS: for rise in RISES_US: seed=BASE_SEED+case t0=time.perf_counter() r,tr=analyze_case(d,n_eq,b,rise,seed) r["case_id"]=case;r["runtime_s"]=time.perf_counter()-t0 rows.append(r) np.savez_compressed(P/"synthetic_trial_summaries"/f"case_{case:03d}.npz",**tr) print(json.dumps({"case_id":case,"khz":khz,"N":n_eq,"b":b,"rise":rise, "runtime_s":r["runtime_s"],"detected":r["detection"]["count"], "accepted_binary":r["fast_slow"]["accepted"]["count"], "wrong_binary":r["fast_slow"]["wrong_among_accepted"]["count"]}),flush=True) case+=1 noise=[] for khz in SAMPLE_KHZ: d=setup(khz) for j,b in enumerate(BACKGROUNDS): seed=BASE_SEED+10000+100*khz+j dr=draw_case(d,0,b,80,seed,1000,signal=False) noise.append({"sample_rate_khz":khz,"background_cps":b,"seed":seed, "trials":1000,"false_detections":int(dr["detected"].sum()), "individual_95pct_Wilson_interval":wilson(dr["detected"].sum(),1000), "qualification":"1000 noise trials cannot empirically establish a 2.87e-7 tail. Exact conditional test supplies the per-event-gate level."}) save_csv(P/"practical_morphology_scenarios.csv",csv_rows(rows)) save_csv(P/"noise_only_controls.csv",[{k:v for k,v in x.items() if k not in ["individual_95pct_Wilson_interval","qualification"]} for x in noise]) out={"schema_version":"1.0","status":"EXECUTED_CONDITIONAL_IDEAL_POISSON_PRACTICAL_MORPHOLOGY_PROTOTYPE", "start_utc":START_UTC,"finished_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(), "runtime_s":time.perf_counter()-start,"seed_base":BASE_SEED, "scenario_count":len(rows),"trials_each_scenario":TRIALS,"total_signal_trials":len(rows)*TRIALS, "parameters":{"receiver_rises_10_90_us":RISES_US,"plateau_equivalent_PE_in_3ms":COUNTS, "background_cps":BACKGROUNDS,"sample_rate_khz":SAMPLE_KHZ, "pulse_end_after_onset_ms":3.,"prewindow_ms":10., "prewindow_position_ms":[-10.5,-.5],"candidate_gate_ms":[-.2,3.2], "morphology_early_window_ms":[-.2,.6],"late_plateau_window_ms":[1.,2.5], "true_onset_distribution_us":[-100,100],"onset_grid_us":[-150,150,5], "alpha_per_predefined_gate":ALPHA_PER_TEST, "acceptance_profile_loglikelihood_margin":MARGIN, "deadtime_diagnostic_ns":17.}, "algorithm":{"label_definition":"Linear receiver-input expected PE-rate morphology after any assumed transport, not intrinsic discharge history or cloud depth.", "waveform":"A linear0-to100% ramp lasting rise_10_90/0.8, then unit plateau, with total pulse ending3ms after onset.", "amplitude":"Peak rate A=N_plateau_equivalent/3ms. Actual integrated source counts=A*(3ms-ramp_duration/2); source and total/background counts are separate.", "sampling":"Independent Poisson hardware-bin counts from exact primitive differences. Onset uniform[-100,100]us spans integer phase cycles at every sample rate.", "background":"Constant unknown rate. One explicit10ms background-only count plus the fit-window counts enter a joint Poisson profile likelihood.", "likelihood":"For each rise and onset, lambda_j=A*q_j+b*d_j, A,b>=0, including q_pre=0,d_pre=10ms. Normalize h=q/sum(q),w=d/sum(d). Their MLE total expectation equals total observed counts N; profile the concave mixture p=(1-theta)w+theta*h for theta in[0,1].", "onset":"Global exhaustive5us grid[-150,150]us for all three morphologies. Brightness and background are profiled for every onset; fit uses early rise and late plateau, excludes final cutoff.", "detection":"Exact on/off conditional-binomial test: N_gate|N_gate+M_pre ~Binomial(N_gate+M_pre,T_gate/(T_gate+T_pre)) under constant-background null. Per predefined gate only.", "classification":"Separate3-class and fast30or80-vs-slow200 choices. Accept only if detected AND winning profile score exceeds runner/opposite group by ln19. This fixed margin is not a95% posterior/probability guarantee.", "training_test_separation":"No fitting of the margin, waveform grid or decision rule to simulation labels. Every performance trial uses a fresh seeded Poisson realization. No labeled training set exists.", "response":"No electronic or arrival-delay filter. Analog post-counting RC response would correlate samples and invalidate this independent-Poisson likelihood.", "MC_intervals":"Individual two-sided95% Wilson score intervals for reported binomial rates. Not simultaneous over all scenarios, not uncertainty in the true event population."}, "math_QA":checks,"scenarios":rows,"noise_only_controls":noise, "limitations":[ "Conditional synthetic morphology experiment, not validated flight performance or cloud-depth inference.", "Exact three templates, known constant background law and known late plateau. Model mismatch, real source variation, cloud-kernel variation and actual detector response are absent.", "Unknown onset is restricted to an already available200us coarse bracket; identifying that bracket in a blind stream is outside scope.", "The fixed3.4ms detection gate covers each simulated pulse; constructing such gates blindly and overlapping-window search trials remain operational inputs.", "The nonnegative nuisance fit can favor flexible fast templates; reported Monte Carlo confusion counts, not the score margin, quantify performance under this generator.", "Only linear independent photon registration is simulated. Peak rate times17ns flags an assumed counter's possible occupancy; it is not a deadtime correction or restoration of Poisson statistics.", "Per-gate false alarm alpha is not a mission false-trigger rate. A specified finite test count permits a Bonferroni familywise upper bound; actual streaming trials are unknown.", "PE counts include an unspecified end-to-end photon-detection efficiency. Aperture photons require a separately verified spectral/angular efficiency.", "No molecular source yield, source brightness distribution, population completeness or optical duty/exposure model is supplied. Useful mission event rate is NOT COMPUTABLE.", "A real morphology classifier requires end-to-end hardware waveforms and separately labeled optical/radio control populations.", "Fitted-grid errors and Monte Carlo intervals do not cover systematic model uncertainty. Samples are10/20/100us integrations; sampling rate is not analog bandwidth."] } (P/"practical_morphology_results.json").write_text(json.dumps(out,indent=2),encoding="utf-8") print(json.dumps({"finished":True,"runtime_s":out["runtime_s"],"scenarios":len(rows)}),flush=True) if __name__=="__main__": parser=argparse.ArgumentParser();parser.add_argument("--mode",choices=["qa","benchmark","full"],default="full") run(parser.parse_args().mode)