From 558458149289dbc2b714aade8f34ff0aaff128ab Mon Sep 17 00:00:00 2001 From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com> Date: Sun, 30 Aug 2026 16:02:21 +0800 Subject: [PATCH] Add finite-orbit decision certificates beyond local gauge observability Add exact shared-orbit query bounds, robust candidate-versus-fallback advantage certificates, bounded-anchor support, constructed decision-value evidence, existing-API integration tests, documentation, and focused CI on top of the nonlinear axial-gauge representation. --- .github/workflows/finite-orbit-query.yml | 44 ++ docs/finite-orbit-query-certificates.md | 228 +++++++++ src/prob4d/axial_query_certificate.py | 446 ++++++++++++++++++ src/prob4d/axial_query_study.py | 259 ++++++++++ tests/test_axial_query_certificate.py | 356 ++++++++++++++ ...est_axial_query_certificate_integration.py | 117 +++++ 6 files changed, 1450 insertions(+) create mode 100644 .github/workflows/finite-orbit-query.yml create mode 100644 docs/finite-orbit-query-certificates.md create mode 100644 src/prob4d/axial_query_certificate.py create mode 100644 src/prob4d/axial_query_study.py create mode 100644 tests/test_axial_query_certificate.py create mode 100644 tests/test_axial_query_certificate_integration.py diff --git a/.github/workflows/finite-orbit-query.yml b/.github/workflows/finite-orbit-query.yml new file mode 100644 index 000000000..b3b79e892 --- /dev/null +++ b/.github/workflows/finite-orbit-query.yml @@ -0,0 +1,44 @@ +name: Finite-orbit query mechanism + +on: + pull_request: + branches: [main] + paths: + - src/prob4d/axial_query_certificate.py + - src/prob4d/axial_query_study.py + - tests/test_axial_query_certificate*.py + - docs/finite-orbit-query-certificates.md + - .github/workflows/finite-orbit-query.yml + +permissions: + contents: read + +jobs: + mechanism: + runs-on: ubuntu-latest + timeout-minutes: 10 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0 + with: + python-version: "3.12" + - name: Install lightweight development dependencies + run: python -m pip install -e ".[dev]" + - name: Check new source and tests + run: | + python -m ruff check src/prob4d/axial_query_certificate.py src/prob4d/axial_query_study.py tests/test_axial_query_certificate*.py + # Current NumPy stubs use Python 3.12 syntax. Runtime-floor coverage remains in Tests. + python -m mypy --python-version 3.12 src/prob4d/axial_query_certificate.py src/prob4d/axial_query_study.py + python -m pytest -q tests/test_axial_query_certificate.py tests/test_axial_query_certificate_integration.py tests/test_query_observability.py --junitxml=outputs/finite-orbit-query-tests.xml + - name: Execute constructed mechanism only + run: | + python -m prob4d.axial_query_study --seed 73029 --cases-per-family 512 --output outputs/finite-orbit-query-v1/result.json + - uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7 + if: always() + with: + name: finite-orbit-query-mechanism + path: outputs/ + if-no-files-found: warn + retention-days: 14 diff --git a/docs/finite-orbit-query-certificates.md b/docs/finite-orbit-query-certificates.md new file mode 100644 index 000000000..68ba745a9 --- /dev/null +++ b/docs/finite-orbit-query-certificates.md @@ -0,0 +1,228 @@ +# Finite-orbit query certificates + +## Scientific purpose and scope + +The local diagnostic in `query_observability.py` reports the sensitivity of a +query Jacobian to observable and null directions. It does not assert a global +invariance theorem. This extension asks a different question: can an action +comparison remain valid for **every finite rotation left unresolved by an exact +line-shaped overlap**? + +`prob4d.axial_query_certificate` preserves that ambiguity rather than completing +a deficient information matrix or treating a zero derivative as a finite-angle +certificate. This is an experimental, conditional-model kernel. It is not a +replacement for source support, identity, covariance, or provider-value gates. + +## 1. Exact correspondence geometry + +Assume labeled point correspondences, at least two distinct source points on one +line, exact agreement, positive Sim(3) scale, and proper rotations. Let `g0` be +one exact fit. Any other fit composed with `inverse(g0)` fixes every target +correspondence. The distance between two distinct fixed points forces relative +scale one; their difference forces the relative rotation to fix the line's +unit direction. The remaining transformations are exactly rotations about the +target line. Conversely, every such rotation fixes the target correspondences. + +This identifies an SO(2) stabilizer. It is **not** implied by a generic rank-six +matrix. With noisy or near-collinear correspondences, the exact-fit argument no +longer applies. `maximum_support_displacement` reports twice the maximum radial +distance to the declared line; it does not silently replace a nearly straight +cloud by an exact symmetry using a numerical rank threshold. + +## 2. A local false assurance + +Take the overlap on the z-axis and a probe at `(1, 0, 0)`. Its x-coordinate is + +```text +q(theta) = cos(theta). +``` + +At the identity representative its derivative in the unresolved rotation is +zero, but its full-orbit range is `[-1, 1]`. Inflating variance in that direction +cannot change the first-order query covariance when the corresponding Jacobian +column is zero. + +Under the analytic rank-six factor with information `10 P_observed` and complete +prior `I_7`, the existing scalar query diagnostic returns direct support `1`, +variance reduction `10/11`, and worst variance ratio `1/11`. It passes thresholds +`(0.8, 0.8, 0.5)`. The new integration test calls the existing API and exhibits an +affine fallback-minus-candidate loss `0.25 + cos(theta)`, whose lower bound is +`-0.75`. Local support therefore cannot substitute for a finite-orbit action +certificate. This is a limitation of the local approximation, not evidence that +the existing implementation calculates its declared local quantities wrongly. + +## 3. Exact finite-angle bounds + +For line origin `o`, unit direction `u`, and representative point `x_j`, write + +```text +p_j = u u^T (x_j - o) +r_j = (I - u u^T) (x_j - o) +x_j(theta) = o + p_j + cos(theta) r_j + sin(theta) (u cross r_j). +``` + +Every scalar affine query `q = b + sum_j w_j^T x_j(theta)` has the form + +```text +q(theta) = C + A cos(theta) + B sin(theta), +C = b + sum_j w_j^T (o + p_j), +A = sum_j w_j^T r_j, +B = sum_j w_j^T (u cross r_j). +``` + +Its full-circle extrema are exactly `C +/- hypot(A, B)`. On a closed circular +arc, it suffices to evaluate the two endpoints and the maximum/minimum stationary +angles that belong to the arc. The implementation handles wraparound and +singleton arcs. A zero derivative at zero establishes `B = 0`, not `A = 0`. + +This also covers isotropic squared-distance objectives to fixed target points: +`||x_j(theta) - y_j||^2` is a first harmonic because the rotating radial component +has constant norm. Its coefficients are + +```text +C_j = ||o + p_j - y_j||^2 + ||r_j||^2, +A_j = 2 (o + p_j - y_j)^T r_j, +B_j = 2 (o + p_j - y_j)^T (u cross r_j). +``` + +Nonnegative weighted sums retain this form. General simulator rollouts, +anisotropic quadratic objectives, contact changes, and arbitrary nonlinear +queries do not automatically have this representation. + +## 4. Compare actions under the SAME unresolved gauge + +Let `L_f` and `L_c` denote fallback and candidate losses with a shared latent +angle and the same representative convention. Form `D = L_f - L_c` **before** +optimizing over the angle. Bounding the two losses independently loses their +shared dependence. For example, + +```text +L_f(theta) = 4 + 2 cos(theta) + 3 sin(theta), +L_c(theta) = 3.75 + 2 cos(theta) + 3 sin(theta) +``` + +has advantage exactly `0.25`, although independent intervals overlap widely. + +An explicit `shared_gauge_id` is mandatory. Matching origins and axes does not +prove shared uncertainty: two independent gauge variables with identical +geometry must not cancel. The key also binds the geometric convention. A +source/protocol owner must justify the identity and representative; the kernel +does not infer them from coincident numbers. + +**Conditional dominance proposition.** Suppose the true shared angle belongs to +the supplied nonempty arc, and the true advantage differs from the represented +`D(theta)` by at most `epsilon`, uniformly on that support. If + +```text +min_arc D(theta) - epsilon > required_margin, +``` + +then the candidate improves the fallback loss by more than `required_margin` +for every admitted state. Proof: subtract the uniform error bound pointwise and +then use the minimum. The real-arithmetic full-circle criterion is +`C_D - hypot(A_D, B_D) - epsilon > required_margin`. + +The implementation also requires a caller-owned `scope_admitted` flag and a +strict numerical slack. A false scope or empty support always rejects. The flag +is not evidence on its own. `epsilon` must cover all omitted effects on the +advantage: observed gauge coordinates, model discrepancy, unmodeled nonlinear +responses, and relevant action uncertainty. It is not interchangeable with a +standard deviation. A Gaussian prior with unbounded support does not imply a +finite deterministic envelope. A statistically constructed support would need +its own coverage argument, including selection and dependence. + +Floating-point formulas are not outward-rounded interval arithmetic. The +explicit numerical slack is an engineering tolerance in loss units, not a +proved rounding-error enclosure or a deployment-safety theorem. + +## 5. One bounded-error metric anchor + +A separately admitted point anchor restricts angles by +`||x_a(theta) - y_a|| <= radius`. Decomposing the points into axial and radial +components gives + +```text +squared_residual(theta) = K - 2 (a cos(theta) + b sin(theta)). +``` + +Thus its feasible support is a circular arc obtained from +`cos(theta - atan2(b,a)) >= (K - radius^2)/(2 hypot(a,b))`. The implementation +returns a full circle for an uninformative feasible anchor and `None` for +inconsistent support; it never accepts vacuously from an empty set. An on-axis +anchor cannot resolve rotation. The radius is a supplied bound, not a learned +confidence interval. Reusing a visual point as an independent new sensor is not +justified. + +## Usage and downstream ownership + +```python +import numpy as np +from prob4d.axial_query_certificate import ( + AxialRotationOrbit, + certify_shared_orbit_advantage, +) + +orbit = AxialRotationOrbit( + origin=np.zeros(3), + axis=np.array([0.0, 0.0, 1.0]), + shared_gauge_id="source-bound-factor-and-representative-id", +) +fallback = orbit.affine_query([[1, 0, 0]], [[0, 0, 0]], offset=4.0) +candidate = orbit.affine_query([[1, 0, 0]], [[-1, 0, 0]], offset=3.75) +certificate = certify_shared_orbit_advantage( + fallback_loss=fallback, + candidate_loss=candidate, + scope_admitted=True, # justified analytically for this exact constructed example +) +assert not certificate.admitted +assert certificate.lower_advantage == -0.75 +``` + +In a real integration, Prob4D supplies an admitted geometric factor and its +lineage; the application supplies query coefficients and a justified envelope. +The certificate is only an additional guard input. BayesianPhysTwin's +`inference.v2` policy must still construct its own `CompleteBeliefGuardDecisionV1` +and route a rejection to the exact original complete baseline belief. This +module does not replace that router or authorize a latent-state correction. +Causal4D can consume a selected BayesianPhysTwin belief, but this extension does +not establish counterfactual or intervention validity. + +## Reproduction and evidence boundary + +```bash +python -m pytest -q tests/test_axial_query_certificate.py \ + tests/test_axial_query_certificate_integration.py +python -m prob4d.axial_query_study --seed 73029 --cases-per-family 512 \ + --output outputs/finite-orbit-query-v1/result.json +``` + +The four constructed families deliberately distinguish a stationary but +ambiguous query, a positive whole-orbit margin, shared-gauge cancellation, and +an omitted-effect envelope. They are failure controls, not a sampled deployment +population. The separate anchor arm has one additional observation and is not +an equal-information comparator. The study's standalone object-identity check +is not execution of BayesianPhysTwin's full-belief router. Paper-facing JSON, +run identities, and tables belong in `FlorianPfaff/BayesianPhysTwin-Paper`. + +The next empirical step is an additional source-developed, frozen method arm +for the PointWorld/Flat'n'Fold path tracked in #333, followed by a fresh +object/garment-disjoint evaluation only when its existing source gates authorize +it. Freeze one query, correct shared identities, an admitted support/envelope, +a coverage/width trade-off, fallback, and all decision thresholds before target +access. Report accepted-update harm, acceptance, proper scores where a genuine +predictive distribution is available, and clustered decision value. Preserve +all unsupported cases and negatives. Do not reopen any terminal MotionCrafter, +Deform360, or CUT3R cohort or call this mechanism result provider competence. + +## Related work and defensible positioning + +Symmetry-aware pose sets are established prior work: Bregier et al., *Defining +the Pose of any 3D Rigid Object and an Associated Distance*, arXiv:1612.04631. +Distributional pose estimates for uncertainty-aware manipulation are also +studied, including Jin et al., *SE(3)-PoseFlow*, arXiv:2511.01501. The proposed +contribution is not the invention of pose ambiguity, SO(2) geometry, harmonic +optimization, or robust dominance. It is a finite-orbit, shared-dependence query +certificate and explicit local-to-finite failure control for this probabilistic +4D-to-physical-query pipeline. Establishing broad novelty requires a fuller +related-work comparison, and establishing practical value requires fresh +provider evidence. diff --git a/src/prob4d/axial_query_certificate.py b/src/prob4d/axial_query_certificate.py new file mode 100644 index 000000000..6f8790651 --- /dev/null +++ b/src/prob4d/axial_query_certificate.py @@ -0,0 +1,446 @@ +"""Finite-angle query bounds for a shared, unresolved axial gauge. + +This is an analytic, conditional-model certificate, not a covariance estimate +or a provider-admission rule. An affine query of points rotated about one +fixed line has the form ``c + a*cos(theta) + b*sin(theta)``. Keeping that +shared angle permits exact action comparisons without completing a deficient +Sim(3) information matrix. A zero derivative at one representative does not +imply that the query is constant on the orbit. + +The caller must establish that the admitted uncertainty is contained in the +specified orbit, angle arc, and uniform advantage-error envelope. A generic +rank-six Jacobian, a near-collinear cloud, or an arbitrary nonlinear physical +query does not establish those assumptions. Rejected certificates must be +routed to the caller's original complete fallback belief by BayesianPhysTwin. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from numbers import Real +from typing import TypeAlias + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +FloatArray: TypeAlias = NDArray[np.float64] +OrbitKey: TypeAlias = tuple[str, tuple[float, ...]] +_SCOPE = "conditional-shared-axial-orbit-affine-query-v1" + + +def _scalar(value: object, name: str, *, nonnegative: bool = False) -> float: + if isinstance(value, (bool, np.bool_)) or not isinstance(value, Real): + raise TypeError(f"{name} must be a real scalar") + number = float(value) + if not math.isfinite(number): + raise ValueError(f"{name} must be finite") + if nonnegative and number < 0.0: + raise ValueError(f"{name} must be nonnegative") + return number + + +def _gauge_id(value: object) -> str: + if not isinstance(value, str) or not value.strip(): + raise ValueError("shared_gauge_id must be an explicit nonempty string") + return value + + +def _points(value: ArrayLike, name: str) -> FloatArray: + result = np.asarray(value, dtype=np.float64) + if result.ndim != 2 or result.shape[0] == 0 or result.shape[1] != 3: + raise ValueError(f"{name} must have shape (N, 3), N positive") + if not np.all(np.isfinite(result)): + raise ValueError(f"{name} must be finite") + return result + + +def _vector(value: ArrayLike, name: str) -> FloatArray: + result = np.asarray(value, dtype=np.float64) + if result.shape != (3,) or not np.all(np.isfinite(result)): + raise ValueError(f"{name} must be a finite three-vector") + # A bytes-backed array cannot have its write flag re-enabled by a caller. + return np.frombuffer(result.tobytes(), dtype=np.float64) + + +@dataclass(frozen=True) +class AngleArc: + """Closed circular arc, parameterized as center +/- half_width radians. + + The full circle is the default. An arc is a supplied support assumption, + not a confidence level. ``None`` represents an infeasible/empty arc in + the functions below; it must never yield vacuous acceptance. + """ + + center: float = 0.0 + half_width: float = math.pi + + def __post_init__(self) -> None: + center = _scalar(self.center, "center") + width = _scalar(self.half_width, "half_width", nonnegative=True) + if width > math.pi: + raise ValueError("half_width must not exceed pi") + object.__setattr__(self, "center", math.remainder(center, 2.0 * math.pi)) + object.__setattr__(self, "half_width", width) + + def contains(self, angle: float, *, atol: float = 0.0) -> bool: + angle = _scalar(angle, "angle") + tolerance = _scalar(atol, "atol", nonnegative=True) + displacement = math.remainder(angle - self.center, 2.0 * math.pi) + return abs(displacement) <= self.half_width + tolerance + + +_FULL_CIRCLE = AngleArc() + + +@dataclass(frozen=True) +class HarmonicQuery: + """Scalar affine query on one declared, shared axial orbit. + + ``orbit_key`` binds both the shared latent identity and coordinate convention. + Matching axis geometry alone does not establish shared uncertainty. + """ + + constant: float + cosine: float + sine: float + orbit_key: OrbitKey + + def __post_init__(self) -> None: + for name in ("constant", "cosine", "sine"): + object.__setattr__(self, name, _scalar(getattr(self, name), name)) + if not isinstance(self.orbit_key, tuple) or len(self.orbit_key) != 2: + raise ValueError("orbit_key must bind a shared identity and geometry") + identity, geometry = self.orbit_key + identity = _gauge_id(identity) + if len(geometry) != 6: + raise ValueError("orbit geometry must contain an origin and an axis") + key = tuple(_scalar(v, "orbit_key entry") for v in geometry) + if not math.isclose(math.hypot(*key[3:]), 1.0, abs_tol=1e-12): + raise ValueError("orbit_key must contain a unit axis") + object.__setattr__(self, "orbit_key", (identity, key)) + + @property + def amplitude(self) -> float: + value = math.hypot(self.cosine, self.sine) + if not math.isfinite(value): + raise ValueError("query amplitude overflowed") + return value + + @property + def derivative_at_zero(self) -> float: + return self.sine + + def evaluate(self, angle: float) -> float: + angle = _scalar(angle, "angle") + value = math.fsum( + (self.constant, self.cosine * math.cos(angle), self.sine * math.sin(angle)) + ) + if not math.isfinite(value): + raise ValueError("query evaluation overflowed") + return value + + def minus(self, other: HarmonicQuery) -> HarmonicQuery: + """Subtract before bounding, preserving shared-angle cancellation.""" + if not isinstance(other, HarmonicQuery): + raise TypeError("other must be a HarmonicQuery") + if self.orbit_key != other.orbit_key: + raise ValueError("queries must use the same shared orbit") + return HarmonicQuery( + self.constant - other.constant, + self.cosine - other.cosine, + self.sine - other.sine, + self.orbit_key, + ) + + def bounds(self, arc: AngleArc = _FULL_CIRCLE) -> QueryBounds: + """Analytic extrema, including interior stationary angles. + + Endpoints alone are insufficient when an arc crosses an extremum. + Floating-point results are numerical evaluations of the exact formula, + not an interval-arithmetic verification of every rounding operation. + """ + if not isinstance(arc, AngleArc): + raise TypeError("arc must be a nonempty AngleArc") + radius = self.amplitude + if radius == 0.0: + return QueryBounds(self.constant, self.constant, arc.center, arc.center) + if arc.half_width == math.pi: + maximum_angle = math.atan2(self.sine, self.cosine) + minimum_angle = math.remainder(maximum_angle + math.pi, 2.0 * math.pi) + return QueryBounds( + self.constant - radius, + self.constant + radius, + minimum_angle, + maximum_angle, + ) + candidates = [arc.center - arc.half_width, arc.center + arc.half_width] + maximum_angle = math.atan2(self.sine, self.cosine) + minimum_angle = math.remainder(maximum_angle + math.pi, 2.0 * math.pi) + for angle in (minimum_angle, maximum_angle): + if arc.contains(angle): + candidates.append(angle) + values = [self.evaluate(angle) for angle in candidates] + lower_index = int(np.argmin(values)) + upper_index = int(np.argmax(values)) + return QueryBounds( + values[lower_index], + values[upper_index], + candidates[lower_index], + candidates[upper_index], + ) + + +@dataclass(frozen=True) +class QueryBounds: + """Numerical analytic extrema and angles attaining them.""" + + lower: float + upper: float + lower_angle: float + upper_angle: float + + def __post_init__(self) -> None: + for name in ("lower", "upper", "lower_angle", "upper_angle"): + object.__setattr__(self, name, _scalar(getattr(self, name), name)) + if self.lower > self.upper: + raise ValueError("lower must not exceed upper") + + +@dataclass(frozen=True, eq=False) +class AxialRotationOrbit: + """Rotations of representative metric points about a fixed reference line. + + Points must already be transformed by the fitted representative Sim(3). + Constructing this object does not assert that a visual factor identifies + the line or that the orbit exhausts a noisy seven-dimensional posterior. + ``shared_gauge_id`` must identify the same latent angle and representative + convention in both queries; it is never inferred from matching geometry. + """ + + origin: FloatArray + axis: FloatArray + shared_gauge_id: str + + def __post_init__(self) -> None: + origin = _vector(self.origin, "origin") + raw_axis = _vector(self.axis, "axis") + norm = math.hypot(*(float(x) for x in raw_axis)) + if norm == 0.0 or not math.isfinite(norm): + raise ValueError("axis must have finite positive norm") + axis = _vector(raw_axis / norm, "normalized axis") + object.__setattr__(self, "origin", origin) + object.__setattr__(self, "axis", axis) + object.__setattr__(self, "shared_gauge_id", _gauge_id(self.shared_gauge_id)) + + @property + def key(self) -> OrbitKey: + geometry = tuple(float(x) for x in np.concatenate((self.origin, self.axis))) + return self.shared_gauge_id, geometry + + def transform(self, points: ArrayLike, angle: float) -> FloatArray: + points_array = _points(points, "points") + angle = _scalar(angle, "angle") + centered = points_array - self.origin + parallel = np.outer(centered @ self.axis, self.axis) + perpendicular = centered - parallel + transformed = ( + self.origin + + parallel + + math.cos(angle) * perpendicular + + math.sin(angle) * np.cross(self.axis, perpendicular) + ) + if not np.all(np.isfinite(transformed)): + raise ValueError("transformed points overflowed") + return transformed + + def maximum_support_displacement(self, points: ArrayLike) -> float: + """Exact full-orbit maximum point displacement from the representative. + + This is twice the largest distance to the axis. It is zero for exact + line support, but explicitly nonzero for near-collinear input. There + is no rank threshold that silently promotes a near-line to an exact + stabilizer. + """ + centered = _points(points, "points") - self.origin + perpendicular = centered - np.outer(centered @ self.axis, self.axis) + maximum = 2.0 * max(math.hypot(*(float(x) for x in row)) for row in perpendicular) + if not math.isfinite(maximum): + raise ValueError("support displacement overflowed") + return maximum + + def affine_query( + self, + points: ArrayLike, + weights: ArrayLike, + *, + offset: float = 0.0, + ) -> HarmonicQuery: + """Represent ``offset + sum_j weights[j] @ rotated_points[j]``.""" + points_array = _points(points, "points") + weights_array = _points(weights, "weights") + if points_array.shape != weights_array.shape: + raise ValueError("points and weights must have the same shape") + offset = _scalar(offset, "offset") + centered = points_array - self.origin + parallel = np.outer(centered @ self.axis, self.axis) + perpendicular = centered - parallel + return HarmonicQuery( + constant=offset + float(np.sum(weights_array * (self.origin + parallel))), + cosine=float(np.sum(weights_array * perpendicular)), + sine=float(np.sum(weights_array * np.cross(self.axis, perpendicular))), + orbit_key=self.key, + ) + + def bounded_anchor_arc( + self, + representative_point: ArrayLike, + observed_point: ArrayLike, + *, + error_radius: float, + ) -> AngleArc | None: + """Angles consistent with one Euclidean bounded-error metric anchor. + + Returns ``None`` for inconsistent support and a full circle for an + uninformative feasible anchor. The supplied radius is a deterministic + bound; this method does not calibrate it or assign a confidence level. + Reusing the visual observation as an independent anchor is not justified. + """ + point = _vector(representative_point, "representative_point") - self.origin + observed = _vector(observed_point, "observed_point") - self.origin + radius = _scalar(error_radius, "error_radius", nonnegative=True) + point_parallel = float(point @ self.axis) + observed_parallel = float(observed @ self.axis) + point_perpendicular = point - point_parallel * self.axis + observed_perpendicular = observed - observed_parallel * self.axis + cosine = float(point_perpendicular @ observed_perpendicular) + sine = float(np.cross(self.axis, point_perpendicular) @ observed_perpendicular) + amplitude = math.hypot(cosine, sine) + constant = math.fsum( + ( + (point_parallel - observed_parallel) ** 2, + float(point_perpendicular @ point_perpendicular), + float(observed_perpendicular @ observed_perpendicular), + ) + ) + radius_squared = radius * radius + if not all(math.isfinite(v) for v in (amplitude, constant, radius_squared)): + raise ValueError("anchor geometry overflowed") + if amplitude == 0.0: + return AngleArc() if constant <= radius_squared else None + threshold = (constant - radius_squared) / (2.0 * amplitude) + if threshold > 1.0: + return None + if threshold <= -1.0: + return AngleArc() + return AngleArc( + center=math.atan2(sine, cosine), + half_width=math.acos(threshold), + ) + + +@dataclass(frozen=True) +class OrbitAdvantageCertificate: + """A conditional advantage certificate; not a deployment authorization.""" + + admitted: bool + reason_codes: tuple[str, ...] + lower_advantage: float | None + upper_advantage: float | None + required_margin: float + advantage_error_bound: float + numerical_slack: float + scope: str = _SCOPE + + def __post_init__(self) -> None: + if type(self.admitted) is not bool: + raise TypeError("admitted must be a bool") + if self.scope != _SCOPE: + raise ValueError("certificate scope changed") + reasons = tuple(self.reason_codes) + allowed = { + "orbit-model-scope-not-admitted", + "infeasible-anchor-support", + "nonpositive-robust-advantage", + } + if len(set(reasons)) != len(reasons) or any(r not in allowed for r in reasons): + raise ValueError("invalid or duplicate rejection reasons") + if self.admitted == bool(reasons): + raise ValueError("admission and rejection reasons disagree") + for name in ("required_margin", "advantage_error_bound", "numerical_slack"): + object.__setattr__(self, name, _scalar(getattr(self, name), name, nonnegative=True)) + if (self.lower_advantage is None) != (self.upper_advantage is None): + raise ValueError("both advantage bounds must be present or absent") + if self.lower_advantage is None: + if "infeasible-anchor-support" not in reasons: + raise ValueError("missing bounds require infeasible support") + else: + lower = _scalar(self.lower_advantage, "lower_advantage") + upper = _scalar(self.upper_advantage, "upper_advantage") + if lower > upper: + raise ValueError("advantage bounds are reversed") + if self.admitted and not lower > self.required_margin + self.numerical_slack: + raise ValueError("admitted certificate has no positive robust advantage") + object.__setattr__(self, "reason_codes", reasons) + + +def certify_shared_orbit_advantage( + *, + fallback_loss: HarmonicQuery, + candidate_loss: HarmonicQuery, + scope_admitted: bool, + arc: AngleArc | None = _FULL_CIRCLE, + advantage_error_bound: float = 0.0, + required_margin: float = 0.0, + numerical_slack: float = 1e-12, +) -> OrbitAdvantageCertificate: + """Admit only uniformly positive fallback-minus-candidate advantage. + + Both losses must use the same orbit. Subtraction precedes optimization, + preserving shared geometric uncertainty. The uniform error bound covers + all omitted effects on their difference, not just a fitted standard error. + ``scope_admitted`` must come from the caller's independently justified + source/model gate. A false scope or empty anchor support always rejects. + + ``numerical_slack`` is an explicit absolute tolerance in loss units, not a + statistical margin or a formal interval-arithmetic error bound. Strict + comparison rejects equality and numerically indistinguishable improvement. + """ + if type(scope_admitted) is not bool: + raise TypeError("scope_admitted must be a bool") + error = _scalar(advantage_error_bound, "advantage_error_bound", nonnegative=True) + margin = _scalar(required_margin, "required_margin", nonnegative=True) + slack = _scalar(numerical_slack, "numerical_slack", nonnegative=True) + difference = fallback_loss.minus(candidate_loss) + reasons: list[str] = [] + if not scope_admitted: + reasons.append("orbit-model-scope-not-admitted") + lower: float | None = None + upper: float | None = None + if arc is None: + reasons.append("infeasible-anchor-support") + else: + bounds = difference.bounds(arc) + lower = _scalar(bounds.lower - error, "lower advantage") + upper = _scalar(bounds.upper + error, "upper advantage") + if not lower > margin + slack: + reasons.append("nonpositive-robust-advantage") + return OrbitAdvantageCertificate( + admitted=not reasons, + reason_codes=tuple(reasons), + lower_advantage=lower, + upper_advantage=upper, + required_margin=margin, + advantage_error_bound=error, + numerical_slack=slack, + ) + + +__all__ = [ + "AngleArc", + "AxialRotationOrbit", + "HarmonicQuery", + "OrbitAdvantageCertificate", + "QueryBounds", + "certify_shared_orbit_advantage", +] diff --git a/src/prob4d/axial_query_study.py b/src/prob4d/axial_query_study.py new file mode 100644 index 000000000..3f86a59b0 --- /dev/null +++ b/src/prob4d/axial_query_study.py @@ -0,0 +1,259 @@ +"""Controlled finite-orbit failure and decision-value study; no provider data. + +Run with ``python -m prob4d.axial_query_study --output result.json``. +The study deliberately constructs its model assumptions. It is not a +calibration experiment, a fresh cohort, or a comparison of real providers. +""" + +from __future__ import annotations + +import argparse +import json +import math +from pathlib import Path +from typing import Any + +import numpy as np + +from .axial_query_certificate import ( + AxialRotationOrbit, + HarmonicQuery, + certify_shared_orbit_advantage, +) + +FAMILIES = ( + "stationary-but-ambiguous", + "positive-whole-orbit-margin", + "shared-gauge-cancellation", + "bounded-omitted-discrepancy", +) +ARMS = ( + "local-query-gate-then-plugin", + "reject-all-deficient", + "independent-query-intervals", + "shared-orbit-certificate", + "shared-orbit-plus-one-bounded-anchor", +) + + +def _summary() -> dict[str, int | float]: + return { + "cases": 0, + "accepted": 0, + "sampled_harmful_accepts": 0, + "admitted_with_possible_harm": 0, + "exact_fallback_identity_failures": 0, + "deployed_advantage_sum": 0.0, + "maximum_possible_harm": 0.0, + } + + +def analytic_local_control(orbit: AxialRotationOrbit, radial: np.ndarray) -> dict[str, float]: + """Reference algebra for the existing scalar local query gate. + + The complete prior is I_7 and visual information is ten times the rank-six + projector excluding rotation about the line. The integration test checks + this control against the repository's actual query_observability API. + """ + normal = radial / np.linalg.norm(radial) + jacobian = np.concatenate(([float(normal @ radial)], np.cross(radial, normal), normal)) + nullspace = np.concatenate(([0.0], orbit.axis, np.zeros(3))) + projector = np.eye(7) - np.outer(nullspace, nullspace) + posterior = np.linalg.solve(np.eye(7) + 10.0 * projector, np.eye(7)) + energy = float(jacobian @ jacobian) + direct = float(jacobian @ projector @ jacobian) / energy + ratio = float(jacobian @ posterior @ jacobian) / energy + return { + "direct_observability_fraction": direct, + "metric_variance_reduction_fraction": 1.0 - ratio, + "worst_supported_variance_ratio": ratio, + } + + +def run_axial_query_study(*, seed: int = 73029, cases_per_family: int = 512) -> dict[str, Any]: + if isinstance(cases_per_family, bool) or not isinstance(cases_per_family, int): + raise TypeError("cases_per_family must be an integer") + if cases_per_family < 1: + raise ValueError("cases_per_family must be positive") + if isinstance(seed, bool) or not isinstance(seed, int) or seed < 0: + raise ValueError("seed must be a nonnegative integer") + rng = np.random.default_rng(seed) + results = {family: {arm: _summary() for arm in ARMS} for family in FAMILIES} + max_support_drift = 0.0 + max_local_nullspace_fraction = 0.0 + anchor_truth_exclusion_count = 0 + for family in FAMILIES: + for case_index in range(cases_per_family): + axis = rng.normal(size=3) + axis /= np.linalg.norm(axis) + origin = rng.normal(scale=0.2, size=3) + direction = rng.normal(size=3) + direction -= axis * float(direction @ axis) + direction /= np.linalg.norm(direction) + radius = float(rng.uniform(0.01, 0.10)) + radial = radius * direction + orbit = AxialRotationOrbit(origin, axis, f"{family}/{case_index}") + support = origin + np.array([-0.4, -0.2, 0.0, 0.2, 0.4])[:, None] * axis + max_support_drift = max(max_support_drift, orbit.maximum_support_displacement(support)) + q = orbit.affine_query( + (origin + radial)[None, :], + direction[None, :], + offset=-float(direction @ origin), + ) + base = 4.0 * radius + error = 0.0 + if family == "shared-gauge-cancellation": + fallback = HarmonicQuery( + base + 2.0 * q.constant, 2.0 * q.cosine, 2.0 * q.sine, orbit.key + ) + candidate = HarmonicQuery( + base - 0.25 * radius + 2.0 * q.constant, + 2.0 * q.cosine, + 2.0 * q.sine, + orbit.key, + ) + else: + gap = (0.25 if family == "stationary-but-ambiguous" else 1.25) * radius + fallback = HarmonicQuery(base, 0.0, 0.0, orbit.key) + candidate = HarmonicQuery(base - gap - q.constant, -q.cosine, -q.sine, orbit.key) + if family == "bounded-omitted-discrepancy": + error = 0.5 * radius + difference = fallback.minus(candidate) + true_angle = float(rng.uniform(-math.pi, math.pi)) + # Deliberately use the worst signed omitted effect. This belongs to + # the stated uniform envelope, not to an estimated noise model. + true_advantage = difference.evaluate(true_angle) - error + full_bounds = difference.bounds() + possible_harm = max(0.0, error - full_bounds.lower) + local = analytic_local_control(orbit, radial) + max_local_nullspace_fraction = max( + max_local_nullspace_fraction, 1.0 - local["direct_observability_fraction"] + ) + local_admitted = ( + local["direct_observability_fraction"] >= 0.8 + and local["metric_variance_reduction_fraction"] >= 0.8 + and local["worst_supported_variance_ratio"] <= 0.5 + and difference.evaluate(0.0) > 1e-12 + ) + certificate = certify_shared_orbit_advantage( + fallback_loss=fallback, + candidate_loss=candidate, + scope_admitted=True, + advantage_error_bound=error, + ) + independent_lower = fallback.bounds().lower - candidate.bounds().upper - error + # Separate information-budget arm: one synthetic bounded-error + # metric anchor, with known radius and no fitted calibration. + anchor_reference = origin + 0.15 * direction + anchor_observation = orbit.transform(anchor_reference[None, :], true_angle)[0] + anchor_arc = orbit.bounded_anchor_arc( + anchor_reference, anchor_observation, error_radius=0.02 + ) + if anchor_arc is None or not anchor_arc.contains(true_angle, atol=1e-12): + anchor_truth_exclusion_count += 1 + anchored = certify_shared_orbit_advantage( + fallback_loss=fallback, + candidate_loss=candidate, + scope_admitted=True, + arc=anchor_arc, + advantage_error_bound=error, + ) + decisions = { + ARMS[0]: local_admitted, + ARMS[1]: False, + ARMS[2]: independent_lower > 1e-12, + ARMS[3]: certificate.admitted, + ARMS[4]: anchored.admitted, + } + fallback_object = object() + candidate_object = object() + for arm, admitted in decisions.items(): + row = results[family][arm] + row["cases"] += 1 + row["accepted"] += int(admitted) + row["sampled_harmful_accepts"] += int(admitted and true_advantage < -1e-12) + bound_harm = possible_harm + if arm == ARMS[4]: + bound_harm = max(0.0, -(anchored.lower_advantage or 0.0)) + row["admitted_with_possible_harm"] += int(admitted and bound_harm > 1e-12) + if admitted: + row["deployed_advantage_sum"] += true_advantage + row["maximum_possible_harm"] = max(row["maximum_possible_harm"], bound_harm) + selected = candidate_object if admitted else fallback_object + row["exact_fallback_identity_failures"] += int( + not admitted and selected is not fallback_object + ) + totals = {arm: _summary() for arm in ARMS} + for family_rows in results.values(): + for arm, row in family_rows.items(): + for key, value in row.items(): + if key == "maximum_possible_harm": + totals[arm][key] = max(totals[arm][key], value) + else: + totals[arm][key] += value + for group in (totals, *results.values()): + for row in group.values(): + row["acceptance_fraction"] = row["accepted"] / row["cases"] + row["mean_deployed_advantage"] = row.pop("deployed_advantage_sum") / row["cases"] + return { + "schema": "prob4d.finite-orbit-query-mechanism-study", + "schema_version": 1, + "seed": seed, + "cases_per_family": cases_per_family, + "total_cases": len(FAMILIES) * cases_per_family, + "evidence_class": "constructed-controlled-mechanism-not-real-provider", + "uncertainty_scope": "one-shared-axial-orbit-plus-declared-uniform-advantage-error", + "local_gate_thresholds": {"direct_min": 0.8, "reduction_min": 0.8, "worst_ratio_max": 0.5}, + "additional_anchor_arm": { + "anchor_radius_m": 0.15, + "error_bound_m": 0.02, + "actual_sensor_error_m": 0.0, + "information_budget": "one-extra-synthetic-metric-anchor", + }, + "maximum_support_orbit_drift_m": max_support_drift, + "maximum_local_nullspace_fraction": max_local_nullspace_fraction, + "anchor_truth_exclusion_count": anchor_truth_exclusion_count, + "totals": totals, + "families": results, + "claim_boundary": [ + "The constructed orbit and error envelope contain the simulated truth by design.", + ( + "Zero harmful accepts is a conditional algebraic mechanism result, " + "not measured real-world safety or calibration." + ), + ( + "The local arm is an analytic reference control checked against the " + "existing API by a separate integration test." + ), + ( + "The anchor arm has one additional observation " + "and is not an equal-information comparison." + ), + ( + "The object-identity check is a standalone routing control, not execution " + "of BayesianPhysTwin's complete-belief router." + ), + ( + "No real provider, protected source, target, physical simulator, " + "or Causal4D outcome was accessed." + ), + ], + } + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--seed", type=int, default=73029) + parser.add_argument("--cases-per-family", type=int, default=512) + parser.add_argument("--output", type=Path, required=True) + args = parser.parse_args() + result = run_axial_query_study(seed=args.seed, cases_per_family=args.cases_per_family) + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text( + json.dumps(result, indent=2, sort_keys=True, allow_nan=False) + "\n", + encoding="utf-8", + ) + + +if __name__ == "__main__": + main() diff --git a/tests/test_axial_query_certificate.py b/tests/test_axial_query_certificate.py new file mode 100644 index 000000000..24e77b63b --- /dev/null +++ b/tests/test_axial_query_certificate.py @@ -0,0 +1,356 @@ +from __future__ import annotations + +import math +from dataclasses import replace + +import numpy as np +import pytest + +from prob4d.axial_query_certificate import ( + AngleArc, + AxialRotationOrbit, + HarmonicQuery, + certify_shared_orbit_advantage, +) +from prob4d.axial_query_study import run_axial_query_study + + +def orbit() -> AxialRotationOrbit: + return AxialRotationOrbit(np.zeros(3), np.array([0.0, 0.0, 1.0]), "shared-test-gauge") + + +def query(c: float, a: float, b: float) -> HarmonicQuery: + return HarmonicQuery(c, a, b, orbit().key) + + +def test_stationary_query_is_not_constant_on_finite_orbit() -> None: + model = orbit() + support = np.array([[0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]) + q = model.affine_query([[1.0, 0.0, 0.0]], [[1.0, 0.0, 0.0]]) + assert model.maximum_support_displacement(support) == 0.0 + assert q.derivative_at_zero == 0.0 + assert q.amplitude == 1.0 + assert q.bounds().lower == -1.0 + assert q.bounds().upper == 1.0 + for angle in np.linspace(-math.pi, math.pi, 31): + np.testing.assert_array_equal(model.transform(support, float(angle)), support) + assert q.evaluate(0.0) == 1.0 + assert q.evaluate(math.pi) == -1.0 + + +def test_shared_gauge_cancels_before_interval_optimization() -> None: + fallback = query(4.0, 2.0, 3.0) + candidate = query(3.75, 2.0, 3.0) + assert fallback.bounds().lower - candidate.bounds().upper < 0.0 + result = certify_shared_orbit_advantage( + fallback_loss=fallback, candidate_loss=candidate, scope_admitted=True + ) + assert result.admitted + assert result.lower_advantage == 0.25 + assert result.upper_advantage == 0.25 + + +def test_ambiguous_query_can_still_have_a_certified_decision() -> None: + result = certify_shared_orbit_advantage( + fallback_loss=query(4.0, 0.0, 0.0), + candidate_loss=query(2.75, -1.0, 0.0), + scope_admitted=True, + ) + assert result.admitted + assert result.lower_advantage == 0.25 + + +def test_omitted_effect_envelope_can_reverse_admission() -> None: + kwargs = dict( + fallback_loss=query(4.0, 0.0, 0.0), + candidate_loss=query(2.75, -1.0, 0.0), + scope_admitted=True, + ) + assert certify_shared_orbit_advantage(**kwargs).admitted + result = certify_shared_orbit_advantage(**kwargs, advantage_error_bound=0.5) + assert not result.admitted + assert result.lower_advantage == -0.25 + assert result.reason_codes == ("nonpositive-robust-advantage",) + + +def test_unadmitted_scope_and_empty_support_never_accept_vacuously() -> None: + result = certify_shared_orbit_advantage( + fallback_loss=query(10.0, 0.0, 0.0), + candidate_loss=query(0.0, 0.0, 0.0), + scope_admitted=False, + arc=None, + ) + assert not result.admitted + assert result.lower_advantage is None + assert result.upper_advantage is None + assert result.reason_codes == ( + "orbit-model-scope-not-admitted", + "infeasible-anchor-support", + ) + + +def test_different_orbits_cannot_be_treated_as_shared_uncertainty() -> None: + first = query(2.0, 1.0, 0.0) + other_orbit = AxialRotationOrbit(np.ones(3), np.array([0.0, 1.0, 0.0]), "shared-test-gauge") + other = HarmonicQuery(1.0, 1.0, 0.0, other_orbit.key) + with pytest.raises(ValueError, match="same shared orbit"): + first.minus(other) + + +def test_arc_crossing_pi_and_interior_extrema_are_handled() -> None: + q = query(0.0, 1.0, 0.0) + wrapped = AngleArc(math.pi - 0.05, 0.20) + bounds = q.bounds(wrapped) + assert bounds.lower == pytest.approx(-1.0) + assert bounds.upper == pytest.approx(math.cos(math.pi - 0.25)) + assert wrapped.contains(-math.pi) + assert wrapped.contains(bounds.lower_angle) + assert wrapped.contains(bounds.upper_angle, atol=1e-14) + maximum = q.bounds(AngleArc(0.1, 0.4)) + assert maximum.upper == 1.0 + + +def test_zero_width_and_constant_queries() -> None: + q = query(0.3, 1.0, 2.0) + result = q.bounds(AngleArc(0.45, 0.0)) + assert result.lower == result.upper == q.evaluate(0.45) + constant = query(3.0, 0.0, 0.0).bounds(AngleArc(2.0, 0.1)) + assert constant.lower == constant.upper == 3.0 + + +def test_random_affine_queries_match_direct_transformed_points() -> None: + rng = np.random.default_rng(11) + for _ in range(100): + model = AxialRotationOrbit(rng.normal(size=3), rng.normal(size=3), "shared-test-gauge") + points = rng.normal(size=(12, 3)) + weights = rng.normal(size=(12, 3)) + offset = float(rng.normal()) + q = model.affine_query(points, weights, offset=offset) + for angle in rng.uniform(-8.0, 8.0, 8): + direct = offset + np.sum(weights * model.transform(points, float(angle))) + assert q.evaluate(float(angle)) == pytest.approx(direct, abs=1e-12) + + +def test_random_arc_bounds_contain_dense_grid_and_attain_their_extrema() -> None: + rng = np.random.default_rng(22) + for _ in range(400): + q = query(*(float(v) for v in rng.normal(size=3))) + arc = AngleArc(float(rng.uniform(-5, 5)), float(rng.uniform(0, math.pi))) + bounds = q.bounds(arc) + grid = np.linspace(arc.center - arc.half_width, arc.center + arc.half_width, 2049) + values = q.constant + q.cosine * np.cos(grid) + q.sine * np.sin(grid) + assert np.min(values) >= bounds.lower - 2e-12 + assert np.max(values) <= bounds.upper + 2e-12 + assert q.evaluate(bounds.lower_angle) == pytest.approx(bounds.lower, abs=2e-12) + assert q.evaluate(bounds.upper_angle) == pytest.approx(bounds.upper, abs=2e-12) + assert arc.contains(bounds.lower_angle, atol=2e-12) + assert arc.contains(bounds.upper_angle, atol=2e-12) + + +def test_rigid_coordinate_change_preserves_scalar_query_ranges() -> None: + rng = np.random.default_rng(33) + for _ in range(50): + matrix, _ = np.linalg.qr(rng.normal(size=(3, 3))) + matrix[:, -1] *= np.linalg.det(matrix) + shift = rng.normal(size=3) + source = AxialRotationOrbit(rng.normal(size=3), rng.normal(size=3), "shared-test-gauge") + target = AxialRotationOrbit( + matrix @ source.origin + shift, matrix @ source.axis, "shared-test-gauge" + ) + points = rng.normal(size=(8, 3)) + weights = rng.normal(size=(8, 3)) + moved_weights = weights @ matrix.T + first = source.affine_query(points, weights, offset=0.7) + second = target.affine_query( + points @ matrix.T + shift, + moved_weights, + offset=0.7 - float(np.sum(moved_weights @ shift)), + ) + np.testing.assert_allclose( + [first.constant, first.cosine, first.sine], + [second.constant, second.cosine, second.sine], + atol=2e-12, + ) + + +def test_near_line_support_is_not_silently_promoted_to_exact_symmetry() -> None: + points = np.array([[0.003, 0.0, 0.2], [-0.001, 0.0, -0.2]]) + model = orbit() + assert model.maximum_support_displacement(points) == pytest.approx(0.006) + motion = np.linalg.norm(model.transform(points, math.pi) - points, axis=1) + assert np.max(motion) == pytest.approx(model.maximum_support_displacement(points)) + + +def test_off_axis_bounded_anchor_has_exact_feasible_arc() -> None: + model = orbit() + arc = model.bounded_anchor_arc([1.0, 0.0, 0.0], [1.0, 0.0, 0.0], error_radius=1.0) + assert arc is not None + assert arc.center == 0.0 + assert arc.half_width == pytest.approx(math.pi / 3) + bounds = query(0.25, 1.0, 0.0).bounds(arc) + assert bounds.lower == pytest.approx(0.75) + assert bounds.upper == pytest.approx(1.25) + + +def test_on_axis_anchor_is_either_uninformative_or_inconsistent() -> None: + model = orbit() + arc = model.bounded_anchor_arc([0.0, 0.0, 1.0], [0.0, 0.0, 1.0], error_radius=0.0) + assert arc == AngleArc() + assert model.bounded_anchor_arc( + [0.0, 0.0, 1.0], [0.0, 0.0, 2.0], error_radius=0.5 + ) is None + + +def test_anchor_full_circle_singleton_and_infeasible_boundaries() -> None: + model = orbit() + assert model.bounded_anchor_arc( + [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], error_radius=2.0 + ) == AngleArc() + singleton = model.bounded_anchor_arc( + [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], error_radius=0.0 + ) + assert singleton is not None + assert singleton.half_width == 0.0 + assert singleton.center == pytest.approx(math.pi / 2) + assert model.bounded_anchor_arc( + [1.0, 0.0, 0.0], [2.0, 0.0, 0.0], error_radius=0.5 + ) is None + + +def test_anchor_arc_matches_direct_residual_and_contains_bounded_sensor_truth() -> None: + rng = np.random.default_rng(44) + for _ in range(200): + model = AxialRotationOrbit(rng.normal(size=3), rng.normal(size=3), "shared-test-gauge") + point = rng.normal(size=3) + theta = float(rng.uniform(-math.pi, math.pi)) + radius = float(rng.uniform(0.001, 0.3)) + noise = rng.normal(size=3) + noise *= radius * float(rng.uniform(0.0, 0.9)) / np.linalg.norm(noise) + observed = model.transform(point[None, :], theta)[0] + noise + arc = model.bounded_anchor_arc(point, observed, error_radius=radius) + assert arc is not None and arc.contains(theta, atol=5e-11) + for candidate in rng.uniform(-math.pi, math.pi, 30): + residual = np.linalg.norm( + model.transform(point[None, :], float(candidate))[0] - observed + ) + if abs(residual - radius) > 1e-10: + assert arc.contains(float(candidate)) == (residual <= radius) + + +def test_axes_are_normalized_and_arrays_irreversibly_readonly() -> None: + original = np.array([0.0, 0.0, 8.0]) + model = AxialRotationOrbit(np.zeros(3), original, "shared-test-gauge") + original[:] = 1.0 + np.testing.assert_array_equal(model.axis, [0.0, 0.0, 1.0]) + for array in (model.origin, model.axis): + with pytest.raises(ValueError): + array.setflags(write=True) + + +@pytest.mark.parametrize("value", [float("nan"), float("inf"), -float("inf")]) +def test_nonfinite_coefficients_and_angles_fail_closed(value: float) -> None: + with pytest.raises(ValueError, match="finite"): + query(value, 0.0, 0.0) + with pytest.raises(ValueError, match="finite"): + AngleArc(value, 0.1) + with pytest.raises(ValueError, match="finite"): + query(0.0, 1.0, 1.0).evaluate(value) + + +@pytest.mark.parametrize("field", ["advantage_error_bound", "required_margin", "numerical_slack"]) +def test_negative_gate_controls_fail_closed(field: str) -> None: + kwargs = { + "fallback_loss": query(1.0, 0.0, 0.0), + "candidate_loss": query(0.0, 0.0, 0.0), + "scope_admitted": True, + field: -0.1, + } + with pytest.raises(ValueError, match="nonnegative"): + certify_shared_orbit_advantage(**kwargs) + + +def test_invalid_geometry_scope_arc_and_orbit_key_fail_closed() -> None: + with pytest.raises(ValueError, match="positive norm"): + AxialRotationOrbit(np.zeros(3), np.zeros(3), "shared-test-gauge") + with pytest.raises(ValueError, match="three-vector"): + AxialRotationOrbit(np.ones(2), np.ones(3), "shared-test-gauge") + with pytest.raises(ValueError, match="same shape"): + orbit().affine_query(np.ones((2, 3)), np.ones((3, 3))) + with pytest.raises(ValueError, match="shape"): + orbit().affine_query([], []) + with pytest.raises(TypeError, match="scope_admitted"): + certify_shared_orbit_advantage( + fallback_loss=query(1.0, 0.0, 0.0), + candidate_loss=query(0.0, 0.0, 0.0), + scope_admitted=1, + ) + with pytest.raises(ValueError, match="exceed pi"): + AngleArc(0.0, 4.0) + with pytest.raises(TypeError, match="real scalar"): + AngleArc(True, 0.5) + with pytest.raises(ValueError, match="unit axis"): + replace(query(0.0, 0.0, 0.0), orbit_key=("shared-test-gauge", (0.0,) * 6)) + + +def test_exact_ties_and_numerically_indistinguishable_improvements_reject() -> None: + for gain in (0.0, 1e-13, 1e-12): + result = certify_shared_orbit_advantage( + fallback_loss=query(gain, 0.0, 0.0), + candidate_loss=query(0.0, 0.0, 0.0), + scope_admitted=True, + ) + assert not result.admitted + + +def test_controlled_study_is_reproducible_and_preserves_claim_boundaries() -> None: + first = run_axial_query_study(seed=7, cases_per_family=16) + assert first == run_axial_query_study(seed=7, cases_per_family=16) + assert first["total_cases"] == 64 + assert first["evidence_class"] == "constructed-controlled-mechanism-not-real-provider" + assert first["anchor_truth_exclusion_count"] == 0 + totals = first["totals"] + assert totals["local-query-gate-then-plugin"]["accepted"] == 64 + assert totals["local-query-gate-then-plugin"]["admitted_with_possible_harm"] == 32 + assert totals["independent-query-intervals"]["accepted"] == 16 + assert totals["shared-orbit-certificate"]["accepted"] == 32 + for arm in ("shared-orbit-certificate", "shared-orbit-plus-one-bounded-anchor"): + assert totals[arm]["sampled_harmful_accepts"] == 0 + assert totals[arm]["admitted_with_possible_harm"] == 0 + assert totals[arm]["exact_fallback_identity_failures"] == 0 + + +@pytest.mark.parametrize("size", [0, -1, True, 1.5]) +def test_invalid_study_size_rejected(size: object) -> None: + with pytest.raises((ValueError, TypeError)): + run_axial_query_study(cases_per_family=size) + + +def test_matching_axes_do_not_authorize_cancellation_of_independent_gauges() -> None: + first = orbit() + independent = AxialRotationOrbit(first.origin, first.axis, "independent-test-gauge") + fallback = first.affine_query([[1, 0, 0]], [[1, 0, 0]], offset=4.0) + candidate = independent.affine_query([[1, 0, 0]], [[1, 0, 0]], offset=3.75) + with pytest.raises(ValueError, match="same shared orbit"): + certify_shared_orbit_advantage( + fallback_loss=fallback, candidate_loss=candidate, scope_admitted=True + ) + + +@pytest.mark.parametrize("identity", ["", " ", None]) +def test_gauge_identity_is_required(identity: object) -> None: + with pytest.raises(ValueError, match="shared_gauge_id"): + AxialRotationOrbit(np.zeros(3), np.array([0.0, 0.0, 1.0]), identity) + + +def test_certificate_record_rejects_contradictory_acceptance() -> None: + accepted = certify_shared_orbit_advantage( + fallback_loss=query(1.0, 0.0, 0.0), + candidate_loss=query(0.0, 0.0, 0.0), + scope_admitted=True, + ) + with pytest.raises(ValueError, match="no positive robust advantage"): + replace(accepted, lower_advantage=-1.0) + with pytest.raises(ValueError, match="reasons disagree"): + replace(accepted, admitted=False) + with pytest.raises(ValueError, match="scope changed"): + replace(accepted, scope="deployment-safety") diff --git a/tests/test_axial_query_certificate_integration.py b/tests/test_axial_query_certificate_integration.py new file mode 100644 index 000000000..678c2a878 --- /dev/null +++ b/tests/test_axial_query_certificate_integration.py @@ -0,0 +1,117 @@ +"""Integration with the existing local query API, not a recreated estimator.""" + +from __future__ import annotations + +from types import SimpleNamespace +from typing import cast + +import numpy as np +import pytest + +from prob4d.axial_query_certificate import ( + AxialRotationOrbit, + certify_shared_orbit_advantage, +) +from prob4d.axial_query_study import analytic_local_control +from prob4d.observable_gauge import ObservableGaugeFactor +from prob4d.query_observability import ( + QueryObservabilityGate, + evaluate_query_observability, + point_position_query_jacobian, +) +from prob4d.sim3 import Sim3 + + +def partial_factor(origin: np.ndarray, axis: np.ndarray) -> ObservableGaugeFactor: + null = np.concatenate(([0.0], axis, np.zeros(3)))[:, None] + projection = np.eye(7) - null @ null.T + values, vectors = np.linalg.eigh(projection) + observed = vectors[:, values > 0.5] + return cast( + ObservableGaugeFactor, + SimpleNamespace( + rank=6, + observable_basis=observed, + nullspace_basis=null, + information_matrix=10.0 * observed @ observed.T, + chart=SimpleNamespace( + linearization=Sim3.identity(), + reference_centroid=origin, + cloud_scale=1.0, + ), + ), + ) + + +def test_actual_local_gate_accepts_a_finite_angle_ambiguous_scalar_query() -> None: + model = AxialRotationOrbit(np.zeros(3), np.array([0.0, 0.0, 1.0]), "shared-test-gauge") + factor = partial_factor(model.origin, model.axis) + point = np.array([1.0, 0.0, 0.0]) + point_jacobian = point_position_query_jacobian(factor, point) + scalar_jacobian = point_jacobian[0:1, :] + report = evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=scalar_jacobian, + ) + gate = QueryObservabilityGate( + minimum_direct_observability_fraction=0.8, + minimum_metric_variance_reduction_fraction=0.8, + maximum_worst_supported_variance_ratio=0.5, + ) + assert report.direct_observability_fraction == pytest.approx(1.0) + assert report.metric_variance_reduction_fraction == pytest.approx(10.0 / 11.0) + assert report.worst_supported_variance_ratio == pytest.approx(1.0 / 11.0) + assert gate.evaluate(report).admitted + # An affine loss, not an arbitrary unbounded nonlinear readout: + # fallback - candidate = 0.25 + cos(theta), whose minimum is -0.75. + fallback = model.affine_query(point[None], np.zeros((1, 3)), offset=4.0) + candidate = model.affine_query(point[None], [[-1.0, 0.0, 0.0]], offset=3.75) + certificate = certify_shared_orbit_advantage( + fallback_loss=fallback, + candidate_loss=candidate, + scope_admitted=True, + ) + assert not certificate.admitted + assert certificate.lower_advantage == pytest.approx(-0.75) + + +def test_study_reference_algebra_matches_existing_api_in_random_frames() -> None: + rng = np.random.default_rng(55) + for _ in range(40): + axis = rng.normal(size=3) + axis /= np.linalg.norm(axis) + model = AxialRotationOrbit(rng.normal(size=3), axis, "shared-test-gauge") + radial = rng.normal(size=3) + radial -= axis * float(axis @ radial) + factor = partial_factor(model.origin, model.axis) + normal = radial / np.linalg.norm(radial) + jacobian = normal[None, :] @ point_position_query_jacobian( + factor, model.origin + radial + ) + actual = evaluate_query_observability( + factor, + prior_covariance_local=np.eye(7), + query_jacobian_local=jacobian, + ) + reference = analytic_local_control(model, radial) + for name, value in reference.items(): + assert getattr(actual, name) == pytest.approx(value, abs=1e-12) + + +@pytest.mark.parametrize("unobserved_variance", [1.0, 1e6, 1e12]) +def test_nullspace_variance_inflation_does_not_fix_zero_query_derivative( + unobserved_variance: float, +) -> None: + model = AxialRotationOrbit(np.zeros(3), np.array([0.0, 0.0, 1.0]), "shared-test-gauge") + factor = partial_factor(model.origin, model.axis) + prior = np.eye(7) + prior[3, 3] = unobserved_variance + jacobian = point_position_query_jacobian(factor, np.array([1.0, 0.0, 0.0]))[0:1] + report = evaluate_query_observability( + factor, prior_covariance_local=prior, query_jacobian_local=jacobian + ) + assert jacobian[0, 3] == 0.0 + assert report.direct_observability_fraction == pytest.approx(1.0) + assert report.prior_query_covariance[0, 0] == pytest.approx(2.0) + assert report.posterior_query_covariance[0, 0] == pytest.approx(2.0 / 11.0)