Building the benchmark harness — it forces the u/p/i/g operationalizations anyway, so it's both. Design before code, so the commitments are explicit: Three scenarios, one deliberately unfair to residual-based methods. **A: silent sensor freeze** at t=1000 (channel holds last value — numerically plausible, partially redundant, so prediction residuals barely move). **B: geometry gap** — the target gains dependence on an unobserved latent factor; residuals rise then plateau, the learner converges within its class. **C: stationary control** for false-alarm rates. Baseline is Page-Hinkley drift detection on residuals. Thresholds calibrated as percentiles of the control burn-in — a declared procedure, answering the free-threshold critique rather than dodging it. 10/10 on every scenario. Writing up the results honestly, then packaging: Done. Ten seeds, three scenarios, and the manuscript came out changed: **Headline result:** Scenario A is e≈0 ⇏ R≈0 as data. A frozen sensor with a redundant twin — residual-invisible by construction — caught 10/10 by the Remainder's instrumentation component, missed 4/10 by Page-Hinkley. Scenario B is the differentiation argument §10.3 was missing: PH detects faster but can only say *drift*; GEO-01b renders the disposition — tune during adaptation, promote to GEOMETRY_EXPANSION_CANDIDATE once improvement stalls. Zero false alarms and zero false expansions in control. **The correction that matters most:** §6.1's gradient-norm criterion failed structurally — 0/10 on the geometry-gap scenario it exists for, because SGD gradient norm tracks residual magnitude and never looks "low" under a persistent unmodeled factor. The working replacement is a loss-slope estimator with the low-gain threshold set at the *upper* range of stationary jitter. That's a §6.1 rewrite, not a footnote. Worth noting without ceremony: the disposition chatter was fixed by your own τ-gate — dwell on every transition, promotion and demotion. The criterion needed the runtime's primitive inside it. Open items are in the results doc — B's tail still demotes intermittently on some seeds, p_t is unexercised, and everything so far survives only a linear world. But the manuscript is no longer future-tense on its central claim. # GEO-01 Benchmark v0.1 — Results **Question:** Does the geometry-evaluation criterion (manuscript Appendix B / §6.1) detect faults that residual-based drift detection cannot — and does it distinguish *tune* from *revise*? **Setup:** Linear-Gaussian world, 3 latents, 6 sensors, online linear model (SGD). Baseline: Page-Hinkley on residuals. Thresholds calibrated as percentiles of a control burn-in, with a declared finite-calibration margin. 10 seeds. ## Scenarios - **A — Silent sensor freeze.** One channel holds its last value at t=1000. An exact-duplicate channel makes the freeze residual-invisible. - **B — Geometry gap.** Target gains dependence on an unobserved latent factor at t=1000. Residuals rise, then plateau at the model class's floor. - **C — Control.** Stationary. Measures false alarms. ## Results (10 seeds) | | GEO-01b detect | Page-Hinkley detect | |---|---|---| | **A** (freeze) | 10/10, delay 79–120 | **6/10**, delay 26–194 | | **B** (geometry gap) | 10/10 GEOMETRY_EXPANSION_CANDIDATE, delay 62–106 | 10/10 alarm, delay 9–66, **no disposition** | | **C** (control) | 0 false alarms, **0 false expansions** | 0 false alarms | Scenario A is the e≈0 ⇏ R≈0 claim as data: a fault invisible to the residual channel, caught every time by the Remainder's instrumentation component, missed 4/10 times by the residual-only baseline. Scenario B is the differentiation argument: PH detects change faster but can only say *drift*. GEO-01b renders the disposition the manuscript exists for — parameter updates during the adaptation phase (+30), promotion to expansion candidate once improvement stalls (+79), sustained thereafter. ## Findings that change the manuscript 1. **§6.1's ||∇_θL|| < τ_θ fails structurally.** SGD gradient norm is proportional to residual magnitude; under a persistent unmodeled factor it never looks "low," so the literal criterion never reaches GEOMETRY_EXPANSION_CANDIDATE on the one scenario it exists for (0/10). Replace with an improvement-slope estimator: relative decrease of rolling loss over a horizon, low gain = improving no faster than the upper range (75th pct) of stationary jitter. 2. **The τ-gate belongs in the criterion, not just the runtime.** Undamped, the disposition chattered 58/42 between tune and revise. Dwell-time gating on *every* disposition (promotion and demotion) stabilized it. Instrumentation needed the same gate: real freezes hold i_t high indefinitely; natural quiet stretches don't. 3. **i_t requires a declared statistic, and linear variance ratio is the wrong one.** Control-regime dips reach ratio ~0.2; a freeze reaches ~1e-8. Linear space gives poor separation and, combined with a threshold margin, made detection impossible on 6/10 seeds. Log-space variance ratio separates the regimes by ~7 orders of magnitude. 4. **Percentile calibration needs an effective-sample-size correction.** A 99th percentile from ~6 effective independent windows underestimates the tail; uncorrected, false-alarm rates hit 10% (including false expansions). A declared 1.25× margin plus dwell gating brought all false alarms to zero across seeds. The paper should state calibration length requirements relative to the autocorrelation time of the monitored statistic. ## Open items - Scenario B tail disposition holds GEOMETRY_EXPANSION 66–100% of post-plateau steps depending on seed. Demotion chatter is reduced, not eliminated. Candidate fix: asymmetric dwell (fast promote, slow demote) or Bayesian slope estimation. - p_t (provenance) not exercised — fixed at 0. Needs its own scenario (numerically plausible corrupted-lineage data). - u_t computed but not wired into any disposition. Either give it a decision role or drop it from v0.2's trigger. - Linear world, linear model. The claims survive one model class. Nonlinear replication is the obvious next tier. ## Files - `geo01_benchmark.py` — full harness, numpy-only, seeded, single command - `geo01_results.png` — traces and dispositions per scenario """ GEO-01 Benchmark v0.1 ===================== Tests the QOTE/Resona geometry-evaluation criterion (Appendix B) against a standard residual-based drift detector (Page-Hinkley) on synthetic scenarios where observable residuals are deceptive. Scenarios --------- A Sensor freeze: one input channel silently holds its last value at t=FAULT_T. Channel is partially redundant -> residuals barely move. B Geometry gap: target gains dependence on an unobserved latent factor at t=FAULT_T. Residuals rise then plateau; learner converges within its (now inadequate) model class. C Control: stationary. Measures false-alarm behavior. Remainder component operationalizations (declared, minimal) ------------------------------------------------------------ e_t : |z - z_hat| / rolling sigma of target (normalized residual) u_t : |log(rolling residual var / model noise var)| (calibration mismatch) p_t : provenance deficiency flag (not exercised here; fixed 0 — declared limitation) i_t : max over channels of variance-collapse score (1 - var_ratio, clipped) g_t : plateau proxy = elevated mean residual AND low gradient norm over window (inferred, not observed — consistent with manuscript §4.4 concession) Threshold calibration --------------------- All tau values set as the 99th percentile of the corresponding statistic over the control burn-in window. This is the declared calibration procedure. Dependencies: numpy only. """ import numpy as np RNG_SEED = 42 T_TOTAL = 3000 FAULT_T = 1000 BURN_IN = 800 # calibration window (pre-fault, also used for control) WIN = 100 # rolling window for monitors N_SENSORS = 6 LATENT_D = 3 LR = 0.02 # SGD learning rate for online linear model # ---------------------------------------------------------------------- # World + scenarios # ---------------------------------------------------------------------- def simulate(scenario: str, seed: int = RNG_SEED): """Returns sensors (T x k), target z (T,), and ground-truth fault time.""" rng = np.random.default_rng(seed) A = 0.95 x = np.zeros((T_TOTAL, LATENT_D)) for t in range(1, T_TOTAL): x[t] = A * x[t - 1] + rng.normal(0, 0.3, LATENT_D) # Sensors: linear reads of latents + noise. Sensors 0 and 1 both read # latent 0 (redundancy), so freezing sensor 0 is survivable for prediction. H = np.array([ [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], # exact duplicate of sensor 0 -> freeze is residual-invisible [0.0, 1.0, 0.0], [0.0, 0.8, 0.2], [0.0, 0.0, 1.0], [0.2, 0.0, 0.9], ]) sensors = x @ H.T + rng.normal(0, 0.1, (T_TOTAL, N_SENSORS)) # Target depends on all three latents w_true = np.array([1.0, -0.7, 0.5]) z = x @ w_true + rng.normal(0, 0.15, T_TOTAL) if scenario == "A": # silent sensor freeze sensors[FAULT_T:, 0] = sensors[FAULT_T - 1, 0] fault = FAULT_T elif scenario == "B": # geometry gap: new unobserved latent enters target x_new = np.zeros(T_TOTAL) for t in range(FAULT_T + 1, T_TOTAL): x_new[t] = 0.95 * x_new[t - 1] + rng.normal(0, 0.3) z = z + 0.9 * x_new # no sensor reads x_new fault = FAULT_T elif scenario == "C": # control fault = None else: raise ValueError(scenario) return sensors, z, fault # ---------------------------------------------------------------------- # Online model: linear regression via SGD (fixed geometry = linear in sensors) # ---------------------------------------------------------------------- def run_model(sensors, z): T, k = sensors.shape w = np.zeros(k) b = 0.0 preds = np.zeros(T) grads = np.zeros(T) for t in range(T): zh = sensors[t] @ w + b preds[t] = zh err = zh - z[t] g = np.concatenate([err * sensors[t], [err]]) grads[t] = np.linalg.norm(g) w -= LR * err * sensors[t] b -= LR * err return preds, grads # ---------------------------------------------------------------------- # Monitors # ---------------------------------------------------------------------- def rolling(arr, win, fn): out = np.full(len(arr), np.nan) for t in range(win, len(arr)): out[t] = fn(arr[t - win:t]) return out def compute_remainder(sensors, z, preds, grads): T = len(z) resid = z - preds abs_resid = np.abs(resid) sigma_z = rolling(z, WIN, np.std) e = abs_resid / np.where(sigma_z > 1e-8, sigma_z, 1.0) # u_t: calibration mismatch — rolling residual variance vs. burn-in residual variance base_var = np.var(resid[WIN:BURN_IN]) + 1e-8 rvar = rolling(resid, WIN, np.var) u = np.abs(np.log((rvar + 1e-8) / base_var)) # i_t: per-channel variance collapse vs burn-in, in log space. # A frozen channel drives its variance ratio toward machine epsilon # (-log10 ~ 8+), while natural quiet stretches of an AR(1) latent # bottom out around ratio ~ 0.2 (-log10 ~ 0.7). Log space separates # these by an order of magnitude; the linear ratio does not. i = np.zeros(T) base_ch_var = np.var(sensors[WIN:BURN_IN], axis=0) + 1e-12 for t in range(WIN, T): ch_var = np.var(sensors[t - WIN:t], axis=0) + 1e-12 ratio = (ch_var / base_ch_var).min() i[t] = np.clip(-np.log10(ratio) / 6.0, 0.0, 1.0) # g_t proxy: elevated rolling mean |resid| AND low rolling gradient norm mean_r = rolling(abs_resid, WIN, np.mean) mean_g = rolling(grads, WIN, np.mean) p = np.zeros(T) # provenance not exercised # improvement potential, slope version: relative decrease of rolling mean # residual over the last WIN steps. Near-zero or positive slope while # residual is elevated => parameter updates are no longer buying anything. H = 3 * WIN # longer horizon: slope over 300 steps, not 100 improve = np.full(T, np.nan) for t in range(WIN + H, T): prev, cur = mean_r[t - H], mean_r[t] if prev and prev > 1e-8: improve[t] = (prev - cur) / prev return {"e": e, "u": u, "i": i, "p": p, "mean_resid": mean_r, "mean_grad": mean_g, "improve": improve} def geo01(components, taus): """Appendix B evaluator, per timestep. Returns disposition string array.""" T = len(components["e"]) out = np.array(["OK"] * T, dtype=object) for t in range(BURN_IN, T): if components["i"][t] > taus["i"]: out[t] = "INSTRUMENTATION_FAILURE" continue elevated = components["mean_resid"][t] > taus["resid"] low_gain = components["mean_grad"][t] < taus["grad"] if elevated and not low_gain: out[t] = "PARAMETER_UPDATE_SUFFICIENT" elif elevated and low_gain: out[t] = "GEOMETRY_EXPANSION_CANDIDATE" return out def geo01b(components, taus, dwell=50): """Variant: improvement potential = loss slope, with dwell-time gating. Promotion to GEOMETRY_EXPANSION_CANDIDATE requires the low-gain+elevated condition to hold for `dwell` consecutive steps; demotion likewise requires sustained counter-evidence. (The manuscript's tau-gate, applied here.)""" T = len(components["e"]) out = np.array(["OK"] * T, dtype=object) pro, dem = 0, 0 instr_run = 0 in_geom = False for t in range(BURN_IN, T): instr_run = instr_run + 1 if components["i"][t] > taus["i"] else 0 if instr_run >= 30: # dwell: a real freeze holds i_t high indefinitely out[t] = "INSTRUMENTATION_FAILURE" continue elevated = components["mean_resid"][t] > taus["resid"] imp = components["improve"][t] low_gain = (not np.isnan(imp)) and imp < taus["improve"] cond = elevated and low_gain pro = pro + 1 if cond else 0 dem = dem + 1 if not cond else 0 if not in_geom and pro >= dwell: in_geom = True elif in_geom and dem >= dwell: in_geom = False if in_geom: out[t] = "GEOMETRY_EXPANSION_CANDIDATE" elif elevated: out[t] = "PARAMETER_UPDATE_SUFFICIENT" return out def calibrate(components, margin=1.25): """tau = extreme percentile of control-region statistics. `margin` is a finite-calibration correction: with WIN-step rolling statistics over a ~600-step calibration slice, there are only ~6 effective independent samples, so raw 99th percentiles underestimate the true tail. The multiplicative margin compensates; a longer calibration window would allow margin -> 1.""" sl = slice(WIN * 2, BURN_IN) return { "i": np.nanpercentile(components["i"][sl], 99) * margin, "resid": np.nanpercentile(components["mean_resid"][sl], 99) * margin, "grad": np.nanpercentile(components["mean_grad"][sl], 1), # slope threshold: 75th pct of control improvement jitter. # "low gain" = not improving faster than the upper range of # stationary noise; a true plateau (slope ~ 0) satisfies this # persistently, while active adaptation (steeply falling loss) # does not. "improve": np.nanpercentile(components["improve"][sl], 75), } # ---------------------------------------------------------------------- # Baseline: Page-Hinkley on absolute residuals # ---------------------------------------------------------------------- def page_hinkley(resid, delta=0.005, lam=None, calib_end=BURN_IN): x = np.abs(resid) # calibrate lambda: max PH statistic over control region * 1.5 def ph_stat(xs): mean = 0.0 mT = 0.0 m_min = 0.0 stats = np.zeros(len(xs)) for t, v in enumerate(xs): mean += (v - mean) / (t + 1) mT += v - mean - delta m_min = min(m_min, mT) stats[t] = mT - m_min return stats if lam is None: lam = ph_stat(x[:calib_end]).max() * 1.5 stats = ph_stat(x) alarms = stats > lam return alarms, lam # ---------------------------------------------------------------------- # Evaluation # ---------------------------------------------------------------------- def first_sustained(mask, t0, k=20): """First t >= t0 where mask holds for k consecutive steps.""" run = 0 for t in range(t0, len(mask)): run = run + 1 if mask[t] else 0 if run >= k: return t - k + 1 return None def run_scenario(name): sensors, z, fault = simulate(name) preds, grads = run_model(sensors, z) comps = compute_remainder(sensors, z, preds, grads) taus = calibrate(comps) dispositions = geo01(comps, taus) dispositions_b = geo01b(comps, taus) ph_alarms, lam = page_hinkley(z - preds) res = {"scenario": name, "fault_t": fault, "taus": taus, "ph_lambda": lam} if fault is not None: for tag, d in [("geo01", dispositions), ("geo01b", dispositions_b)]: res[f"{tag}_instr_detect"] = first_sustained(d == "INSTRUMENTATION_FAILURE", fault) res[f"{tag}_geom_detect"] = first_sustained(d == "GEOMETRY_EXPANSION_CANDIDATE", fault) res[f"{tag}_tune_detect"] = first_sustained(d == "PARAMETER_UPDATE_SUFFICIENT", fault) pre = slice(BURN_IN, fault) res[f"{tag}_pre_fault_alarm_rate"] = float(np.mean(d[pre] != "OK")) res["ph_detect"] = first_sustained(ph_alarms, fault) res["ph_pre_fault_alarm_rate"] = float(np.mean(ph_alarms[slice(BURN_IN, fault)])) else: post = slice(BURN_IN, T_TOTAL) res["geo01_false_alarm_rate"] = float(np.mean(dispositions[post] != "OK")) res["geo01b_false_alarm_rate"] = float(np.mean(dispositions_b[post] != "OK")) res["ph_false_alarm_rate"] = float(np.mean(ph_alarms[post])) # keep traces for plotting res["_traces"] = { "resid": np.abs(z - preds), "i": comps["i"], "mean_resid": comps["mean_resid"], "mean_grad": comps["mean_grad"], "dispositions": dispositions_b, "ph_alarms": ph_alarms, } return res def main(): results = [run_scenario(s) for s in ["A", "B", "C"]] print("=" * 72) print("GEO-01 vs Page-Hinkley — synthetic benchmark v0.1") print(f"T={T_TOTAL}, fault at t={FAULT_T}, burn-in={BURN_IN}, window={WIN}") print("=" * 72) for r in results: print(f"\nScenario {r['scenario']}" + (f" (fault @ {r['fault_t']})" if r["fault_t"] else " (control)")) if r["fault_t"]: def fmt(v): return f"t={v} (+{v - r['fault_t']})" if v else "never" for tag, label in [("geo01", "GEO-01 (grad-norm) "), ("geo01b", "GEO-01b (loss-slope)")]: print(f" {label} INSTR_FAIL : {fmt(r[f'{tag}_instr_detect'])}") print(f" {label} GEOM_EXP : {fmt(r[f'{tag}_geom_detect'])}") print(f" {label} PARAM_OK : {fmt(r[f'{tag}_tune_detect'])}") print(f" Page-Hinkley alarm : {fmt(r['ph_detect'])}") print(f" Pre-fault alarm rate GEO-01={r['geo01_pre_fault_alarm_rate']:.3f}" f" GEO-01b={r['geo01b_pre_fault_alarm_rate']:.3f}" f" PH={r['ph_pre_fault_alarm_rate']:.3f}") else: print(f" False alarm rate GEO-01={r['geo01_false_alarm_rate']:.3f}" f" GEO-01b={r['geo01b_false_alarm_rate']:.3f}" f" PH={r['ph_false_alarm_rate']:.3f}") # plot try: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True) for ax, r in zip(axes, results): tr = r["_traces"] ax.plot(tr["mean_resid"], lw=1, label="rolling |resid|", color="steelblue") ax.plot(tr["i"], lw=1, label="i_t (instr.)", color="darkorange") d = tr["dispositions"] for lbl, col in [("INSTRUMENTATION_FAILURE", "red"), ("GEOMETRY_EXPANSION_CANDIDATE", "purple")]: idx = np.where(d == lbl)[0] if len(idx): ax.scatter(idx, np.full(len(idx), -0.05), s=2, color=col, label=f"GEO-01: {lbl}") ph_idx = np.where(tr["ph_alarms"])[0] if len(ph_idx): ax.scatter(ph_idx, np.full(len(ph_idx), -0.12), s=2, color="green", label="PH alarm") if r["fault_t"]: ax.axvline(r["fault_t"], color="k", ls="--", lw=1) ax.set_title(f"Scenario {r['scenario']}") ax.legend(fontsize=7, loc="upper left", ncol=2) axes[-1].set_xlabel("t") plt.tight_layout() plt.savefig("/home/claude/geo01_results.png", dpi=110) print("\nPlot saved: geo01_results.png") except Exception as ex: print(f"\n(plotting skipped: {ex})") if __name__ == "__main__": main()