chore(sync): mirror backend .py + Dockerfile to Space (hf-sync-backend)

Automated backend sync from szl-holdings/a11oy main via hf-sync-backend.
Updated (differed from the Space): Dockerfile, serve.py, szl_governed_ipinn.py, szl_pinn_inverse.py
Deleted (gone from the repo + Dockerfile COPY set): (none)

Keeps the Space-built backend (serve.py + the Dockerfile-COPY'd .py
modules) identical to GitHub main so the Space never rebuilds from a
stale backend, new endpoints don't 404 there, and orphaned modules
removed from the repo don't linger in the Space tree.

Dockerfile

@@ -321,6 +321,10 @@ COPY a11oy_uds_portability_nav.py ./
 # module honestly serves a SAMPLE certificate until Forge writes real ones on the box.
 COPY szl_pinn_bounds.py ./
 COPY physical_bounds_certificate.json agentic_decision_trail.json physical_bounds_certificate.dsse.json ./
+# Governed Inverse-PINN engine (governed-inverse-pinn) — adds POST /api/a11oy/v1/pinn/identify
+# (+ GET demo, GET /pinn/health). Both modules MUST be COPY'd or serve.py's guarded import
+# falls back (merged-but-not-live) in the HF image. NumPy-only (no torch/DeepXDE/scipy added).
+COPY szl_pinn_inverse.py szl_governed_ipinn.py ./
 # PNT / quantum-sensing mesh (pure-stdlib closed-form web path; serves /api/a11oy/v1/pnt/*).
 # szl_pnt_mesh.py loads the 4 engine modules dynamically via importlib, so ALL FIVE MUST be
 # COPY'd or serve.py's guarded import falls back to a stub (merged-but-not-live) in the HF

serve.py

@@ -643,6 +643,28 @@ try:
 except Exception as _szl_pinn_e:  # pragma: no cover
     print(f"[a11oy] Agentic-PINN + physical-bounds mesh NOT registered: {_szl_pinn_e!r}", file=__import__("sys").stderr)
 
+# ── Governed Inverse-PINN engine (governed-inverse-pinn) — adds the INVERSE
+# discovery surface POST /api/a11oy/v1/pinn/identify (+ GET demo, GET /pinn/health,
+# alias prefix /v1/pinn). Discovers unknown PHYSICAL parameters of an ODE/PDE from
+# data with an HONEST self-doubt gate: a parameter the data cannot identify is
+# labelled RED/UNIDENTIFIABLE and the engine REFUSES to assert a value. NumPy-only
+# (no torch/DeepXDE/scipy): spectral surrogate with exact analytic derivatives,
+# ridge-LS data fit, exact LS for linear params, Adam GD on the physics residual,
+# FIM identifiability, bootstrap 95% CI, three-state GREEN/YELLOW/RED convergence.
+# Values are MODELED (a fit to data, never MEASURED); F19/Bekenstein is a
+# locked-proven inequality APPLIED (not re-claimed); Λ=Conjecture 1 (advisory ≤0.99);
+# DSSE receipt is honest-UNSIGNED until the on-metal cosign key signs it (never faked).
+# Additive, try/except-guarded, registered BEFORE the /api/a11oy/{path:path} Node
+# proxy + SPA catch-all (defined at the file tail) so it wins ordered matching. The
+# guard is HARD: any import/register failure logs and continues — a11oy boots even if
+# this engine is broken, and the engine NEVER raises into app startup.
+try:
+    import szl_governed_ipinn as _szl_governed_ipinn
+    _szl_ipinn_routes = _szl_governed_ipinn.register(app, ns="a11oy")
+    print(f"[a11oy] Governed Inverse-PINN registered: POST /api/a11oy/v1/pinn/identify (+ /pinn/health) {_szl_ipinn_routes}", file=__import__("sys").stderr)
+except Exception as _szl_ipinn_e:  # pragma: no cover
+    print(f"[a11oy] Governed Inverse-PINN NOT registered (a11oy continues): {_szl_ipinn_e!r}", file=__import__("sys").stderr)
+
 # ── Compliance crosswalk MESH (compliance-mesh) — closes the audited gap where the
 # doctrine-v11 → NIST AI RMF / ISO 42001 / EU AI Act crosswalk module existed
 # (szl_compliance_mesh.py + compliance_crosswalk.py, REAL honest data with

szl_governed_ipinn.py

@@ -0,0 +1,412 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: Apache-2.0
+# (c) 2026 Lutar, Stephen P. - SZL Holdings - ORCID 0009-0001-0110-4173
+#
+# szl_governed_ipinn.py — Governed wrapper + HTTP surface for the SZL Inverse-PINN
+# engine (szl_pinn_inverse). Taxonomy: services (frontier discovery) + provenance.
+#
+# Doctrine v11 LOCKED 749/14/163 @ c7c0ba17 · Lambda = Conjecture 1 (advisory).
+#
+# WHAT THIS ADDS over the bare engine
+#   1. governed_discover(spec): runs the engine and returns, PER discovered
+#      parameter: value, 95% CI, a three-state convergence label (GREEN/YELLOW/
+#      RED) with the EXACT numeric criteria, the physics residual, and a
+#      Bekenstein/F19 information-cost ratio (PHYSICALLY_PLAUSIBLE / IMPLAUSIBLE).
+#      A parameter the data cannot identify is RED/UNIDENTIFIABLE and the engine
+#      REFUSES to assert a value for it (value withheld, reason given).
+#   2. A Lambda advisory (Conjecture 1) in [0, 0.99] — NEVER 1.0.
+#   3. A DSSE-signable receipt dict (organ="a11oy-pinn"); signed with the real
+#      cosign key when present, otherwise an HONEST UNSIGNED envelope (never a
+#      fabricated signature). The ledger write itself is NOT done here — we call
+#      record_pinn_receipt(receipt), a no-op-safe hook that Dev C wires to
+#      szl_lake_ingest.record_receipt.
+#   4. POST /api/a11oy/v1/pinn/identify — the endpoint. A built-in demo
+#      ("demo":"duffing") returns a real GREEN alpha ~ 1.0 out of the box.
+#
+# HONEST LABELS: all discovered values are MODELED (a fit to data), never
+#   MEASURED. F19 (Bekenstein bound) is one of the 8 locked-proven inequalities
+#   {F1,F4,F7,F11,F12,F18,F19,F22}@c7c0ba17 — a PROVEN inequality, never an
+#   assertion; its APPLICATION to an information-cost ratio here is MODELED.
+#   The locked-proven count is 8. Lambda = Conjecture 1. No user-visible codenames.
+
+import math
+import time
+import json
+
+import numpy as np
+
+from szl_pinn_inverse import (
+    SZLInversePINN, SZLInversePINNTrainer, SZLSpectralSurrogate,
+    duffing_residual, integrate_duffing,
+    CAUSAL_GREEN, CAUSAL_RED, GRAD_GREEN, KAPPA_IDENT, KAPPA_RED, FISHER_FLOOR,
+    MIN_DATA_POINTS,
+)
+
+RECEIPT_SCHEMA = "szl.lake.receipt/v1"
+RECEIPT_ORGAN = "a11oy-pinn"
+RECEIPT_PAYLOAD_TYPE = "application/vnd.szl.ipinn+json"
+LOCKED_PROVEN = ("F1", "F4", "F7", "F11", "F12", "F18", "F19", "F22")
+LOCKED_PROVEN_AT = "c7c0ba17"
+
+# Built-in, code-safe systems. We NEVER eval user-supplied residual source — a
+# request selects a named built-in physics, or passes data to identify against
+# one. (Security rule: no arbitrary code-as-action on this path.)
+_BUILTIN_SYSTEMS = {
+    "duffing": {
+        "residual": duffing_residual,
+        "unknowns": ["alpha"],
+        "linear": ["alpha"],
+        "bounds": {"alpha": (-10.0, 10.0)},
+        "inits": {"alpha": 0.4},
+        "truth": {"alpha": 1.0},
+        "desc": "Duffing oscillator m x'' + c x' + delta x + alpha x^3 = F cos(omega t)",
+    },
+}
+
+
+# ---------------------------------------------------------------------------
+# Bekenstein / F19 information-cost ratio.
+#   I_eta   = log2(sigma_prior / sigma_posterior)         bits of info gained
+#   I_max   = 2*pi*R*E / (hbar*c) / ln(2)                 Bekenstein bound (bits)
+#   ratio   = I_eta / I_max                               > 1 => IMPLAUSIBLE
+# F19 is the PROVEN Bekenstein inequality (locked-8). The numbers below are a
+# MODELED application with SAMPLE R, E unless the caller supplies real ones.
+# ---------------------------------------------------------------------------
+_HBAR = 1.054571817e-34   # J*s
+_C = 2.99792458e8         # m/s
+
+
+def bekenstein_ratio(sigma_prior, sigma_posterior, radius_m=1.0, energy_j=1.0):
+    sp = max(float(sigma_prior), 1e-30)
+    sq = max(float(sigma_posterior), 1e-30)
+    info_bits = math.log2(sp / sq) if sp > sq else 0.0
+    i_max = (2.0 * math.pi * float(radius_m) * float(energy_j)) / (_HBAR * _C) / math.log(2.0)
+    ratio = info_bits / i_max if i_max > 0 else float("inf")
+    return {
+        "info_bits": info_bits,
+        "bekenstein_max_bits": i_max,
+        "ratio": ratio,
+        "label": "PHYSICALLY_PLAUSIBLE" if ratio <= 1.0 else "PHYSICALLY_IMPLAUSIBLE",
+        "radius_m": float(radius_m),
+        "energy_j": float(energy_j),
+        "basis": ("F19 Bekenstein bound = PROVEN inequality (locked-8 @ %s); "
+                  "this application is MODELED with SAMPLE R,E unless supplied"
+                  % LOCKED_PROVEN_AT),
+    }
+
+
+# ---------------------------------------------------------------------------
+# Lambda advisory (Conjecture 1) — weighted geometric mean of honest factors,
+# HARD-capped at 0.99. Never 1.0, never presented as proven.
+# ---------------------------------------------------------------------------
+def compute_lambda(label, frac_asserted, data_rms, delta_param_rel):
+    f_label = {"GREEN": 0.9, "YELLOW": 0.6, "RED": 0.2}.get(label, 0.2)
+    f_ident = 0.05 + 0.95 * float(np.clip(frac_asserted, 0.0, 1.0))
+    f_data = float(np.clip(math.exp(-5.0 * max(data_rms, 0.0)), 0.05, 1.0))
+    f_stab = float(np.clip(1.0 / (1.0 + max(delta_param_rel, 0.0)), 0.05, 1.0))
+    geom = (f_label * f_ident * f_data * f_stab) ** 0.25
+    return {
+        "value": round(min(geom, 0.99), 4),
+        "status": "ADVISORY",
+        "basis": "Lambda = Conjecture 1 (advisory, capped <= 0.99; NEVER a proof)",
+        "factors": {"label": f_label, "identifiable": round(f_ident, 4),
+                    "data_fit": round(f_data, 4), "stability": round(f_stab, 4)},
+    }
+
+
+# ---------------------------------------------------------------------------
+# JSON-safe coercion. Starlette's JSONResponse uses allow_nan=False, so a
+# non-finite float (e.g. kappa(FIM)=inf — the honest non-identifiable signal)
+# would 500 the response. We convert non-finite floats to honest string tokens
+# and numpy scalars to plain Python, so the receipt and the wire agree.
+# ---------------------------------------------------------------------------
+def _json_safe(obj):
+    if isinstance(obj, dict):
+        return {k: _json_safe(v) for k, v in obj.items()}
+    if isinstance(obj, (list, tuple)):
+        return [_json_safe(v) for v in obj]
+    if isinstance(obj, np.generic):
+        obj = obj.item()
+    if isinstance(obj, float):
+        if math.isinf(obj):
+            return "Infinity" if obj > 0 else "-Infinity"
+        if math.isnan(obj):
+            return "NaN"
+    return obj
+
+
+# ---------------------------------------------------------------------------
+# The no-op-safe ledger hook. Dev C wires this to szl_lake_ingest.record_receipt;
+# until then (and when running standalone) it degrades honestly to not-recorded.
+# ---------------------------------------------------------------------------
+def record_pinn_receipt(receipt):
+    """Hook for Dev C. Signature: record_pinn_receipt(receipt: dict) -> dict.
+    Attempts an in-process ledger append via szl_lake_ingest.record_receipt with
+    organ="a11oy-pinn"; never raises, returns a status dict."""
+    try:
+        import szl_lake_ingest  # type: ignore
+        res = szl_lake_ingest.record_receipt(receipt, organ=RECEIPT_ORGAN)
+        return {"recorded": True, "backend": "szl_lake_ingest.record_receipt",
+                "result": res if isinstance(res, dict) else str(res)}
+    except Exception as e:  # noqa: BLE001 — honest degrade, never fatal
+        return {"recorded": False, "reason": "ledger hook not wired (%r)" % e,
+                "note": "Dev C wires record_pinn_receipt -> szl_lake_ingest"}
+
+
+def build_ipinn_receipt(system, method, params_block, convergence, lambda_adv, sign=True):
+    payload = {
+        "schema": RECEIPT_SCHEMA,
+        "organ": RECEIPT_ORGAN,
+        "kind": "inverse_pinn_identify",
+        "ts": time.time(),
+        "system": system,
+        "method": method,
+        "label_provenance": "MODELED (fit to data; not MEASURED)",
+        "discovered": params_block,
+        "convergence": convergence,
+        "lambda_advisory": lambda_adv,
+        "doctrine": {
+            "locked_proven_count": 8,
+            "locked_proven": list(LOCKED_PROVEN),
+            "locked_at": LOCKED_PROVEN_AT,
+            "lambda": "Conjecture 1",
+            "f19": "Bekenstein bound = PROVEN inequality (locked-8); application MODELED",
+        },
+    }
+    receipt = {"payload": payload}
+    if sign:
+        try:
+            import szl_dsse  # type: ignore
+            env = szl_dsse.sign_payload(payload, RECEIPT_PAYLOAD_TYPE)
+            receipt["dsse"] = env
+            receipt["signed"] = bool(env.get("signatures"))
+        except Exception as e:  # noqa: BLE001
+            receipt["dsse"] = {"signed": False,
+                               "reason": "szl_dsse unavailable (%r)" % e}
+            receipt["signed"] = False
+    else:
+        receipt["signed"] = False
+    return receipt
+
+
+# ---------------------------------------------------------------------------
+# The governed discovery orchestrator.
+# ---------------------------------------------------------------------------
+def _coerce_data(spec):
+    """Return (t, y, noise_sigma, system_key, used_demo). Either a built-in demo
+    ('demo':'duffing') generating synthetic data, or caller-supplied
+    {'data': {'t': [...], 'x': [...]}, 'system': 'duffing'}."""
+    demo = spec.get("demo")
+    system_key = (spec.get("system") or demo or "duffing")
+    if isinstance(system_key, str):
+        system_key = system_key.strip().lower()
+    if system_key not in _BUILTIN_SYSTEMS:
+        raise ValueError("unsupported system %r; supported: %s (or pass demo='duffing')"
+                         % (system_key, list(_BUILTIN_SYSTEMS)))
+    data = spec.get("data")
+    if demo or not data:
+        opts = spec.get("options") or {}
+        n = int(opts.get("n_points", 160))
+        n = max(MIN_DATA_POINTS, min(n, 600))
+        t1 = float(opts.get("t_max", 10.0))
+        t = np.linspace(0.0, t1, n)
+        truth = _BUILTIN_SYSTEMS[system_key]["truth"]
+        x = integrate_duffing(t, alpha=truth.get("alpha", 1.0))
+        noise = float(opts.get("noise", 0.0))
+        if noise > 0:
+            x = x + np.random.default_rng(0).normal(0.0, noise, size=x.shape)
+        return t, x, max(noise, 1e-3), system_key, True
+    # caller-supplied data
+    t = np.asarray(data["t"], float).reshape(-1)
+    y = np.asarray(data.get("x", data.get("y")), float).reshape(-1)
+    if t.shape[0] != y.shape[0]:
+        raise ValueError("data.t and data.x must have equal length")
+    return t, y, 1e-3, system_key, False
+
+
+def governed_discover(spec):
+    spec = dict(spec or {})
+    try:
+        t, y, noise_sigma, system_key, used_demo = _coerce_data(spec)
+    except (ValueError, KeyError, TypeError) as ce:
+        return {"ok": False, "error": str(ce), "honesty": _honesty()}
+    if t.shape[0] < MIN_DATA_POINTS:
+        return {"ok": False,
+                "error": "need >= %d data points (got %d)" % (MIN_DATA_POINTS, t.shape[0]),
+                "honesty": _honesty()}
+    sysdef = _BUILTIN_SYSTEMS[system_key]
+    requested = spec.get("unknowns") or list(sysdef["unknowns"])
+    # the residual only knows about its own params; unknowns not in the residual
+    # (e.g. "ghost") are admitted on purpose so the self-doubt gate can REFUSE them.
+    inits = dict(sysdef["inits"])
+    bounds = dict(sysdef["bounds"])
+    linear = list(sysdef["linear"])
+    for nm in requested:
+        if nm not in inits:
+            inits[nm] = 0.4
+            bounds.setdefault(nm, (-10.0, 10.0))
+            linear.append(nm)   # treat unknowns as linear unless engine proves otherwise
+
+    opts = spec.get("options") or {}
+    n_modes = int(opts.get("n_modes", 24))
+    epochs = int(opts.get("epochs", 50))
+    restarts = int(opts.get("restarts", 6))
+
+    surrogate = SZLSpectralSurrogate(n_modes=n_modes, poly_deg=3)
+    model = SZLInversePINN(sysdef["residual"], {k: inits[k] for k in requested},
+                           surrogate=surrogate, param_bounds=bounds,
+                           linear_params=[k for k in linear if k in requested])
+    trainer = SZLInversePINNTrainer(model, t, y, t, noise_sigma=noise_sigma, seed=1)
+    try:
+        record = trainer.fit(epochs=epochs)
+    except ValueError as ve:
+        return {"ok": False, "error": str(ve), "honesty": _honesty()}
+    results = trainer.param_results(n_restarts=restarts, epochs=max(20, epochs // 2))
+
+    # prior std per param (for the Bekenstein info-gain): uniform-prior std over bounds.
+    discovered = []
+    n_asserted = 0
+    for pr in results:
+        lo, hi = bounds.get(pr.name, (-10.0, 10.0))
+        sigma_prior = (hi - lo) / math.sqrt(12.0) if math.isfinite(hi - lo) else 10.0
+        sigma_post = pr.std if pr.std > 0 else max(abs(pr.ci_high - pr.ci_low) / 3.92, 1e-6)
+        bek = bekenstein_ratio(sigma_prior, sigma_post,
+                               radius_m=float(opts.get("radius_m", 1.0)),
+                               energy_j=float(opts.get("energy_j", 1.0)))
+        block = {
+            "name": pr.name,
+            "asserted": pr.asserted,
+            "value": (round(pr.value, 6) if pr.asserted else None),
+            "ci95": ([round(pr.ci_low, 6), round(pr.ci_high, 6)] if pr.asserted else None),
+            "std": round(pr.std, 6),
+            "fisher_information": pr.fisher,
+            "identifiable": pr.identifiable,
+            "convergence_label": record.label if pr.asserted else "RED",
+            "bekenstein": bek,
+            "label": "MODELED",
+        }
+        if not pr.asserted:
+            if pr.fisher < FISHER_FLOOR:
+                why = ("Fisher information %.2e is below the floor %.0e — the data carry "
+                       "no information about this parameter" % (pr.fisher, FISHER_FLOOR))
+            else:
+                why = ("the FIM is ill-conditioned (kappa=%.2e >= %.0e) — the parameters "
+                       "are jointly non-identifiable" % (record.kappa_fim, KAPPA_RED))
+            block["refusal"] = ("UNIDENTIFIABLE: %s. The engine REFUSES to assert this "
+                                "parameter." % why)
+        else:
+            n_asserted += 1
+        discovered.append(block)
+
+    frac_asserted = n_asserted / max(1, len(results))
+    convergence = {
+        "label": record.label,
+        "criteria": record.criteria,
+        "min_causal_weight": round(record.min_causal_weight, 6),
+        "grad_norm": record.grad_norm,
+        "kappa_fim": record.kappa_fim,
+        "residual_rms": round(record.residual_rms, 6),
+        "data_rms": round(record.data_rms, 6),
+        "epochs_run": record.epochs_run,
+        "thresholds": {
+            "causal_green": CAUSAL_GREEN, "causal_red": CAUSAL_RED,
+            "grad_green": GRAD_GREEN, "kappa_ident": KAPPA_IDENT, "kappa_red": KAPPA_RED,
+        },
+    }
+    convergence = _json_safe(convergence)
+    discovered = _json_safe(discovered)
+    lambda_adv = compute_lambda(record.label, frac_asserted,
+                                record.data_rms, record.delta_param_rel)
+    method = {
+        "engine": "szl_pinn_inverse.SZLInversePINN",
+        "surrogate": "spectral basis (Fourier %d modes + poly deg 3), exact analytic "
+                     "derivatives; NumPy-only (no torch/DeepXDE)" % n_modes,
+        "param_solve": "exact least-squares for linear params; Adam GD on physics "
+                       "residual for nonlinear; FIM identifiability self-doubt gate",
+        "data_label": "MODELED synthetic (demo)" if used_demo else "caller-supplied",
+        "system": sysdef["desc"],
+    }
+    receipt = build_ipinn_receipt({"key": system_key, "desc": sysdef["desc"]},
+                                  method, discovered, convergence, lambda_adv,
+                                  sign=bool(opts.get("sign", True)))
+    ledger = record_pinn_receipt(receipt)
+
+    return {
+        "ok": True,
+        "system": system_key,
+        "convergence": convergence,
+        "discovered": discovered,
+        "lambda_advisory": lambda_adv,
+        "receipt": receipt,
+        "ledger": ledger,
+        "method": method,
+        "honesty": _honesty(),
+    }
+
+
+def _honesty():
+    return {
+        "values": "MODELED (fit to data; not MEASURED)",
+        "locked_proven_count": 8,
+        "locked_proven": list(LOCKED_PROVEN),
+        "lambda": "Conjecture 1 (advisory, <= 0.99)",
+        "f19": "Bekenstein bound = PROVEN inequality (locked-8); application MODELED",
+        "self_doubt_gate": ("a non-identifiable parameter (Fisher < %.0e or kappa(FIM) "
+                            ">= %.0e) is labelled RED/UNIDENTIFIABLE and NOT asserted"
+                            % (FISHER_FLOOR, KAPPA_RED)),
+    }
+
+
+# ---------------------------------------------------------------------------
+# HTTP surface — POST /api/a11oy/v1/pinn/identify (+ GET health/info).
+# Registered BEFORE the SPA catch-all (front-inserted by serve.py).
+# ---------------------------------------------------------------------------
+def register(app, ns="a11oy"):
+    from fastapi.responses import JSONResponse
+    from fastapi import Request
+
+    async def _identify(request: Request):
+        try:
+            try:
+                spec = await request.json()
+            except Exception:
+                spec = {}
+            if not isinstance(spec, dict):
+                spec = {}
+            if not spec:
+                spec = {"demo": "duffing"}
+            out = governed_discover(spec)
+            code = 200 if out.get("ok") else 400
+            label = out.get("convergence", {}).get("label", "NA")
+            return JSONResponse(out, status_code=code,
+                                headers={"x-szl-pinn-label": str(label),
+                                         "x-szl-organ": RECEIPT_ORGAN})
+        except Exception as e:  # noqa: BLE001
+            return JSONResponse({"ok": False, "error": "%r" % e, "honesty": _honesty()},
+                                status_code=500)
+
+    async def _health():
+        return JSONResponse({
+            "ok": True, "organ": RECEIPT_ORGAN,
+            "endpoint": "POST /api/%s/v1/pinn/identify" % ns,
+            "supported_systems": list(_BUILTIN_SYSTEMS),
+            "demo": "POST {\"demo\":\"duffing\"} -> GREEN alpha ~ 1.0",
+            "honesty": _honesty(),
+        })
+
+    prefixes = ["/api/%s/v1/pinn" % ns, "/v1/pinn"]
+    routes = []
+    for p in prefixes:
+        app.add_api_route("%s/identify" % p, _identify, methods=["POST", "GET"],
+                          include_in_schema=True)
+        app.add_api_route("%s/health" % p, _health, methods=["GET"],
+                          include_in_schema=True)
+        routes += ["%s/identify" % p, "%s/health" % p]
+    return routes
+
+
+if __name__ == "__main__":
+    out = governed_discover({"demo": "duffing"})
+    print(json.dumps({k: out[k] for k in ("ok", "system", "convergence", "discovered",
+                                          "lambda_advisory", "ledger")},
+                     indent=2, default=float)[:2000])

szl_pinn_inverse.py

@@ -0,0 +1,578 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: Apache-2.0
+# (c) 2026 Lutar, Stephen P. - SZL Holdings - ORCID 0009-0001-0110-4173
+#
+# szl_pinn_inverse.py — SZL Governed Inverse-PINN engine (a11oy frontier).
+# Taxonomy home: services (frontier discovery surface) + provenance (signed receipt).
+#
+# Doctrine v11 LOCKED 749/14/163 @ c7c0ba17 · Lambda = Conjecture 1 (advisory only).
+#
+# WHAT THIS IS
+#   A self-contained INVERSE physics-informed solver that discovers unknown
+#   PHYSICAL PARAMETERS of an ODE/PDE from data, with an HONEST self-doubt gate:
+#   a parameter the data cannot identify is labelled RED / UNIDENTIFIABLE and the
+#   engine REFUSES to assert a value for it (Fisher-information gate, below).
+#
+# OWN CODE / PERMISSIVE DEPS ONLY (NumPy BSD-3)
+#   No torch, no DeepXDE (LGPL), nothing proprietary. Every equation is
+#   RE-IMPLEMENTED from the public literature (per-equation citation map in
+#   team/frontier/PINN_BACKEND.md), not copied from any GPL/LGPL package.
+#
+#   The surrogate that represents x(t) is a LINEAR spectral basis (Fourier modes
+#   + a low-order polynomial trend). A linear basis is the pragmatic, robust
+#   choice for an autograd-free NumPy build: its first/second time-derivatives
+#   are EXACT and analytic (no fragile finite differencing, no second-order
+#   backprop), the data fit is a single regularised least-squares solve (fast,
+#   CPU-only, seconds), and the physics residual is then well-conditioned. A
+#   tanh-MLP surrogate (SZLPinnNet) with exact analytic input-derivatives is also
+#   provided for callers who prefer it, but the governed endpoint defaults to the
+#   spectral basis because it is the one that converges reliably on a cpu-basic
+#   Space without a heavy autodiff dependency.
+#
+#   Parameters that enter the residual LINEARLY (e.g. Duffing alpha) are solved by
+#   exact least squares; any others are refined by gradient descent on the
+#   physics residual. Identifiability is then checked via the Fisher Information
+#   Matrix BEFORE any value is asserted.
+#
+# HONEST LABELS: every numeric result here is MODELED (a fit to data), never
+#   MEASURED. The convergence label is GREEN / YELLOW / RED with the EXACT numeric
+#   criteria from team/frontier/ARXIV_LEADERS.md — see _classify_convergence().
+#   The half-state ("looks done but isn't") is unacceptable.
+
+import math
+from dataclasses import dataclass, field
+from typing import Callable, Dict, List, Optional, Sequence, Tuple
+
+import numpy as np
+
+# ---------------------------------------------------------------------------
+# EXACT convergence thresholds — from ARXIV_LEADERS.md (three-state GREEN/YELLOW/
+# RED rules). These are the contract the front-end + receipt rely on; do NOT
+# loosen them without updating the spec.
+# ---------------------------------------------------------------------------
+CAUSAL_GREEN = 0.99          # min(w_causal) > 0.99  -> temporally converged
+CAUSAL_RED = 0.50            # min(w_causal) <= 0.50 -> diverged
+GRAD_GREEN = 1e-5            # ||grad L_A (params)|| < 1e-5 -> stationary
+KAPPA_IDENT = 1e6            # kappa(FIM) < 1e6  -> identifiable
+KAPPA_RED = 1e8              # kappa(FIM) >= 1e8 -> non-identifiable (RED)
+FISHER_FLOOR = 1e-8          # per-param Fisher information floor (self-doubt gate)
+EPSILON_CAUSAL = 0.01        # epsilon in w_causal = exp(-eps * cumsum r^2)
+MIN_DATA_POINTS = 10         # below this the engine refuses to assert anything
+
+__all__ = [
+    "SZLSpectralSurrogate",
+    "SZLPinnNet",
+    "SZLInversePINN",
+    "SZLInversePINNTrainer",
+    "ConvergenceRecord",
+    "ParamResult",
+    "duffing_residual",
+    "integrate_duffing",
+]
+
+
+# ===========================================================================
+# 1a. Spectral-basis surrogate (DEFAULT) — exact analytic time-derivatives.
+# ===========================================================================
+class SZLSpectralSurrogate:
+    """x(t) ~ sum_p a_p t^p  +  sum_k [ b_k sin(w_k tau) + c_k cos(w_k tau) ],
+    tau = t - t0, w_k = 2*pi*k / span. A LINEAR-in-coefficients model: the
+    design matrices for x, dx/dt, d2x/dt2 are exact and analytic, so the physics
+    residual needs no finite differencing of the surrogate. Re-implemented from
+    standard Fourier/Chebyshev collocation theory (own code)."""
+
+    def __init__(self, n_modes: int = 24, poly_deg: int = 3, ridge: float = 1e-6):
+        self.K = int(n_modes)
+        self.D = int(poly_deg)
+        self.ridge = float(ridge)
+        self.coef: Optional[np.ndarray] = None
+        self.t0 = 0.0
+        self.span = 1.0
+
+    def n_features(self) -> int:
+        return (self.D + 1) + 2 * self.K
+
+    def design(self, t: np.ndarray, order: int = 0) -> np.ndarray:
+        t = np.asarray(t, float).reshape(-1)
+        tau = t - self.t0
+        cols: List[np.ndarray] = []
+        for p in range(self.D + 1):
+            if order == 0:
+                cols.append(t ** p)
+            elif order == 1:
+                cols.append(p * t ** (p - 1) if p >= 1 else np.zeros_like(t))
+            else:
+                cols.append(p * (p - 1) * t ** (p - 2) if p >= 2 else np.zeros_like(t))
+        for k in range(1, self.K + 1):
+            w = 2.0 * math.pi * k / self.span
+            if order == 0:
+                cols += [np.sin(w * tau), np.cos(w * tau)]
+            elif order == 1:
+                cols += [w * np.cos(w * tau), -w * np.sin(w * tau)]
+            else:
+                cols += [-w * w * np.sin(w * tau), -w * w * np.cos(w * tau)]
+        return np.stack(cols, axis=1)
+
+    def fit(self, t: np.ndarray, y: np.ndarray):
+        t = np.asarray(t, float).reshape(-1)
+        y = np.asarray(y, float).reshape(-1)
+        self.t0 = float(t.min())
+        self.span = float(t.max() - t.min()) or 1.0
+        A = self.design(t, 0)
+        G = A.T @ A + self.ridge * np.eye(A.shape[1])
+        self.coef = np.linalg.solve(G, A.T @ y)
+        return self
+
+    def derivatives(self, t: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        if self.coef is None:
+            raise RuntimeError("surrogate not fitted")
+        x = self.design(t, 0) @ self.coef
+        dx = self.design(t, 1) @ self.coef
+        ddx = self.design(t, 2) @ self.coef
+        return x, dx, ddx
+
+    def predict(self, t: np.ndarray) -> np.ndarray:
+        return self.design(t, 0) @ self.coef
+
+
+# ===========================================================================
+# 1b. Optional tanh-MLP surrogate — exact analytic input-derivatives.
+#     Provided for spec-completeness/generality; the governed endpoint defaults
+#     to SZLSpectralSurrogate (more robust autograd-free convergence on CPU).
+# ===========================================================================
+class SZLPinnNet:
+    """1-D -> 1-D tanh MLP with EXACT first/second analytic derivatives of the
+    output w.r.t. the scalar input, plus manual reverse-mode backprop of the data
+    MSE. Re-implemented from first principles (chain rule for tanh layers;
+    Rumelhart et al. backprop). NumPy only."""
+
+    def __init__(self, layers: Sequence[int] = (1, 32, 32, 1), seed: int = 0):
+        self.layers = list(layers)
+        rng = np.random.default_rng(seed)
+        self.W: List[np.ndarray] = []
+        self.b: List[np.ndarray] = []
+        for nin, nout in zip(self.layers[:-1], self.layers[1:]):
+            self.W.append(rng.normal(0.0, math.sqrt(1.0 / nin), size=(nin, nout)))
+            self.b.append(np.zeros((1, nout)))
+        self.in_mean = 0.0
+        self.in_scale = 1.0
+        self.out_mean = 0.0
+        self.out_scale = 1.0
+        self._cache: dict = {}
+
+    def set_norm(self, t, y):
+        self.in_mean = float(np.mean(t)); self.in_scale = float(np.std(t)) or 1.0
+        self.out_mean = float(np.mean(y)); self.out_scale = float(np.std(y)) or 1.0
+
+    def _raw_forward(self, tn: np.ndarray) -> np.ndarray:
+        a = tn.reshape(-1, 1)
+        zs, acts = [], [a]
+        for i, (W, b) in enumerate(zip(self.W, self.b)):
+            z = a @ W + b
+            zs.append(z)
+            a = np.tanh(z) if i < len(self.W) - 1 else z
+            acts.append(a)
+        self._cache = {"zs": zs, "acts": acts}
+        return a.reshape(-1)
+
+    def predict(self, t: np.ndarray) -> np.ndarray:
+        tn = (np.asarray(t, float) - self.in_mean) / self.in_scale
+        return self.out_mean + self.out_scale * self._raw_forward(tn)
+
+    def derivatives(self, t: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        tn = (np.asarray(t, float) - self.in_mean) / self.in_scale
+        a = tn.reshape(-1, 1)
+        da = np.ones_like(a)
+        dda = np.zeros_like(a)
+        for i, (W, b) in enumerate(zip(self.W, self.b)):
+            z = a @ W + b
+            dz = da @ W
+            ddz = dda @ W
+            if i < len(self.W) - 1:
+                th = np.tanh(z)
+                tp = 1.0 - th * th
+                tpp = -2.0 * th * tp
+                a = th
+                dda = tpp * dz * dz + tp * ddz
+                da = tp * dz
+            else:
+                a = z; da = dz; dda = ddz
+        s = self.out_scale
+        x = self.out_mean + s * a.reshape(-1)
+        dx = s * da.reshape(-1) / self.in_scale
+        ddx = s * dda.reshape(-1) / (self.in_scale ** 2)
+        return x, dx, ddx
+
+    def data_grads(self, t, y):
+        tn = (np.asarray(t, float) - self.in_mean) / self.in_scale
+        yn = (np.asarray(y, float) - self.out_mean) / self.out_scale
+        x = self._raw_forward(tn)
+        n = x.shape[0]
+        resid = (x - yn)
+        loss = float(0.5 * np.mean(resid * resid))
+        g = (resid / n).reshape(-1, 1)
+        acts = self._cache["acts"]; zs = self._cache["zs"]
+        gW = [None] * len(self.W); gb = [None] * len(self.b)
+        delta = g
+        for i in reversed(range(len(self.W))):
+            gW[i] = acts[i].T @ delta
+            gb[i] = delta.sum(axis=0, keepdims=True)
+            if i > 0:
+                delta = (delta @ self.W[i].T) * (1.0 - np.tanh(zs[i - 1]) ** 2)
+        return gW, gb, loss
+
+    def fit(self, t, y, epochs: int = 1500, lr: float = 5e-3):
+        self.set_norm(t, y)
+        mW = [np.zeros_like(w) for w in self.W]; vW = [np.zeros_like(w) for w in self.W]
+        mb = [np.zeros_like(b) for b in self.b]; vb = [np.zeros_like(b) for b in self.b]
+        b1, b2, e = 0.9, 0.999, 1e-8
+        for it in range(int(epochs)):
+            gW, gb, _ = self.data_grads(t, y)
+            i1 = it + 1
+            for j in range(len(self.W)):
+                mW[j] = b1 * mW[j] + (1 - b1) * gW[j]; vW[j] = b2 * vW[j] + (1 - b2) * gW[j] ** 2
+                self.W[j] -= lr * (mW[j] / (1 - b1 ** i1)) / (np.sqrt(vW[j] / (1 - b2 ** i1)) + e)
+                mb[j] = b1 * mb[j] + (1 - b1) * gb[j]; vb[j] = b2 * vb[j] + (1 - b2) * gb[j] ** 2
+                self.b[j] -= lr * (mb[j] / (1 - b1 ** i1)) / (np.sqrt(vb[j] / (1 - b2 ** i1)) + e)
+        return self
+
+
+# ===========================================================================
+# 2. The inverse model — surrogate + learnable physical parameters + residual.
+# ===========================================================================
+class SZLInversePINN:
+    """Bundles a surrogate (spectral by default), a dict of learnable physical
+    parameters, and a user-supplied residual callable r = f(t,x,dx,ddx,params)."""
+
+    def __init__(self, residual_fn: Callable[..., np.ndarray],
+                 param_inits: Dict[str, float],
+                 surrogate=None,
+                 param_bounds: Optional[Dict[str, Tuple[float, float]]] = None,
+                 linear_params: Optional[Sequence[str]] = None):
+        self.surrogate = surrogate if surrogate is not None else SZLSpectralSurrogate()
+        self.residual_fn = residual_fn
+        self.params: Dict[str, float] = dict(param_inits)
+        self.param_bounds = dict(param_bounds or {})
+        self.linear_params = list(linear_params or [])
+
+    def param_values(self) -> Dict[str, float]:
+        return dict(self.params)
+
+    def derivatives(self, t: np.ndarray):
+        return self.surrogate.derivatives(t)
+
+    def residual(self, t: np.ndarray) -> np.ndarray:
+        x, dx, ddx = self.surrogate.derivatives(t)
+        return np.asarray(self.residual_fn(t, x, dx, ddx, self.params), dtype=float)
+
+
+# ===========================================================================
+# 3. Records.
+# ===========================================================================
+@dataclass
+class ParamResult:
+    name: str
+    value: float
+    ci_low: float
+    ci_high: float
+    std: float
+    fisher: float
+    identifiable: bool
+    asserted: bool          # False -> engine REFUSES (RED / UNIDENTIFIABLE)
+
+
+@dataclass
+class ConvergenceRecord:
+    label: str                      # GREEN | YELLOW | RED
+    min_causal_weight: float
+    grad_norm: float
+    kappa_fim: float
+    delta_param_rel: float
+    residual_rms: float
+    data_rms: float
+    epochs_run: int
+    criteria: Dict[str, str] = field(default_factory=dict)
+    note: str = ""
+
+
+# ===========================================================================
+# 4. The trainer — surrogate fit + param solve (LS for linear, GD for nonlinear)
+#    + SA-PINN weights + causal weights + FIM self-doubt gate + ensemble CI.
+# ===========================================================================
+class SZLInversePINNTrainer:
+    def __init__(self, model: SZLInversePINN,
+                 t_data: np.ndarray, y_data: np.ndarray,
+                 t_colloc: Optional[np.ndarray] = None,
+                 w_phys: float = 1.0,
+                 lr_param: float = 5e-2,
+                 noise_sigma: Optional[float] = None,
+                 seed: int = 0):
+        self.m = model
+        self.t_data = np.asarray(t_data, float).reshape(-1)
+        self.y_data = np.asarray(y_data, float).reshape(-1)
+        self.t_colloc = (np.asarray(t_colloc, float).reshape(-1)
+                         if t_colloc is not None else self.t_data.copy())
+        self.w_phys = float(w_phys)
+        self.lr_param = float(lr_param)
+        self.noise_sigma = noise_sigma
+        self.rng = np.random.default_rng(seed)
+        self.sa = np.ones_like(self.t_colloc)       # SA-PINN self-adaptive weights
+        self.param_history: List[Dict[str, float]] = []
+
+    # ---- residual on the collocation grid ----
+    def _phys_residual(self) -> np.ndarray:
+        return self.m.residual(self.t_colloc)
+
+    # ---- exact least squares for params that enter r LINEARLY ----
+    def _ls_solve_linear(self):
+        names = list(self.m.linear_params)
+        if not names:
+            return
+        x, dx, ddx = self.m.derivatives(self.t_colloc)
+        base = dict(self.m.params)
+        G = np.zeros((self.t_colloc.shape[0], len(names)))
+        for j, nm in enumerate(names):
+            p = dict(base); p[nm] = base[nm] + 1.0
+            r1 = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p), float)
+            p[nm] = base[nm]
+            r0 = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p), float)
+            G[:, j] = r1 - r0
+        p0 = dict(base)
+        for nm in names:
+            p0[nm] = 0.0
+        r_at_zero = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p0), float)
+        w = np.sqrt(self.sa)
+        try:
+            sol, *_ = np.linalg.lstsq(G * w[:, None], -r_at_zero * w, rcond=None)
+            for j, nm in enumerate(names):
+                v = float(sol[j])
+                lo, hi = self.m.param_bounds.get(nm, (-np.inf, np.inf))
+                self.m.params[nm] = float(min(max(v, lo), hi))
+        except np.linalg.LinAlgError:
+            pass
+
+    # ---- gradient descent for NONLINEAR params (Adam on physics residual) ----
+    def _gd_nonlinear(self, epochs: int):
+        names = [k for k in self.m.params if k not in self.m.linear_params]
+        if not names:
+            return
+        mom = {k: 0.0 for k in names}; vel = {k: 0.0 for k in names}
+        b1, b2, e = 0.9, 0.999, 1e-8
+        for it in range(int(epochs)):
+            r = self._phys_residual()
+            self.sa = np.clip(self.sa + 0.01 * np.abs(r), 1.0, 50.0)
+            x, dx, ddx = self.m.derivatives(self.t_colloc)
+            i1 = it + 1
+            for nm in names:
+                d = max(1e-6, 1e-4 * (abs(self.m.params[nm]) + 1.0))
+                p = dict(self.m.params); p[nm] += d
+                rp = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p), float)
+                grad = float(np.mean(self.sa * r * (rp - r) / d)) * self.w_phys
+                mom[nm] = b1 * mom[nm] + (1 - b1) * grad
+                vel[nm] = b2 * vel[nm] + (1 - b2) * grad * grad
+                step = self.lr_param * (mom[nm] / (1 - b1 ** i1)) / (math.sqrt(vel[nm] / (1 - b2 ** i1)) + e)
+                v = self.m.params[nm] - step
+                lo, hi = self.m.param_bounds.get(nm, (-np.inf, np.inf))
+                self.m.params[nm] = float(min(max(v, lo), hi))
+            self.param_history.append(dict(self.m.params))
+
+    # ---- causal temporal weights: w_i = exp(-eps * sum_{k<i} r_k^2) ----
+    def _causal_weights(self) -> np.ndarray:
+        order = np.argsort(self.t_colloc)
+        r2 = self._phys_residual()[order] ** 2
+        cum = np.concatenate([[0.0], np.cumsum(r2)[:-1]])   # strictly earlier pts
+        w = np.exp(-EPSILON_CAUSAL * cum)
+        out = np.empty_like(w)
+        out[order] = w
+        return out
+
+    # ---- parameter gradient norm of the physics loss ----
+    def _param_grad_norm(self) -> float:
+        r = self._phys_residual()
+        x, dx, ddx = self.m.derivatives(self.t_colloc)
+        gs = []
+        for nm in self.m.params:
+            d = max(1e-6, 1e-4 * (abs(self.m.params[nm]) + 1.0))
+            p = dict(self.m.params); p[nm] += d
+            rp = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p), float)
+            gs.append(float(np.mean(r * (rp - r) / d)) * self.w_phys)
+        return float(np.linalg.norm(gs))
+
+    # ---- FIM identifiability — the SELF-DOUBT GATE ----
+    def fisher_information(self) -> Tuple[np.ndarray, np.ndarray, float]:
+        """FIM = J^T J / (Nc * sigma^2), J_{i,j} = d r_i / d eta_j at the solution.
+        A near-zero column => the data carries no information about that parameter
+        => UNIDENTIFIABLE. Returns (FIM, per-param Fisher diag, kappa(FIM))."""
+        x, dx, ddx = self.m.derivatives(self.t_colloc)
+        names = list(self.m.params)
+        nc = self.t_colloc.shape[0]
+        J = np.zeros((nc, len(names)))
+        for j, nm in enumerate(names):
+            d = max(1e-6, 1e-4 * (abs(self.m.params[nm]) + 1.0))
+            p = dict(self.m.params); p[nm] += d
+            rp = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p), float)
+            p[nm] = self.m.params[nm] - d
+            rm = np.asarray(self.m.residual_fn(self.t_colloc, x, dx, ddx, p), float)
+            J[:, j] = (rp - rm) / (2 * d)
+        sig2 = (self.noise_sigma ** 2) if self.noise_sigma else max(
+            float(np.var(self.m.surrogate.predict(self.t_data) - self.y_data)), 1e-8)
+        fim = (J.T @ J) / (nc * sig2)
+        diag = np.diag(fim).copy()
+        try:
+            s = np.linalg.svd(fim, compute_uv=False)
+            smax = float(s[0]); smin = float(s[-1])
+            kappa = (smax / smin) if smin > 0 else float("inf")
+        except np.linalg.LinAlgError:
+            kappa = float("inf")
+        return fim, diag, kappa
+
+    # ---- full fit ----
+    def fit(self, epochs: int = 600) -> ConvergenceRecord:
+        if self.t_data.shape[0] < MIN_DATA_POINTS:
+            raise ValueError(
+                f"need >= {MIN_DATA_POINTS} data points to assert anything "
+                f"(got {self.t_data.shape[0]})")
+        self.m.surrogate.fit(self.t_data, self.y_data)
+        # nonlinear params first (uses fixed surrogate derivatives), then exact LS
+        self._gd_nonlinear(epochs)
+        self._ls_solve_linear()
+        return self._build_record(epochs)
+
+    def _delta_param_rel(self, window: int = 20) -> float:
+        if len(self.param_history) < window + 1:
+            return 0.0 if self.m.linear_params else float("inf")
+        recent = self.param_history[-window:]
+        rels = []
+        for k in self.m.params:
+            if k in self.m.linear_params:
+                continue
+            vals = np.array([h[k] for h in recent])
+            denom = max(abs(np.mean(vals)), 1e-9)
+            rels.append(float(np.max(np.abs(np.diff(vals))) / denom))
+        return float(max(rels)) if rels else 0.0
+
+    def _build_record(self, epochs: int) -> ConvergenceRecord:
+        wc = self._causal_weights()
+        min_wc = float(np.min(wc)) if wc.size else 0.0
+        gnorm = self._param_grad_norm()
+        _, diag, kappa = self.fisher_information()
+        dpr = self._delta_param_rel()
+        r = self._phys_residual()
+        rms = float(np.sqrt(np.mean(r * r)))
+        d_rms = float(np.sqrt(np.mean((self.m.surrogate.predict(self.t_data) - self.y_data) ** 2)))
+        # self-doubt: any per-param Fisher below the floor forces RED.
+        min_fisher = float(np.min(diag)) if diag.size else 0.0
+        label, crit = _classify_convergence(min_wc, gnorm, kappa, dpr, min_fisher)
+        return ConvergenceRecord(
+            label=label, min_causal_weight=min_wc, grad_norm=gnorm,
+            kappa_fim=kappa, delta_param_rel=dpr, residual_rms=rms,
+            data_rms=d_rms, epochs_run=int(epochs), criteria=crit)
+
+    # ---- per-parameter governed results (value, CI, identifiability, assert) ----
+    def param_results(self, n_restarts: int = 6, epochs: int = 400) -> List[ParamResult]:
+        ci = self.ensemble_ci(n_restarts=n_restarts, epochs=epochs)
+        _, diag, kappa = self.fisher_information()
+        names = list(self.m.params)
+        out: List[ParamResult] = []
+        for j, nm in enumerate(names):
+            mean, std, lo, hi = ci[nm]
+            fisher = float(diag[j])
+            identifiable = (fisher >= FISHER_FLOOR) and (kappa < KAPPA_RED) and math.isfinite(kappa)
+            out.append(ParamResult(
+                name=nm, value=float(self.m.params[nm]),
+                ci_low=lo, ci_high=hi, std=std, fisher=fisher,
+                identifiable=identifiable, asserted=identifiable))
+        return out
+
+    # ---- ensemble CI (E-PINN): bootstrap-resample data -> refit -> resolve ----
+    def ensemble_ci(self, n_restarts: int = 6, epochs: int = 400
+                    ) -> Dict[str, Tuple[float, float, float, float]]:
+        names = list(self.m.params)
+        samples: Dict[str, List[float]] = {k: [self.m.params[k]] for k in names}
+        n = self.t_data.shape[0]
+        for s in range(1, max(1, n_restarts)):
+            idx = self.rng.integers(0, n, size=n)         # bootstrap resample
+            mdl = SZLInversePINN(
+                self.m.residual_fn,
+                {k: (self.m.params[k] if k in self.m.linear_params
+                     else self.rng.normal(self.m.params[k], 0.25 * abs(self.m.params[k]) + 0.1))
+                 for k in names},
+                surrogate=type(self.m.surrogate)(
+                    getattr(self.m.surrogate, "K", 24), getattr(self.m.surrogate, "D", 3),
+                    getattr(self.m.surrogate, "ridge", 1e-6))
+                if isinstance(self.m.surrogate, SZLSpectralSurrogate) else SZLSpectralSurrogate(),
+                param_bounds=self.m.param_bounds, linear_params=self.m.linear_params)
+            tr = SZLInversePINNTrainer(
+                mdl, self.t_data[idx], self.y_data[idx], self.t_colloc,
+                self.w_phys, self.lr_param, self.noise_sigma, seed=s + 13)
+            try:
+                tr.fit(epochs=epochs)
+                for k in names:
+                    samples[k].append(mdl.params[k])
+            except Exception:
+                pass
+        out: Dict[str, Tuple[float, float, float, float]] = {}
+        for k in names:
+            arr = np.array(samples[k], float)
+            mean = float(np.mean(arr))
+            std = float(np.std(arr, ddof=1)) if arr.size > 1 else 0.0
+            out[k] = (mean, std, mean - 1.96 * std, mean + 1.96 * std)
+        return out
+
+
+# ===========================================================================
+# 5. Three-state convergence classifier — EXACT criteria from ARXIV_LEADERS.
+# ===========================================================================
+def _classify_convergence(min_wc: float, grad_norm: float, kappa: float,
+                          delta_param_rel: float, min_fisher: float
+                          ) -> Tuple[str, Dict[str, str]]:
+    crit = {
+        "min_causal_weight": f"{min_wc:.4f} (GREEN>{CAUSAL_GREEN}, RED<={CAUSAL_RED})",
+        "grad_norm": f"{grad_norm:.2e} (GREEN<{GRAD_GREEN:.0e})",
+        "kappa_fim": f"{kappa:.2e} (IDENT<{KAPPA_IDENT:.0e}, RED>={KAPPA_RED:.0e})",
+        "min_fisher": f"{min_fisher:.2e} (floor {FISHER_FLOOR:.0e}; below=UNIDENTIFIABLE)",
+        "delta_param_rel": f"{delta_param_rel:.2e}",
+    }
+    # RED dominates: divergence, non-identifiability, or a below-floor Fisher.
+    if (min_wc <= CAUSAL_RED) or (kappa >= KAPPA_RED) or (not math.isfinite(kappa)) \
+            or (min_fisher < FISHER_FLOOR):
+        return "RED", crit
+    # GREEN requires ALL exact gates.
+    if (min_wc > CAUSAL_GREEN) and (grad_norm < GRAD_GREEN) and (kappa < KAPPA_IDENT):
+        return "GREEN", crit
+    return "YELLOW", crit
+
+
+# ===========================================================================
+# 6. Built-in Duffing system — the runnable demo physics.
+#    m x'' + c x' + delta x + alpha x^3 = F cos(omega t); alpha is the unknown.
+# ===========================================================================
+def duffing_residual(t, x, dx, ddx, params,
+                     m=1.0, c=0.2, delta=1.0, F=0.5, omega=1.0):
+    alpha = params.get("alpha", 0.0)
+    ghost = params.get("ghost", 0.0)   # enters with coefficient 0 -> UNIDENTIFIABLE
+    return (m * ddx + c * dx + delta * x + alpha * (x ** 3)
+            - F * np.cos(omega * t) + 0.0 * ghost)
+
+
+def integrate_duffing(t, m=1.0, c=0.2, delta=1.0, alpha=1.0, F=0.5, omega=1.0,
+                      x0=0.0, v0=0.0):
+    """RK4 integration of the Duffing oscillator (own NumPy integrator; no scipy
+    dependency). Returns x(t) on the given grid."""
+    t = np.asarray(t, float)
+    def deriv(state, tt):
+        x, v = state
+        a = (F * math.cos(omega * tt) - c * v - delta * x - alpha * x ** 3) / m
+        return np.array([v, a])
+    xs = np.empty_like(t)
+    state = np.array([x0, v0], float)
+    xs[0] = state[0]
+    for i in range(1, len(t)):
+        h = t[i] - t[i - 1]
+        k1 = deriv(state, t[i - 1])
+        k2 = deriv(state + 0.5 * h * k1, t[i - 1] + 0.5 * h)
+        k3 = deriv(state + 0.5 * h * k2, t[i - 1] + 0.5 * h)
+        k4 = deriv(state + h * k3, t[i])
+        state = state + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
+        xs[i] = state[0]
+    return xs