fix(inference): preflight observed-information work budgets

CHANGELOG.md

@@ -428,6 +428,11 @@
   tau update, and the covariate update while preserving the Rust production
   backend and public model contracts.
 
+#### S-X2 and person-fit Rust ownership fail-closed
+
+- Public `s_x2()` and `person_fit()` require the compiled Rust core entrypoints and no longer fall back to Python/NumPy numerical implementations when the core or symbols are missing.
+- `s_x2()` always dispatches trait `prior_mean` through the native S-X² entrypoint instead of selecting the Python reference path whenever a prior is supplied.
+
 #### Parallel-analysis input and workspace bounds
 
 - `parallel_analysis()` now rejects booleans, floats, strings, and caller-defined integer-conversion hooks for integer controls instead of silently coercing them before Rust dispatch.
@@ -535,6 +540,24 @@
 - Cap raw LLM-judge response JSON nesting at 32 levels before `json.loads`, failing closed with `JudgeFormatError` so hostile recursive objects cannot expand into parser resource exhaustion.
 - Keep valid shallow judge payloads accepted with the existing criterion/score contracts.
 
+#### Fit-statistics infit/outfit and M2 fail closed
+
+- Public `infit_outfit()` and ordinary `m2()` fail closed when the compiled Rust
+  core or required entrypoints are missing, completing the residual ownership
+  gaps from issue #627 after S-X² and person-fit hardening.
+
+#### Documentation coverage vocabulary and shipped-capability matrix
+
+- Align architecture documentation contracts with the protected-main maturity
+  vocabulary and mark parallel-analysis control bounds and essay-report native
+  dark-mode accents as ancestral after their integration.
+
+#### Observed-information work-budget preflight
+
+- Public observed-information finite-difference assembly preflights dense
+  workspace ceilings before stencil evaluation so hostile dimensions cannot
+  request unbounded Hessian work.
+
 #### Diagnostics-report focus and contrast preservation
 
 - Revealed the visually hidden diagnostics-report skip link for every actual `:focus` state while retaining the explicit `:focus-visible` treatment and strong outline.

docs/changelog.d/776-fitstats-infit-m2-failclosed.md

@@ -0,0 +1,7 @@
+# Fit-statistics infit/outfit and M2 fail closed
+
+## Fixed
+
+- Public `infit_outfit()` and ordinary `m2()` fail closed when the compiled Rust
+  core or required entrypoints are missing, completing the residual ownership
+  gaps from issue #627 after S-X² and person-fit hardening.

docs/changelog.d/777-docs-coverage-vocabulary-sync.md

@@ -0,0 +1,7 @@
+# Documentation coverage vocabulary and shipped-capability matrix
+
+## Fixed
+
+- Align architecture documentation contracts with the protected-main maturity
+  vocabulary and mark parallel-analysis control bounds and essay-report native
+  dark-mode accents as ancestral after their integration.

docs/changelog.d/781-observed-information-work-budget.md

@@ -0,0 +1,7 @@
+# Observed-information work-budget preflight
+
+## Fixed
+
+- Public observed-information finite-difference assembly preflights dense
+  workspace ceilings before stencil evaluation so hostile dimensions cannot
+  request unbounded Hessian work.

docs/documentation_coverage.md

@@ -82,10 +82,10 @@ The table below records product truth, not documentation-file presence. “Imple
 | Fixed-anchor parameter linking arithmetic | IMPLEMENTED_ON_PROTECTED_MAIN | protected main owns scale/shift estimation and theta/alpha/b transformation in Rust/PyO3 |
 | Observed-information Hessian and second-order diagnostics | IMPLEMENTED_ON_PROTECTED_MAIN | protected main owns finite-difference coefficients/symmetric Hessian assembly and eigenvalue/positive-definiteness diagnostics in Rust/PyO3; Python only evaluates objective samples and transports results |
 | JMLE Adam/L-BFGS optimizer arithmetic | IMPLEMENTED_ON_PROTECTED_MAIN | PR #760 is ancestral to current protected main; `backend="rust"` delegates Adam/L-BFGS/combined optimizer control to compiled Rust while recovery evidence remains governed separately by issue #626 |
-| Parallel-analysis public control/resource hardening | IMPLEMENTED_ON_ACTIVE_PR | current fail-first/implementation PR owns strict integer/control and bounded-workspace hardening; it remains non-shipped until exact-head integration |
+| Parallel-analysis public control/resource hardening | IMPLEMENTED_ON_PROTECTED_MAIN | strict integer/control validation and bounded random-benchmark workspace ceilings are ancestral to protected main |
 | Hourly review-repair caller | IMPLEMENTED_ON_PROTECTED_MAIN / PARTIAL | PR #763 integrated the product-side bounded caller; operational scheduler/control-plane acceptance remains external evidence rather than a library capability |
 | LLM-judge raw JSON depth hardening | IMPLEMENTED_ON_PROTECTED_MAIN | PR #764 is ancestral to current protected main and bounds recursive JSON nesting before parser materialization |
-| Essay-report native dark-mode status accents | IMPLEMENTED_ON_ACTIVE_PR | current accessibility PR owns the CSS-variable/media-query change; do not treat it as protected-main until integration |
+| Essay-report native dark-mode status accents | IMPLEMENTED_ON_PROTECTED_MAIN | CSS-variable and prefers-color-scheme dark-mode status accents are ancestral to protected main |
 | Canonical PyO3/public-export governance | ACCEPTED_ARCHITECTURE / PARTIAL | ADR-0011 governs convergence; feature-by-feature hardening continues |
 | Purpose-limited sensitive-data handling | IMPLEMENTED_ON_PROTECTED_MAIN / DOWNSTREAM | reusable contracts prefer purpose limitation/minimization/separated identities; hosted authorization/retention execution remains downstream |
 | LLM orchestration/model credentials | IMPLEMENTED_ON_PROTECTED_MAIN | provider execution and independent reviewer identity/credential boundaries are governed; provider calls remain outside psychometric numerical core |
@@ -95,9 +95,8 @@ The table below records product truth, not documentation-file presence. “Imple
 
 At this review, material open work includes:
 
-- strict/bounded public controls and Rust allocation preflight for parallel analysis;
-- the documentation-fitness refresh itself, which may describe current protected truth but is not authoritative until merged; and
-- native dark-mode report status accents.
+- Rust allocation preflight for parallel analysis, if it remains separate;
+- the documentation-fitness refresh itself, which may describe current protected truth but is not authoritative until merged.
 
 These remain active-PR evidence, not protected-main capability. Their source heads, checks, reviews, writer leases, and mergeability are operational evidence and must be re-fetched rather than copied into timeless architecture prose.
 

python/fast_mlsirm/fitstats.py

@@ -635,61 +635,39 @@ def infit_outfit(
     mask: np.ndarray | None = None,
     eps_distance: float = 1e-8,
 ) -> dict[str, np.ndarray]:
-    """Per-item infit/outfit mean squares at the EAP estimates."""
+    """Per-item infit/outfit mean squares at the EAP estimates.
+
+    Production numerical ownership is the compiled Rust core
+    (``infit_outfit_stat``). Missing or incomplete cores fail closed.
+    """
     model = model.upper()
-    free_alpha = model not in {"MLSRM", "ULSRM"}
-    uses_space = model != "MIRT"
     y, observed, d_of_i = _prepare_dichotomous_diagnostic_inputs(
         responses, factor_id, mask
     )
     core = _core_module()
-    if core is not None and hasattr(core, "infit_outfit_stat"):
-        n_persons = y.shape[0]
-        n_dims = int(d_of_i.max()) + 1
-        bank = _bank_args(params, d_of_i, model, n_dims, eps_distance)
-        res = core.infit_outfit_stat(
-            y.ravel(),
-            observed.ravel(),
-            int(n_persons),
-            bank["alpha"],
-            bank["b"],
-            bank["zeta"],
-            bank["tau"],
-            bank["factor_id"],
-            bank["model"],
-            bank["n_dims"],
-            bank["latent_dim"],
-            bank["eps_distance"],
-            np.asarray(params.theta, dtype=np.float64).ravel(),
-            np.asarray(params.xi, dtype=np.float64).ravel(),
-        )
-        return {"infit": np.asarray(res["infit"]), "outfit": np.asarray(res["outfit"])}
-    a = np.exp(params.alpha) if free_alpha else np.ones(len(params.b))
-    eta = a[None, :] * np.asarray(params.theta)[:, d_of_i] + params.b[None, :]
-    if uses_space:
-        # Optimized distance computation: replace O(N*J*D) 3D broadcast with O(N*J) 2D dot product
-        xi = np.asarray(params.xi)
-        zeta = np.asarray(params.zeta)
-        x_sq = np.einsum("ij,ij->i", xi, xi)
-        z_sq = np.einsum("ij,ij->i", zeta, zeta)
-        dist_sq = x_sq[:, None] + z_sq[None, :] - 2 * np.dot(xi, zeta.T)
-        dist = np.sqrt(eps_distance + np.maximum(dist_sq, 0.0))
-        eta = eta - math.exp(params.tau) * dist
-    p = np.clip(1.0 / (1.0 + np.exp(-np.clip(eta, -700, 700))), 1e-12, 1 - 1e-12)
-    v = p * (1.0 - p)
-    resid2 = np.subtract(y, p)
-    np.square(resid2, out=resid2)
-    np.multiply(resid2, observed, out=resid2)
-    n_obs = np.maximum(observed.sum(axis=0), 1)
-
-    # Preserve the masked squared-residual numerator, then reuse its owned
-    # float64 buffer for the outfit division without a numeric mask copy.
-    resid2_sum = resid2.sum(axis=0)
-    infit_denominator = np.sum(v, axis=0, where=observed)
-    np.divide(resid2, v, out=resid2)
-    outfit = resid2.sum(axis=0) / n_obs
-    infit = resid2_sum / np.maximum(infit_denominator, 1e-12)
-    return {"infit": infit, "outfit": outfit}
+    if core is None or not hasattr(core, "infit_outfit_stat"):
+        raise RuntimeError("fit statistics require the compiled Rust core")
+    n_persons = y.shape[0]
+    n_dims = int(d_of_i.max()) + 1
+    bank = _bank_args(params, d_of_i, model, n_dims, eps_distance)
+    res = core.infit_outfit_stat(
+        y.ravel(),
+        observed.ravel(),
+        int(n_persons),
+        bank["alpha"],
+        bank["b"],
+        bank["zeta"],
+        bank["tau"],
+        bank["factor_id"],
+        bank["model"],
+        bank["n_dims"],
+        bank["latent_dim"],
+        bank["eps_distance"],
+        np.asarray(params.theta, dtype=np.float64).ravel(),
+        np.asarray(params.xi, dtype=np.float64).ravel(),
+    )
+    return {"infit": np.asarray(res["infit"]), "outfit": np.asarray(res["outfit"])}
+
 
 
 # --------------------------------------------------------------------------
@@ -1851,56 +1829,44 @@ def m2(
         )
 
     core = _core_module()
-    if core is not None and hasattr(core, "m2_stat"):
-        bank = _bank_args(params, d_of_i, model, n_dims, eps_distance)
-        res = core.m2_stat(
-            np.where(observed0, y0, 0.0).ravel(),
-            observed0.ravel(),
-            int(y0.shape[0]),
-            bank["alpha"],
-            bank["b"],
-            bank["zeta"],
-            bank["tau"],
-            bank["factor_id"],
-            bank["model"],
-            bank["n_dims"],
-            bank["latent_dim"],
-            bank["eps_distance"],
-            prior_mean,
-            prior_sd,
-            q_theta=int(q_theta),
-            xi_rule="gh",
-            q_xi=int(q_xi),
-        )
-        result = M2Result(
-            m2=float(res["m2"]),
-            df=float(res["df"]),
-            p_value=float(res["p_value"]),
-            rmsea2=float(res["rmsea2"]),
-            rmsea2_ci_lower=float(res["rmsea2_ci_lower"]),
-            rmsea2_ci_upper=float(res["rmsea2_ci_upper"]),
-            srmsr=float(res["srmsr"]),
-            null_m2=float(res["null_m2"]),
-            null_df=float(res["null_df"]),
-            cfi=float(res["cfi"]),
-            tli=float(res["tli"]),
-            n_moments=int(res["n_moments"]),
-            n_parameters=int(res["n_parameters"]),
-            n_complete=int(res["n_complete"]),
-        )
-    else:
-        result = _m2_numpy(
-            y0,
-            observed0,
-            d_of_i,
-            params,
-            model,
-            q_theta,
-            q_xi,
-            eps_distance,
-            prior_mean,
-            prior_sd,
-        )
+    if core is None or not hasattr(core, "m2_stat"):
+        raise RuntimeError("fit statistics require the compiled Rust core")
+    bank = _bank_args(params, d_of_i, model, n_dims, eps_distance)
+    res = core.m2_stat(
+        np.where(observed0, y0, 0.0).ravel(),
+        observed0.ravel(),
+        int(y0.shape[0]),
+        bank["alpha"],
+        bank["b"],
+        bank["zeta"],
+        bank["tau"],
+        bank["factor_id"],
+        bank["model"],
+        bank["n_dims"],
+        bank["latent_dim"],
+        bank["eps_distance"],
+        prior_mean,
+        prior_sd,
+        q_theta=int(q_theta),
+        xi_rule="gh",
+        q_xi=int(q_xi),
+    )
+    result = M2Result(
+        m2=float(res["m2"]),
+        df=float(res["df"]),
+        p_value=float(res["p_value"]),
+        rmsea2=float(res["rmsea2"]),
+        rmsea2_ci_lower=float(res["rmsea2_ci_lower"]),
+        rmsea2_ci_upper=float(res["rmsea2_ci_upper"]),
+        srmsr=float(res["srmsr"]),
+        null_m2=float(res["null_m2"]),
+        null_df=float(res["null_df"]),
+        cfi=float(res["cfi"]),
+        tli=float(res["tli"]),
+        n_moments=int(res["n_moments"]),
+        n_parameters=int(res["n_parameters"]),
+        n_complete=int(res["n_complete"]),
+    )
     if estimator == "mmle":
         return result
     result.estimator = estimator

python/fast_mlsirm/inference.py

@@ -7,6 +7,48 @@
 from .objective import neg_loglik_and_grad
 from .types import MLSIRMParams
 
+# Package-owned support ceilings for the O(n^2) finite-difference observed
+# information path. These are implementation safety limits, not psychometric
+# recommendations or universal hardware-capacity claims.
+_MAX_OBSERVED_INFORMATION_OBJECTIVE_CALLS = 250_001
+_MAX_OBSERVED_INFORMATION_WORKSPACE_BYTES = 128 * 1024 * 1024
+
+
+def _observed_information_work(n: int) -> tuple[int, int]:
+    """Return exact objective-call and fixed-width workspace requirements.
+
+    The central-difference stencil evaluates the objective once at the base,
+    twice per diagonal, and four times per off-diagonal pair, for exactly
+    ``1 + 2*n**2`` calls. Fixed-width workspace accounts conservatively for the
+    finite-difference value arrays, one reusable perturbation vector, and the
+    dense Rust-owned result while those inputs remain live.
+    """
+    n = int(n)
+    if n < 0:
+        raise ValueError("observed-information parameter count must be non-negative")
+    objective_calls = 1 + 2 * n * n
+    float64_bytes = np.dtype(np.float64).itemsize
+    workspace_values = 3 * n * n + n
+    workspace_bytes = workspace_values * float64_bytes
+    return objective_calls, workspace_bytes
+
+
+def _preflight_observed_information(n: int) -> None:
+    """Reject unsupported finite-difference work before objective evaluation."""
+    objective_calls, workspace_bytes = _observed_information_work(n)
+    if objective_calls > _MAX_OBSERVED_INFORMATION_OBJECTIVE_CALLS:
+        raise ValueError(
+            "observed_information objective-call budget exceeded: "
+            f"requires {objective_calls} calls, at most "
+            f"{_MAX_OBSERVED_INFORMATION_OBJECTIVE_CALLS} are supported"
+        )
+    if workspace_bytes > _MAX_OBSERVED_INFORMATION_WORKSPACE_BYTES:
+        raise ValueError(
+            "observed_information workspace budget exceeded: "
+            f"requires {workspace_bytes} bytes, limit is "
+            f"{_MAX_OBSERVED_INFORMATION_WORKSPACE_BYTES}"
+        )
+
 
 def observed_information(
     responses: np.ndarray,
@@ -23,6 +65,10 @@ def observed_information(
     The default Rust device is CPU so finite-difference curvature uses the f64
     path even when model fitting defaults to ``rust_device="auto"`` on GPU hosts.
     Pass ``device=None`` to honor ``config.rust_device`` instead.
+
+    The dense finite-difference path is preflighted against package-owned
+    objective-call and fixed-width workspace budgets before the first objective
+    evaluation. Exceeding either support ceiling fails closed with ``ValueError``.
     """
     config = config or FitConfig()
     model = config.normalized_model()
@@ -31,6 +77,9 @@ def observed_information(
     if not np.isfinite(step) or step <= 0:
         raise ValueError("step must be > 0 and finite")
 
+    n = x0.size
+    _preflight_observed_information(n)
+
     def objective(x: np.ndarray) -> float:
         """Return the penalized negative log-likelihood at packed parameter vector ``x``."""
         value, _, _ = neg_loglik_and_grad(
@@ -46,35 +95,42 @@ def objective(x: np.ndarray) -> float:
             raise ValueError("objective must be finite for Hessian calculation")
         return float(value)
 
-    n = x0.size
-    MAX_HESSIAN_DIM = 5_000
-    if n > MAX_HESSIAN_DIM:
-        raise ValueError(
-            f"observed_information supports at most {MAX_HESSIAN_DIM} parameters (got {n}); "
-            "the dense finite-difference Hessian is O(n^2) memory and O(n^2) objective calls"
-        )
     base = objective(x0)
-    eye = np.eye(n, dtype=np.float64)
     h = float(step)
 
     # Python evaluates the scalar objective at FD offsets; Rust owns the
-    # finite-difference coefficients and symmetrised matrix assembly.
+    # finite-difference coefficients and symmetrised matrix assembly. One
+    # reusable trial vector replaces the former dense n x n identity matrix.
     diag_plus = np.empty(n, dtype=np.float64)
     diag_minus = np.empty(n, dtype=np.float64)
     off_n = n * (n - 1) // 2
     off_pp = np.empty(off_n, dtype=np.float64)
     off_pm = np.empty(off_n, dtype=np.float64)
     off_mp = np.empty(off_n, dtype=np.float64)
     off_mm = np.empty(off_n, dtype=np.float64)
+    trial = np.array(x0, dtype=np.float64, copy=True)
     k = 0
     for i in range(n):
-        diag_plus[i] = objective(x0 + h * eye[i])
-        diag_minus[i] = objective(x0 - h * eye[i])
+        base_i = float(x0[i])
+        trial[i] = base_i + h
+        diag_plus[i] = objective(trial)
+        trial[i] = base_i - h
+        diag_minus[i] = objective(trial)
+        trial[i] = base_i
         for j in range(i + 1, n):
-            off_pp[k] = objective(x0 + h * eye[i] + h * eye[j])
-            off_pm[k] = objective(x0 + h * eye[i] - h * eye[j])
-            off_mp[k] = objective(x0 - h * eye[i] + h * eye[j])
-            off_mm[k] = objective(x0 - h * eye[i] - h * eye[j])
+            base_j = float(x0[j])
+            trial[i] = base_i + h
+            trial[j] = base_j + h
+            off_pp[k] = objective(trial)
+            trial[j] = base_j - h
+            off_pm[k] = objective(trial)
+            trial[i] = base_i - h
+            trial[j] = base_j + h
+            off_mp[k] = objective(trial)
+            trial[j] = base_j - h
+            off_mm[k] = objective(trial)
+            trial[i] = base_i
+            trial[j] = base_j
             k += 1
 
     from . import _core as core