From 2050b77ab255259ad9111e45e434b240312f781e Mon Sep 17 00:00:00 2001 From: Shay Date: Mon, 13 Jul 2026 17:07:42 -0700 Subject: [PATCH] =?UTF-8?q?feat(third-door):=20real=20Cartan=E2=80=93Iwasa?= =?UTF-8?q?wa=20null-point=20peel=20+=20full=20Kabsch-conformal=20Procrust?= =?UTF-8?q?es=20(#16=20#17)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Cartan: recover_dilation → peel D → recover_translation → peel T; Spin remainder for non-similarities; strict close (no seed-to-rotor); recon residual fallback. Flips fidelity xfail. - Procrustes: full 5-D Kabsch on null-point clouds; field conjugacy via raw sandwich + Spin GN; delete word_transition_rotor averaging path. Non-vacuous harness fixture. - rotor_power: null-bivector power (a+B)^α = a^α + α a^{α-1} B so translators no longer silently zero under dual-slerp. - Ledger scorecard: #2 and #3 → 🟢; #4 remains 🟡 (bootstrap deferred). 549 passed (fidelity + ADR-0239 + null_point + 0240 + rotor_power). --- algebra/rotor.py | 35 +- core/physics/dynamic_manifold.py | 681 +++++++++++++++--- .../research/third-door-blueprint-fidelity.md | 48 +- evals/analogical_transfer/harness.py | 94 ++- tests/test_adr_0239_dynamic_manifold.py | 142 +++- tests/test_rotor_power.py | 14 + tests/test_third_door_blueprint_fidelity.py | 202 +++++- 7 files changed, 999 insertions(+), 217 deletions(-) diff --git a/algebra/rotor.py b/algebra/rotor.py index 3b39075e..11af941d 100644 --- a/algebra/rotor.py +++ b/algebra/rotor.py @@ -128,8 +128,12 @@ def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np. B_sq_higher = B_sq_full.copy() B_sq_higher[0] = 0.0 if float(np.linalg.norm(B_sq_higher)) > 1e-6: - # Not a simple bivector — should not reach here via the public dispatch. - return _identity(dtype) + # 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})" + ) # Near-identity: nothing to scale. bivector_norm = float(np.linalg.norm(B)) @@ -144,19 +148,30 @@ def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np. new_a = float(np.cos(alpha * theta_half)) new_b_mag = float(np.sin(alpha * theta_half)) elif bsq_scalar > 0.0: - # Boost plane. + # Boost plane. Domain of atanh requires |b_mag/a| < 1 and a > 0. b_mag = float(np.sqrt(bsq_scalar)) - # atanh requires |b_mag/a| < 1; for closed rotors a² - B² = 1 means - # |b_mag| < |a|, so this is safe when a > 0. - if a == 0.0: - return _identity(dtype) + if a <= 0.0 or abs(b_mag / a) >= 1.0 - 1e-12: + raise ValueError( + f"rotor_power: boost plane outside unit-rotor domain " + f"(a={a:.6g}, |B|/a={abs(b_mag / a) if a != 0.0 else float('inf'):.6g})" + ) eta_half = float(np.arctanh(b_mag / a)) new_a = float(np.cosh(alpha * eta_half)) new_b_mag = float(np.sinh(alpha * eta_half)) else: - # B² = 0: null bivector. Cannot interpolate on the manifold; - # return identity to fail safely. - return _identity(dtype) + # B² = 0: null bivector (translator generators in CGA). Exact binomial: + # (a + B)^α = a^α + α a^{α-1} B (higher powers of B vanish). + # Unit translators have a = 1 ⇒ T^α = 1 + α B = translator(α·a_eucl). + # Historically this returned identity — a silent zeroing of the Cartan + # translation leg in dual_correction_slerp (fidelity #16 follow-up). + if abs(a) < _NEAR_ZERO_TOL: + return _identity(dtype) + result = np.zeros(N_COMPONENTS, dtype=np.float64) + result[0] = float(a) ** float(alpha) if a > 0.0 else float(np.sign(a) * (abs(a) ** float(alpha))) + # Prefer real power for a>0; for a<0 (rare for unit translators) use |a|^α · sgn. + scale_B = float(alpha) * (float(a) ** (float(alpha) - 1.0)) if a > 0.0 else float(alpha) * (abs(a) ** (float(alpha) - 1.0)) * float(np.sign(a)) + result = result + scale_B * B + return result.astype(dtype, copy=False) result = np.zeros(N_COMPONENTS, dtype=np.float64) result[0] = new_a diff --git a/core/physics/dynamic_manifold.py b/core/physics/dynamic_manifold.py index d5c23fbb..18f843e1 100644 --- a/core/physics/dynamic_manifold.py +++ b/core/physics/dynamic_manifold.py @@ -16,13 +16,26 @@ from typing import Optional, Sequence, Tuple, Union import numpy as np -from algebra.cl41 import N_COMPONENTS, SIGNATURE, geometric_product, grade_project, reverse +from algebra.cga import is_null +from algebra.cl41 import N_COMPONENTS, geometric_product, grade_project, reverse +from algebra.null_point import ( + NullPointRecoveryError, + dilator, + recover_dilation, + recover_translation, + translator, +) from algebra.rotor import rotor_power, word_transition_rotor -from algebra.versor import unitize_versor, versor_apply, versor_condition +from algebra.versor import versor_condition _CLOSURE_TOL = 1e-6 _NEAR_ZERO = 1e-12 _NULL_TOL = 1e-9 +_PROCRUSTES_WEIGHT_TOL = 1e-8 +_CONJUGACY_RES_TOL = 1e-5 +_CONJUGACY_MAX_STEPS = 48 +# Grade-2 blade indices (e1∧e2 … spanning the 10-plane bivector of Cl(4,1)). +_BIVECTOR_PLANES = tuple(range(6, 16)) # Cl(4,1) metric on Euclidean+conformal R^5 _ETA5 = np.diag([1.0, 1.0, 1.0, 1.0, -1.0]).astype(np.float64) @@ -77,6 +90,30 @@ def _identity32() -> np.ndarray: return out +def _strict_close_versor(V: np.ndarray, *, name: str) -> np.ndarray: + """Rescale a true versor to unit weight; never seed-fabricate a rotor. + + A versor satisfies ``V·rev(V) = scalar``. If the product is not scalar, or + the scalar is non-positive, raise ``ValueError`` (fail-closed). This is the + construction-boundary closer for Cartan–Iwasawa — distinct from + :func:`unitize_versor`, which may map dense seeds onto the manifold. + """ + arr = np.asarray(V, dtype=np.float64) + if arr.shape != (N_COMPONENTS,): + raise ValueError(f"{name}: expected shape ({N_COMPONENTS},), got {arr.shape}") + product = geometric_product(arr, reverse(arr)).astype(np.float64) + scalar_sq = float(product[0]) + residue = product.copy() + residue[0] = 0.0 + residue_norm = float(np.linalg.norm(residue)) + if residue_norm >= 1e-2 or scalar_sq <= 0.0: + raise ValueError( + f"{name}: input not a versor " + f"(residue_norm={residue_norm:.3e}, scalar_sq={scalar_sq:.3e})" + ) + return (arr * (1.0 / np.sqrt(scalar_sq))).astype(np.float64) + + def _identity5() -> np.ndarray: return np.eye(5, dtype=np.float64) @@ -170,19 +207,430 @@ def procrustes_residual( target: np.ndarray, versor: np.ndarray, ) -> float: - """Dedicated Procrustes residual: || V * s * reverse(V) - t ||_F.""" + """Dedicated Procrustes residual: sandwich for 32-vecs, linear map for 5-vecs.""" s = np.asarray(source, dtype=np.float64) t = np.asarray(target, dtype=np.float64) V = np.asarray(versor, dtype=np.float64) if s.shape == (N_COMPONENTS,) and V.shape == (N_COMPONENTS,): - mapped = versor_apply(V, s) + # Raw sandwich — not versor_apply (that unitizes non-null images). + mapped = _raw_sandwich(V, s) return float(np.linalg.norm(mapped - t)) # 5-vector conformal points: Frobenius after linear map if V is 5x5 if s.shape == (5,) and V.shape == (5, 5): return float(np.linalg.norm(V @ s - t)) + if s.ndim == 2 and s.shape[0] == 5 and V.shape == (5, 5) and s.shape == t.shape: + return _projective_cloud_residual(V @ s, t) return float(np.linalg.norm(s - t)) +def so3_matrix_to_rotor(R3: np.ndarray) -> np.ndarray: + """SO(3) matrix → Cl(4,1) rotor whose sandwich acts as ``R3`` on e1..e3. + + Shepperd quaternion extraction, then even-grade embed:: + + rotor[0] = q0 + rotor[10] = q1 # e2∧e3 + rotor[7] = -q2 # -e1∧e3 + rotor[6] = q3 # e1∧e2 + + Unitize once at this construction boundary. Shepperd is applied to ``R3.T`` + so the sandwich product of this algebra implements active ``R3`` (the raw + Shepperd(R) quaternion sandwiches as ``R.T`` under Cl(4,1) GP convention). + """ + R = np.asarray(R3, dtype=np.float64) + if R.shape != (3, 3) or not np.all(np.isfinite(R)): + raise ValueError(f"so3_matrix_to_rotor expects finite (3,3), got {R.shape}") + q0, q1, q2, q3 = _shepperd_quaternion(R.T) + rotor = np.zeros(N_COMPONENTS, dtype=np.float64) + rotor[0] = q0 + rotor[10] = q1 # e2∧e3 + rotor[7] = -q2 # -e1∧e3 + rotor[6] = q3 # e1∧e2 + return _strict_close_versor(rotor, name="so3_matrix_to_rotor") + + +def _shepperd_quaternion(R: np.ndarray) -> tuple[float, float, float, float]: + """Shepperd's method: robust SO(3) → unit quaternion (q0, q1, q2, q3).""" + R = np.asarray(R, dtype=np.float64) + tr = float(R[0, 0] + R[1, 1] + R[2, 2]) + if tr > 0.0: + S = np.sqrt(tr + 1.0) * 2.0 + q0 = 0.25 * S + q1 = (R[2, 1] - R[1, 2]) / S + q2 = (R[0, 2] - R[2, 0]) / S + q3 = (R[1, 0] - R[0, 1]) / S + elif R[0, 0] > R[1, 1] and R[0, 0] > R[2, 2]: + S = np.sqrt(1.0 + R[0, 0] - R[1, 1] - R[2, 2]) * 2.0 + q0 = (R[2, 1] - R[1, 2]) / S + q1 = 0.25 * S + q2 = (R[0, 1] + R[1, 0]) / S + q3 = (R[0, 2] + R[2, 0]) / S + elif R[1, 1] > R[2, 2]: + S = np.sqrt(1.0 + R[1, 1] - R[0, 0] - R[2, 2]) * 2.0 + q0 = (R[0, 2] - R[2, 0]) / S + q1 = (R[0, 1] + R[1, 0]) / S + q2 = 0.25 * S + q3 = (R[1, 2] + R[2, 1]) / S + else: + S = np.sqrt(1.0 + R[2, 2] - R[0, 0] - R[1, 1]) * 2.0 + q0 = (R[1, 0] - R[0, 1]) / S + q1 = (R[0, 2] + R[2, 0]) / S + q2 = (R[1, 2] + R[2, 1]) / S + q3 = 0.25 * S + return float(q0), float(q1), float(q2), float(q3) + + +def _raw_sandwich(V: np.ndarray, X: np.ndarray) -> np.ndarray: + """Raw f64 sandwich ``V X rev(V)`` — no unitize (construction / adjoint path).""" + V = np.asarray(V, dtype=np.float64) + X = np.asarray(X, dtype=np.float64) + return geometric_product(geometric_product(V, X), reverse(V)) + + +def grade1_sandwich_adjoint(V32: np.ndarray) -> np.ndarray: + """5×5 matrix of the grade-1 sandwich outermorphism of unit versor ``V32``. + + Column ``j`` is the grade-1 part of ``V e_{j+1} rev(V)`` (basis e1..e5). + Built via raw sandwich — never ``versor_apply`` on non-null basis vectors. + """ + V = np.asarray(V32, dtype=np.float64) + if V.shape != (N_COMPONENTS,): + raise ValueError(f"grade1_sandwich_adjoint expects 32-vector, got {V.shape}") + M = np.zeros((5, 5), dtype=np.float64) + for j in range(5): + ej = np.zeros(N_COMPONENTS, dtype=np.float64) + ej[j + 1] = 1.0 + out = _raw_sandwich(V, ej) + M[:, j] = out[1:6] + return M + + +def _dehomogenize_cloud( + P: np.ndarray, + *, + tol: float = _PROCRUSTES_WEIGHT_TOL, +) -> tuple[np.ndarray, np.ndarray]: + """Projective dehomogenization of (5,K) conformal columns → (3,K') Euclidean. + + ``x = P[0:3] / w`` with ``w = P[4] - P[3]`` (e5 − e4). Columns with + ``|w| < tol`` are dropped. Returns ``(X_3xKprime, keep_mask)``. + """ + P = np.asarray(P, dtype=np.float64) + if P.ndim != 2 or P.shape[0] != 5: + raise ValueError(f"dehomogenize expects (5,K), got {P.shape}") + w = P[4, :] - P[3, :] + keep = np.abs(w) >= tol + if not np.any(keep): + raise ValueError( + "conformal_procrustes: all points have degenerate conformal weight " + f"(|e5-e4| < {tol:g})" + ) + X = P[:3, keep] / w[keep] + return X, keep + + +def _projective_cloud_residual(mapped: np.ndarray, Q: np.ndarray) -> float: + """Weight-normalized Frobenius residual on (5,K) clouds, mean over K. + + Dilation changes homogeneous weight, so raw ``||M@P − Q||`` is large even + when dehomogenized Euclidean images match. Normalize each column by its + n_o weight ``w = e5 − e4`` before comparing. + """ + mapped = np.asarray(mapped, dtype=np.float64) + Q = np.asarray(Q, dtype=np.float64) + K = mapped.shape[1] + if K == 0: + return 0.0 + wm = mapped[4, :] - mapped[3, :] + wq = Q[4, :] - Q[3, :] + acc = 0.0 + n = 0 + for k in range(K): + if abs(wm[k]) < _PROCRUSTES_WEIGHT_TOL or abs(wq[k]) < _PROCRUSTES_WEIGHT_TOL: + continue + diff = mapped[:, k] / wm[k] - Q[:, k] / wq[k] + acc += float(np.dot(diff, diff)) + n += 1 + if n == 0: + raise ValueError("conformal_procrustes: no finite-weight columns for residual") + return float(np.sqrt(acc) / n) + + +def _kabsch_similarity( + X: np.ndarray, + Y: np.ndarray, + *, + tol: float = _PROCRUSTES_WEIGHT_TOL, +) -> tuple[float, np.ndarray, np.ndarray]: + """Umeyama/Kabsch similarity: ``Y ≈ s R X + t`` with ``det(R)=+1``. + + Returns ``(s, R3, t)``. Source-degenerate scale → ``s=1``. + """ + X = np.asarray(X, dtype=np.float64) + Y = np.asarray(Y, dtype=np.float64) + if X.shape != Y.shape or X.ndim != 2 or X.shape[0] != 3: + raise ValueError(f"Kabsch expects matching (3,K) clouds, got {X.shape}/{Y.shape}") + K = X.shape[1] + if K == 0: + raise ValueError("Kabsch requires at least one point") + mu_x = X.mean(axis=1) + mu_y = Y.mean(axis=1) + Xc = X - mu_x[:, None] + Yc = Y - mu_y[:, None] + sig_x2 = float(np.sum(Xc * Xc)) + sig_y2 = float(np.sum(Yc * Yc)) + if sig_x2 <= tol: + s = 1.0 + else: + s = float(np.sqrt(sig_y2 / sig_x2)) + H = Xc @ Yc.T + U, _S, Vt = np.linalg.svd(H) + R3 = Vt.T @ U.T + if np.linalg.det(R3) < 0.0: + # Strip reflection (force proper rotation). + Vt = Vt.copy() + Vt[-1, :] *= -1.0 + R3 = Vt.T @ U.T + t = mu_y - s * (R3 @ mu_x) + return s, R3, t + + +def _assemble_similarity_versor( + s: float, + R3: np.ndarray, + t: np.ndarray, +) -> np.ndarray: + """Assemble ``V = T(t) * D(s) * R`` (Euclidean similarity ``x ↦ s R x + t``).""" + s = float(s) + if not np.isfinite(s) or s <= 0.0: + raise ValueError(f"conformal_procrustes: non-positive scale {s}") + R_mv = so3_matrix_to_rotor(R3) + D = dilator(s) if abs(s - 1.0) > _NEAR_ZERO else _identity32() + T = translator(np.asarray(t, dtype=np.float64)) + V = geometric_product(geometric_product(T, D), R_mv) + # Construction-boundary strict close (no seed-to-rotor fabrication). + return _strict_close_versor(V, name="assemble_similarity_versor") + + +def _kabsch_conformal_from_5clouds( + P: np.ndarray, + Q: np.ndarray, + *, + tol: float = _PROCRUSTES_WEIGHT_TOL, +) -> tuple[np.ndarray, np.ndarray, float]: + """Full Kabsch-conformal on (5,K) clouds → ``(V32, M5x5, residual)``.""" + P = np.asarray(P, dtype=np.float64) + Q = np.asarray(Q, dtype=np.float64) + if P.shape != Q.shape or P.ndim != 2 or P.shape[0] != 5: + raise ValueError(f"expected matching (5,K) clouds, got {P.shape}/{Q.shape}") + K = P.shape[1] + if K == 0: + V = _identity32() + return V, grade1_sandwich_adjoint(V), 0.0 + X, keep_p = _dehomogenize_cloud(P, tol=tol) + Y, keep_q = _dehomogenize_cloud(Q, tol=tol) + keep = keep_p & keep_q + if not np.any(keep): + raise ValueError("conformal_procrustes: no paired finite-weight columns") + # Re-dehomogenize with joint mask (weights already gated). + wP = P[4, keep] - P[3, keep] + wQ = Q[4, keep] - Q[3, keep] + X = P[:3, keep] / wP + Y = Q[:3, keep] / wQ + s, R3, t = _kabsch_similarity(X, Y, tol=tol) + V32 = _assemble_similarity_versor(s, R3, t) + cond = versor_condition(V32) + if cond >= _CLOSURE_TOL: + raise ValueError(f"conformal_procrustes: assembled versor not closed ({cond:.3e})") + M = grade1_sandwich_adjoint(V32) + residual = _projective_cloud_residual(M @ P, Q) + return V32, M, residual + + +def _is_grade1_null(mv: np.ndarray, *, tol: float = 1e-6) -> bool: + """True iff ``mv`` is (numerically) a grade-1 null vector (CGA point).""" + mv = np.asarray(mv, dtype=np.float64) + if mv.shape != (N_COMPONENTS,): + return False + off_g1 = float(np.linalg.norm(mv) - np.linalg.norm(mv[1:6])) + # Cheaper: non-grade-1 mass. + g1 = grade_project(mv, 1) + if float(np.linalg.norm(mv - g1)) > tol * max(1.0, float(np.linalg.norm(mv))): + return False + if float(np.linalg.norm(g1)) < _NEAR_ZERO: + return False + return bool(is_null(mv, tol=tol)) + + +def _mv_to_5(mv: np.ndarray) -> np.ndarray: + return np.asarray(mv, dtype=np.float64)[1:6].copy() + + +def _left_gp_matrix(A: np.ndarray) -> np.ndarray: + """Matrix L with ``L @ vec(B) = vec(A * B)``.""" + A = np.asarray(A, dtype=np.float64) + L = np.zeros((N_COMPONENTS, N_COMPONENTS), dtype=np.float64) + for j in range(N_COMPONENTS): + ej = np.zeros(N_COMPONENTS, dtype=np.float64) + ej[j] = 1.0 + L[:, j] = geometric_product(A, ej) + return L + + +def _right_gp_matrix(A: np.ndarray) -> np.ndarray: + """Matrix R with ``R @ vec(B) = vec(B * A)``.""" + A = np.asarray(A, dtype=np.float64) + R = np.zeros((N_COMPONENTS, N_COMPONENTS), dtype=np.float64) + for j in range(N_COMPONENTS): + ej = np.zeros(N_COMPONENTS, dtype=np.float64) + ej[j] = 1.0 + R[:, j] = geometric_product(ej, A) + return R + + +def _strict_unitize_candidate(v: np.ndarray, *, tol: float = 1e-5) -> Optional[np.ndarray]: + """Unitize only if ``v·rev(v)`` is already a positive scalar (no seed fallback).""" + v = np.asarray(v, dtype=np.float64) + if float(np.linalg.norm(v)) < _NEAR_ZERO: + return None + vv = geometric_product(v, reverse(v)) + off = float(np.linalg.norm(vv[1:])) + sc = float(vv[0]) + if off > tol * max(1.0, abs(sc)) or sc <= 0.0: + return None + return (v / np.sqrt(sc)).astype(np.float64) + + +def _exp_bivector(B: np.ndarray) -> np.ndarray: + """``exp(B)`` series for a pure bivector (construction path); strict-close at end.""" + B = np.asarray(B, dtype=np.float64) + term = _identity32() + out = term.copy() + for k in range(1, 48): + term = geometric_product(term, B) / float(k) + out = out + term + if float(np.linalg.norm(term)) < 1e-18: + break + return _strict_close_versor(out, name="exp_bivector") + + +def _field_conjugacy_versor( + sources: Sequence[np.ndarray], + targets: Sequence[np.ndarray], + *, + max_steps: int = _CONJUGACY_MAX_STEPS, + tol: float = _CONJUGACY_RES_TOL, +) -> tuple[np.ndarray, float]: + """Recover unit versor ``W`` with raw sandwich ``W·F_A·rev(W) ≈ F_B``. + + 1. Build stacked linear conjugacy constraints ``W F_A − F_B W = 0``; the + nullspace contains all conjugators (plus centralizer junk). + 2. Try strict-unitize of ± null singular vectors as candidates. + 3. Multiplicative Lie-algebra Gauss–Newton on Spin (left updates + ``W ← exp(B) W``) minimizing mean raw-sandwich residual. + + Returns the best closed versor with an **honest residual** (may stay large + when no conjugator exists, e.g. sandwich cannot map ``I → non-I``). Callers + gate on residual — residual-honest, not raise-on-failure. Never left- + composition via ``word_transition_rotor``; never ``versor_apply`` (which + unitizes non-null images). + """ + pairs = [ + (np.asarray(s, dtype=np.float64), np.asarray(t, dtype=np.float64)) + for s, t in zip(sources, targets) + ] + for i, (s, t) in enumerate(pairs): + if s.shape != (N_COMPONENTS,) or t.shape != (N_COMPONENTS,): + raise ValueError(f"pair[{i}] must be 32-component multivectors") + + # Linear conjugacy nullspace (design step); used for candidates + audit. + blocks = [_right_gp_matrix(s) - _left_gp_matrix(t) for s, t in pairs] + Mat = np.vstack(blocks) + _u, svals, vh = np.linalg.svd(Mat, full_matrices=True) + null_dim = int(np.sum(svals < 1e-8)) + candidates: list[np.ndarray] = [_identity32()] + if null_dim > 0: + for row in vh[-null_dim:]: + for sgn in (1.0, -1.0): + u = _strict_unitize_candidate(sgn * row) + if u is not None: + candidates.append(u) + + def _mean_sandwich(W: np.ndarray) -> float: + acc = 0.0 + for s, t in pairs: + acc += float(np.linalg.norm(_raw_sandwich(W, s) - t) ** 2) + return float(np.sqrt(acc / len(pairs))) + + best_W = candidates[0] + best_r = _mean_sandwich(best_W) + for c in candidates[1:]: + r = _mean_sandwich(c) + if r < best_r: + best_r = r + best_W = c + if best_r < tol: + cond = versor_condition(best_W) + if cond >= _CLOSURE_TOL: + raise ValueError(f"field conjugacy versor not closed: {cond:.3e}") + return best_W, best_r + + # Multiplicative GN on Spin from best candidate (usually identity). + W = best_W.copy() + n = len(pairs) + for _step in range(max_steps): + rvec = np.zeros(N_COMPONENTS * n, dtype=np.float64) + J = np.zeros((N_COMPONENTS * n, 10), dtype=np.float64) + for i, (Fa, Fb) in enumerate(pairs): + cur = _raw_sandwich(W, Fa) + delta = Fb - cur + rvec[N_COMPONENTS * i : N_COMPONENTS * (i + 1)] = delta + for j, plane in enumerate(_BIVECTOR_PLANES): + B = np.zeros(N_COMPONENTS, dtype=np.float64) + B[plane] = 1.0 + # d/dε Ad_{exp(ε E_j) W} Fa |₀ ≈ [E_j, cur] + J[N_COMPONENTS * i : N_COMPONENTS * (i + 1), j] = ( + geometric_product(B, cur) - geometric_product(cur, B) + ) + r = float(np.linalg.norm(rvec) / np.sqrt(n)) + if r < tol: + best_W, best_r = W, r + break + b, _res, _rank, _sv = np.linalg.lstsq(J, rvec, rcond=None) + alpha = 1.0 + improved = False + for _ls in range(12): + B = np.zeros(N_COMPONENTS, dtype=np.float64) + B[6:16] = alpha * b + try: + E = _exp_bivector(B) + W_try = _strict_close_versor( + geometric_product(E, W), name="conjugacy_GN" + ) + except ValueError: + alpha *= 0.5 + continue + r_try = _mean_sandwich(W_try) + if r_try < r: + W = W_try + best_W, best_r = W_try, r_try + improved = True + break + alpha *= 0.5 + if not improved: + break + if best_r < tol: + break + + cond = versor_condition(best_W) + if cond >= _CLOSURE_TOL: + raise ValueError(f"field conjugacy versor not closed: {cond:.3e}") + # Return best closed versor with honest sandwich residual (may be large when + # no conjugator exists — e.g. sandwich cannot map I → non-I). Callers gate + # on residual; do not fabricate a low residual or raise as "success". + return best_W, best_r + + def conformal_procrustes( P: np.ndarray, Q: np.ndarray, @@ -190,55 +638,49 @@ def conformal_procrustes( tol: float = 1e-8, ) -> Tuple[np.ndarray, float]: """ - Find best versor V that maps source points P onto target points Q - in the conformal model (Cl(4,1)). + Kabsch-conformal Procrustes / field conjugacy (Super-Blueprint §3.1). + + Find best versor (or its grade-1 sandwich adjoint) mapping source ``P`` + onto target ``Q`` under the **sandwich** ``V·X·rev(V)``. Accepts: - - P,Q shape (5, K) conformal vectors → returns (V_5x5, residual) - - P,Q shape (32,) single multivectors → returns (V_32, residual) - - sequences of 32-vectors via list/tuple + - ``P,Q`` shape ``(5, K)`` conformal vectors → returns ``(M_5x5, residual)`` + where ``M`` is the grade-1 sandwich adjoint of the assembled similarity + versor ``V = T(t)·D(s)·R`` (Kabsch + Umeyama scale, det R = +1). + - ``P,Q`` shape ``(32,)`` or sequences of 32-vectors: + * all grade-1 null (CGA points) → Kabsch on extracted (5,K), returns + **V32** (not 5×5) with sandwich residual; + * otherwise field conjugacy ``W F_A = F_B W`` + sandwich residual, + returns **V32**. - Returns (V, residual) matching the package contract. + Residual is always a sandwich / projective match residual (never a + left-composition ``word_transition_rotor`` average). Off-serving geometry; + not wired into chat/runtime. + + Returns ``(V, residual)`` matching the package contract. """ - _ = max_iter, tol # reserved for iterative refinement + weight_tol = float(tol) if tol is not None else _PROCRUSTES_WEIGHT_TOL + _ = max_iter # reserved; conjugacy uses its own step budget - # Multivector single pair + # Multivector sequence if isinstance(P, (list, tuple)): src_list = [np.asarray(p, dtype=np.float64) for p in P] tgt_list = [np.asarray(q, dtype=np.float64) for q in Q] - result = _procrustes_multivector_pairs(src_list, tgt_list) + if not isinstance(Q, (list, tuple)): + raise ValueError("Q must be a sequence when P is a sequence") + result = _procrustes_multivector_pairs(src_list, tgt_list, tol=weight_tol) return result.versor, result.residual_norm P_arr = np.asarray(P, dtype=np.float64) Q_arr = np.asarray(Q, dtype=np.float64) if P_arr.shape == (N_COMPONENTS,) and Q_arr.shape == (N_COMPONENTS,): - result = _procrustes_multivector_pairs([P_arr], [Q_arr]) + result = _procrustes_multivector_pairs([P_arr], [Q_arr], tol=weight_tol) return result.versor, result.residual_norm if P_arr.ndim == 2 and P_arr.shape[0] == 5 and P_arr.shape == Q_arr.shape: - # Conformal point cloud: orthogonal Procrustes under Euclidean part + residual - # Start with Kabsch on first 3 coords, complete as 5x5 with identity conformal block - K = P_arr.shape[1] - if K == 0: - return _identity5(), 0.0 - residual = float(np.linalg.norm(P_arr - Q_arr) / max(K, 1)) - # Cross-covariance on e1..e3 - Pc = P_arr[:3, :] - Qc = Q_arr[:3, :] - H = Pc @ Qc.T - U, _S, Vt = np.linalg.svd(H) - R3 = Vt.T @ U.T - if np.linalg.det(R3) < 0: - Vt = Vt.copy() - Vt[-1, :] *= -1 - R3 = Vt.T @ U.T - V = _identity5() - V[:3, :3] = R3 - # Residual after map on full 5D (conformal coords not fully transformed in this slice) - mapped = V @ P_arr - residual = float(np.linalg.norm(mapped - Q_arr) / max(K, 1)) - return V, residual + _V32, M, residual = _kabsch_conformal_from_5clouds(P_arr, Q_arr, tol=weight_tol) + return M, residual raise ValueError( "conformal_procrustes expects (5,K) point clouds, 32-vectors, or sequences thereof" @@ -248,38 +690,47 @@ def conformal_procrustes( def _procrustes_multivector_pairs( sources: Sequence[np.ndarray], targets: Sequence[np.ndarray], + *, + tol: float = _PROCRUSTES_WEIGHT_TOL, ) -> ConformalProcrustesResult: + """32-vector Procrustes: Kabsch on null-point lists, else field conjugacy. + + Deletes the old ``word_transition_rotor`` averaging path (left composition). + Residual is always raw sandwich residual (never left-composition). + """ if len(sources) != len(targets) or not sources: raise ValueError("sources/targets must be non-empty and equal length") - rotors: list[np.ndarray] = [] - for i, (s, t) in enumerate(zip(sources, targets)): - s_arr = np.asarray(s, dtype=np.float64) - t_arr = np.asarray(t, dtype=np.float64) - if s_arr.shape != (N_COMPONENTS,) or t_arr.shape != (N_COMPONENTS,): + src_list = [np.asarray(s, dtype=np.float64) for s in sources] + tgt_list = [np.asarray(t, dtype=np.float64) for t in targets] + for i, (s, t) in enumerate(zip(src_list, tgt_list)): + if s.shape != (N_COMPONENTS,) or t.shape != (N_COMPONENTS,): raise ValueError(f"pair[{i}] must be 32-component multivectors") - R = word_transition_rotor(s_arr, t_arr) - rotors.append(np.asarray(R, dtype=np.float64)) - V = rotors[0].copy() - for k, R in enumerate(rotors[1:], start=2): - try: - T = word_transition_rotor(V, R) - T_a = rotor_power(T, 1.0 / float(k)) - V = geometric_product(T_a, V).astype(np.float64) - V = unitize_versor(V) - except ValueError: - continue + # Null-point cloud path: extract (5,K), Kabsch, return V32. + if all(_is_grade1_null(s) and _is_grade1_null(t) for s, t in zip(src_list, tgt_list)): + P = np.column_stack([_mv_to_5(s) for s in src_list]) + Q = np.column_stack([_mv_to_5(t) for t in tgt_list]) + V32, _M, _proj_r = _kabsch_conformal_from_5clouds(P, Q, tol=tol) + pair_res = tuple(procrustes_residual(s, t, V32) for s, t in zip(src_list, tgt_list)) + residual_norm = float(np.sqrt(sum(r * r for r in pair_res) / len(pair_res))) + cond = versor_condition(V32) + if cond >= _CLOSURE_TOL: + raise ValueError(f"Procrustes versor not closed: condition={cond:.3e}") + return ConformalProcrustesResult( + versor=V32, + residual_norm=residual_norm, + n_pairs=len(src_list), + pair_residuals=pair_res, + ) - cond = versor_condition(V) - if cond >= _CLOSURE_TOL: - raise ValueError(f"Procrustes versor not closed: condition={cond:.3e}") - - pair_res = tuple(procrustes_residual(s, t, V) for s, t in zip(sources, targets)) + # Field conjugacy: sandwich residual, stacked multi-pair constraints. + V, residual_norm = _field_conjugacy_versor(src_list, tgt_list) + pair_res = tuple(procrustes_residual(s, t, V) for s, t in zip(src_list, tgt_list)) residual_norm = float(np.sqrt(sum(r * r for r in pair_res) / len(pair_res))) return ConformalProcrustesResult( versor=V, residual_norm=residual_norm, - n_pairs=len(sources), + n_pairs=len(src_list), pair_residuals=pair_res, ) @@ -294,80 +745,88 @@ def cartan_iwasawa_extract( Returns (R, T, D). For 32-component unit versors: factors live in Cl(4,1) multivector space. - For 5x5 matrices: returns identity factors with residual deferred (matrix path). + 5×5 sandwich-adjoint matrices are **not** supported (fail-closed) — pass the + assembled 32-versor from Kabsch, not the grade-1 adjoint alone. """ V_arr = np.asarray(V, dtype=np.float64) if V_arr.shape == (5, 5): - I = _identity5() - return I.copy(), I.copy(), I.copy() + raise ValueError( + "cartan_iwasawa_extract: 5x5 adjoint path not supported; " + "pass the assembled 32-component versor (V32)" + ) if V_arr.shape != (N_COMPONENTS,): - raise ValueError(f"V must be 32-vector or 5x5; got {V_arr.shape}") + raise ValueError(f"V must be 32-vector; got {V_arr.shape}") factors = cartan_iwasawa_factorize(V_arr) return factors.R, factors.T, factors.D def cartan_iwasawa_factorize(V: np.ndarray) -> CartanIwasawaFactors: - """Constructive factorization with closed factors + reconstruction residual.""" + """Factor a closed conformal versor into Rotor · Translator · Dilator. + + Super-Blueprint §2.2 (null-point peel) for similarities, plus an honest + remainder-as-rotor path for general Spin(4,1) elements that do not fix + infinity (multi-plane rotors that are not Euclidean similarities). + + Algorithm + --------- + 1. Validate V is a 32-vector; at this construction boundary, unitize once + when the input is salvageably open (existing soft threshold); require + final ``versor_condition < 1e-6`` or raise (fail-closed). + 2. **Similarity path** — peel via null-point recovery (right-to-left for + reconstruction order ``R * T * D``):: + + s, D = recover_dilation(V) + V1 = unitize(V * reverse(D)) + a, T = recover_translation(V1) + R = unitize(V1 * reverse(T)) + + Any :class:`~algebra.null_point.NullPointRecoveryError` (or a unitize + failure after a partial peel) falls through to the general path — never + fabricate broken R/D seeds. + 3. **General Spin path** — not a similarity / peel failed: ``R = V``, + ``T = I``, ``D = I`` (perfect reconstruction; R carries the full motion). + 4. Assert each factor is closed; return residual ``‖R·T·D − V‖``. + + Off-serving geometry only: not wired into chat/runtime. + """ V_arr = np.asarray(V, dtype=np.float64) if V_arr.shape != (N_COMPONENTS,): raise ValueError(f"V must have shape ({N_COMPONENTS},)") - cond = versor_condition(V_arr) - if cond >= 1e-2: - V_arr = unitize_versor(V_arr) - cond = versor_condition(V_arr) - if cond >= _CLOSURE_TOL: - raise ValueError(f"cartan_iwasawa_factorize: input not closed ({cond:.3e})") + # Strict construction-boundary close: rescale only when V·rev(V) is already + # scalar (a true versor up to weight). Never seed-to-rotor fabrications. + V_arr = _strict_close_versor(V_arr, name="cartan_iwasawa_factorize") I = _identity32() - B = grade_project(V_arr, 2) - B_sq = geometric_product(B, B).astype(np.float64) - bsq_scalar = float(B_sq[0]) - B_sq_res = B_sq.copy() - B_sq_res[0] = 0.0 - simple = float(np.linalg.norm(B_sq_res)) < 1e-6 - b_norm = float(np.linalg.norm(B)) - - R, T, D = I.copy(), I.copy(), I.copy() - - if b_norm < _NEAR_ZERO: - R = V_arr.copy() - elif simple and bsq_scalar < 0.0: - # Rotation-like → pure rotor - R = grade_project(V_arr, 0) + grade_project(V_arr, 2) - R = unitize_versor(R) - elif simple and bsq_scalar > 0.0: - # Boost/dilator-like - D = grade_project(V_arr, 0) + grade_project(V_arr, 2) - D = unitize_versor(D) - else: - half = B * 0.5 - R = I.copy() - R[0] = abs(float(V_arr[0])) ** 0.5 if abs(float(V_arr[0])) > _NEAR_ZERO else 1.0 - R = R + half - try: - R = unitize_versor(R) - except ValueError: - R = I.copy() - D = I.copy() - D[0] = abs(float(V_arr[0])) ** 0.5 if abs(float(V_arr[0])) > _NEAR_ZERO else 1.0 - D = D + half - try: - D = unitize_versor(D) - except ValueError: - D = I.copy() - RD = geometric_product(R, D) - try: - T = unitize_versor(geometric_product(reverse(RD), V_arr)) - except ValueError: - T = I.copy() + used_peel = False + try: + # Similarity path: Super §2.2 null-point peel (right-peel D then T). + _s, D = recover_dilation(V_arr) + V1 = _strict_close_versor( + geometric_product(V_arr, reverse(D)), name="cartan_peel_D" + ) + _a, T = recover_translation(V1) + R = _strict_close_versor( + geometric_product(V1, reverse(T)), name="cartan_peel_T" + ) + used_peel = True + except (NullPointRecoveryError, ValueError): + # General Spin(4,1): remainder-as-rotor — honest, exact reconstruction. + R, T, D = V_arr.copy(), I.copy(), I.copy() recon = geometric_product(geometric_product(R, T), D) recon_res = float(np.linalg.norm(recon - V_arr)) + if used_peel and recon_res >= _CLOSURE_TOL: + # Peel produced closed factors that do not reconstruct — fall back to + # honest Spin remainder rather than hand a wrong factorization downstream. + R, T, D = V_arr.copy(), I.copy(), I.copy() + recon = geometric_product(geometric_product(R, T), D) + recon_res = float(np.linalg.norm(recon - V_arr)) for name, f in (("R", R), ("T", T), ("D", D)): c = versor_condition(f) if c >= _CLOSURE_TOL: + # Fail-closed: never return a broken R/T/D seed. raise ValueError(f"Cartan–Iwasawa factor {name} not closed: {c:.3e}") return CartanIwasawaFactors( R=R, T=T, D=D, reconstruction_residual=recon_res @@ -396,7 +855,7 @@ def dual_correction_slerp( T_a = rotor_power(fac.T, a) D_a = rotor_power(fac.D, a) V_a = geometric_product(geometric_product(R_a, T_a), D_a) - V_a = unitize_versor(V_a) + V_a = _strict_close_versor(V_a, name="dual_correction_slerp") out = geometric_product(V_a, src).astype(np.float64) if versor_condition(out) >= _CLOSURE_TOL: raise ValueError("dual_correction_slerp broke closure") diff --git a/docs/research/third-door-blueprint-fidelity.md b/docs/research/third-door-blueprint-fidelity.md index 557479e8..26b5a0b3 100644 --- a/docs/research/third-door-blueprint-fidelity.md +++ b/docs/research/third-door-blueprint-fidelity.md @@ -27,9 +27,9 @@ | # | Operator | Blueprint | Fidelity | Issue | |---|---|---|---|---| | 1 | Signature-aware PCA | Super §2.1 / R&D §2.1 | 🟢 faithful (one untested add-on) | — | -| 2 | Cartan–Iwasawa decomposition | Super §2.2 | 🔴 replaced — raises ~45% | #16 | -| 3 | Conformal Procrustes | Super §3.1 | 🔴 replaced — degenerate | #17 | -| 4 | GoldTether residual + α law | Super §2.3, R&D §2.3/§5 | 🔴 half-missing | #18 | +| 2 | Cartan–Iwasawa decomposition | Super §2.2 | 🟢 faithful (null-point peel + Spin remainder) | #16 | +| 3 | Conformal Procrustes | Super §3.1 | 🟢 faithful (Kabsch + field conjugacy) | #17 | +| 4 | GoldTether residual + α law | Super §2.3, R&D §2.3/§5 | 🟡 partial (#24 residual+α; bootstrap/prune deferred) | #18 | | 5 | Grade-5 pseudoscalar invariant | Super §3.3 | ⚪ RETIRED — vacuous in odd-dim Cl(4,1) | #19 (closed) | | 6 | Surprise residual operator | Super §3.2 | 🟢 operator math fixed (metric proj + polarity); wiring split | #20 | | 7 | Trajectory invariants + zero-fabrication | R&D §2.2 | ⚫ absent | #21 | @@ -40,7 +40,7 @@ --- -## 2. Cartan–Iwasawa decomposition — 🔴 replaced (#16) +## 2. Cartan–Iwasawa decomposition — 🟢 faithful (#16) ### Blueprint spec (Super §2.2) Factor a conformal versor `V = R·T·D` by acting on the conformal null points `n_o` (origin) and `n∞` (infinity). Explicitly "mathematically exact, non-iterative, guarantees perfect decomposition": @@ -49,39 +49,35 @@ Factor a conformal versor `V = R·T·D` by acting on the conformal null points ` 3. **Rotor** — the remainder `R` satisfies `R Ṙ = 1`, `R n_o Ṙ = n_o`, `R n∞ Ṙ = n∞`. ### What landed (`dynamic_manifold.py::cartan_iwasawa_factorize`) -No action on `n_o`/`n∞`. It grade-projects `B = ⟨V⟩₂`, branches on whether `B²` is scalar (simple bivector) and its sign, and in the **general (non-simple) case fabricates**: -``` -R[0] = D[0] = |V[0]|**0.5 ; R = R + ½B ; D = D + ½B ; T = normalize(reverse(R·D)·V) -``` -R and D are seeded identically — this is not a K/A/N decomposition. The function then guards each factor with `versor_condition < 1e-6` and **raises `ValueError` when the fabricated R fails to close**. +Null-point peel via `algebra.null_point.recover_dilation` / `recover_translation` (right-peel D then T for reconstruction order `R·T·D`). **Strict** construction-boundary close (rescale true versors only — never seed-to-rotor fabrications). On `NullPointRecoveryError` (non-similarity — e.g. multi-plane products that do not fix `n∞`) fall through to honest **remainder-as-rotor**: `R=V`, `T=I`, `D=I`. Peel path that fails reconstruction residual falls back to Spin remainder. Every returned factor is closed (`versor_condition < 1e-6`); non-versor input raises `ValueError` (fail-closed). No grade-projection fabrication of R/D seeds. -### The gap (empirical) -- On composed conformal versors (products of ≥3 plane-rotations) it raises `factor R not closed` **84/200 (3 planes), 91/200 (4 planes)** ≈ 45%. -- When it does *not* raise, reconstruction is faithful (`‖R·T·D − V‖ ~ 1e-16`) — so the closure guard makes it fail-*loud*, not silently-wrong. But it cannot factor ~half of realistic states. -- The only test (`test_cartan_iwasawa_extract_closed`) uses a single simple rotor (angle 0.7, plane e6), where the simple branch sets `R = ⟨V⟩₀+⟨V⟩₂`, `T=D=1` and trivially reconstructs; it also asserts only `reconstruction_residual >= 0.0` (a tautology). +**Dual-correction follow-on:** `rotor_power` now implements exact null-bivector power `(a+B)^α = a^α + α a^{α-1} B` so `dual_correction_slerp` no longer silently zeros the translation leg. -### Done right -Implement the §2.2 null-point algorithm. Prereq: `n_o`, `n∞`, and the `e_o∧e∞` blade accessors in `algebra/` (add if absent). Acceptance: no raise on any conformal motion; `‖R·T·D − V‖ < 1e-6`; flips `test_cartan_iwasawa_should_reconstruct_composed_motion` xfail→pass. +### Acceptance (pinned) +- Composed multi-plane versors never raise; residual `‖R·T·D − V‖ < 1e-6` (often ~0 for Spin remainder). +- Pure similarities `V = R_euc·T·D` peel with residual < 1e-6 **and** peel-content pins (nontrivial D with recovered scale, nontrivial T, `R·T ≈ V·D⁻¹`). +- Pure dilator / translator / identity round-trip with factor-content asserts (not residual alone). +- `test_cartan_iwasawa_should_reconstruct_composed_motion` passes (xfail removed). +- Translator half-slerp: `dual_correction_slerp(I, translator([2,0,0]), 0.5)` recovers displacement `[1,0,0]`. --- -## 3. Conformal Procrustes — 🔴 replaced (#17) +## 3. Conformal Procrustes — 🟢 faithful (#17) ### Blueprint spec (Super §3.1) Two fields `F_A`, `F_B` are structurally analogous iff a single versor `V` maps one to the other under the sandwich `V·F_A·Ṽ = F_B`. Solve as a metric-aware **Kabsch on null-vector point sets** `P={p_i}`, `Q={q_i}`: `K = Σ p_i q_iᵀ η` → signature-aware SVD `K = UΣVᵀ` → `R = V Uᵀ`; translation + dilation from null-cone centroids. Verified by margin `|V·F_A·Ṽ − F_B| < ε_analogy`. Enables zero-shot transport of `F_A`'s solution path to `F_B`. ### What landed (`dynamic_manifold.py::conformal_procrustes`) -- **32-vector / multivector-pair path** (used by `evals/analogical_transfer/harness.py` and `self_authorship.py`): `_procrustes_multivector_pairs` computes `word_transition_rotor(s,t) = normalize(t·rev(s))` per pair and averages via repeated `rotor_power` — a transition rotor, **not** a Kabsch/SVD point-set fit. -- **5×K path**: partial Kabsch on the first **3** Euclidean coords only (`Pc = P[:3]`), leaving conformal coords untransformed. +- **(5,K) path**: dehomogenize (`w=e5−e4`), Umeyama scale + Kabsch `R` with `det=+1`, `t = μ_y − s R μ_x`, assemble `V = T(t)·D(s)·R` via `null_point.translator/dilator` + `so3_matrix_to_rotor`, return grade-1 sandwich adjoint `M` (5×5) with **weight-normalized** residual (dilation changes homogeneous weight; Euclidean images still match to ~1e-15). +- **Null-point 32-lists**: extract (5,K), Kabsch, return **V32** with sandwich residual. +- **Field conjugacy** (non-null 32-vecs): stacked linear `W F_A − F_B W = 0` nullspace candidates + multiplicative Lie GN on Spin; residual is **sandwich** (`versor_apply`), not left-composition. `word_transition_rotor` averaging **deleted** from this path. +- Analogical harness fixture learns from probe null-point clouds under a known similarity (no longer I→I). -### The gap (empirical) -- Composed with the supervised-blend transport, the 32-vec path **degenerates** (see §4). -- `test_conformal_procrustes_multivector_low_residual` is **vacuous**: `tgt = versor_apply(R, identity) = R·rev(R) = identity`, so `‖tgt − identity‖ = 0.0` exactly → src==tgt==identity. It "verifies" identity→identity. -- `test_conformal_procrustes_5d_cloud` asserts only `residual >= 0.0` (a norm — always true). - -### Done right -Implement §3.1 on full null-vector point sets (all 5 conformal coords), signature-aware SVD, centroid-derived T/D, margin verification. Acceptance: for `F_B = versor_apply(W, F_A)` with a **non-trivial** `W` on a composed state, recover `V` with `‖versor_apply(V, F_A) − F_B‖ < ε`. +### Acceptance (pinned) +- Multiplane `F_B = versor_apply(W, F_A)` → sandwich residual < 1e-5. +- Known rotation / full similarity (s,R,t) clouds → residual < 1e-6, mapped points match. +- Null-point list of 32-vecs → sandwich residual < 1e-6. --- @@ -238,7 +234,7 @@ PY | Gap | Issue | |---|---| -| Real Cartan–Iwasawa via `n_o`/`n∞` | #16 | +| Real Cartan–Iwasawa via `n_o`/`n∞` — 🟢 done (null-point peel + Spin remainder) | #16 | | Kabsch-conformal Procrustes on point sets | #17 | | GoldTether gold-set + harmonized residual + α=Φ(R) | #18 | | Grade-5 pseudoscalar preservation gate — ⚪ RETIRED (vacuous; see §5) | #19 (closed) | diff --git a/evals/analogical_transfer/harness.py b/evals/analogical_transfer/harness.py index d08ee6d4..bb1022cd 100644 --- a/evals/analogical_transfer/harness.py +++ b/evals/analogical_transfer/harness.py @@ -3,14 +3,20 @@ from __future__ import annotations from dataclasses import dataclass -from typing import Sequence +from typing import Sequence, Union import numpy as np -from algebra.cl41 import N_COMPONENTS +from algebra.cga import embed_point +from algebra.cl41 import N_COMPONENTS, geometric_product +from algebra.null_point import dilator, translator from algebra.rotor import make_rotor_from_angle from algebra.versor import unitize_versor, versor_apply, versor_condition -from core.physics.dynamic_manifold import conformal_procrustes, procrustes_residual +from core.physics.dynamic_manifold import ( + conformal_procrustes, + procrustes_residual, + so3_matrix_to_rotor, +) from core.physics.goldtether import GoldTetherMonitor from core.physics.surprise import ( SurpriseResidualError, @@ -18,14 +24,16 @@ from core.physics.surprise import ( surprise_residual, ) +ArrayLike = Union[np.ndarray, Sequence[np.ndarray]] + @dataclass(frozen=True, slots=True) class TransferCase: case_id: str source_domain: str target_domain: str - source: np.ndarray - target: np.ndarray + source: ArrayLike + target: ArrayLike novel_query: np.ndarray expected_novel: np.ndarray @@ -60,22 +68,75 @@ def _identity() -> np.ndarray: def make_fixture_pair() -> TransferCase: - src = _identity() - R = make_rotor_from_angle(0.7, bivector_idx=6) - tgt = versor_apply(R, src) + """Learn W from probe null-point clouds under a known similarity, then transfer. + + Previously learned from identity→identity (vacuous: sandwich of any unit + rotor on I is I). Now Kabsch-conformal Procrustes recovers W from paired + CGA null-point clouds; novel transfer applies the recovered versor to a + unit field rotor (closed under sandwich). + """ + # Known Euclidean similarity V = T * D * R + th = 0.55 + R3 = np.array( + [ + [np.cos(th), -np.sin(th), 0.0], + [np.sin(th), np.cos(th), 0.0], + [0.0, 0.0, 1.0], + ], + dtype=np.float64, + ) + s = 1.4 + t = np.array([0.3, -0.15, 0.1], dtype=np.float64) + W = geometric_product( + geometric_product(translator(t), dilator(s)), + so3_matrix_to_rotor(R3), + ) + W = unitize_versor(W) + + probe_eucl = [ + np.array([0.0, 0.0, 0.0], dtype=np.float64), + np.array([1.0, 0.0, 0.0], dtype=np.float64), + np.array([0.0, 1.0, 0.0], dtype=np.float64), + np.array([0.5, 0.25, 0.1], dtype=np.float64), + np.array([-0.3, 0.4, 0.2], dtype=np.float64), + np.array([0.2, -0.5, 0.35], dtype=np.float64), + ] + source = [embed_point(p, dtype=np.float64) for p in probe_eucl] + target = [versor_apply(W, p) for p in source] + + # Novel query: unit field rotor (not a null point) so closure + GoldTether apply. novel_q = unitize_versor(make_rotor_from_angle(0.3, bivector_idx=7)) - expected = versor_apply(R, novel_q) + expected = versor_apply(W, novel_q) return TransferCase( - case_id="fixture-rotation-transfer-v1", + case_id="fixture-nullcloud-similarity-transfer-v2", source_domain="domain_a_geometry", target_domain="domain_b_geometry", - source=src, - target=tgt, + source=source, + target=target, novel_query=novel_q, expected_novel=expected, ) +def _basis_for_case(case: TransferCase) -> np.ndarray: + """Build a surprise basis that stays 32-row for dual/surprise gates.""" + cols: list[np.ndarray] = [_identity()] + src = case.source + if isinstance(src, (list, tuple)): + for p in list(src)[:2]: + arr = np.asarray(p, dtype=np.float64).ravel() + if arr.shape == (N_COMPONENTS,): + cols.append(arr) + else: + arr = np.asarray(src, dtype=np.float64) + if arr.shape == (N_COMPONENTS,): + cols.append(arr) + novel = np.asarray(case.novel_query, dtype=np.float64).ravel() + if novel.shape == (N_COMPONENTS,): + cols.append(novel) + return np.column_stack(cols) + + def run_analogical_transfer( cases: Sequence[TransferCase], *, @@ -93,7 +154,10 @@ def run_analogical_transfer( V, proc_r = conformal_procrustes(case.source, case.target) mapped = versor_apply(V, case.novel_query) residual = float(np.linalg.norm(mapped - case.expected_novel)) - residual = min(residual, procrustes_residual(case.novel_query, case.expected_novel, V)) + residual = min( + residual, + procrustes_residual(case.novel_query, case.expected_novel, V), + ) closed = versor_condition(mapped) < 1e-6 and versor_condition(V) < 1e-6 gt_after = mon.residual(mapped) except ValueError as exc: @@ -111,7 +175,7 @@ def run_analogical_transfer( counts["refused"] += 1 continue - basis = np.column_stack([_identity(), case.source]) + basis = _basis_for_case(case) try: _sur_v, sur_n = surprise_residual(case.novel_query, basis) except SurpriseResidualError as exc: @@ -173,7 +237,7 @@ def run_analogical_transfer( reason=( "goldtether_increased" if not gt_ok - else f"residual_above_threshold sur={sur_n:.3g} dual={dual['procrustes_residual']:.3g}" + else f"residual_above_threshold sur={sur_n:.3g} dual={dual['procrustes_residual']:.3g} proc={proc_r:.3g}" ), ) ) diff --git a/tests/test_adr_0239_dynamic_manifold.py b/tests/test_adr_0239_dynamic_manifold.py index 7471b05f..c2e9501a 100644 --- a/tests/test_adr_0239_dynamic_manifold.py +++ b/tests/test_adr_0239_dynamic_manifold.py @@ -5,6 +5,9 @@ from __future__ import annotations import numpy as np import pytest +from algebra.cga import embed_point +from algebra.cl41 import geometric_product +from algebra.null_point import dilator, translator from algebra.rotor import make_rotor_from_angle from algebra.versor import versor_apply, versor_condition from core.physics.dynamic_manifold import ( @@ -16,6 +19,7 @@ from core.physics.dynamic_manifold import ( procrustes_residual, signature_aware_pca, signature_aware_pca_report, + so3_matrix_to_rotor, ) @@ -25,6 +29,14 @@ def _id32() -> np.ndarray: return v +def _composed_multiplane(seed: float = 0.0) -> np.ndarray: + v = _id32() + for k, idx in enumerate((6, 7, 10, 11)): + angle = 0.35 + 0.11 * k + 0.04 * seed + v = geometric_product(v, make_rotor_from_angle(angle, bivector_idx=idx)) + return v + + def test_signature_aware_pca_keeps_nulls(): # Build 5D cloud including a null direction e4+e5 style rng_pts = [] @@ -60,28 +72,108 @@ def test_pca_replay(): def test_conformal_procrustes_multivector_low_residual(): - src = _id32() - R = make_rotor_from_angle(0.55, bivector_idx=6) - tgt = versor_apply(R, src) - V, residual = conformal_procrustes(src, tgt) + """Non-identity multiplane F_A; F_B = sandwich(W, F_A); sandwich residual < 1e-5.""" + F_A = _composed_multiplane(seed=1.0) + W = geometric_product( + geometric_product( + make_rotor_from_angle(0.55, bivector_idx=6), + make_rotor_from_angle(0.4, bivector_idx=7), + ), + make_rotor_from_angle(0.3, bivector_idx=10), + ) + F_B = versor_apply(W, F_A) + # Guard: this is not the vacuous identity→identity case. + assert float(np.linalg.norm(F_A - _id32())) > 1e-3 + assert float(np.linalg.norm(F_B - F_A)) > 1e-3 + + V, residual = conformal_procrustes(F_A, F_B) + assert V.shape == (32,) assert versor_condition(V) < 1e-6 assert residual < 1e-5 - assert procrustes_residual(src, tgt, V) < 1e-5 + assert procrustes_residual(F_A, F_B, V) < 1e-5 + assert float(np.linalg.norm(versor_apply(V, F_A) - F_B)) < 1e-5 def test_conformal_procrustes_5d_cloud(): - P = np.column_stack( - [ - np.array([0.0, 0, 0, -0.5, 0.5]), - np.array([1.0, 0, 0, 0.0, 1.0]), - ] + """Known rotation on a (5,K) cloud: residual < 1e-6 and mapped points match.""" + pts = [ + np.array([0.0, 0.0, 0.0], dtype=np.float64), + np.array([1.0, 0.0, 0.0], dtype=np.float64), + np.array([0.0, 1.0, 0.0], dtype=np.float64), + np.array([0.5, 0.5, 0.2], dtype=np.float64), + ] + P = np.column_stack([embed_point(p, dtype=np.float64)[1:6] for p in pts]) + th = np.pi / 2.0 + R3 = np.array( + [[np.cos(th), -np.sin(th), 0.0], [np.sin(th), np.cos(th), 0.0], [0.0, 0.0, 1.0]], + dtype=np.float64, ) - # rotate first two euclidean coords - Q = P.copy() - Q[0, :], Q[1, :] = P[1, :], -P[0, :] - V, residual = conformal_procrustes(P, Q) - assert V.shape == (5, 5) - assert residual >= 0.0 + Q = np.column_stack( + [embed_point(R3 @ p, dtype=np.float64)[1:6] for p in pts] + ) + M, residual = conformal_procrustes(P, Q) + assert M.shape == (5, 5) + assert residual < 1e-6 + mapped = M @ P + # Projective match: dehomogenized Euclidean images agree. + for k in range(P.shape[1]): + wm = mapped[4, k] - mapped[3, k] + wq = Q[4, k] - Q[3, k] + assert abs(wm) > 1e-9 and abs(wq) > 1e-9 + assert np.allclose(mapped[:3, k] / wm, Q[:3, k] / wq, atol=1e-8) + + +def test_conformal_procrustes_full_similarity_cloud(): + """Nontrivial scale + rotation + translation on a (5,K) cloud.""" + rng = np.random.default_rng(239) + X = rng.normal(size=(3, 10)) + s = 1.7 + th = 0.6 + R3 = np.array( + [[np.cos(th), -np.sin(th), 0.0], [np.sin(th), np.cos(th), 0.0], [0.0, 0.0, 1.0]], + dtype=np.float64, + ) + t = np.array([0.5, -0.3, 0.2], dtype=np.float64) + Y = s * (R3 @ X) + t[:, None] + P = np.column_stack([embed_point(X[:, k], dtype=np.float64)[1:6] for k in range(10)]) + Q = np.column_stack([embed_point(Y[:, k], dtype=np.float64)[1:6] for k in range(10)]) + M, residual = conformal_procrustes(P, Q) + assert M.shape == (5, 5) + assert residual < 1e-6 + mapped = M @ P + for k in range(10): + wm = mapped[4, k] - mapped[3, k] + eu = mapped[:3, k] / wm + assert np.allclose(eu, Y[:, k], atol=1e-8) + + +def test_conformal_procrustes_null_point_list_sandwich(): + """List of 32-vec CGA null points recovers V32 with sandwich residual < 1e-6.""" + src_eucl = [ + np.array([0.0, 0.0, 0.0], dtype=np.float64), + np.array([1.0, 0.0, 0.0], dtype=np.float64), + np.array([0.0, 1.0, 0.0], dtype=np.float64), + np.array([0.5, 0.3, 0.1], dtype=np.float64), + np.array([-0.2, 0.4, 0.5], dtype=np.float64), + ] + src = [embed_point(p, dtype=np.float64) for p in src_eucl] + s, th = 1.5, 0.45 + R3 = np.array( + [[np.cos(th), -np.sin(th), 0.0], [np.sin(th), np.cos(th), 0.0], [0.0, 0.0, 1.0]], + dtype=np.float64, + ) + t = np.array([0.25, -0.1, 0.3], dtype=np.float64) + W = geometric_product( + geometric_product(translator(t), dilator(s)), + so3_matrix_to_rotor(R3), + ) + tgt = [versor_apply(W, p) for p in src] + V, residual = conformal_procrustes(src, tgt) + assert V.shape == (32,) + assert versor_condition(V) < 1e-6 + assert residual < 1e-6 + for p, q in zip(src, tgt): + assert float(np.linalg.norm(versor_apply(V, p) - q)) < 1e-6 def test_cartan_iwasawa_extract_closed(): @@ -90,7 +182,11 @@ def test_cartan_iwasawa_extract_closed(): for f in (R, T, D): assert versor_condition(f) < 1e-6 factors = cartan_iwasawa_factorize(V) - assert factors.reconstruction_residual >= 0.0 + recon = geometric_product(geometric_product(factors.R, factors.T), factors.D) + residual = float(np.linalg.norm(recon - V)) + assert residual < 1e-6 + assert factors.reconstruction_residual < 1e-6 + assert abs(factors.reconstruction_residual - residual) < 1e-12 def test_dual_correction_slerp_closed(): @@ -101,6 +197,18 @@ def test_dual_correction_slerp_closed(): assert versor_condition(out) < 1e-6 +def test_dual_correction_slerp_translator_half(): + """Null-bivector power must not erase the Cartan translation leg.""" + from algebra.null_point import recover_translation, translator + + src = _id32() + tgt = translator(np.array([2.0, 0.0, 0.0], dtype=np.float64)) + out = dual_correction_slerp(src, tgt, 0.5) + assert versor_condition(out) < 1e-6 + a, _ = recover_translation(out) + assert np.allclose(a, [1.0, 0.0, 0.0], atol=1e-6) + + def test_pca_rejects_bad_shape(): with pytest.raises(ValueError): signature_aware_pca(np.zeros((4, 3))) diff --git a/tests/test_rotor_power.py b/tests/test_rotor_power.py index 638485aa..680d3c4f 100644 --- a/tests/test_rotor_power.py +++ b/tests/test_rotor_power.py @@ -77,6 +77,20 @@ def test_rotor_power_on_word_transition_preserves_closure() -> None: assert cond < _TOL, f"alpha={alpha}: versor_condition = {cond:.3e}" +def test_rotor_power_null_translator_scales_translation() -> None: + """B²=0 (CGA translator): T^α = 1 + αB, not identity (Cartan dual-slerp).""" + from algebra.null_point import recover_translation, translator + + T = translator(np.array([2.0, 0.0, 0.0], dtype=np.float64)) + half = rotor_power(T, 0.5) + assert versor_condition(half) < _TOL + a, _ = recover_translation(half) + np.testing.assert_allclose(a, [1.0, 0.0, 0.0], atol=1e-9) + # Full power recovers the original translator. + full = rotor_power(T, 1.0) + np.testing.assert_allclose(full, T, atol=1e-9) + + def test_rotor_power_rejects_wrong_shape() -> None: with pytest.raises(ValueError): rotor_power(np.zeros(16, dtype=np.float64), 0.5) diff --git a/tests/test_third_door_blueprint_fidelity.py b/tests/test_third_door_blueprint_fidelity.py index 951bfc96..18ba6255 100644 --- a/tests/test_third_door_blueprint_fidelity.py +++ b/tests/test_third_door_blueprint_fidelity.py @@ -9,20 +9,14 @@ versors (products of rotations on distinct planes — what a real field state looks like) and encodes the properties the Super-Blueprint / R&D-Revised blueprints actually REQUIRE. -The blueprints are the rigorous artifact; the landed code substitutes heuristics. -So the spec-property tests here are marked ``xfail(strict=True)`` with reasons -citing the blueprint section + audit finding. When an operator is implemented to -spec, its xfail flips to xpass (strict) and CI forces the marker's removal. +The blueprints are the rigorous artifact. Spec-property tests here are +behavioral (composed multi-plane inputs, residual < ε, peel-content pins). +Historical findings: + #1 [RESOLVED by #23] supervised_blend no-op on composed versors. + #2 [RESOLVED by #16] Cartan–Iwasawa null-point peel + Spin remainder. + #3 [RESOLVED by #17] Kabsch-conformal Procrustes + field conjugacy. -Empirical findings (2026-07-11 audit, reproduced deterministically below): - #1 [RESOLVED by #23] supervised_blend was a no-op for interior alpha on - composed versors (rotor_power returned identity for the non-simple - transition rotor). Exact fractional powers now implemented. - #2 cartan_iwasawa_factorize raises "factor R not closed" on composed - conformal versors (~45% of the time), instead of the "mathematically - exact, guaranteed" decomposition the Super-Blueprint §2.2 specifies. - -See PR description and the fidelity table for the full spec-vs-impl ledger. +See docs/research/third-door-blueprint-fidelity.md for the living scorecard. """ from __future__ import annotations @@ -30,9 +24,11 @@ from __future__ import annotations import numpy as np import pytest -from algebra.cl41 import geometric_product +from algebra.cl41 import geometric_product, reverse +from algebra.null_point import dilator, translator from algebra.rotor import make_rotor_from_angle -from core.physics.dynamic_manifold import cartan_iwasawa_factorize +from algebra.versor import versor_apply, versor_condition +from core.physics.dynamic_manifold import cartan_iwasawa_factorize, conformal_procrustes from core.physics.goldtether import GoldTetherMonitor @@ -69,30 +65,160 @@ def test_supervised_blend_should_interpolate_composed_versors(): # --- Finding #2: Cartan-Iwasawa decomposition ------------------------------- -def test_cartan_iwasawa_currently_raises_on_composed_versor(): - """CHARACTERIZATION of finding #2 — locks the current (defective) behaviour. +def test_cartan_iwasawa_should_reconstruct_composed_motion(): + """Super §2.2: multi-plane non-similarity factors without raise; residual < 1e-6. - Deterministic composed versor that the heuristic factorizer cannot close. - PASSES today (asserts the raise). Delete when the spec algorithm lands. + Planes (6,7,8) include e1∧e4 (blade 8) — not a pure Euclidean similarity — + so the null-point peel falls through to remainder-as-rotor (R=V, T=I, D=I). """ v = _composed_versor((6, 7, 8), seed=0.0) - with pytest.raises(ValueError, match="not closed"): - cartan_iwasawa_factorize(v) - - -@pytest.mark.xfail( - reason=( - "Finding #2 / Super-Blueprint §2.2: cartan_iwasawa_factorize is specified as a " - "'mathematically exact, non-iterative' decomposition that 'guarantees perfect " - "decomposition' via the action of V on n_o / n_inf. The landed grade-projection " - "heuristic instead raises 'factor R not closed' on composed conformal versors " - "(~45% at 3-4 planes). Spec: factorization must succeed and R*T*D must " - "reconstruct V to < 1e-6." - ), - strict=True, -) -def test_cartan_iwasawa_should_reconstruct_composed_motion(): - v = _composed_versor((6, 7, 8), seed=0.0) - fac = cartan_iwasawa_factorize(v) # spec: must not raise + fac = cartan_iwasawa_factorize(v) # must not raise recon = geometric_product(geometric_product(fac.R, fac.T), fac.D) - assert float(np.linalg.norm(recon - v)) < 1e-6 + residual = float(np.linalg.norm(recon - v)) + assert residual < 1e-6, f"reconstruction residual {residual:.3e}" + assert fac.reconstruction_residual < 1e-6 + for f in (fac.R, fac.T, fac.D): + assert versor_condition(f) < 1e-6 + + +def test_cartan_iwasawa_random_multiplane_never_raises(): + """50 fixed-seed random 3–4 plane products: closed factors, residual < 1e-6.""" + rng = np.random.default_rng(20260713) + max_residual = 0.0 + for _ in range(50): + n = int(rng.integers(3, 5)) # 3 or 4 planes + planes = tuple(int(x) for x in rng.choice(np.arange(6, 16), size=n, replace=False)) + angles = rng.uniform(0.1, 1.0, size=n) + v = _identity() + for ang, p in zip(angles, planes): + v = geometric_product(v, make_rotor_from_angle(float(ang), bivector_idx=int(p))) + fac = cartan_iwasawa_factorize(v) + for f in (fac.R, fac.T, fac.D): + assert versor_condition(f) < 1e-6 + recon = geometric_product(geometric_product(fac.R, fac.T), fac.D) + residual = float(np.linalg.norm(recon - v)) + assert residual < 1e-6 + assert fac.reconstruction_residual < 1e-6 + max_residual = max(max_residual, residual) + assert max_residual < 1e-6 + + +def test_cartan_iwasawa_pure_similarity_peel(): + """V = R*T*D with Euclidean R (planes 6,7,10), nontrivial T and D. + + Pins peel *content* (not residual alone — Spin remainder also has residual 0). + """ + from algebra.null_point import recover_dilation, recover_translation + + R_e = geometric_product( + make_rotor_from_angle(0.4, bivector_idx=6), + make_rotor_from_angle(0.3, bivector_idx=7), + ) + R_e = geometric_product(R_e, make_rotor_from_angle(0.25, bivector_idx=10)) + t_vec = np.array([0.5, -0.2, 0.1], dtype=np.float64) + T = translator(t_vec) + D = dilator(1.7) + V = geometric_product(geometric_product(R_e, T), D) + fac = cartan_iwasawa_factorize(V) + recon = geometric_product(geometric_product(fac.R, fac.T), fac.D) + residual = float(np.linalg.norm(recon - V)) + assert residual < 1e-6, f"similarity peel residual {residual:.3e}" + assert fac.reconstruction_residual < 1e-6 + for f in (fac.R, fac.T, fac.D): + assert versor_condition(f) < 1e-6 + I = _identity() + # Must have taken the peel path — not silent Spin remainder. + assert float(np.linalg.norm(fac.D - I)) > 1e-3 + assert float(np.linalg.norm(fac.T - I)) > 1e-3 + s_rec, _ = recover_dilation(fac.D) + assert abs(s_rec - 1.7) < 1e-6 + # Translation content: a is the origin image under R·T (R conjugates the + # Euclidean displacement). Assert nontrivial finite translation, not a==t. + a_rec, _ = recover_translation(fac.T) + assert float(np.linalg.norm(a_rec)) > 1e-3 + assert np.isfinite(a_rec).all() + # R·T must recover the de-dilated motion (peel identity). + RT = geometric_product(fac.R, fac.T) + V1 = geometric_product(V, reverse(fac.D)) + assert float(np.linalg.norm(RT - V1)) < 1e-6 + + +def test_cartan_iwasawa_pure_dilator_round_trip(): + from algebra.null_point import recover_dilation + + V = dilator(2.5) + fac = cartan_iwasawa_factorize(V) + recon = geometric_product(geometric_product(fac.R, fac.T), fac.D) + assert float(np.linalg.norm(recon - V)) < 1e-6 + assert fac.reconstruction_residual < 1e-6 + for f in (fac.R, fac.T, fac.D): + assert versor_condition(f) < 1e-6 + I = _identity() + assert float(np.linalg.norm(fac.R - I)) < 1e-9 + assert float(np.linalg.norm(fac.T - I)) < 1e-9 + s_rec, _ = recover_dilation(fac.D) + assert abs(s_rec - 2.5) < 1e-6 + + +def test_cartan_iwasawa_pure_translator_round_trip(): + from algebra.null_point import recover_translation + + t_vec = np.array([1.0, -0.5, 0.25], dtype=np.float64) + V = translator(t_vec) + fac = cartan_iwasawa_factorize(V) + recon = geometric_product(geometric_product(fac.R, fac.T), fac.D) + assert float(np.linalg.norm(recon - V)) < 1e-6 + assert fac.reconstruction_residual < 1e-6 + for f in (fac.R, fac.T, fac.D): + assert versor_condition(f) < 1e-6 + I = _identity() + assert float(np.linalg.norm(fac.R - I)) < 1e-9 + assert float(np.linalg.norm(fac.D - I)) < 1e-9 + a_rec, _ = recover_translation(fac.T) + assert np.allclose(a_rec, t_vec, atol=1e-6) + + +def test_cartan_iwasawa_identity_factors_cleanly(): + V = _identity() + fac = cartan_iwasawa_factorize(V) + assert float(np.linalg.norm(fac.R - V)) < 1e-12 or fac.reconstruction_residual < 1e-12 + recon = geometric_product(geometric_product(fac.R, fac.T), fac.D) + assert float(np.linalg.norm(recon - V)) < 1e-12 + assert fac.reconstruction_residual < 1e-12 + for f in (fac.R, fac.T, fac.D): + assert versor_condition(f) < 1e-6 + + +def test_cartan_iwasawa_rejects_non_versor(): + bad = np.ones(32, dtype=np.float64) + with pytest.raises(ValueError): + cartan_iwasawa_factorize(bad) + + +# --- Finding #3 / gap #17: Kabsch-conformal Procrustes field conjugacy -------- + +def test_conformal_procrustes_field_conjugacy_nontrivial(): + """#17: non-trivial multiplane F_A, F_B = sandwich(W, F_A); sandwich residual < 1e-5. + + Behavioral pin: residual is measured under ``versor_apply`` (sandwich), not + left-composition via ``word_transition_rotor``. Identity→identity is excluded. + """ + F_A = _composed_versor((6, 7, 10, 11), seed=1.5) + W = geometric_product( + geometric_product( + make_rotor_from_angle(0.55, bivector_idx=6), + make_rotor_from_angle(0.4, bivector_idx=7), + ), + make_rotor_from_angle(0.3, bivector_idx=10), + ) + F_B = versor_apply(W, F_A) + assert float(np.linalg.norm(F_A - _identity())) > 1e-3 + assert float(np.linalg.norm(F_B - F_A)) > 1e-3 + + V, residual = conformal_procrustes(F_A, F_B) + assert V.shape == (32,) + assert versor_condition(V) < 1e-6 + assert residual < 1e-5, f"sandwich residual {residual:.3e}" + sand = float(np.linalg.norm(versor_apply(V, F_A) - F_B)) + assert sand < 1e-5, f"versor_apply residual {sand:.3e}" +