Review index · Download raw source
reference/statistics/practical_morphology.py
1from pathlib import Path2import os3for _k in ("OMP_NUM_THREADS","OPENBLAS_NUM_THREADS","MKL_NUM_THREADS"):4 os.environ[_k]="1"5import sys, math, json, csv, time, datetime, hashlib, gzip, argparse6import numpy as np7P=Path(__file__).resolve().parent8OLD=Path("D:/Photometer_research_20261022/research")9sys.path.insert(0,str(OLD/"cloud_transport/private_readers"))10from scipy.stats import binom11from scipy.optimize import minimize_scalar12T=0.00313T_PRE=0.01014RISES_US=[30,80,200]15COUNTS=[100,1000,10000]16BACKGROUNDS=[50.,1e4,1e6]17SAMPLE_KHZ=[10,50,100]18TRIALS=100019BASE_SEED=2026102208100120ALPHA_PER_TEST=2.866515718791933e-721MARGIN=math.log(19.)22ONSET_GRID=np.arange(-150.,150.0001,5.)*1e-623START_UTC=datetime.datetime.now(datetime.timezone.utc).isoformat()2425def primitive(t,d):26 z=np.clip(np.asarray(t,dtype=float),0,T)27 return np.where(z<d,z*z/(2*d),z-d/2)2829def source_integrals(edges,onsets,rise_us):30 return np.diff(primitive(edges[None,:]-np.atleast_1d(onsets)[:,None],rise_us*1e-6/.8),axis=1)3132def setup(sample_khz):33 dt=1/(1000*sample_khz)34 # The full fixed candidate-event gate is [-0.2,3.2] ms.35 edges=np.linspace(-.0002,.0032,round(.0034/dt)+1)36 centers=(edges[:-1]+edges[1:])/237 early=np.where(centers<.0006-1e-12)[0]38 late=np.where((centers>.001-1e-12)&(centers<.0025-1e-12))[0]39 assert len(early)==round(.0008/dt) and len(late)==round(.0015/dt)40 durations=np.r_[np.full(len(early),dt),len(late)*dt,T_PRE]41 models=[]42 for rise in RISES_US:43 q_all=source_integrals(edges,ONSET_GRID,rise)44 for k,onset in enumerate(ONSET_GRID):45 q=np.r_[q_all[k,early],q_all[k,late].sum(),0.]46 f=q/durations47 full=f>=1-1e-1148 zero=f<=1e-1149 active=~(full|zero)50 Q=float(q.sum()); D=float(durations.sum())51 K=D/Q52 a=f[active]*K-1.53 models.append(dict(rise_us=rise,onset_s=float(onset),q=q,f=f,Q=Q,D=D,54 full=full,zero=zero,active=active,a=a,K=K))55 return dict(dt=dt,edges=edges,early=early,late=late,durations=durations,models=models)5657def fit_one(counts,model,max_iter=40,tol=2e-10):58 """Exact joint Poisson profile over nonnegative peak rate and background.59 Normalized source/background weights sum to1. At the MLE their total60 expectation equals observed N, leaving a concave 1-D mixture fraction theta.61 """62 n=counts.sum(axis=1).astype(float)63 n0=counts[:,model["zero"]].sum(axis=1).astype(float)64 np_=counts[:,model["full"]].sum(axis=1).astype(float)65 nr=counts[:,model["active"]].astype(float)66 a=model["a"]; c=model["K"]-1.67 g0=-n0+np_*c+(nr*a).sum(axis=1)68 theta=np.zeros(len(n))69 signal=g0>1e-1070 if np.any(signal):71 ix=np.flatnonzero(signal)72 N0=n0[ix]; NP=np_[ix]; NR=nr[ix]73 # Boundary theta=1 is possible only with no counts where q=0.74 g1=np.full(len(ix),-np.inf)75 support=N0==076 if np.any(support):77 g1[support]=NP[support]*c/(1+c)+(NR[support]*a/(1+a)).sum(axis=1)78 upper=g1>=079 theta[ix[upper]]=1.80 inner=~upper81 if np.any(inner):82 jx=ix[inner]; N0=N0[inner]; NP=NP[inner]; NR=NR[inner]83 lo=np.zeros(len(jx)); hi=np.ones(len(jx))84 # One Newton step from theta=0 is a stable nonnegative initializer.85 curvature=N0+NP*c*c+(NR*a*a).sum(axis=1)86 x=np.clip(g0[jx]/np.maximum(curvature,1e-300),1e-12,.99)87 done=np.zeros(len(jx),dtype=bool)88 for iteration in range(max_iter):89 den=1+x[:,None]*a90 g=-N0/(1-x)+NP*c/(1+x*c)+(NR*a/den).sum(axis=1)91 h=-N0/(1-x)**2-NP*c*c/(1+x*c)**2-(NR*a*a/(den*den)).sum(axis=1)92 lo=np.where(g>0,x,lo); hi=np.where(g<=0,x,hi)93 proposal=x-g/np.minimum(h,-1e-300)94 bad=(~np.isfinite(proposal))|(proposal<=lo)|(proposal>=hi)95 proposal=np.where(bad,(lo+hi)/2,proposal)96 step=np.abs(proposal-x)97 converged=(np.abs(g)/(n[jx]+1)<tol)|(hi-lo<tol)|(step<tol)98 done|=converged99 x=np.where(done,x,proposal)100 if np.all(done):101 break102 theta[jx]=x103 # Zero count bins contribute0 even if fitted probability is0.104 v0=np.zeros(len(n)); nz=n0>0105 v0[nz]=n0[nz]*np.log1p(-theta[nz])106 llr=v0+np_*np.log1p(theta*c)+(nr*np.log1p(theta[:,None]*a)).sum(axis=1)107 peak=theta*n/model["Q"]108 b=(1-theta)*n/model["D"]109 return llr,theta,peak,b110111def profile(counts,design):112 ntrial=len(counts)113 profiles=np.full((ntrial,3),-np.inf)114 onsets=np.full((ntrial,3),np.nan)115 peaks=np.full((ntrial,3),np.nan)116 bs=np.full((ntrial,3),np.nan)117 for model in design["models"]:118 k=RISES_US.index(model["rise_us"])119 llr,theta,peak,b=fit_one(counts,model)120 update=llr>profiles[:,k]121 profiles[update,k]=llr[update]122 onsets[update,k]=np.where(theta[update]>0,model["onset_s"],np.nan)123 peaks[update,k]=peak[update];bs[update,k]=b[update]124 return profiles,onsets,peaks,bs125126def wilson(k,n,z=1.959963984540054):127 if n==0:128 return [None,None]129 p=k/n; den=1+z*z/n130 middle=(p+z*z/(2*n))/den131 half=z*math.sqrt(p*(1-p)/n+z*z/(4*n*n))/den132 return [max(0.,middle-half),min(1.,middle+half)]133134def diagnostic_metric(k,n):135 return {"count":int(k),"denominator":int(n),"fraction":float(k/n) if n else None,136 "individual_95pct_Wilson_interval":wilson(k,n)}137138def draw_case(design,n_eq,b,rise_us,seed,ntrial=TRIALS,signal=True):139 rng=np.random.default_rng(seed)140 onset=rng.uniform(-100e-6,100e-6,size=ntrial)141 q=source_integrals(design["edges"],onset,rise_us)142 source=(n_eq/T)*q if signal else np.zeros_like(q)143 means=b*design["dt"]+source144 full=rng.poisson(means)145 pre=rng.poisson(b*T_PRE,size=ntrial)146 fitcounts=np.column_stack([full[:,design["early"]],full[:,design["late"]].sum(axis=1),pre])147 # Exact conditional on/off test of equal unknown constant background rates.148 total=full.sum(axis=1)149 p=binom.sf(total-1,total+pre,.0034/(.0034+T_PRE))150 detected=p<=ALPHA_PER_TEST151 return dict(onset=onset,full=full,pre=pre,fitcounts=fitcounts,total=total,152 pvalue=p,detected=detected,source=source)153154def analyze_case(design,n_eq,b,rise_us,seed,ntrial=TRIALS):155 draw=draw_case(design,n_eq,b,rise_us,seed,ntrial)156 scores,onsets,peaks,bs=profile(draw["fitcounts"],design)157 n=len(scores); ar=np.arange(n)158 best=scores.argmax(axis=1)159 sorted_scores=np.sort(scores,axis=1)160 margin3=sorted_scores[:,-1]-sorted_scores[:,-2]161 accept3=draw["detected"]&(margin3>=MARGIN)162 fastscore=scores[:,:2].max(axis=1)163 slow=scores[:,2]164 bestgroup=(slow>fastscore).astype(int) #0 fast morphology,1 slow morphology165 margin2=np.abs(slow-fastscore)166 accept2=draw["detected"]&(margin2>=MARGIN)167 truth=RISES_US.index(rise_us)168 truthgroup=0 if truth<2 else 1169 pred3=np.where(~draw["detected"],4,np.where(~accept3,3,best))170 pred2=np.where(~draw["detected"],3,np.where(~accept2,2,bestgroup))171 actual=n_eq*(1-rise_us*1e-6/.8/(2*T))172 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]}173 row={"N_plateau_equivalent_3ms_PE":n_eq,"background_cps":b,"sample_rate_khz":round(1/design["dt"]/1000),174 "truth_receiver_rise_10_90_us":rise_us,"trials":ntrial,"seed":seed,175 "actual_total_source_PE":actual,"peak_source_cps":n_eq/T,**early,176 "expected_background_first_200us_PE":b*.0002,177 "expected_source_in_gate_PE_mean":float(draw["source"].sum(axis=1).mean()),178 "expected_background_in_gate_PE":b*.0034,179 "expected_background_prewindow_PE":b*T_PRE,180 "peak_total_rate_times_17ns":(b+n_eq/T)*17e-9,181 "count_loss_over_1pct_stationary_mean_flag":(b+n_eq/T)*17e-9>.01/ .99,182 "detection":diagnostic_metric(draw["detected"].sum(),ntrial),183 "three_class":{"columns":["rise30","rise80","rise200","rejected_detected","undetected"],184 "confusion_counts":np.bincount(pred3,minlength=5).tolist(),185 "accepted":diagnostic_metric(accept3.sum(),ntrial),186 "wrong_among_accepted":diagnostic_metric((accept3&(best!=truth)).sum(),accept3.sum()),187 "rejected_among_detected":diagnostic_metric((draw["detected"]&~accept3).sum(),draw["detected"].sum())},188 "fast_slow":{"columns":["fast30or80","slow200","rejected_detected","undetected"],189 "confusion_counts":np.bincount(pred2,minlength=4).tolist(),190 "accepted":diagnostic_metric(accept2.sum(),ntrial),191 "wrong_among_accepted":diagnostic_metric((accept2&(bestgroup!=truthgroup)).sum(),accept2.sum()),192 "rejected_among_detected":diagnostic_metric((draw["detected"]&~accept2).sum(),draw["detected"].sum())},193 "oracle_same_data_three_class":{"error":diagnostic_metric((best!=truth).sum(),ntrial),194 "qualification":"Name retained only for compatibility; this is unthresholded practical profile choice, NOT an oracle."},195 "fitted_diagnostics":{"median_peak_rate_cps":float(np.median(peaks[ar,best])),196 "median_background_cps":float(np.median(bs[ar,best])),197 "median_onset_error_us":float(np.nanmedian((onsets[ar,best]-draw["onset"])*1e6)),198 "median_abs_onset_error_us":float(np.nanmedian(abs(onsets[ar,best]-draw["onset"])*1e6)),199 "onset_grid_boundary_fraction":float(np.mean(abs(onsets[ar,best])>=149.99e-6)),200 "median_three_class_score_margin":float(np.median(margin3)),201 "median_fast_slow_score_margin":float(np.median(margin2))}}202 del row["oracle_same_data_three_class"]203 row["unthresholded_practical_three_class_error"]=diagnostic_metric((best!=truth).sum(),ntrial)204 trials={"onset_s":draw["onset"],"pre_count":draw["pre"],"gate_count":draw["total"],205 "detection_pvalue":draw["pvalue"],"detected":draw["detected"],206 "scores":scores,"best_rise_index":best,"margin3":margin3,"accept3":accept3,207 "best_group":bestgroup,"margin2":margin2,"accept2":accept2,208 "fitted_onset_s":onsets[ar,best],"fitted_peak_cps":peaks[ar,best],209 "fitted_background_cps":bs[ar,best]}210 return row,trials211212def check_math():213 checks=[]214 for khz in SAMPLE_KHZ:215 design=setup(khz)216 for rise in RISES_US:217 onsets=np.array([-.0001,-.000037,0,.000099])218 q=source_integrals(design["edges"],onsets,rise)219 expected=T-rise*1e-6/.8/2220 assert np.max(abs(q.sum(axis=1)-expected))<1e-14221 checks.append({"test":"full_gate_exact_source_integral","sample_rate_khz":khz,"passed":True})222 design=setup(50)223 rng=np.random.default_rng(5821031)224 test_counts=np.vstack([np.zeros(len(design["durations"])),225 rng.poisson(np.r_[np.full(40,3.),200.,500.],size=(8,42))])226 # Independent scalar bounded optimization uses the full normalized mixture.227 errors=[]228 for model in design["models"][::13]:229 ll,theta,_,_=fit_one(test_counts,model)230 h=model["q"]/model["Q"]; w=design["durations"]/model["D"]231 for j,n in enumerate(test_counts):232 a=h/w-1233 def objective(x):234 valid=n>0235 return -float((n[valid]*np.log1p(x*a[valid])).sum())236 opt=minimize_scalar(objective,bounds=(0,1-1e-13),method="bounded",237 options={"xatol":1e-12})238 best=min(0.,opt.fun,objective(1-1e-13))239 errors.append(abs(ll[j]+best))240 assert max(errors)<2e-6241 checks.append({"test":"profile_MLE_vs_independent_scalar_optimizer","comparisons":len(errors),242 "max_absolute_loglikelihood_difference":max(errors),"passed":True})243 # Same fitted parameters imply expected count sum equals actual total.244 model=design["models"][60]245 ll,theta,A,b=fit_one(test_counts,model)246 lam=A[:,None]*model["q"]+b[:,None]*design["durations"]247 assert np.max(abs(lam.sum(axis=1)-test_counts.sum(axis=1)))<1e-8248 assert np.all(ll>=-1e-8)249 checks.append({"test":"joint_MLE_count_scale_identity_and_LLR_nonnegative","passed":True})250 # Time shifts wholly inside fit windows leave late plateau flat.251 assert all(np.allclose(m["f"][-2],1) and m["f"][-1]==0 for m in design["models"])252 checks.append({"test":"late_plateau_and_prewindow_source_exclusion","passed":True})253 return checks254255def csv_rows(rows):256 for r in rows:257 yield {k:r[k] for k in ["N_plateau_equivalent_3ms_PE","background_cps","sample_rate_khz",258 "truth_receiver_rise_10_90_us","trials","seed","actual_total_source_PE",259 "expected_source_first_30us_PE","expected_source_first_80us_PE",260 "expected_source_first_200us_PE","expected_background_first_200us_PE",261 "peak_total_rate_times_17ns"]} | {262 "detected":r["detection"]["count"],"detection_efficiency":r["detection"]["fraction"],263 "accepted3":r["three_class"]["accepted"]["count"],264 "wrong3":r["three_class"]["wrong_among_accepted"]["count"],265 "error_among_accepted3":r["three_class"]["wrong_among_accepted"]["fraction"],266 "accepted_fast_slow":r["fast_slow"]["accepted"]["count"],267 "wrong_fast_slow":r["fast_slow"]["wrong_among_accepted"]["count"],268 "error_among_accepted_fast_slow":r["fast_slow"]["wrong_among_accepted"]["fraction"],269 "rejected_detected_fast_slow":r["fast_slow"]["confusion_counts"][2],270 "undetected":r["fast_slow"]["confusion_counts"][3]}271272def save_csv(path,rows):273 rows=list(rows)274 with path.open("w",newline="",encoding="utf-8") as f:275 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)276277def run(mode):278 P.mkdir(parents=True,exist_ok=True)279 start=time.perf_counter()280 checks=check_math()281 print(json.dumps({"QA":checks}),flush=True)282 if mode=="qa":283 (P/"math_QA.json").write_text(json.dumps({"checks":checks,"passed":True},indent=2))284 return285 if mode=="benchmark":286 d=setup(100)287 r,tr=analyze_case(d,1000,1e6,80,BASE_SEED,1000)288 print(json.dumps({"runtime_s":time.perf_counter()-start,"example":r}),flush=True)289 return290 rows=[];case=0291 (P/"synthetic_trial_summaries").mkdir(exist_ok=True)292 for khz in SAMPLE_KHZ:293 d=setup(khz)294 for n_eq in COUNTS:295 for b in BACKGROUNDS:296 for rise in RISES_US:297 seed=BASE_SEED+case298 t0=time.perf_counter()299 r,tr=analyze_case(d,n_eq,b,rise,seed)300 r["case_id"]=case;r["runtime_s"]=time.perf_counter()-t0301 rows.append(r)302 np.savez_compressed(P/"synthetic_trial_summaries"/f"case_{case:03d}.npz",**tr)303 print(json.dumps({"case_id":case,"khz":khz,"N":n_eq,"b":b,"rise":rise,304 "runtime_s":r["runtime_s"],"detected":r["detection"]["count"],305 "accepted_binary":r["fast_slow"]["accepted"]["count"],306 "wrong_binary":r["fast_slow"]["wrong_among_accepted"]["count"]}),flush=True)307 case+=1308 noise=[]309 for khz in SAMPLE_KHZ:310 d=setup(khz)311 for j,b in enumerate(BACKGROUNDS):312 seed=BASE_SEED+10000+100*khz+j313 dr=draw_case(d,0,b,80,seed,1000,signal=False)314 noise.append({"sample_rate_khz":khz,"background_cps":b,"seed":seed,315 "trials":1000,"false_detections":int(dr["detected"].sum()),316 "individual_95pct_Wilson_interval":wilson(dr["detected"].sum(),1000),317 "qualification":"1000 noise trials cannot empirically establish a 2.87e-7 tail. Exact conditional test supplies the per-event-gate level."})318 save_csv(P/"practical_morphology_scenarios.csv",csv_rows(rows))319 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])320 out={"schema_version":"1.0","status":"EXECUTED_CONDITIONAL_IDEAL_POISSON_PRACTICAL_MORPHOLOGY_PROTOTYPE",321 "start_utc":START_UTC,"finished_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),322 "runtime_s":time.perf_counter()-start,"seed_base":BASE_SEED,323 "scenario_count":len(rows),"trials_each_scenario":TRIALS,"total_signal_trials":len(rows)*TRIALS,324 "parameters":{"receiver_rises_10_90_us":RISES_US,"plateau_equivalent_PE_in_3ms":COUNTS,325 "background_cps":BACKGROUNDS,"sample_rate_khz":SAMPLE_KHZ,326 "pulse_end_after_onset_ms":3.,"prewindow_ms":10.,327 "prewindow_position_ms":[-10.5,-.5],"candidate_gate_ms":[-.2,3.2],328 "morphology_early_window_ms":[-.2,.6],"late_plateau_window_ms":[1.,2.5],329 "true_onset_distribution_us":[-100,100],"onset_grid_us":[-150,150,5],330 "alpha_per_predefined_gate":ALPHA_PER_TEST,331 "acceptance_profile_loglikelihood_margin":MARGIN,332 "deadtime_diagnostic_ns":17.},333 "algorithm":{"label_definition":"Linear receiver-input expected PE-rate morphology after any assumed transport, not intrinsic discharge history or cloud depth.",334 "waveform":"A linear0-to100% ramp lasting rise_10_90/0.8, then unit plateau, with total pulse ending3ms after onset.",335 "amplitude":"Peak rate A=N_plateau_equivalent/3ms. Actual integrated source counts=A*(3ms-ramp_duration/2); source and total/background counts are separate.",336 "sampling":"Independent Poisson hardware-bin counts from exact primitive differences. Onset uniform[-100,100]us spans integer phase cycles at every sample rate.",337 "background":"Constant unknown rate. One explicit10ms background-only count plus the fit-window counts enter a joint Poisson profile likelihood.",338 "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].",339 "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.",340 "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.",341 "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.",342 "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.",343 "response":"No electronic or arrival-delay filter. Analog post-counting RC response would correlate samples and invalidate this independent-Poisson likelihood.",344 "MC_intervals":"Individual two-sided95% Wilson score intervals for reported binomial rates. Not simultaneous over all scenarios, not uncertainty in the true event population."},345 "math_QA":checks,"scenarios":rows,"noise_only_controls":noise,346 "limitations":[347 "Conditional synthetic morphology experiment, not validated flight performance or cloud-depth inference.",348 "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.",349 "Unknown onset is restricted to an already available200us coarse bracket; identifying that bracket in a blind stream is outside scope.",350 "The fixed3.4ms detection gate covers each simulated pulse; constructing such gates blindly and overlapping-window search trials remain operational inputs.",351 "The nonnegative nuisance fit can favor flexible fast templates; reported Monte Carlo confusion counts, not the score margin, quantify performance under this generator.",352 "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.",353 "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.",354 "PE counts include an unspecified end-to-end photon-detection efficiency. Aperture photons require a separately verified spectral/angular efficiency.",355 "No molecular source yield, source brightness distribution, population completeness or optical duty/exposure model is supplied. Useful mission event rate is NOT COMPUTABLE.",356 "A real morphology classifier requires end-to-end hardware waveforms and separately labeled optical/radio control populations.",357 "Fitted-grid errors and Monte Carlo intervals do not cover systematic model uncertainty. Samples are10/20/100us integrations; sampling rate is not analog bandwidth."] }358 (P/"practical_morphology_results.json").write_text(json.dumps(out,indent=2),encoding="utf-8")359 print(json.dumps({"finished":True,"runtime_s":out["runtime_s"],"scenarios":len(rows)}),flush=True)360361if __name__=="__main__":362 parser=argparse.ArgumentParser();parser.add_argument("--mode",choices=["qa","benchmark","full"],default="full")363 run(parser.parse_args().mode)364365