Merge pull request 'fix(algebra): rotor_power smoke — α≈0 identity + simple B² float-dust tol' (#33) from feat/fix-rotor-power-smoke-simple-dispatch into main
Unblock smoke: early R^0/R^1 exits; simple B² higher-grade float-dust tol 1e-3 after invariant split.
This commit is contained in:
commit
9ab4cca821
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
|
||||
# single invariant plane). Above this, the rotor needs the invariant split.
|
||||
_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).
|
||||
_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}."
|
||||
)
|
||||
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.
|
||||
if float(np.linalg.norm(grade_project(R_arr, 4))) >= _SIMPLE_GRADE4_TOL:
|
||||
return _general_rotor_power(R_arr, alpha, dtype)
|
||||
return _simple_rotor_power(R_arr, alpha, dtype)
|
||||
return _general_rotor_power(R_arr, a, dtype)
|
||||
return _simple_rotor_power(R_arr, a, dtype)
|
||||
|
||||
|
||||
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
|
||||
|
||||
# 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)
|
||||
bsq_scalar = float(B_sq_full[0])
|
||||
B_sq_higher = B_sq_full.copy()
|
||||
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
|
||||
# silently return identity (that zeros motion without a signal).
|
||||
raise ValueError(
|
||||
"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.
|
||||
|
|
|
|||
|
|
@ -94,3 +94,44 @@ def test_rotor_power_null_translator_scales_translation() -> None:
|
|||
def test_rotor_power_rejects_wrong_shape() -> None:
|
||||
with pytest.raises(ValueError):
|
||||
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