From 5f038243aec0aa2f3f3fc284f89552a2b47c86e9 Mon Sep 17 00:00:00 2001 From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com> Date: Sun, 30 Aug 2026 06:46:28 +0800 Subject: [PATCH 1/3] Add conditional query information for correlated partial gauges --- .../conditional-query-information.yml | 66 +++ docs/conditional-query-information.md | 197 +++++++++ .../conditional-query-design-study-v1.json | 25 ++ src/prob4d/conditional_gauge_design.py | 400 ++++++++++++++++++ src/prob4d/conditional_gauge_study.py | 255 +++++++++++ tests/test_conditional_gauge_design.py | 245 +++++++++++ tests/test_conditional_gauge_integration.py | 51 +++ 7 files changed, 1239 insertions(+) create mode 100644 .github/workflows/conditional-query-information.yml create mode 100644 docs/conditional-query-information.md create mode 100644 protocols/conditional-query-design-study-v1.json create mode 100644 src/prob4d/conditional_gauge_design.py create mode 100644 src/prob4d/conditional_gauge_study.py create mode 100644 tests/test_conditional_gauge_design.py create mode 100644 tests/test_conditional_gauge_integration.py diff --git a/.github/workflows/conditional-query-information.yml b/.github/workflows/conditional-query-information.yml new file mode 100644 index 000000000..79c80097d --- /dev/null +++ b/.github/workflows/conditional-query-information.yml @@ -0,0 +1,66 @@ +name: Conditional query information + +on: + pull_request: + branches: [main] + paths: + - ".github/workflows/conditional-query-information.yml" + - "src/prob4d/conditional_gauge_design.py" + - "src/prob4d/conditional_gauge_study.py" + - "src/prob4d/observable_gauge.py" + - "src/prob4d/query_observability.py" + - "tests/test_conditional_gauge_design.py" + - "tests/test_conditional_gauge_integration.py" + - "protocols/conditional-query-design-study-v1.json" + - "docs/conditional-query-information.md" + workflow_dispatch: + +permissions: + contents: read + +concurrency: + group: conditional-query-information-${{ github.ref }} + cancel-in-progress: true + +jobs: + conditional-information: + name: Conditional likelihood, query design, and observable-factor parity + runs-on: ubuntu-latest + timeout-minutes: 15 + steps: + - name: Check out exact revision + uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - name: Set up Python + uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0 + with: + python-version: "3.12" + - name: Install development environment + run: python -m pip install -e ".[dev]" + - name: Check source formatting + run: python -m ruff format --diff src/prob4d/conditional_gauge_{design,study}.py tests/test_conditional_gauge_{design,integration}.py + - name: Check source lint + if: ${{ !cancelled() }} + run: python -m ruff check src/prob4d/conditional_gauge_{design,study}.py tests/test_conditional_gauge_{design,integration}.py + - name: Check source typing + if: ${{ !cancelled() }} + run: python -m mypy src/prob4d/conditional_gauge_design.py src/prob4d/conditional_gauge_study.py + - name: Run conditional and existing observability tests + if: ${{ !cancelled() }} + run: | + python -m pytest -q tests/test_conditional_gauge_design.py tests/test_conditional_gauge_integration.py tests/test_observable_gauge.py tests/test_query_observability.py + - name: Reproduce complete controlled study + if: ${{ !cancelled() }} + run: | + python -m prob4d.conditional_gauge_study --output "$RUNNER_TEMP/conditional-query-result.json" + python - <<'PY' + import json + import os + from pathlib import Path + from prob4d.conditional_gauge_study import PROTOCOL + assert json.loads(Path("protocols/conditional-query-design-study-v1.json").read_text()) == PROTOCOL + result = json.loads((Path(os.environ["RUNNER_TEMP"]) / "conditional-query-result.json").read_text()) + assert result["maximum_kernel_vs_independent_dense_reference_error"] < 1e-10 + print(json.dumps(result, indent=2, sort_keys=True)) + PY diff --git a/docs/conditional-query-information.md b/docs/conditional-query-information.md new file mode 100644 index 000000000..52fd2c49a --- /dev/null +++ b/docs/conditional-query-information.md @@ -0,0 +1,197 @@ +# Conditional query information for correlated partial gauges + +## Scientific question + +Which additional 4-D prediction window is useful **after** the already consumed +windows, for one declared downstream physical query? + +An overlap's standalone covariance is not its incremental information. Repeated +source evidence can look precise without adding information; a geometrically +complementary window can be more useful; reducing total gauge uncertainty can be +irrelevant to the query. This experimental module joins the existing observable +Sim(3) factor and query-observability work with explicit cross-window dependence. +It changes no stable exporter, production guard, evidence gate, or target access. + +## Model and conditional likelihood + +Work in one fixed, centroid-normalized, seven-dimensional local gauge chart: + +\[ +x\sim\mathcal N(m_0,P_0),\qquad +\begin{bmatrix}y_H\\y_c\end{bmatrix} +=\begin{bmatrix}H_H\\H_c\end{bmatrix}x+ +\begin{bmatrix}\epsilon_H\\\epsilon_c\end{bmatrix},\qquad +\operatorname{Cov}(\epsilon)= +\begin{bmatrix}R_{HH}&R_{Hc}\\R_{cH}&R_{cc}\end{bmatrix}. +\] + +The complete prior is independent of this observation noise. The joint noise +covariance is known/frozen for this calculation; it is **not** a predictive +covariance that already contains state uncertainty. For nonsingular history +noise, define + +\[ +L=R_{cH}R_{HH}^{-1},\quad +\widetilde y_c=y_c-Ly_H,\quad +\widetilde H_c=H_c-LH_H,\quad +\widetilde R_c=R_{cc}-LR_{Hc}. +\] + +Then the new likelihood is + +\[ +p(y_c\mid x,y_H) +=\mathcal N(\widetilde y_c;\widetilde H_cx,\widetilde R_c). +\] + +The implementation whitens only the supported noise subspace. It neither +inverts a geometric nullspace nor fills it with a ridge. Singular noise is +accepted only when its zero-noise directions carry no state information and the +observations satisfy their deterministic identities. Informative noiseless +constraints require a constrained solver and are explicitly rejected, rather +than silently discarded by a pseudoinverse. Numerical support is relative to +the declared covariance scale, with `rtol=1e-10` by default. + +### Proposition: exact incremental query value + +Let the actual history posterior be `(m_H, P_H)`, and let `q=Jx` be the fixed +local query with positive-definite output metric `W`. For positive-definite +conditional innovation covariance +`C = H_tilde P_H H_tilde.T + R_tilde`, + +\[ +\Delta_c=\operatorname{tr}\!\left[ + WJ P_H\widetilde H_c^T C^{-1}\widetilde H_c P_H J^T\right]\geq0. +\] + +This is the expected reduction of squared query error when using the posterior +mean, conditional on the existing history and under the stated Gaussian model. +It is zero exactly when `J P_H H_tilde.T = 0`. It can be computed before reading +any candidate outcome. The supported-subspace implementation extends the +calculation to compatible, redundant singular observations. + +**Proof.** Gaussian conditioning gives +`P_new = P_H - P_H H_tilde.T C^-1 H_tilde P_H`. The posterior mean's conditional +Bayes risk for the metric-squared query loss is `trace(W J P J.T)`. Subtracting +risks gives the displayed trace, equivalently the squared Frobenius norm of +`W^(1/2) J P_H H_tilde.T C^(-1/2)`. Positive definiteness proves the zero-gain +criterion. This is an application of standard Gaussian conditioning and +quantity-of-interest design, not a claim to have invented those identities. + +### Corollary: source-replay invariance + +If a candidate is a deterministic replay `y_c=T y_H`, with matching designs and +noise covariance blocks, both its conditional design and conditional noise are +zero. Its incremental information and query value are zero. The session returns +the identical `GaussianGaugeBelief` object, not a numerically reconstructed copy. +An inconsistent replay is rejected before session state is changed. + +Information nullspaces and prior correlations remain distinct. A factor does +not add precision along its nullspace, but a correlated prior can legitimately +transmit information to that direction. The implementation does not incorrectly +force nullspace marginal variances to stay constant for every prior. + +### Counterexample: standalone filtering and greedy guarantees + +Consider `x ~ N(0,1)`, `y_A=x+b`, and `y_B=b+eta`, with independent +`b ~ N(0,1)` and `eta ~ N(0,0.1)`. Window B has no standalone state information. +Yet after A, it changes posterior variance from `1/2` to `1/12`: a gain of +`5/12`. It estimates the shared nuisance rather than directly measuring x. + +Thus discarding every candidate with zero standalone query information is not +valid under dependence. Also, variance-reduction set utility is not generally +submodular: B's gain increases after A. The implemented selector is an exact +**one-step** choice among the supplied candidates, not a globally optimal +multi-window scheduler or a greedy approximation guarantee. This analytic +reference-channel control does not claim that a particular provider supplies +such a channel. + +## API and integration + +`CorrelatedGaugeDesign` holds window row partitions, common-chart design rows, +and the full externally specified joint noise covariance. `ConditionalGaugeSession` +owns the actually consumed history. `preview_query` has no candidate-outcome +argument. `select_query_window` ranks expected reduction per positive cost and +returns `None` if no candidate clears the frozen minimum. + +For an existing `ObservableGaugeFactor` in **the same chart**, use +`factor.observable_basis.T` as the observation design, +`factor.observable_covariance` as its marginal noise covariance, and zero as its +factor-centered observation. The integration test verifies exact parity with +`fuse_local_gaussian` and `evaluate_query_observability`, then verifies that +replaying the same source evidence changes nothing. Cross-window covariance +blocks still need independent justification; diagonal blocks alone are not +sufficient. Factors from different linearization charts must first be transported +to the declared common chart; merely assigning the same text identifier is not a +coordinate transformation. + +A session must begin at a prior that has not consumed any of its windows. Replay +its admitted history through the session rather than wrapping an already updated +BayesianPhysTwin belief as an untouched prior. Downstream BayesianPhysTwin still +owns complete physical-belief construction and its exact caller-owned fallback. +The identity guarantee here applies to the local gauge-belief object only. +Causal4D query Jacobians can define a future endpoint or contrast objective, but +this module does not execute Causal4D or establish causal benefit. + +The kernel is dense in the compressed window rows. It is a small research +implementation, not a demonstrated large-stream sparse solver. Known chart, +state-independent Gaussian covariance, first-order query validity, provider +support, identity and lineage, calibration, physical consistency, and downstream +guards remain separate requirements. A positive utility does not authorize a +real update or a target evaluation. + +## Reproducible controlled study + +```bash +python -m pytest -q tests/test_conditional_gauge_design.py \ + tests/test_conditional_gauge_integration.py +python -m prob4d.conditional_gauge_study --output /tmp/conditional-query-result.json +``` + +The frozen numerical design is in +`protocols/conditional-query-design-study-v1.json`. The study generates 10,000 +independent Gaussian episodes, a rank-six line-like history, a 0.999-correlated +repeat, complementary rank-six support, and a precise but query-irrelevant +rotation measurement. Each policy receives one additional window at equal cost. +Selections are made before any episode is drawn. The point-position query uses +a fixed local Sim(3) Jacobian, not a nonlinear scene renderer or physical solver. + +Five arms separate selection from covariance accounting: history only, marginal +query selection with independent updates, marginal selection with correct +updates, global gauge-variance selection with correct updates, and conditional +query selection with correct updates. Report Euclidean query RMSE, Gaussian NLL, +normalized query NEES, 90% ellipsoid coverage, harmful-episode fraction, and paired +episode-bootstrap squared-loss differences. An exact-replay sweep and the +reference-channel counterexample are analytic controls. An independently computed +dense joint posterior checks the conditional implementation. + +Paper-facing generated results and exact source manifests belong in +`FlorianPfaff/BayesianPhysTwin-Paper`, not in the public release claim table. +This is a constructed mechanism study, not a fresh real-provider evaluation, +independent calibration, safety result, or proof of state of the art. Individual +episodes may become worse despite a positive expected value; those fractions +must be reported rather than suppressed. + +## Relation to prior work and the larger paper + +Partial-factor/localizability treatment already exists in registration and +SLAM. Quantity-of-interest optimal design and sensor selection with correlated +noise also predate this work. The proposed Prob4D contribution is the narrower +combination: partial learned-window gauge likelihoods, conditional source-evidence +accounting, decision-relevant query value, and replay/selection failure controls +at the physical-twin interface. A literature review and a fresh real-provider +study are still needed before claiming a broadly novel or empirically validated +system. This work does not reopen any terminal MotionCrafter, CUT3R, or Deform360 +cohort and does not require collecting new physical data. + +Primary references: + +- Tuna et al., *X-ICP: Localizability-Aware LiDAR Registration for Robust + Localization in Extreme Environments*, IEEE T-RO, 2024; + arXiv:2211.16335, https://arxiv.org/abs/2211.16335. +- Attia, Alexanderian, and Saibaba, *Goal-Oriented Optimal Design of Experiments + for Large-Scale Bayesian Linear Inverse Problems*, Inverse Problems, 2018; + arXiv:1802.06517, https://arxiv.org/abs/1802.06517. +- Liu et al., *Sensor Selection for Estimation with Correlated Measurement + Noise*, IEEE T-SP 64(13):3509--3522, 2016; + DOI:10.1109/TSP.2016.2550005, https://arxiv.org/abs/1508.03690. diff --git a/protocols/conditional-query-design-study-v1.json b/protocols/conditional-query-design-study-v1.json new file mode 100644 index 000000000..8bb265319 --- /dev/null +++ b/protocols/conditional-query-design-study-v1.json @@ -0,0 +1,25 @@ +{ + "bootstrap_replicates": 2000, + "bootstrap_seed": 8302026, + "candidate_cost": 1.0, + "classification": "constructed-linear-Gaussian-mechanism-study", + "cloud_scale_m": 0.1, + "complement_noise_variance": 0.00025, + "ellipsoid_90_chi_square_df3": 6.251388631170325, + "episodes": 10000, + "global_only_noise_variance": 1e-06, + "history_noise_standard_deviation": 0.01, + "independent_unit": "synthetic episode, not coordinate or point", + "new_window_budget": 1, + "prior_standard_deviation": 0.05, + "provider_executed": false, + "query_point_in_cloud_radii": [ + 0.0, + 0.02, + 0.01 + ], + "repeat_noise_correlation": 0.999, + "schema": "prob4d.conditional-query-design-study-v1", + "seed": 20260830, + "targets_opened": false +} diff --git a/src/prob4d/conditional_gauge_design.py b/src/prob4d/conditional_gauge_design.py new file mode 100644 index 000000000..07d0288a1 --- /dev/null +++ b/src/prob4d/conditional_gauge_design.py @@ -0,0 +1,400 @@ +"""Conditional query information for correlated, partial 4-D gauge windows. + +Experimental, linear-Gaussian research API. All rows must refer to one declared +seven-dimensional gauge chart. A full, externally calibrated joint *noise* +covariance is required: absence of cross-covariance is not independence. + +The implementation conditions candidate noise on the actually assimilated +history. It never substitutes the candidate's marginal information, inserts a +ridge into a geometric nullspace, or uses candidate outcomes to select a window. +A utility is an expected squared-query-loss reduction, not a safety certificate. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any + +import numpy as np +from numpy.typing import NDArray + +FloatArray = NDArray[np.floating[Any]] + + +class UnsupportedDeterministicConstraint(ValueError): + """A zero-noise direction constrains the state; use a constrained solver.""" + + +def _array(value: object, *, name: str, ndim: int) -> FloatArray: + result = np.asarray(value, dtype=np.float64).copy() + if result.ndim != ndim or not np.all(np.isfinite(result)): + raise ValueError(f"{name} must be a finite {ndim}-dimensional array") + result.setflags(write=False) + return result + + +def _symmetric(value: object, *, name: str, size: int) -> FloatArray: + result = _array(value, name=name, ndim=2) + if result.shape != (size, size): + raise ValueError(f"{name} must have shape ({size}, {size})") + scale = max(float(np.max(np.abs(result), initial=0.0)), np.finfo(float).tiny) + if float(np.max(np.abs(result - result.T), initial=0.0)) > 1e-10 * scale: + raise ValueError(f"{name} must be symmetric") + return _array((result + result.T) / 2, name=name, ndim=2) + + +def _positive_definite(value: object, *, name: str, size: int) -> FloatArray: + result = _symmetric(value, name=name, size=size) + try: + np.linalg.cholesky(result) + except np.linalg.LinAlgError as exc: + raise ValueError(f"{name} must be positive definite") from exc + return result + + +def _text(value: object, name: str) -> str: + if not isinstance(value, str) or not value.strip(): + raise ValueError(f"{name} must be a nonempty string") + return value + + +def _noise_support( + covariance: FloatArray, + design: FloatArray, + *, + rtol: float, + reference_noise_scale: float = 0.0, + reference_design_scale: float = 0.0, +) -> tuple[FloatArray, FloatArray]: + """Return a supported whitener and zero-noise basis, without a ridge.""" + if covariance.shape[0] == 0: + return np.empty((0, 0)), np.empty((0, 0)) + eigenvalues, eigenvectors = np.linalg.eigh(covariance) + scale = max(float(np.max(np.abs(eigenvalues))), reference_noise_scale) + threshold = rtol * scale + if float(eigenvalues[0]) < -threshold: + raise ValueError("joint or conditional noise covariance is not positive semidefinite") + supported = eigenvalues > threshold + zero_basis = eigenvectors[:, ~supported] + design_scale = max(float(np.linalg.norm(design)), reference_design_scale) + if float(np.linalg.norm(zero_basis.T @ design)) > rtol * design_scale: + raise UnsupportedDeterministicConstraint( + "zero-noise direction contains state information; no pseudoinverse rescue" + ) + whitener = (eigenvectors[:, supported] / np.sqrt(eigenvalues[supported])).T + return whitener, zero_basis + + +def _covariance_update(covariance: FloatArray, design: FloatArray) -> FloatArray: + """Joseph-form update for unit-noise supported measurements.""" + if not np.any(design): + return covariance + innovation = np.eye(design.shape[0]) + design @ covariance @ design.T + gain = np.linalg.solve(innovation, design @ covariance).T + residual = np.eye(7) - gain @ design + updated = residual @ covariance @ residual.T + gain @ gain.T + return _positive_definite(updated, name="posterior covariance", size=7) + + +@dataclass(frozen=True) +class GaussianGaugeBelief: + """Complete Gaussian belief in the caller's declared local gauge chart.""" + + chart_id: str + mean: FloatArray + covariance: FloatArray + + def __post_init__(self) -> None: + _text(self.chart_id, "chart_id") + mean = _array(self.mean, name="mean", ndim=1) + if mean.shape != (7,): + raise ValueError("mean must have shape (7,)") + covariance = _positive_definite(self.covariance, name="covariance", size=7) + object.__setattr__(self, "mean", mean) + object.__setattr__(self, "covariance", covariance) + + +@dataclass(frozen=True) +class ConditionalWindowFactor: + """An outcome-free conditional likelihood, including singular-noise support.""" + + candidate_id: str + history_ids: tuple[str, ...] + regression: FloatArray + design: FloatArray + noise_covariance: FloatArray + whitener: FloatArray + zero_noise_basis: FloatArray + history_zero_noise_basis: FloatArray + noise_scale: float + rtol: float + + @property + def whitened_design(self) -> FloatArray: + return self.whitener @ self.design + + @property + def information_matrix(self) -> FloatArray: + matrix = self.whitened_design + return matrix.T @ matrix + + @property + def information_rank(self) -> int: + singular = np.linalg.svd(self.whitened_design, compute_uv=False) + if singular.size == 0 or singular[0] == 0: + return 0 + return int(np.count_nonzero(singular > self.rtol * singular[0])) + + def whitened_value(self, history: FloatArray, candidate: FloatArray) -> FloatArray: + """Validate zero-noise identities before using supported residuals.""" + residual = candidate - self.regression @ history + scale = max( + float(np.linalg.norm(history)), + float(np.linalg.norm(candidate)), + float(np.linalg.norm(self.regression @ history)), + float(np.sqrt(self.noise_scale)), + ) + if ( + float(np.linalg.norm(self.history_zero_noise_basis.T @ history)) + > self.rtol * scale + or float(np.linalg.norm(self.zero_noise_basis.T @ residual)) > self.rtol * scale + ): + raise ValueError("observations violate a deterministic noise-support identity") + return self.whitener @ residual + + +@dataclass(frozen=True) +class CorrelatedGaugeDesign: + """Frozen window designs and their full, known joint observation noise. + + ``design_matrix`` maps one common local gauge vector to stacked measurements. + ``noise_covariance`` is conditional on that vector, NOT a predictive covariance + containing state uncertainty. Window IDs and sizes partition its row axis. + ``dependence_id`` names the externally fitted covariance model; it is metadata, + not independent evidence that the covariance is calibrated. + + Singular covariance is supported only when its zero-noise directions contain + no state information (e.g. exact repeated source evidence). Informative exact + constraints fail closed rather than being dropped by a pseudoinverse. + """ + + chart_id: str + dependence_id: str + window_ids: tuple[str, ...] + window_sizes: tuple[int, ...] + design_matrix: FloatArray + noise_covariance: FloatArray + rtol: float = 1e-10 + + def __post_init__(self) -> None: + _text(self.chart_id, "chart_id") + _text(self.dependence_id, "dependence_id") + ids = tuple(self.window_ids) + sizes = tuple(self.window_sizes) + if not ids or len(ids) != len(sizes) or len(set(ids)) != len(ids): + raise ValueError("window IDs must be nonempty, unique, and match window sizes") + for window_id in ids: + _text(window_id, "window ID") + if any( + isinstance(size, bool) + or not isinstance(size, (int, np.integer)) + or size < 1 + for size in sizes + ): + raise ValueError("window sizes must be positive integers") + if isinstance(self.rtol, bool) or not np.isfinite(self.rtol) or not 0 < self.rtol < 1: + raise ValueError("rtol must lie strictly between zero and one") + design = _array(self.design_matrix, name="design_matrix", ndim=2) + if design.shape != (sum(sizes), 7): + raise ValueError("design_matrix must have shape (sum(window_sizes), 7)") + covariance = _symmetric(self.noise_covariance, name="noise_covariance", size=sum(sizes)) + _noise_support(covariance, design, rtol=self.rtol) + object.__setattr__(self, "window_ids", ids) + object.__setattr__(self, "window_sizes", sizes) + object.__setattr__(self, "design_matrix", design) + object.__setattr__(self, "noise_covariance", covariance) + + def rows(self, window_ids: tuple[str, ...]) -> tuple[int, ...]: + if len(set(window_ids)) != len(window_ids): + raise ValueError("a window cannot appear twice in a history") + result: list[int] = [] + for window_id in window_ids: + if window_id not in self.window_ids: + raise ValueError(f"unknown window: {window_id}") + index = self.window_ids.index(window_id) + start = sum(self.window_sizes[:index]) + result.extend(range(start, start + self.window_sizes[index])) + return tuple(result) + + def conditional_factor( + self, history_ids: tuple[str, ...], candidate_id: str + ) -> ConditionalWindowFactor: + """Construct p(y_candidate | x, y_history) without reading any y.""" + if candidate_id in history_ids: + raise ValueError("candidate has already been assimilated") + old_rows = self.rows(history_ids) + new_rows = self.rows((candidate_id,)) + old_design = self.design_matrix[list(old_rows)] + new_design = self.design_matrix[list(new_rows)] + old_noise = self.noise_covariance[np.ix_(old_rows, old_rows)] + cross_noise = self.noise_covariance[np.ix_(new_rows, old_rows)] + new_noise = self.noise_covariance[np.ix_(new_rows, new_rows)] + old_whitener, old_zero = _noise_support(old_noise, old_design, rtol=self.rtol) + old_inverse = old_whitener.T @ old_whitener + regression = cross_noise @ old_inverse + design = new_design - regression @ old_design + design_scale = max( + float(np.linalg.norm(new_design)), + float(np.linalg.norm(regression @ old_design)), + ) + if float(np.linalg.norm(design)) <= self.rtol * design_scale: + design = np.zeros_like(design) + noise = new_noise - regression @ cross_noise.T + noise = (noise + noise.T) / 2 + noise_scale = float(np.max(np.linalg.eigvalsh(new_noise), initial=0.0)) + whitener, zero = _noise_support( + noise, + design, + rtol=self.rtol, + reference_noise_scale=noise_scale, + reference_design_scale=design_scale, + ) + return ConditionalWindowFactor( + candidate_id=candidate_id, + history_ids=tuple(history_ids), + regression=_array(regression, name="regression", ndim=2), + design=_array(design, name="conditional design", ndim=2), + noise_covariance=_array(noise, name="conditional noise", ndim=2), + whitener=_array(whitener, name="whitener", ndim=2), + zero_noise_basis=_array(zero, name="zero_noise_basis", ndim=2), + history_zero_noise_basis=_array(old_zero, name="history_zero_noise_basis", ndim=2), + noise_scale=noise_scale, + rtol=self.rtol, + ) + + +@dataclass(frozen=True) +class QueryWindowUtility: + """Expected metric-squared query-loss reduction, not realized benefit.""" + + candidate_id: str + conditional_information_rank: int + prior_metric_variance: float + posterior_metric_variance: float + cost: float + + def __post_init__(self) -> None: + _text(self.candidate_id, "candidate_id") + rank = self.conditional_information_rank + if isinstance(rank, bool) or not isinstance(rank, (int, np.integer)) or not 0 <= rank <= 7: + raise ValueError("conditional_information_rank must be an integer in [0, 7]") + for name in ("prior_metric_variance", "posterior_metric_variance"): + value = getattr(self, name) + if isinstance(value, bool) or not np.isfinite(value) or value < 0: + raise ValueError(f"{name} must be finite and nonnegative") + if isinstance(self.cost, bool) or not np.isfinite(self.cost) or self.cost <= 0: + raise ValueError("cost must be finite and positive") + scale = max(self.prior_metric_variance, np.finfo(float).tiny) + if self.posterior_metric_variance > self.prior_metric_variance + 1e-10 * scale: + raise ValueError("posterior query variance must not exceed prior variance") + + @property + def variance_reduction(self) -> float: + return max(0.0, self.prior_metric_variance - self.posterior_metric_variance) + + @property + def reduction_per_cost(self) -> float: + return self.variance_reduction / self.cost + + +class ConditionalGaugeSession: + """Track actual evidence consumption so conditional updates use the right prior. + + The initial prior must be independent of this model's noise and must not have + consumed any model window. Start a new session from that prior and replay the + admitted history; do not wrap an already-updated belief as an untouched prior. + """ + + def __init__(self, model: CorrelatedGaugeDesign, prior: GaussianGaugeBelief) -> None: + if model.chart_id != prior.chart_id: + raise ValueError("model and prior must use the same declared gauge chart") + self.model = model + self._belief = prior + self._history: dict[str, FloatArray] = {} + + @property + def belief(self) -> GaussianGaugeBelief: + return self._belief + + @property + def history_ids(self) -> tuple[str, ...]: + return tuple(self._history) + + def preview_query( + self, + candidate_id: str, + query_jacobian: FloatArray, + *, + query_metric: FloatArray | None = None, + cost: float = 1.0, + ) -> QueryWindowUtility: + query = _array(query_jacobian, name="query_jacobian", ndim=2) + if query.shape[0] < 1 or query.shape[1] != 7: + raise ValueError("query_jacobian must have shape (Q, 7), Q positive") + metric = ( + np.eye(query.shape[0]) + if query_metric is None + else _positive_definite(query_metric, name="query_metric", size=query.shape[0]) + ) + if isinstance(cost, bool) or not np.isfinite(cost) or cost <= 0: + raise ValueError("cost must be finite and positive") + factor = self.model.conditional_factor(self.history_ids, candidate_id) + before = self.belief.covariance + after = _covariance_update(before, factor.whitened_design) + prior_variance = float(np.trace(metric @ query @ before @ query.T)) + posterior_variance = float(np.trace(metric @ query @ after @ query.T)) + if posterior_variance > prior_variance + 1e-10 * max(prior_variance, np.finfo(float).tiny): + raise ValueError("conditional update increased query covariance") + return QueryWindowUtility( + candidate_id, factor.information_rank, prior_variance, posterior_variance, float(cost) + ) + + def assimilate(self, candidate_id: str, value: FloatArray) -> GaussianGaugeBelief: + factor = self.model.conditional_factor(self.history_ids, candidate_id) + candidate = _array(value, name="candidate value", ndim=1) + if candidate.shape != (factor.design.shape[0],): + raise ValueError("candidate value does not match the window size") + history = np.concatenate(tuple(self._history.values())) if self._history else np.empty(0) + whitened_value = factor.whitened_value(history, candidate) + design = factor.whitened_design + prior = self.belief + posterior = prior + if np.any(design): + innovation = np.eye(design.shape[0]) + design @ prior.covariance @ design.T + gain = np.linalg.solve(innovation, design @ prior.covariance).T + mean = prior.mean + gain @ (whitened_value - design @ prior.mean) + posterior = GaussianGaugeBelief( + prior.chart_id, mean, _covariance_update(prior.covariance, design) + ) + self._history[candidate_id] = candidate + self._belief = posterior + return posterior + + +def select_query_window( + utilities: tuple[QueryWindowUtility, ...], *, minimum_gain_per_cost: float = 0.0 +) -> str | None: + """One-step, outcome-free selection; lexical ties; None means no new update. + + This is not a globally optimal multi-window schedule. Correlated observation + utilities need not be submodular; a multi-step approximation guarantee is not + asserted. Thresholds and costs must be frozen before candidate outcomes. + """ + if not np.isfinite(minimum_gain_per_cost) or minimum_gain_per_cost < 0: + raise ValueError("minimum_gain_per_cost must be finite and nonnegative") + if len({item.candidate_id for item in utilities}) != len(utilities): + raise ValueError("candidate IDs must be unique") + if not utilities: + return None + best = min(utilities, key=lambda item: (-item.reduction_per_cost, item.candidate_id)) + return best.candidate_id if best.reduction_per_cost > minimum_gain_per_cost else None diff --git a/src/prob4d/conditional_gauge_study.py b/src/prob4d/conditional_gauge_study.py new file mode 100644 index 000000000..2b4c76f87 --- /dev/null +++ b/src/prob4d/conditional_gauge_study.py @@ -0,0 +1,255 @@ +"""Reproducible, constructed mechanism study; not a real-provider experiment.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +from pathlib import Path +from typing import Any + +import numpy as np + +from .conditional_gauge_design import ( + ConditionalGaugeSession, + CorrelatedGaugeDesign, + GaussianGaugeBelief, + select_query_window, +) + +PROTOCOL: dict[str, Any] = { + "schema": "prob4d.conditional-query-design-study-v1", + "classification": "constructed-linear-Gaussian-mechanism-study", + "seed": 20260830, + "episodes": 10000, + "bootstrap_seed": 8302026, + "bootstrap_replicates": 2000, + "prior_standard_deviation": 0.05, + "history_noise_standard_deviation": 0.01, + "repeat_noise_correlation": 0.999, + "complement_noise_variance": 0.00025, + "global_only_noise_variance": 0.000001, + "cloud_scale_m": 0.1, + "query_point_in_cloud_radii": [0.0, 0.02, 0.01], + "new_window_budget": 1, + "candidate_cost": 1.0, + "ellipsoid_90_chi_square_df3": 6.251388631170325, + "independent_unit": "synthetic episode, not coordinate or point", + "targets_opened": False, + "provider_executed": False, +} + + +def make_design(correlation: float = 0.999) -> CorrelatedGaugeDesign: + """Line-like rank-six history, repeat, complementary support, and twist-only.""" + history = np.eye(7)[[0, 2, 3, 4, 5, 6]] + complement = np.eye(7)[[0, 1, 3, 4, 5, 6]] + global_only = np.eye(7)[[1]] + design = np.vstack((history, history, complement, global_only)) + noise = np.diag([0.0001] * 12 + [0.00025] * 6 + [0.000001]) + noise[:6, 6:12] = correlation * 0.0001 * np.eye(6) + noise[6:12, :6] = correlation * 0.0001 * np.eye(6) + return CorrelatedGaugeDesign( + "centroid-normalized-local-sim3-study-v1", + "known-generative-noise-v1", + ("history", "near_repeat", "complement", "global_only"), + (6, 6, 6, 1), + design, + noise, + ) + + +def query_jacobian() -> np.ndarray: + point = np.array(PROTOCOL["query_point_in_cloud_radii"], dtype=float) + x, y, z = point + skew = np.array([[0.0, -z, y], [z, 0.0, -x], [-y, x, 0.0]]) + return 0.1 * np.column_stack((point, -skew, np.eye(3))) + + +def prior_belief(model: CorrelatedGaugeDesign) -> GaussianGaugeBelief: + return GaussianGaugeBelief(model.chart_id, np.zeros(7), 0.0025 * np.eye(7)) + + +def _history_session(model: CorrelatedGaugeDesign) -> ConditionalGaugeSession: + session = ConditionalGaugeSession(model, prior_belief(model)) + session.assimilate("history", np.zeros(6)) + return session + + +def _vectorized_prediction( + model: CorrelatedGaugeDesign, observations: np.ndarray, selected: str | None +) -> tuple[np.ndarray, np.ndarray]: + """Apply exactly the session's conditional factor to all independent episodes.""" + session = _history_session(model) + old = session.belief + h = model.design_matrix[:6] + prior = prior_belief(model) + k0 = np.linalg.solve( + model.noise_covariance[:6, :6] + h @ prior.covariance @ h.T, + h @ prior.covariance, + ).T + means = observations[:, :6] @ k0.T + if selected is not None: + factor = model.conditional_factor(("history",), selected) + rows = list(model.rows((selected,))) + values = observations[:, rows] - observations[:, :6] @ factor.regression.T + a = factor.whitened_design + k = np.linalg.solve(np.eye(a.shape[0]) + a @ old.covariance @ a.T, a @ old.covariance).T + means += (values @ factor.whitener.T - means @ a.T) @ k.T + session.assimilate(selected, np.zeros(len(rows))) + return means, session.belief.covariance + + +def _paired_bootstrap_interval( + differences: np.ndarray, replicates: int, seed: int +) -> list[float]: + rng = np.random.default_rng(seed) + means = [] + for start in range(0, replicates, 100): + size = min(100, replicates - start) + indices = rng.integers(0, differences.size, size=(size, differences.size)) + means.extend(np.mean(differences[indices], axis=1).tolist()) + return np.quantile(means, [0.025, 0.975]).tolist() + + +def _reference_channel_control() -> dict[str, float]: + design = np.vstack((np.eye(7)[0], np.zeros(7))) + model = CorrelatedGaugeDesign( + "scalar-reference-embedded-in-gauge", "known-reference-noise", + ("signal", "noise_reference"), (1, 1), design, + np.array([[1.0, 1.0], [1.0, 1.1]]), + ) + session = ConditionalGaugeSession( + model, GaussianGaugeBelief(model.chart_id, np.zeros(7), np.eye(7)) + ) + query = np.eye(7)[[0]] + standalone = session.preview_query("noise_reference", query) + session.assimilate("signal", np.zeros(1)) + conditional = session.preview_query("noise_reference", query) + return { + "standalone_reference_gain": standalone.variance_reduction, + "conditional_reference_gain": conditional.variance_reduction, + "variance_before_reference": conditional.prior_metric_variance, + "variance_after_reference": conditional.posterior_metric_variance, + } + + +def run_study(*, episodes: int = 10000, bootstrap_replicates: int = 2000) -> dict[str, Any]: + if episodes < 100 or bootstrap_replicates < 20: + raise ValueError("at least 100 episodes and 20 bootstrap replicates are required") + protocol = dict(PROTOCOL, episodes=episodes, bootstrap_replicates=bootstrap_replicates) + model = make_design() + independent_noise = np.diag(np.diag(model.noise_covariance)) + independence_model = CorrelatedGaugeDesign( + model.chart_id, "deliberately-invalid-independent-noise-control", + model.window_ids, model.window_sizes, model.design_matrix, independent_noise, + ) + query = query_jacobian() + candidates = ("near_repeat", "complement", "global_only") + correct = _history_session(model) + independent = _history_session(independence_model) + correct_utilities = tuple(correct.preview_query(w, query) for w in candidates) + marginal_utilities = tuple(independent.preview_query(w, query) for w in candidates) + global_utilities = tuple(correct.preview_query(w, np.eye(7)) for w in candidates) + query_choice = select_query_window(correct_utilities) + marginal_choice = select_query_window(marginal_utilities) + global_choice = select_query_window(global_utilities) + # Selection is complete before the first observation or latent episode is drawn. + rng = np.random.default_rng(PROTOCOL["seed"]) + truth = 0.05 * rng.normal(size=(episodes, 7)) + noise = rng.normal(size=(episodes, 19)) @ np.linalg.cholesky(model.noise_covariance).T + observations = truth @ model.design_matrix.T + noise + arms: dict[str, tuple[CorrelatedGaugeDesign, str | None]] = { + "history_only": (model, None), + "marginal_query_selection_independent_update": (independence_model, marginal_choice), + "marginal_query_selection_correct_update": (model, marginal_choice), + "global_variance_selection_correct_update": (model, global_choice), + "conditional_query_selection_correct_update": (model, query_choice), + } + errors: dict[str, np.ndarray] = {} + covariances: dict[str, np.ndarray] = {} + parity_maximum = 0.0 + for name, (assumed_model, selected) in arms.items(): + means, covariance = _vectorized_prediction(assumed_model, observations, selected) + selected_ids = ("history",) if selected is None else ("history", selected) + rows = list(assumed_model.rows(selected_ids)) + h = assumed_model.design_matrix[rows] + r = assumed_model.noise_covariance[np.ix_(rows, rows)] + p = prior_belief(assumed_model).covariance + batch_gain = np.linalg.solve(r + h @ p @ h.T, h @ p).T + batch_means = observations[:, rows] @ batch_gain.T + batch_covariance = np.linalg.inv(np.linalg.inv(p) + h.T @ np.linalg.solve(r, h)) + parity_maximum = max(parity_maximum, float(np.max(np.abs(means - batch_means))), + float(np.max(np.abs(covariance - batch_covariance)))) + errors[name] = (means - truth) @ query.T + covariances[name] = query @ covariance @ query.T + baseline_loss = np.sum(errors["history_only"] ** 2, axis=1) + metrics: dict[str, Any] = {} + for name, error in errors.items(): + covariance = covariances[name] + loss = np.sum(error**2, axis=1) + nees = np.sum(error * np.linalg.solve(covariance, error.T).T, axis=1) + sign, logdet = np.linalg.slogdet(covariance) + if sign <= 0: + raise ValueError("query covariance must be positive definite for scoring") + improvement_mm2 = 1e6 * (baseline_loss - loss) + interval = _paired_bootstrap_interval(improvement_mm2, bootstrap_replicates, + PROTOCOL["bootstrap_seed"]) + metrics[name] = { + "selected_window": arms[name][1], + "query_euclidean_rmse_mm": float(1000 * np.sqrt(np.mean(loss))), + "query_expected_rmse_mm": float(1000 * np.sqrt(np.trace(covariance))), + "normalized_query_nees": float(np.mean(nees) / 3), + "query_ellipsoid_90_coverage": float(np.mean(nees <= PROTOCOL["ellipsoid_90_chi_square_df3"])), + "gaussian_query_nll_nats": float(0.5 * (3 * np.log(2 * np.pi) + logdet + np.mean(nees))), + "harmful_episode_fraction_vs_history": float(np.mean(loss > baseline_loss + 1e-18)), + "paired_mean_squared_query_error_improvement_mm2": float(np.mean(improvement_mm2)), + "paired_95_percentile_bootstrap_interval_mm2": interval, + } + correlation_sweep = [] + for correlation in (0.0, 0.5, 0.9, 0.999, 1.0): + sweep_session = _history_session(make_design(correlation)) + utilities = tuple(sweep_session.preview_query(w, query) for w in candidates) + correlation_sweep.append({ + "correlation": correlation, + "selected_window": select_query_window(utilities), + "query_variance_reduction_mm2": {u.candidate_id: 1e6 * u.variance_reduction for u in utilities}, + }) + return { + "schema": "prob4d.conditional-query-design-result-v1", + "protocol": protocol, + "protocol_sha256": hashlib.sha256(json.dumps(protocol, sort_keys=True, separators=(",", ":")).encode()).hexdigest(), + "query_variance_reduction_mm2": { + "correct_conditional": {u.candidate_id: 1e6 * u.variance_reduction for u in correct_utilities}, + "incorrect_marginal": {u.candidate_id: 1e6 * u.variance_reduction for u in marginal_utilities}, + }, + "arms": metrics, + "correlation_sweep": correlation_sweep, + "noise_reference_non_submodularity_control": _reference_channel_control(), + "maximum_kernel_vs_independent_dense_reference_error": parity_maximum, + "boundaries": [ + "Known, correctly specified Gaussian noise except the explicitly invalid control.", + "Constructed local-linear Sim(3) chart, not nonlinear perception or physical simulation.", + "Independent synthetic episodes, not real objects, frames, or target evaluations.", + "Utility guarantees expected squared-query-loss reduction only under the assumed model.", + "Nonzero harmful-episode fractions are retained; no per-update safety claim.", + "No PointWorld, BayesianPhysTwin, or Causal4D runtime or empirical benefit is established.", + ], + } + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument("--episodes", type=int, default=10000) + parser.add_argument("--bootstrap-replicates", type=int, default=2000) + args = parser.parse_args() + result = run_study(episodes=args.episodes, bootstrap_replicates=args.bootstrap_replicates) + args.output.parent.mkdir(parents=True, exist_ok=True) + with args.output.open("x", encoding="utf-8") as stream: + json.dump(result, stream, indent=2, sort_keys=True, allow_nan=False) + stream.write("\n") + + +if __name__ == "__main__": + main() diff --git a/tests/test_conditional_gauge_design.py b/tests/test_conditional_gauge_design.py new file mode 100644 index 000000000..f950ea797 --- /dev/null +++ b/tests/test_conditional_gauge_design.py @@ -0,0 +1,245 @@ +"""Numerical, causal-input, and exact-replay controls for conditional windows.""" + +from __future__ import annotations + +import itertools + +import numpy as np +import pytest + +from prob4d.conditional_gauge_design import ( + ConditionalGaugeSession, + CorrelatedGaugeDesign, + GaussianGaugeBelief, + QueryWindowUtility, + UnsupportedDeterministicConstraint, + select_query_window, +) + + +def _model(design, noise, sizes, ids=None): + return CorrelatedGaugeDesign( + "normalized-test-chart", "known-synthetic-noise", tuple(ids or [f"w{i}" for i in range(len(sizes))]), + tuple(sizes), design, noise, + ) + + +def _prior(mean=None, covariance=None): + return GaussianGaugeBelief( + "normalized-test-chart", np.zeros(7) if mean is None else mean, + np.eye(7) if covariance is None else covariance, + ) + + +@pytest.mark.parametrize("seed", range(5)) +def test_sequential_correlated_updates_equal_dense_joint_posterior(seed): + rng = np.random.default_rng(seed) + design = rng.normal(size=(9, 7)) + root = rng.normal(size=(9, 9)) + noise = root @ root.T + 0.3 * np.eye(9) + p_root = rng.normal(size=(7, 7)) + prior = _prior(rng.normal(size=7), p_root @ p_root.T + np.eye(7)) + values = rng.normal(size=9) + model = _model(design, noise, [2, 3, 4]) + prior_information = np.linalg.solve(prior.covariance, np.eye(7)) + expected_cov = np.linalg.inv(prior_information + design.T @ np.linalg.solve(noise, design)) + expected_mean = expected_cov @ (prior_information @ prior.mean + design.T @ np.linalg.solve(noise, values)) + for order in itertools.permutations(model.window_ids): + session = ConditionalGaugeSession(model, prior) + for window in order: + session.assimilate(window, values[list(model.rows((window,)))]) + np.testing.assert_allclose(session.belief.covariance, expected_cov, rtol=2e-10, atol=2e-11) + np.testing.assert_allclose(session.belief.mean, expected_mean, rtol=2e-10, atol=2e-11) + + +def test_exact_replay_has_zero_utility_and_preserves_complete_belief_identity(): + design = np.eye(7)[[0, 2, 3, 4, 5, 6]] + covariance = 0.1 * np.eye(6) + model = _model(np.vstack((design, design)), np.tile(covariance, (2, 2)), [6, 6]) + session = ConditionalGaugeSession(model, _prior()) + first = session.assimilate("w0", np.arange(6.0)) + utility = session.preview_query("w1", np.eye(7)) + assert utility.conditional_information_rank == 0 + assert utility.variance_reduction == 0 + assert select_query_window((utility,)) is None + replay = session.assimilate("w1", np.arange(6.0)) + assert replay is first + assert session.history_ids == ("w0", "w1") + + +def test_inconsistent_duplicate_fails_before_mutating_session(): + h = np.eye(7)[[0]] + session = ConditionalGaugeSession(_model(np.vstack((h, h)), np.ones((2, 2)), [1, 1]), _prior()) + before = session.assimilate("w0", np.array([1.0])) + with pytest.raises(ValueError, match="noise-support identity"): + session.assimilate("w1", np.array([2.0])) + assert session.belief is before + assert session.history_ids == ("w0",) + + +def test_rank_deficiency_is_preserved_without_fabricated_precision(): + h = np.eye(7)[[0, 2, 3, 4, 5, 6]] + session = ConditionalGaugeSession(_model(h, 0.1 * np.eye(6), [6]), _prior()) + posterior = session.assimilate("w0", np.ones(6)) + assert posterior.mean[1] == 0 + assert posterior.covariance[1, 1] == 1 + assert session.model.conditional_factor((), "w0").information_rank == 6 + + +def test_prior_mediated_gain_is_not_falsely_forbidden(): + p = np.eye(7) + p[0, 1] = p[1, 0] = 0.7 + session = ConditionalGaugeSession(_model(np.eye(7)[[0]], np.eye(1), [1]), _prior(covariance=p)) + assert session.preview_query("w0", np.eye(7)[[1]]).variance_reduction > 0 + assert session.assimilate("w0", np.array([1.0])).covariance[1, 1] < 1 + + +def test_zero_standalone_state_information_can_have_conditional_query_value(): + h = np.vstack((np.eye(7)[0], np.zeros(7))) + model = _model(h, np.array([[1.0, 1.0], [1.0, 1.1]]), [1, 1]) + untouched = ConditionalGaugeSession(model, _prior()) + assert untouched.preview_query("w1", np.eye(7)[[0]]).variance_reduction == 0 + untouched.assimilate("w0", np.array([1.0])) + conditional = untouched.preview_query("w1", np.eye(7)[[0]]) + assert conditional.conditional_information_rank == 1 + assert conditional.variance_reduction > 0.4 + posterior = untouched.assimilate("w1", np.array([0.5])) + assert posterior.covariance[0, 0] < 0.1 + + +def test_candidate_selection_does_not_depend_on_observed_history_values(): + h = np.eye(7)[[0, 0, 1]] + noise = np.array([[1.0, 0.9, 0], [0.9, 1.0, 0], [0, 0, 2.0]]) + model = _model(h, noise, [1, 1, 1]) + utilities = [] + for value in [-1e6, 0.0, 1e6]: + session = ConditionalGaugeSession(model, _prior()) + session.assimilate("w0", np.array([value])) + utilities.append(tuple(session.preview_query(w, np.eye(7)) for w in ["w1", "w2"])) + assert utilities[0] == utilities[1] == utilities[2] + + +def test_global_chart_change_and_query_units_preserve_utility(): + rng = np.random.default_rng(80) + h = rng.normal(size=(5, 7)) + n = rng.normal(size=(5, 5)) + r = n @ n.T + np.eye(5) + j = rng.normal(size=(3, 7)) + transform = rng.normal(size=(7, 7)) + 4 * np.eye(7) + inverse = np.linalg.inv(transform) + original = ConditionalGaugeSession(_model(h, r, [2, 3]), _prior()) + changed = ConditionalGaugeSession(_model(h @ inverse, r, [2, 3]), _prior(covariance=transform @ transform.T)) + original.assimilate("w0", np.ones(2)) + changed.assimilate("w0", np.ones(2)) + a = original.preview_query("w1", j) + b = changed.preview_query("w1", 1000 * j @ inverse, query_metric=1e-6 * np.eye(3)) + np.testing.assert_allclose(a.variance_reduction, b.variance_reduction, rtol=1e-10) + + +@pytest.mark.parametrize("scale", [1e-7, 1.0, 1e7]) +def test_measurement_unit_scaling_preserves_utility(scale): + h = np.eye(7)[[0, 1]] + r = np.array([[1.0, 0.5], [0.5, 1.0]]) + session = ConditionalGaugeSession(_model(scale * h, scale**2 * r, [1, 1]), _prior()) + session.assimilate("w0", np.array([scale])) + reference = ConditionalGaugeSession(_model(h, r, [1, 1]), _prior()) + reference.assimilate("w0", np.ones(1)) + np.testing.assert_allclose(session.preview_query("w1", np.eye(7)).variance_reduction, + reference.preview_query("w1", np.eye(7)).variance_reduction, rtol=1e-10) + + +def test_informative_zero_noise_is_rejected_not_silently_discarded(): + with pytest.raises(UnsupportedDeterministicConstraint): + _model(np.eye(7)[[0]], np.zeros((1, 1)), [1]) + model = _model(np.zeros((1, 7)), np.zeros((1, 1)), [1]) + session = ConditionalGaugeSession(model, _prior()) + prior = session.belief + assert session.assimilate("w0", np.zeros(1)) is prior + + +def test_arrays_are_defensive_copies_and_immutable(): + h = np.eye(7) + r = np.eye(7) + model = _model(h, r, [7]) + h[:] = 3 + r[:] = 4 + np.testing.assert_array_equal(model.design_matrix, np.eye(7)) + with pytest.raises(ValueError): + model.noise_covariance[0, 0] = 7 + prior = _prior() + with pytest.raises(ValueError): + prior.covariance[0, 0] = 7 + + +@pytest.mark.parametrize("bad_noise", [np.array([[1.0, 2.0], [2.0, 1.0]]), np.array([[1.0, 0.2], [0.0, 1.0]]), np.full((2, 2), np.nan)]) +def test_invalid_joint_covariance_is_rejected(bad_noise): + with pytest.raises(ValueError): + _model(np.eye(7)[[0, 1]], bad_noise, [1, 1]) + + +def test_validation_and_no_reconsumption(): + model = _model(np.eye(7)[[0]], np.eye(1), [1]) + with pytest.raises(ValueError, match="same declared"): + ConditionalGaugeSession(model, GaussianGaugeBelief("different", np.zeros(7), np.eye(7))) + session = ConditionalGaugeSession(model, _prior()) + for bad in [np.zeros((1, 1)), np.zeros(2), np.array([np.nan])]: + with pytest.raises(ValueError): + session.assimilate("w0", bad) + assert session.history_ids == () + with pytest.raises(ValueError): + session.preview_query("w0", np.ones((1, 6))) + with pytest.raises(ValueError): + session.preview_query("w0", np.ones((2, 7)), query_metric=np.diag([1.0, -1.0])) + assert session.preview_query("w0", np.zeros((1, 7))).variance_reduction == 0 + session.assimilate("w0", np.zeros(1)) + with pytest.raises(ValueError, match="already been assimilated"): + session.assimilate("w0", np.zeros(1)) + + +def test_selection_threshold_ties_cost_and_empty_cases(): + a = QueryWindowUtility("a", 1, 10.0, 8.0, 2.0) + b = QueryWindowUtility("b", 1, 10.0, 9.0, 1.0) + assert select_query_window((b, a)) == "a" + assert select_query_window((a, b), minimum_gain_per_cost=1) is None + assert select_query_window(()) is None + with pytest.raises(ValueError): + select_query_window((a, a)) + with pytest.raises(ValueError): + select_query_window((a,), minimum_gain_per_cost=-1) + + +@pytest.mark.parametrize("cost", [0, -1, np.inf, np.nan, True]) +def test_invalid_utility_cannot_enter_selector(cost): + with pytest.raises(ValueError): + QueryWindowUtility("bad", 1, 1.0, 0.5, cost) + + +def test_study_reproducibility_and_predeclared_choice(): + from prob4d.conditional_gauge_study import run_study + + first = run_study(episodes=200, bootstrap_replicates=20) + second = run_study(episodes=200, bootstrap_replicates=20) + assert first == second + arms = first["arms"] + assert arms["conditional_query_selection_correct_update"]["selected_window"] == "complement" + assert arms["marginal_query_selection_independent_update"]["selected_window"] == "near_repeat" + assert arms["global_variance_selection_correct_update"]["selected_window"] == "global_only" + assert first["maximum_kernel_vs_independent_dense_reference_error"] < 1e-10 + assert first["correlation_sweep"][-1]["query_variance_reduction_mm2"]["near_repeat"] == 0 + assert first["protocol"]["targets_opened"] is False + + +def test_singular_repeated_history_still_conditions_future_noise_correctly(): + h = np.eye(7)[[0, 0, 1]] + r = np.array([[1.0, 1.0, 0.5], [1.0, 1.0, 0.5], [0.5, 0.5, 1.0]]) + model = _model(h, r, [1, 1, 1]) + a = ConditionalGaugeSession(model, _prior()) + b = ConditionalGaugeSession(model, _prior()) + a.assimilate("w0", np.array([1.0])) + a.assimilate("w1", np.array([1.0])) + a.assimilate("w2", np.array([2.0])) + b.assimilate("w0", np.array([1.0])) + before_replay = b.assimilate("w2", np.array([2.0])) + assert b.assimilate("w1", np.array([1.0])) is before_replay + np.testing.assert_allclose(a.belief.mean, b.belief.mean, atol=1e-12) + np.testing.assert_allclose(a.belief.covariance, b.belief.covariance, atol=1e-12) diff --git a/tests/test_conditional_gauge_integration.py b/tests/test_conditional_gauge_integration.py new file mode 100644 index 000000000..a044bdc6f --- /dev/null +++ b/tests/test_conditional_gauge_integration.py @@ -0,0 +1,51 @@ +"""Parity with the existing observable factor and query projection APIs.""" + +import numpy as np + +from prob4d.conditional_gauge_design import ( + ConditionalGaugeSession, + CorrelatedGaugeDesign, + GaussianGaugeBelief, +) +from prob4d.observable_gauge import CentroidGaugeChart, ObservableGaugeFactor +from prob4d.query_observability import evaluate_query_observability, point_position_query_jacobian +from prob4d.sim3 import Sim3 + + +def test_existing_partial_factor_and_query_api_parity(): + chart = CentroidGaugeChart( + Sim3(scale=1.0, rotation=np.eye(3), translation=np.zeros(3)), np.zeros(3), 0.1 + ) + factor = ObservableGaugeFactor( + chart=chart, + observable_basis=np.eye(7)[:, [0, 2, 3, 4, 5, 6]], + nullspace_basis=np.eye(7)[:, [1]], + observable_information=10 * np.eye(6), + normalized_geometry_spectrum=np.array([1.0] * 6 + [0.0]), + rank_threshold=1e-8, + residual_rms=0.01, + residual_variance=0.0001, + inlier_fraction=1.0, + num_correspondences=16, + covariance_method="iid_observable_information_v1", + ) + model = CorrelatedGaugeDesign( + "one-exact-centroid-chart", "known-observable-factor-covariance", + ("original", "exact-source-replay"), (6, 6), + np.vstack((factor.observable_basis.T, factor.observable_basis.T)), + np.tile(factor.observable_covariance, (2, 2)), + ) + prior = GaussianGaugeBelief(model.chart_id, np.arange(7) / 10, np.eye(7)) + session = ConditionalGaugeSession(model, prior) + query = point_position_query_jacobian(factor, np.array([0.01, 0.02, 0.03])) + existing_report = evaluate_query_observability( + factor, prior_covariance_local=prior.covariance, query_jacobian_local=query + ) + preview = session.preview_query("original", query) + np.testing.assert_allclose(preview.posterior_metric_variance, + np.trace(existing_report.posterior_query_covariance), rtol=1e-12) + expected = factor.fuse_local_gaussian(prior.mean, prior.covariance) + actual = session.assimilate("original", np.zeros(6)) + np.testing.assert_allclose(actual.mean, expected.mean_local, rtol=1e-12) + np.testing.assert_allclose(actual.covariance, expected.covariance_local, rtol=1e-12) + assert session.assimilate("exact-source-replay", np.zeros(6)) is actual From 9743c2e40e0f2ec1c8b4c0be5bf8d240c574cd64 Mon Sep 17 00:00:00 2001 From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com> Date: Sun, 30 Aug 2026 06:53:41 +0800 Subject: [PATCH 2/3] Apply hosted quality fixes and compatible NumPy typing environment --- .../conditional-query-information.yml | 4 +- docs/conditional-query-information.md | 2 +- src/prob4d/conditional_gauge_design.py | 7 +- src/prob4d/conditional_gauge_study.py | 80 +++++++++++++------ tests/test_conditional_gauge_design.py | 35 ++++++-- tests/test_conditional_gauge_integration.py | 13 ++- 6 files changed, 98 insertions(+), 43 deletions(-) diff --git a/.github/workflows/conditional-query-information.yml b/.github/workflows/conditional-query-information.yml index 79c80097d..3ef2842ad 100644 --- a/.github/workflows/conditional-query-information.yml +++ b/.github/workflows/conditional-query-information.yml @@ -37,7 +37,9 @@ jobs: with: python-version: "3.12" - name: Install development environment - run: python -m pip install -e ".[dev]" + run: | + python -m pip install -e ".[dev]" + python -m pip install "numpy>=1.24,<2.3" - name: Check source formatting run: python -m ruff format --diff src/prob4d/conditional_gauge_{design,study}.py tests/test_conditional_gauge_{design,integration}.py - name: Check source lint diff --git a/docs/conditional-query-information.md b/docs/conditional-query-information.md index 52fd2c49a..9db3bd5a7 100644 --- a/docs/conditional-query-information.md +++ b/docs/conditional-query-information.md @@ -151,7 +151,7 @@ python -m prob4d.conditional_gauge_study --output /tmp/conditional-query-result. The frozen numerical design is in `protocols/conditional-query-design-study-v1.json`. The study generates 10,000 independent Gaussian episodes, a rank-six line-like history, a 0.999-correlated -repeat, complementary rank-six support, and a precise but query-irrelevant +repeat, complementary rank-six support, and a precise but weakly query-relevant rotation measurement. Each policy receives one additional window at equal cost. Selections are made before any episode is drawn. The point-position query uses a fixed local Sim(3) Jacobian, not a nonlinear scene renderer or physical solver. diff --git a/src/prob4d/conditional_gauge_design.py b/src/prob4d/conditional_gauge_design.py index 07d0288a1..bf50d2780 100644 --- a/src/prob4d/conditional_gauge_design.py +++ b/src/prob4d/conditional_gauge_design.py @@ -155,8 +155,7 @@ def whitened_value(self, history: FloatArray, candidate: FloatArray) -> FloatArr float(np.sqrt(self.noise_scale)), ) if ( - float(np.linalg.norm(self.history_zero_noise_basis.T @ history)) - > self.rtol * scale + float(np.linalg.norm(self.history_zero_noise_basis.T @ history)) > self.rtol * scale or float(np.linalg.norm(self.zero_noise_basis.T @ residual)) > self.rtol * scale ): raise ValueError("observations violate a deterministic noise-support identity") @@ -196,9 +195,7 @@ def __post_init__(self) -> None: for window_id in ids: _text(window_id, "window ID") if any( - isinstance(size, bool) - or not isinstance(size, (int, np.integer)) - or size < 1 + isinstance(size, bool) or not isinstance(size, (int, np.integer)) or size < 1 for size in sizes ): raise ValueError("window sizes must be positive integers") diff --git a/src/prob4d/conditional_gauge_study.py b/src/prob4d/conditional_gauge_study.py index 2b4c76f87..d32432e4f 100644 --- a/src/prob4d/conditional_gauge_study.py +++ b/src/prob4d/conditional_gauge_study.py @@ -100,9 +100,7 @@ def _vectorized_prediction( return means, session.belief.covariance -def _paired_bootstrap_interval( - differences: np.ndarray, replicates: int, seed: int -) -> list[float]: +def _paired_bootstrap_interval(differences: np.ndarray, replicates: int, seed: int) -> list[float]: rng = np.random.default_rng(seed) means = [] for start in range(0, replicates, 100): @@ -115,8 +113,11 @@ def _paired_bootstrap_interval( def _reference_channel_control() -> dict[str, float]: design = np.vstack((np.eye(7)[0], np.zeros(7))) model = CorrelatedGaugeDesign( - "scalar-reference-embedded-in-gauge", "known-reference-noise", - ("signal", "noise_reference"), (1, 1), design, + "scalar-reference-embedded-in-gauge", + "known-reference-noise", + ("signal", "noise_reference"), + (1, 1), + design, np.array([[1.0, 1.0], [1.0, 1.1]]), ) session = ConditionalGaugeSession( @@ -141,8 +142,12 @@ def run_study(*, episodes: int = 10000, bootstrap_replicates: int = 2000) -> dic model = make_design() independent_noise = np.diag(np.diag(model.noise_covariance)) independence_model = CorrelatedGaugeDesign( - model.chart_id, "deliberately-invalid-independent-noise-control", - model.window_ids, model.window_sizes, model.design_matrix, independent_noise, + model.chart_id, + "deliberately-invalid-independent-noise-control", + model.window_ids, + model.window_sizes, + model.design_matrix, + independent_noise, ) query = query_jacobian() candidates = ("near_repeat", "complement", "global_only") @@ -179,8 +184,11 @@ def run_study(*, episodes: int = 10000, bootstrap_replicates: int = 2000) -> dic batch_gain = np.linalg.solve(r + h @ p @ h.T, h @ p).T batch_means = observations[:, rows] @ batch_gain.T batch_covariance = np.linalg.inv(np.linalg.inv(p) + h.T @ np.linalg.solve(r, h)) - parity_maximum = max(parity_maximum, float(np.max(np.abs(means - batch_means))), - float(np.max(np.abs(covariance - batch_covariance)))) + parity_maximum = max( + parity_maximum, + float(np.max(np.abs(means - batch_means))), + float(np.max(np.abs(covariance - batch_covariance))), + ) errors[name] = (means - truth) @ query.T covariances[name] = query @ covariance @ query.T baseline_loss = np.sum(errors["history_only"] ** 2, axis=1) @@ -193,15 +201,20 @@ def run_study(*, episodes: int = 10000, bootstrap_replicates: int = 2000) -> dic if sign <= 0: raise ValueError("query covariance must be positive definite for scoring") improvement_mm2 = 1e6 * (baseline_loss - loss) - interval = _paired_bootstrap_interval(improvement_mm2, bootstrap_replicates, - PROTOCOL["bootstrap_seed"]) + interval = _paired_bootstrap_interval( + improvement_mm2, bootstrap_replicates, PROTOCOL["bootstrap_seed"] + ) metrics[name] = { "selected_window": arms[name][1], "query_euclidean_rmse_mm": float(1000 * np.sqrt(np.mean(loss))), "query_expected_rmse_mm": float(1000 * np.sqrt(np.trace(covariance))), "normalized_query_nees": float(np.mean(nees) / 3), - "query_ellipsoid_90_coverage": float(np.mean(nees <= PROTOCOL["ellipsoid_90_chi_square_df3"])), - "gaussian_query_nll_nats": float(0.5 * (3 * np.log(2 * np.pi) + logdet + np.mean(nees))), + "query_ellipsoid_90_coverage": float( + np.mean(nees <= PROTOCOL["ellipsoid_90_chi_square_df3"]) + ), + "gaussian_query_nll_nats": float( + 0.5 * (3 * np.log(2 * np.pi) + logdet + np.mean(nees)) + ), "harmful_episode_fraction_vs_history": float(np.mean(loss > baseline_loss + 1e-18)), "paired_mean_squared_query_error_improvement_mm2": float(np.mean(improvement_mm2)), "paired_95_percentile_bootstrap_interval_mm2": interval, @@ -210,18 +223,28 @@ def run_study(*, episodes: int = 10000, bootstrap_replicates: int = 2000) -> dic for correlation in (0.0, 0.5, 0.9, 0.999, 1.0): sweep_session = _history_session(make_design(correlation)) utilities = tuple(sweep_session.preview_query(w, query) for w in candidates) - correlation_sweep.append({ - "correlation": correlation, - "selected_window": select_query_window(utilities), - "query_variance_reduction_mm2": {u.candidate_id: 1e6 * u.variance_reduction for u in utilities}, - }) + correlation_sweep.append( + { + "correlation": correlation, + "selected_window": select_query_window(utilities), + "query_variance_reduction_mm2": { + u.candidate_id: 1e6 * u.variance_reduction for u in utilities + }, + } + ) return { "schema": "prob4d.conditional-query-design-result-v1", "protocol": protocol, - "protocol_sha256": hashlib.sha256(json.dumps(protocol, sort_keys=True, separators=(",", ":")).encode()).hexdigest(), + "protocol_sha256": hashlib.sha256( + json.dumps(protocol, sort_keys=True, separators=(",", ":")).encode() + ).hexdigest(), "query_variance_reduction_mm2": { - "correct_conditional": {u.candidate_id: 1e6 * u.variance_reduction for u in correct_utilities}, - "incorrect_marginal": {u.candidate_id: 1e6 * u.variance_reduction for u in marginal_utilities}, + "correct_conditional": { + u.candidate_id: 1e6 * u.variance_reduction for u in correct_utilities + }, + "incorrect_marginal": { + u.candidate_id: 1e6 * u.variance_reduction for u in marginal_utilities + }, }, "arms": metrics, "correlation_sweep": correlation_sweep, @@ -229,11 +252,20 @@ def run_study(*, episodes: int = 10000, bootstrap_replicates: int = 2000) -> dic "maximum_kernel_vs_independent_dense_reference_error": parity_maximum, "boundaries": [ "Known, correctly specified Gaussian noise except the explicitly invalid control.", - "Constructed local-linear Sim(3) chart, not nonlinear perception or physical simulation.", + ( + "Constructed local-linear Sim(3) chart, " + "not nonlinear perception or physical simulation." + ), "Independent synthetic episodes, not real objects, frames, or target evaluations.", - "Utility guarantees expected squared-query-loss reduction only under the assumed model.", + ( + "Utility guarantees expected squared-query-loss reduction " + "only under the assumed model." + ), "Nonzero harmful-episode fractions are retained; no per-update safety claim.", - "No PointWorld, BayesianPhysTwin, or Causal4D runtime or empirical benefit is established.", + ( + "No PointWorld, BayesianPhysTwin, or Causal4D runtime " + "or empirical benefit is established." + ), ], } diff --git a/tests/test_conditional_gauge_design.py b/tests/test_conditional_gauge_design.py index f950ea797..a3329a396 100644 --- a/tests/test_conditional_gauge_design.py +++ b/tests/test_conditional_gauge_design.py @@ -19,14 +19,19 @@ def _model(design, noise, sizes, ids=None): return CorrelatedGaugeDesign( - "normalized-test-chart", "known-synthetic-noise", tuple(ids or [f"w{i}" for i in range(len(sizes))]), - tuple(sizes), design, noise, + "normalized-test-chart", + "known-synthetic-noise", + tuple(ids or [f"w{i}" for i in range(len(sizes))]), + tuple(sizes), + design, + noise, ) def _prior(mean=None, covariance=None): return GaussianGaugeBelief( - "normalized-test-chart", np.zeros(7) if mean is None else mean, + "normalized-test-chart", + np.zeros(7) if mean is None else mean, np.eye(7) if covariance is None else covariance, ) @@ -43,7 +48,9 @@ def test_sequential_correlated_updates_equal_dense_joint_posterior(seed): model = _model(design, noise, [2, 3, 4]) prior_information = np.linalg.solve(prior.covariance, np.eye(7)) expected_cov = np.linalg.inv(prior_information + design.T @ np.linalg.solve(noise, design)) - expected_mean = expected_cov @ (prior_information @ prior.mean + design.T @ np.linalg.solve(noise, values)) + expected_mean = expected_cov @ ( + prior_information @ prior.mean + design.T @ np.linalg.solve(noise, values) + ) for order in itertools.permutations(model.window_ids): session = ConditionalGaugeSession(model, prior) for window in order: @@ -128,7 +135,9 @@ def test_global_chart_change_and_query_units_preserve_utility(): transform = rng.normal(size=(7, 7)) + 4 * np.eye(7) inverse = np.linalg.inv(transform) original = ConditionalGaugeSession(_model(h, r, [2, 3]), _prior()) - changed = ConditionalGaugeSession(_model(h @ inverse, r, [2, 3]), _prior(covariance=transform @ transform.T)) + changed = ConditionalGaugeSession( + _model(h @ inverse, r, [2, 3]), _prior(covariance=transform @ transform.T) + ) original.assimilate("w0", np.ones(2)) changed.assimilate("w0", np.ones(2)) a = original.preview_query("w1", j) @@ -144,8 +153,11 @@ def test_measurement_unit_scaling_preserves_utility(scale): session.assimilate("w0", np.array([scale])) reference = ConditionalGaugeSession(_model(h, r, [1, 1]), _prior()) reference.assimilate("w0", np.ones(1)) - np.testing.assert_allclose(session.preview_query("w1", np.eye(7)).variance_reduction, - reference.preview_query("w1", np.eye(7)).variance_reduction, rtol=1e-10) + np.testing.assert_allclose( + session.preview_query("w1", np.eye(7)).variance_reduction, + reference.preview_query("w1", np.eye(7)).variance_reduction, + rtol=1e-10, + ) def test_informative_zero_noise_is_rejected_not_silently_discarded(): @@ -171,7 +183,14 @@ def test_arrays_are_defensive_copies_and_immutable(): prior.covariance[0, 0] = 7 -@pytest.mark.parametrize("bad_noise", [np.array([[1.0, 2.0], [2.0, 1.0]]), np.array([[1.0, 0.2], [0.0, 1.0]]), np.full((2, 2), np.nan)]) +@pytest.mark.parametrize( + "bad_noise", + [ + np.array([[1.0, 2.0], [2.0, 1.0]]), + np.array([[1.0, 0.2], [0.0, 1.0]]), + np.full((2, 2), np.nan), + ], +) def test_invalid_joint_covariance_is_rejected(bad_noise): with pytest.raises(ValueError): _model(np.eye(7)[[0, 1]], bad_noise, [1, 1]) diff --git a/tests/test_conditional_gauge_integration.py b/tests/test_conditional_gauge_integration.py index a044bdc6f..746c6696a 100644 --- a/tests/test_conditional_gauge_integration.py +++ b/tests/test_conditional_gauge_integration.py @@ -30,8 +30,10 @@ def test_existing_partial_factor_and_query_api_parity(): covariance_method="iid_observable_information_v1", ) model = CorrelatedGaugeDesign( - "one-exact-centroid-chart", "known-observable-factor-covariance", - ("original", "exact-source-replay"), (6, 6), + "one-exact-centroid-chart", + "known-observable-factor-covariance", + ("original", "exact-source-replay"), + (6, 6), np.vstack((factor.observable_basis.T, factor.observable_basis.T)), np.tile(factor.observable_covariance, (2, 2)), ) @@ -42,8 +44,11 @@ def test_existing_partial_factor_and_query_api_parity(): factor, prior_covariance_local=prior.covariance, query_jacobian_local=query ) preview = session.preview_query("original", query) - np.testing.assert_allclose(preview.posterior_metric_variance, - np.trace(existing_report.posterior_query_covariance), rtol=1e-12) + np.testing.assert_allclose( + preview.posterior_metric_variance, + np.trace(existing_report.posterior_query_covariance), + rtol=1e-12, + ) expected = factor.fuse_local_gaussian(prior.mean, prior.covariance) actual = session.assimilate("original", np.zeros(6)) np.testing.assert_allclose(actual.mean, expected.mean_local, rtol=1e-12) From ee75b49acd7a70ebdc89c066246e189832fdb02c Mon Sep 17 00:00:00 2001 From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com> Date: Sun, 30 Aug 2026 07:01:02 +0800 Subject: [PATCH 3/3] Fix explicit NumPy array alias and scalar return annotations --- src/prob4d/conditional_gauge_design.py | 11 ++++++----- src/prob4d/conditional_gauge_study.py | 5 +++-- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/prob4d/conditional_gauge_design.py b/src/prob4d/conditional_gauge_design.py index bf50d2780..0a8dd1e46 100644 --- a/src/prob4d/conditional_gauge_design.py +++ b/src/prob4d/conditional_gauge_design.py @@ -13,12 +13,12 @@ from __future__ import annotations from dataclasses import dataclass -from typing import Any +from typing import Any, TypeAlias import numpy as np from numpy.typing import NDArray -FloatArray = NDArray[np.floating[Any]] +FloatArray: TypeAlias = NDArray[np.floating[Any]] class UnsupportedDeterministicConstraint(ValueError): @@ -37,7 +37,7 @@ def _symmetric(value: object, *, name: str, size: int) -> FloatArray: result = _array(value, name=name, ndim=2) if result.shape != (size, size): raise ValueError(f"{name} must have shape ({size}, {size})") - scale = max(float(np.max(np.abs(result), initial=0.0)), np.finfo(float).tiny) + scale = max(float(np.max(np.abs(result), initial=0.0)), float(np.finfo(float).tiny)) if float(np.max(np.abs(result - result.T), initial=0.0)) > 1e-10 * scale: raise ValueError(f"{name} must be symmetric") return _array((result + result.T) / 2, name=name, ndim=2) @@ -291,7 +291,7 @@ def __post_init__(self) -> None: raise ValueError(f"{name} must be finite and nonnegative") if isinstance(self.cost, bool) or not np.isfinite(self.cost) or self.cost <= 0: raise ValueError("cost must be finite and positive") - scale = max(self.prior_metric_variance, np.finfo(float).tiny) + scale = max(self.prior_metric_variance, float(np.finfo(float).tiny)) if self.posterior_metric_variance > self.prior_metric_variance + 1e-10 * scale: raise ValueError("posterior query variance must not exceed prior variance") @@ -350,7 +350,8 @@ def preview_query( after = _covariance_update(before, factor.whitened_design) prior_variance = float(np.trace(metric @ query @ before @ query.T)) posterior_variance = float(np.trace(metric @ query @ after @ query.T)) - if posterior_variance > prior_variance + 1e-10 * max(prior_variance, np.finfo(float).tiny): + variance_scale = max(prior_variance, float(np.finfo(float).tiny)) + if posterior_variance > prior_variance + 1e-10 * variance_scale: raise ValueError("conditional update increased query covariance") return QueryWindowUtility( candidate_id, factor.information_rank, prior_variance, posterior_variance, float(cost) diff --git a/src/prob4d/conditional_gauge_study.py b/src/prob4d/conditional_gauge_study.py index d32432e4f..248570a2e 100644 --- a/src/prob4d/conditional_gauge_study.py +++ b/src/prob4d/conditional_gauge_study.py @@ -102,12 +102,13 @@ def _vectorized_prediction( def _paired_bootstrap_interval(differences: np.ndarray, replicates: int, seed: int) -> list[float]: rng = np.random.default_rng(seed) - means = [] + means: list[float] = [] for start in range(0, replicates, 100): size = min(100, replicates - start) indices = rng.integers(0, differences.size, size=(size, differences.size)) means.extend(np.mean(differences[indices], axis=1).tolist()) - return np.quantile(means, [0.025, 0.975]).tolist() + quantiles = np.asarray(np.quantile(means, [0.025, 0.975]), dtype=float) + return [float(quantiles[0]), float(quantiles[1])] def _reference_channel_control() -> dict[str, float]: