Review index · Download raw source
reference/cloud/cloudscat_control.py
1"""Matched homogeneous-HG cross-code control using retained CloudScat calculations."""2from pathlib import Path3import json,csv,math,datetime,hashlib4import numpy as np5H=Path(__file__).resolve().parent;O=H/"validation"/"cloudscat_control";O.mkdir(exist_ok=True)6# Separate controlled extension: default/main cylinder engine remains conservative.7s=(H/"cylinder_mc.py").read_text(encoding="utf-8")8s=s.replace('importance_cone_deg=None):','importance_cone_deg=None,albedo=1.):')9s=s.replace('first_hist=np.zeros(nbins+1)','first_hist=np.zeros(nbins+1)\n absorbed_weight=0.')10s=s.replace('toward=-pos.copy();toward[:,2]+=geo.observer_z_m','absorbed_weight+=float(w.sum())*(1-albedo)\n w=w*albedo\n toward=-pos.copy();toward[:,2]+=geo.observer_z_m')11s=s.replace('residual=launch_sum-escaped_weight-truncated_weight','residual=launch_sum-escaped_weight-truncated_weight-absorbed_weight')12s=s.replace('"weighted_launch_mean":launch_sum/n,','"weighted_launch_mean":launch_sum/n,"absorbed_weight_fraction":absorbed_weight/n,"albedo":albedo,')13(O/"matched_control_engine.py").write_text(s,encoding="utf-8")14import sys15sys.path.insert(0,str(O));from matched_control_engine import Phase,Geometry,simulate,sha16import timing_metrics as tm17OLD=Path("C:/Users/user/Documents/Codex/2026-10-14/given-a-source-at-15-18/outputs/CloudScat/calculations")18def dump(p,x):p.write_text(json.dumps(x,indent=2,allow_nan=False),encoding="utf-8")19def sem(x):return float(np.std(x,ddof=1)/np.sqrt(len(x)))20def shape_metrics(mass,direct):21 # Place scatteredbincenters at1,3,..us; directatom at0 on a1us grid.22 h=np.zeros(12001);h[0]=direct;h[1::2]=mass23 t,y=tm.propagate_waveform(h/1e-6,1e-6,30.,None)24 return tm.waveform_metrics(t,y)25def main():26 start=datetime.datetime.now(datetime.timezone.utc).isoformat()27 rows=list(csv.DictReader((OLD/"results"/"batch_summary.csv").open(encoding="utf-8")))28 result=[]29 for oldid in ["noray_3","noray_10"]:30 oldrows=[x for x in rows if x["case_id"]==oldid];r=oldrows[0]31 source=float(r["zsrc_m"]);base=float(r["zbottom_m"]);top=float(r["ztop_m"]);radius=float(r["cloudradius_m"])32 omega=float(r["albedo"]);g=float(r["g"]);tau=float(r["tau_cloud_total"]);ph=Phase("hg",g=g)33 geo=Geometry(base_m=base,top_m=top,radius_m=radius)34 vac=1/(4*math.pi*(550000-source)**2);direct=float(r["direct_fluence_m-2"])/vac35 ref=[];ref_total=[];refs=[]36 for rr in oldrows:37 p=OLD/"results"/f'{oldid}_batch{rr["batch"]}.csv'38 a=np.genfromtxt(p,delimiter=",",names=True)39 mass=a[a.dtype.names[1]]*2e-6/vac40 ref_total.append(float(mass.sum()));mass[0]-=direct41 if mass[0]<-1e-10:raise ValueError("Reference direct subtraction failed")42 mass=np.maximum(mass,0);ref.append(mass)43 refs.append({"file":p.name,"sha256":sha(p),"seed":int(rr["seed"]),"n_packets":int(rr["N"])})44 records=[];new=[]45 for b in range(6):46 p=O/f"{oldid}_new_batch{b}.json";q=O/f"{oldid}_new_batch{b}.npz"47 if p.exists()and q.exists():info=json.loads(p.read_text(encoding="utf-8"));h=np.load(q)["mass"]48 else:49 info,h=simulate(ph,source,tau,50000,93200100+(0 if oldid=="noray_3" else 100)+b,geo=geo,albedo=omega)50 dump(p,info);np.savez_compressed(q,mass=h,dt_s=2e-6)51 records.append(info);new.append(h)52 new=np.asarray(new);ref=np.asarray(ref)53 np.savez_compressed(O/f"{oldid}_reference_kernels.npz",mass=ref,dt_s=2e-6,direct=direct)54 new_total=np.array([x["total_relative_to_vacuum"]-x["response_overflow_relative_to_vacuum"]for x in records])55 delta=float(new_total.mean()-np.mean(ref_total));sigma=math.hypot(sem(new_total),sem(ref_total))56 row={"reference_case":oldid,"physical_parameters":{"source_m":source,"base_m":base,"top_m":top,"radius_m":radius,"total_tau":tau,"overlying_tau":float(r["tau_over"]),"g":g,"SSA":omega,"observer_m":550000,"full_fov_deg":60,"dt_us":2,"record_ms":12},57 "reference_physical_method":"CloudScat HG with runtime plane-aperture observer, exponential cloud optical path quadrature, vacuum atmosphere. Original package/source and configuration retained read-only.",58 "new_method":"Independent exact cylinder intersections and next-event estimator; HG and identical SSA via implicit absorption weights in this control only.",59 "reference_files":refs,"reference_batch_gain":ref_total,"reference_gain_mean":float(np.mean(ref_total)),"reference_gain_batch_SE":sem(ref_total),60 "new_batch_gain":new_total.tolist(),"new_gain_mean":float(new_total.mean()),"new_gain_batch_SE":sem(new_total),61 "difference_new_minus_reference":delta,"combined_MC_standard_error":sigma,"difference_in_combined_SE":delta/sigma,"within_3SE":abs(delta)<3*sigma,62 "new_escaped_plus_absorbed_weight_balance_max_abs":max(abs(x["energy_balance_residual_per_packet"])for x in records),63 "new_truncated_packets":sum(x["truncated_packet_count"]for x in records),64 "new_overflow_mean":float(np.mean([x["response_overflow_relative_to_vacuum"]for x in records])),65 "reference_30us_source_ideal_metrics":shape_metrics(ref.mean(axis=0),direct),66 "new_30us_source_ideal_metrics":shape_metrics(new.mean(axis=0)[:-1],direct),67 "qualification":"Independent-code comparison for these ideal HG cases only. It does not validate tabulated GHM, real ice, atmosphere or source-depth classification."}68 result.append(row);print(json.dumps({k:row[k]for k in ["reference_case","reference_gain_mean","new_gain_mean","difference_in_combined_SE","within_3SE"]}),flush=True)69 out={"started_utc":start,"ended_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),"cases":result,70 "original_model_sha256":sha(OLD/"photometer_model.jl"),"original_runner_sha256":sha(OLD/"run_grid.jl"),"original_CloudScat_source_sha256":sha(OLD.parent/"source"/"src"/"CloudScat.jl"),"control_engine_sha256":sha(O/"matched_control_engine.py"),71 "status":"EXECUTED_INDEPENDENT_CODE_CONTROL_FOR_MATCHED_IDEAL_HG_ONLY"}72 dump(O/"comparison.json",out)73if __name__=="__main__":main()74