Review index · Download raw source
reference/statistics/count_inference/review_morphology_profile.py
1"""Read-only independent boundary check of parent morphology profile likelihood."""2import sys3sys.dont_write_bytecode = True4import pathlib, json, hashlib, importlib.util, datetime5import numpy as np6sys.path.insert(0,"D:/Photometer_research_20261022/research/cloud_transport/private_readers")7from scipy.optimize import minimize_scalar8from scipy.special import xlogy910P=pathlib.Path(__file__).resolve().parent11SOURCE=P.parent/"practical_morphology.py"1213def main():14 raw=SOURCE.read_bytes()15 spec=importlib.util.spec_from_file_location("profile_review_target",SOURCE)16 module=importlib.util.module_from_spec(spec)17 spec.loader.exec_module(module) # __main__ is not called; no target writes.18 rng=np.random.default_rng(202610220817)19 reviews=[]20 for khz in [10,50,100]:21 d=module.setup(khz)22 for rise in module.RISES_US:23 selected=[m for m in d["models"] if m["rise_us"]==rise]24 for model in [selected[0],selected[len(selected)//2],selected[-1]]:25 q=model["q"]; dur=d["durations"]; Q=q.sum(); D=dur.sum()26 h=q/Q; w=dur/D27 cases=[("zero_counts",np.zeros(len(q),dtype=int))]28 c=np.zeros(len(q),dtype=int);c[-1]=1729 cases.append(("off_only_theta0",c))30 c=np.zeros(len(q),dtype=int);c[-2]=1731 cases.append(("plateau_only_theta1",c))32 for A,b in [(0,50),(0,1e6),(1e3,0),(3.333333333e4,50),33 (3.333333333e6,0),(3.333333333e6,1e6)]:34 cases.append((f"A{A}_b{b}",rng.poisson(A*q+b*dur)))35 counts=np.stack([z[1] for z in cases])36 ll,theta,Ahat,bhat=module.fit_one(counts,model)37 for j,(name,n) in enumerate(cases):38 N=int(n.sum())39 # Direct original-bin mixture objective, no zero/full/active compression.40 def direct_score(t):41 return float(np.sum(xlogy(n,(1-t)*w+t*h)-xlogy(n,w)))42 opt=minimize_scalar(lambda t:-direct_score(t),bounds=(0,1),method="bounded",43 options={"xatol":1e-13,"maxiter":1000})44 candidates=[(0.,direct_score(0.)),(1.,direct_score(1.)),(float(opt.x),-float(opt.fun))]45 td,ld=max(candidates,key=lambda x:x[1])46 fitted=Ahat[j]*q+bhat[j]*dur47 scale_error=abs(float(fitted.sum())-N)48 score_error=abs(float(ll[j])-ld)49 if score_error>3e-6 or scale_error>1e-7:50 raise AssertionError((khz,rise,model["onset_s"],name,score_error,scale_error))51 if name=="off_only_theta0" and theta[j]!=0:52 raise AssertionError("off-only boundary should theta=0")53 if name=="plateau_only_theta1" and theta[j]!=1:54 raise AssertionError("plateau-only boundary should theta=1")55 reviews.append({"sample_rate_khz":khz,"rise_us":rise,56 "onset_s":model["onset_s"],"case":name,57 "observed_count":N,"theta_target":float(theta[j]),58 "theta_direct":td,"absolute_LLR_error":score_error,59 "absolute_fitted_total_error":scale_error})60 out={61 "status":"INDEPENDENT_READ_ONLY_PROFILE_DERIVATION_AND_BOUNDARY_QA_PASSED",62 "finished_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),63 "target_path":str(SOURCE),"target_sha256":hashlib.sha256(raw).hexdigest(),64 "target_unchanged_during_check":raw==SOURCE.read_bytes(),65 "seed":202610220817,66 "derivation":[67 "lambda_j=A*q_j+b*d_j, Q=sum(q),D=sum(d). Set eta=AQ+bD and theta=AQ/eta.",68 "Normalized h=q/Q,w=d/D sum to1. logL=Nlog(eta)-eta+sum nlog((1-theta)w+theta*h)+constant.",69 "eta_hat=N for N>0. At N=0 theta is unidentifiable and A=b=0 suffices.",70 "Relative LLR=sum n*log1p(theta*(h/w-1)); this is concave in theta on[0,1].",71 "A_hat=theta*N/Q,b_hat=(1-theta)*N/D. Root zero/full/active grouping gives this same objective.",72 "Boundary theta=1 requires no counts in bins with h=0. Off-only theta0 and plateau-only theta1 controls were exercised."],73 "comparison":"Direct original-bin objective using scipy.special.xlogy, separate scalar bounded optimization and explicit endpoints; no grouped derivative implementation reused.",74 "comparison_count":len(reviews),75 "maximum_absolute_LLR_error":max(z["absolute_LLR_error"] for z in reviews),76 "maximum_absolute_fitted_total_error":max(z["absolute_fitted_total_error"] for z in reviews),77 "rows":reviews,78 "qualifications":[79 "Read-only target import; its __main__/run functions were never called, no target outputs changed.",80 "Checks establish algebra/numerical implementation for exercised cases, not source-population or flight validation.",81 "Morphology model fits only designated early bins, aggregate late plateau and prewindow. Detection uses the full gate separately.",82 "ln19 profile-score margin has no automatic95% confidence/posterior interpretation; target output states that limit.",83 "Data-dependent onset/gate searches, mismatch, analog correlated response and deadtime need separate treatment.",84 "Zero source MLE leaves onset/morphology unidentified; target uses NaN onset at theta=0."85 ],86 "writes_only_own_subfolder":True87 }88 (P/"morphology_profile_independent_review.json").write_text(json.dumps(out,indent=2)+"\n",encoding="utf-8")89 print(json.dumps({k:out[k] for k in ["status","comparison_count","maximum_absolute_LLR_error","maximum_absolute_fitted_total_error","target_unchanged_during_check"]}))90if __name__=="__main__":91 main()92