init: algebra layer — cl41, versor, cga, holonomy

algebra/__init__.py

@@ -0,0 +1,4 @@
+from .cl41 import geometric_product, reverse, grade_project, scalar_part, norm_squared, basis_vector
+from .versor import versor_apply, normalize_to_versor, versor_condition
+from .cga import cga_inner, outer_product, is_null, null_project, embed_point
+from .holonomy import holonomy_encode, holonomy_similarity

algebra/cga.py

@@ -0,0 +1,70 @@
+"""
+Conformal Geometric Algebra geometry on Cl(4,1).
+
+Key identity: for null vectors X, Y on the horosphere,
+    cga_inner(X, Y) = -d(X, Y)^2 / 2
+where d is Euclidean distance.
+
+This is the ONLY distance metric in CORE-AI.
+No cosine similarity. No L2 norm. No approximate indexing.
+"""
+
+import numpy as np
+from .cl41 import geometric_product, reverse, scalar_part, basis_vector, N_COMPONENTS
+
+
+def cga_inner(X: np.ndarray, Y: np.ndarray) -> float:
+    """
+    Symmetric inner product: 0.5 * scalar_part(X*Y + Y*X).
+    For null vectors representing conformal points: equals -d^2 / 2.
+    """
+    XY = geometric_product(X, Y)
+    YX = geometric_product(Y, X)
+    return 0.5 * scalar_part(XY + YX)
+
+
+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.
+    A real (non-imaginary) result means the dialogue is coherent.
+    """
+    XY = geometric_product(X, Y)
+    YX = geometric_product(Y, X)
+    return 0.5 * (XY - YX)
+
+
+def is_null(X: np.ndarray, tol: float = 1e-6) -> bool:
+    """Check if X lies on the null cone: X*X = 0."""
+    return abs(cga_inner(X, X)) < tol
+
+
+def null_project(X: np.ndarray) -> np.ndarray:
+    """
+    Re-project X onto the null cone.
+    Call on vault entries periodically to correct floating-point null-cone drift.
+    This is numerical maintenance, not a heat shield.
+    Method: extract Euclidean part, re-embed via standard CGA point map.
+    """
+    euclidean = X[1:4].copy().astype(np.float32)
+    x_sq = float(np.dot(euclidean, euclidean))
+    result = np.zeros(N_COMPONENTS, dtype=np.float32)
+    result[1:4] = euclidean
+    result[4] = 0.5 * x_sq   # e+ coefficient
+    result[5] = 1.0           # e- coefficient
+    return result
+
+
+def embed_point(x: np.ndarray) -> np.ndarray:
+    """
+    Embed a Euclidean point x in R^3 into the CGA null cone.
+    Standard map: X = x + (1/2)|x|^2 * e+ + e-
+    """
+    x = np.asarray(x, dtype=np.float32)
+    assert len(x) == 3, "embed_point expects a 3D vector"
+    result = np.zeros(N_COMPONENTS, dtype=np.float32)
+    result[1:4] = x
+    result[4] = 0.5 * float(np.dot(x, x))
+    result[5] = 1.0
+    return result

algebra/cl41.py

@@ -0,0 +1,174 @@
+"""
+Cl(4,1) multivector arithmetic.
+
+Signature: (+,+,+,+,-). Basis e1..e5.
+Multivectors are float32 arrays of shape (32,) ordered by grade:
+  grade-0: index 0        (1 component)
+  grade-1: indices 1-5    (5 components)
+  grade-2: indices 6-15   (10 components)
+  grade-3: indices 16-25  (10 components)
+  grade-4: indices 26-30  (5 components)
+  grade-5: index 31       (1 component)
+"""
+
+from itertools import combinations
+from math import comb
+import numpy as np
+
+N_DIMS = 5
+N_COMPONENTS = 32
+SIGNATURE = np.array([1, 1, 1, 1, -1], dtype=np.float64)
+
+# --- Grade offset table ---
+
+def _grade_offsets():
+    offsets = []
+    start = 0
+    for k in range(N_DIMS + 1):
+        count = comb(N_DIMS, k)
+        offsets.append((start, count))
+        start += count
+    return offsets
+
+_GRADE_OFFSETS = _grade_offsets()
+
+def grade_start(k: int) -> int:
+    return _GRADE_OFFSETS[k][0]
+
+def grade_count(k: int) -> int:
+    return _GRADE_OFFSETS[k][1]
+
+# --- Blade index maps ---
+
+def _all_blades():
+    """Return ordered list of blade tuples (one per component, ordered by grade)."""
+    blades = []
+    for k in range(N_DIMS + 1):
+        for combo in combinations(range(N_DIMS), k):
+            blades.append(combo)
+    return blades
+
+_BLADES = _all_blades()  # index -> tuple of basis vector indices
+_BLADE_TO_IDX = {b: i for i, b in enumerate(_BLADES)}
+
+def _blade_product(blade_a, blade_b):
+    """Compute geometric product of two basis blades. Returns (sign, result_blade_tuple)."""
+    # Concatenate and bubble-sort, tracking sign flips and metric contractions
+    seq = list(blade_a) + list(blade_b)
+    sign = 1
+    # Bubble sort to canonical order, tracking swaps
+    n = len(seq)
+    for i in range(n):
+        for j in range(n - i - 1):
+            if seq[j] > seq[j + 1]:
+                seq[j], seq[j + 1] = seq[j + 1], seq[j]
+                sign *= -1
+            elif seq[j] == seq[j + 1]:
+                # Metric contraction: e_i^2 = signature[i]
+                metric = int(SIGNATURE[seq[j]])
+                sign *= metric
+                seq.pop(j)
+                seq.pop(j)  # second element now at same index after first pop
+                n -= 2
+                break
+        else:
+            continue
+        break
+    # After contraction there may still be duplicates — recurse
+    result = tuple(sorted(set(seq)))  # this is wrong for multi-contraction; use proper loop
+    return sign, tuple(seq)
+
+def _build_multiplication_table():
+    """Precompute full 32x32 geometric product table."""
+    table_idx = np.zeros((N_COMPONENTS, N_COMPONENTS), dtype=np.int32)
+    table_sign = np.zeros((N_COMPONENTS, N_COMPONENTS), dtype=np.float32)
+
+    for i, blade_a in enumerate(_BLADES):
+        for j, blade_b in enumerate(_BLADES):
+            sign, result_blade = _compute_blade_product(blade_a, blade_b)
+            result_idx = _BLADE_TO_IDX.get(result_blade, 0)
+            table_idx[i, j] = result_idx
+            table_sign[i, j] = sign
+
+    return table_idx, table_sign
+
+def _compute_blade_product(blade_a, blade_b):
+    """Compute geometric product of two basis blades via bubble sort + metric."""
+    seq = list(blade_a) + list(blade_b)
+    sign = 1
+    i = 0
+    while i < len(seq) - 1:
+        j = i
+        while j < len(seq) - 1:
+            if seq[j] == seq[j + 1]:
+                # Contract: e_k^2 = signature[k]
+                sign *= int(SIGNATURE[seq[j]])
+                seq.pop(j)
+                seq.pop(j)
+                if j > 0:
+                    i = max(0, j - 1)
+                break
+            elif seq[j] > seq[j + 1]:
+                seq[j], seq[j + 1] = seq[j + 1], seq[j]
+                sign *= -1
+                j += 1
+            else:
+                j += 1
+        else:
+            i += 1
+    result_blade = tuple(seq)
+    if result_blade not in _BLADE_TO_IDX:
+        return 0, ()
+    return sign, result_blade
+
+_TABLE_IDX, _TABLE_SIGN = _build_multiplication_table()
+
+# --- Core operations ---
+
+def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
+    """Full geometric product in Cl(4,1)."""
+    A = np.asarray(A, dtype=np.float32)
+    B = np.asarray(B, dtype=np.float32)
+    result = np.zeros(N_COMPONENTS, dtype=np.float32)
+    for i in range(N_COMPONENTS):
+        if A[i] == 0.0:
+            continue
+        for j in range(N_COMPONENTS):
+            if B[j] == 0.0:
+                continue
+            result[_TABLE_IDX[i, j]] += _TABLE_SIGN[i, j] * A[i] * B[j]
+    return result
+
+def reverse(A: np.ndarray) -> np.ndarray:
+    """
+    Reverse (main anti-automorphism).
+    Grade-k blades pick up sign (-1)^(k*(k-1)/2).
+    Grade 0,1: +1.  Grade 2,3: -1.  Grade 4,5: +1.
+    """
+    A = np.asarray(A, dtype=np.float32).copy()
+    # Grade 2: indices 6-15
+    A[6:16] *= -1.0
+    # Grade 3: indices 16-25
+    A[16:26] *= -1.0
+    return A
+
+def grade_project(A: np.ndarray, k: int) -> np.ndarray:
+    """Extract grade-k part of A."""
+    result = np.zeros(N_COMPONENTS, dtype=np.float32)
+    start, count = _GRADE_OFFSETS[k]
+    result[start:start + count] = A[start:start + count]
+    return result
+
+def scalar_part(A: np.ndarray) -> float:
+    """Return grade-0 component."""
+    return float(A[0])
+
+def norm_squared(A: np.ndarray) -> float:
+    """||A||^2 = scalar_part(A * reverse(A))."""
+    return scalar_part(geometric_product(A, reverse(A)))
+
+def basis_vector(i: int) -> np.ndarray:
+    """Return the i-th basis vector (0-indexed) as a 32-component multivector."""
+    v = np.zeros(N_COMPONENTS, dtype=np.float32)
+    v[1 + i] = 1.0
+    return v

algebra/holonomy.py

@@ -0,0 +1,71 @@
+"""
+Holonomy prompt encoding.
+
+A prompt w1, w2, ..., wn is encoded as the geometric holonomy of its
+forward+reverse versor walk. The walk closes, producing a versor that
+is bounded by construction and invariant to global phase.
+
+The holonomy IS a versor — it drops directly into versor_apply with
+no bridging code. The fuel and the engine are the same substance.
+"""
+
+import numpy as np
+from .cl41 import geometric_product, reverse as cl_reverse
+from .versor import normalize_to_versor
+from .cga import cga_inner
+
+
+def holonomy_encode(
+    word_versors: list,
+    alpha: float = 0.5,
+    weights: list = None,
+) -> np.ndarray:
+    """
+    Compute the holonomy of the word versor sequence.
+
+    Forward walk:  F = w1 * w2 * ... * wn  (weighted by word frequency inverse)
+    Reverse walk:  R = (1-alpha) * reverse(wn) * ... * reverse(w1)
+    Holonomy:      H = geometric_product(F, R)
+
+    H is a versor. For alpha=0.5, the holonomy captures the geometric
+    curvature of the prompt path. Prompts with different semantic content
+    produce geometrically distinct holonomies even at the same length.
+
+    weights: optional list of float scalars (e.g. inverse token frequency).
+             Rare content words rotate more than common function words.
+             If None, uniform weights are used.
+    """
+    if not word_versors:
+        raise ValueError("Cannot encode empty prompt.")
+
+    n = len(word_versors)
+    if weights is None:
+        weights = [1.0] * n
+    assert len(weights) == n
+
+    # Forward accumulation
+    F = word_versors[0].copy() * weights[0]
+    F = normalize_to_versor(F)
+    for k in range(1, n):
+        w = word_versors[k] * weights[k]
+        w = normalize_to_versor(w)
+        F = geometric_product(F, w)
+
+    # Reverse accumulation with alpha damping
+    R = cl_reverse(word_versors[-1]) * (1.0 - alpha)
+    R = normalize_to_versor(R)
+    for k in range(n - 2, -1, -1):
+        r = cl_reverse(word_versors[k])
+        r = normalize_to_versor(r)
+        R = geometric_product(r, R)
+
+    H = geometric_product(F, R)
+    return normalize_to_versor(H)
+
+
+def holonomy_similarity(H1: np.ndarray, H2: np.ndarray) -> float:
+    """
+    Compare two holonomies via CGA inner product.
+    Used for prompt-level semantic similarity without embedding lookup.
+    """
+    return cga_inner(H1, H2)

algebra/versor.py

@@ -0,0 +1,50 @@
+"""
+The three versor primitives.
+
+These are the ONLY normalization/transition/check functions in the system.
+Do not add correction, monitoring, or grade-guard functions here.
+If you think you need something else, you have an unclosed operation upstream.
+"""
+
+import numpy as np
+from .cl41 import geometric_product, reverse, scalar_part, norm_squared
+
+
+def versor_apply(V: np.ndarray, F: np.ndarray) -> np.ndarray:
+    """
+    Sandwich product: V * F * reverse(V).
+
+    The ONLY allowed field transition in the system.
+    Algebraically closed on the versor manifold:
+    if V and F are versors, V*F*reverse(V) is a versor.
+    No pre/post normalization. No grade projection. No guards.
+    """
+    return geometric_product(V, geometric_product(F, reverse(V)))
+
+
+def normalize_to_versor(F: np.ndarray) -> np.ndarray:
+    """
+    Project F onto the versor manifold: F / sqrt(|F * reverse(F)|).
+
+    Call this ONCE per input at the injection gate (ingest/gate.py).
+    Never call mid-propagation, mid-generation, or in the vault.
+    If you feel the urge to call this elsewhere, fix the upstream operation.
+    """
+    n2 = norm_squared(F)
+    if abs(n2) < 1e-12:
+        raise ValueError("Cannot normalize a null multivector to a versor.")
+    return F / np.sqrt(abs(n2))
+
+
+def versor_condition(F: np.ndarray) -> float:
+    """
+    Returns ||F * reverse(F) - 1||_F.
+
+    Zero means F is on the versor manifold.
+    Use in tests and at the injection gate only.
+    Never call in the generation hot path.
+    """
+    product = geometric_product(F, reverse(F))
+    product = product.copy()
+    product[0] -= 1.0
+    return float(np.linalg.norm(product))