Review index · Download raw source
reference/cloud/validate_model.py
1"""Bounded validation and convergence diagnostics; no observational validation."""2from pathlib import Path3import json,math,time,datetime4import numpy as np5from cylinder_mc import Phase,Geometry,simulate,intervals,segment_inside,load_phases,sha6import timing_metrics as tm7H=Path(__file__).resolve().parent;O=H/"validation";O.mkdir(exist_ok=True)8LUT=Path("D:/Photometer_research_20261022/research/cloud_transport/HYGEOS_ghm_337_bracketing_RAW.private.npz")9def dump(p,x):p.write_text(json.dumps(x,indent=2,allow_nan=False),encoding="utf-8")10def sem(x):return float(np.std(x,ddof=1)/np.sqrt(len(x)))11def run(label,phase,source,tau,n,batches=6,geo=None,importance=.5,dt=2e-6,cap=20000):12 d=O/label;d.mkdir(exist_ok=True);records=[];kernels=[]13 for b in range(batches):14 p=d/f"batch{b}.json";q=d/f"batch{b}.npz"15 if p.exists()and q.exists():16 info=json.loads(p.read_text(encoding="utf-8"));h=np.load(q)["mass"]17 else:18 info,h=simulate(phase,source,tau,n,910000+sum(ord(x)for x in label)*100+b,geo=geo,importance=importance,dt_s=dt,max_scatter=cap)19 dump(p,info);np.savez_compressed(q,mass=h,dt_s=dt)20 records.append(info);kernels.append(h)21 gains=[x["total_relative_to_vacuum"]for x in records]22 r={"label":label,"n_packets_per_batch":n,"batches":batches,"mean_gain":float(np.mean(gains)),"gain_batch_standard_error":sem(gains)if batches>1 else 0.,23 "mean_first_scatter_gain":float(np.mean([x["first_scatter_relative_to_vacuum"]for x in records])),24 "first_scatter_batch_standard_error":sem([x["first_scatter_relative_to_vacuum"]for x in records])if batches>1 else 0.,25 "truncated_packets":sum(x["truncated_packet_count"]for x in records),26 "elapsed_s":sum(x["elapsed_s"]for x in records)}27 print(json.dumps(r),flush=True);return r,records,np.asarray(kernels)28def single_scatter_quad(source,tau,order):29 """Axisymmetric angular/path quadrature; deterministic first-collision integral."""30 geo=Geometry();beta=tau/geo.thickness;r0=geo.observer_z_m-source31 bounds=[-1.,-(source-geo.base_m)/math.hypot(geo.radius_m,source-geo.base_m),0.,(geo.top_m-source)/math.hypot(geo.radius_m,geo.top_m-source),1.]32 q,w=np.polynomial.legendre.leggauss(order);total=0.33 for a,b in zip(bounds[:-1],bounds[1:]):34 mu=(a+b)/2+(b-a)*q/2;wm=w*(b-a)/235 u=np.column_stack((np.sqrt(1-mu*mu),np.zeros(len(mu)),mu));p=np.tile([0.,0.,source],(len(mu),1))36 lo,hi,hit=intervals(p,u,geo)37 s=(lo[:,None]+hi[:,None])/2+(hi-lo)[:,None]*q/2;ws=(hi-lo)[:,None]*w/238 at=p[:,None,:]+s[:,:,None]*u[:,None,:];flat=at.reshape(-1,3)39 toward=-flat;toward[:,2]+=geo.observer_z_m;rr=np.linalg.norm(toward,axis=1);uu=toward/rr[:,None]40 length=segment_inside(flat,uu,rr,geo).reshape(s.shape);rr=rr.reshape(s.shape);cos=uu[:,2].reshape(s.shape)41 f=.5*beta*np.exp(-beta*(s-lo[:,None]))*r0*r0*cos/(rr*rr)*np.exp(-beta*length)42 total+=float(np.sum(wm[:,None]*ws*f))43 return total44def main():45 started=datetime.datetime.now(datetime.timezone.utc).isoformat();t0=time.perf_counter()46 phases,_=load_phases(LUT);out={"started_utc":started,"scope":"Own analytic/numerical tests and sensitivity checks, not an independent-code or real-cloud validation.","checks":[]}47 tests=[]48 for source in [15000.,16000.,17000.,18000.]:49 i,h=simulate(phases["hg"],source,3.,1000,447+int(source),importance=0)50 target=math.exp(-3*max(0,17000-source)/3000)51 tests.append({"source_m":source,"computed_direct":i["direct_relative_to_vacuum"],"expected_direct":target,"pass":abs(i["direct_relative_to_vacuum"]-target)<1e-14})52 vac=[]53 for obs in [550000.,800000.]:54 i,h=simulate(phases["hg"],15000.,0.,1000,int(obs),geo=Geometry(observer_z_m=obs),importance=0)55 vac.append(i)56 ratio=vac[0]["vacuum_transfer_m_inverse2"]/vac[1]["vacuum_transfer_m_inverse2"]57 out["direct_tests"]=tests;out["vacuum_inverse_square"]={"ratio":ratio,"expected":((800000-15000)/(550000-15000))**2,"zero_scattered":all(x["scattered_relative_to_vacuum"]==0 for x in vac),"zero_truncation":all(x["truncated_packet_count"]==0 for x in vac)}58 quad=[{"order":n,"scattered_first_order_gain":single_scatter_quad(15000.,.3,n)}for n in [24,48,96,192]]59 control,_,_=run("isotropic_single_scatter",Phase("hg",g=0.),15000.,.3,30000,importance=0)60 err=control["mean_first_scatter_gain"]-quad[-1]["scattered_first_order_gain"]61 out["single_scatter_integral"]={"quadrature":quad,"MC":control,"difference":err,"difference_in_MC_standard_errors":err/control["first_scatter_batch_standard_error"],"pass_5SE":abs(err)<5*control["first_scatter_batch_standard_error"],"qualification":"Deterministic equation-based reference; no independent RT code."}62 out["representative_controls"]=[]63 for name,phase in phases.items():64 for imp in [0.,.5]:65 r,_,_=run(f"{name}_importance{imp}",phase,15000.,10.,20000,importance=imp)66 out["representative_controls"].append(r)67 for radius in [15000.,60000.]:68 r,_,_=run(f"{name}_radius{int(radius)}",phase,15000.,10.,20000,geo=Geometry(radius_m=radius))69 out["representative_controls"].append(r)70 # Higher packet count is a diagnostic; it does not silently replace primary grid data.71 for name in ["ghm","hg"]:72 r,_,_=run(f"{name}_larger_packets",phases[name],15000.,10.,200000)73 out["representative_controls"].append(r)74 # Same physical packet histories,1us vs2us binning; coarse4us by exact mass regrouping.75 r,recs,h=run("ghm_timebin_1us",phases["ghm"],16000.,10.,30000,dt=1e-6)76 mean=h.mean(axis=0);direct=recs[0]["direct_relative_to_vacuum"];tim=[]77 for width in [1,2,4]:78 masses=np.zeros(int(.012/(.5e-6))+1);masses[0]=direct79 rebin=mean[:-1].reshape(-1,width).sum(axis=1);idx=((np.arange(len(rebin))+.5)*width/.5).astype(int)80 masses[idx]=rebin81 for rise in [30.,80.]:82 tt,yy=tm.propagate_waveform(masses/.5e-6,.5e-6,rise,None)83 z=tm.waveform_metrics(tt,yy)84 tim.append({"scatter_bin_width_us":width,"intrinsic_rise_10_90_us":rise,"rise_10_90_us":z["rise_10_90_us"],"onset_to_peak_us":z["onset_to_peak_us"]})85 out["time_bin_sensitivity"]={"records":tim,"qualification":"Identical trajectories with native1us bins regrouped to2/4us. This isolates bin placement, not packet uncertainty."}86 truncated,_,_=run("deliberate_truncation_cap10",phases["hg"],15000.,45.,2000,batches=1,cap=10)87 out["negative_control"]={"deliberate_truncation_cap10":truncated,"expected":"Nonzero truncated packets must remain visible; this control is excluded from production."}88 out["ended_utc"]=datetime.datetime.now(datetime.timezone.utc).isoformat();out["wall_s"]=time.perf_counter()-t089 dump(O/"validation_results.json",out);print(json.dumps({"validation_status":"DONE","seconds":out["wall_s"]}),flush=True)90if __name__=="__main__":main()91