diff --git a/CHANGELOG.d/query-conditioned-observability.md b/CHANGELOG.d/query-conditioned-observability.md new file mode 100644 index 000000000..99bba221b --- /dev/null +++ b/CHANGELOG.d/query-conditioned-observability.md @@ -0,0 +1,6 @@ +Add an experimental query-conditioned observability kernel for rank-deficient +Sim(3) gauge factors. The new diagnostics separate direct geometric query +support from complete-prior covariance reduction, expose a source-frozen +three-criterion gate, provide point-position Jacobians, and retain a +deterministic on-support/off-support failure control. No real-provider or +downstream benefit claim is introduced. diff --git a/docs/query-conditioned-observability.md b/docs/query-conditioned-observability.md new file mode 100644 index 000000000..d49769976 --- /dev/null +++ b/docs/query-conditioned-observability.md @@ -0,0 +1,253 @@ +# Query-conditioned observability for partial 4D gauge factors + +Status: **experimental scientific kernel and deterministic mechanism evidence**. +This path is not yet admitted to the claim-bearing provider-v2 export. + +## Why this is needed + +A rank-deficient overlap can be informative without identifying every component +of its relative Sim(3) gauge. The existing observable-subspace factor retains +that valid information and leaves the missing directions to a complete prior. +The remaining deployment question is query-specific: + +> Does the retained information actually constrain the physical quantity that a +> downstream Bayesian twin is about to use? + +A centerline overlap, for example, does not observe twist about the line. It can +still locate points on that line accurately. The same factor may be unsafe for +an off-axis contact point, gripper pose, or counterfactual probe whose position +depends strongly on the missing twist. Pooled reconstruction error does not +express this distinction. + +`prob4d.query_observability` projects one +`ObservableGaugeFactor` through a declared downstream query. It reports direct +geometric support separately from covariance reduction mediated by the prior, +then exposes a source-frozen gate. The caller retains ownership of exact physical +fallback. + +## Formulation + +Let the local gauge coordinate in the factor's centroid-normalized chart be + +\[ +\boldsymbol\zeta = +[\delta\ell,\delta\boldsymbol\phi^\top, + \delta\boldsymbol\tau^\top]^\top \in \mathbb R^7 . +\] + +An observable-subspace factor has information + +\[ +\boldsymbol\Lambda_{\mathrm{obs}} += +\boldsymbol U_r \boldsymbol\Lambda_r \boldsymbol U_r^\top , +\] + +where \(\boldsymbol U_r\) spans the retained observable directions and +\(\boldsymbol N\) spans its orthogonal nullspace. For a downstream query +\(\boldsymbol q(\boldsymbol\zeta)\), freeze its local Jacobian + +\[ +\boldsymbol J_q = +\left.\frac{\partial\boldsymbol q}{\partial\boldsymbol\zeta}\right|_{\boldsymbol 0} +\] + +and a positive-definite output metric \(\boldsymbol M_q\). The metric is +important for mixed-unit or differently weighted query outputs; omitting it +selects the identity metric explicitly. + +### Direct query observability + +The direct fraction is + +\[ +s_{\mathrm{dir}} = +\frac{\lVert\boldsymbol M_q^{1/2}\boldsymbol J_q\boldsymbol U_r\rVert_F^2} +{\lVert\boldsymbol M_q^{1/2}\boldsymbol J_q\boldsymbol U_r\rVert_F^2+ + \lVert\boldsymbol M_q^{1/2}\boldsymbol J_q\boldsymbol N\rVert_F^2}. +\] + +This is a geometric statement. It does not increase merely because the prior +correlates an unobserved direction with an observed one. The complementary +nullspace-sensitivity fraction is \(1-s_{\mathrm{dir}}\). + +### Complete-prior covariance reduction + +For a full-rank local prior +\(\boldsymbol P^-\), the posterior covariance is + +\[ +\boldsymbol P^+ = +[(\boldsymbol P^-)^{-1}+\boldsymbol\Lambda_{\mathrm{obs}}]^{-1}. +\] + +The corresponding query covariances are + +\[ +\boldsymbol\Sigma_q^- = +\boldsymbol J_q\boldsymbol P^-\boldsymbol J_q^\top,\qquad +\boldsymbol\Sigma_q^+ = +\boldsymbol J_q\boldsymbol P^+\boldsymbol J_q^\top. +\] + +The metric-weighted variance reduction is + +\[ +s_{\mathrm{red}} = +1- +\frac{\operatorname{tr}(\boldsymbol M_q\boldsymbol\Sigma_q^+)} + {\operatorname{tr}(\boldsymbol M_q\boldsymbol\Sigma_q^-)}. +\] + +This quantity may be positive even when \(s_{\mathrm{dir}}=0\), because a +correlated complete prior can transmit information between gauge directions. +The implementation deliberately reports both quantities instead of calling +prior-mediated reduction direct visual observability. + +The worst supported variance ratio is the largest eigenvalue of + +\[ +(\boldsymbol\Sigma_q^-)^{-1/2} +\boldsymbol\Sigma_q^+ +(\boldsymbol\Sigma_q^-)^{-1/2} +\] + +on the positive-variance support of \(\boldsymbol\Sigma_q^-\). It catches a +single unresolved query direction that can be hidden by a favorable average +trace reduction. + +## Query gate + +`QueryObservabilityGate` applies three independently frozen thresholds: + +1. minimum direct observability fraction; +2. minimum metric-weighted variance reduction; and +3. maximum worst-direction variance ratio. + +A rejection returns stable reason codes. It does **not** mutate the candidate, +replace covariance, invent a ridge, or implement fallback itself. BayesianPhysTwin +must either select a complete admitted candidate or return its exact caller-owned +physical belief. + +Thresholds, the query Jacobian, and the query metric must be selected on +development/source or calibration groups before any protected target outcome is +opened. A target-side threshold adjustment would be a new experiment. + +## Point-position helper + +For a source point \(\boldsymbol x\), the transformed-point query has Jacobian + +\[ +\boldsymbol J_x = +\begin{bmatrix} +\boldsymbol q & -[\boldsymbol q]_\times & \rho\boldsymbol I_3 +\end{bmatrix}, +\] + +where \(\boldsymbol q\) is the fitted point relative to the transformed centroid +and \(\rho\) is the chart cloud scale. Use: + +```python +from prob4d.query_observability import ( + QueryObservabilityGate, + evaluate_query_observability, + point_position_query_jacobian, +) + +jacobian = point_position_query_jacobian(factor, source_point) +report = evaluate_query_observability( + factor, + prior_covariance_local=complete_prior_covariance, + query_jacobian_local=jacobian, + query_metric=point_metric, +) +decision = frozen_gate.evaluate(report) +``` + +Arbitrary nonlinear physical queries remain supported by supplying their own +local Jacobian. + +## Deterministic control + +The checked-in control isolates one missing rotation-about-line direction. It +uses the same complete identity prior and one frozen gate for three cases. + +| Case | Rank | Direct fraction | Variance reduction | Worst ratio | Decision | +|---|---:|---:|---:|---:|---| +| Point on observed line | 6 | 1.000 | 0.909 | 0.091 | admit | +| Distant off-axis probe | 6 | 0.679 | 0.618 | 0.965 | reject | +| Invalid full-rank completion, same probe | 7 | 1.000 | 0.909 | 0.091 | admit | + +The invalid completion is intentionally included as a failure control. It +demonstrates that query conditioning cannot repair a factor that already +fabricated information in the geometric nullspace. Preserving rank deficiency +is therefore a prerequisite for a meaningful downstream gate. + +Reproduce the result with: + +```bash +PYTHONPATH=src python -m prob4d.query_observability_study \ + --output evidence/query-observability-control-v1/result.json +``` + +## Paper contribution enabled by this kernel + +The defensible Prob4D companion-paper claim is not that geometric degeneracy was +newly discovered. Localizability-aware registration already studies weak pose +directions. The contribution is the complete chain: + +1. retain a rank-deficient **Sim(3) likelihood** from overlapping learned 4D + predictions rather than reject it or numerically complete it; +2. fuse it with a complete correlated prior without pretending that the visual + provider measured its nullspace; +3. project the resulting belief into a declared action-conditioned physical + query; +4. distinguish direct geometric support from prior-mediated uncertainty + reduction; +5. admit or reject the complete candidate prospectively; and +6. expose exact fallback, harmful accepted updates, proper scores, and + worst-group regret as primary endpoints. + +This moves the scientific object from “better aligned 4D reconstruction” to +“which physical decisions are identified by an uncertain 4D observation?” + +## Required real-provider promotion + +The highest-value promotion path remains PointWorld on a fresh, +garment-disjoint Flat'n'Fold robot cohort after issue #333 passes every +source-support gate. + +Before target access, freeze: + +- exact PointWorld checkpoint, runtime, normalization statistics, and input + lineage; +- persistent sparse point identities and cross-window association; +- camera/action/metric-frame support; +- source-only covariance and reliability calibration; +- the physical query Jacobian and output metric; +- query-observability thresholds; +- the BayesianPhysTwin regret guard and exact fallback; and +- garment/session statistical units and one-shot target scoring. + +The held-out comparison should include: + +1. unchanged physical fallback; +2. direct PointWorld or strongest simple deterministic comparator; +3. full-rank/ridge alignment failure control; +4. observable-subspace Prob4D candidate without query admission; +5. query-conditioned Prob4D candidate; and +6. guarded BayesianPhysTwin deployment with exact fallback. + +Provider competence and downstream physical-query value must be reported +separately. Primary endpoints are object/session-clustered proper score, physical +query error, accepted-update harm, worst-group regret, calibration and width, +admission/fallback rate, and off-support query behavior. Causal4D may consume +only the admitted belief and remains an optional downstream experiment. + +## Claim boundary + +The current implementation and deterministic control establish mechanism and +software evidence only. They do not establish real-provider competence, +held-out calibration, BayesianPhysTwin benefit, Causal4D benefit, deployment +safety, or state of the art. A material paper claim requires the prospective +fresh-provider experiment above. diff --git a/evidence/query-observability-control-v1/result.json b/evidence/query-observability-control-v1/result.json new file mode 100644 index 000000000..85bb4b92f --- /dev/null +++ b/evidence/query-observability-control-v1/result.json @@ -0,0 +1,69 @@ +{ + "chart_coordinates": "[log-scale, left-rotation(3), centroid-translation/cloud-scale]", + "claim_boundary": "Deterministic mechanism evidence only; no real-provider accuracy, calibration, BayesianPhysTwin benefit, Causal4D benefit, or deployment claim.", + "controlled_missing_direction": "rotation-x", + "factor_information_precision": 10.0, + "gate": { + "maximum_worst_supported_variance_ratio": 0.5, + "minimum_direct_observability_fraction": 0.8, + "minimum_metric_variance_reduction_fraction": 0.8 + }, + "prior_covariance": "identity-7", + "queries": { + "off_support_source_point": [ + 0.0, + 5.0, + 0.0 + ], + "output_metric": "identity-3", + "supported_source_point": [ + 1.0, + 0.0, + 0.0 + ] + }, + "results": { + "invalid_full_rank_off_support_query": { + "admitted": true, + "direct_observability_fraction": 1.0, + "factor_rank": 7, + "metric_variance_reduction_fraction": 0.9090909090909091, + "nullspace_sensitivity_fraction": 0.0, + "posterior_metric_variance": 7.090909090909092, + "prior_metric_variance": 78.0, + "query_dimension": 3, + "reason_codes": [], + "worst_supported_variance_ratio": 0.09090909090909093 + }, + "rank_six_off_support_query": { + "admitted": false, + "direct_observability_fraction": 0.6794871794871795, + "factor_rank": 6, + "metric_variance_reduction_fraction": 0.6177156177156178, + "nullspace_sensitivity_fraction": 0.32051282051282054, + "posterior_metric_variance": 29.818181818181817, + "prior_metric_variance": 78.0, + "query_dimension": 3, + "reason_codes": [ + "insufficient-direct-query-observability", + "insufficient-query-variance-reduction", + "excessive-worst-direction-variance-ratio" + ], + "worst_supported_variance_ratio": 0.9650349650349651 + }, + "rank_six_supported_query": { + "admitted": true, + "direct_observability_fraction": 1.0, + "factor_rank": 6, + "metric_variance_reduction_fraction": 0.9090909090909092, + "nullspace_sensitivity_fraction": 0.0, + "posterior_metric_variance": 0.5454545454545454, + "prior_metric_variance": 6.0, + "query_dimension": 3, + "reason_codes": [], + "worst_supported_variance_ratio": 0.09090909090909088 + } + }, + "schema_name": "prob4d.query-conditioned-observability-control", + "schema_version": 1 +} diff --git a/src/prob4d/query_observability.py b/src/prob4d/query_observability.py new file mode 100644 index 000000000..2b2a7ad02 --- /dev/null +++ b/src/prob4d/query_observability.py @@ -0,0 +1,431 @@ +"""Decision-relevant observability diagnostics for partial Sim(3) gauge factors. + +A rank-deficient visual factor can be useful for one physical query and +uninformative for another. This module projects an +:class:`~prob4d.observable_gauge.ObservableGaugeFactor` through a user-supplied +query Jacobian in the factor's centroid-normalized local chart. It reports +direct geometric support separately from prior-mediated posterior variance +reduction and provides a threshold gate whose values must be frozen on source +or calibration groups. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any, TypeAlias + +import numpy as np +from numpy.typing import NDArray + +from .observable_gauge import ObservableGaugeFactor +from .sim3 import skew + +FloatArray: TypeAlias = NDArray[np.floating[Any]] + +_DIRECT_SUPPORT_REASON = "insufficient-direct-query-observability" +_VARIANCE_REDUCTION_REASON = "insufficient-query-variance-reduction" +_WORST_DIRECTION_REASON = "excessive-worst-direction-variance-ratio" + + +def _readonly_matrix( + value: object, + *, + name: str, + shape: tuple[int, int] | None = None, + columns: int | None = None, +) -> FloatArray: + array = np.asarray(value, dtype=np.float64).copy() + if array.ndim != 2: + raise ValueError(f"{name} must be a matrix") + if shape is not None and array.shape != shape: + raise ValueError(f"{name} must have shape {shape}") + if columns is not None and (array.shape[0] < 1 or array.shape[1] != columns): + raise ValueError(f"{name} must have shape (Q, {columns}) with Q positive") + if not np.all(np.isfinite(array)): + raise ValueError(f"{name} must be finite") + array.setflags(write=False) + return array + + +def _symmetric_positive_definite_square( + value: object, + *, + name: str, +) -> FloatArray: + array = np.asarray(value, dtype=np.float64).copy() + if array.ndim != 2 or array.shape[0] < 1 or array.shape[0] != array.shape[1]: + raise ValueError(f"{name} must be a nonempty square matrix") + if not np.all(np.isfinite(array)): + raise ValueError(f"{name} must be finite") + symmetric = 0.5 * (array + array.T) + if not np.allclose(array, symmetric, atol=1e-12, rtol=1e-10): + raise ValueError(f"{name} must be symmetric") + if float(np.min(np.linalg.eigvalsh(symmetric))) <= 0.0: + raise ValueError(f"{name} must be positive definite") + symmetric.setflags(write=False) + return symmetric + + +def _symmetric_positive_definite(value: object, *, name: str) -> FloatArray: + matrix = _readonly_matrix(value, name=name, shape=(7, 7)) + return _symmetric_positive_definite_square(matrix, name=name) + + +def _symmetric_positive_semidefinite( + value: FloatArray, + *, + name: str, +) -> FloatArray: + array = np.asarray(value, dtype=np.float64) + symmetric = 0.5 * (array + array.T) + scale = max(float(np.max(np.abs(symmetric), initial=0.0)), 1.0) + if float(np.min(np.linalg.eigvalsh(symmetric))) < -1e-10 * scale: + raise ValueError(f"{name} must be positive semidefinite") + symmetric = symmetric.copy() + symmetric.setflags(write=False) + return symmetric + + +def _unit_interval(value: object, *, name: str) -> float: + numeric = float(value) + if not np.isfinite(numeric) or not 0.0 <= numeric <= 1.0: + raise ValueError(f"{name} must lie in [0, 1]") + return numeric + + +def _supported_variance_ratio( + prior_covariance: FloatArray, + posterior_covariance: FloatArray, +) -> float: + eigenvalues, eigenvectors = np.linalg.eigh(prior_covariance) + maximum = float(np.max(eigenvalues, initial=0.0)) + threshold = max(maximum * 1e-10, np.finfo(np.float64).eps) + supported = eigenvalues > threshold + if not np.any(supported): + return 0.0 + whitener = ( + eigenvectors[:, supported] / np.sqrt(eigenvalues[supported]) + ).T + relative = whitener @ posterior_covariance @ whitener.T + relative = 0.5 * (relative + relative.T) + maximum_ratio = float(np.max(np.linalg.eigvalsh(relative), initial=0.0)) + if maximum_ratio > 1.0 + 1e-8: + raise ValueError("posterior query covariance exceeds the prior covariance") + return float(np.clip(maximum_ratio, 0.0, 1.0)) + + +@dataclass(frozen=True) +class QueryObservabilityReport: + """Projection of one partial gauge factor into a downstream query.""" + + factor_rank: int + query_dimension: int + direct_observability_fraction: float + nullspace_sensitivity_fraction: float + query_metric: FloatArray + prior_query_covariance: FloatArray + posterior_query_covariance: FloatArray + metric_variance_reduction_fraction: float + worst_supported_variance_ratio: float + gauge_invariant_query: bool + + def __post_init__(self) -> None: + if isinstance(self.factor_rank, bool) or not isinstance( + self.factor_rank, (int, np.integer) + ): + raise TypeError("factor_rank must be an integer") + factor_rank = int(self.factor_rank) + if not 1 <= factor_rank <= 7: + raise ValueError("factor_rank must lie in [1, 7]") + if isinstance(self.query_dimension, bool) or not isinstance( + self.query_dimension, (int, np.integer) + ): + raise TypeError("query_dimension must be an integer") + query_dimension = int(self.query_dimension) + if query_dimension < 1: + raise ValueError("query_dimension must be positive") + direct = _unit_interval( + self.direct_observability_fraction, + name="direct_observability_fraction", + ) + nullspace = _unit_interval( + self.nullspace_sensitivity_fraction, + name="nullspace_sensitivity_fraction", + ) + if not np.isclose(direct + nullspace, 1.0, atol=1e-10, rtol=1e-10): + raise ValueError( + "direct and nullspace sensitivity fractions must sum to one" + ) + query_metric = _readonly_matrix( + self.query_metric, + name="query_metric", + shape=(query_dimension, query_dimension), + ) + query_metric = _symmetric_positive_definite_square( + query_metric, + name="query_metric", + ) + prior = _readonly_matrix( + self.prior_query_covariance, + name="prior_query_covariance", + shape=(query_dimension, query_dimension), + ) + prior = _symmetric_positive_semidefinite( + prior, + name="prior_query_covariance", + ) + posterior = _readonly_matrix( + self.posterior_query_covariance, + name="posterior_query_covariance", + shape=(query_dimension, query_dimension), + ) + posterior = _symmetric_positive_semidefinite( + posterior, + name="posterior_query_covariance", + ) + reduction = _unit_interval( + self.metric_variance_reduction_fraction, + name="metric_variance_reduction_fraction", + ) + worst_ratio = _unit_interval( + self.worst_supported_variance_ratio, + name="worst_supported_variance_ratio", + ) + if type(self.gauge_invariant_query) is not bool: + raise TypeError("gauge_invariant_query must be a bool") + object.__setattr__(self, "factor_rank", factor_rank) + object.__setattr__(self, "query_dimension", query_dimension) + object.__setattr__(self, "direct_observability_fraction", direct) + object.__setattr__(self, "nullspace_sensitivity_fraction", nullspace) + object.__setattr__(self, "query_metric", query_metric) + object.__setattr__(self, "prior_query_covariance", prior) + object.__setattr__(self, "posterior_query_covariance", posterior) + object.__setattr__( + self, + "metric_variance_reduction_fraction", + reduction, + ) + object.__setattr__( + self, + "worst_supported_variance_ratio", + worst_ratio, + ) + + @property + def prior_metric_variance(self) -> float: + return float(np.trace(self.query_metric @ self.prior_query_covariance)) + + @property + def posterior_metric_variance(self) -> float: + return float(np.trace(self.query_metric @ self.posterior_query_covariance)) + + +@dataclass(frozen=True) +class QueryObservabilityGate: + """Source-frozen thresholds for admitting a query-specific factor.""" + + minimum_direct_observability_fraction: float = 0.0 + minimum_metric_variance_reduction_fraction: float = 0.0 + maximum_worst_supported_variance_ratio: float = 1.0 + + def __post_init__(self) -> None: + object.__setattr__( + self, + "minimum_direct_observability_fraction", + _unit_interval( + self.minimum_direct_observability_fraction, + name="minimum_direct_observability_fraction", + ), + ) + object.__setattr__( + self, + "minimum_metric_variance_reduction_fraction", + _unit_interval( + self.minimum_metric_variance_reduction_fraction, + name="minimum_metric_variance_reduction_fraction", + ), + ) + object.__setattr__( + self, + "maximum_worst_supported_variance_ratio", + _unit_interval( + self.maximum_worst_supported_variance_ratio, + name="maximum_worst_supported_variance_ratio", + ), + ) + + def evaluate( + self, + report: QueryObservabilityReport, + ) -> QueryObservabilityDecision: + """Return a deterministic admission decision for one report.""" + + reasons: list[str] = [] + if ( + report.direct_observability_fraction + < self.minimum_direct_observability_fraction + ): + reasons.append(_DIRECT_SUPPORT_REASON) + if ( + report.metric_variance_reduction_fraction + < self.minimum_metric_variance_reduction_fraction + ): + reasons.append(_VARIANCE_REDUCTION_REASON) + if ( + report.worst_supported_variance_ratio + > self.maximum_worst_supported_variance_ratio + ): + reasons.append(_WORST_DIRECTION_REASON) + return QueryObservabilityDecision( + admitted=not reasons, + reason_codes=tuple(reasons), + ) + + +@dataclass(frozen=True) +class QueryObservabilityDecision: + """Admission metadata; the caller owns exact fallback semantics.""" + + admitted: bool + reason_codes: tuple[str, ...] + + def __post_init__(self) -> None: + if type(self.admitted) is not bool: + raise TypeError("admitted must be a bool") + reason_codes = tuple(str(reason) for reason in self.reason_codes) + if any(not reason for reason in reason_codes): + raise ValueError("reason_codes must contain only nonempty strings") + if len(set(reason_codes)) != len(reason_codes): + raise ValueError("reason_codes must not contain duplicates") + if self.admitted and reason_codes: + raise ValueError("an admitted decision cannot contain rejection reasons") + if not self.admitted and not reason_codes: + raise ValueError("a rejected decision must contain at least one reason") + object.__setattr__(self, "reason_codes", reason_codes) + + +def point_position_query_jacobian( + factor: ObservableGaugeFactor, + source_point: FloatArray, +) -> FloatArray: + """Linearize one gauge-transformed source-point position. + + The returned ``3 x 7`` Jacobian uses the factor's intrinsic coordinates. + It is suitable for endpoint, marker, or probe-position queries whose point + is expressed in the moving/source window coordinates. + """ + + point = np.asarray(source_point, dtype=np.float64).copy() + if point.shape != (3,): + raise ValueError("source_point must have shape (3,)") + if not np.all(np.isfinite(point)): + raise ValueError("source_point must be finite") + transformed = factor.chart.linearization.transform_points(point) + centered = transformed - factor.chart.reference_centroid + jacobian = np.empty((3, 7), dtype=np.float64) + jacobian[:, 0] = centered + jacobian[:, 1:4] = -skew(centered) + jacobian[:, 4:7] = factor.chart.cloud_scale * np.eye(3) + jacobian.setflags(write=False) + return jacobian + + +def evaluate_query_observability( + factor: ObservableGaugeFactor, + *, + prior_covariance_local: FloatArray, + query_jacobian_local: FloatArray, + query_metric: FloatArray | None = None, +) -> QueryObservabilityReport: + """Project factor information and a complete prior into one local query. + + ``query_jacobian_local`` linearizes a downstream query with respect to the + factor chart coordinates + ``[log scale, left rotation(3), centroid translation / cloud scale]``. + Direct observability depends only on the query, its declared output metric, + and the factor subspaces. Variance reduction additionally depends on the + supplied complete prior and therefore remains explicitly separate. The + identity metric is used only when ``query_metric`` is omitted. + """ + + prior = _symmetric_positive_definite( + prior_covariance_local, + name="prior_covariance_local", + ) + query_jacobian = _readonly_matrix( + query_jacobian_local, + name="query_jacobian_local", + columns=7, + ) + metric = ( + np.eye(query_jacobian.shape[0], dtype=np.float64) + if query_metric is None + else _symmetric_positive_definite_square( + query_metric, + name="query_metric", + ) + ) + prior_information = np.linalg.solve(prior, np.eye(7)) + posterior_information = prior_information + factor.information_matrix + posterior = np.linalg.solve(posterior_information, np.eye(7)) + posterior = 0.5 * (posterior + posterior.T) + + prior_query = _symmetric_positive_semidefinite( + query_jacobian @ prior @ query_jacobian.T, + name="prior_query_covariance", + ) + posterior_query = _symmetric_positive_semidefinite( + query_jacobian @ posterior @ query_jacobian.T, + name="posterior_query_covariance", + ) + + metric_sqrt = np.linalg.cholesky(metric).T + metric_jacobian = metric_sqrt @ query_jacobian + observable_energy = float( + np.sum((metric_jacobian @ factor.observable_basis) ** 2) + ) + nullspace_energy = float( + np.sum((metric_jacobian @ factor.nullspace_basis) ** 2) + ) + total_energy = observable_energy + nullspace_energy + gauge_invariant = bool(total_energy <= np.finfo(np.float64).eps) + if gauge_invariant: + direct_fraction = 1.0 + nullspace_fraction = 0.0 + else: + direct_fraction = observable_energy / total_energy + nullspace_fraction = nullspace_energy / total_energy + + prior_trace = float(np.trace(metric @ prior_query)) + posterior_trace = float(np.trace(metric @ posterior_query)) + if prior_trace <= np.finfo(np.float64).eps: + trace_reduction = 0.0 + else: + trace_reduction = float( + np.clip((prior_trace - posterior_trace) / prior_trace, 0.0, 1.0) + ) + + return QueryObservabilityReport( + factor_rank=factor.rank, + query_dimension=int(query_jacobian.shape[0]), + direct_observability_fraction=direct_fraction, + nullspace_sensitivity_fraction=nullspace_fraction, + query_metric=metric, + prior_query_covariance=prior_query, + posterior_query_covariance=posterior_query, + metric_variance_reduction_fraction=trace_reduction, + worst_supported_variance_ratio=_supported_variance_ratio( + prior_query, + posterior_query, + ), + gauge_invariant_query=gauge_invariant, + ) + + +__all__ = [ + "QueryObservabilityDecision", + "QueryObservabilityGate", + "QueryObservabilityReport", + "evaluate_query_observability", + "point_position_query_jacobian", +] diff --git a/src/prob4d/query_observability_study.py b/src/prob4d/query_observability_study.py new file mode 100644 index 000000000..c419b70e5 --- /dev/null +++ b/src/prob4d/query_observability_study.py @@ -0,0 +1,200 @@ +"""Analytic control for query-conditioned observability. + +The control isolates one unobservable twist direction. A rank-six factor +correctly supports a point on the observed line and rejects a distant off-axis +probe under a frozen query gate. An intentionally invalid full-rank completion +passes the same off-axis query because it fabricates information in the missing +direction. +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path +from typing import Any + +import numpy as np + +from .observable_gauge import ( + IID_OBSERVABLE_INFORMATION, + CentroidGaugeChart, + ObservableGaugeFactor, +) +from .query_observability import ( + QueryObservabilityDecision, + QueryObservabilityGate, + QueryObservabilityReport, + evaluate_query_observability, + point_position_query_jacobian, +) +from .sim3 import Sim3 + +SCHEMA_NAME = "prob4d.query-conditioned-observability-control" +SCHEMA_VERSION = 1 + + +def _controlled_factor(*, complete_nullspace: bool) -> ObservableGaugeFactor: + identity = np.eye(7) + if complete_nullspace: + observable_basis = identity + nullspace_basis = np.empty((7, 0)) + observable_information = 10.0 * identity + spectrum = np.ones(7) + else: + observable_basis = np.delete(identity, 1, axis=1) + nullspace_basis = identity[:, 1:2] + observable_information = 10.0 * np.eye(6) + spectrum = np.concatenate((np.ones(6), np.zeros(1))) + return ObservableGaugeFactor( + chart=CentroidGaugeChart( + linearization=Sim3.identity(), + source_centroid=np.zeros(3), + cloud_scale=1.0, + ), + observable_basis=observable_basis, + nullspace_basis=nullspace_basis, + observable_information=observable_information, + normalized_geometry_spectrum=spectrum, + rank_threshold=1e-8, + residual_rms=0.01, + residual_variance=0.01, + inlier_fraction=1.0, + num_correspondences=48, + covariance_method=IID_OBSERVABLE_INFORMATION, + ) + + +def _report_payload( + report: QueryObservabilityReport, + decision: QueryObservabilityDecision, +) -> dict[str, Any]: + return { + "factor_rank": report.factor_rank, + "query_dimension": report.query_dimension, + "direct_observability_fraction": report.direct_observability_fraction, + "nullspace_sensitivity_fraction": report.nullspace_sensitivity_fraction, + "metric_variance_reduction_fraction": ( + report.metric_variance_reduction_fraction + ), + "worst_supported_variance_ratio": ( + report.worst_supported_variance_ratio + ), + "prior_metric_variance": report.prior_metric_variance, + "posterior_metric_variance": report.posterior_metric_variance, + "admitted": decision.admitted, + "reason_codes": list(decision.reason_codes), + } + + +def run_query_observability_study() -> dict[str, Any]: + """Run the deterministic partial-observability control.""" + + partial = _controlled_factor(complete_nullspace=False) + invalid_completion = _controlled_factor(complete_nullspace=True) + prior_covariance = np.eye(7) + supported_point = np.array([1.0, 0.0, 0.0]) + off_support_point = np.array([0.0, 5.0, 0.0]) + gate = QueryObservabilityGate( + minimum_direct_observability_fraction=0.80, + minimum_metric_variance_reduction_fraction=0.80, + maximum_worst_supported_variance_ratio=0.50, + ) + + partial_supported = evaluate_query_observability( + partial, + prior_covariance_local=prior_covariance, + query_jacobian_local=point_position_query_jacobian( + partial, + supported_point, + ), + ) + partial_off_support = evaluate_query_observability( + partial, + prior_covariance_local=prior_covariance, + query_jacobian_local=point_position_query_jacobian( + partial, + off_support_point, + ), + ) + completion_off_support = evaluate_query_observability( + invalid_completion, + prior_covariance_local=prior_covariance, + query_jacobian_local=point_position_query_jacobian( + invalid_completion, + off_support_point, + ), + ) + + return { + "schema_name": SCHEMA_NAME, + "schema_version": SCHEMA_VERSION, + "claim_boundary": ( + "Deterministic mechanism evidence only; no real-provider accuracy, " + "calibration, BayesianPhysTwin benefit, Causal4D benefit, or " + "deployment claim." + ), + "chart_coordinates": ( + "[log-scale, left-rotation(3), centroid-translation/cloud-scale]" + ), + "controlled_missing_direction": "rotation-x", + "factor_information_precision": 10.0, + "prior_covariance": "identity-7", + "queries": { + "supported_source_point": supported_point.tolist(), + "off_support_source_point": off_support_point.tolist(), + "output_metric": "identity-3", + }, + "gate": { + "minimum_direct_observability_fraction": ( + gate.minimum_direct_observability_fraction + ), + "minimum_metric_variance_reduction_fraction": ( + gate.minimum_metric_variance_reduction_fraction + ), + "maximum_worst_supported_variance_ratio": ( + gate.maximum_worst_supported_variance_ratio + ), + }, + "results": { + "rank_six_supported_query": _report_payload( + partial_supported, + gate.evaluate(partial_supported), + ), + "rank_six_off_support_query": _report_payload( + partial_off_support, + gate.evaluate(partial_off_support), + ), + "invalid_full_rank_off_support_query": _report_payload( + completion_off_support, + gate.evaluate(completion_off_support), + ), + }, + } + + +def _parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Run the query-conditioned observability analytic control." + ) + parser.add_argument( + "--output", + type=Path, + help="Optional JSON output path.", + ) + return parser.parse_args() + + +def main() -> None: + args = _parse_args() + result = run_query_observability_study() + serialized = json.dumps(result, indent=2, sort_keys=True) + "\n" + if args.output is None: + print(serialized, end="") + return + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(serialized, encoding="utf-8") + + +if __name__ == "__main__": + main() diff --git a/tests/test_query_observability.py b/tests/test_query_observability.py new file mode 100644 index 000000000..5c17eaed3 --- /dev/null +++ b/tests/test_query_observability.py @@ -0,0 +1,258 @@ +from __future__ import annotations + +import json +from dataclasses import dataclass +from pathlib import Path +from typing import cast + +import numpy as np +import pytest + +from prob4d.observable_gauge import ObservableGaugeFactor +from prob4d.query_observability import ( + QueryObservabilityGate, + evaluate_query_observability, + point_position_query_jacobian, +) +from prob4d.query_observability_study import run_query_observability_study +from prob4d.sim3 import Sim3 + + +@dataclass(frozen=True) +class _Chart: + linearization: Sim3 + reference_centroid: np.ndarray + cloud_scale: float + + +@dataclass(frozen=True) +class _PartialFactor: + observable_basis: np.ndarray + nullspace_basis: np.ndarray + observable_information: np.ndarray + chart: _Chart + + @property + def rank(self) -> int: + return int(self.observable_basis.shape[1]) + + @property + def information_matrix(self) -> np.ndarray: + return ( + self.observable_basis + @ self.observable_information + @ self.observable_basis.T + ) + + +def _twist_ambiguous_factor() -> ObservableGaugeFactor: + identity = np.eye(7) + observable = np.delete(identity, 1, axis=1) + nullspace = identity[:, 1:2] + factor = _PartialFactor( + observable_basis=observable, + nullspace_basis=nullspace, + observable_information=10.0 * np.eye(6), + chart=_Chart( + linearization=Sim3.identity(), + reference_centroid=np.zeros(3), + cloud_scale=1.0, + ), + ) + return cast(ObservableGaugeFactor, factor) + + +def test_query_conditioning_distinguishes_support_from_off_support() -> None: + factor = _twist_ambiguous_factor() + prior = np.eye(7) + supported = evaluate_query_observability( + factor, + prior_covariance_local=prior, + query_jacobian_local=point_position_query_jacobian( + factor, + np.array([1.0, 0.0, 0.0]), + ), + ) + off_support = evaluate_query_observability( + factor, + prior_covariance_local=prior, + query_jacobian_local=point_position_query_jacobian( + factor, + np.array([0.0, 5.0, 0.0]), + ), + ) + + assert supported.direct_observability_fraction == pytest.approx(1.0) + assert supported.nullspace_sensitivity_fraction == pytest.approx(0.0) + assert supported.metric_variance_reduction_fraction > 0.90 + assert supported.worst_supported_variance_ratio < 0.10 + + assert off_support.direct_observability_fraction < 0.70 + assert off_support.nullspace_sensitivity_fraction > 0.30 + assert off_support.metric_variance_reduction_fraction < 0.65 + assert off_support.worst_supported_variance_ratio > 0.95 + + gate = QueryObservabilityGate( + minimum_direct_observability_fraction=0.80, + minimum_metric_variance_reduction_fraction=0.80, + maximum_worst_supported_variance_ratio=0.50, + ) + assert gate.evaluate(supported).admitted + decision = gate.evaluate(off_support) + assert not decision.admitted + assert decision.reason_codes == ( + "insufficient-direct-query-observability", + "insufficient-query-variance-reduction", + "excessive-worst-direction-variance-ratio", + ) + + +def test_prior_mediated_reduction_is_not_called_direct_observability() -> None: + factor = _twist_ambiguous_factor() + prior = np.eye(7) + prior[1, 2] = 0.8 + prior[2, 1] = 0.8 + nullspace_query = factor.nullspace_basis.T + + report = evaluate_query_observability( + factor, + prior_covariance_local=prior, + query_jacobian_local=nullspace_query, + ) + + assert report.direct_observability_fraction == pytest.approx(0.0) + assert report.nullspace_sensitivity_fraction == pytest.approx(1.0) + assert report.metric_variance_reduction_fraction > 0.50 + assert report.worst_supported_variance_ratio < 0.50 + + +def test_declared_query_metric_controls_multioutput_weighting() -> None: + factor = _twist_ambiguous_factor() + observable_query = factor.observable_basis[:, 0] + nullspace_query = factor.nullspace_basis[:, 0] + jacobian = np.vstack((observable_query, nullspace_query)) + + identity_metric = evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=jacobian, + ) + nullspace_weighted = evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=jacobian, + query_metric=np.diag([1.0, 4.0]), + ) + + assert identity_metric.direct_observability_fraction == pytest.approx(0.5) + assert nullspace_weighted.direct_observability_fraction == pytest.approx(0.2) + assert ( + nullspace_weighted.metric_variance_reduction_fraction + < identity_metric.metric_variance_reduction_fraction + ) + + +def test_full_rank_factor_has_no_direct_nullspace_sensitivity() -> None: + identity = np.eye(7) + factor = cast( + ObservableGaugeFactor, + _PartialFactor( + observable_basis=identity, + nullspace_basis=np.empty((7, 0)), + observable_information=4.0 * identity, + chart=_Chart( + linearization=Sim3.identity(), + reference_centroid=np.zeros(3), + cloud_scale=1.0, + ), + ), + ) + jacobian = np.array( + [ + [1.0, 0.2, 0.0, -0.1, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.3, 0.1, 0.0, 1.0, 0.0], + ] + ) + + report = evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=jacobian, + ) + + assert report.direct_observability_fraction == pytest.approx(1.0) + assert report.nullspace_sensitivity_fraction == pytest.approx(0.0) + assert report.metric_variance_reduction_fraction == pytest.approx(0.8) + assert report.worst_supported_variance_ratio == pytest.approx(0.2) + + +def test_point_query_jacobian_is_read_only_and_validated() -> None: + factor = _twist_ambiguous_factor() + + jacobian = point_position_query_jacobian( + factor, + np.array([1.0, 2.0, 3.0]), + ) + + assert jacobian.shape == (3, 7) + assert not jacobian.flags.writeable + with pytest.raises(ValueError, match="shape"): + point_position_query_jacobian(factor, np.ones(2)) + with pytest.raises(ValueError, match="finite"): + point_position_query_jacobian( + factor, + np.array([0.0, np.nan, 0.0]), + ) + + +def test_invalid_query_metric_and_prior_fail_closed() -> None: + factor = _twist_ambiguous_factor() + jacobian = np.eye(7)[0:1] + + with pytest.raises(ValueError, match="positive definite"): + evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=jacobian, + query_metric=np.zeros((1, 1)), + ) + with pytest.raises(ValueError, match="positive definite"): + evaluate_query_observability( + factor, + prior_covariance_local=np.zeros((7, 7)), + query_jacobian_local=jacobian, + ) + with pytest.raises(ValueError, match=r"shape \(Q, 7\)"): + evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=np.ones((2, 6)), + ) + + +def test_analytic_control_exposes_fabricated_full_rank_information() -> None: + result = run_query_observability_study() + supported = result["results"]["rank_six_supported_query"] + off_support = result["results"]["rank_six_off_support_query"] + invalid_completion = result["results"][ + "invalid_full_rank_off_support_query" + ] + + assert supported["admitted"] + assert not off_support["admitted"] + assert invalid_completion["admitted"] + assert off_support["worst_supported_variance_ratio"] > 0.95 + assert invalid_completion["worst_supported_variance_ratio"] < 0.10 + + +def test_checked_control_evidence_matches_implementation() -> None: + path = ( + Path(__file__).parents[1] + / "evidence" + / "query-observability-control-v1" + / "result.json" + ) + + checked = json.loads(path.read_text(encoding="utf-8")) + + assert checked == run_query_observability_study()