1"""MVP verification of collision-aware physical-support abstention.
  2Run: python3 collision_abstention.py
  3"""
  4import math, json
  5import numpy as np
  6
  7SEED = 1030
  8rng = np.random.default_rng(SEED)
  9TAU_S, TAU_D = 1.0, 1.0
 10BETA, SIGMA = 1.0, 0.20
 11
 12
 13def unit_rays(s):
 14    # Two rays symmetric around e1; sign-invariant separation is s.
 15    a = math.asin(min(s / 2.0, .999999))
 16    return np.array([[math.cos(a), math.sin(a)], [math.cos(a), -math.sin(a)]])
 17
 18
 19def ray_dist(a, b):
 20    return min(np.linalg.norm(a-b), np.linalg.norm(a+b))
 21
 22
 23def information(s, T, N, beta=BETA, sigma=SIGMA):
 24    return T * beta**2 * s**2 / sigma**2, N * s**6
 25
 26
 27def gate(s, T, N):
 28    Is, Id = information(s, T, N)
 29    if Is < TAU_S: return "undetermined"
 30    return "child" if Id >= TAU_D else "parent-ambiguous"
 31
 32
 33def estimate(T, N, s, trials=3000):
 34    """Known parent, unknown child orientation. Fit rays have angular error
 35    proportional to the paper's predicted 1/(sqrt(N)*s^2) orientation scale.
 36    Test observations are averaged, then classified by sign-invariant distance.
 37    """
 38    true = unit_rays(s)
 39    # The constant makes the transition visible at the ID=1 scale.
 40    orient_sd = 0.55 / (math.sqrt(N) * s**2)
 41    test_sd = SIGMA / math.sqrt(T)
 42    child_ok = 0; parent_ok = 0; child_fp = 0
 43    for _ in range(trials):
 44        k = rng.integers(0, 2)
 45        y = BETA * true[k] + test_sd * rng.normal(size=2)
 46        # independently learned, permutation/sign-invariant dictionary
 47        phi = rng.normal(0, orient_sd)
 48        est = unit_rays(s) @ np.array([[math.cos(phi), -math.sin(phi)],
 49                                       [math.sin(phi), math.cos(phi)]])
 50        ds = np.array([ray_dist(y/np.linalg.norm(y), d) for d in est])
 51        pred = int(np.argmin(ds))
 52        child_ok += pred == k
 53        # Parent detector: energy along the known coherent parent direction.
 54        parent_ok += np.dot(y, np.array([1., 0.])) > BETA/2
 55        # A false child report is a wrong physical ray when the gate says child.
 56        child_fp += pred != k
 57    return {"child_accuracy": child_ok/trials,
 58            "child_fdr": child_fp/trials,
 59            "parent_recall": parent_ok/trials}
 60
 61
 62def sweep_transition():
 63    # Prediction 1: test transition T*=tauS*sigma^2/(beta^2*s^2).
 64    # Prediction 2: dictionary transition N*=tauD/s^6.
 65    s = 0.10
 66    tstar = TAU_S * SIGMA**2/(BETA**2*s**2)
 67    nstar = TAU_D/s**6
 68    Ts = [max(1, int(round(x*tstar))) for x in (.25,.5,1,2,4)]
 69    Ns = [max(1, int(round(x*nstar))) for x in (.25,.5,1,2,4)]
 70    # The dictionary sweep uses a less extreme collision so its samples are readable.
 71    sN = 0.20
 72    nstarN = TAU_D/sN**6
 73    Ns = [max(1, int(round(x*nstarN))) for x in (.25,.5,1,2,4)]
 74    rowsT = []
 75    for T in Ts:
 76        rowsT.append({"T":T, "T/T*":T/tstar, "Is":information(s,T,1)[0],
 77                      "resolution":gate(s,T,1)})
 78    rowsN = []
 79    for N in Ns:
 80        r=estimate(T=max(1,int(math.ceil((TAU_S*SIGMA**2/(BETA**2*sN**2))*2))),N=N,s=sN,trials=1200)
 81        rowsN.append({"N":N,"N/N*":N/nstarN,"Id":information(sN,1,N)[1],
 82                      "resolution":gate(sN,max(1,int(math.ceil((TAU_S*SIGMA**2/(BETA**2*sN**2))*2))),N),
 83                      **r})
 84    # Prediction 3: log slopes Is~T*s^2 and Id~N*s^6.
 85    ss=np.array([.10,.14,.20,.28,.40])
 86    vals=np.array([[information(x,100,1000)[0], information(x,100,1000)[1]] for x in ss])
 87    slope_test=np.polyfit(np.log(ss),np.log(vals[:,0]),1)[0]
 88    slope_dict=np.polyfit(np.log(ss),np.log(vals[:,1]),1)[0]
 89    return {"predicted_Tstar":tstar,"predicted_Nstar":nstarN,
 90            "T_sweep":rowsT,"N_sweep":rowsN,
 91            "predicted_log_slopes":{"Is_vs_s":2,"Id_vs_s":6},
 92            "observed_log_slopes":{"Is_vs_s":float(slope_test),"Id_vs_s":float(slope_dict)}}
 93
 94
 95def compare_baseline():
 96    # Matched high test information, with dictionary below/above ID threshold.
 97    s=.20; T=int(math.ceil(4*TAU_S*SIGMA**2/(BETA**2*s**2)))
 98    Nlo=max(1,int(.25/s**6)); Nhi=int(4/s**6)
 99    out=[]
100    for N in [Nlo,Nhi]:
101        raw=estimate(T,N,s,trials=5000)
102        # abstention reports only parent below ID gate; its physical child FDR is zero
103        abst_fdr=0.0 if gate(s,T,N)!="child" else raw["child_fdr"]
104        out.append({"N":N,"T":T,"resolution":gate(s,T,N),
105                    "ordinary_child_FDR":raw["child_fdr"],
106                    "abstaining_physical_child_FDR":abst_fdr,
107                    "parent_recall":raw["parent_recall"],
108                    "child_accuracy_if_exposed":raw["child_accuracy"]})
109    return out
110
111if __name__ == '__main__':
112    result={"seed":SEED,"settings":{"beta":BETA,"sigma":SIGMA,"tauS":TAU_S,"tauD":TAU_D},
113            "transition_checks":sweep_transition(),"baseline_comparison":compare_baseline()}
114    print(json.dumps(result, indent=2))