Review index · Download raw source
reference/cloud/postprocess.py
1"""Postprocess completed independent MC; retain batch statistics and definitions."""2from pathlib import Path3import json,csv,math,hashlib,datetime,sys4import numpy as np5from cylinder_mc import sha6import timing_metrics as tm7HERE=Path(__file__).resolve().parent;RUN=HERE/"production"8def se(x):9 x=np.asarray(x,float)10 return float(x.std(ddof=1)/math.sqrt(len(x)))if len(x)>1 else 0.11def tidy(x):12 if isinstance(x,dict):return {k:tidy(v)for k,v in x.items()}13 if isinstance(x,list):return [tidy(v)for v in x]14 if isinstance(x,(np.floating,np.integer)):x=x.item()15 if isinstance(x,float)and not math.isfinite(x):return None16 return x17def main():18 m=json.loads((RUN/"run_manifest.json").read_text(encoding="utf-8"))19 if m["status"]!="COMPLETE":raise RuntimeError("Grid is not complete; no final results written")20 dt=m["dt_s"];fine=dt/2;nb=int(m["response_window_s"]/dt)21 # Interactive bins preserve2us onset information and coarsen only late tails.22 edges_us=np.unique(np.r_[np.arange(0,202,2),np.arange(205,1005,5),np.arange(1020,12001,20)])23 # Fine sampled pulse times are common to everycase; first12ms only.24 curve_idx=np.unique(np.r_[np.arange(0,501),np.arange(520,12001,20)]).astype(int)25 curve_t_us=curve_idx*fine*1e626 cases=[];summary=[];timing_rows=[]27 for cid in m["case_ids"]:28 entry=json.loads((RUN/"batch_info"/(cid+".json")).read_text(encoding="utf-8"))29 bs=entry["batches"];raw=np.load(RUN/"batch_kernels"/(cid+".npz"))["batch_mass_relative"]30 avg=raw.mean(axis=0);std=raw.std(axis=0,ddof=1)/math.sqrt(len(raw))if len(raw)>1 else np.zeros_like(avg)31 b0=bs[0];direct=b0["direct_relative_to_vacuum"];vac=b0["vacuum_transfer_m_inverse2"]32 scattered=np.array([b["scattered_relative_to_vacuum"]for b in bs]);gain=direct+scattered33 relative_se=se(gain)/float(gain.mean())34 coarse=[]35 for row in raw:36 centers=(np.arange(nb)+.5)*dt*1e637 coarse.append(np.histogram(centers,bins=edges_us,weights=row[:-1])[0])38 coarse=np.asarray(coarse)39 small={"case_id":cid,"phase_id":b0["phase_id"],40 "source_altitude_km":b0["source_altitude_m"]/1000,41 "total_vertical_cloud_tau":b0["total_tau"],"overlying_cloud_tau":b0["overlying_tau"],42 "source_location":"below_top"if b0["source_altitude_m"]<m["geometry"]["top_m"]else("at_top"if b0["source_altitude_m"]==m["geometry"]["top_m"]else"above_top"),43 "direct_relative_to_vacuum":direct,44 "scattered_relative_to_vacuum":float(scattered.mean()),45 "scattered_mc_standard_error":se(scattered),46 "total_relative_to_vacuum":float(gain.mean()),47 "total_mc_standard_error":se(gain),48 "relative_mc_standard_error_total":relative_se,49 "vacuum_transfer_m_inverse2":vac,"unscattered_flight_time_us":(m["geometry"]["observer_z_m"]-b0["source_altitude_m"])/299792458.*1e6,50 "direct_transfer_m_inverse2":direct*vac,51 "scattered_transfer_m_inverse2":float(scattered.mean())*vac,52 "total_transfer_m_inverse2":float(gain.mean())*vac,53 "total_transfer_mc_standard_error_m_inverse2":se(gain)*vac,54 "n_batches":len(bs) if b0["total_tau"] else 0,"packets_per_batch":b0["n_packets"] if b0["total_tau"] else 0,55 "analytic_vacuum_case":b0["total_tau"]==0,56 "packet_count":sum(b["n_packets"]for b in bs)if b0["total_tau"]else 0,57 "truncated_packet_count":sum(b["truncated_packet_count"]for b in bs),58 "weighted_truncation_fraction":float(np.mean([b["weighted_truncation_fraction"]for b in bs])),59 "energy_balance_max_abs_per_packet":max(abs(b["energy_balance_residual_per_packet"])for b in bs),60 "weighted_launch_mean":float(np.mean([b["weighted_launch_mean"]for b in bs])),61 "weighted_launch_mc_standard_error":se([b["weighted_launch_mean"]for b in bs]),62 "weighted_escape_fraction_mean":float(np.mean([b["weighted_escape_fraction"]for b in bs])),63 "weighted_upward_escape_fraction_mean":float(np.mean([b["weighted_upward_escape_fraction"]for b in bs])) if b0["total_tau"] else .5,64 "weighted_upward_escape_fraction_mc_se":se([b["weighted_upward_escape_fraction"]for b in bs]) if b0["total_tau"] else 0.,65 "escape_statistics_method":"weighted_MC" if b0["total_tau"] else "analytic_isotropic_vacuum",66 "mean_scattering_events_per_packet":float(np.mean([b["mean_collisions_per_packet"]for b in bs])),67 "min_scattered_delay_us":min([b["min_scattered_delay_s"]*1e6 for b in bs if b["min_scattered_delay_s"]is not None],default=None),68 "response_overflow_relative_to_vacuum":float(avg[-1]),69 "run_seconds":sum(b["elapsed_s"]for b in bs),70 "quality_flags":(["GHM_rare_path_convergence_not_uniformly_established"] if b0["phase_id"]=="ghm" and b0["total_tau"] else [])+(["high_point_nadir_MC_uncertainty"]if relative_se>.05 else [])+(["max_scatter_truncation"]if any(b["truncated_packet_count"]for b in bs)else[]),71 "batch_info_file":"production/batch_info/"+cid+".json",72 "batch_kernel_file":"production/batch_kernels/"+cid+".npz"}73 if b0["total_tau"]:74 small["finite_receiver_scattered_gains"]={deg:{"mean":float(np.mean([b["finite_receiver_scattered_gain"][deg]for b in bs])),"batch_standard_error":se([b["finite_receiver_scattered_gain"][deg]for b in bs])}for deg in ["3.0","5.0","10.0"]}75 # Explicitcenter rule: directatom att0, MCbins at2n+1us.76 masses=np.zeros((len(raw),2*nb+1));masses[:,0]=direct;masses[:,1:2*nb:2]=raw[:,:-1]77 mass=masses.mean(axis=0)78 pulse_curves={};timing=[]79 for rise in [30.,80.,200.]:80 for fc in [0.,10000.,50000.]:81 t,y=tm.propagate_waveform(mass/fine,fine,rise,None if not fc else fc)82 met=tm.waveform_metrics(t,y)83 batchmets=[]84 for bmass in masses:85 tb,yb=tm.propagate_waveform(bmass/fine,fine,rise,None if not fc else fc)86 batchmets.append(tm.waveform_metrics(tb,yb))87 item={"intrinsic_rise_10_90_us":rise,"analog_one_pole_fc_hz":fc,88 "rise_10_90_us":met["rise_10_90_us"],"onset_to_peak_us":met["onset_to_peak_us"],89 "rise_batch_standard_error_us":se([x["rise_10_90_us"]for x in batchmets]),90 "onset_peak_batch_standard_error_us":se([x["onset_to_peak_us"]for x in batchmets]),91 "batch_rise_10_90_us":[x["rise_10_90_us"]for x in batchmets] if b0["total_tau"] else [],92 "batch_onset_to_peak_us":[x["onset_to_peak_us"]for x in batchmets] if b0["total_tau"] else [],93 "peak_gain_relative_to_vacuum_unit_peak_source":met["peak_amplitude"],94 "peak_at_window_end":met["peak_at_window_end"]}95 timing.append(item);timing_rows.append({"case_id":cid,**{k:v for k,v in item.items() if not isinstance(v,list)}})96 pulse_curves[f"rise{int(rise)}_fc{int(fc)}"]=(y[curve_idx]/y.max()).tolist()97 row=dict(small);row["timing"]=timing98 row["scattered_bin_mass_relative"]=coarse.mean(axis=0).tolist()99 row["scattered_bin_mass_mc_se"]=((coarse.std(axis=0,ddof=1)/math.sqrt(len(coarse)))if len(coarse)>1 else np.zeros(coarse.shape[1])).tolist()100 row["pulse_curves_unit_peak"]=pulse_curves101 row["batch_total_gain"]=[b["total_relative_to_vacuum"]for b in bs] if b0["total_tau"] else []102 row["batch_upward_escape_fraction"]=[b["weighted_upward_escape_fraction"]for b in bs] if b0["total_tau"] else []103 cases.append(row);summary.append(small)104 output={"schema_version":"1.0","status":"EXECUTED_INDEPENDENT_SCALAR_CONDITIONAL_MC_NOT_REAL_CLOUD_VALIDATION",105 "generated_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),106 "model":{"geometry":m["geometry"],"source_altitude_km":m["source_altitude_km"],107 "total_tau_grid":m["total_tau"],"phases":["ghm","hg"],"physical_source_emission":"Isotropic; exact source-only importance weights used internally below the cloud top.",108 "medium":"Finite homogeneous cylinder; vacuum exterior; scalar conservative SSA=1. No Rayleigh scattering, gas absorption, aerosol, ground reflection, polarization or real-cloud heterogeneity.",109 "wavelength_label":"337.1nm numerical input coordinate; linear interpolation of native330/340nm data. N2 band integration and air/vacuum wavelength registration are not modeled.",110 "phase_input":m["phase_input"],111 "run_elapsed_s":m["elapsed_s"],112 "tau_sweep_scope":"Changing total tau changes extinction throughout the whole cylinder, including the cloud below and above an internal source. Overlying tau alone does not specify this family.",113 "earth_geometry":"Locally flat cloud boundaries and straight vacuum rays; no Earth curvature, atmospheric refraction or finite source size.",114 "boundary_conditions":"Vacuum after exit from the convex finite cylinder; no reflecting surface or ground. A photon outside cannot re-enter without exterior scattering, which is omitted.",115 "receiver":"Small flat horizontal aperture at nadir; collecting area enters once and FOV is a60degree full circular cone. Transfer is per square metre before any QE, PDE or throughput.","nonvacuum_packets":m["nonvacuum_total_packets"],116 "independent_batches_per_nonvacuum_case":m["independent_batches"],117 "engine_sha256_at_run":m["engine_sha256"],"runner_sha256_at_run":m["runner_sha256"],118 "engine_file":"cylinder_mc.py","runner_file":"run_grid.py","raw_LUT_included":False},119 "definitions":{"transfer_m_inverse2":"Expected photon fluence at the aperture per emitted photon,m^-2. Multiply by small aperture area in m² for the collected photon fraction. QE/PDE is not included.",120 "relative_to_vacuum":"Transfer divided by1/(4pi R0²),where R0 is the source-observer distance. This directional ratio is not a probability or an event rate.",121 "tau_total":"Extinction coefficient times total3km cloud thickness.",122 "tau_overlying":"Extinction along the straight vertical source-observer line; sources at/above the top have zero overlying cloud optical depth.",123 "MC_error":"Sample SD of independent batch means divided by sqrt(B). Rare forward paths may remain underresolved, so SEM is not a guaranteed confidence interval. Timing SEM is a diagnostic from nonlinear metrics of individual batch waveforms.",124 "time":"Excess physical path delay relative to unscattered source-observer flight. The direct atom at delay0 is supplied separately.",125 "response_bins":"Scattered bin MASS relative to vacuum. Divide by bin width in seconds for response density in s^-1. Bin errors cannot be summed in quadrature because bins share photon paths.",126 "rise":"Rising10–90% of the waveform peak and KNOWN MODEL ONSET-to-peak are reported separately. Noisy measured onset or population classification is not validated.",127 "source_waveform":"Assumed exp(-t/3ms)*(1-exp(-t/tau_r));tau_r is tuned to the intrinsic10–90% rise. These intrinsic rises and the3ms decay are not fitted NBE observations.",128 "electronics":"Ideal first-order low-pass only; not an actual analog or antialias chain. Cutoff frequency is distinct from sampling rate.",129 "finite_receiver_control":"Escaped-ray scattered fluence averaged over a large horizontal receiver disk subtending3/5/10degrees. This diagnostic is not the exact point-nadir value or an independent-code validation."},130 "time_grid":{"native_scattered_bin_width_us":dt*1e6,"response_window_us":m["response_window_s"]*1e6,131 "scattered_plot_bin_edges_us":edges_us.tolist(),"pulse_plot_time_us":curve_t_us.tolist(),132 "pulse_plot_window":"First12ms plotted; source and convolution metrics use an approximately48ms record.",133 "convolution_sample_dt_us":fine*1e6,"scatter_mass_placement":"At native bin centers; direct atom at delay0. Unresolved within-bin spread is limited by the native bin width."},134 "source_parameters":[tm.source_parameters(x)for x in [30.,80.,200.]],135 "numerical_warning":"The tabulated GHM point-nadir rare-path estimator is not uniformly converged. Representative increased-packet and source-proposal controls can differ more than naive finite-batch SEM suggests. Do not infer a physical GHM/HG rise-time bias from unconverged ratios. All original batches and controls are retained.",136 "qa_files":{"phase_and_geometry":"production/phase_and_geometry_QA.json","full_audit":"QA.json","representative_controls":"validation/validation_results.json","cone_proposal_control":"validation/cone_control/summary.json","matched_HG_CloudScat_control":"validation/cloudscat_control/comparison.json","compact_validation_summary":"validation_summary.json","readme":"README.md","vacuum_metadata_repair":"vacuum_metadata_repair.json"},137 "cases":cases,138 "limitations":["Conditional model comparison,not a universal transmission curve or inversion.","Equal g does not fix higher phase-function moments. This grid supplies conditional estimates with MC errors; the current GHM rare-path controls do not establish a physical timing bias.","Two matched homogeneous-HG CloudScat controls were executed and agree within their combined MC errors. This independent-code check does not validate tabulated GHM, atmospheric transport or a real-cloud classifier; own analytic tests are labeled separately.","The native processed phase grid is assumed. Interpolation, ice habit and size-distribution uncertainties are not represented by MC error.","No actual event-brightness population,absolute rate,current calibration,flight noise or rise-time classifier is provided."],139 "primary_urls":["https://gmd.copernicus.org/articles/13/5549/2020/","https://www.ssec.wisc.edu/ice_models/polarization.html","https://docs.hygeos.com/public.php/webdav/ic_baum_ghm_sol.nc"]}140 (HERE/"cloud_results.json").write_text(json.dumps(tidy(output),separators=(",",":"),allow_nan=False),encoding="utf-8")141 # CompactscalarCSV; nested diagnosticfields remain inJSON.142 scalar=[{k:v for k,v in x.items()if not isinstance(v,(dict,list))}for x in summary]143 with(HERE/"cloud_summary.csv").open("w",newline="",encoding="utf-8")as f:144 w=csv.DictWriter(f,fieldnames=list(scalar[0]));w.writeheader();w.writerows(scalar)145 with(HERE/"cloud_timing.csv").open("w",newline="",encoding="utf-8")as f:146 w=csv.DictWriter(f,fieldnames=list(timing_rows[0]));w.writeheader();w.writerows(timing_rows)147 print(json.dumps({"cases":len(cases),"timing_rows":len(timing_rows),"cloud_results_bytes":(HERE/"cloud_results.json").stat().st_size,"high_SE_cases":sum("high_point_nadir_MC_uncertainty" in x["quality_flags"] for x in cases),"GHM_nonuniform_convergence_flags":sum("GHM_rare_path_convergence_not_uniformly_established" in x["quality_flags"] for x in cases)}))148if __name__=="__main__":main()149