From 26270ed846efa5266cbec0963c71d73c28718acf Mon Sep 17 00:00:00 2001 From: Shay Date: Sun, 12 Jul 2026 16:28:19 -0700 Subject: [PATCH] feat(algebra): null-point recovery primitives + frozen CGA null constants Shared CGA substrate for the #17 conformal-Procrustes/Kabsch and #16 Cartan-Iwasawa decompositions. Adds algebra/null_point.py and hoists the two conformal null directions to frozen module constants. algebra/cga.py - Add frozen read-only f64 N_O / N_INF constants: the same vectors embed_point builds inline (origin embeds to N_O; N_INF is fixed by every Euclidean isometry), so the null-point primitives share one exact sign definition instead of re-deriving it per call site. - Fix header-docstring sign typo: n_o = 0.5*(e5 - e4), not 0.5*(e4 - e5). embed_point was already correct; only the module header disagreed. algebra/null_point.py (new) - dilator(scale), translator(a): CGA similarity constructors; both round-trip through the recoverers. - recover_dilation(V) -> (scale, D): reads V n_inf rev(V), weight-normalised so recovery is invariant to a non-unit versor weight (verified vs V -> kV). - recover_translation(V) -> (a, T): reads V n_o rev(V), projective dehomogenisation. - NullPointRecoveryError carries machine-readable reason codes. - Fail-closed symmetric similarity gate (_require_similarity): BOTH recoverers now reject non-versors (not_a_versor) and non-similarities (not_similarity, e.g. transversions). Closes an asymmetry where recover_translation silently accepted a transversion / non-versor and returned a plausible translation, violating the module's own wrong=0 contract. - Orientation-reversing (reflection / det=-1) versors are refused by recover_dilation with a distinct reason improper_versor, kept separate from degenerate_scale; recover_translation still accepts them (the origin image is well defined). conformal_procrustes strips reflections upstream, so this is a documented boundary, not a silent one. - Default tol=1e-9 documented: matches f64-exact recovery of a cleanly assembled versor (~1e-14 round-trip); noisy/SVD callers must pass a wider tol. tests/test_null_point_primitives.py (new): 33 tests - null-cone/pairing invariants, constant immutability, constructor round-trips, composed T.D.R peel, versor-weight invariance, and the full fail-closed matrix (transversion, non-versor, inversion, reflection asymmetry, non-positive scale, bad vector). Invariant protected: wrong=0 - no recovery returns a silently wrong value on a degenerate / non-versor / non-similarity input. Validation: 33/33 new pass; 88 passed / 1 xfailed across the CGA substrate + physics Procrustes consumers (dynamic_manifold, surprise, versor closure, rotor, holonomy). Hardened via a 3-lens adversarial verification (soundness / sign-convention / consumer-contract, each executing counterexample versors, every finding skeptic-verified): 2 CONFIRMED findings fixed (asymmetric validation gap; reflection reason conflation); tol-tightness resolved by documentation rather than a guard-weakening default change. --- algebra/cga.py | 20 +- algebra/null_point.py | 275 ++++++++++++++++++++++++++++ tests/test_null_point_primitives.py | 252 +++++++++++++++++++++++++ 3 files changed, 546 insertions(+), 1 deletion(-) create mode 100644 algebra/null_point.py create mode 100644 tests/test_null_point_primitives.py diff --git a/algebra/cga.py b/algebra/cga.py index a4cdf9a8..4c1e5f78 100644 --- a/algebra/cga.py +++ b/algebra/cga.py @@ -4,7 +4,7 @@ Conformal Geometric Algebra geometry on Cl(4,1). Signature: (+,+,+,+,-), with Euclidean coordinates on e1,e2,e3. The two conformal null directions are built from e4 and e5: - n_o = 0.5 * (e4 - e5) # origin, n_o^2 = 0 + n_o = 0.5 * (e5 - e4) # origin, n_o^2 = 0 n_inf = e4 + e5 # infinity, n_inf^2 = 0 n_o · n_inf = -1 @@ -39,6 +39,24 @@ _I5[_PSEUDOSCALAR_INDEX] = 1.0 _E4_IDX = 4 _E5_IDX = 5 +# The two conformal null directions, frozen as f64 32-vectors — the canonical +# origin/infinity of the CGA point map. These are the SAME vectors ``embed_point`` +# builds inline (origin embeds to N_O; N_INF is fixed by every Euclidean isometry), +# hoisted to module constants so the null-point recovery primitives (dilation / +# translation peel) and any incidence code share one exact definition instead of +# re-deriving the signs. Invariants (pinned in tests/test_null_point_primitives.py): +# N_O · N_O = 0, N_INF · N_INF = 0, N_O · N_INF = -1. +# Never mutated; callers that need a scratch copy must ``.copy()``. +N_O = np.zeros(N_COMPONENTS, dtype=np.float64) +N_O[_E4_IDX] = -0.5 # n_o = 0.5 * (e5 - e4) +N_O[_E5_IDX] = 0.5 +N_O.setflags(write=False) + +N_INF = np.zeros(N_COMPONENTS, dtype=np.float64) +N_INF[_E4_IDX] = 1.0 # n_inf = e4 + e5 +N_INF[_E5_IDX] = 1.0 +N_INF.setflags(write=False) + # Pinned magnitude ceiling for f64-exact embedding + read-back (Phase 0A). # Below this bound, ``embed_point(..., dtype=np.float64)`` round-trips integer # coordinates exactly through ``read_scalar_e1`` and the conformal distance metric diff --git a/algebra/null_point.py b/algebra/null_point.py new file mode 100644 index 00000000..c10006a8 --- /dev/null +++ b/algebra/null_point.py @@ -0,0 +1,275 @@ +"""Null-point recovery primitives for CGA conformal versors. + +Shared substrate for the conformal-Procrustes (#17) and Cartan–Iwasawa (#16) +decompositions. Given a *similarity* versor V (rotation · dilation · translation, +in any order), these peel off the translation it applies to the origin and the +uniform dilation it applies to lengths, using only the exact CGA sandwich +``V·X·rev(V)`` on the two null directions ``N_O`` / ``N_INF`` (see algebra/cga.py: +``n_o = 0.5(e5 - e4)``, ``n_inf = e4 + e5``). + +Empirically pinned (f64-exact; probes reproduced in the test module): + + * ``V n_inf rev(V) = scale · n_inf`` — a similarity FIXES the point at + infinity, so its n_inf image is a *pure* positive multiple of n_inf whose + coefficient is the dilation factor. Anything else — a transversion / special + conformal versor — leaves an off-n_inf residual and is REFUSED. + * ``V n_o rev(V) = w_o·n_o + scale^-1·a + …`` — the origin's image is a + conformal point; ``a = euclidean_part / w_o`` recovers the translation by + projective dehomogenization (the weight divides out the dilation, and + rotation fixes the origin, so ``a`` is exact regardless of V's rotation or + scale content — the same trick as :func:`algebra.cga.read_scalar_e1`). + +Conventions — both constructors round-trip through the recoverers: + ``dilator(scale)`` scales Euclidean lengths by ``scale`` (> 0); + ``recover_dilation(dilator(s)) == s``. + ``translator(a)`` maps the origin to Euclidean point ``a`` (3-vector); + ``recover_translation(translator(a)) == a``. + +Fail-closed discipline (the wrong=0 rule): every recovery raises +:class:`NullPointRecoveryError` on a degenerate, non-versor, or non-similarity +input rather than returning a silently wrong value — ``recover_dilation`` and +``recover_translation`` share one versor+similarity gate +(:func:`_require_similarity`), so neither accepts what the other refuses. Guards +are scale-relative so a versor with non-unit weight (e.g. one assembled from a +Kabsch/SVD point cloud) is judged by its *shape*, not its magnitude. + +Tolerance: the default ``tol=1e-9`` matches the f64-exact recovery of a cleanly +assembled versor (an SVD-orthogonal rotation composed with an exact +dilator/translator round-trips to ~1e-14). A caller whose versor carries larger +numerical noise — e.g. an iteratively refined Procrustes fit — must pass a ``tol`` +at least as large as that residual, or a valid similarity may be refused as +``not_a_versor`` / ``not_similarity`` (fail-closed: it is never *accepted* with a +wrong value). ``core.physics.conformal_procrustes`` uses ``tol=1e-8`` by convention. +""" + +from __future__ import annotations + +import numpy as np + +from .cga import N_INF, N_O, cga_inner, graded_wedge +from .cl41 import N_COMPONENTS, geometric_product, reverse + +# e4 / e5 component indices inside the grade-1 block (mirror of algebra.cga; kept +# local to avoid importing a private name across modules). +_E4_IDX = 4 +_E5_IDX = 5 + +# The dilation bivector E = n_o ^ n_inf. E^2 = +1 (boost-like), so the dilator is +# a hyperbolic exponential cosh + sinh·E. Frozen f64; never mutated. +_E_DILATION = graded_wedge(N_O, N_INF).astype(np.float64) +_E_DILATION.setflags(write=False) + + +class NullPointRecoveryError(ValueError): + """A versor is degenerate or not a similarity transform. + + Carries a machine-readable ``reason`` for callers that route on the failure + mode (e.g. #17 margin reporting) rather than only surfacing the message. + """ + + def __init__(self, message: str, *, reason: str) -> None: + super().__init__(message) + self.reason = reason + + +def _sandwich(V: np.ndarray, X: np.ndarray) -> np.ndarray: + """The raw f64 sandwich ``V X rev(V)`` — no closure, no unitisation. + + (Deliberately not :func:`algebra.versor.versor_apply`: that path unitises + non-null inputs and coerces to the runtime field dtype. Null-point recovery + needs the exact algebraic image in f64.) + """ + V = np.asarray(V, dtype=np.float64) + X = np.asarray(X, dtype=np.float64) + return geometric_product(geometric_product(V, X), reverse(V)) + + +def dilator(scale: float) -> np.ndarray: + """Uniform-scale versor that scales Euclidean lengths by ``scale`` (> 0). + + ``D = exp(0.5·ln(scale)·E) = cosh(h) + sinh(h)·E`` with ``h = 0.5·ln(scale)`` + and ``E = n_o ^ n_inf`` (``E^2 = +1``). Acts as + ``D n_inf rev(D) = scale·n_inf`` and ``D n_o rev(D) = scale^-1·n_o``. + """ + scale = float(scale) + if not np.isfinite(scale) or scale <= 0.0: + raise NullPointRecoveryError( + f"dilator scale must be finite and positive, got {scale}", + reason="nonpositive_scale", + ) + half = 0.5 * np.log(scale) + D = np.zeros(N_COMPONENTS, dtype=np.float64) + D[0] = np.cosh(half) + D = D + np.sinh(half) * _E_DILATION + return D + + +def translator(a: np.ndarray) -> np.ndarray: + """Translator versor that maps the origin to Euclidean point ``a`` (3-vector). + + ``T = 1 - 0.5·a·n_inf`` (a embedded on e1..e3). ``T n_o rev(T)`` equals the + conformal embedding of ``a`` (== :func:`algebra.cga.embed_point`). + """ + a = np.asarray(a, dtype=np.float64) + if a.shape != (3,) or not np.all(np.isfinite(a)): + raise NullPointRecoveryError( + f"translator expects a finite 3-vector, got shape {a.shape}", + reason="bad_translation_vector", + ) + a_mv = np.zeros(N_COMPONENTS, dtype=np.float64) + a_mv[1:4] = a + T = np.zeros(N_COMPONENTS, dtype=np.float64) + T[0] = 1.0 + T = T - 0.5 * geometric_product(a_mv, N_INF) + return T + + +def _versor_scalar_weight(V: np.ndarray, tol: float) -> float: + """Return ``scalar_part(V·rev(V))`` after checking ``V`` is a versor. + + A versor satisfies ``V·rev(V) = scalar``; a non-versor multivector leaves an + off-scalar residual. Raises :class:`NullPointRecoveryError` (``not_a_versor`` + / ``degenerate_weight``) otherwise. The weight is what makes + :func:`recover_dilation` weight-invariant — the raw ``n_inf`` coefficient + scales with this weight, so the true dilation is the coefficient divided by it. + """ + V = np.asarray(V, dtype=np.float64) + vv = geometric_product(V, reverse(V)) + w = float(vv[0]) + off_scalar = float(np.linalg.norm(vv[1:])) + ref = max(1.0, abs(w)) + if off_scalar > tol * ref: + raise NullPointRecoveryError( + f"V·rev(V) is not scalar (off-scalar residual {off_scalar / ref:.3e}); " + "not a versor", + reason="not_a_versor", + ) + if abs(w) <= tol: + raise NullPointRecoveryError( + f"degenerate versor weight {w:.3e}", reason="degenerate_weight", + ) + return w + + +def _require_similarity(V: np.ndarray, tol: float) -> tuple[float, float]: + """Gate ``V`` as a similarity versor; return ``(weight, signed_scale)``. + + A similarity (rotation · dilation · translation, in any order) is the only + class both recoverers accept: it is a versor (``V·rev(V)`` scalar) *and* it + fixes infinity (``V n_inf rev(V)`` is a pure multiple of ``n_inf``). The + returned ``signed_scale = c_inf / weight`` is positive for a proper similarity + and negative for an orientation-reversing (improper / reflection) one; sign + and degeneracy classification is left to the caller, so + :func:`recover_translation` can accept a reflection — whose origin image is + still well defined — while :func:`recover_dilation` refuses it. + + Raises :class:`NullPointRecoveryError` with ``not_a_versor`` / + ``degenerate_weight`` (from :func:`_versor_scalar_weight`) or ``not_similarity``. + """ + weight = _versor_scalar_weight(V, tol) + W = _sandwich(V, N_INF) + c_inf = 0.5 * (float(W[_E4_IDX]) + float(W[_E5_IDX])) + + resid = W.copy() + resid[_E4_IDX] -= c_inf + resid[_E5_IDX] -= c_inf + resid_norm = float(np.linalg.norm(resid)) + ref = max(1.0, float(np.linalg.norm(W))) + if resid_norm > tol * ref: + raise NullPointRecoveryError( + f"versor does not fix infinity (off-n_inf residual " + f"{resid_norm / ref:.3e} > {tol:.1e}); not a similarity transform", + reason="not_similarity", + ) + return weight, c_inf / weight + + +def recover_dilation(V: np.ndarray, *, tol: float = 1e-9) -> tuple[float, np.ndarray]: + """Recover the uniform scale a similarity versor ``V`` applies to lengths. + + Returns ``(scale, D)`` with ``D == dilator(scale)`` and ``scale > 0``. Reads + the image of the point at infinity ``W = V n_inf rev(V)`` (for a similarity a + pure multiple of ``n_inf``) and normalises its coefficient by the versor weight + ``V·rev(V)`` — the sandwich scales with that weight, so a non-unit versor still + yields the true scale (verified against ``V -> kV``). + + Raises :class:`NullPointRecoveryError` when + * ``V`` is not a versor (``not_a_versor`` / ``degenerate_weight``) or does + not fix infinity, i.e. is not a similarity (``not_similarity`` — e.g. a + transversion); + * ``V`` is orientation-reversing — a reflection / improper rotation, the + ``det = -1`` case a raw Kabsch/SVD fit produces before it strips the + reflection (``core.physics.conformal_procrustes`` does strip it). Its + signed scale is a clean negative, refused as ``improper_versor``, kept + distinct from true degeneracy so a caller can tell "flip a singular + vector" from "numerically broken". :func:`recover_translation` still + accepts such a versor — only the *dilation* is ill-defined for an improper + map here; or + * the recovered scale is non-finite or collapses to zero (``degenerate_scale``). + """ + _, scale = _require_similarity(V, tol) + # Preserve the original accept-set exactly (finite *positive* scale, any + # magnitude); split the negative case out to a distinct, honest reason. + if not np.isfinite(scale): + raise NullPointRecoveryError( + f"degenerate dilation coefficient {scale}", + reason="degenerate_scale", + ) + if scale < 0.0: + raise NullPointRecoveryError( + f"orientation-reversing versor (signed scale {scale:.6g}); an improper " + "similarity has no positive dilation — strip the reflection first", + reason="improper_versor", + ) + if scale == 0.0: + raise NullPointRecoveryError( + "degenerate dilation coefficient 0.0 (versor collapses n_inf)", + reason="degenerate_scale", + ) + return scale, dilator(scale) + + +def recover_translation(V: np.ndarray, *, tol: float = 1e-9) -> tuple[np.ndarray, np.ndarray]: + """Recover the translation a similarity versor ``V`` applies to the origin. + + Returns ``(a, T)`` with ``a`` the Euclidean image of the origin (3-vector) + and ``T == translator(a)``. Reads ``W = V n_o rev(V)`` and dehomogenizes + projectively: ``a = W[e1:e3+1] / w_o`` where ``w_o = W[e5] - W[e4]``. The + weight divides out any dilation, and rotation — proper *or* a reflection — + fixes the origin, so ``a`` is exact regardless of ``V``'s rotation/scale + content. An improper (reflection) similarity is therefore accepted here even + though :func:`recover_dilation` refuses it: the origin image is well defined, + only the positive dilation is not. + + Gates ``V`` as a similarity versor first (the same :func:`_require_similarity` + gate as :func:`recover_dilation`), so a non-versor or a non-similarity — e.g. a + transversion, which fixes the origin and would otherwise return a plausible + ``a`` silently — fails closed rather than returning a wrong value. + + Raises :class:`NullPointRecoveryError` when + * ``V`` is not a versor (``not_a_versor`` / ``degenerate_weight``) or does + not fix infinity (``not_similarity``); + * the origin maps to infinity (``origin_at_infinity`` — ``|w_o|`` at/below + ``tol``; guards the projective division, subsumed by the similarity gate + for genuine inversions); or + * the origin image leaves the null cone (``non_null_image`` — scale-relative + defect > ``tol``), so ``W`` is not a conformal point. + """ + _require_similarity(V, tol) + W = _sandwich(V, N_O) + w_o = float(W[_E5_IDX] - W[_E4_IDX]) + if abs(w_o) <= tol: + raise NullPointRecoveryError( + f"origin maps to infinity (n_o weight {w_o:.3e}); no finite translation", + reason="origin_at_infinity", + ) + null_defect = abs(cga_inner(W, W)) + ref = max(1.0, float(np.dot(W, W))) + if null_defect > tol * ref: + raise NullPointRecoveryError( + f"origin image leaves the null cone (defect {null_defect / ref:.3e}); " + "not a conformal point", + reason="non_null_image", + ) + a = np.asarray(W[1:4], dtype=np.float64) / w_o + return a, translator(a) diff --git a/tests/test_null_point_primitives.py b/tests/test_null_point_primitives.py new file mode 100644 index 00000000..0f130cba --- /dev/null +++ b/tests/test_null_point_primitives.py @@ -0,0 +1,252 @@ +"""Pin tests for the conformal null-point primitives. + +These lock the CGA null geometry that the shared #17 (Kabsch / conformal +Procrustes) and #16 (Cartan–Iwasawa) recovery helpers stand on: + + * the frozen ``N_O`` / ``N_INF`` constants agree exactly with the vectors + ``embed_point`` builds inline, and + * the sign convention is ``n_o = 0.5 * (e5 - e4)`` (NOT ``e4 - e5`` — the old + module-header docstring had that backwards). + +The inner-product identities are exact (0.5 / 1.0 are representable), so the +tolerances are tight. +""" + +import numpy as np +import pytest + +from algebra.cga import N_INF, N_O, cga_inner, embed_point, read_scalar_e1 +from algebra.cl41 import basis_vector, geometric_product, reverse, scalar_part +from algebra.null_point import ( + NullPointRecoveryError, + dilator, + recover_dilation, + recover_translation, + translator, + _E_DILATION, +) +from algebra.rotor import make_rotor_from_angle + + +def _sandwich(V, X): + V = np.asarray(V, dtype=np.float64) + X = np.asarray(X, dtype=np.float64) + return geometric_product(geometric_product(V, X), reverse(V)) + + +def test_n_o_sign_convention_matches_basis_vectors(): + """n_o = 0.5 * (e5 - e4); basis_vector is 0-indexed so e4=bv(3), e5=bv(4).""" + n_o = 0.5 * (basis_vector(4) - basis_vector(3)) + assert np.allclose(N_O, n_o), "frozen N_O disagrees with 0.5*(e5 - e4)" + + +def test_n_inf_matches_basis_vectors(): + n_inf = basis_vector(3) + basis_vector(4) # e4 + e5 + assert np.allclose(N_INF, n_inf), "frozen N_INF disagrees with e4 + e5" + + +def test_null_cone_invariants(): + """N_O and N_INF both lie on the null cone: X . X = 0.""" + assert abs(cga_inner(N_O, N_O)) < 1e-12 + assert abs(cga_inner(N_INF, N_INF)) < 1e-12 + + +def test_no_ninf_pairing_is_minus_one(): + """N_O . N_INF = -1 exactly. Parenthesis is INSIDE abs: abs(x + 1), not abs(x) + 1.""" + assert abs(cga_inner(N_O, N_INF) + 1.0) < 1e-12 + + +def test_embed_origin_is_n_o(): + """The Euclidean origin embeds to n_o: e4 coeff = -0.5, e5 coeff = +0.5.""" + x0 = embed_point(np.zeros(3), dtype=np.float64) + assert np.allclose(x0[4:6], [-0.5, 0.5]) + # ...and the whole embedding equals N_O (origin has zero Euclidean part). + assert np.allclose(x0, N_O) + + +def test_constants_are_read_only(): + """The module constants must not be mutable in place.""" + for const in (N_O, N_INF): + assert const.flags.writeable is False + + +# --------------------------------------------------------------------------- +# Recovery primitives: constructors, round-trips, composed peel, fail-closed. +# --------------------------------------------------------------------------- + + +def test_E_dilation_squares_to_one(): + """The dilation bivector E = n_o ^ n_inf lives at index 15 and E^2 = +1.""" + assert _E_DILATION[15] == -1.0 + assert np.count_nonzero(_E_DILATION) == 1 + e_sq = geometric_product(_E_DILATION, _E_DILATION) + assert abs(scalar_part(e_sq) - 1.0) < 1e-12 + assert np.linalg.norm(e_sq[1:]) < 1e-12 # pure scalar + + +def test_dilator_scales_euclidean_lengths(): + """dilator(s) scales the Euclidean coordinate of a point by s.""" + X = embed_point(np.array([3.0, 0.0, 0.0]), dtype=np.float64) + for s in (2.0, 0.5, 4.0): + Y = _sandwich(dilator(s), X) + assert abs(read_scalar_e1(Y) - s * 3.0) < 1e-9 + + +def test_translator_maps_origin_to_point(): + """translator(a) carries the origin exactly to embed_point(a).""" + for a in ([1.0, 0.0, 0.0], [2.0, -1.0, 0.5]): + a = np.array(a) + image = _sandwich(translator(a), N_O) + assert np.allclose(image, embed_point(a, dtype=np.float64), atol=1e-9) + + +@pytest.mark.parametrize("scale", [2.5, 0.4, 1.0, 7.0, 0.125]) +def test_recover_dilation_round_trip(scale): + rec_scale, D = recover_dilation(dilator(scale)) + assert abs(rec_scale - scale) < 1e-9 + assert np.allclose(D, dilator(scale), atol=1e-12) + + +@pytest.mark.parametrize("a", [[1.5, -0.5, 2.0], [-3.0, 1.0, 0.0], [0.0, 0.0, 0.0]]) +def test_recover_translation_round_trip(a): + a = np.array(a) + rec_a, T = recover_translation(translator(a)) + assert np.allclose(rec_a, a, atol=1e-9) + assert np.allclose(T, translator(a), atol=1e-12) + + +@pytest.mark.parametrize( + "scale,a,angle", + [(2.5, [1.5, -0.5, 2.0], 0.7), (0.4, [-3.0, 1.0, 0.0], 1.9), (3.0, [0.2, 0.2, 0.2], -1.1)], +) +def test_recover_from_composed_similarity(scale, a, angle): + """V = T . D . R : dilation and translation peel out exactly, rotation and + each other's presence notwithstanding.""" + a = np.array(a) + R = make_rotor_from_angle(angle, bivector_idx=6).astype(np.float64) + V = geometric_product(geometric_product(translator(a), dilator(scale)), R) + rec_scale, _ = recover_dilation(V) + rec_a, _ = recover_translation(V) + assert abs(rec_scale - scale) < 1e-8 + assert np.allclose(rec_a, a, atol=1e-8) + + +@pytest.mark.parametrize("k", [3.0, 0.5, 10.0]) +def test_recover_dilation_is_versor_weight_invariant(k): + """Regression: the raw n_inf coefficient scales with the versor weight; the + recovered scale must NOT. recover_dilation(k*V) == recover_dilation(V).""" + V = geometric_product(translator(np.array([1.0, 2.0, -1.0])), dilator(2.5)) + base, _ = recover_dilation(V) + scaled, _ = recover_dilation(k * V) + assert abs(base - 2.5) < 1e-9 + assert abs(scaled - base) < 1e-9 + + +def test_recover_translation_is_weight_invariant(): + V = geometric_product(translator(np.array([1.0, 2.0, -1.0])), dilator(2.5)) + a0, _ = recover_translation(V) + a1, _ = recover_translation(3.0 * V) + assert np.allclose(a0, [1.0, 2.0, -1.0], atol=1e-9) + assert np.allclose(a1, a0, atol=1e-9) + + +def test_recover_dilation_refuses_transversion(): + """A transversion (special conformal) does NOT fix infinity -> not_similarity.""" + b = np.zeros(32) + b[1] = 0.3 + K = np.zeros(32) + K[0] = 1.0 + K = K - 0.5 * geometric_product(b, N_O) # transversion = 1 - 0.5 b n_o + with pytest.raises(NullPointRecoveryError) as exc: + recover_dilation(K) + assert exc.value.reason == "not_similarity" + + +def test_recover_dilation_refuses_non_versor(): + """A mixed-grade multivector is not a versor -> not_a_versor.""" + bad = np.zeros(32) + bad[0] = 1.0 + bad[1] = 1.0 + bad[6] = 1.0 # scalar + vector + bivector: V rev(V) not scalar + with pytest.raises(NullPointRecoveryError) as exc: + recover_dilation(bad) + assert exc.value.reason == "not_a_versor" + + +def test_recover_translation_refuses_inversion_as_not_similarity(): + """Unit-sphere inversion sigma = n_o - 0.5 n_inf swaps the null directions, so + it is not a similarity (does not fix infinity). The shared similarity gate + refuses it as not_similarity — the fundamental cause. (It also sends the origin + to infinity; origin_at_infinity remains a defensive division guard, subsumed + here for genuine inversions.)""" + sigma = N_O - 0.5 * N_INF # sigma^2 = 1, an honest inversion reflector + with pytest.raises(NullPointRecoveryError) as exc: + recover_translation(sigma) + assert exc.value.reason == "not_similarity" + + +def test_recover_translation_refuses_transversion(): + """Symmetric with recover_dilation: a transversion IS a versor and fixes the + origin, so without the similarity gate it silently returned a plausible + a=[0,0,0]. The gate must refuse it (it does not fix infinity) -> not_similarity.""" + b = np.zeros(32) + b[1] = 0.3 + K = np.zeros(32) + K[0] = 1.0 + K = K - 0.5 * geometric_product(b, N_O) # transversion = 1 - 0.5 b n_o + with pytest.raises(NullPointRecoveryError) as exc: + recover_translation(K) + assert exc.value.reason == "not_similarity" + + +def test_recover_translation_refuses_non_versor(): + """A mixed-grade multivector is not a versor -> not_a_versor (symmetric with + recover_dilation; previously recover_translation accepted it and returned a + silent value).""" + bad = np.zeros(32) + bad[0] = 1.0 + bad[1] = 1.0 + bad[6] = 1.0 # scalar + vector + bivector: V rev(V) not scalar + with pytest.raises(NullPointRecoveryError) as exc: + recover_translation(bad) + assert exc.value.reason == "not_a_versor" + + +def _reflection_similarity(a, scale): + """T . D . (e1-reflection): an orientation-reversing (det=-1) similarity, what a + raw Kabsch/SVD fit yields before it strips the reflection.""" + return geometric_product( + geometric_product(translator(np.array(a)), dilator(scale)), basis_vector(0) + ) + + +def test_recover_dilation_refuses_reflection_as_improper(): + """A reflection (improper rotation, det=-1) is refused as improper_versor, NOT + degenerate_scale: its scale magnitude is a clean, well-conditioned number, and + the distinct reason lets a consumer route 'strip the reflection' vs 'broken'.""" + V = _reflection_similarity([1.0, 0.5, -0.3], 2.0) + with pytest.raises(NullPointRecoveryError) as exc: + recover_dilation(V) + assert exc.value.reason == "improper_versor" + + +def test_recover_translation_accepts_reflection(): + """Asymmetry by design: the origin image is well defined under a reflection, so + recover_translation SUCCEEDS on the very versor recover_dilation refuses.""" + V = _reflection_similarity([1.0, 0.5, -0.3], 2.0) + rec_a, _ = recover_translation(V) + assert np.allclose(rec_a, [1.0, 0.5, -0.3], atol=1e-9) + + +def test_dilator_rejects_nonpositive_scale(): + for bad in (0.0, -1.0, float("inf"), float("nan")): + with pytest.raises(NullPointRecoveryError) as exc: + dilator(bad) + assert exc.value.reason == "nonpositive_scale" + + +def test_translator_rejects_bad_vector(): + with pytest.raises(NullPointRecoveryError): + translator(np.array([1.0, 2.0])) # wrong shape + with pytest.raises(NullPointRecoveryError): + translator(np.array([1.0, np.nan, 0.0])) # non-finite