"""Postprocess completed independent MC; retain batch statistics and definitions."""
from pathlib import Path
import json,csv,math,hashlib,datetime,sys
import numpy as np
from cylinder_mc import sha
import timing_metrics as tm
HERE=Path(__file__).resolve().parent;RUN=HERE/"production"
def se(x):
    x=np.asarray(x,float)
    return float(x.std(ddof=1)/math.sqrt(len(x)))if len(x)>1 else 0.
def tidy(x):
    if isinstance(x,dict):return {k:tidy(v)for k,v in x.items()}
    if isinstance(x,list):return [tidy(v)for v in x]
    if isinstance(x,(np.floating,np.integer)):x=x.item()
    if isinstance(x,float)and not math.isfinite(x):return None
    return x
def main():
    m=json.loads((RUN/"run_manifest.json").read_text(encoding="utf-8"))
    if m["status"]!="COMPLETE":raise RuntimeError("Grid is not complete; no final results written")
    dt=m["dt_s"];fine=dt/2;nb=int(m["response_window_s"]/dt)
    # Interactive bins preserve2us onset information and coarsen only late tails.
    edges_us=np.unique(np.r_[np.arange(0,202,2),np.arange(205,1005,5),np.arange(1020,12001,20)])
    # Fine sampled pulse times are common to everycase; first12ms only.
    curve_idx=np.unique(np.r_[np.arange(0,501),np.arange(520,12001,20)]).astype(int)
    curve_t_us=curve_idx*fine*1e6
    cases=[];summary=[];timing_rows=[]
    for cid in m["case_ids"]:
        entry=json.loads((RUN/"batch_info"/(cid+".json")).read_text(encoding="utf-8"))
        bs=entry["batches"];raw=np.load(RUN/"batch_kernels"/(cid+".npz"))["batch_mass_relative"]
        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)
        b0=bs[0];direct=b0["direct_relative_to_vacuum"];vac=b0["vacuum_transfer_m_inverse2"]
        scattered=np.array([b["scattered_relative_to_vacuum"]for b in bs]);gain=direct+scattered
        relative_se=se(gain)/float(gain.mean())
        coarse=[]
        for row in raw:
            centers=(np.arange(nb)+.5)*dt*1e6
            coarse.append(np.histogram(centers,bins=edges_us,weights=row[:-1])[0])
        coarse=np.asarray(coarse)
        small={"case_id":cid,"phase_id":b0["phase_id"],
          "source_altitude_km":b0["source_altitude_m"]/1000,
          "total_vertical_cloud_tau":b0["total_tau"],"overlying_cloud_tau":b0["overlying_tau"],
          "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"),
          "direct_relative_to_vacuum":direct,
          "scattered_relative_to_vacuum":float(scattered.mean()),
          "scattered_mc_standard_error":se(scattered),
          "total_relative_to_vacuum":float(gain.mean()),
          "total_mc_standard_error":se(gain),
          "relative_mc_standard_error_total":relative_se,
          "vacuum_transfer_m_inverse2":vac,"unscattered_flight_time_us":(m["geometry"]["observer_z_m"]-b0["source_altitude_m"])/299792458.*1e6,
          "direct_transfer_m_inverse2":direct*vac,
          "scattered_transfer_m_inverse2":float(scattered.mean())*vac,
          "total_transfer_m_inverse2":float(gain.mean())*vac,
          "total_transfer_mc_standard_error_m_inverse2":se(gain)*vac,
          "n_batches":len(bs) if b0["total_tau"] else 0,"packets_per_batch":b0["n_packets"] if b0["total_tau"] else 0,
          "analytic_vacuum_case":b0["total_tau"]==0,
          "packet_count":sum(b["n_packets"]for b in bs)if b0["total_tau"]else 0,
          "truncated_packet_count":sum(b["truncated_packet_count"]for b in bs),
          "weighted_truncation_fraction":float(np.mean([b["weighted_truncation_fraction"]for b in bs])),
          "energy_balance_max_abs_per_packet":max(abs(b["energy_balance_residual_per_packet"])for b in bs),
          "weighted_launch_mean":float(np.mean([b["weighted_launch_mean"]for b in bs])),
          "weighted_launch_mc_standard_error":se([b["weighted_launch_mean"]for b in bs]),
          "weighted_escape_fraction_mean":float(np.mean([b["weighted_escape_fraction"]for b in bs])),
          "weighted_upward_escape_fraction_mean":float(np.mean([b["weighted_upward_escape_fraction"]for b in bs])) if b0["total_tau"] else .5,
          "weighted_upward_escape_fraction_mc_se":se([b["weighted_upward_escape_fraction"]for b in bs]) if b0["total_tau"] else 0.,
          "escape_statistics_method":"weighted_MC" if b0["total_tau"] else "analytic_isotropic_vacuum",
          "mean_scattering_events_per_packet":float(np.mean([b["mean_collisions_per_packet"]for b in bs])),
          "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),
          "response_overflow_relative_to_vacuum":float(avg[-1]),
          "run_seconds":sum(b["elapsed_s"]for b in bs),
          "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[]),
          "batch_info_file":"production/batch_info/"+cid+".json",
          "batch_kernel_file":"production/batch_kernels/"+cid+".npz"}
        if b0["total_tau"]:
            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"]}
        # Explicitcenter rule: directatom att0, MCbins at2n+1us.
        masses=np.zeros((len(raw),2*nb+1));masses[:,0]=direct;masses[:,1:2*nb:2]=raw[:,:-1]
        mass=masses.mean(axis=0)
        pulse_curves={};timing=[]
        for rise in [30.,80.,200.]:
          for fc in [0.,10000.,50000.]:
            t,y=tm.propagate_waveform(mass/fine,fine,rise,None if not fc else fc)
            met=tm.waveform_metrics(t,y)
            batchmets=[]
            for bmass in masses:
                tb,yb=tm.propagate_waveform(bmass/fine,fine,rise,None if not fc else fc)
                batchmets.append(tm.waveform_metrics(tb,yb))
            item={"intrinsic_rise_10_90_us":rise,"analog_one_pole_fc_hz":fc,
              "rise_10_90_us":met["rise_10_90_us"],"onset_to_peak_us":met["onset_to_peak_us"],
              "rise_batch_standard_error_us":se([x["rise_10_90_us"]for x in batchmets]),
              "onset_peak_batch_standard_error_us":se([x["onset_to_peak_us"]for x in batchmets]),
              "batch_rise_10_90_us":[x["rise_10_90_us"]for x in batchmets] if b0["total_tau"] else [],
              "batch_onset_to_peak_us":[x["onset_to_peak_us"]for x in batchmets] if b0["total_tau"] else [],
              "peak_gain_relative_to_vacuum_unit_peak_source":met["peak_amplitude"],
              "peak_at_window_end":met["peak_at_window_end"]}
            timing.append(item);timing_rows.append({"case_id":cid,**{k:v for k,v in item.items() if not isinstance(v,list)}})
            pulse_curves[f"rise{int(rise)}_fc{int(fc)}"]=(y[curve_idx]/y.max()).tolist()
        row=dict(small);row["timing"]=timing
        row["scattered_bin_mass_relative"]=coarse.mean(axis=0).tolist()
        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()
        row["pulse_curves_unit_peak"]=pulse_curves
        row["batch_total_gain"]=[b["total_relative_to_vacuum"]for b in bs] if b0["total_tau"] else []
        row["batch_upward_escape_fraction"]=[b["weighted_upward_escape_fraction"]for b in bs] if b0["total_tau"] else []
        cases.append(row);summary.append(small)
    output={"schema_version":"1.0","status":"EXECUTED_INDEPENDENT_SCALAR_CONDITIONAL_MC_NOT_REAL_CLOUD_VALIDATION",
      "generated_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),
      "model":{"geometry":m["geometry"],"source_altitude_km":m["source_altitude_km"],
        "total_tau_grid":m["total_tau"],"phases":["ghm","hg"],"physical_source_emission":"Isotropic; exact source-only importance weights used internally below the cloud top.",
        "medium":"Finite homogeneous cylinder; vacuum exterior; scalar conservative SSA=1. No Rayleigh scattering, gas absorption, aerosol, ground reflection, polarization or real-cloud heterogeneity.",
        "wavelength_label":"337.1nm numerical input coordinate; linear interpolation of native330/340nm data. N2 band integration and air/vacuum wavelength registration are not modeled.",
        "phase_input":m["phase_input"],
        "run_elapsed_s":m["elapsed_s"],
        "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.",
        "earth_geometry":"Locally flat cloud boundaries and straight vacuum rays; no Earth curvature, atmospheric refraction or finite source size.",
        "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.",
        "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"],
        "independent_batches_per_nonvacuum_case":m["independent_batches"],
        "engine_sha256_at_run":m["engine_sha256"],"runner_sha256_at_run":m["runner_sha256"],
        "engine_file":"cylinder_mc.py","runner_file":"run_grid.py","raw_LUT_included":False},
      "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.",
        "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.",
        "tau_total":"Extinction coefficient times total3km cloud thickness.",
        "tau_overlying":"Extinction along the straight vertical source-observer line; sources at/above the top have zero overlying cloud optical depth.",
        "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.",
        "time":"Excess physical path delay relative to unscattered source-observer flight. The direct atom at delay0 is supplied separately.",
        "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.",
        "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.",
        "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.",
        "electronics":"Ideal first-order low-pass only; not an actual analog or antialias chain. Cutoff frequency is distinct from sampling rate.",
        "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."},
      "time_grid":{"native_scattered_bin_width_us":dt*1e6,"response_window_us":m["response_window_s"]*1e6,
        "scattered_plot_bin_edges_us":edges_us.tolist(),"pulse_plot_time_us":curve_t_us.tolist(),
        "pulse_plot_window":"First12ms plotted; source and convolution metrics use an approximately48ms record.",
        "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."},
      "source_parameters":[tm.source_parameters(x)for x in [30.,80.,200.]],
      "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.",
      "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"},
      "cases":cases,
      "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."],
      "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"]}
    (HERE/"cloud_results.json").write_text(json.dumps(tidy(output),separators=(",",":"),allow_nan=False),encoding="utf-8")
    # CompactscalarCSV; nested diagnosticfields remain inJSON.
    scalar=[{k:v for k,v in x.items()if not isinstance(v,(dict,list))}for x in summary]
    with(HERE/"cloud_summary.csv").open("w",newline="",encoding="utf-8")as f:
        w=csv.DictWriter(f,fieldnames=list(scalar[0]));w.writeheader();w.writerows(scalar)
    with(HERE/"cloud_timing.csv").open("w",newline="",encoding="utf-8")as f:
        w=csv.DictWriter(f,fieldnames=list(timing_rows[0]));w.writeheader();w.writerows(timing_rows)
    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)}))
if __name__=="__main__":main()