"""Independent light audit of scalar transport outputs; NumPy only. No raw ice table and no full radiative-transfer run are needed. All derived checks are recomputed from native kernels/metadata where available. """ from pathlib import Path import ast,argparse,csv,hashlib,json,math,datetime import numpy as np C=299792458. def sha(path):return hashlib.sha256(Path(path).read_bytes()).hexdigest() def find_root(start=None): starts=[Path(start or Path.cwd()).resolve(),Path(__file__).resolve().parent] for s in starts: for p in [s,*s.parents]: if(p/"modules").is_dir()and(p/"notebooks").is_dir()and(p/"data").is_dir():return p raise FileNotFoundError("Open the verification package or a subdirectory; expected modules/,notebooks/,data/.") def _first(paths): for p in paths: if p.is_file():return p raise FileNotFoundError("Missing packaged input; checked: "+", ".join(str(p)for p in paths)) def load_inputs(root=None,data_dir=None): root=Path(root)if root else find_root();d=Path(data_dir)if data_dir else root/"data" def file(*names):return _first([d/x for x in names]) paths={"results":file("cloud_results.json"),"qa":file("cloud_QA.json","QA.json"), "validation":file("cloud_validation_summary.json","validation_summary.json")} return root,d,{k:json.loads(p.read_text(encoding="utf-8"))for k,p in paths.items()},paths def independent_direct(case,geometry): """Equations D1–D3; on-axis point source and flat cylinder.""" z=case["source_altitude_km"]*1000;z0=geometry["base_m"];zt=geometry["top_m"] R=geometry["observer_z_m"]-z beta=case["total_vertical_cloud_tau"]/(zt-z0) length=max(0.,min(zt-z,zt-z0));tau=beta*length vacuum=1/(4*math.pi*R*R) return dict(beta_ext_m_inverse=beta,R0_m=R,overlying_tau=tau, direct_gain=math.exp(-tau),vacuum_transfer_m_inverse2=vacuum, direct_transfer_m_inverse2=vacuum*math.exp(-tau),flight_time_us=R/C*1e6) def native_paths(d,cid): kernel=_first([d/"native_kernels"/(cid+".npz"),d.parent/"native_kernels"/(cid+".npz"),d/"production/batch_kernels"/(cid+".npz")]) info=_first([d/"batch_info"/(cid+".json"),d.parent/"batch_info"/(cid+".json"),d/"production/batch_info"/(cid+".json")]) return kernel,info def batch_sem(x): x=np.asarray(x,float) return float(x.std(ddof=1)/math.sqrt(len(x)))if len(x)>1 else 0. def audit_native(root=None,data_dir=None): root,d,obj,paths=load_inputs(root,data_dir);result=obj["results"];g=result["model"]["geometry"];rows=[] for case in result["cases"]: cid=case["case_id"];kp,ip=native_paths(d,cid);meta=json.loads(ip.read_text(encoding="utf-8"));bs=meta["batches"] with np.load(kp,allow_pickle=False)as z: mass=np.asarray(z["batch_mass_relative"],float);dt=float(z["dt_s"]) direct=independent_direct(case,g);scattered=mass.sum(axis=1);gains=direct["direct_gain"]+scattered # BxJ center matrix: compute sum of covariance without allocating JxJ. B=len(mass);centered=mass-mass.mean(axis=0) cov_sum=float(np.sum(centered.sum(axis=1)**2)/(B*(B-1)))if B>1 else 0. diagonal_sum=float(np.sum(centered**2)/(B*(B-1)))if B>1 else 0. launch=np.array([b["weighted_launch_mean"]for b in bs]) escape=np.array([b["weighted_escape_fraction"]for b in bs]) trunc=np.array([b["weighted_truncation_fraction"]for b in bs]) boundary=np.array([sum(b["boundary_escape_weights_per_packet"].values())for b in bs]) rows.append({"case_id":cid,"phase":case["phase_id"],"source_km":case["source_altitude_km"], "tau_total":case["total_vertical_cloud_tau"],"tau_overlying_recomputed":direct["overlying_tau"], "direct_gain_recomputed":direct["direct_gain"],"direct_gain_error":direct["direct_gain"]-case["direct_relative_to_vacuum"], "direct_transfer_relative_error":direct["direct_transfer_m_inverse2"]/case["direct_transfer_m_inverse2"]-1, "native_batches":B,"statistical_batches":0 if case["analytic_vacuum_case"]else B, "native_dt_us":dt*1e6,"gain_recomputed":float(gains.mean()),"gain_error":float(gains.mean()-case["total_relative_to_vacuum"]), "SEM_recomputed":batch_sem(gains),"SEM_error":batch_sem(gains)-case["total_mc_standard_error"], "variance_of_sum_from_covariance":cov_sum,"variance_from_independent_bins_incorrect":diagonal_sum, "covariance_identity_error":cov_sum-batch_sem(gains)**2, "offdiagonal_covariance_sum":cov_sum-diagonal_sum, "naive_independent_bin_SEM":math.sqrt(diagonal_sum),"overflow_gain":float(mass[:,-1].mean()), "minimum_native_bin_mass":float(mass.min()),"energy_residual_max_abs":float(np.max(abs(launch-escape-trunc))), "boundary_accounting_max_abs":float(np.max(abs(boundary-escape))), "packets":0 if case["analytic_vacuum_case"]else sum(b["n_packets"]for b in bs), "truncated_packets":sum(b["truncated_packet_count"]for b in bs), "packet_count_balance_pass":all(b["escaped_packet_count"]+b["truncated_packet_count"]==b["n_packets"]for b in bs), "noncausal_delay_violations":sum(b["negative_delay_beyond_tolerance_count"]for b in bs), "quality":"EXPECTED_WARNING_GHM_CONVERGENCE_INCOMPLETE"if case["phase_id"]=="ghm"and case["total_vertical_cloud_tau"]else"conditional_numerical_check", "kernel_sha256":sha(kp),"metadata_sha256":sha(ip)}) summary={"case_count":len(rows),"nonvacuum_packets":sum(r["packets"]for r in rows), "max_abs_direct_error":max(abs(r["direct_gain_error"])for r in rows), "max_abs_direct_transfer_relative_error":max(abs(r["direct_transfer_relative_error"])for r in rows), "max_abs_gain_error":max(abs(r["gain_error"])for r in rows), "max_abs_SEM_error":max(abs(r["SEM_error"])for r in rows), "max_abs_covariance_identity_error":max(abs(r["covariance_identity_error"])for r in rows), "max_abs_energy_residual":max(r["energy_residual_max_abs"]for r in rows), "max_abs_boundary_accounting_error":max(r["boundary_accounting_max_abs"]for r in rows), "truncated_packets":sum(r["truncated_packets"]for r in rows), "noncausal_delay_violations":sum(r["noncausal_delay_violations"]for r in rows), "maximum_observed_overflow_gain":max(r["overflow_gain"]for r in rows), "GHM_warning_case_count":sum(r["quality"].startswith("EXPECTED_WARNING")for r in rows)} assert summary["case_count"]==56 and summary["nonvacuum_packets"]==14400000 assert summary["max_abs_direct_error"]<1e-13 and summary["max_abs_gain_error"]<1e-12 assert summary["max_abs_SEM_error"]<1e-12 assert all(r["packet_count_balance_pass"]and r["minimum_native_bin_mass"]>=0 for r in rows) return rows,summary,obj,paths def hg_check(g,n=100000): """Independent Gauss-Legendre integration and inverse-CDF roundtrip.""" mu,w=np.polynomial.legendre.leggauss(256) f=(1-g*g)/(2*(1+g*g-2*g*mu)**1.5) norm=float(w@f);m1=float(w@(f*mu));m2=float(w@(f*(3*mu*mu-1)/2)) u=np.linspace(0,1,10001) inverse=(1+g*g-((1-g*g)/(1-g+2*g*u))**2)/(2*g)if abs(g)>1e-12 else 2*u-1 back=(1-g*g)/(2*g)*(1/np.sqrt(1+g*g-2*g*inverse)-1/(1+g))if abs(g)>1e-12 else(inverse+1)/2 rng=np.random.default_rng(411337);v=rng.random(n) sample=(1+g*g-((1-g*g)/(1-g+2*g*v))**2)/(2*g)if abs(g)>1e-12 else 2*v-1 se=float(sample.std(ddof=1)/math.sqrt(n)) report={"g":g,"quadrature_nodes":256,"integral_f_dmu":norm,"first_moment":m1,"second_Legendre":m2, "expected":[1,g,g*g],"CDF_roundtrip_max_error":float(abs(back-u).max()), "sample_n":n,"sample_mean_mu":float(sample.mean()),"sample_mean_SE":se,"sample_g_z":float((sample.mean()-g)/se), "scope":"Independent HG consistency check. It does not reconstruct or validate absent GHM phase tables."} assert max(abs(norm-1),abs(m1-g),abs(m2-g*g))<1e-10 and report["CDF_roundtrip_max_error"]<1e-10 return report def _exit_distance(x,z,ux,uz,base,top,radius): # All positions are inside cylinder, y=0; observer rays stay in x-z plane. a=ux*ux;b=x*ux;cc=x*x-radius*radius radial=np.divide(-b+np.sqrt(np.maximum(0,b*b-a*cc)),a,out=np.full_like(x,np.inf),where=a>1e-28) vertical=np.divide(np.where(uz>0,top-z,base-z),uz,out=np.full_like(z,np.inf),where=abs(uz)>1e-15) return np.minimum(radial,vertical) def isotropic_first_scatter_gain(order=64,tau=.3,source=15000.,base=14000.,top=17000.,radius=30000.,observer=550000.): """Independent axisymmetric deterministic first-collision integral, Eq S1.""" if tau==0:return 0. beta=tau/(top-base);r0=observer-source a=(source-base)/math.hypot(radius,source-base);b=(top-source)/math.hypot(radius,top-source) bounds=[-1.,-a,0.,b,1.];q,w=np.polynomial.legendre.leggauss(order);total=0. for low,high in zip(bounds[:-1],bounds[1:]): mu=(low+high)/2+(high-low)*q/2;wm=w*(high-low)/2 ux=np.sqrt(1-mu*mu);exit_s=_exit_distance(np.zeros_like(mu),np.full_like(mu,source),ux,mu,base,top,radius) s=exit_s[:,None]*(q+1)/2;ws=exit_s[:,None]*w/2 x=s*ux[:,None];z=source+s*mu[:,None];rr=np.sqrt(x*x+(observer-z)**2) ox=-x/rr;oz=(observer-z)/rr exit_obs=_exit_distance(x,z,ox,oz,base,top,radius) integrand=.5*beta*np.exp(-beta*(s+exit_obs))*r0*r0*oz/(rr*rr) total+=float(np.sum(wm[:,None]*ws*integrand)) return total def convergence_comparison(validation): out=[] for r in validation["packet_count_controls"]: delta=r["larger_gain"]-r["primary_gain"];den=math.hypot(r["primary_SE"],r["larger_SE"]) out.append({**r,"difference_larger_minus_primary":delta,"difference_in_naive_combined_SE":delta/den, "interpretation":"EXPECTED_WARNING: rare GHM paths remain unresolved; combined-SE discrepancy is diagnostic,not a significance test."if r["phase"]=="ghm"else"Representative ideal-HG estimates agree; not a real-cloud validation."}) return out def code_map(root=None): root=Path(root)if root else find_root();out=[] targets={"cylinder_mc.py":["load_phases","intervals","segment_inside","launch","simulate"], "run_grid.py":["main"],"postprocess.py":["main"],"engine_snapshot.py":["load_phases","intervals","segment_inside","launch","simulate"],"runner_snapshot.py":["main"], "timing_metrics.py":["source_parameters","propagate_waveform","waveform_metrics"]} for name,names in targets.items(): p=root/"reference/cloud"/name if not p.exists():continue tree=ast.parse(p.read_text(encoding="utf-8")) for node in tree.body: if isinstance(node,(ast.FunctionDef,ast.ClassDef))and node.name in names: out.append({"file":p.relative_to(root).as_posix(),"function":node.name,"line":node.lineno,"end_line":node.end_lineno,"sha256":sha(p)}) return out def continuous_source_tau(target_rise_us,decay=.003): """Independent continuous-time bisections,not the reference sampled calibration.""" def peak(tau):return tau*math.log1p(decay/tau) def shape(t,tau):return math.exp(-t/decay)*(-math.expm1(-t/tau)) def rise(tau): tp=peak(tau);pk=shape(tp,tau);times=[] for level in [.1,.9]: lo=0.;hi=tp for _ in range(55): mid=(lo+hi)/2 if shape(mid,tau)=level*height)[0]) cross.append(float(t[j-1]+(level*height-y[j-1])/(y[j]-y[j-1])*(t[j]-t[j-1]))if j else 0.) metrics={"rise_10_90_us":(cross[1]-cross[0])*1e6,"onset_to_peak_us":float(t[k]*1e6), "peak_gain":height,"source_tau_r_us":tau*1e6,"direct_atom_delay_us":0.,"scatter_mass_center_offset_us":dt*1e6} return t,y,ts,source,metrics def worked_hg_packet(g,seed=2026102202,max_collisions=10): """Fixed-seed illustrative history; not an original production path or fluence estimate.""" rng=np.random.default_rng(seed);beta=10/3000.;base=14000.;top=17000.;radius=30000. source=np.array([0.,0.,16000.]);observer=np.array([0.,0.,550000.]);r0=float(np.linalg.norm(observer-source)) def inv(u):return(1+g*g-((1-g*g)/(1-g+2*g*u))**2)/(2*g) def phase(mu):return(1-g*g)/(4*math.pi*(1+g*g-2*g*mu)**1.5) def exit_distance(pos,u): a=u[0]**2+u[1]**2;b=pos[0]*u[0]+pos[1]*u[1];cc=pos[0]**2+pos[1]**2-radius*radius radial=(-b+math.sqrt(max(0,b*b-a*cc)))/a if a>1e-28 else math.inf vertical=((top-pos[2])if u[2]>0 else(base-pos[2]))/u[2]if abs(u[2])>1e-15 else math.inf return min(radial,vertical) chosen=bool(rng.random()<.5);mu=inv(rng.random())if chosen else rng.uniform(-1,1);phi=rng.uniform(0,2*math.pi) u=np.array([math.sqrt(1-mu*mu)*math.cos(phi),math.sqrt(1-mu*mu)*math.sin(phi),mu]) weight=1/(.5+.5*4*math.pi*phase(mu));initial=u.copy();pos=source.copy();path=0.;rows=[] status="illustration_stopped_at_declared_collision_limit" for k in range(max_collisions): optical_free=float(rng.exponential());free=optical_free/beta;edge=exit_distance(pos,u) if free>=edge: status="escaped_before_next_collision";break pos=pos+free*u;path+=free;toward=observer-pos;r=float(np.linalg.norm(toward));v=toward/r cos=float(u@v);psr=phase(cos);lout=min(r,exit_distance(pos,v));tauout=beta*lout survival=math.exp(-tauout);cA=float(v[2]);delay=(path+r-r0)/C per_area=weight*psr*cA/(r*r)*survival sample_cos=inv(rng.random());sample_phi=rng.uniform(0,2*math.pi) axis=np.array([0.,0.,1.])if abs(u[2])<.9 else np.array([1.,0.,0.]) e1=np.cross(u,axis);e1/=np.linalg.norm(e1);e2=np.cross(u,e1) new_u=sample_cos*u+math.sqrt(max(0,1-sample_cos**2))*(math.cos(sample_phi)*e1+math.sin(sample_phi)*e2) rows.append({"collision":k+1,"free_optical_depth":optical_free,"free_path_m":free,"prior_exit_distance_m":edge, "x_m":float(pos[0]),"y_m":float(pos[1]),"z_m":float(pos[2]),"incoming_u":u.tolist(), "source_weight":weight,"accumulated_source_path_m":path,"observer_range_m":r, "cos_to_observer":cos,"p_per_sr":psr,"aperture_cosine":cA,"observer_cloud_length_m":lout, "observer_tau":tauout,"observer_survival":survival,"delay_us":delay*1e6, "NEE_summand_per_m2":per_area,"NEE_summand_relative_vacuum":per_area*4*math.pi*r0*r0, "sampled_next_cos":sample_cos,"outgoing_u":new_u.tolist()}) u=new_u return {"seed":seed,"g":g,"beta_ext_m_inverse":beta,"source_m":source.tolist(),"observer_m":observer.tolist(), "initial_direction":initial.tolist(),"initial_proposal":"0.5isotropic+0.5HG about zenith","sampled_phase_branch":chosen, "initial_source_weight":weight,"max_recorded_collisions":max_collisions,"end_status":status,"collisions":rows, "scope":"Worked mechanical ledger only. Each NEE value is a raw summand before division by the number of emitted histories. A partial single history is not an unbiased complete fluence estimate."} def run_review(root=None,data_dir=None,write=True): root=Path(root)if root else find_root();rows,summary,obj,paths=audit_native(root,data_dir) hg=hg_check(obj["results"]["model"]["phase_input"]["derived_g"]) single=[{"Gauss_order":n,"first_scatter_gain":isotropic_first_scatter_gain(n)}for n in [24,48,96]] mc=obj["validation"]["single_scatter_reference"]["MC"] delta=mc["mean_first_scatter_gain"]-single[-1]["first_scatter_gain"] analytic={"quadrature":single,"retained_MC_mean":mc["mean_first_scatter_gain"],"retained_MC_SE":mc["first_scatter_batch_standard_error"], "difference_in_retained_MC_SE":delta/mc["first_scatter_batch_standard_error"], "qualification":"Own deterministic first-collision equation check,not independent full RT and not a GHM validation."} report={"status":"LIGHT_RECOMPUTATION_PASS_WITH_EXPECTED_GHM_CONVERGENCE_WARNING","generated_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(), "summary":summary,"HG_phase":hg,"single_scatter":analytic,"convergence":convergence_comparison(obj["validation"]), "inputs":[{"name":k,"path":str(p),"sha256":sha(p)}for k,p in paths.items()], "code_map":code_map(root),"worked_packet":worked_hg_packet(hg["g"]),"not_recomputed":["GHM phase normalization without raw LUT","Physical cloud/atmosphere/absorption validity","Noisy onset or population classifier","Per-packet geometry,not retained in kernels"]} if write: out=root/"checks";out.mkdir(exist_ok=True) (out/"transport_review.json").write_text(json.dumps(report,indent=2,allow_nan=False),encoding="utf-8") with(out/"transport_native_reconciliation.csv").open("w",newline="",encoding="utf-8")as f: writer=csv.DictWriter(f,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows) (out/"transport_code_map.json").write_text(json.dumps(report["code_map"],indent=2),encoding="utf-8") (out/"transport_packet_ledger.json").write_text(json.dumps(report["worked_packet"],indent=2),encoding="utf-8") return report,rows,obj if __name__=="__main__": p=argparse.ArgumentParser(description=__doc__);p.add_argument("--root");p.add_argument("--data-dir") a=p.parse_args();report,_,_=run_review(a.root,a.data_dir);print(json.dumps(report["summary"],indent=2));print(report["status"])