From de645055ea7a3f53cd4c75a38366d409db9903b0 Mon Sep 17 00:00:00 2001 From: Shay Date: Thu, 4 Jun 2026 21:43:35 -0700 Subject: [PATCH] =?UTF-8?q?feat(algebra):=20incidence=20algebra=20?= =?UTF-8?q?=E2=80=94=20graded=5Fwedge,=20dual,=20meet=20+=20honest=20outer?= =?UTF-8?q?=5Fproduct?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the correct grade-raising "wire" the field substrate was missing — so cga_inner can operate on RELATIONS among entities (lines/planes/incidence), not just pairwise point distance. Built only from existing Cl(4,1) primitives (geometric_product, grade_project) + the pseudoscalar; no normalization, no approximation, versor_condition path untouched (flats are null-cone wedges, not unit versors). - outer_product: DOCSTRING-ONLY honesty fix (behavior byte-identical, every caller unchanged). It is the commutator 0.5*(XY-YX) = the wedge ONLY for grade-1 vectors; for higher grades it is the Lie bracket, NOT the wedge, and does NOT build a k-blade by repetition. Existing callers consume it as an opaque cga_inner-reduced feature (none read it by grade), so the relabel is safe. Points to graded_wedge for the real exterior product. - graded_wedge(X,Y) = _{grade(X)+grade(Y)} — the true wedge; agrees with outer_product on grade-1, differs above (pinned by test). Builds lines/planes. - is_incident(point, flat): EXACT zero-test (point^flat == 0, no float tolerance to admit — near-incident is refused, per wrong=0). Exact at scale in f64. - dual(X) = X*I5^{-1} (I5^2=-1 confirmed); involutive up to sign. - meet(A,B) = dual(dual(A)^dual(B)): correct for spanning operands (two planes -> their line, incidence verified). HONEST ENVELOPE: degenerates for non-spanning operands (coplanar lines) — returns the ZERO multivector (detectable, documented, tested), never a silent wrong value. The general coplanar intersection needs the join-relative meet, deliberately NOT faked here. Green: smoke 87, algebra 82, incidence 8, outer_product consumers + invariants 109; zero regressions (outer_product behavior unchanged). --- algebra/__init__.py | 6 ++ algebra/cga.py | 114 ++++++++++++++++++++++++++++++-- tests/test_cga_incidence.py | 128 ++++++++++++++++++++++++++++++++++++ 3 files changed, 243 insertions(+), 5 deletions(-) create mode 100644 tests/test_cga_incidence.py diff --git a/algebra/__init__.py b/algebra/__init__.py index 6d6c3abf..bd06f292 100644 --- a/algebra/__init__.py +++ b/algebra/__init__.py @@ -8,6 +8,12 @@ from .cga import ( null_project, embed_point, read_scalar_e1, + blade_grade, + blade_norm, + graded_wedge, + is_incident, + dual, + meet, ) from .holonomy import holonomy_encode, holonomy_similarity from .rotor import word_transition_rotor diff --git a/algebra/cga.py b/algebra/cga.py index 42ff6b23..a4cdf9a8 100644 --- a/algebra/cga.py +++ b/algebra/cga.py @@ -19,7 +19,20 @@ No cosine similarity. No L2 norm. No approximate indexing. """ import numpy as np -from .cl41 import geometric_product, scalar_part, basis_vector, N_COMPONENTS +from .cl41 import ( + geometric_product, + grade_project, + reverse, + scalar_part, + N_COMPONENTS, +) + +# The unit pseudoscalar I5 = e1 e2 e3 e4 e5 (the grade-5 blade, component 31). +# In Cl(4,1) with signature (+,+,+,+,-), I5^2 = -1, so I5^{-1} = -I5. Used by +# ``dual`` / ``meet``. Module-level singleton; never mutated. +_PSEUDOSCALAR_INDEX = 31 +_I5 = np.zeros(N_COMPONENTS, dtype=np.float64) +_I5[_PSEUDOSCALAR_INDEX] = 1.0 # Basis-vector component indices for e4/e5 inside the grade-1 block. # component 1=e1, 2=e2, 3=e3, 4=e4, 5=e5. @@ -46,10 +59,19 @@ def cga_inner(X: np.ndarray, Y: np.ndarray) -> float: def outer_product(X: np.ndarray, Y: np.ndarray) -> np.ndarray: - """ - Outer (wedge) product: X ^ Y. - For a prompt versor X_p and response versor X_r, - X_p ^ X_r is a grade-2 object encoding their geometric relationship. + """The antisymmetric (commutator) product ``0.5 * (XY - YX)``. + + HONEST CONTRACT: this equals the grade-raising wedge ``X ^ Y`` **only when both + operands are grade 1** (vectors). For higher-grade operands it is the *commutator* + (Lie bracket), which is NOT the wedge — in particular it does NOT build a k-blade + by repeated application (a bivector commuted with a vector collapses the grade-3 + part to grade 1). Existing callers use the result as an opaque, deterministic + relationship feature (folded into a scalar via :func:`cga_inner`), where the + commutator is well-defined regardless; none read it by grade. + + For the true grade-raising exterior product (lines/planes/incidence) use + :func:`graded_wedge`. (Renamed contract only — behaviour is unchanged, so every + current caller is byte-identical.) """ XY = geometric_product(X, Y) YX = geometric_product(Y, X) @@ -120,3 +142,85 @@ def read_scalar_e1(X: np.ndarray) -> float: "read_scalar_e1: degenerate n_o weight (point at infinity or f32 collapse)" ) return float(X[1]) / no_weight + + +# --------------------------------------------------------------------------- +# Incidence algebra — the corrected grade-raising wedge, dual, and meet. +# These let the inner product operate on RELATIONS among entities (lines, planes, +# incidence) rather than only pairwise point distance. Built only from the existing +# Cl(4,1) primitives (geometric_product, grade_project) + the pseudoscalar; they add +# no normalization, no approximation, and leave the versor_condition path untouched +# (flats are null-cone outer products, not unit versors). +# --------------------------------------------------------------------------- + +_MAX_GRADE = 5 # Cl(4,1): grades 0..5 + + +def blade_grade(X: np.ndarray) -> int: + """The single grade of a homogeneous blade. Raises if X is zero or grade-mixed. + + Grade is detected by EXACT nonzero (no tolerance): integer-coordinate embeddings + produce exact integer blades in float64, so a grade block is exactly 0 or not. + """ + grades = [k for k in range(_MAX_GRADE + 1) if np.any(grade_project(X, k))] + if len(grades) != 1: + raise ValueError(f"not a homogeneous blade: nonzero grades {grades}") + return grades[0] + + +def graded_wedge(X: np.ndarray, Y: np.ndarray) -> np.ndarray: + """The true grade-raising exterior product ``X ^ Y`` for homogeneous blades. + + ``X ^ Y = _{grade(X)+grade(Y)}`` — the top-grade part of the geometric + product. Unlike :func:`outer_product` (the commutator) this composes correctly: + ``graded_wedge(graded_wedge(P, Q), n_inf)`` builds the grade-3 line P^Q^n_inf, + and so on. If the grades sum past the pseudoscalar (>5) the wedge is identically + zero. For two grade-1 vectors it agrees with :func:`outer_product` exactly. + """ + gx, gy = blade_grade(X), blade_grade(Y) + if gx + gy > _MAX_GRADE: + return np.zeros(N_COMPONENTS, dtype=geometric_product(X, Y).dtype) + return grade_project(geometric_product(X, Y), gx + gy) + + +def blade_norm(X: np.ndarray) -> float: + """Reversion norm ``sqrt(|_0|)`` — zero iff X is the zero blade.""" + return float(np.sqrt(abs(scalar_part(geometric_product(X, reverse(X)))))) + + +def is_incident(point: np.ndarray, flat: np.ndarray) -> bool: + """Exact incidence test: is ``point`` on ``flat`` (a line/plane OPNS blade)? + + True iff ``point ^ flat == 0`` EXACTLY (every component zero) — no float + tolerance to admit (the wrong=0 discipline: a near-incident point is REFUSED, + not admitted). Exact for integer-coordinate points within ``EMBED_EXACT_MAX``. + """ + return not bool(np.any(graded_wedge(point, flat))) + + +def dual(X: np.ndarray) -> np.ndarray: + """Pseudoscalar dual ``X * I5^{-1}`` (``I5^{-1} = -I5`` since ``I5^2 = -1``). + + Maps a grade-k blade to grade ``5-k``. Involutive up to sign: + ``dual(dual(X)) == -X``. + """ + return geometric_product(X, -_I5) + + +def meet(A: np.ndarray, B: np.ndarray) -> np.ndarray: + """The meet (intersection) ``dual(dual(A) ^ dual(B))`` of two homogeneous blades. + + Correct for operands in GENERAL POSITION whose join spans the space — e.g. two + non-parallel planes meet in their intersection line. The grade of the result is + ``grade(A)+grade(B)-5``. + + HONEST ENVELOPE: this full-pseudoscalar meet DEGENERATES for operands that share + a proper subspace (e.g. two coplanar lines, two parallel planes): the inner wedge + ``dual(A) ^ dual(B)`` is then identically zero, so ``meet`` returns the **zero + multivector** — a detectable signal of "no transversal meet", never a silently + wrong value. The general intersection of such operands (e.g. the point where two + coplanar lines cross) requires the *join-relative* meet, which is deliberately + NOT implemented here; the caller MUST check ``blade_norm(result) == 0`` and treat + zero as degenerate/refuse rather than as a geometric object. + """ + return dual(graded_wedge(dual(A), dual(B))) diff --git a/tests/test_cga_incidence.py b/tests/test_cga_incidence.py new file mode 100644 index 00000000..2fbf95d8 --- /dev/null +++ b/tests/test_cga_incidence.py @@ -0,0 +1,128 @@ +"""CGA incidence algebra — the corrected grade-raising wedge, dual, and meet. + +These primitives let `cga_inner` operate on RELATIONS among entities (lines, planes, +incidence) rather than only pairwise point distance — the missing "wire" for the +section-level relational layer. Every test pins EXACT behaviour (integer-coordinate +embeddings → exact integer blades in float64), no float tolerance to admit. + +It also pins the honest distinction the `outer_product` docstring now states: the +corrected `graded_wedge` agrees with `outer_product` for grade-1 vectors but DIFFERS +for higher grades (where `outer_product` is the commutator, not the wedge) — and the +honest envelope of `meet` (degenerate operands return zero, never a silent wrong). +""" + +from __future__ import annotations + +import numpy as np + +from algebra.cga import ( + EMBED_EXACT_MAX, + blade_grade, + blade_norm, + dual, + embed_point, + graded_wedge, + is_incident, + meet, + outer_product, +) +from algebra.cl41 import N_COMPONENTS, geometric_product, scalar_part + +_F64 = np.float64 + + +def _pt(x: float, y: float = 0.0, z: float = 0.0) -> np.ndarray: + return embed_point(np.array([x, y, z], dtype=_F64), dtype=_F64) + + +def _n_inf() -> np.ndarray: + v = np.zeros(N_COMPONENTS, dtype=_F64) + v[4] = 1.0 + v[5] = 1.0 + return v + + +def _line(p: np.ndarray, q: np.ndarray) -> np.ndarray: + """OPNS line p ^ q ^ n_inf (grade 3).""" + return graded_wedge(graded_wedge(p, q), _n_inf()) + + +def _plane(p: np.ndarray, q: np.ndarray, r: np.ndarray) -> np.ndarray: + """OPNS plane p ^ q ^ r ^ n_inf (grade 4).""" + return graded_wedge(graded_wedge(graded_wedge(p, q), r), _n_inf()) + + +# --- graded_wedge agrees with outer_product on grade-1, differs above ------ + + +def test_graded_wedge_agrees_with_outer_product_for_vectors(): + a = np.zeros(N_COMPONENTS, dtype=_F64); a[1] = 1.0 # e1 + b = np.zeros(N_COMPONENTS, dtype=_F64); b[2] = 1.0 # e2 + np.testing.assert_array_equal(graded_wedge(a, b), outer_product(a, b)) + assert blade_grade(graded_wedge(a, b)) == 2 + + +def test_graded_wedge_differs_from_commutator_above_grade_1(): + """The honest distinction outer_product's docstring now states: building a + 3-blade by repeated wedge works for graded_wedge but COLLAPSES for the commutator.""" + p, q, ninf = _pt(0, 0), _pt(2, 0), _n_inf() + pq = graded_wedge(p, q) # grade 2 + line_true = graded_wedge(pq, ninf) # grade 3 — the real line + line_commutator = outer_product(pq, ninf) + assert blade_grade(line_true) == 3 + # the commutator does NOT yield the grade-3 line (it collapses the top grade) + assert not np.array_equal(line_true, line_commutator) + + +# --- incidence: exact, and exact at scale (f64) ---------------------------- + + +def test_incidence_collinear_exact(): + line = _line(_pt(0, 0), _pt(2, 0)) + assert is_incident(_pt(5, 0), line) # collinear beyond the segment + assert is_incident(_pt(1, 0), line) # on the segment + assert not is_incident(_pt(0, 1), line) # off the line + assert not is_incident(_pt(3, 2), line) + + +def test_incidence_exact_at_scale(): + """f64 keeps incidence exact for large integer coordinates (within the ceiling).""" + v = EMBED_EXACT_MAX // 2 + line = _line(_pt(0, 0), _pt(2, 0)) # the x-axis + assert is_incident(_pt(float(v), 0), line) + assert not is_incident(_pt(float(v), 1), line) + + +# --- dual ------------------------------------------------------------------ + + +def test_pseudoscalar_squares_to_minus_one(): + i5 = np.zeros(N_COMPONENTS, dtype=_F64); i5[31] = 1.0 + assert scalar_part(geometric_product(i5, i5)) == -1.0 + + +def test_dual_is_involution_up_to_sign(): + x = _pt(3, 0) + np.testing.assert_allclose(dual(dual(x)), -x, atol=0.0) + + +# --- meet: correct for spanning operands, honest-zero for degenerate ------- + + +def test_meet_of_two_planes_is_their_line(): + p_z0 = _plane(_pt(0, 0, 0), _pt(1, 0, 0), _pt(0, 1, 0)) # z = 0 + p_y0 = _plane(_pt(0, 0, 0), _pt(1, 0, 0), _pt(0, 0, 1)) # y = 0 + line = meet(p_z0, p_y0) # expect the x-axis + assert blade_norm(line) != 0.0 + assert blade_grade(line) == 3 + assert is_incident(_pt(5, 0, 0), line) # x-axis point is on it + assert not is_incident(_pt(0, 5, 0), line) # off it + + +def test_meet_degenerate_operands_return_zero_not_silent_wrong(): + """The honest envelope: coplanar lines do not span, so the full-pseudoscalar meet + DEGENERATES — it returns the zero multivector (detectable), never a wrong object.""" + l1 = _line(_pt(0, 0), _pt(2, 0)) # x-axis (z=0 plane) + l2 = _line(_pt(1, -1), _pt(1, 1)) # x=1 vertical (z=0 plane) — coplanar with l1 + result = meet(l1, l2) + assert blade_norm(result) == 0.0 # degenerate → zero, caller must refuse