Review index · Download raw source
reference/statistics/assemble_statistics.py
1from pathlib import Path2import json, csv, math, hashlib, platform, datetime, shutil, sys3import numpy as np4P=Path(__file__).resolve().parent5sys.path.insert(0,"D:/Photometer_research_20261022/research/cloud_transport/private_readers")6import scipy7from scipy.stats import beta8import practical_morphology as pm9def read(name):10 return json.loads((P/name).read_text(encoding="utf-8-sig"))11pr=read("practical_morphology_results.json")12oc=read("oracle_current_generator_results.json")13orr=read("oracle_reference/calculation_summary.json")14counts=read("count_inference/count_inference_results.json")15review=read("count_inference/morphology_profile_independent_review.json")16checks=[]17assert len(pr["scenarios"])==8118expected={(n,b,k,r) for n in [100,1000,10000] for b in [50.,1e4,1e6] for k in [10,50,100] for r in [30,80,200]}19actual={(r["N_plateau_equivalent_3ms_PE"],r["background_cps"],r["sample_rate_khz"],r["truth_receiver_rise_10_90_us"]) for r in pr["scenarios"]}20assert expected==actual21checks.append({"test":"complete81_scenario_cartesian_grid","passed":True})22synthetic_total=023for r in pr["scenarios"]:24 n=r["trials"]25 assert n>=100026 assert sum(r["three_class"]["confusion_counts"])==n27 assert sum(r["fast_slow"]["confusion_counts"])==n28 assert r["detection"]["count"]==n-r["three_class"]["confusion_counts"][-1]29 assert r["detection"]["count"]==n-r["fast_slow"]["confusion_counts"][-1]30 assert abs(r["expected_source_in_gate_PE_mean"]-r["actual_total_source_PE"])<1e-831 tr=np.load(P/"synthetic_trial_summaries"/f"case_{r['case_id']:03d}.npz")32 assert len(tr["detected"])==n33 assert int(tr["detected"].sum())==r["detection"]["count"]34 assert np.all(tr["detected"]==(tr["detection_pvalue"]<=pm.ALPHA_PER_TEST))35 assert np.all(tr["accept3"]==(tr["detected"]&(tr["margin3"]>=pm.MARGIN)))36 assert np.all(tr["accept2"]==(tr["detected"]&(tr["margin2"]>=pm.MARGIN)))37 assert np.all(np.isfinite(tr["scores"]))38 assert np.all(tr["scores"]>=-1e-7)39 if r["case_id"] in [0,40,80]:40 d=pm.setup(r["sample_rate_khz"])41 draw=pm.draw_case(d,r["N_plateau_equivalent_3ms_PE"],r["background_cps"],42 r["truth_receiver_rise_10_90_us"],r["seed"],n)43 assert np.array_equal(draw["pre"],tr["pre_count"])44 assert np.array_equal(draw["total"],tr["gate_count"])45 assert np.array_equal(draw["detected"],tr["detected"])46 assert np.array_equal(draw["onset"],tr["onset_s"])47 synthetic_total+=n48checks.append({"test":"all81_confusion_and_acceptance_matrices_and_source_integrals","passed":True})49checks.append({"test":"all81_saved_synthetic_trial_summaries_match_reported_rules","passed":True})50checks.append({"test":"three_scenarios_reproduced_exactly_from_seeds","cases":[0,40,80],"passed":True})51assert synthetic_total==8100052assert all(c["passed"] for c in counts["deterministic_tests"])53assert counts["coverage_summary"]["MC_family99_all_overlap_analytic"]54assert review["target_sha256"]==hashlib.sha256((P/"practical_morphology.py").read_bytes()).hexdigest()55assert all(r["byte_identical"] for r in oc["oracle_rerun"]["csv_matches"])56checks.append({"test":"exact_count_controls_and_family_MC_audit_pass","passed":True})57checks.append({"test":"independent_profile_review_targets_current_unchanged_script","passed":True})58checks.append({"test":"five_original_oracle_CSVs_match_revision_rerun","passed":True})59# Add exact individual binomial MC intervals to the consolidated JSON only.60# These supplement the original Wilson score intervals; not a simultaneous claim.61def add_cp(obj):62 if isinstance(obj,dict):63 if {"count","denominator","fraction","individual_95pct_Wilson_interval"}<=set(obj):64 k=obj["count"];n=obj["denominator"]65 obj["individual_95pct_Clopper_Pearson_interval"] = (66 [float(beta.ppf(.025,k,n-k+1)) if k else 0.,67 float(beta.ppf(.975,k+1,n-k)) if k<n else 1.] if n else [None,None])68 for v in list(obj.values()):add_cp(v)69 elif isinstance(obj,list):70 for v in obj:add_cp(v)71add_cp(pr)72refs=[73 {"id":"DERIVATION","title":"Own derivations, exact controls and independent audit","url":None,74 "locator":"README.md; count_inference/README.md; count_inference/morphology_profile_independent_review.json",75 "role":"Mathematical basis, not a new experimental source."},76 {"id":"SCIPY_POISSON","title":"SciPy official Poisson distribution API","url":"https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.poisson.html",77 "locator":"Methods sf/logsf; Notes PMF","checked":"2026-10-22",78 "role":"Numerical implementation documentation; web manual1.18.0, installed library1.18.1."},79 {"id":"SCIPY_BINOM","title":"SciPy official binomial distribution API","url":"https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.html",80 "locator":"Methods sf; Notes binomial PMF","checked":"2026-10-22","role":"On/off conditional test implementation."},81 {"id":"SCIPY_CHI2","title":"SciPy official chi-square distribution API","url":"https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html",82 "locator":"Methods ppf/isf","checked":"2026-10-22","role":"Poisson interval inversion."},83 {"id":"NUMPY_POISSON","title":"NumPy official Generator.poisson API","url":"https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.poisson.html",84 "locator":"Generator.poisson; Notes PMF","checked":"2026-10-22",85 "role":"Synthetic exact Poisson variate generator; current web manual differs from installed2.3.5."},86 {"id":"MMIA_HERITAGE","title":"The Modular Multispectral Imaging Array (MMIA) of the ASIM payload","url":"https://link.springer.com/article/10.1007/s11214-019-0593-y",87 "locator":"Earlier primary heritage audit in research/detectors_optics/mmia_optical_heritage/optical_heritage_review.json",88 "role":"Instrument context only. This prototype is not a reanalysis of its proprietary event waveforms."}]89headline=[r for r in pr["scenarios"] if r["sample_rate_khz"]==50 and r["background_cps"] in [50.,1e6]]90limitations=[91 "This is a conditional synthetic receiver-morphology study, not a flight sensitivity/completeness validation.",92 "The classes mean30/80/200us10-to90% rises at a hypothetical linear receiver input after transport; they do not label source depth.",93 "Source brightness, cloud transfer, line spectrum, detector PDE and source occurrence population are outside this count-level calculation.",94 "The practical code has an external coarse onset bracket and a specified event gate. Blind stream triggering and a mission false-alarm rate are unvalidated.",95 "A fixed ln19 profile-score margin is not calibrated95% confidence. Always show accepted fraction, rejected fraction and the error interval together.",96 "The3-class fast-subtype distinction is harder than the grouped fast-vs-slow test; these are separate decisions.",97 "Only exact model-matched templates and constant backgrounds are generated. Real independently labeled optical/radio controls are still needed.",98 "Post-counting electronic filtering, detector excess noise, afterpulses, saturation and deadtime are not simulated.17ns occupancy is an illustrative hardware flag, not Poisson restoration.",99 "Known-background red limits are invalid if a noisy finite baseline is silently substituted. The provided alternative uses a conservative Bonferroni construction.",100 "Red limits constrain expected registered counts in a specified window, not thermal-source absence. Selection, red-window timing, throughput and population ratios remain physical inputs.",101 "A useful mission event rate is NOT COMPUTABLE without an absolute received-fluence distribution, completeness and exposure definition.",102 "No empirical wide-FOV event rate is inferred from the conditional geometric/population proxy; source geography and coverage remain unknown."]103out={104 "schema_version":"1.0","revision":"20261022-r4","title":"Photon statistics: detection, practical morphology, and red-channel upper limits",105 "status":"EXECUTED_CONDITIONAL_STATISTICAL_PROTOTYPES_NOT_MISSION_VALIDATION",106 "assembled_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),107 "question_tags":["A2","C8","C9","E16","E17","G21","G22"],108 "software":{"python":platform.python_version(),"numpy":np.__version__,"scipy":scipy.__version__,109 "random_generator":"NumPy default_rng PCG64","no_new_installation":True},110 "summary":{"signal_scenarios":81,"signal_trials_each":1000,"independent_signal_trials":81000,111 "noise_only_trials":9000,"practical_runtime_after_imports_s":pr["runtime_s"],112 "old_oracle_result_CSVs_byte_identical":5,113 "same_generator_oracle_balanced_pairs":54,114 "independent_profile_objective_comparisons":review["comparison_count"],115 "exact_count_deterministic_controls":len(counts["deterministic_tests"]),116 "interpretation":"Detecting the3ms pulse does not ensure enough early photons to classify its rise.",117 "mission_useful_rate":"NOT_COMPUTABLE_FROM_AVAILABLE_INPUTS"},118 "practical_morphology":pr,"oracle_reference":orr,"oracle_current_generator":oc,119 "exact_count_inference":counts,120 "independent_profile_QA":{"status":review["status"],"comparison_count":review["comparison_count"],121 "maximum_absolute_LLR_error":review["maximum_absolute_LLR_error"],122 "maximum_absolute_fitted_total_error":review["maximum_absolute_fitted_total_error"],123 "file":"count_inference/morphology_profile_independent_review.json"},124 "QA":{"passed":True,"checks":checks,"note":"Mathematical and software QA only; not physical validation."},125 "headline50kHz_examples":headline,126 "limitations":limitations,"source_refs":refs,127 "files":{"practical_scenarios_csv":"practical_morphology_scenarios.csv",128 "compact_alias_csv":"statistics_scenarios.csv",129 "practical_full_json":"practical_morphology_results.json",130 "current_oracle_json":"oracle_current_generator_results.json",131 "count_inference_json":"count_inference/count_inference_results.json",132 "synthetic_per_trial_arrays":"synthetic_trial_summaries/case_000.npz through case_080.npz",133 "scripts":["practical_morphology.py","oracle_same_generator.py","oracle_reference/calculate_ideal_poisson_rise.py",134 "count_inference/count_inference.py","count_inference/review_morphology_profile.py",135 "assemble_statistics.py"]}}136(P/"statistics_results.json").write_text(json.dumps(out,indent=2,allow_nan=False),encoding="utf-8")137(P/"statistics_QA.json").write_text(json.dumps(out["QA"],indent=2),encoding="utf-8")138shutil.copyfile(P/"practical_morphology_scenarios.csv",P/"statistics_scenarios.csv")139(P/"source_refs.json").write_text(json.dumps(refs,indent=2),encoding="utf-8")140manifest=[]141for f in sorted(P.rglob("*")):142 if f.is_file() and "__pycache__" not in f.parts and f.name not in ["CHECKPOINT.json","manifest.json"]:143 manifest.append({"file":f.relative_to(P).as_posix(),"bytes":f.stat().st_size,144 "sha256":hashlib.sha256(f.read_bytes()).hexdigest()})145(P/"manifest.json").write_text(json.dumps({"files":manifest},indent=2),encoding="utf-8")146checkpoint={"status":out["status"],"finished_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),147 "results":"statistics_results.json","QA":"statistics_QA.json",148 "manifest":"manifest.json","manifest_sha256":hashlib.sha256((P/"manifest.json").read_bytes()).hexdigest(),149 "file_count":len(manifest),"total_bytes":sum(f["bytes"] for f in manifest),150 "preserved_originals":"No prior research, public site, PDFs or baseline outputs changed.",151 "numerical_scope":"Synthetic statistical calculations only; no licensed or proprietary experimental trace republished.",152 "mission_rate":"NOT_COMPUTABLE","QA_passed":True}153(P/"CHECKPOINT.json").write_text(json.dumps(checkpoint,indent=2),encoding="utf-8")154print(json.dumps(checkpoint,indent=2),flush=True)155156