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.
This commit is contained in:
parent
d376339049
commit
26270ed846
3 changed files with 546 additions and 1 deletions
|
|
@ -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
|
||||
|
|
|
|||
275
algebra/null_point.py
Normal file
275
algebra/null_point.py
Normal file
|
|
@ -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)
|
||||
252
tests/test_null_point_primitives.py
Normal file
252
tests/test_null_point_primitives.py
Normal file
|
|
@ -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
|
||||
Loading…
Reference in a new issue