Review index · Download raw source
modules/transport_review.py
1"""Independent light audit of scalar transport outputs; NumPy only.2No raw ice table and no full radiative-transfer run are needed.3All derived checks are recomputed from native kernels/metadata where available.4"""5from pathlib import Path6import ast,argparse,csv,hashlib,json,math,datetime7import numpy as np8C=299792458.9def sha(path):return hashlib.sha256(Path(path).read_bytes()).hexdigest()10def find_root(start=None):11 starts=[Path(start or Path.cwd()).resolve(),Path(__file__).resolve().parent]12 for s in starts:13 for p in [s,*s.parents]:14 if(p/"modules").is_dir()and(p/"notebooks").is_dir()and(p/"data").is_dir():return p15 raise FileNotFoundError("Open the verification package or a subdirectory; expected modules/,notebooks/,data/.")16def _first(paths):17 for p in paths:18 if p.is_file():return p19 raise FileNotFoundError("Missing packaged input; checked: "+", ".join(str(p)for p in paths))20def load_inputs(root=None,data_dir=None):21 root=Path(root)if root else find_root();d=Path(data_dir)if data_dir else root/"data"22 def file(*names):return _first([d/x for x in names])23 paths={"results":file("cloud_results.json"),"qa":file("cloud_QA.json","QA.json"),24 "validation":file("cloud_validation_summary.json","validation_summary.json")}25 return root,d,{k:json.loads(p.read_text(encoding="utf-8"))for k,p in paths.items()},paths26def independent_direct(case,geometry):27 """Equations D1–D3; on-axis point source and flat cylinder."""28 z=case["source_altitude_km"]*1000;z0=geometry["base_m"];zt=geometry["top_m"]29 R=geometry["observer_z_m"]-z30 beta=case["total_vertical_cloud_tau"]/(zt-z0)31 length=max(0.,min(zt-z,zt-z0));tau=beta*length32 vacuum=1/(4*math.pi*R*R)33 return dict(beta_ext_m_inverse=beta,R0_m=R,overlying_tau=tau,34 direct_gain=math.exp(-tau),vacuum_transfer_m_inverse2=vacuum,35 direct_transfer_m_inverse2=vacuum*math.exp(-tau),flight_time_us=R/C*1e6)36def native_paths(d,cid):37 kernel=_first([d/"native_kernels"/(cid+".npz"),d.parent/"native_kernels"/(cid+".npz"),d/"production/batch_kernels"/(cid+".npz")])38 info=_first([d/"batch_info"/(cid+".json"),d.parent/"batch_info"/(cid+".json"),d/"production/batch_info"/(cid+".json")])39 return kernel,info40def batch_sem(x):41 x=np.asarray(x,float)42 return float(x.std(ddof=1)/math.sqrt(len(x)))if len(x)>1 else 0.43def audit_native(root=None,data_dir=None):44 root,d,obj,paths=load_inputs(root,data_dir);result=obj["results"];g=result["model"]["geometry"];rows=[]45 for case in result["cases"]:46 cid=case["case_id"];kp,ip=native_paths(d,cid);meta=json.loads(ip.read_text(encoding="utf-8"));bs=meta["batches"]47 with np.load(kp,allow_pickle=False)as z:48 mass=np.asarray(z["batch_mass_relative"],float);dt=float(z["dt_s"])49 direct=independent_direct(case,g);scattered=mass.sum(axis=1);gains=direct["direct_gain"]+scattered50 # BxJ center matrix: compute sum of covariance without allocating JxJ.51 B=len(mass);centered=mass-mass.mean(axis=0)52 cov_sum=float(np.sum(centered.sum(axis=1)**2)/(B*(B-1)))if B>1 else 0.53 diagonal_sum=float(np.sum(centered**2)/(B*(B-1)))if B>1 else 0.54 launch=np.array([b["weighted_launch_mean"]for b in bs])55 escape=np.array([b["weighted_escape_fraction"]for b in bs])56 trunc=np.array([b["weighted_truncation_fraction"]for b in bs])57 boundary=np.array([sum(b["boundary_escape_weights_per_packet"].values())for b in bs])58 rows.append({"case_id":cid,"phase":case["phase_id"],"source_km":case["source_altitude_km"],59 "tau_total":case["total_vertical_cloud_tau"],"tau_overlying_recomputed":direct["overlying_tau"],60 "direct_gain_recomputed":direct["direct_gain"],"direct_gain_error":direct["direct_gain"]-case["direct_relative_to_vacuum"],61 "direct_transfer_relative_error":direct["direct_transfer_m_inverse2"]/case["direct_transfer_m_inverse2"]-1,62 "native_batches":B,"statistical_batches":0 if case["analytic_vacuum_case"]else B,63 "native_dt_us":dt*1e6,"gain_recomputed":float(gains.mean()),"gain_error":float(gains.mean()-case["total_relative_to_vacuum"]),64 "SEM_recomputed":batch_sem(gains),"SEM_error":batch_sem(gains)-case["total_mc_standard_error"],65 "variance_of_sum_from_covariance":cov_sum,"variance_from_independent_bins_incorrect":diagonal_sum,66 "covariance_identity_error":cov_sum-batch_sem(gains)**2,67 "offdiagonal_covariance_sum":cov_sum-diagonal_sum,68 "naive_independent_bin_SEM":math.sqrt(diagonal_sum),"overflow_gain":float(mass[:,-1].mean()),69 "minimum_native_bin_mass":float(mass.min()),"energy_residual_max_abs":float(np.max(abs(launch-escape-trunc))),70 "boundary_accounting_max_abs":float(np.max(abs(boundary-escape))),71 "packets":0 if case["analytic_vacuum_case"]else sum(b["n_packets"]for b in bs),72 "truncated_packets":sum(b["truncated_packet_count"]for b in bs),73 "packet_count_balance_pass":all(b["escaped_packet_count"]+b["truncated_packet_count"]==b["n_packets"]for b in bs),74 "noncausal_delay_violations":sum(b["negative_delay_beyond_tolerance_count"]for b in bs),75 "quality":"EXPECTED_WARNING_GHM_CONVERGENCE_INCOMPLETE"if case["phase_id"]=="ghm"and case["total_vertical_cloud_tau"]else"conditional_numerical_check",76 "kernel_sha256":sha(kp),"metadata_sha256":sha(ip)})77 summary={"case_count":len(rows),"nonvacuum_packets":sum(r["packets"]for r in rows),78 "max_abs_direct_error":max(abs(r["direct_gain_error"])for r in rows),79 "max_abs_direct_transfer_relative_error":max(abs(r["direct_transfer_relative_error"])for r in rows),80 "max_abs_gain_error":max(abs(r["gain_error"])for r in rows),81 "max_abs_SEM_error":max(abs(r["SEM_error"])for r in rows),82 "max_abs_covariance_identity_error":max(abs(r["covariance_identity_error"])for r in rows),83 "max_abs_energy_residual":max(r["energy_residual_max_abs"]for r in rows),84 "max_abs_boundary_accounting_error":max(r["boundary_accounting_max_abs"]for r in rows),85 "truncated_packets":sum(r["truncated_packets"]for r in rows),86 "noncausal_delay_violations":sum(r["noncausal_delay_violations"]for r in rows),87 "maximum_observed_overflow_gain":max(r["overflow_gain"]for r in rows),88 "GHM_warning_case_count":sum(r["quality"].startswith("EXPECTED_WARNING")for r in rows)}89 assert summary["case_count"]==56 and summary["nonvacuum_packets"]==1440000090 assert summary["max_abs_direct_error"]<1e-13 and summary["max_abs_gain_error"]<1e-1291 assert summary["max_abs_SEM_error"]<1e-1292 assert all(r["packet_count_balance_pass"]and r["minimum_native_bin_mass"]>=0 for r in rows)93 return rows,summary,obj,paths94def hg_check(g,n=100000):95 """Independent Gauss-Legendre integration and inverse-CDF roundtrip."""96 mu,w=np.polynomial.legendre.leggauss(256)97 f=(1-g*g)/(2*(1+g*g-2*g*mu)**1.5)98 norm=float(w@f);m1=float(w@(f*mu));m2=float(w@(f*(3*mu*mu-1)/2))99 u=np.linspace(0,1,10001)100 inverse=(1+g*g-((1-g*g)/(1-g+2*g*u))**2)/(2*g)if abs(g)>1e-12 else 2*u-1101 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)/2102 rng=np.random.default_rng(411337);v=rng.random(n)103 sample=(1+g*g-((1-g*g)/(1-g+2*g*v))**2)/(2*g)if abs(g)>1e-12 else 2*v-1104 se=float(sample.std(ddof=1)/math.sqrt(n))105 report={"g":g,"quadrature_nodes":256,"integral_f_dmu":norm,"first_moment":m1,"second_Legendre":m2,106 "expected":[1,g,g*g],"CDF_roundtrip_max_error":float(abs(back-u).max()),107 "sample_n":n,"sample_mean_mu":float(sample.mean()),"sample_mean_SE":se,"sample_g_z":float((sample.mean()-g)/se),108 "scope":"Independent HG consistency check. It does not reconstruct or validate absent GHM phase tables."}109 assert max(abs(norm-1),abs(m1-g),abs(m2-g*g))<1e-10 and report["CDF_roundtrip_max_error"]<1e-10110 return report111def _exit_distance(x,z,ux,uz,base,top,radius):112 # All positions are inside cylinder, y=0; observer rays stay in x-z plane.113 a=ux*ux;b=x*ux;cc=x*x-radius*radius114 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)115 vertical=np.divide(np.where(uz>0,top-z,base-z),uz,out=np.full_like(z,np.inf),where=abs(uz)>1e-15)116 return np.minimum(radial,vertical)117def isotropic_first_scatter_gain(order=64,tau=.3,source=15000.,base=14000.,top=17000.,radius=30000.,observer=550000.):118 """Independent axisymmetric deterministic first-collision integral, Eq S1."""119 if tau==0:return 0.120 beta=tau/(top-base);r0=observer-source121 a=(source-base)/math.hypot(radius,source-base);b=(top-source)/math.hypot(radius,top-source)122 bounds=[-1.,-a,0.,b,1.];q,w=np.polynomial.legendre.leggauss(order);total=0.123 for low,high in zip(bounds[:-1],bounds[1:]):124 mu=(low+high)/2+(high-low)*q/2;wm=w*(high-low)/2125 ux=np.sqrt(1-mu*mu);exit_s=_exit_distance(np.zeros_like(mu),np.full_like(mu,source),ux,mu,base,top,radius)126 s=exit_s[:,None]*(q+1)/2;ws=exit_s[:,None]*w/2127 x=s*ux[:,None];z=source+s*mu[:,None];rr=np.sqrt(x*x+(observer-z)**2)128 ox=-x/rr;oz=(observer-z)/rr129 exit_obs=_exit_distance(x,z,ox,oz,base,top,radius)130 integrand=.5*beta*np.exp(-beta*(s+exit_obs))*r0*r0*oz/(rr*rr)131 total+=float(np.sum(wm[:,None]*ws*integrand))132 return total133def convergence_comparison(validation):134 out=[]135 for r in validation["packet_count_controls"]:136 delta=r["larger_gain"]-r["primary_gain"];den=math.hypot(r["primary_SE"],r["larger_SE"])137 out.append({**r,"difference_larger_minus_primary":delta,"difference_in_naive_combined_SE":delta/den,138 "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."})139 return out140def code_map(root=None):141 root=Path(root)if root else find_root();out=[]142 targets={"cylinder_mc.py":["load_phases","intervals","segment_inside","launch","simulate"],143 "run_grid.py":["main"],"postprocess.py":["main"],"engine_snapshot.py":["load_phases","intervals","segment_inside","launch","simulate"],"runner_snapshot.py":["main"],144 "timing_metrics.py":["source_parameters","propagate_waveform","waveform_metrics"]}145 for name,names in targets.items():146 p=root/"reference/cloud"/name147 if not p.exists():continue148 tree=ast.parse(p.read_text(encoding="utf-8"))149 for node in tree.body:150 if isinstance(node,(ast.FunctionDef,ast.ClassDef))and node.name in names:151 out.append({"file":p.relative_to(root).as_posix(),"function":node.name,"line":node.lineno,"end_line":node.end_lineno,"sha256":sha(p)})152 return out153def continuous_source_tau(target_rise_us,decay=.003):154 """Independent continuous-time bisections,not the reference sampled calibration."""155 def peak(tau):return tau*math.log1p(decay/tau)156 def shape(t,tau):return math.exp(-t/decay)*(-math.expm1(-t/tau))157 def rise(tau):158 tp=peak(tau);pk=shape(tp,tau);times=[]159 for level in [.1,.9]:160 lo=0.;hi=tp161 for _ in range(55):162 mid=(lo+hi)/2163 if shape(mid,tau)<level*pk:lo=mid164 else:hi=mid165 times.append((lo+hi)/2)166 return times[1]-times[0]167 lo=1e-10;hi=.003168 for _ in range(60):169 mid=(lo+hi)/2170 if rise(mid)<target_rise_us*1e-6:lo=mid171 else:hi=mid172 return(lo+hi)/2173def independent_pulse(native_mass,direct_gain,target_rise_us=80.,native_dt=2e-6):174 """Independent source calibration and discrete convolution of response MASS."""175 tau=continuous_source_tau(target_rise_us);dt=native_dt/2176 tp=tau*math.log1p(.003/tau);pk=math.exp(-tp/.003)*(-math.expm1(-tp/tau))177 ts=np.arange(0,.036+dt/2,dt);source=np.exp(-ts/.003)*(-np.expm1(-ts/tau))/pk178 response=np.zeros(2*(len(native_mass)-1)+1);response[0]=direct_gain;response[1::2]=native_mass[:-1]179 n=len(source)+len(response)-1;nfft=1<<(n-1).bit_length()180 y=np.fft.irfft(np.fft.rfft(source,nfft)*np.fft.rfft(response,nfft),nfft)[:n]181 y=np.maximum(y,0);t=np.arange(n)*dt182 k=int(y.argmax());height=float(y[k]);cross=[]183 for level in [.1,.9]:184 j=int(np.flatnonzero(y[:k+1]>=level*height)[0])185 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.)186 metrics={"rise_10_90_us":(cross[1]-cross[0])*1e6,"onset_to_peak_us":float(t[k]*1e6),187 "peak_gain":height,"source_tau_r_us":tau*1e6,"direct_atom_delay_us":0.,"scatter_mass_center_offset_us":dt*1e6}188 return t,y,ts,source,metrics189190def worked_hg_packet(g,seed=2026102202,max_collisions=10):191 """Fixed-seed illustrative history; not an original production path or fluence estimate."""192 rng=np.random.default_rng(seed);beta=10/3000.;base=14000.;top=17000.;radius=30000.193 source=np.array([0.,0.,16000.]);observer=np.array([0.,0.,550000.]);r0=float(np.linalg.norm(observer-source))194 def inv(u):return(1+g*g-((1-g*g)/(1-g+2*g*u))**2)/(2*g)195 def phase(mu):return(1-g*g)/(4*math.pi*(1+g*g-2*g*mu)**1.5)196 def exit_distance(pos,u):197 a=u[0]**2+u[1]**2;b=pos[0]*u[0]+pos[1]*u[1];cc=pos[0]**2+pos[1]**2-radius*radius198 radial=(-b+math.sqrt(max(0,b*b-a*cc)))/a if a>1e-28 else math.inf199 vertical=((top-pos[2])if u[2]>0 else(base-pos[2]))/u[2]if abs(u[2])>1e-15 else math.inf200 return min(radial,vertical)201 chosen=bool(rng.random()<.5);mu=inv(rng.random())if chosen else rng.uniform(-1,1);phi=rng.uniform(0,2*math.pi)202 u=np.array([math.sqrt(1-mu*mu)*math.cos(phi),math.sqrt(1-mu*mu)*math.sin(phi),mu])203 weight=1/(.5+.5*4*math.pi*phase(mu));initial=u.copy();pos=source.copy();path=0.;rows=[]204 status="illustration_stopped_at_declared_collision_limit"205 for k in range(max_collisions):206 optical_free=float(rng.exponential());free=optical_free/beta;edge=exit_distance(pos,u)207 if free>=edge:208 status="escaped_before_next_collision";break209 pos=pos+free*u;path+=free;toward=observer-pos;r=float(np.linalg.norm(toward));v=toward/r210 cos=float(u@v);psr=phase(cos);lout=min(r,exit_distance(pos,v));tauout=beta*lout211 survival=math.exp(-tauout);cA=float(v[2]);delay=(path+r-r0)/C212 per_area=weight*psr*cA/(r*r)*survival213 sample_cos=inv(rng.random());sample_phi=rng.uniform(0,2*math.pi)214 axis=np.array([0.,0.,1.])if abs(u[2])<.9 else np.array([1.,0.,0.])215 e1=np.cross(u,axis);e1/=np.linalg.norm(e1);e2=np.cross(u,e1)216 new_u=sample_cos*u+math.sqrt(max(0,1-sample_cos**2))*(math.cos(sample_phi)*e1+math.sin(sample_phi)*e2)217 rows.append({"collision":k+1,"free_optical_depth":optical_free,"free_path_m":free,"prior_exit_distance_m":edge,218 "x_m":float(pos[0]),"y_m":float(pos[1]),"z_m":float(pos[2]),"incoming_u":u.tolist(),219 "source_weight":weight,"accumulated_source_path_m":path,"observer_range_m":r,220 "cos_to_observer":cos,"p_per_sr":psr,"aperture_cosine":cA,"observer_cloud_length_m":lout,221 "observer_tau":tauout,"observer_survival":survival,"delay_us":delay*1e6,222 "NEE_summand_per_m2":per_area,"NEE_summand_relative_vacuum":per_area*4*math.pi*r0*r0,223 "sampled_next_cos":sample_cos,"outgoing_u":new_u.tolist()})224 u=new_u225 return {"seed":seed,"g":g,"beta_ext_m_inverse":beta,"source_m":source.tolist(),"observer_m":observer.tolist(),226 "initial_direction":initial.tolist(),"initial_proposal":"0.5isotropic+0.5HG about zenith","sampled_phase_branch":chosen,227 "initial_source_weight":weight,"max_recorded_collisions":max_collisions,"end_status":status,"collisions":rows,228 "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."}229230def run_review(root=None,data_dir=None,write=True):231 root=Path(root)if root else find_root();rows,summary,obj,paths=audit_native(root,data_dir)232 hg=hg_check(obj["results"]["model"]["phase_input"]["derived_g"])233 single=[{"Gauss_order":n,"first_scatter_gain":isotropic_first_scatter_gain(n)}for n in [24,48,96]]234 mc=obj["validation"]["single_scatter_reference"]["MC"]235 delta=mc["mean_first_scatter_gain"]-single[-1]["first_scatter_gain"]236 analytic={"quadrature":single,"retained_MC_mean":mc["mean_first_scatter_gain"],"retained_MC_SE":mc["first_scatter_batch_standard_error"],237 "difference_in_retained_MC_SE":delta/mc["first_scatter_batch_standard_error"],238 "qualification":"Own deterministic first-collision equation check,not independent full RT and not a GHM validation."}239 report={"status":"LIGHT_RECOMPUTATION_PASS_WITH_EXPECTED_GHM_CONVERGENCE_WARNING","generated_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),240 "summary":summary,"HG_phase":hg,"single_scatter":analytic,"convergence":convergence_comparison(obj["validation"]),241 "inputs":[{"name":k,"path":str(p),"sha256":sha(p)}for k,p in paths.items()],242 "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"]}243 if write:244 out=root/"checks";out.mkdir(exist_ok=True)245 (out/"transport_review.json").write_text(json.dumps(report,indent=2,allow_nan=False),encoding="utf-8")246 with(out/"transport_native_reconciliation.csv").open("w",newline="",encoding="utf-8")as f:247 writer=csv.DictWriter(f,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows)248 (out/"transport_code_map.json").write_text(json.dumps(report["code_map"],indent=2),encoding="utf-8")249 (out/"transport_packet_ledger.json").write_text(json.dumps(report["worked_packet"],indent=2),encoding="utf-8")250 return report,rows,obj251if __name__=="__main__":252 p=argparse.ArgumentParser(description=__doc__);p.add_argument("--root");p.add_argument("--data-dir")253 a=p.parse_args();report,_,_=run_review(a.root,a.data_dir);print(json.dumps(report["summary"],indent=2));print(report["status"])254