Review index · Download raw source

reference/cloud/cone_control.py

1"""Bounded source-direction importance control; no change to physical phase."""2from pathlib import Path3import json,math,datetime,time4import numpy as np5from cylinder_mc import load_phases,launch,simulate,sha6H=Path(__file__).resolve().parent;O=H/"validation"/"cone_control";O.mkdir(exist_ok=True)7def dump(p,x):p.write_text(json.dumps(x,indent=2,allow_nan=False),encoding="utf-8")8def main():9    phases,_=load_phases("D:/Photometer_research_20261022/research/cloud_transport/HYGEOS_ghm_337_bracketing_RAW.private.npz");ph=phases["ghm"]10    start=datetime.datetime.now(datetime.timezone.utc).isoformat();t0=time.perf_counter()11    u,w=launch(np.random.default_rng(991144),1000000,ph,.5,1.)12    checks=[]13    for label,v,target in [("weighted_normalization",w,1.),("weighted_mu",w*u[:,2],0.),("weighted_P2",w*(3*u[:,2]**2-1)/2,0.)]:14        mean=float(v.mean());se=float(v.std(ddof=1)/np.sqrt(len(v)))15        checks.append({"metric":label,"mean":mean,"expected":target,"standard_error":se,"z":(mean-target)/se,"pass_5SE":abs(mean-target)<5*se})16    records=[]17    for b in range(6):18        p=O/f"batch{b}.json";q=O/f"batch{b}.npz"19        if p.exists()and q.exists():info=json.loads(p.read_text(encoding="utf-8"))20        else:21            info,h=simulate(ph,15000.,10.,200000,92100400+b,importance=.5,importance_cone_deg=1.)22            dump(p,info);np.savez_compressed(q,mass=h,dt_s=2e-6)23        records.append(info);print(json.dumps({"batch":b,"gain":info["total_relative_to_vacuum"],"max_packet":info["max_single_packet_scattered_gain"]}),flush=True)24    gains=np.array([x["total_relative_to_vacuum"]for x in records])25    out={"description":"Physical isotropic source. Proposal q=(1-a)/(4pi)+a/Omega_cone in1degree cap about zenith; a=.5. Exact likelihood weight (1/4pi)/q. All collision directions and estimator use unchanged GHM phase.",26      "weight_bound_max":2.,"cone_full_solid_angle_sr":2*math.pi*(1-math.cos(math.radians(1.))),27      "source_sampling_checks":checks,"sample_check_packets":1000000,"geometry":"Same as primarygrid;source15km,totaltau10.",28      "mean_gain":float(gains.mean()),"batch_standard_error":float(gains.std(ddof=1)/np.sqrt(6)),"batches":6,"packets_per_batch":200000,29      "batch_gains":gains.tolist(),"max_top10_contribution_fraction":max(x["fraction_scattered_gain_from_top10_packets"]for x in records),30      "status":"CONTROL_NOT_AUTOMATIC_REPLACEMENT_FOR_PRIMARY_GRID","started_utc":start,"ended_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),"wall_s":time.perf_counter()-t0,"engine_sha256":sha(H/"cylinder_mc.py")}31    dump(O/"summary.json",out);print(json.dumps(out),flush=True)32if __name__=="__main__":main()33