From 9f3cab5db1dcb09841ce53d448b9ce657379be69 Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Sun, 30 Aug 2026 06:40:53 +0800
Subject: [PATCH] Add finite axial-gauge query bounds and shared-loss regret
 diagnostics

---
 docs/finite-gauge-query-ambiguity.md  | 210 ++++++++++++++++++++
 src/prob4d/axial_gauge_query.py       | 276 ++++++++++++++++++++++++++
 src/prob4d/axial_gauge_query_study.py | 171 ++++++++++++++++
 tests/test_axial_gauge_query.py       | 269 +++++++++++++++++++++++++
 4 files changed, 926 insertions(+)
 create mode 100644 docs/finite-gauge-query-ambiguity.md
 create mode 100644 src/prob4d/axial_gauge_query.py
 create mode 100644 src/prob4d/axial_gauge_query_study.py
 create mode 100644 tests/test_axial_gauge_query.py

diff --git a/docs/finite-gauge-query-ambiguity.md b/docs/finite-gauge-query-ambiguity.md
new file mode 100644
index 000000000..2655ced5b
--- /dev/null
+++ b/docs/finite-gauge-query-ambiguity.md
@@ -0,0 +1,210 @@
+# Finite axial-gauge ambiguity for physical queries
+
+## Purpose and status
+
+This experimental, NumPy-only analysis extends the **local**
+[query-observability diagnostic](query-conditioned-observability.md) with exact
+finite-orbit bounds for a restricted but useful geometry: one free rotation
+about a registered line per shared gauge group, with all other transformation
+parameters fixed. It also supplies sharp worst-case regret for actual affine
+action losses that share those same gauges.
+
+The existing local diagnostic is not mathematically incorrect. The new control
+shows why its direct-observability fraction must not be interpreted as a global
+invariance or decision certificate. A small first derivative does not establish
+that a query is constant on the finite observation-equivalence class.
+
+This is analytic development evidence, not real-provider competence,
+calibration, BayesianPhysTwin benefit, Causal4D benefit, or deployment safety.
+No provider, dataset, protected target, or physical execution is used. No
+existing result, production exporter, stable API, or admission rule is changed.
+
+## A falsifying control with a positive counterpart
+
+An exactly collinear overlap lies on the x-axis. Rotation about that axis leaves
+all its point coordinates unchanged. A probe initially at `(0, 0.1, 0)` has
+
+$$q(\phi)=0.1\cos\phi,\qquad q'(0)=0,\qquad q''(0)=-0.1.$$
+
+For the existing analytic rank-six factor, the scalar y-position Jacobian has
+zero projection onto the missing rotation-x direction. With its complete
+identity prior and observable precision 10, the local report gives direct
+support 1, variance reduction 10/11, and worst variance ratio 1/11. The existing
+illustrative `(0.8, 0.8, 0.5)` gate admits this scalar query. Nevertheless, the
+finite y-query range over a full rotation is `[-100, 100] mm`.
+
+These comparisons concern **different guarantees**. A full-circle ambiguity
+bound is deliberately stronger than a prior-weighted local uncertainty
+calculation. The control does not establish that the local posterior's actual
+credible interval undercovers or that its Bayesian decision is wrong.
+
+There are also useful positive controls. The axial coordinate of an off-axis
+probe remains fixed. Two points with the same transverse displacement have an
+exactly invariant y-difference when they share one angle. Treating their angles
+as separately variable destroys that cancellation; falsely merging unrelated
+angles invents it.
+
+## Domain assumptions
+
+The caller registers metric-space points, their persistent identities, a pivot
+and unit axis for each group, and which points share exactly the same angle.
+Each group angle ranges over a full circle; the joint domain is the Cartesian
+product of those circles. This is a set-valued domain, **not** an assertion of
+statistical independence. A smaller or coupled domain is conservatively
+contained by this product, but the bounds need not be sharp on that smaller
+domain.
+
+A rank-six information matrix alone does not prove the presence of an exact
+axial symmetry. The caller must establish that the finite transformations are
+admissible for the relevant observation model. An exactly collinear support
+has zero geometric motion under its axial orbit. A near-collinear support does
+not: `maximum_support_motion` returns twice its maximum distance to the axis.
+That is a geometric motion bound, not a bound on likelihood ratio, calibration,
+or posterior probability. Additional observations, metric anchors, covariance
+orientation, priors, or physical constraints can restrict or break the proposed
+symmetry and must not be ignored.
+
+Points are already expressed in the aligned metric frame. Uncertainty in scale,
+axis location, other gauge directions, identity, deformable dynamics, and
+nonlinear physical queries is outside this exact special case. The helper does
+not automatically turn a provider dependence-group label into a shared angle.
+
+## Exact query bounds
+
+Let `a_g` be a unit axis through `c_g`, and decompose each point in group `g` as
+
+$$d_i=p_i-c_g,\quad d_i^\parallel=a_ga_g^\top d_i,
+\quad d_i^\perp=d_i-d_i^\parallel.$$
+
+Rodrigues' formula gives
+
+$$p_i(\phi_g)=c_g+d_i^\parallel
+ +d_i^\perp\cos\phi_g+(a_g\times d_i^\perp)\sin\phi_g.$$
+
+For affine queries `q_k = b_k + sum_i w_ki^T p_i`, collect contributions
+**inside each shared group before taking an amplitude**:
+
+$$q_k(\phi)=C_k+\sum_g[A_{kg}\cos\phi_g+B_{kg}\sin\phi_g].$$
+
+Consequently the sharp componentwise interval is
+
+$$q_k\in\left[C_k-\sum_g\sqrt{A_{kg}^2+B_{kg}^2},
+                 C_k+\sum_g\sqrt{A_{kg}^2+B_{kg}^2}\right].$$
+
+Proof: each harmonic has extrema plus or minus its amplitude, attained at
+`atan2(B_kg, A_kg)` and that angle plus pi. Product-domain freedom lets the
+extrema be attained simultaneously for one query. Different query endpoints
+need not be jointly attainable: the interval collection is not an exact
+rectangular joint feasible set.
+
+Global invariance holds exactly when every `A_kg` and `B_kg` is zero. Stationarity
+at the reference only requires every `B_kg` to be zero. This isolates the
+first-order blind spot without a Monte Carlo tolerance or an arbitrary grid.
+
+## Sharp action regret preserves common uncertainty
+
+When the query family contains actual lower-is-better affine action losses
+`L_k`, define regret for a caller-selected action `k` by
+
+$$R_k=\sup_\phi\left[L_k(\phi)-\min_j L_j(\phi)\right].$$
+
+For finitely many actions this has the closed form
+
+$$R_k=\max_j\left\{C_k-C_j+
+ \sum_g\sqrt{(A_{kg}-A_{jg})^2+(B_{kg}-B_{jg})^2}\right\}.$$
+
+Proof: interchange the supremum and the finite maximum over `j`, and apply the
+harmonic bound to each loss difference. Including `j=k` ensures nonnegative
+regret. `regret_witness(k)` returns a competing action and an angle vector
+attaining this bound. Query contrasts share the same angles, so amplitudes are
+computed **after subtraction**, not from separate marginal interval endpoints.
+
+For `L_0=0.1-y` and `L_1=0.1+y`, the reference selects action 0 with zero loss.
+At half a turn, action 0 has loss 200 mm while action 1 has zero loss. Both
+worst-case regrets are 200 mm, so an illustrative 50 mm regret budget rejects
+that reference-selected action. This is a constructed loss, not robot-control
+or physical-state estimation performance.
+
+Conversely, `L_0=0.05+y` and `L_1=0.10+y` have overlapping marginal loss
+intervals but action 0 is uniformly preferable: its exact regret is zero.
+The common uncertain term cancels. The method therefore does not simply
+reject everything with a broad marginal interval.
+
+## Reproduce and use
+
+From an installed source checkout:
+
+```bash
+python -m pytest tests/test_axial_gauge_query.py
+python -m prob4d.axial_gauge_query_study \
+  --output outputs/axial-gauge-query-development-v1.json
+```
+
+The output is labelled deterministic development evidence. Repeating the command
+at the same path is allowed only for identical bytes; different existing output
+is never overwritten. The tests include an independent quaternion-rotation
+reference, attaining extrema, common-nuisance cancellation, frame and axis-sign
+invariance, malformed inputs, and parity with the actual existing local gate.
+The integration test must run in a complete checkout; it does not silently skip
+when the existing modules cannot be imported. Random algebraic test cases are
+development controls, not independent empirical trials.
+
+```python
+import numpy as np
+from prob4d.axial_gauge_query import AxialGaugeOrbit, affine_axial_queries
+
+orbit = AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+losses = affine_axial_queries(
+    np.array([[0.0, 0.1, 0.0]]),
+    np.array([[[0.0, -1.0, 0.0]], [[0.0, 1.0, 0.0]]]),
+    offsets=np.array([0.1, 0.1]),
+    point_group_ids=("window-gauge-0",),
+    orbits={"window-gauge-0": orbit},
+)
+regret = losses.worst_case_regrets()  # [0.2, 0.2], in these loss units
+within_budget = losses.within_regret_budget(0, maximum_regret=0.05)  # False
+```
+
+`within_regret_budget` executes no action and changes no belief. A downstream caller
+must retain every existing provider, calibration, identifiability, and support
+gate, then route rejection through its existing complete-belief fallback.
+Passing this additional check cannot rescue an upstream failure. Applying a
+position bound to a nonlinear simulator loss without a valid reduction does not
+produce an action-regret guarantee.
+
+The formulas are analytic, but the implementation evaluates them in float64,
+not certified interval arithmetic. A caller-frozen nonnegative numerical margin
+can make a threshold decision more conservative; it is not a proof of a
+floating-point enclosure. The implementation rejects nonfinite values rather
+than manufacturing a bounded result.
+
+## Relation to existing work and paper ownership
+
+Weak-direction-aware registration is established prior art: see Tuna et al.,
+[X-ICP](https://arxiv.org/abs/2211.16335). Group-invariant observable/unobservable
+decompositions are also established; see Shen and Leok,
+[Geometric Symmetry Reduction of the Unobservable Subspace for Kalman
+Filtering](https://arxiv.org/abs/1901.03474). Rodrigues' formula, harmonic support
+functions, and worst-case regret are not claimed as new generic mathematics.
+This note does not establish an exhaustive novelty claim.
+
+The bounded contribution candidate is the finite-gauge, shared-lineage query
+and decision audit for partially observable learned 4D windows, together with
+its stationary-derivative counterexample and cancellation-preserving positive
+controls. The existing [Gaussian linearization-closure
+diagnostic](gauge-linearization-closure.md) instead checks numerical moment
+closure under a declared joint Gaussian. Neither result subsumes the other's
+assumptions or gives empirical calibration.
+
+This belongs to the Prob4D geometric observation boundary. It does not duplicate
+the general query-quotient, Jeffrey-lift, or cross-intervention theorems owned by
+the existing theory companion in the paper repository. No fourth manuscript or
+new physical acquisition is proposed. Paper-facing development results and
+interpretation belong in `FlorianPfaff/BayesianPhysTwin-Paper`; this public
+repository owns the method, tests, and reproducer.
+
+Real-provider value remains a separate empirical question. The existing
+PointWorld/Flat'n'Fold source-qualification route in issue #333 remains governed
+by its own unresolved source requirements and target-access rules. This method
+neither authorizes its execution nor makes it a dependency for the already
+bounded BayesianPhysTwin and Causal4D manuscripts.
diff --git a/src/prob4d/axial_gauge_query.py b/src/prob4d/axial_gauge_query.py
new file mode 100644
index 000000000..00d38c17c
--- /dev/null
+++ b/src/prob4d/axial_gauge_query.py
@@ -0,0 +1,276 @@
+"""Exact finite-orbit bounds for affine queries of axial-gauge point families.
+
+Experimental, NumPy-only analysis. This is not a posterior or an admission
+certificate for a provider. The caller must register the axial symmetry,
+point identities, shared gauge groups, and the domain of allowed angles.
+A small local Jacobian or a rank-six factor does not establish that domain.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+
+import numpy as np
+from numpy.typing import NDArray
+
+FloatArray = NDArray[np.float64]
+
+
+def _finite_array(value: object, *, name: str) -> FloatArray:
+    result = np.asarray(value, dtype=np.float64).copy()
+    if not np.all(np.isfinite(result)):
+        raise ValueError(f"{name} must contain only finite values")
+    result.setflags(write=False)
+    return result
+
+
+def _vector(value: object, *, name: str) -> FloatArray:
+    result = _finite_array(value, name=name)
+    if result.shape != (3,):
+        raise ValueError(f"{name} must have shape (3,)")
+    return result
+
+
+def _points(value: object, *, name: str) -> FloatArray:
+    result = _finite_array(value, name=name)
+    if result.ndim != 2 or result.shape[0] < 1 or result.shape[1] != 3:
+        raise ValueError(f"{name} must have shape (N, 3), N >= 1")
+    return result
+
+
+def _group_id(value: object) -> str:
+    if not isinstance(value, str) or not value or value != value.strip():
+        raise ValueError("gauge group IDs must be nonempty, unpadded strings")
+    return value
+
+
+def _nonnegative(value: float, *, name: str) -> float:
+    result = float(value)
+    if not np.isfinite(result) or result < 0.0:
+        raise ValueError(f"{name} must be finite and nonnegative")
+    return result
+
+
+@dataclass(frozen=True)
+class AxialGaugeOrbit:
+    """All rotations about a registered metric-space line, with fixed scale.
+
+    ``axis`` is normalized on construction. ``pivot`` lies on the line.
+    The class defines an orbit; it does not infer a symmetry from observations.
+    """
+
+    pivot: FloatArray
+    axis: FloatArray
+
+    def __post_init__(self) -> None:
+        pivot = _vector(self.pivot, name="pivot")
+        axis = _vector(self.axis, name="axis")
+        scale = float(np.max(np.abs(axis)))
+        if scale == 0.0:
+            raise ValueError("axis must be nonzero")
+        normalized = axis / scale
+        normalized = normalized / np.linalg.norm(normalized)
+        normalized.setflags(write=False)
+        object.__setattr__(self, "pivot", pivot)
+        object.__setattr__(self, "axis", normalized)
+
+    def components(self, points: FloatArray) -> tuple[FloatArray, FloatArray, FloatArray]:
+        """Return constant, cosine, and sine point components (Rodrigues)."""
+        values = _points(points, name="points")
+        relative = values - self.pivot
+        parallel = (relative @ self.axis)[:, None] * self.axis
+        constant = self.pivot + parallel
+        cosine = relative - parallel
+        sine = np.cross(self.axis, cosine)
+        return (
+            _finite_array(constant, name="constant orbit component"),
+            _finite_array(cosine, name="cosine orbit component"),
+            _finite_array(sine, name="sine orbit component"),
+        )
+
+    def transform(self, points: FloatArray, angle: float) -> FloatArray:
+        """Apply a finite rotation, without a local Gaussian approximation."""
+        angle = float(angle)
+        if not np.isfinite(angle):
+            raise ValueError("angle must be finite")
+        constant, cosine, sine = self.components(points)
+        return _finite_array(
+            constant + np.cos(angle) * cosine + np.sin(angle) * sine,
+            name="transformed points",
+        )
+
+    def maximum_support_motion(self, support_points: FloatArray) -> float:
+        """Sharp maximum point displacement over the full orbit, in input units.
+
+        This is twice the largest distance to the axis. Zero means geometric
+        invariance of the supplied support. A small nonzero value is only a
+        geometric bound: it is NOT equal likelihood or statistical equivalence.
+        """
+        _, cosine, _ = self.components(support_points)
+        bound = float(2.0 * np.max(np.linalg.norm(cosine, axis=1)))
+        if not np.isfinite(bound):
+            raise ValueError("support motion is not finite")
+        return bound
+
+
+@dataclass(frozen=True)
+class AxialQueryFamily:
+    """q_k(phi) = constant_k + sum_g (cosine_kg cos phi_g + sine_kg sin phi_g).
+
+    All queries share the SAME angle for each group. Bounds are sharp when
+    group angles range freely over a product of full circles. For a smaller
+    coupled angle domain they remain conservative, not necessarily sharp.
+    """
+
+    group_ids: tuple[str, ...]
+    constant: FloatArray
+    cosine: FloatArray
+    sine: FloatArray
+
+    def __post_init__(self) -> None:
+        if isinstance(self.group_ids, (str, bytes)):
+            raise ValueError("group_ids must be a sequence of group IDs")
+        groups = tuple(_group_id(group) for group in self.group_ids)
+        if len(groups) != len(set(groups)):
+            raise ValueError("group_ids must be unique")
+        constant = _finite_array(self.constant, name="constant")
+        if constant.ndim != 1 or constant.size < 1:
+            raise ValueError("constant must have shape (Q,), Q >= 1")
+        cosine = _finite_array(self.cosine, name="cosine")
+        sine = _finite_array(self.sine, name="sine")
+        expected = (constant.size, len(groups))
+        if cosine.shape != expected or sine.shape != expected:
+            raise ValueError("cosine and sine must have shape (Q, number of groups)")
+        object.__setattr__(self, "group_ids", groups)
+        object.__setattr__(self, "constant", constant)
+        object.__setattr__(self, "cosine", cosine)
+        object.__setattr__(self, "sine", sine)
+
+    def evaluate(self, angles: FloatArray) -> FloatArray:
+        angles = _finite_array(angles, name="angles")
+        if angles.shape != (len(self.group_ids),):
+            raise ValueError("angles must have shape (number of groups,)")
+        return _finite_array(
+            self.constant + self.cosine @ np.cos(angles) + self.sine @ np.sin(angles),
+            name="query values",
+        )
+
+    def bounds(self) -> tuple[FloatArray, FloatArray]:
+        """Sharp componentwise intervals; not a joint rectangular feasible set."""
+        radius = np.sum(np.hypot(self.cosine, self.sine), axis=1)
+        return (
+            _finite_array(self.constant - radius, name="lower bounds"),
+            _finite_array(self.constant + radius, name="upper bounds"),
+        )
+
+    def contrast_bounds(self) -> tuple[FloatArray, FloatArray]:
+        """Sharp bounds on q_i - q_j, preserving shared-angle cancellation."""
+        center = self.constant[:, None] - self.constant[None, :]
+        cosine = self.cosine[:, None, :] - self.cosine[None, :, :]
+        sine = self.sine[:, None, :] - self.sine[None, :, :]
+        radius = np.sum(np.hypot(cosine, sine), axis=2)
+        return (
+            _finite_array(center - radius, name="lower contrast bounds"),
+            _finite_array(center + radius, name="upper contrast bounds"),
+        )
+
+    def worst_case_regrets(self) -> FloatArray:
+        """Sharp regrets when queries are actual lower-is-better action losses.
+
+        The caller must supply those loss functions. This does not convert
+        position uncertainty into a simulator-loss or robot-safety guarantee.
+        """
+        _, upper = self.contrast_bounds()
+        return _finite_array(np.max(upper, axis=1), name="worst-case regrets")
+
+    def _action_index(self, action: int) -> int:
+        if isinstance(action, bool) or not isinstance(action, (int, np.integer)):
+            raise TypeError("action must be an integer index")
+        action = int(action)
+        if not 0 <= action < self.constant.size:
+            raise ValueError("action is outside the query family")
+        return action
+
+    def regret_witness(self, action: int) -> tuple[int, FloatArray]:
+        """Return a competing action and angle vector attaining the worst regret."""
+        action = self._action_index(action)
+        _, upper = self.contrast_bounds()
+        competitor = int(np.argmax(upper[action]))
+        angles = np.arctan2(
+            self.sine[action] - self.sine[competitor],
+            self.cosine[action] - self.cosine[competitor],
+        )
+        return competitor, _finite_array(angles, name="witness angles")
+
+    def within_regret_budget(
+        self,
+        action: int,
+        *,
+        maximum_regret: float,
+        numerical_margin: float = 0.0,
+    ) -> bool:
+        """Check a caller-selected action against a caller-frozen regret budget.
+
+        No action is executed and no belief is modified. Rejection must be
+        routed to the caller's existing complete-belief fallback. These are
+        float64 evaluations of analytic bounds, not interval-arithmetic proofs.
+        """
+        action = self._action_index(action)
+        budget = _nonnegative(maximum_regret, name="maximum_regret")
+        margin = _nonnegative(numerical_margin, name="numerical_margin")
+        return bool(float(self.worst_case_regrets()[action]) + margin <= budget)
+
+
+def affine_axial_queries(
+    points: FloatArray,
+    weights: FloatArray,
+    *,
+    point_group_ids: Sequence[str],
+    orbits: Mapping[str, AxialGaugeOrbit],
+    offsets: FloatArray | None = None,
+) -> AxialQueryFamily:
+    """Build affine queries ``offsets[k] + sum_i weights[k,i] dot points_i``.
+
+    Points are already in the aligned metric frame. Their group labels encode
+    a declared shared finite rotation, not merely similar geometry or an
+    arbitrary dependence label. Contributions are summed INSIDE each group
+    before an amplitude is computed. Incorrect group merging can fabricate
+    cancellation, just as incorrect splitting can discard real cancellation.
+    """
+    values = _points(points, name="points")
+    weights = _finite_array(weights, name="weights")
+    if weights.ndim != 3 or weights.shape[0] < 1 or weights.shape[1:] != values.shape:
+        raise ValueError("weights must have shape (Q, N, 3), Q >= 1")
+    if isinstance(point_group_ids, (str, bytes)):
+        raise ValueError("point_group_ids must be a sequence of group IDs")
+    labels = tuple(_group_id(group) for group in point_group_ids)
+    if len(labels) != values.shape[0]:
+        raise ValueError("one gauge group ID is required per point")
+    keys = tuple(_group_id(group) for group in orbits)
+    if set(keys) != set(labels):
+        raise ValueError("orbits must contain exactly the referenced gauge groups")
+    if any(not isinstance(orbit, AxialGaugeOrbit) for orbit in orbits.values()):
+        raise TypeError("every orbit must be an AxialGaugeOrbit")
+    groups = tuple(sorted(keys))
+    query_count = weights.shape[0]
+    constant = (
+        np.zeros(query_count, dtype=np.float64)
+        if offsets is None
+        else _finite_array(offsets, name="offsets").copy()
+    )
+    if constant.shape != (query_count,):
+        raise ValueError("offsets must have shape (Q,)")
+    cosine = np.zeros((query_count, len(groups)), dtype=np.float64)
+    sine = np.zeros_like(cosine)
+    for index, group in enumerate(groups):
+        mask = np.array([label == group for label in labels])
+        base, cos_part, sin_part = orbits[group].components(values[mask])
+        group_weights = weights[:, mask, :]
+        constant += np.einsum("qnc,nc->q", group_weights, base)
+        cosine[:, index] = np.einsum("qnc,nc->q", group_weights, cos_part)
+        sine[:, index] = np.einsum("qnc,nc->q", group_weights, sin_part)
+    return AxialQueryFamily(groups, constant, cosine, sine)
+
+
+__all__ = ["AxialGaugeOrbit", "AxialQueryFamily", "affine_axial_queries"]
diff --git a/src/prob4d/axial_gauge_query_study.py b/src/prob4d/axial_gauge_query_study.py
new file mode 100644
index 000000000..f96bab2fb
--- /dev/null
+++ b/src/prob4d/axial_gauge_query_study.py
@@ -0,0 +1,171 @@
+"""Deterministic development controls for finite axial-gauge query ambiguity.
+
+No provider, dataset, posterior, held-out target, or physical execution is used.
+The local reference matches the rank-six analytic factor already used by the
+query-observability study; a separate integration test checks that parity.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from .axial_gauge_query import AxialGaugeOrbit, AxialQueryFamily, affine_axial_queries
+
+
+def run_axial_gauge_query_study() -> dict[str, Any]:
+    orbit = AxialGaugeOrbit(pivot=np.zeros(3), axis=np.array([1.0, 0.0, 0.0]))
+    support = np.column_stack((np.linspace(-1.0, 1.0, 48), np.zeros((48, 2))))
+    point = np.array([[0.0, 0.1, 0.0]])
+    radial = affine_axial_queries(
+        point,
+        np.array([[[0.0, 1.0, 0.0]]]),
+        point_group_ids=("shared",),
+        orbits={"shared": orbit},
+    )
+    axial = affine_axial_queries(
+        np.array([[0.04, 0.1, 0.0]]),
+        np.array([[[1.0, 0.0, 0.0]]]),
+        point_group_ids=("shared",),
+        orbits={"shared": orbit},
+    )
+    points = np.array([[0.0, 0.1, 0.0], [1.0, 0.1, 0.0]])
+    weights = np.array([[[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]]])
+    common = affine_axial_queries(
+        points,
+        weights,
+        point_group_ids=("shared", "shared"),
+        orbits={"shared": orbit},
+    )
+    separate = affine_axial_queries(
+        points,
+        weights,
+        point_group_ids=("first", "second"),
+        orbits={"first": orbit, "second": orbit},
+    )
+    losses = affine_axial_queries(
+        point,
+        np.array([[[0.0, -1.0, 0.0]], [[0.0, 1.0, 0.0]]]),
+        offsets=np.array([0.1, 0.1]),
+        point_group_ids=("shared",),
+        orbits={"shared": orbit},
+    )
+    common_losses = affine_axial_queries(
+        point,
+        np.array([[[0.0, 1.0, 0.0]], [[0.0, 1.0, 0.0]]]),
+        offsets=np.array([0.05, 0.10]),
+        point_group_ids=("shared",),
+        orbits={"shared": orbit},
+    )
+    # Independent closed-form local reference: J_y=[.1,0,0,0,0,1,0].
+    # Only rotation-x has zero precision. The complete prior is I_7.
+    jacobian = np.array([[0.1, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]])
+    precision = np.diag([10.0, 0.0, 10.0, 10.0, 10.0, 10.0, 10.0])
+    posterior = np.linalg.solve(np.eye(7) + precision, np.eye(7))
+    prior_variance = float((jacobian @ jacobian.T)[0, 0])
+    posterior_variance = float((jacobian @ posterior @ jacobian.T)[0, 0])
+    variance_reduction = 1.0 - posterior_variance / prior_variance
+    worst_ratio = posterior_variance / prior_variance
+    local_admitted = variance_reduction >= 0.8 and worst_ratio <= 0.5
+    competitor, angles = losses.regret_witness(0)
+    at_witness = losses.evaluate(angles)
+
+    def interval(family: AxialQueryFamily) -> list[float]:
+        lower, upper = family.bounds()
+        return [float(lower[0]), float(upper[0])]
+
+    return {
+        "schema_name": "prob4d.axial-gauge-query-development-control",
+        "schema_version": 1,
+        "classification": "deterministic analytic development evidence",
+        "claim_boundary": (
+            "Exact only over the declared axial-orbit domain and affine queries; "
+            "no provider competence, posterior calibration, physical validation, "
+            "BayesianPhysTwin benefit, Causal4D benefit, or deployment safety."
+        ),
+        "geometry": {
+            "axis": [1.0, 0.0, 0.0],
+            "pivot": [0.0, 0.0, 0.0],
+            "support_points": 48,
+            "probe_radius_m": 0.1,
+            "maximum_support_motion_m": orbit.maximum_support_motion(support),
+        },
+        "stationary_derivative_counterexample": {
+            "query": "y-coordinate of a 100 mm off-axis point",
+            "query_at_reference_m": float(radial.evaluate(np.zeros(1))[0]),
+            "first_twist_derivative_m_per_rad": float(radial.sine[0, 0]),
+            "second_twist_derivative_m_per_rad2": -float(radial.cosine[0, 0]),
+            "finite_orbit_interval_m": interval(radial),
+            "local_linear_reference": {
+                "factor_rank": 6,
+                "direct_observability_fraction": 1.0,
+                "metric_variance_reduction_fraction": variance_reduction,
+                "worst_supported_variance_ratio": worst_ratio,
+                "gate_thresholds": [0.8, 0.8, 0.5],
+                "local_gate_admits": local_admitted,
+            },
+        },
+        "positive_control": {
+            "query": "axial coordinate of an off-axis point",
+            "finite_orbit_interval_m": interval(axial),
+        },
+        "lineage_control": {
+            "query": "difference of equal-radius y-coordinates",
+            "shared_angle_interval_m": interval(common),
+            "separately_variable_angles_interval_m": interval(separate),
+            "boundary": "Group equality is declared, never inferred from coincident axes.",
+        },
+        "action_control": {
+            "losses": ["0.1 - y", "0.1 + y"],
+            "nominal_action": 0,
+            "nominal_losses_m": losses.evaluate(np.zeros(1)).tolist(),
+            "worst_case_regrets_m": losses.worst_case_regrets().tolist(),
+            "illustrative_regret_budget_m": 0.05,
+            "within_budget": losses.within_regret_budget(0, maximum_regret=0.05),
+            "witness_competitor": competitor,
+            "witness_angle_rad": angles.tolist(),
+            "witness_losses_m": at_witness.tolist(),
+        },
+        "shared_action_nuisance_control": {
+            "losses": ["0.05 + y", "0.10 + y"],
+            "marginal_intervals_m": np.column_stack(common_losses.bounds()).tolist(),
+            "worst_case_regrets_m": common_losses.worst_case_regrets().tolist(),
+            "action_zero_within_zero_budget": common_losses.within_regret_budget(
+                0, maximum_regret=0.0
+            ),
+        },
+        "information_boundary": {
+            "provider_forward_calls": 0,
+            "dataset_records_accessed": 0,
+            "protected_targets_opened": 0,
+            "physical_executions": 0,
+        },
+    }
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--output", type=Path)
+    args = parser.parse_args()
+    text = json.dumps(
+        run_axial_gauge_query_study(), indent=2, sort_keys=True, allow_nan=False
+    ) + "\n"
+    if args.output is None:
+        print(text, end="")
+    else:
+        args.output.parent.mkdir(parents=True, exist_ok=True)
+        # Never overwrite a different retained result.
+        if args.output.exists():
+            if args.output.read_text(encoding="utf-8") != text:
+                raise FileExistsError(f"refusing to replace different evidence: {args.output}")
+        else:
+            with args.output.open("x", encoding="utf-8") as stream:
+                stream.write(text)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/test_axial_gauge_query.py b/tests/test_axial_gauge_query.py
new file mode 100644
index 000000000..6f4cdca7a
--- /dev/null
+++ b/tests/test_axial_gauge_query.py
@@ -0,0 +1,269 @@
+"""Analytic and adversarial controls, with no provider or data access."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from prob4d.axial_gauge_query import AxialGaugeOrbit, AxialQueryFamily, affine_axial_queries
+from prob4d.axial_gauge_query_study import run_axial_gauge_query_study
+
+
+def _orbit() -> AxialGaugeOrbit:
+    return AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+
+
+def _family(points: np.ndarray, weights: np.ndarray) -> AxialQueryFamily:
+    return affine_axial_queries(
+        points, weights, point_group_ids=("g",) * len(points), orbits={"g": _orbit()}
+    )
+
+
+def _quaternion_rotation(axis: np.ndarray, angle: float) -> np.ndarray:
+    """Independent finite-rotation reference, not the production decomposition."""
+    axis = axis / np.linalg.norm(axis)
+    x, y, z = np.sin(angle / 2.0) * axis
+    w = np.cos(angle / 2.0)
+    return np.array(
+        [
+            [1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w)],
+            [2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w)],
+            [2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)],
+        ]
+    )
+
+
+def test_zero_first_derivative_does_not_mean_finite_invariance() -> None:
+    family = _family(np.array([[0.0, 0.1, 0.0]]), np.array([[[0.0, 1.0, 0.0]]]))
+    assert family.sine[0, 0] == 0.0
+    assert family.cosine[0, 0] == 0.1
+    np.testing.assert_allclose(family.bounds(), [[-0.1], [0.1]])
+    np.testing.assert_allclose(family.evaluate(np.array([np.pi])), [-0.1])
+    support = np.array([[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+    assert _orbit().maximum_support_motion(support) == 0.0
+    np.testing.assert_array_equal(_orbit().transform(support, np.pi), support)
+
+
+@pytest.mark.parametrize("point", [[0.04, 0.1, 0.0], [0.04, 0.0, 0.0]])
+def test_axial_query_is_globally_invariant(point: list[float]) -> None:
+    family = _family(np.array([point]), np.array([[[1.0, 0.0, 0.0]]]))
+    np.testing.assert_allclose(family.bounds(), [[0.04], [0.04]])
+
+
+def test_shared_gauge_cancellation_and_invalid_merging_control() -> None:
+    points = np.array([[0.0, 0.1, 0.0], [1.0, 0.1, 0.0]])
+    weights = np.array([[[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]]])
+    shared = _family(points, weights)
+    separate = affine_axial_queries(
+        points,
+        weights,
+        point_group_ids=("a", "b"),
+        orbits={"a": _orbit(), "b": _orbit()},
+    )
+    np.testing.assert_array_equal(shared.bounds(), [[0.0], [0.0]])
+    np.testing.assert_allclose(separate.bounds(), [[-0.2], [0.2]])
+    # Coincident axes do not justify merging two separately variable gauges.
+    np.testing.assert_allclose(separate.evaluate(np.array([0.0, np.pi])), [0.2])
+
+
+def test_action_regret_has_an_attaining_witness() -> None:
+    family = affine_axial_queries(
+        np.array([[0.0, 0.1, 0.0]]),
+        np.array([[[0.0, -1.0, 0.0]], [[0.0, 1.0, 0.0]]]),
+        offsets=np.array([0.1, 0.1]),
+        point_group_ids=("g",),
+        orbits={"g": _orbit()},
+    )
+    np.testing.assert_allclose(family.evaluate(np.zeros(1)), [0.0, 0.2])
+    np.testing.assert_allclose(family.worst_case_regrets(), [0.2, 0.2])
+    for action in range(2):
+        competitor, angles = family.regret_witness(action)
+        losses = family.evaluate(angles)
+        assert losses[action] - np.min(losses) == pytest.approx(0.2)
+        assert losses[action] - losses[competitor] == pytest.approx(0.2)
+    assert not family.within_regret_budget(0, maximum_regret=0.05)
+    assert family.within_regret_budget(0, maximum_regret=0.2)
+    assert not family.within_regret_budget(0, maximum_regret=0.2, numerical_margin=1e-9)
+
+
+def test_shared_action_nuisance_cancels_before_bounding() -> None:
+    family = AxialQueryFamily(
+        ("g",), np.array([0.05, 0.10]), np.array([[0.1], [0.1]]), np.zeros((2, 1))
+    )
+    lower, upper = family.bounds()
+    assert upper[0] > lower[1]  # Marginal intervals overlap.
+    np.testing.assert_allclose(family.contrast_bounds()[1][0, 1], -0.05)
+    np.testing.assert_allclose(family.worst_case_regrets(), [0.0, 0.05])
+    assert family.within_regret_budget(0, maximum_regret=0.0)
+
+
+def test_multigroup_affine_bounds_and_regrets_against_finite_rotations() -> None:
+    rng = np.random.default_rng(8302026)
+    for _ in range(20):
+        points = rng.normal(size=(6, 3))
+        weights = rng.normal(size=(4, 6, 3))
+        offsets = rng.normal(size=4)
+        labels = ("a", "b", "a", "b", "a", "b")
+        orbits = {
+            "a": AxialGaugeOrbit(rng.normal(size=3), rng.normal(size=3)),
+            "b": AxialGaugeOrbit(rng.normal(size=3), rng.normal(size=3)),
+        }
+        family = affine_axial_queries(
+            points, weights, offsets=offsets, point_group_ids=labels, orbits=orbits
+        )
+        angles = rng.uniform(-np.pi, np.pi, size=2)
+        transformed = np.empty_like(points)
+        for index, group in enumerate(labels):
+            orbit = orbits[group]
+            rotation = _quaternion_rotation(orbit.axis, angles[family.group_ids.index(group)])
+            transformed[index] = orbit.pivot + rotation @ (points[index] - orbit.pivot)
+        reference = offsets + np.einsum("qnc,nc->q", weights, transformed)
+        np.testing.assert_allclose(family.evaluate(angles), reference, atol=1e-12)
+        lower, upper = family.bounds()
+        assert np.all(reference >= lower - 1e-12)
+        assert np.all(reference <= upper + 1e-12)
+        for query in range(4):
+            maximum_angles = np.arctan2(family.sine[query], family.cosine[query])
+            assert family.evaluate(maximum_angles)[query] == pytest.approx(upper[query])
+            assert family.evaluate(maximum_angles + np.pi)[query] == pytest.approx(lower[query])
+            competitor, witness = family.regret_witness(query)
+            losses = family.evaluate(witness)
+            regret = family.worst_case_regrets()[query]
+            assert losses[query] - losses[competitor] == pytest.approx(regret, abs=1e-12)
+            assert losses[query] - min(losses) == pytest.approx(regret, abs=1e-12)
+
+
+def test_rigid_frame_and_axis_sign_invariance() -> None:
+    points = np.array([[0.5, 0.2, -0.1], [-0.4, 0.8, 0.3]])
+    weights = np.array([[[0.1, 1.0, -0.3], [0.2, -0.2, 0.7]]])
+    original = _family(points, weights)
+    rotation = _quaternion_rotation(np.array([1.0, 2.0, -1.0]), 1.2)
+    shift = np.array([2.0, -0.4, 0.7])
+    moved_points = points @ rotation.T + shift
+    moved_weights = weights @ rotation.T
+    offsets = -np.einsum("qnc,c->q", moved_weights, shift)
+    moved = affine_axial_queries(
+        moved_points,
+        moved_weights,
+        offsets=offsets,
+        point_group_ids=("g", "g"),
+        orbits={"g": AxialGaugeOrbit(shift, rotation @ np.array([1.0, 0.0, 0.0]))},
+    )
+    np.testing.assert_allclose(original.bounds(), moved.bounds(), atol=1e-12)
+    flipped = affine_axial_queries(
+        points,
+        weights,
+        point_group_ids=("g", "g"),
+        orbits={"g": AxialGaugeOrbit(np.zeros(3), np.array([-1.0, 0.0, 0.0]))},
+    )
+    np.testing.assert_allclose(original.bounds(), flipped.bounds())
+    np.testing.assert_allclose(
+        original.evaluate(np.array([0.7])), flipped.evaluate(np.array([-0.7]))
+    )
+
+
+def test_near_line_is_not_silently_called_an_exact_symmetry() -> None:
+    support = np.array([[-1.0, 0.0, 0.0], [1.0, 0.002, 0.0]])
+    assert _orbit().maximum_support_motion(support) == pytest.approx(0.004)
+    assert _orbit().maximum_support_motion(support) > 0.001
+
+
+def test_constant_family_without_free_angles() -> None:
+    family = AxialQueryFamily((), np.array([1.0, 2.0]), np.empty((2, 0)), np.empty((2, 0)))
+    np.testing.assert_allclose(family.bounds(), [[1.0, 2.0], [1.0, 2.0]])
+    np.testing.assert_allclose(family.worst_case_regrets(), [0.0, 1.0])
+    assert family.regret_witness(1)[1].size == 0
+
+
+def test_readonly_copies_and_determinism() -> None:
+    axis = np.array([10.0, 0.0, 0.0])
+    orbit = AxialGaugeOrbit(np.zeros(3), axis)
+    axis[0] = 0.0
+    np.testing.assert_array_equal(orbit.axis, [1.0, 0.0, 0.0])
+    with pytest.raises(ValueError):
+        orbit.axis[0] = 2.0
+    result = run_axial_gauge_query_study()
+    assert result == run_axial_gauge_query_study()
+    assert result["stationary_derivative_counterexample"]["finite_orbit_interval_m"] == [-0.1, 0.1]
+    assert result["action_control"]["within_budget"] is False
+    assert result["shared_action_nuisance_control"]["action_zero_within_zero_budget"] is True
+    assert all(value == 0 for value in result["information_boundary"].values())
+
+
+@pytest.mark.parametrize("axis", [[0.0, 0.0, 0.0], [np.nan, 0.0, 0.0], [1.0, 2.0]])
+def test_invalid_axis_rejected(axis: list[float]) -> None:
+    with pytest.raises(ValueError):
+        AxialGaugeOrbit(np.zeros(3), np.array(axis))
+
+
+@pytest.mark.parametrize("value", [-1.0, np.nan, np.inf])
+def test_invalid_regret_budget_rejected(value: float) -> None:
+    family = AxialQueryFamily((), np.ones(1), np.empty((1, 0)), np.empty((1, 0)))
+    with pytest.raises(ValueError):
+        family.within_regret_budget(0, maximum_regret=value)
+    with pytest.raises(ValueError):
+        family.within_regret_budget(0, maximum_regret=0.0, numerical_margin=value)
+
+
+@pytest.mark.parametrize("action", [True, 0.5, -1, 2])
+def test_invalid_action_rejected(action: object) -> None:
+    family = AxialQueryFamily((), np.ones(1), np.empty((1, 0)), np.empty((1, 0)))
+    with pytest.raises((ValueError, TypeError)):
+        family.regret_witness(action)
+
+
+@pytest.mark.parametrize("labels", [("missing",), ("g", "g"), ("",), (" g",), "g"])
+def test_invalid_group_lineage_rejected(labels: object) -> None:
+    with pytest.raises(ValueError):
+        affine_axial_queries(
+            np.zeros((1, 3)), np.ones((1, 1, 3)), point_group_ids=labels, orbits={"g": _orbit()}
+        )
+
+
+def test_bad_shapes_and_nonfinite_values_rejected() -> None:
+    with pytest.raises(ValueError):
+        _family(np.zeros((0, 3)), np.ones((1, 0, 3)))
+    with pytest.raises(ValueError):
+        _family(np.zeros((1, 3)), np.ones((1, 3)))
+    with pytest.raises(ValueError):
+        _family(np.full((1, 3), np.inf), np.ones((1, 1, 3)))
+    with pytest.raises(ValueError):
+        AxialQueryFamily("g", np.ones(1), np.ones((1, 1)), np.ones((1, 1)))
+    with pytest.raises(ValueError):
+        AxialQueryFamily(("g", "g"), np.ones(1), np.ones((1, 2)), np.ones((1, 2)))
+    with pytest.raises(ValueError):
+        AxialQueryFamily(("g",), np.ones(1), np.ones((2, 1)), np.ones((1, 1)))
+    with pytest.raises(ValueError):
+        affine_axial_queries(
+            np.zeros((1, 3)),
+            np.ones((1, 1, 3)),
+            offsets=np.ones(2),
+            point_group_ids=("g",),
+            orbits={"g": _orbit()},
+        )
+    with pytest.raises(ValueError):
+        _orbit().transform(np.zeros((1, 3)), np.nan)
+
+
+def test_existing_local_gate_parity() -> None:
+    # Imports must succeed in the complete repository: no silent integration skip.
+    from prob4d.query_observability import (
+        QueryObservabilityGate,
+        evaluate_query_observability,
+        point_position_query_jacobian,
+    )
+    from prob4d.query_observability_study import _controlled_factor
+
+    factor = _controlled_factor(complete_nullspace=False)
+    jacobian = point_position_query_jacobian(factor, np.array([0.0, 0.1, 0.0]))[1:2]
+    report = evaluate_query_observability(
+        factor, prior_covariance_local=np.eye(7), query_jacobian_local=jacobian
+    )
+    gate = QueryObservabilityGate(0.8, 0.8, 0.5)
+    reference = run_axial_gauge_query_study()["stationary_derivative_counterexample"]
+    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(10.0 / 11.0)
+    assert report.worst_supported_variance_ratio == pytest.approx(1.0 / 11.0)
+    expected = reference["local_linear_reference"]["local_gate_admits"]
+    assert gate.evaluate(report).admitted is expected