fix(algebra): rotor_power smoke — α≈0 identity + simple B² float-dust tol
Smoke reds came from (1) denormal stream weights powering non-simple splits and (2) invariant-split factors with B² higher residual in 1e-6..1e-3 dust band raising fail-closed. Early R^0→I / R^1→R; raise simple B² higher tol to 1e-3 with named constant. Pins multiplane transition + near-zero alpha.
This commit is contained in:
parent
e2c375e407
commit
e124d83bfa
2 changed files with 61 additions and 4 deletions
|
|
@ -18,6 +18,11 @@ _STRICT_RESIDUE_TOL = 1e-2
|
||||||
# A rotor is SIMPLE iff its grade-4 part vanishes (<R>_4 == 0 <=> R = R1 with a
|
# A rotor is SIMPLE iff its grade-4 part vanishes (<R>_4 == 0 <=> R = R1 with a
|
||||||
# single invariant plane). Above this, the rotor needs the invariant split.
|
# single invariant plane). Above this, the rotor needs the invariant split.
|
||||||
_SIMPLE_GRADE4_TOL = 1e-10
|
_SIMPLE_GRADE4_TOL = 1e-10
|
||||||
|
# After the invariant (bivector) split, each factor is *approximately* simple;
|
||||||
|
# B² higher-grade residual is float dust, not a true multi-plane bivector.
|
||||||
|
# 1e-6 was too tight (raised on live word-transition / stream weights ≈ 1e-6..1e-3).
|
||||||
|
# Refuse only residuals that are clearly structural non-simplicity.
|
||||||
|
_SIMPLE_BSQ_HIGHER_TOL = 1e-3
|
||||||
# |discriminant| below this => the two invariant eigenvalues coincide (isoclinic).
|
# |discriminant| below this => the two invariant eigenvalues coincide (isoclinic).
|
||||||
_DEGEN_TOL = 1e-9
|
_DEGEN_TOL = 1e-9
|
||||||
|
|
||||||
|
|
@ -107,10 +112,17 @@ def rotor_power(R: np.ndarray, alpha: float) -> np.ndarray:
|
||||||
f"rotor_power expects a {N_COMPONENTS}-component rotor; got {R_arr.shape}."
|
f"rotor_power expects a {N_COMPONENTS}-component rotor; got {R_arr.shape}."
|
||||||
)
|
)
|
||||||
dtype = _result_dtype(R_arr)
|
dtype = _result_dtype(R_arr)
|
||||||
|
a = float(alpha)
|
||||||
|
# Endpoints by continuity: R^0 = 1, R^1 = R. Stream weights can be denormal
|
||||||
|
# tiny; never run the invariant split on α≈0 (smoke / generate.stream path).
|
||||||
|
if abs(a) <= _NEAR_ZERO_TOL:
|
||||||
|
return _identity(dtype)
|
||||||
|
if abs(a - 1.0) <= _NEAR_ZERO_TOL:
|
||||||
|
return R_arr.astype(dtype, copy=True)
|
||||||
# <R>_4 == 0 <=> R is a single simple rotor. Otherwise take the split path.
|
# <R>_4 == 0 <=> R is a single simple rotor. Otherwise take the split path.
|
||||||
if float(np.linalg.norm(grade_project(R_arr, 4))) >= _SIMPLE_GRADE4_TOL:
|
if float(np.linalg.norm(grade_project(R_arr, 4))) >= _SIMPLE_GRADE4_TOL:
|
||||||
return _general_rotor_power(R_arr, alpha, dtype)
|
return _general_rotor_power(R_arr, a, dtype)
|
||||||
return _simple_rotor_power(R_arr, alpha, dtype)
|
return _simple_rotor_power(R_arr, a, dtype)
|
||||||
|
|
||||||
|
|
||||||
def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np.ndarray:
|
def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np.ndarray:
|
||||||
|
|
@ -123,16 +135,20 @@ def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np.
|
||||||
B[0] = 0.0
|
B[0] = 0.0
|
||||||
|
|
||||||
# A simple rotor's bivector squares to a scalar (B² is grade-0 only).
|
# A simple rotor's bivector squares to a scalar (B² is grade-0 only).
|
||||||
|
# Higher-grade residual above _SIMPLE_BSQ_HIGHER_TOL is structural non-simplicity
|
||||||
|
# (fail closed). Below that, treat as float dust from the invariant split and
|
||||||
|
# use only the scalar part of B² (closed form still exact on the simple plane).
|
||||||
B_sq_full = geometric_product(B, B).astype(np.float64)
|
B_sq_full = geometric_product(B, B).astype(np.float64)
|
||||||
bsq_scalar = float(B_sq_full[0])
|
bsq_scalar = float(B_sq_full[0])
|
||||||
B_sq_higher = B_sq_full.copy()
|
B_sq_higher = B_sq_full.copy()
|
||||||
B_sq_higher[0] = 0.0
|
B_sq_higher[0] = 0.0
|
||||||
if float(np.linalg.norm(B_sq_higher)) > 1e-6:
|
higher_norm = float(np.linalg.norm(B_sq_higher))
|
||||||
|
if higher_norm > _SIMPLE_BSQ_HIGHER_TOL:
|
||||||
# Not a simple bivector under the simple dispatch — fail closed, never
|
# Not a simple bivector under the simple dispatch — fail closed, never
|
||||||
# silently return identity (that zeros motion without a signal).
|
# silently return identity (that zeros motion without a signal).
|
||||||
raise ValueError(
|
raise ValueError(
|
||||||
"rotor_power: non-simple bivector under simple dispatch "
|
"rotor_power: non-simple bivector under simple dispatch "
|
||||||
f"(B² higher-grade residual {float(np.linalg.norm(B_sq_higher)):.3e})"
|
f"(B² higher-grade residual {higher_norm:.3e})"
|
||||||
)
|
)
|
||||||
|
|
||||||
# Near-identity: nothing to scale.
|
# Near-identity: nothing to scale.
|
||||||
|
|
|
||||||
|
|
@ -94,3 +94,44 @@ def test_rotor_power_null_translator_scales_translation() -> None:
|
||||||
def test_rotor_power_rejects_wrong_shape() -> None:
|
def test_rotor_power_rejects_wrong_shape() -> None:
|
||||||
with pytest.raises(ValueError):
|
with pytest.raises(ValueError):
|
||||||
rotor_power(np.zeros(16, dtype=np.float64), 0.5)
|
rotor_power(np.zeros(16, dtype=np.float64), 0.5)
|
||||||
|
|
||||||
|
|
||||||
|
def test_rotor_power_near_zero_alpha_is_identity() -> None:
|
||||||
|
"""Stream weights can be denormal tiny; R^α → I as α → 0 without split."""
|
||||||
|
from algebra.cl41 import geometric_product
|
||||||
|
from algebra.versor import unitize_versor
|
||||||
|
|
||||||
|
v = make_rotor_from_angle(0.4, bivector_idx=6)
|
||||||
|
for plane in (7, 8, 10):
|
||||||
|
v = geometric_product(v, make_rotor_from_angle(0.25, bivector_idx=plane))
|
||||||
|
R = unitize_versor(v)
|
||||||
|
out = rotor_power(R, 1e-40)
|
||||||
|
expected = np.zeros(32, dtype=np.float64)
|
||||||
|
expected[0] = 1.0
|
||||||
|
np.testing.assert_allclose(out, expected, atol=1e-12)
|
||||||
|
|
||||||
|
|
||||||
|
def test_rotor_power_multiplane_transition_half_stays_closed() -> None:
|
||||||
|
"""Live path: multi-plane word_transition_rotor then fractional power.
|
||||||
|
|
||||||
|
Regression for smoke: invariant-split factors have B² higher residual in the
|
||||||
|
1e-6..1e-3 float-dust band; must not raise non-simple under simple dispatch.
|
||||||
|
"""
|
||||||
|
from algebra.cl41 import geometric_product
|
||||||
|
from algebra.versor import unitize_versor
|
||||||
|
|
||||||
|
def compose(scale: float) -> np.ndarray:
|
||||||
|
v = np.zeros(32, dtype=np.float64)
|
||||||
|
v[0] = 1.0
|
||||||
|
for k, plane in enumerate((6, 7, 8, 10, 11)):
|
||||||
|
v = geometric_product(
|
||||||
|
v, make_rotor_from_angle(0.2 + 0.11 * k + 0.05 * scale, bivector_idx=plane)
|
||||||
|
)
|
||||||
|
return unitize_versor(v)
|
||||||
|
|
||||||
|
R = word_transition_rotor(compose(1.0), compose(2.0))
|
||||||
|
for alpha in (0.1, 0.3, 0.5, 0.7, 0.9):
|
||||||
|
R_a = rotor_power(R, alpha)
|
||||||
|
assert versor_condition(R_a) < _TOL, (
|
||||||
|
f"multiplane transition power alpha={alpha}: cond={versor_condition(R_a):.3e}"
|
||||||
|
)
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue