Merge pull request #562 from AssetOverflow/feat/cga-incidence-primitives

algebra: incidence algebra — graded_wedge, dual, meet + honest outer_product

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

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 = <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(|<X reverse(X)>_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)))

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