"""Run pinned independent finite-cylinder grid; resumable per-case outputs."""
from pathlib import Path
import sys,json,hashlib,time,datetime,argparse
import numpy as np
from cylinder_mc import Geometry,load_phases,simulate,phase_qa,intervals,segment_inside,sha
def main():
    ap=argparse.ArgumentParser()
    ap.add_argument("--lut",required=True);ap.add_argument("--out",required=True)
    ap.add_argument("--packets",type=int,default=50000);ap.add_argument("--batches",type=int,default=6)
    args=ap.parse_args();out=Path(args.out);out.mkdir(exist_ok=True)
    for name in ["batch_info","batch_kernels"]: (out/name).mkdir(exist_ok=True)
    phases,meta=load_phases(args.lut);geo=Geometry()
    q={"phase_input":meta,"phases":{k:phase_qa(v)for k,v in phases.items()}}
    q["geometry_tests"]=[]
    tests=[([0,0,15000],[0,0,1],0,2000,True),
      ([0,0,15000],[0,0,-1],0,1000,True),
      ([0,0,15000],[1,0,0],0,30000,True),
      ([0,0,18000],[0,0,-1],1000,4000,True),
      ([0,0,18000],[0,0,1],None,None,False),
      ([0,0,17000],[0,0,1],None,None,False),
      ([0,0,17000],[0,0,-1],0,3000,True),
      ([40000,0,15000],[-1,0,0],10000,70000,True)]
    for pos,vec,exlo,exhi,hit in tests:
        lo,hi,h=intervals(np.array([pos],float),np.array([vec],float),geo)
        passed=bool(h[0])==hit and(not hit or(abs(lo[0]-exlo)<1e-8 and abs(hi[0]-exhi)<1e-8))
        q["geometry_tests"].append({"position_m":pos,"direction":vec,"passed":bool(passed)})
    q["all_geometry_tests_pass"]=all(x["passed"]for x in q["geometry_tests"])
    if not q["all_geometry_tests_pass"] or not all(x["pass"]for x in q["phases"].values()):raise RuntimeError("Preflight QA failed")
    (out/"phase_and_geometry_QA.json").write_text(json.dumps(q,indent=2),encoding="utf-8")
    manifest={"status":"RUNNING","started_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),
      "packets_per_batch":args.packets,"independent_batches":args.batches,
      "total_tau":[0,.3,1,3,10,30,45],"source_altitude_km":[15,16,17,18],
      "geometry":geo.__dict__,"dt_s":2e-6,"response_window_s":.012,"max_scatter":20000,
      "phase_input":meta,"engine_sha256":sha(Path(__file__).parent/"cylinder_mc.py"),
      "runner_sha256":sha(__file__),"case_ids":[]}
    manifest_path=out/"run_manifest.json"
    manifest_path.write_text(json.dumps(manifest,indent=2),encoding="utf-8")
    start=time.perf_counter()
    for phase_i,(phase_name,phase)in enumerate(phases.items()):
      for alt_i,source in enumerate([15000.,16000.,17000.,18000.]):
       for ti,tau in enumerate([0.,.3,1.,3.,10.,30.,45.]):
        case=f"{phase_name}_z{int(source/1000)}_tau{str(tau).replace('.','p')}"
        info_path=out/"batch_info"/(case+".json");kernel_path=out/"batch_kernels"/(case+".npz")
        if info_path.exists()and kernel_path.exists():
            print("REUSE",case,flush=True)
        else:
            infos=[];hists=[]
            for bi in range(args.batches if tau else 1):
                seed=2026102200+phase_i*100000+alt_i*10000+ti*100+bi
                n=args.packets if tau else 1
                info,h=simulate(phase,source,tau,n,seed,geo,importance=.5)
                info["phase_id"]=phase_name
                infos.append(info);hists.append(h)
            info_path.write_text(json.dumps({"case_id":case,"geometry":geo.__dict__,"batches":infos},indent=2),encoding="utf-8")
            np.savez_compressed(kernel_path,batch_mass_relative=np.array(hists),dt_s=np.array(2e-6),window_s=np.array(.012))
            a=np.array([x["total_relative_to_vacuum"]for x in infos])
            se=float(a.std(ddof=1)/np.sqrt(len(a)))if len(a)>1 else 0.
            print(json.dumps({"case":case,"mean_gain":float(a.mean()),"batch_SE":se,"seconds":sum(x["elapsed_s"]for x in infos)}),flush=True)
        manifest["case_ids"].append(case)
        manifest["elapsed_s"]=time.perf_counter()-start
        manifest_path.write_text(json.dumps(manifest,indent=2),encoding="utf-8")
    manifest["status"]="COMPLETE";manifest["ended_utc"]=datetime.datetime.now(datetime.timezone.utc).isoformat()
    manifest["elapsed_s"]=time.perf_counter()-start
    manifest["nonvacuum_total_packets"]=2*4*6*args.batches*args.packets
    manifest_path.write_text(json.dumps(manifest,indent=2),encoding="utf-8")
    print("COMPLETE",manifest["elapsed_s"],flush=True)
if __name__=="__main__":main()
