Core

The mathematical kernel of Versor.

Algebra

CliffordAlgebra

Differentiable Clifford algebra kernel with memory-optimized blocked accumulation.

Handles geometric product, grade projection, and rotor operations.

Supports degenerate (null) dimensions via the r parameter: Cl(p, q, r) has p positive, q negative, and r null basis vectors (e_i^2 = 0).

Attributes:

Name	Type	Description
p	int	Positive signature dimensions.
q	int	Negative signature dimensions.
r	int	Degenerate (null) dimensions.
n	int	Total dimensions (p + q + r).
dim	int	Total basis elements (2^n).
device	str	Computation device.

Source code in core/algebra.py

class CliffordAlgebra:
    """Differentiable Clifford algebra kernel with memory-optimized blocked accumulation.

    Handles geometric product, grade projection, and rotor operations.

    Supports degenerate (null) dimensions via the ``r`` parameter:
    ``Cl(p, q, r)`` has ``p`` positive, ``q`` negative, and ``r`` null
    basis vectors (``e_i^2 = 0``).

    Attributes:
        p (int): Positive signature dimensions.
        q (int): Negative signature dimensions.
        r (int): Degenerate (null) dimensions.
        n (int): Total dimensions (p + q + r).
        dim (int): Total basis elements (2^n).
        device (str): Computation device.
    """
    _CACHED_TABLES = {}

    def __init__(self, p: int, q: int = 0, r: int = 0, device='cuda'):
        """Initialize the algebra and cache the Cayley table.

        Args:
            p (int): Positive dimensions (+1).
            q (int, optional): Negative dimensions (-1). Defaults to 0.
            r (int, optional): Degenerate dimensions (0). Defaults to 0.
            device (str, optional): The device on which computations are performed. Defaults to 'cuda'.
        """
        assert p >= 0, f"p must be non-negative, got {p}"
        assert q >= 0, f"q must be non-negative, got {q}"
        assert r >= 0, f"r must be non-negative, got {r}"
        assert p + q + r <= 12, f"p + q + r must be <= 12, got {p + q + r}"

        self.p, self.q, self.r = p, q, r
        self.n = p + q + r
        self.dim = 2 ** self.n
        self.device = device

        # Cache Cayley tables to avoid recomputation
        cache_key = (p, q, r, str(device))
        if cache_key not in CliffordAlgebra._CACHED_TABLES:
            CliffordAlgebra._CACHED_TABLES[cache_key] = self._generate_cayley_table()

        (
            self.cayley_indices,
            self.cayley_signs,
            self.gp_signs,
            self.grade_masks,
            self.rev_signs,
            self.bv_sq_scalar,
            self.wedge_gp_signs,
            self.inner_gp_signs,
            self.grade_index,
            self.rc_action,
        ) = CliffordAlgebra._CACHED_TABLES[cache_key]

    @property
    def num_grades(self) -> int:
        """Counts the number of grades (n + 1).

        Returns:
            int: Number of grades.
        """
        return self.n + 1

    def embed_vector(self, vectors: torch.Tensor) -> torch.Tensor:
        """Injects vectors into the Grade-1 subspace.

        Args:
            vectors (torch.Tensor): Raw vectors [..., n].

        Returns:
            torch.Tensor: Multivector coefficients [..., dim].
        """
        batch_shape = vectors.shape[:-1]
        mv = torch.zeros(*batch_shape, self.dim, device=vectors.device, dtype=vectors.dtype)
        for i in range(self.n):
            mv[..., 1 << i] = vectors[..., i]
        return mv

    def get_grade_norms(self, mv: torch.Tensor) -> torch.Tensor:
        """Calculates norms per grade. Useful for invariant features.

        Vectorized via scatter_add -- no Python loops over grades.

        Args:
            mv (torch.Tensor): Input multivector [..., dim].

        Returns:
            torch.Tensor: Grade norms [..., num_grades].
        """
        gi = self.grade_index
        if gi.device != mv.device:
            gi = gi.to(mv.device)
        batch_shape = mv.shape[:-1]
        sq = mv.pow(2)
        flat = sq.reshape(-1, self.dim)
        idx = gi.unsqueeze(0).expand_as(flat)
        result = torch.zeros(flat.shape[0], self.num_grades, device=mv.device, dtype=mv.dtype)
        result.scatter_add_(1, idx, flat)
        return result.reshape(*batch_shape, self.num_grades).clamp(min=1e-6).sqrt()

    def _generate_cayley_table(self):
        """Precompute the Cayley table, grade masks, and reversion signs.

        Returns:
            tuple: Cached tensors for algebra operations.
        """
        indices = torch.arange(self.dim, device=self.device)

        # Result index = A XOR B
        cayley_indices = indices.unsqueeze(0) ^ indices.unsqueeze(1)
        cayley_signs = self._compute_signs(indices)

        # Precompute signs for geometric_product accumulation
        gp_signs = torch.gather(cayley_signs, 1, cayley_indices)

        # Grade masks: one bool tensor per grade (cached to avoid per-call Python loop)
        grade_masks = []
        for k in range(self.n + 1):
            mask = torch.tensor(
                [bin(i).count('1') == k for i in range(self.dim)],
                dtype=torch.bool, device=self.device,
            )
            grade_masks.append(mask)

        # Reverse signs: blade i gets sign (-1)^(k(k-1)/2) where k = grade(i)
        rev_signs = torch.zeros(self.dim, dtype=cayley_signs.dtype, device=self.device)
        for i in range(self.dim):
            k = bin(i).count('1')
            rev_signs[i] = (-1) ** (k * (k - 1) // 2)

        # Bivector squared scalars: for each basis bivector e_ab,
        # (e_ab)^2 = -s_a * s_b where s_i = +1 if i < p, -1 if p <= i < p+q, 0 if i >= p+q.
        # Used by closed-form exp for arbitrary signature.
        if self.n >= 2:
            bv_mask = grade_masks[2]
            bv_indices = bv_mask.nonzero(as_tuple=False).squeeze(-1)
            bv_sq_scalar = torch.zeros(len(bv_indices), dtype=cayley_signs.dtype,
                                       device=self.device)
            for idx_pos, blade_idx in enumerate(bv_indices.tolist()):
                bits = []
                for bit in range(self.n):
                    if blade_idx & (1 << bit):
                        bits.append(bit)
                if len(bits) == 2:
                    a, b = bits
                    s_a = 1.0 if a < self.p else (-1.0 if a < self.p + self.q else 0.0)
                    s_b = 1.0 if b < self.p else (-1.0 if b < self.p + self.q else 0.0)
                    bv_sq_scalar[idx_pos] = -s_a * s_b
        else:
            bv_sq_scalar = torch.zeros(0, dtype=cayley_signs.dtype,
                                       device=self.device)

        # Grade index: maps each basis element to its grade (popcount)
        grade_index = torch.tensor([bin(i).count('1') for i in range(self.dim)],
                                   dtype=torch.long, device=self.device)

        # Precomputed signs for single-pass wedge and inner product
        # wedge(A,B) = (AB - BA)/2 uses antisymmetric part of signs
        # inner(A,B) = (AB + BA)/2 uses symmetric part of signs
        wedge_cayley_signs = (cayley_signs - cayley_signs.T) / 2.0
        inner_cayley_signs = (cayley_signs + cayley_signs.T) / 2.0
        wedge_gp_signs = torch.gather(wedge_cayley_signs, 1, cayley_indices)
        inner_gp_signs = torch.gather(inner_cayley_signs, 1, cayley_indices)

        # Precomputed right-contraction action matrices for bivector-vector case
        # rc_action[b, i, j] encodes how basis bivector b acts on grade-1 vectors
        if self.n >= 2:
            bv_mask_idx = bv_indices.tolist()
            n = self.n
            rc_action = torch.zeros(len(bv_mask_idx), n, n,
                                    dtype=cayley_signs.dtype, device=self.device)
            for idx_pos, blade_idx in enumerate(bv_mask_idx):
                bits = []
                for bit in range(n):
                    if blade_idx & (1 << bit):
                        bits.append(bit)
                if len(bits) == 2:
                    a, b_bit = bits
                    s_a = 1.0 if a < self.p else (-1.0 if a < self.p + self.q else 0.0)
                    s_b = 1.0 if b_bit < self.p else (-1.0 if b_bit < self.p + self.q else 0.0)
                    # e_{ab} . e_b = s_b * e_a  (right contraction)
                    rc_action[idx_pos, a, b_bit] = s_b
                    # e_{ab} . e_a = -s_a * e_b
                    rc_action[idx_pos, b_bit, a] = -s_a
        else:
            rc_action = torch.zeros(0, self.n, self.n,
                                    dtype=cayley_signs.dtype, device=self.device)

        return (cayley_indices, cayley_signs, gp_signs, grade_masks,
                rev_signs, bv_sq_scalar, wedge_gp_signs, inner_gp_signs,
                grade_index, rc_action)

    def _compute_signs(self, indices: torch.Tensor) -> torch.Tensor:
        """Compute the sign matrix from commutation parity and metric signature.

        Handles three signature types:
        - Positive (i < p): e_i^2 = +1
        - Negative (p <= i < p+q): e_i^2 = -1
        - Null (i >= p+q): e_i^2 = 0

        Args:
            indices (torch.Tensor): Basis indices.

        Returns:
            torch.Tensor: Sign matrix.
        """
        # 1. Commutation Sign: Count swaps needed to reorder basis vectors
        # A bit-wise comparison counts inversions
        A = indices.unsqueeze(1) # Row
        B = indices.unsqueeze(0) # Col

        swap_counts = torch.zeros((self.dim, self.dim), dtype=torch.long, device=self.device)
        for i in range(self.n):
            a_i = (A >> i) & 1
            # Count set bits in B strictly lower than bit i
            lower_mask = (1 << i) - 1
            b_lower = (B & lower_mask)

            # Count bits in b_lower
            b_lower_cnt = torch.zeros_like(B)
            temp_b = b_lower
            for _ in range(self.n):
                b_lower_cnt += (temp_b & 1)
                temp_b = temp_b >> 1

            swap_counts += a_i * b_lower_cnt

        commutator_sign = (-1) ** swap_counts

        # 2. Metric Sign: e_i^2 = -1 if p <= i < p+q, 0 if i >= p+q
        intersection = A & B

        # Mask for negative signature dimensions (p <= i < p+q)
        q_mask = 0
        for i in range(self.p, self.p + self.q):
            q_mask |= (1 << i)

        neg_intersection = intersection & q_mask

        # Count set bits in negative intersection
        neg_cnt = torch.zeros_like(neg_intersection)
        temp_neg = neg_intersection
        for _ in range(self.n):
            neg_cnt += (temp_neg & 1)
            temp_neg = temp_neg >> 1

        metric_sign = (-1) ** neg_cnt

        # 3. Null dimensions: if any null basis vector appears in the intersection
        # (i.e., e_i^2 = 0 for i >= p+q), the entire product is killed.
        if self.r > 0:
            r_mask = 0
            for i in range(self.p + self.q, self.n):
                r_mask |= (1 << i)
            null_intersection = intersection & r_mask
            # Any set bit in null_intersection means a null vector is squared -> 0
            has_null = (null_intersection != 0).to(metric_sign.dtype)
            metric_sign = metric_sign * (1 - has_null)

        return (commutator_sign * metric_sign).to(dtype=torch.float32)

    def geometric_product(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
        """Computes the Geometric Product.

        Uses vectorized gather + broadcast multiply + sum. No Python loops.

        Args:
            A (torch.Tensor): Left operand [..., Dim].
            B (torch.Tensor): Right operand [..., Dim].

        Returns:
            torch.Tensor: The product AB [..., Dim].
        """
        from core.validation import check_multivector
        check_multivector(A, self, "geometric_product(A)")
        check_multivector(B, self, "geometric_product(B)")

        # Gather B coefficients according to Cayley table: B_gathered[..., i, k] = B[..., cayley[i,k]]
        idx = self.cayley_indices  # [D, D]
        if idx.device != A.device:
            self.ensure_device(A.device)
            idx = self.cayley_indices

        # Expand B for gather: [..., D] -> [..., D, D] via advanced indexing
        B_gathered = B[..., idx]  # [..., D, D]

        # result[..., k] = sum_i A[..., i] * B[..., cayley[i,k]] * signs[i,k]
        # = sum_i (A[..., i, None] * B_gathered[..., i, :] * signs[i, :])  summed over i
        return (A.unsqueeze(-1) * B_gathered * self.gp_signs).sum(dim=-2)

    def ensure_device(self, device) -> None:
        """Move cached tables to the given device if not already there.

        Args:
            device (torch.device): Target device.
        """
        if self.cayley_indices.device == device:
            return
        self.cayley_indices = self.cayley_indices.to(device)
        self.cayley_signs = self.cayley_signs.to(device)
        self.gp_signs = self.gp_signs.to(device)
        self.grade_masks = [m.to(device) for m in self.grade_masks]
        self.rev_signs = self.rev_signs.to(device)
        self.bv_sq_scalar = self.bv_sq_scalar.to(device)
        self.wedge_gp_signs = self.wedge_gp_signs.to(device)
        self.inner_gp_signs = self.inner_gp_signs.to(device)
        self.grade_index = self.grade_index.to(device)
        self.rc_action = self.rc_action.to(device)

        cache_key = (self.p, self.q, self.r, str(self.device))
        CliffordAlgebra._CACHED_TABLES[cache_key] = (
            self.cayley_indices, self.cayley_signs, self.gp_signs,
            self.grade_masks, self.rev_signs, self.bv_sq_scalar,
            self.wedge_gp_signs, self.inner_gp_signs,
            self.grade_index, self.rc_action,
        )

    def grade_projection(self, mv: torch.Tensor, grade: int) -> torch.Tensor:
        """Isolates a specific grade.

        Args:
            mv (torch.Tensor): Multivector [..., Dim].
            grade (int): Target grade.

        Returns:
            torch.Tensor: Projected multivector [..., Dim].
        """
        mask = self.grade_masks[grade]
        if mask.device != mv.device:
            mask = mask.to(mv.device)
        result = torch.zeros_like(mv)
        result[..., mask] = mv[..., mask]
        return result

    def reverse(self, mv: torch.Tensor) -> torch.Tensor:
        """Computes the reversion. The Clifford conjugate.

        Args:
            mv (torch.Tensor): Input multivector [..., Dim].

        Returns:
            torch.Tensor: Reversed multivector [..., Dim].
        """
        rev = self.rev_signs
        if rev.device != mv.device or rev.dtype != mv.dtype:
            rev = rev.to(dtype=mv.dtype, device=mv.device)
        return mv * rev

    def wedge(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
        """Computes the wedge (outer) product: A ^ B = (AB - BA)/2.

        Single-pass implementation using precomputed antisymmetric signs.

        Reference:
            Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers
            from Irreducibles." arXiv:2507.11688v1 [cs.LG]

        Args:
            A (torch.Tensor): Left operand [..., dim].
            B (torch.Tensor): Right operand [..., dim].

        Returns:
            torch.Tensor: Wedge product A ^ B [..., dim].
        """
        idx = self.cayley_indices
        if idx.device != A.device:
            self.ensure_device(A.device)
            idx = self.cayley_indices
        B_gathered = B[..., idx]
        return (A.unsqueeze(-1) * B_gathered * self.wedge_gp_signs).sum(dim=-2)

    def right_contraction(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
        """Computes the right contraction: A _| B.

        Fast path for bivector-vector case using precomputed skew-symmetric
        action matrices (avoids full geometric product + grade projection).

        Reference:
            Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers
            from Irreducibles." arXiv:2507.11688v1 [cs.LG], Algorithm 2

        Args:
            A (torch.Tensor): Left operand (bivector) [..., dim].
            B (torch.Tensor): Right operand (vector) [..., dim].

        Returns:
            torch.Tensor: Right contraction A _| B [..., dim].
        """
        bv_mask = self.grade_masks[2]
        if bv_mask.device != A.device:
            self.ensure_device(A.device)
            bv_mask = self.grade_masks[2]

        bv_coeffs = A[..., bv_mask]                          # [..., num_bv]
        g1_idx = self.grade_masks[1].nonzero(as_tuple=False).squeeze(-1)
        v_coeffs = B[..., g1_idx]                            # [..., n]

        rc = self.rc_action
        if rc.device != A.device or rc.dtype != A.dtype:
            rc = rc.to(device=A.device, dtype=A.dtype)

        # M[..., i, j] = sum_b bv_coeffs[..., b] * rc_action[b, i, j]
        M = torch.einsum('...b, bij -> ...ij', bv_coeffs, rc)  # [..., n, n]
        result_v = torch.matmul(M, v_coeffs.unsqueeze(-1)).squeeze(-1)  # [..., n]

        result = torch.zeros_like(A)
        result[..., g1_idx] = result_v
        return result

    def inner_product(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
        """Computes the inner product: A . B = (AB + BA)/2.

        Single-pass implementation using precomputed symmetric signs.

        Reference:
            Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers
            from Irreducibles." arXiv:2507.11688v1 [cs.LG]

        Args:
            A (torch.Tensor): Left operand [..., dim].
            B (torch.Tensor): Right operand [..., dim].

        Returns:
            torch.Tensor: Inner product A . B [..., dim].
        """
        idx = self.cayley_indices
        if idx.device != A.device:
            self.ensure_device(A.device)
            idx = self.cayley_indices
        B_gathered = B[..., idx]
        return (A.unsqueeze(-1) * B_gathered * self.inner_gp_signs).sum(dim=-2)

    def blade_inverse(self, blade: torch.Tensor) -> torch.Tensor:
        """Compute the inverse of a blade: B^{-1} = B_rev / <B * B_rev>_0.

        Works for any simple blade (non-degenerate). For null blades the
        scalar denominator is clamped to avoid division by zero.

        Args:
            blade (torch.Tensor): Blade multivector [..., dim].

        Returns:
            torch.Tensor: Inverse blade [..., dim].
        """
        blade_rev = self.reverse(blade)
        blade_sq = self.geometric_product(blade, blade_rev)
        scalar = blade_sq[..., 0:1].clamp(min=1e-12)
        return blade_rev / scalar

    def blade_project(self, mv: torch.Tensor, blade: torch.Tensor) -> torch.Tensor:
        """Project multivector onto blade subspace: (mv . B) B^{-1}.

        Args:
            mv (torch.Tensor): Multivector to project [..., dim].
            blade (torch.Tensor): Blade defining the subspace [..., dim].

        Returns:
            torch.Tensor: Projected multivector [..., dim].
        """
        inner = self.inner_product(mv, blade)
        return self.geometric_product(inner, self.blade_inverse(blade))

    def blade_reject(self, mv: torch.Tensor, blade: torch.Tensor) -> torch.Tensor:
        """Reject multivector from blade subspace: mv - proj_B(mv).

        The orthogonal complement of the projection onto blade.

        Args:
            mv (torch.Tensor): Multivector to reject [..., dim].
            blade (torch.Tensor): Blade defining the subspace [..., dim].

        Returns:
            torch.Tensor: Rejected multivector [..., dim].
        """
        return mv - self.blade_project(mv, blade)

    def grade_involution(self, mv: torch.Tensor) -> torch.Tensor:
        """Grade involution (main involution): x_hat = sum (-1)^k <x>_k.

        Flips sign of all odd-grade components, preserves even-grade.
        This is an algebra automorphism: (AB)^ = A_hat B_hat.

        Args:
            mv (torch.Tensor): Input multivector [..., dim].

        Returns:
            torch.Tensor: Involuted multivector [..., dim].
        """
        gi = self.grade_index
        if gi.device != mv.device:
            gi = gi.to(mv.device)
        signs = (-1.0) ** gi.float()
        return mv * signs.to(dtype=mv.dtype)

    def exp(self, mv: torch.Tensor) -> torch.Tensor:
        """Exponentiates a bivector to produce a rotor via closed-form.

        Three signature regimes:
            - Elliptic  (B^2 < 0): exp(B) = cos(th) + sin(th)/th * B
            - Hyperbolic (B^2 > 0): exp(B) = cosh(th) + sinh(th)/th * B
            - Parabolic  (B^2 ~ 0): exp(B) ~ 1 + B

        For n <= 3 every bivector is simple, so the closed-form is exact.
        For n >= 4 the closed-form is exact for simple bivectors and a
        first-order approximation for non-simple ones.  Use
        ``exp_decomposed()`` when exact non-simple handling is needed
        (inference only -- see ``core.decomposition``).

        Args:
            mv (torch.Tensor): Pure bivector [..., dim].

        Returns:
            torch.Tensor: Rotor exp(mv) [..., dim].
        """
        return self._exp_bivector_closed(mv)

    def _exp_bivector_closed(self, B: torch.Tensor) -> torch.Tensor:
        """Closed-form exponential for bivectors in arbitrary signature.

        For a bivector B, computes B^2 (scalar part) using the metric:
            B^2_scalar = Sum_k b_k^2 . (e_k)^2   where (e_ab)^2 = -s_a.s_b

        Three regimes:
            - B^2 < 0 (elliptic): exp(B) = cos(theta) + sin(theta)/theta . B,  theta = Sqrt(-B^2)
            - B^2 > 0 (hyperbolic): exp(B) = cosh(theta) + sinh(theta)/theta . B,  theta = Sqrt(B^2)
            - B^2 ~= 0 (parabolic): exp(B) ~= 1 + B

        Uses zero geometric products. Exact for simple bivectors in any
        Clifford algebra Cl(p,q,r).

        Args:
            B (torch.Tensor): Pure bivector [..., dim].

        Returns:
            torch.Tensor: Rotor exp(B) [..., dim].
        """
        bv_mask = self.grade_masks[2]
        if bv_mask.device != B.device:
            bv_mask = bv_mask.to(B.device)
        bv_coeffs = B[..., bv_mask]  # [..., num_bivectors]

        bv_sq = self.bv_sq_scalar
        if bv_sq.device != B.device:
            bv_sq = bv_sq.to(B.device)

        # Signed squared norm: alpha = Sum_k b_k^2 . (e_k)^2
        # alpha < 0 -> elliptic (Euclidean-like), alpha > 0 -> hyperbolic
        alpha = (bv_coeffs * bv_coeffs * bv_sq).sum(dim=-1, keepdim=True)

        abs_alpha = alpha.abs().clamp(min=1e-12)
        theta = torch.sqrt(abs_alpha)  # [..., 1]

        # Elliptic branch: cos(theta) and sin(theta)/theta
        cos_theta = torch.cos(theta)
        sinc_theta = torch.where(
            theta > 1e-7,
            torch.sin(theta) / theta,
            1.0 - abs_alpha / 6.0,
        )

        # Hyperbolic branch: cosh(theta) and sinh(theta)/theta
        cosh_theta = torch.cosh(theta)
        sinhc_theta = torch.where(
            theta > 1e-7,
            torch.sinh(theta) / theta,
            1.0 + abs_alpha / 6.0,
        )

        # Select branch based on sign of alpha
        is_elliptic = alpha < -1e-12
        is_hyperbolic = alpha > 1e-12
        # Parabolic (null) falls through: scalar=1, coeff=1

        scalar_part = torch.where(
            is_elliptic, cos_theta,
            torch.where(is_hyperbolic, cosh_theta, torch.ones_like(theta))
        )
        coeff_part = torch.where(
            is_elliptic, sinc_theta,
            torch.where(is_hyperbolic, sinhc_theta, torch.ones_like(theta))
        )

        result = coeff_part * B
        result[..., 0] = scalar_part.squeeze(-1)

        return result

    def _exp_taylor(self, mv: torch.Tensor, order: int = 8) -> torch.Tensor:
        """Taylor series exponential with scaling-and-squaring (fallback).

        Args:
            mv (torch.Tensor): General multivector [..., dim].
            order (int, optional): Taylor order. Defaults to 8.

        Returns:
            torch.Tensor: exp(mv) [..., dim].
        """
        norm = mv.norm(dim=-1, keepdim=True)
        k = torch.ceil(torch.log2(torch.clamp(norm, min=1.0))).int()

        max_k = k.max().item()
        if max_k > 0:
            mv_scaled = mv / (2.0 ** max_k)
        else:
            mv_scaled = mv

        res = torch.zeros_like(mv)
        res[..., 0] = 1.0

        term = torch.zeros_like(mv)
        term[..., 0] = 1.0

        for i in range(1, order + 1):
            term = self.geometric_product(term, mv_scaled)
            res = res + term / math.factorial(i)

        if max_k > 0:
            for _ in range(int(max_k)):
                res = self.geometric_product(res, res)

        return res

    def exp_decomposed(self, mv: torch.Tensor, **kwargs) -> torch.Tensor:
        """Exponentiates a bivector via decomposition into simple components.

        Decomposes a general bivector into a sum of simple (blade) bivectors
        using GA power iteration, exponentiates each in closed form, then
        composes via geometric product.  Exact for all bivectors, but the
        power iteration loop is not differentiable through its normalization
        steps.  During training the loop is detached and gradients flow only
        through the final closed-form exp + composition.

        Reference:
            Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear
            Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG]

        Args:
            mv (torch.Tensor): Pure bivector [..., dim].
            **kwargs (dict): Forwarded to ``core.decomposition.exp_decomposed``.
                use_decomposition (bool): Enable decomposition. Default True.
                k (int | None): Number of simple components (auto if None).
                threshold (float): Convergence threshold. Default 1e-6.
                max_iterations (int): Power iteration cap. Default 100.

        Returns:
            torch.Tensor: Rotor exp(mv) [..., dim].
        """
        from core.decomposition import exp_decomposed
        if 'use_decomposition' not in kwargs:
            kwargs['use_decomposition'] = True
        return exp_decomposed(self, mv, **kwargs)

num_grades property

Counts the number of grades (n + 1).

Returns:

Name	Type	Description
int	int	Number of grades.

__init__(p, q=0, r=0, device='cuda')

Initialize the algebra and cache the Cayley table.

Parameters:

Name	Type	Description	Default
p	int	Positive dimensions (+1).	required
q	int	Negative dimensions (-1). Defaults to 0.	0
r	int	Degenerate dimensions (0). Defaults to 0.	0
device	str	The device on which computations are performed. Defaults to 'cuda'.	'cuda'

Source code in core/algebra.py

def __init__(self, p: int, q: int = 0, r: int = 0, device='cuda'):
    """Initialize the algebra and cache the Cayley table.

    Args:
        p (int): Positive dimensions (+1).
        q (int, optional): Negative dimensions (-1). Defaults to 0.
        r (int, optional): Degenerate dimensions (0). Defaults to 0.
        device (str, optional): The device on which computations are performed. Defaults to 'cuda'.
    """
    assert p >= 0, f"p must be non-negative, got {p}"
    assert q >= 0, f"q must be non-negative, got {q}"
    assert r >= 0, f"r must be non-negative, got {r}"
    assert p + q + r <= 12, f"p + q + r must be <= 12, got {p + q + r}"

    self.p, self.q, self.r = p, q, r
    self.n = p + q + r
    self.dim = 2 ** self.n
    self.device = device

    # Cache Cayley tables to avoid recomputation
    cache_key = (p, q, r, str(device))
    if cache_key not in CliffordAlgebra._CACHED_TABLES:
        CliffordAlgebra._CACHED_TABLES[cache_key] = self._generate_cayley_table()

    (
        self.cayley_indices,
        self.cayley_signs,
        self.gp_signs,
        self.grade_masks,
        self.rev_signs,
        self.bv_sq_scalar,
        self.wedge_gp_signs,
        self.inner_gp_signs,
        self.grade_index,
        self.rc_action,
    ) = CliffordAlgebra._CACHED_TABLES[cache_key]

embed_vector(vectors)

Injects vectors into the Grade-1 subspace.

Parameters:

Name	Type	Description	Default
vectors	Tensor	Raw vectors [..., n].	required

Returns:

Type	Description
Tensor	torch.Tensor: Multivector coefficients [..., dim].

Source code in core/algebra.py

def embed_vector(self, vectors: torch.Tensor) -> torch.Tensor:
    """Injects vectors into the Grade-1 subspace.

    Args:
        vectors (torch.Tensor): Raw vectors [..., n].

    Returns:
        torch.Tensor: Multivector coefficients [..., dim].
    """
    batch_shape = vectors.shape[:-1]
    mv = torch.zeros(*batch_shape, self.dim, device=vectors.device, dtype=vectors.dtype)
    for i in range(self.n):
        mv[..., 1 << i] = vectors[..., i]
    return mv

get_grade_norms(mv)

Calculates norms per grade. Useful for invariant features.

Vectorized via scatter_add -- no Python loops over grades.

Parameters:

Name	Type	Description	Default
mv	Tensor	Input multivector [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Grade norms [..., num_grades].

Source code in core/algebra.py

def get_grade_norms(self, mv: torch.Tensor) -> torch.Tensor:
    """Calculates norms per grade. Useful for invariant features.

    Vectorized via scatter_add -- no Python loops over grades.

    Args:
        mv (torch.Tensor): Input multivector [..., dim].

    Returns:
        torch.Tensor: Grade norms [..., num_grades].
    """
    gi = self.grade_index
    if gi.device != mv.device:
        gi = gi.to(mv.device)
    batch_shape = mv.shape[:-1]
    sq = mv.pow(2)
    flat = sq.reshape(-1, self.dim)
    idx = gi.unsqueeze(0).expand_as(flat)
    result = torch.zeros(flat.shape[0], self.num_grades, device=mv.device, dtype=mv.dtype)
    result.scatter_add_(1, idx, flat)
    return result.reshape(*batch_shape, self.num_grades).clamp(min=1e-6).sqrt()

geometric_product(A, B)

Computes the Geometric Product.

Uses vectorized gather + broadcast multiply + sum. No Python loops.

Parameters:

Name	Type	Description	Default
A	Tensor	Left operand [..., Dim].	required
B	Tensor	Right operand [..., Dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: The product AB [..., Dim].

Source code in core/algebra.py

def geometric_product(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Computes the Geometric Product.

    Uses vectorized gather + broadcast multiply + sum. No Python loops.

    Args:
        A (torch.Tensor): Left operand [..., Dim].
        B (torch.Tensor): Right operand [..., Dim].

    Returns:
        torch.Tensor: The product AB [..., Dim].
    """
    from core.validation import check_multivector
    check_multivector(A, self, "geometric_product(A)")
    check_multivector(B, self, "geometric_product(B)")

    # Gather B coefficients according to Cayley table: B_gathered[..., i, k] = B[..., cayley[i,k]]
    idx = self.cayley_indices  # [D, D]
    if idx.device != A.device:
        self.ensure_device(A.device)
        idx = self.cayley_indices

    # Expand B for gather: [..., D] -> [..., D, D] via advanced indexing
    B_gathered = B[..., idx]  # [..., D, D]

    # result[..., k] = sum_i A[..., i] * B[..., cayley[i,k]] * signs[i,k]
    # = sum_i (A[..., i, None] * B_gathered[..., i, :] * signs[i, :])  summed over i
    return (A.unsqueeze(-1) * B_gathered * self.gp_signs).sum(dim=-2)

ensure_device(device)

Move cached tables to the given device if not already there.

Parameters:

Name	Type	Description	Default
device	device	Target device.	required

Source code in core/algebra.py

def ensure_device(self, device) -> None:
    """Move cached tables to the given device if not already there.

    Args:
        device (torch.device): Target device.
    """
    if self.cayley_indices.device == device:
        return
    self.cayley_indices = self.cayley_indices.to(device)
    self.cayley_signs = self.cayley_signs.to(device)
    self.gp_signs = self.gp_signs.to(device)
    self.grade_masks = [m.to(device) for m in self.grade_masks]
    self.rev_signs = self.rev_signs.to(device)
    self.bv_sq_scalar = self.bv_sq_scalar.to(device)
    self.wedge_gp_signs = self.wedge_gp_signs.to(device)
    self.inner_gp_signs = self.inner_gp_signs.to(device)
    self.grade_index = self.grade_index.to(device)
    self.rc_action = self.rc_action.to(device)

    cache_key = (self.p, self.q, self.r, str(self.device))
    CliffordAlgebra._CACHED_TABLES[cache_key] = (
        self.cayley_indices, self.cayley_signs, self.gp_signs,
        self.grade_masks, self.rev_signs, self.bv_sq_scalar,
        self.wedge_gp_signs, self.inner_gp_signs,
        self.grade_index, self.rc_action,
    )

grade_projection(mv, grade)

Isolates a specific grade.

Parameters:

Name	Type	Description	Default
mv	Tensor	Multivector [..., Dim].	required
grade	int	Target grade.	required

Returns:

Type	Description
Tensor	torch.Tensor: Projected multivector [..., Dim].

Source code in core/algebra.py

def grade_projection(self, mv: torch.Tensor, grade: int) -> torch.Tensor:
    """Isolates a specific grade.

    Args:
        mv (torch.Tensor): Multivector [..., Dim].
        grade (int): Target grade.

    Returns:
        torch.Tensor: Projected multivector [..., Dim].
    """
    mask = self.grade_masks[grade]
    if mask.device != mv.device:
        mask = mask.to(mv.device)
    result = torch.zeros_like(mv)
    result[..., mask] = mv[..., mask]
    return result

reverse(mv)

Computes the reversion. The Clifford conjugate.

Parameters:

Name	Type	Description	Default
mv	Tensor	Input multivector [..., Dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Reversed multivector [..., Dim].

Source code in core/algebra.py

def reverse(self, mv: torch.Tensor) -> torch.Tensor:
    """Computes the reversion. The Clifford conjugate.

    Args:
        mv (torch.Tensor): Input multivector [..., Dim].

    Returns:
        torch.Tensor: Reversed multivector [..., Dim].
    """
    rev = self.rev_signs
    if rev.device != mv.device or rev.dtype != mv.dtype:
        rev = rev.to(dtype=mv.dtype, device=mv.device)
    return mv * rev

wedge(A, B)

Computes the wedge (outer) product: A ^ B = (AB - BA)/2.

Single-pass implementation using precomputed antisymmetric signs.

Reference

Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG]

Parameters:

Name	Type	Description	Default
A	Tensor	Left operand [..., dim].	required
B	Tensor	Right operand [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Wedge product A ^ B [..., dim].

Source code in core/algebra.py

def wedge(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Computes the wedge (outer) product: A ^ B = (AB - BA)/2.

    Single-pass implementation using precomputed antisymmetric signs.

    Reference:
        Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers
        from Irreducibles." arXiv:2507.11688v1 [cs.LG]

    Args:
        A (torch.Tensor): Left operand [..., dim].
        B (torch.Tensor): Right operand [..., dim].

    Returns:
        torch.Tensor: Wedge product A ^ B [..., dim].
    """
    idx = self.cayley_indices
    if idx.device != A.device:
        self.ensure_device(A.device)
        idx = self.cayley_indices
    B_gathered = B[..., idx]
    return (A.unsqueeze(-1) * B_gathered * self.wedge_gp_signs).sum(dim=-2)

right_contraction(A, B)

Computes the right contraction: A _| B.

Fast path for bivector-vector case using precomputed skew-symmetric action matrices (avoids full geometric product + grade projection).

Reference

Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG], Algorithm 2

Parameters:

Name	Type	Description	Default
A	Tensor	Left operand (bivector) [..., dim].	required
B	Tensor	Right operand (vector) [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Right contraction A _| B [..., dim].

Source code in core/algebra.py

def right_contraction(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Computes the right contraction: A _| B.

    Fast path for bivector-vector case using precomputed skew-symmetric
    action matrices (avoids full geometric product + grade projection).

    Reference:
        Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers
        from Irreducibles." arXiv:2507.11688v1 [cs.LG], Algorithm 2

    Args:
        A (torch.Tensor): Left operand (bivector) [..., dim].
        B (torch.Tensor): Right operand (vector) [..., dim].

    Returns:
        torch.Tensor: Right contraction A _| B [..., dim].
    """
    bv_mask = self.grade_masks[2]
    if bv_mask.device != A.device:
        self.ensure_device(A.device)
        bv_mask = self.grade_masks[2]

    bv_coeffs = A[..., bv_mask]                          # [..., num_bv]
    g1_idx = self.grade_masks[1].nonzero(as_tuple=False).squeeze(-1)
    v_coeffs = B[..., g1_idx]                            # [..., n]

    rc = self.rc_action
    if rc.device != A.device or rc.dtype != A.dtype:
        rc = rc.to(device=A.device, dtype=A.dtype)

    # M[..., i, j] = sum_b bv_coeffs[..., b] * rc_action[b, i, j]
    M = torch.einsum('...b, bij -> ...ij', bv_coeffs, rc)  # [..., n, n]
    result_v = torch.matmul(M, v_coeffs.unsqueeze(-1)).squeeze(-1)  # [..., n]

    result = torch.zeros_like(A)
    result[..., g1_idx] = result_v
    return result

inner_product(A, B)

Computes the inner product: A . B = (AB + BA)/2.

Single-pass implementation using precomputed symmetric signs.

Reference

Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG]

Parameters:

Name	Type	Description	Default
A	Tensor	Left operand [..., dim].	required
B	Tensor	Right operand [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Inner product A . B [..., dim].

Source code in core/algebra.py

def inner_product(self, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Computes the inner product: A . B = (AB + BA)/2.

    Single-pass implementation using precomputed symmetric signs.

    Reference:
        Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers
        from Irreducibles." arXiv:2507.11688v1 [cs.LG]

    Args:
        A (torch.Tensor): Left operand [..., dim].
        B (torch.Tensor): Right operand [..., dim].

    Returns:
        torch.Tensor: Inner product A . B [..., dim].
    """
    idx = self.cayley_indices
    if idx.device != A.device:
        self.ensure_device(A.device)
        idx = self.cayley_indices
    B_gathered = B[..., idx]
    return (A.unsqueeze(-1) * B_gathered * self.inner_gp_signs).sum(dim=-2)

blade_inverse(blade)

Compute the inverse of a blade: B^{-1} = B_rev / _0.

Works for any simple blade (non-degenerate). For null blades the scalar denominator is clamped to avoid division by zero.

Parameters:

Name	Type	Description	Default
blade	Tensor	Blade multivector [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Inverse blade [..., dim].

Source code in core/algebra.py

def blade_inverse(self, blade: torch.Tensor) -> torch.Tensor:
    """Compute the inverse of a blade: B^{-1} = B_rev / <B * B_rev>_0.

    Works for any simple blade (non-degenerate). For null blades the
    scalar denominator is clamped to avoid division by zero.

    Args:
        blade (torch.Tensor): Blade multivector [..., dim].

    Returns:
        torch.Tensor: Inverse blade [..., dim].
    """
    blade_rev = self.reverse(blade)
    blade_sq = self.geometric_product(blade, blade_rev)
    scalar = blade_sq[..., 0:1].clamp(min=1e-12)
    return blade_rev / scalar

blade_project(mv, blade)

Project multivector onto blade subspace: (mv . B) B^{-1}.

Parameters:

Name	Type	Description	Default
mv	Tensor	Multivector to project [..., dim].	required
blade	Tensor	Blade defining the subspace [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Projected multivector [..., dim].

Source code in core/algebra.py

def blade_project(self, mv: torch.Tensor, blade: torch.Tensor) -> torch.Tensor:
    """Project multivector onto blade subspace: (mv . B) B^{-1}.

    Args:
        mv (torch.Tensor): Multivector to project [..., dim].
        blade (torch.Tensor): Blade defining the subspace [..., dim].

    Returns:
        torch.Tensor: Projected multivector [..., dim].
    """
    inner = self.inner_product(mv, blade)
    return self.geometric_product(inner, self.blade_inverse(blade))

blade_reject(mv, blade)

Reject multivector from blade subspace: mv - proj_B(mv).

The orthogonal complement of the projection onto blade.

Parameters:

Name	Type	Description	Default
mv	Tensor	Multivector to reject [..., dim].	required
blade	Tensor	Blade defining the subspace [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Rejected multivector [..., dim].

Source code in core/algebra.py

def blade_reject(self, mv: torch.Tensor, blade: torch.Tensor) -> torch.Tensor:
    """Reject multivector from blade subspace: mv - proj_B(mv).

    The orthogonal complement of the projection onto blade.

    Args:
        mv (torch.Tensor): Multivector to reject [..., dim].
        blade (torch.Tensor): Blade defining the subspace [..., dim].

    Returns:
        torch.Tensor: Rejected multivector [..., dim].
    """
    return mv - self.blade_project(mv, blade)

grade_involution(mv)

Grade involution (main involution): x_hat = sum (-1)^k _k.

Flips sign of all odd-grade components, preserves even-grade. This is an algebra automorphism: (AB)^ = A_hat B_hat.

Parameters:

Name	Type	Description	Default
mv	Tensor	Input multivector [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Involuted multivector [..., dim].

Source code in core/algebra.py

def grade_involution(self, mv: torch.Tensor) -> torch.Tensor:
    """Grade involution (main involution): x_hat = sum (-1)^k <x>_k.

    Flips sign of all odd-grade components, preserves even-grade.
    This is an algebra automorphism: (AB)^ = A_hat B_hat.

    Args:
        mv (torch.Tensor): Input multivector [..., dim].

    Returns:
        torch.Tensor: Involuted multivector [..., dim].
    """
    gi = self.grade_index
    if gi.device != mv.device:
        gi = gi.to(mv.device)
    signs = (-1.0) ** gi.float()
    return mv * signs.to(dtype=mv.dtype)

exp(mv)

Exponentiates a bivector to produce a rotor via closed-form.

Three signature regimes

• Elliptic (B^2 < 0): exp(B) = cos(th) + sin(th)/th * B
• Hyperbolic (B^2 > 0): exp(B) = cosh(th) + sinh(th)/th * B
• Parabolic (B^2 ~ 0): exp(B) ~ 1 + B

For n <= 3 every bivector is simple, so the closed-form is exact. For n >= 4 the closed-form is exact for simple bivectors and a first-order approximation for non-simple ones. Use exp_decomposed() when exact non-simple handling is needed (inference only -- see core.decomposition).

Parameters:

Name	Type	Description	Default
mv	Tensor	Pure bivector [..., dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Rotor exp(mv) [..., dim].

Source code in core/algebra.py

def exp(self, mv: torch.Tensor) -> torch.Tensor:
    """Exponentiates a bivector to produce a rotor via closed-form.

    Three signature regimes:
        - Elliptic  (B^2 < 0): exp(B) = cos(th) + sin(th)/th * B
        - Hyperbolic (B^2 > 0): exp(B) = cosh(th) + sinh(th)/th * B
        - Parabolic  (B^2 ~ 0): exp(B) ~ 1 + B

    For n <= 3 every bivector is simple, so the closed-form is exact.
    For n >= 4 the closed-form is exact for simple bivectors and a
    first-order approximation for non-simple ones.  Use
    ``exp_decomposed()`` when exact non-simple handling is needed
    (inference only -- see ``core.decomposition``).

    Args:
        mv (torch.Tensor): Pure bivector [..., dim].

    Returns:
        torch.Tensor: Rotor exp(mv) [..., dim].
    """
    return self._exp_bivector_closed(mv)

exp_decomposed(mv, **kwargs)

Exponentiates a bivector via decomposition into simple components.

Decomposes a general bivector into a sum of simple (blade) bivectors using GA power iteration, exponentiates each in closed form, then composes via geometric product. Exact for all bivectors, but the power iteration loop is not differentiable through its normalization steps. During training the loop is detached and gradients flow only through the final closed-form exp + composition.

Reference

Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG]

Parameters:

Name	Type	Description	Default
mv	Tensor	Pure bivector [..., dim].	required
**kwargs	dict	Forwarded to core.decomposition.exp_decomposed. use_decomposition (bool): Enable decomposition. Default True. k (int | None): Number of simple components (auto if None). threshold (float): Convergence threshold. Default 1e-6. max_iterations (int): Power iteration cap. Default 100.	{}

Returns:

Type	Description
Tensor	torch.Tensor: Rotor exp(mv) [..., dim].

Source code in core/algebra.py

def exp_decomposed(self, mv: torch.Tensor, **kwargs) -> torch.Tensor:
    """Exponentiates a bivector via decomposition into simple components.

    Decomposes a general bivector into a sum of simple (blade) bivectors
    using GA power iteration, exponentiates each in closed form, then
    composes via geometric product.  Exact for all bivectors, but the
    power iteration loop is not differentiable through its normalization
    steps.  During training the loop is detached and gradients flow only
    through the final closed-form exp + composition.

    Reference:
        Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear
        Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG]

    Args:
        mv (torch.Tensor): Pure bivector [..., dim].
        **kwargs (dict): Forwarded to ``core.decomposition.exp_decomposed``.
            use_decomposition (bool): Enable decomposition. Default True.
            k (int | None): Number of simple components (auto if None).
            threshold (float): Convergence threshold. Default 1e-6.
            max_iterations (int): Power iteration cap. Default 100.

    Returns:
        torch.Tensor: Rotor exp(mv) [..., dim].
    """
    from core.decomposition import exp_decomposed
    if 'use_decomposition' not in kwargs:
        kwargs['use_decomposition'] = True
    return exp_decomposed(self, mv, **kwargs)

Conformal Algebra

ConformalAlgebra

Helper for CGA. Handles Conformal Geometric Algebra operations.

Maps Euclidean R^n to Null Cone in R^{n+1, 1}.

Attributes:

Name	Type	Description
d	int	Euclidean dimension.
algebra	CliffordAlgebra	The higher dimensional algebra.
e_o	Tensor	Origin basis vector (null).
e_inf	Tensor	Infinity basis vector (null).
device	str	Computation device.

Source code in core/cga.py

class ConformalAlgebra:
    """Helper for CGA. Handles Conformal Geometric Algebra operations.

    Maps Euclidean R^n to Null Cone in R^{n+1, 1}.

    Attributes:
        d (int): Euclidean dimension.
        algebra (CliffordAlgebra): The higher dimensional algebra.
        e_o (torch.Tensor): Origin basis vector (null).
        e_inf (torch.Tensor): Infinity basis vector (null).
        device (str): Computation device.
    """

    def __init__(self, euclidean_dim: int = 3, device='cpu'):
        """Sets up the CGA stage.

        Args:
            euclidean_dim (int): Physical dimension d.
            device (str): Computation device.
        """
        self.d = euclidean_dim
        # Standard CGA signature: p=d+1 (e1..ed, e+), q=1 (e-)
        self.algebra = CliffordAlgebra(p=euclidean_dim + 1, q=1, device=device)
        self.device = device

        # Basis indices
        # e_i corresponds to bit i-1. e+ is bit d, e- is bit d+1.
        self.idx_ep = 1 << euclidean_dim
        self.idx_em = 1 << (euclidean_dim + 1)

        # Construct e_o and e_inf
        # e_inf = e_- + e_+
        # e_o = 0.5 * (e_- - e_+)
        self.e_o = torch.zeros(1, self.algebra.dim, device=device)
        self.e_inf = torch.zeros(1, self.algebra.dim, device=device)

        self.e_inf[0, self.idx_em] = 1.0
        self.e_inf[0, self.idx_ep] = 1.0

        self.e_o[0, self.idx_em] = 0.5
        self.e_o[0, self.idx_ep] = -0.5

    def to_cga(self, x: torch.Tensor) -> torch.Tensor:
        """Embeds Euclidean points into the conformal null cone.

        P(x) = x + 0.5 * x^2 * e_inf + e_o.

        Args:
            x (torch.Tensor): Euclidean points [Batch, d].

        Returns:
            torch.Tensor: Conformal points [Batch, 2^n].
        """
        batch_size = x.shape[0]

        # 1. Embed vector part x
        x_mv = torch.zeros(batch_size, self.algebra.dim, device=x.device, dtype=x.dtype)
        for i in range(self.d):
            basis_idx = 1 << i
            x_mv[:, basis_idx] = x[:, i]

        # 2. Compute x^2 (squared Euclidean norm)
        x_sq = (x ** 2).sum(dim=1, keepdim=True)

        # 3. Assemble P
        e_inf = self.e_inf.to(device=x.device, dtype=x.dtype)
        e_o = self.e_o.to(device=x.device, dtype=x.dtype)
        term2 = 0.5 * x_sq * e_inf
        term3 = e_o.expand(batch_size, -1)

        P = x_mv + term2 + term3
        return P

    def from_cga(self, P: torch.Tensor) -> torch.Tensor:
        """Projects conformal points back to Euclidean space.

        Normalization: P -> P / (-P . e_inf).

        Args:
            P (torch.Tensor): Conformal points [Batch, 2^n].

        Returns:
            torch.Tensor: Euclidean coordinates [Batch, d].
        """
        # 1. Normalize P such that P inner e_inf = -1
        e_inf = self.e_inf.to(device=P.device, dtype=P.dtype)
        P_einf = self.algebra.geometric_product(P, e_inf)
        scale = -P_einf[..., 0:1] # Scalar part

        P_norm = P / (scale + 1e-6)

        # 2. Extract vector components (indices 0 to d-1)
        x = torch.zeros(P.shape[0], self.d, device=P.device, dtype=P.dtype)
        for i in range(self.d):
            basis_idx = 1 << i
            x[:, i] = P_norm[:, basis_idx]

        return x

__init__(euclidean_dim=3, device='cpu')

Sets up the CGA stage.

Parameters:

Name	Type	Description	Default
euclidean_dim	int	Physical dimension d.	3
device	str	Computation device.	'cpu'

Source code in core/cga.py

def __init__(self, euclidean_dim: int = 3, device='cpu'):
    """Sets up the CGA stage.

    Args:
        euclidean_dim (int): Physical dimension d.
        device (str): Computation device.
    """
    self.d = euclidean_dim
    # Standard CGA signature: p=d+1 (e1..ed, e+), q=1 (e-)
    self.algebra = CliffordAlgebra(p=euclidean_dim + 1, q=1, device=device)
    self.device = device

    # Basis indices
    # e_i corresponds to bit i-1. e+ is bit d, e- is bit d+1.
    self.idx_ep = 1 << euclidean_dim
    self.idx_em = 1 << (euclidean_dim + 1)

    # Construct e_o and e_inf
    # e_inf = e_- + e_+
    # e_o = 0.5 * (e_- - e_+)
    self.e_o = torch.zeros(1, self.algebra.dim, device=device)
    self.e_inf = torch.zeros(1, self.algebra.dim, device=device)

    self.e_inf[0, self.idx_em] = 1.0
    self.e_inf[0, self.idx_ep] = 1.0

    self.e_o[0, self.idx_em] = 0.5
    self.e_o[0, self.idx_ep] = -0.5

to_cga(x)

Embeds Euclidean points into the conformal null cone.

P(x) = x + 0.5 * x^2 * e_inf + e_o.

Parameters:

Name	Type	Description	Default
x	Tensor	Euclidean points [Batch, d].	required

Returns:

Type	Description
Tensor	torch.Tensor: Conformal points [Batch, 2^n].

Source code in core/cga.py

def to_cga(self, x: torch.Tensor) -> torch.Tensor:
    """Embeds Euclidean points into the conformal null cone.

    P(x) = x + 0.5 * x^2 * e_inf + e_o.

    Args:
        x (torch.Tensor): Euclidean points [Batch, d].

    Returns:
        torch.Tensor: Conformal points [Batch, 2^n].
    """
    batch_size = x.shape[0]

    # 1. Embed vector part x
    x_mv = torch.zeros(batch_size, self.algebra.dim, device=x.device, dtype=x.dtype)
    for i in range(self.d):
        basis_idx = 1 << i
        x_mv[:, basis_idx] = x[:, i]

    # 2. Compute x^2 (squared Euclidean norm)
    x_sq = (x ** 2).sum(dim=1, keepdim=True)

    # 3. Assemble P
    e_inf = self.e_inf.to(device=x.device, dtype=x.dtype)
    e_o = self.e_o.to(device=x.device, dtype=x.dtype)
    term2 = 0.5 * x_sq * e_inf
    term3 = e_o.expand(batch_size, -1)

    P = x_mv + term2 + term3
    return P

from_cga(P)

Projects conformal points back to Euclidean space.

Normalization: P -> P / (-P . e_inf).

Parameters:

Name	Type	Description	Default
P	Tensor	Conformal points [Batch, 2^n].	required

Returns:

Type	Description
Tensor	torch.Tensor: Euclidean coordinates [Batch, d].

Source code in core/cga.py

def from_cga(self, P: torch.Tensor) -> torch.Tensor:
    """Projects conformal points back to Euclidean space.

    Normalization: P -> P / (-P . e_inf).

    Args:
        P (torch.Tensor): Conformal points [Batch, 2^n].

    Returns:
        torch.Tensor: Euclidean coordinates [Batch, d].
    """
    # 1. Normalize P such that P inner e_inf = -1
    e_inf = self.e_inf.to(device=P.device, dtype=P.dtype)
    P_einf = self.algebra.geometric_product(P, e_inf)
    scale = -P_einf[..., 0:1] # Scalar part

    P_norm = P / (scale + 1e-6)

    # 2. Extract vector components (indices 0 to d-1)
    x = torch.zeros(P.shape[0], self.d, device=P.device, dtype=P.dtype)
    for i in range(self.d):
        basis_idx = 1 << i
        x[:, i] = P_norm[:, basis_idx]

    return x

Multivector

Multivector

Object-oriented multivector wrapper with operator overloading.

Attributes:

Name	Type	Description
algebra	CliffordAlgebra	The backend.
tensor	Tensor	The raw data [..., Dim].

Source code in core/multivector.py

class Multivector:
    """Object-oriented multivector wrapper with operator overloading.

    Attributes:
        algebra (CliffordAlgebra): The backend.
        tensor (torch.Tensor): The raw data [..., Dim].
    """

    def __init__(self, algebra: CliffordAlgebra, tensor: torch.Tensor):
        """Wraps the tensor.

        Args:
            algebra (CliffordAlgebra): The algebra instance.
            tensor (torch.Tensor): Coefficients.
        """
        self.algebra = algebra
        self.tensor = tensor

    @classmethod
    def from_vectors(cls, algebra: CliffordAlgebra, vectors: torch.Tensor) -> "Multivector":
        """Promotes vectors to multivectors (Grade 1).

        Args:
            algebra (CliffordAlgebra): The algebra instance.
            vectors (torch.Tensor): Vectors [Batch, n].

        Returns:
            Multivector: The wrapper.
        """
        # Map vectors to multivector coefficients
        batch_shape = vectors.shape[:-1]
        mv_tensor = torch.zeros(*batch_shape, algebra.dim, device=vectors.device, dtype=vectors.dtype)

        for i in range(algebra.n):
            mv_tensor[..., 1 << i] = vectors[..., i]

        return cls(algebra, mv_tensor)

    def __repr__(self):
        return f"Multivector(shape={self.tensor.shape}, algebra=Cl({self.algebra.p},{self.algebra.q}))"

    def __add__(self, other):
        """Compute the multivector sum."""
        if isinstance(other, Multivector):
            assert self.algebra.n == other.algebra.n, "Algebras must match"
            return Multivector(self.algebra, self.tensor + other.tensor)
        elif isinstance(other, (int, float, torch.Tensor)):
            # Broadcast add
            return Multivector(self.algebra, self.tensor + other)
        else:
            return NotImplemented

    def __sub__(self, other):
        """Compute the multivector difference."""
        if isinstance(other, Multivector):
            return Multivector(self.algebra, self.tensor - other.tensor)
        else:
            return Multivector(self.algebra, self.tensor - other)

    def __mul__(self, other):
        """Compute the geometric product A * B."""
        if isinstance(other, Multivector):
            res = self.algebra.geometric_product(self.tensor, other.tensor)
            return Multivector(self.algebra, res)
        elif isinstance(other, (int, float)):
            return Multivector(self.algebra, self.tensor * other)
        else:
            return NotImplemented

    def __invert__(self):
        """Compute the reversion ~A."""
        return Multivector(self.algebra, self.algebra.reverse(self.tensor))

    def norm(self):
        """Compute the induced metric norm."""
        from core.metric import induced_norm
        return induced_norm(self.algebra, self.tensor)

    def exp(self):
        """Exponentiate via the algebra exp map."""
        return Multivector(self.algebra, self.algebra.exp(self.tensor))

    def grade(self, k: int):
        """Extract the grade-k component."""
        return Multivector(self.algebra, self.algebra.grade_projection(self.tensor, k))

__init__(algebra, tensor)

Wraps the tensor.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
tensor	Tensor	Coefficients.	required

Source code in core/multivector.py

def __init__(self, algebra: CliffordAlgebra, tensor: torch.Tensor):
    """Wraps the tensor.

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        tensor (torch.Tensor): Coefficients.
    """
    self.algebra = algebra
    self.tensor = tensor

from_vectors(algebra, vectors) classmethod

Promotes vectors to multivectors (Grade 1).

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
vectors	Tensor	Vectors [Batch, n].	required

Returns:

Name	Type	Description
Multivector	Multivector	The wrapper.

Source code in core/multivector.py

@classmethod
def from_vectors(cls, algebra: CliffordAlgebra, vectors: torch.Tensor) -> "Multivector":
    """Promotes vectors to multivectors (Grade 1).

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        vectors (torch.Tensor): Vectors [Batch, n].

    Returns:
        Multivector: The wrapper.
    """
    # Map vectors to multivector coefficients
    batch_shape = vectors.shape[:-1]
    mv_tensor = torch.zeros(*batch_shape, algebra.dim, device=vectors.device, dtype=vectors.dtype)

    for i in range(algebra.n):
        mv_tensor[..., 1 << i] = vectors[..., i]

    return cls(algebra, mv_tensor)

norm()

Compute the induced metric norm.

Source code in core/multivector.py

def norm(self):
    """Compute the induced metric norm."""
    from core.metric import induced_norm
    return induced_norm(self.algebra, self.tensor)

exp()

Exponentiate via the algebra exp map.

Source code in core/multivector.py

def exp(self):
    """Exponentiate via the algebra exp map."""
    return Multivector(self.algebra, self.algebra.exp(self.tensor))

grade(k)

Extract the grade-k component.

Source code in core/multivector.py

def grade(self, k: int):
    """Extract the grade-k component."""
    return Multivector(self.algebra, self.algebra.grade_projection(self.tensor, k))

Metric

metric

Metric definitions for Clifford algebras.

Provides distances, norms, and inner products that respect the metric signature.

inner_product(algebra, A, B)

Compute the scalar product via projection onto grade 0.

Computes _0.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	First multivector [Batch, Dim].	required
B	Tensor	Second multivector [Batch, Dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Scalar part [Batch, 1].

Source code in core/metric.py

def inner_product(algebra: CliffordAlgebra, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Compute the scalar product via projection onto grade 0.

    Computes <A B>_0.

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        A (torch.Tensor): First multivector [Batch, Dim].
        B (torch.Tensor): Second multivector [Batch, Dim].

    Returns:
        torch.Tensor: Scalar part [Batch, 1].
    """
    # Optimized: A * B then extract grade 0
    prod = algebra.geometric_product(A, B)
    scalar_part = prod[..., 0:1] # Grade 0
    return scalar_part

induced_norm(algebra, A)

Compute the induced norm respecting the metric signature.

Computes ||A|| = sqrt(|_0|).

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector [Batch, Dim].	required

Returns:

Type	Description
Tensor	torch.Tensor: Norm [Batch, 1].

Source code in core/metric.py

def induced_norm(algebra: CliffordAlgebra, A: torch.Tensor) -> torch.Tensor:
    """Compute the induced norm respecting the metric signature.

    Computes ||A|| = sqrt(|<A ~A>_0|).

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        A (torch.Tensor): Multivector [Batch, Dim].

    Returns:
        torch.Tensor: Norm [Batch, 1].
    """
    A_rev = algebra.reverse(A)
    # Scalar product <A A~>_0
    sq_norm = inner_product(algebra, A, A_rev)

    # In mixed signatures, sq_norm can be negative.
    # We return sqrt(|sq_norm|)
    return torch.sqrt(torch.abs(sq_norm))

geometric_distance(algebra, A, B)

Computes geometric distance.

dist(A, B) = ||A - B||.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	First multivector.	required
B	Tensor	Second multivector.	required

Returns:

Type	Description
Tensor	torch.Tensor: Distance.

Source code in core/metric.py

def geometric_distance(algebra: CliffordAlgebra, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Computes geometric distance.

    dist(A, B) = ||A - B||.

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        A (torch.Tensor): First multivector.
        B (torch.Tensor): Second multivector.

    Returns:
        torch.Tensor: Distance.
    """
    diff = A - B
    return induced_norm(algebra, diff)

grade_purity(algebra, A, grade)

Checks the purity of the grade by examining coefficient energy.

Purity = ||_k||^2 / ||A||^2.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector [..., Dim].	required
grade	int	Target grade.	required

Returns:

Type	Description
Tensor	torch.Tensor: Purity score [0, 1].

Source code in core/metric.py

def grade_purity(algebra: CliffordAlgebra, A: torch.Tensor, grade: int) -> torch.Tensor:
    """Checks the purity of the grade by examining coefficient energy.

    Purity = ||<A>_k||^2 / ||A||^2.

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        A (torch.Tensor): Multivector [..., Dim].
        grade (int): Target grade.

    Returns:
        torch.Tensor: Purity score [0, 1].
    """
    # Project to grade
    A_k = algebra.grade_projection(A, grade)

    # Compute energies (using standard squared norm of coefficients for stability)
    energy_k = (A_k ** 2).sum(dim=-1)
    energy_total = (A ** 2).sum(dim=-1).clamp(min=1e-6)

    return energy_k / energy_total

mean_active_grade(algebra, A)

Average grade. Identifies the grade where the majority of the energy resides.

Mean Grade = Sum(k * ||_k||^2) / ||A||^2.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector.	required

Returns:

Type	Description
Tensor	torch.Tensor: Average grade index.

Source code in core/metric.py

def mean_active_grade(algebra: CliffordAlgebra, A: torch.Tensor) -> torch.Tensor:
    """Average grade. Identifies the grade where the majority of the energy resides.

    Mean Grade = Sum(k * ||<A>_k||^2) / ||A||^2.

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        A (torch.Tensor): Multivector.

    Returns:
        torch.Tensor: Average grade index.
    """
    energy_total = (A ** 2).sum(dim=-1).clamp(min=1e-6)
    weighted_sum = torch.zeros_like(energy_total)

    for k in range(algebra.n + 1):
        A_k = algebra.grade_projection(A, k)
        energy_k = (A_k ** 2).sum(dim=-1)
        weighted_sum += k * energy_k

    return weighted_sum / energy_total

clifford_conjugate(algebra, mv)

Clifford conjugation (bar involution).

Combines reversion with grade involution

A_bar_k = (-1)^k * (-1)^{k(k-1)/2} * A_k

This is the natural *-involution on Cl(p,q). Useful for algebraic computations (e.g., spinor norms, Lipschitz groups).

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
mv	Tensor	Multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Conjugated multivector [..., Dim].

Source code in core/metric.py

def clifford_conjugate(algebra: CliffordAlgebra, mv: torch.Tensor) -> torch.Tensor:
    """Clifford conjugation (bar involution).

    Combines reversion with grade involution:
        A_bar_k = (-1)^k * (-1)^{k(k-1)/2} * A_k

    This is the natural *-involution on Cl(p,q). Useful for
    algebraic computations (e.g., spinor norms, Lipschitz groups).

    Args:
        algebra: The algebra instance.
        mv: Multivector [..., Dim].

    Returns:
        Conjugated multivector [..., Dim].
    """
    result = mv.clone()
    for i in range(algebra.dim):
        k = bin(i).count('1')
        sign = ((-1) ** k) * ((-1) ** (k * (k - 1) // 2))
        if sign == -1:
            result[..., i] = -result[..., i]
    return result

hermitian_inner_product(algebra, A, B)

Hermitian inner product on Cl(p,q): _0.

_H = Sum_I (conj_sign_I * metric_sign_I) * a_I * b_I

Uses precomputed sign arrays so that the result equals the scalar part of the geometric product of the Clifford conjugate of A with B. For Euclidean algebras (q=0), all signs are +1 and this reduces to the simple coefficient inner product Sum a_I b_I.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	First multivector [..., Dim].	required
B	Tensor	Second multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Scalar inner product [..., 1].

Source code in core/metric.py

def hermitian_inner_product(algebra: CliffordAlgebra, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Hermitian inner product on Cl(p,q): <bar{A} B>_0.

    <A, B>_H = Sum_I (conj_sign_I * metric_sign_I) * a_I * b_I

    Uses precomputed sign arrays so that the result equals the scalar
    part of the geometric product of the Clifford conjugate of A with B.
    For Euclidean algebras (q=0), all signs are +1 and this reduces to
    the simple coefficient inner product Sum a_I b_I.

    Args:
        algebra: The algebra instance.
        A: First multivector [..., Dim].
        B: Second multivector [..., Dim].

    Returns:
        Scalar inner product [..., 1].
    """
    signs = _hermitian_signs(algebra).to(device=A.device, dtype=A.dtype)
    return (signs * A * B).sum(dim=-1, keepdim=True)

hermitian_norm(algebra, A)

Hermitian norm: ||A||_H = sqrt(|_H|).

Always real and non-negative for any signature. Uses abs() since the signed inner product can produce negative self-products in mixed-signature algebras.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Norm [..., 1]. Always >= 0.

Source code in core/metric.py

def hermitian_norm(algebra: CliffordAlgebra, A: torch.Tensor) -> torch.Tensor:
    """Hermitian norm: ||A||_H = sqrt(|<A, A>_H|).

    Always real and non-negative for any signature.
    Uses abs() since the signed inner product can produce negative
    self-products in mixed-signature algebras.

    Args:
        algebra: The algebra instance.
        A: Multivector [..., Dim].

    Returns:
        Norm [..., 1]. Always >= 0.
    """
    sq = hermitian_inner_product(algebra, A, A)
    return torch.sqrt(torch.abs(sq))

hermitian_distance(algebra, A, B)

Hermitian distance: d_H(A, B) = ||A - B||_H.

Positive-definite metric distance for any signature. Satisfies: non-negativity, symmetry, triangle inequality, identity.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	First multivector [..., Dim].	required
B	Tensor	Second multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Distance [..., 1]. Always >= 0.

Source code in core/metric.py

def hermitian_distance(algebra: CliffordAlgebra, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Hermitian distance: d_H(A, B) = ||A - B||_H.

    Positive-definite metric distance for any signature.
    Satisfies: non-negativity, symmetry, triangle inequality, identity.

    Args:
        algebra: The algebra instance.
        A: First multivector [..., Dim].
        B: Second multivector [..., Dim].

    Returns:
        Distance [..., 1]. Always >= 0.
    """
    return hermitian_norm(algebra, A - B)

hermitian_angle(algebra, A, B)

Hermitian angle between multivectors.

cos(theta) = _H / (||A||_H * ||B||_H)

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	First multivector [..., Dim].	required
B	Tensor	Second multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Angle in radians [..., 1].

Source code in core/metric.py

def hermitian_angle(algebra: CliffordAlgebra, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Hermitian angle between multivectors.

    cos(theta) = <A, B>_H / (||A||_H * ||B||_H)

    Args:
        algebra: The algebra instance.
        A: First multivector [..., Dim].
        B: Second multivector [..., Dim].

    Returns:
        Angle in radians [..., 1].
    """
    signs = _hermitian_signs(algebra).to(device=A.device, dtype=A.dtype)
    ip = (signs * A * B).sum(dim=-1, keepdim=True)
    sq_a = (signs * A * A).sum(dim=-1, keepdim=True)
    sq_b = (signs * B * B).sum(dim=-1, keepdim=True)
    # Use sqrt(sq_a * sq_b) instead of sqrt(sq_a)*sqrt(sq_b) to avoid
    # float32 precision loss from two separate sqrt operations.
    denom = torch.sqrt(torch.abs(sq_a) * torch.abs(sq_b)).clamp(min=1e-6)
    cos_theta = ip / denom
    cos_theta = torch.clamp(cos_theta, -1.0, 1.0)
    return torch.acos(cos_theta)

grade_hermitian_norm(algebra, A, grade)

Hermitian norm restricted to a single grade.

||k||_H = sqrt(Sum{I: |I|=k} a_I**2)

Measures the energy contribution of a specific grade in a signature-independent way.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector [..., Dim].	required
grade	int	Target grade.	required

Returns:

Type	Description
Tensor	Grade-specific norm [..., 1].

Source code in core/metric.py

def grade_hermitian_norm(algebra: CliffordAlgebra, A: torch.Tensor, grade: int) -> torch.Tensor:
    """Hermitian norm restricted to a single grade.

    ||<A>_k||_H = sqrt(Sum_{I: |I|=k} a_I**2)

    Measures the energy contribution of a specific grade
    in a signature-independent way.

    Args:
        algebra: The algebra instance.
        A: Multivector [..., Dim].
        grade: Target grade.

    Returns:
        Grade-specific norm [..., 1].
    """
    A_k = algebra.grade_projection(A, grade)
    return hermitian_norm(algebra, A_k)

hermitian_grade_spectrum(algebra, A)

Full Hermitian grade spectrum.

Returns |_H| for each grade k = 0, ..., n. Uses abs() to ensure non-negative values in mixed signatures.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Grade energies [..., n+1]. Each entry >= 0.

Source code in core/metric.py

def hermitian_grade_spectrum(algebra: CliffordAlgebra, A: torch.Tensor) -> torch.Tensor:
    """Full Hermitian grade spectrum.

    Returns |<A_k, A_k>_H| for each grade k = 0, ..., n.
    Uses abs() to ensure non-negative values in mixed signatures.

    Args:
        algebra: The algebra instance.
        A: Multivector [..., Dim].

    Returns:
        Grade energies [..., n+1]. Each entry >= 0.
    """
    signs = _hermitian_signs(algebra).to(device=A.device, dtype=A.dtype)
    spectrum = []
    for k in range(algebra.n + 1):
        A_k = algebra.grade_projection(A, k)
        sq = (signs * A_k * A_k).sum(dim=-1, keepdim=True)
        spectrum.append(torch.abs(sq))
    return torch.cat(spectrum, dim=-1)

signature_trace_form(algebra, A, B)

Signature-aware trace form: <~A B>_0.

The standard Clifford algebra scalar product. NOT positive-definite in mixed signatures. Use hermitian_inner_product for optimization.

This form is signature-aware and useful for: - Rotor normalization (R~R = 1) - Versor validation - Spinor norm computation

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	First multivector [..., Dim].	required
B	Tensor	Second multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Scalar trace form [..., 1]. Can be negative in mixed signatures.

Source code in core/metric.py

def signature_trace_form(algebra: CliffordAlgebra, A: torch.Tensor, B: torch.Tensor) -> torch.Tensor:
    """Signature-aware trace form: <~A B>_0.

    The standard Clifford algebra scalar product. NOT positive-definite
    in mixed signatures. Use hermitian_inner_product for optimization.

    This form is signature-aware and useful for:
    - Rotor normalization (R~R = 1)
    - Versor validation
    - Spinor norm computation

    Args:
        algebra: The algebra instance.
        A: First multivector [..., Dim].
        B: Second multivector [..., Dim].

    Returns:
        Scalar trace form [..., 1]. Can be negative in mixed signatures.
    """
    A_rev = algebra.reverse(A)
    prod = algebra.geometric_product(A_rev, B)
    return prod[..., 0:1]

signature_norm_squared(algebra, A)

Signature-aware squared norm: _0.

Can be negative in mixed-signature algebras. Returns the raw value without absolute value, preserving causal structure information.

For Cl(n,0): always non-negative. For Cl(p,q) with q>0: sign encodes causal character.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
A	Tensor	Multivector [..., Dim].	required

Returns:

Type	Description
Tensor	Signed squared norm [..., 1].

Source code in core/metric.py

def signature_norm_squared(algebra: CliffordAlgebra, A: torch.Tensor) -> torch.Tensor:
    """Signature-aware squared norm: <A~A>_0.

    Can be negative in mixed-signature algebras. Returns the raw value
    without absolute value, preserving causal structure information.

    For Cl(n,0): always non-negative.
    For Cl(p,q) with q>0: sign encodes causal character.

    Args:
        algebra: The algebra instance.
        A: Multivector [..., Dim].

    Returns:
        Signed squared norm [..., 1].
    """
    return signature_trace_form(algebra, A, A)

Decomposition

decomposition

Bivector decomposition via GA power iteration.

Decomposes a general bivector into simple (blade) components that can each be exponentiated with the closed-form formula.

Reference

Pence, T., Yamada, D., & Singh, V. (2025). "Composing Linear Layers from Irreducibles." arXiv:2507.11688v1 [cs.LG]

ga_power_iteration(algebra, b, v_init=None, threshold=1e-06, max_iterations=100)

Find the dominant simple bivector component via power iteration.

Implements Algorithm 2 from Pence et al. (2025). Iterates v <- (b _| v) / ||b _| v|| until convergence, then recovers the simple projection b_s = sigma * (u ^ v).

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	CliffordAlgebra instance.	required
b	Tensor	Bivector to decompose [..., dim].	required
v_init	Optional[Tensor]	Initial grade-1 vector (random if None).	None
threshold	float	Convergence tolerance on ||v - v_prev||.	1e-06
max_iterations	int	Iteration cap.	100

Returns:

Type	Description
Tensor	(b_s, v) where b_s is the simple projection and v the converged
Tensor	vector, both shaped [..., dim].

Source code in core/decomposition.py

def ga_power_iteration(
    algebra,
    b: torch.Tensor,
    v_init: Optional[torch.Tensor] = None,
    threshold: float = 1e-6,
    max_iterations: int = 100
) -> Tuple[torch.Tensor, torch.Tensor]:
    """Find the dominant simple bivector component via power iteration.

    Implements Algorithm 2 from Pence et al. (2025).  Iterates
    ``v <- (b _| v) / ||b _| v||`` until convergence, then recovers
    the simple projection ``b_s = sigma * (u ^ v)``.

    Args:
        algebra (CliffordAlgebra): CliffordAlgebra instance.
        b: Bivector to decompose [..., dim].
        v_init: Initial grade-1 vector (random if None).
        threshold: Convergence tolerance on ``||v - v_prev||``.
        max_iterations: Iteration cap.

    Returns:
        (b_s, v) where b_s is the simple projection and v the converged
        vector, both shaped [..., dim].
    """
    batch_shape = b.shape[:-1]
    device = b.device
    dtype = b.dtype

    if v_init is None:
        v_raw = torch.randn(*batch_shape, algebra.n, device=device, dtype=dtype)
        v = algebra.embed_vector(v_raw)
    else:
        v = v_init

    v_norm = v.norm(dim=-1, keepdim=True)
    v = v / v_norm.clamp(min=1e-6)

    for _ in range(max_iterations):
        v_prev = v
        v = algebra.right_contraction(b, v)
        v_norm = v.norm(dim=-1, keepdim=True)
        v = v / v_norm.clamp(min=1e-6)

        if (v - v_prev).norm(dim=-1).max() < threshold:
            break

    u = algebra.right_contraction(b, v)
    u_norm = u.norm(dim=-1, keepdim=True)
    u = u / u_norm.clamp(min=1e-6)

    sigma = b.norm(dim=-1, keepdim=True)
    b_s = sigma * algebra.wedge(u, v)

    return b_s, v

differentiable_invariant_decomposition(algebra, b, k=None, threshold=1e-06, max_iterations=100)

Decompose a bivector into simple components via greedy projection.

Implements Algorithm 1 from Pence et al. (2025). Iteratively extracts the dominant simple component and subtracts it from the residual.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	CliffordAlgebra instance.	required
b	Tensor	Bivector [..., dim].	required
k	Optional[int]	Number of components (auto = n(n-1)/2 if None).	None
threshold	float	Stop when residual norm falls below this.	1e-06
max_iterations	int	Per-component power iteration cap.	100

Returns:

Type	Description
(decomp, vectors)	lists of simple bivectors and their
List[Tensor]	associated vectors.

Source code in core/decomposition.py

def differentiable_invariant_decomposition(
    algebra,
    b: torch.Tensor,
    k: Optional[int] = None,
    threshold: float = 1e-6,
    max_iterations: int = 100
) -> Tuple[List[torch.Tensor], List[torch.Tensor]]:
    """Decompose a bivector into simple components via greedy projection.

    Implements Algorithm 1 from Pence et al. (2025).  Iteratively
    extracts the dominant simple component and subtracts it from the
    residual.

    Args:
        algebra (CliffordAlgebra): CliffordAlgebra instance.
        b: Bivector [..., dim].
        k: Number of components (auto = n(n-1)/2 if None).
        threshold: Stop when residual norm falls below this.
        max_iterations: Per-component power iteration cap.

    Returns:
        (decomp, vectors): lists of simple bivectors and their
        associated vectors.
    """
    n = algebra.n
    k_max = (n * (n - 1)) // 2
    k = min(k, k_max) if k is not None else k_max

    decomp: List[torch.Tensor] = []
    vectors: List[torch.Tensor] = []
    residual = b.clone()

    for _ in range(k):
        if residual.norm(dim=-1).max() < threshold:
            break

        b_i, v_i = ga_power_iteration(
            algebra, residual, threshold=threshold, max_iterations=max_iterations
        )
        decomp.append(b_i)
        vectors.append(v_i)
        residual = residual - b_i

    return decomp, vectors

exp_simple_bivector(algebra, b)

Closed-form exponential of a simple bivector.

Delegates to algebra._exp_bivector_closed which handles all three signature regimes (elliptic, hyperbolic, parabolic).

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	CliffordAlgebra instance.	required
b	Tensor	Simple bivector [..., dim].	required

Returns:

Type	Description
Tensor	Rotor exp(b) [..., dim].

Source code in core/decomposition.py

def exp_simple_bivector(algebra, b: torch.Tensor) -> torch.Tensor:
    """Closed-form exponential of a simple bivector.

    Delegates to ``algebra._exp_bivector_closed`` which handles all
    three signature regimes (elliptic, hyperbolic, parabolic).

    Args:
        algebra (CliffordAlgebra): CliffordAlgebra instance.
        b: Simple bivector [..., dim].

    Returns:
        Rotor exp(b) [..., dim].
    """
    return algebra._exp_bivector_closed(b)

exp_decomposed(algebra, b, use_decomposition=True, k=None, threshold=1e-06, max_iterations=100)

Exponentiate a bivector via decomposition into simple components.

When use_decomposition is True the bivector is decomposed into simple blades (via differentiable_invariant_decomposition), each is exponentiated in closed form, and the rotors are composed via geometric product.

During training the power iteration loop is detached (run in forward-only mode) so that gradients do not flow through the normalization divisions. Gradients instead flow through the closed-form exp of each component and the final GP composition. This is stable for all bivector magnitudes.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	CliffordAlgebra instance.	required
b	Tensor	Bivector [..., dim].	required
use_decomposition	bool	Enable decomposition (False -> algebra.exp).	True
k	Optional[int]	Number of simple components (auto if None).	None
threshold	float	Convergence threshold.	1e-06
max_iterations	int	Power iteration cap.	100

Returns:

Type	Description
Tensor	Rotor exp(b) [..., dim].

Source code in core/decomposition.py

def exp_decomposed(
    algebra,
    b: torch.Tensor,
    use_decomposition: bool = True,
    k: Optional[int] = None,
    threshold: float = 1e-6,
    max_iterations: int = 100
) -> torch.Tensor:
    """Exponentiate a bivector via decomposition into simple components.

    When ``use_decomposition`` is True the bivector is decomposed into
    simple blades (via ``differentiable_invariant_decomposition``), each
    is exponentiated in closed form, and the rotors are composed via
    geometric product.

    During training the power iteration loop is **detached** (run in
    forward-only mode) so that gradients do not flow through the
    normalization divisions.  Gradients instead flow through the
    closed-form exp of each component and the final GP composition.
    This is stable for all bivector magnitudes.

    Args:
        algebra (CliffordAlgebra): CliffordAlgebra instance.
        b: Bivector [..., dim].
        use_decomposition: Enable decomposition (False -> ``algebra.exp``).
        k: Number of simple components (auto if None).
        threshold: Convergence threshold.
        max_iterations: Power iteration cap.

    Returns:
        Rotor exp(b) [..., dim].
    """
    if not use_decomposition:
        return algebra.exp(b)

    # Detach for decomposition (power iteration is not differentiable)
    # then re-project the original bivector onto the discovered planes.
    with torch.no_grad():
        decomp, _ = differentiable_invariant_decomposition(
            algebra, b.detach(), k=k, threshold=threshold,
            max_iterations=max_iterations
        )

    if len(decomp) == 0:
        result = torch.zeros_like(b)
        result[..., 0] = 1.0
        return result

    bv_mask = algebra.grade_masks[2]
    if bv_mask.device != b.device:
        bv_mask = bv_mask.to(b.device)

    rotors = []
    residual = b
    for b_i_detached in decomp:
        # Plane direction (unit simple bivector) -- detached
        plane_norm = b_i_detached.norm(dim=-1, keepdim=True).clamp(min=1e-12)
        plane_dir = b_i_detached / plane_norm  # detached unit plane

        # Project the live bivector onto this plane
        bv_live = residual[..., bv_mask]
        plane_bv = plane_dir[..., bv_mask]
        coeff = (bv_live * plane_bv).sum(dim=-1, keepdim=True)

        b_i_live = coeff * plane_dir
        residual = residual - b_i_live

        rotors.append(exp_simple_bivector(algebra, b_i_live))

    result = rotors[0]
    for R_i in rotors[1:]:
        result = algebra.geometric_product(result, R_i)

    return result

Signature Search

MetricSearch

Learns optimal (p, q, r) signature via GBN probe training and bivector energy analysis.

Trains small single-rotor GBN probes on conformally-lifted data using coherence + curvature as the loss. After training, reads the learned bivector energy distribution to infer the optimal signature.

Multiple probes with biased initialization combat local minima.

Source code in core/search.py

class MetricSearch:
    """Learns optimal (p, q, r) signature via GBN probe training and bivector
    energy analysis.

    Trains small single-rotor GBN probes on conformally-lifted data using
    coherence + curvature as the loss. After training, reads the learned
    bivector energy distribution to infer the optimal signature.

    Multiple probes with biased initialization combat local minima.
    """

    def __init__(
        self,
        device: str = 'cpu',
        num_probes: int = 6,
        probe_epochs: int = 80,
        probe_lr: float = 0.005,
        probe_channels: int = 4,
        k: int = 8,
        energy_threshold: float = 0.05,
        curvature_weight: float = 0.3,
        sparsity_weight: float = 0.01,
        max_workers: Optional[int] = None,
        micro_batch_size: Optional[int] = None,
        early_stop_patience: int = 0,
    ):
        """Initialize Metric Search.

        Args:
            device (str): Computation device.
            num_probes (int): Number of parallel probes to train.
            probe_epochs (int): Epochs per probe training.
            probe_lr (float): Learning rate for probes.
            probe_channels (int): Channels in probe models.
            k (int): Nearest neighbors for flow analysis.
            energy_threshold (float): Energy cutoff for active base vectors.
            curvature_weight (float): Curvature penalty in probe loss.
            sparsity_weight (float): Sparsity penalty in probe loss.
            max_workers (int, optional): Max threads for parallel training.
            micro_batch_size (int, optional): If set, train probes on random
                micro-batches of this size each epoch instead of full data.
            early_stop_patience (int): Stop probe training if best loss has
                not improved for this many epochs. 0 disables early stopping.
        """
        self.device = device
        self.num_probes = num_probes
        self.probe_epochs = probe_epochs
        self.probe_lr = probe_lr
        self.probe_channels = probe_channels
        self.k = k
        self.energy_threshold = energy_threshold
        self.curvature_weight = curvature_weight
        self.sparsity_weight = sparsity_weight
        self.max_workers = max_workers
        self.micro_batch_size = micro_batch_size
        self.early_stop_patience = early_stop_patience

    def _lift_data(self, data: torch.Tensor) -> Tuple[torch.Tensor, CliffordAlgebra]:
        """Lifts [N, X] data into Cl(X+1, 1, 0) via CGA-style embedding.

        Adds 2 extra dimensions (one positive, one negative):
            lifted = [x_1..x_X, 0.5*||x||^2, 1.0]

        Args:
            data (torch.Tensor): Input data [N, X].

        Returns:
            Tuple[torch.Tensor, CliffordAlgebra]: (mv_data [N, 1, dim], algebra).
        """
        data = data.to(self.device)
        N, X = data.shape

        if X + 2 > 8:
            warnings.warn(
                f"Data dimension {X} yields algebra dim 2^{X+2}={2**(X+2)}. "
                f"Consider PCA pre-reduction to X <= 6 for tractable computation."
            )

        # CGA-style: [x, 0.5*||x||^2, 1.0]
        norm_sq = 0.5 * (data ** 2).sum(dim=-1, keepdim=True)
        ones = torch.ones(N, 1, device=self.device, dtype=data.dtype)
        lifted = torch.cat([data, norm_sq, ones], dim=-1)  # [N, X+2]

        algebra = CliffordAlgebra(X + 1, 1, 0, device=self.device)
        mv = algebra.embed_vector(lifted)  # [N, 2^(X+2)]
        mv = mv.unsqueeze(1)  # [N, 1, dim]
        return mv, algebra

    def _train_probe(
        self,
        mv_data: torch.Tensor,
        algebra: CliffordAlgebra,
        bias_type: str = 'random',
    ) -> Dict:
        """Trains a single probe and returns results.

        Args:
            mv_data (torch.Tensor): Conformally lifted data [N, 1, dim].
            algebra (CliffordAlgebra): The conformal algebra.
            bias_type (str): Initialization bias type.

        Returns:
            Dict: Results with 'loss', 'coherence', 'curvature', 'probe'.
        """
        probe = _SignatureProbe(algebra, channels=self.probe_channels)
        probe.to(self.device)
        _apply_biased_init(probe, algebra, bias_type)

        gf = GeodesicFlow(algebra, k=self.k)
        optimizer = torch.optim.Adam(probe.parameters(), lr=self.probe_lr)

        best_loss = float('inf')
        best_state = None
        patience_counter = 0
        N = mv_data.shape[0]

        for _ in range(self.probe_epochs):
            # Micro-batch: sample a random subset each epoch
            if self.micro_batch_size and self.micro_batch_size < N:
                idx = torch.randperm(N, device=mv_data.device)[:self.micro_batch_size]
                batch = mv_data[idx]
            else:
                batch = mv_data

            optimizer.zero_grad()
            output = probe(batch)  # [B, 1, dim]
            output_flat = output.squeeze(1)  # [B, dim]

            coherence_t = gf._coherence_tensor(output_flat)
            curvature_t = gf._curvature_tensor(output_flat)

            sparsity = sum(
                r.sparsity_loss() for r in probe.get_rotor_layers()
            )

            loss = (
                -coherence_t
                + self.curvature_weight * curvature_t
                + self.sparsity_weight * sparsity
            )

            loss.backward()
            optimizer.step()

            loss_val = loss.item()
            if loss_val < best_loss:
                best_loss = loss_val
                best_state = copy.deepcopy(probe.state_dict())
                patience_counter = 0
            elif self.early_stop_patience > 0:
                patience_counter += 1
                if patience_counter >= self.early_stop_patience:
                    break

        # Restore best
        if best_state is not None:
            probe.load_state_dict(best_state)

        # Final metrics
        with torch.no_grad():
            output = probe(mv_data).squeeze(1)
            coh = gf.coherence(output)
            curv = gf.curvature(output)

        return {
            'loss': best_loss,
            'coherence': coh,
            'curvature': curv,
            'probe': probe,
        }

    def _analyze_bivector_energy(
        self,
        probe: _SignatureProbe,
        algebra: CliffordAlgebra,
        original_dim: int,
    ) -> Tuple[Tuple[int, int, int], Dict]:
        """Maps learned bivector energy to (p, q, r) signature.

        Algorithm:
        1. Collect bivector_weights from all RotorLayers -> energy = weight^2
        2. Average across channels
        3. For each basis bivector e_ab, look up bv_sq_scalar:
           -1 -> elliptic, +1 -> hyperbolic, 0 -> null
        4. Collect active base vector indices (energy > threshold) by type
        5. Subtract 2 conformal dimensions from counts -> (p, q, r)

        Args:
            probe (_SignatureProbe): Trained probe model.
            algebra (CliffordAlgebra): The conformal algebra.
            original_dim (int): Original data dimension X (before lifting).

        Returns:
            Tuple[Tuple[int, int, int], Dict]: ((p, q, r), energy_breakdown dict).
        """
        bv_sq = algebra.bv_sq_scalar  # [num_bivectors]
        bv_mask = algebra.grade_masks[2]
        bv_indices = bv_mask.nonzero(as_tuple=False).squeeze(-1)

        # Collect energy from all rotor layers
        total_energy = torch.zeros(len(bv_indices), device=self.device)
        n_layers = 0
        for rotor in probe.get_rotor_layers():
            with torch.no_grad():
                # [channels, num_bivectors] -> energy = weight^2, mean across channels
                energy = (rotor.bivector_weights ** 2).mean(dim=0)
                total_energy += energy
                n_layers += 1

        if n_layers > 0:
            total_energy /= n_layers

        # Normalize
        max_energy = total_energy.max().clamp(min=1e-8)
        normalized_energy = total_energy / max_energy

        # Classify base vectors by energy-weighted dominant type.
        n = algebra.n
        # base_type_energy[b] -> accumulated energy per signature type
        base_type_energy: dict = {}
        base_active: dict = {}

        for bv_idx_pos, blade_idx in enumerate(bv_indices.tolist()):
            energy_val = normalized_energy[bv_idx_pos].item()
            # Decode which base vectors form this bivector
            bits = []
            for bit in range(n):
                if blade_idx & (1 << bit):
                    bits.append(bit)
            if len(bits) != 2:
                continue

            sq_val = bv_sq[bv_idx_pos].item()
            if sq_val < -0.5:
                sig_type = 'elliptic'
            elif sq_val > 0.5:
                sig_type = 'hyperbolic'
            else:
                sig_type = 'null'

            for b in bits:
                if b not in base_type_energy:
                    base_type_energy[b] = {'elliptic': 0.0, 'hyperbolic': 0.0, 'null': 0.0}
                    base_active[b] = 0.0
                base_type_energy[b][sig_type] += energy_val
                base_active[b] = max(base_active[b], energy_val)

        # Count active base vectors by their dominant (highest cumulative energy) type
        active_positive = 0  # base vectors whose dominant bivectors are elliptic
        active_negative = 0  # base vectors whose dominant bivectors are hyperbolic
        active_null = 0      # base vectors whose dominant bivectors are null

        for b_idx in range(n):
            if b_idx not in base_active or base_active[b_idx] < self.energy_threshold:
                continue
            type_energy = base_type_energy[b_idx]
            dominant = max(type_energy, key=type_energy.get)
            if dominant == 'null':
                active_null += 1
            elif dominant == 'hyperbolic':
                active_negative += 1
            else:
                active_positive += 1

        # Subtract 2 conformal dimensions (1 positive + 1 negative from CGA lift)
        p = max(0, active_positive - 1)
        q = max(0, active_negative - 1)
        r = active_null

        # Clamp to original data dimension
        total = p + q + r
        if total > original_dim:
            # Scale down proportionally
            scale = original_dim / max(total, 1)
            p = max(1, round(p * scale))
            q = round(q * scale)
            r = round(r * scale)
            # Ensure they sum to <= original_dim
            while p + q + r > original_dim:
                if r > 0:
                    r -= 1
                elif q > 0:
                    q -= 1
                else:
                    p -= 1
        elif total == 0:
            p = original_dim  # Default to Euclidean

        energy_breakdown = {
            'per_bivector_energy': normalized_energy.tolist(),
            'active_positive': active_positive,
            'active_negative': active_negative,
            'active_null': active_null,
            'bv_sq_scalar': bv_sq.tolist(),
        }

        return (p, q, r), energy_breakdown

    def search(self, data: torch.Tensor) -> Tuple[int, int, int]:
        """Returns optimal (p, q, r) signature for the data.

        Args:
            data (torch.Tensor): Input data [N, D].

        Returns:
            Tuple[int, int, int]: Optimal signature (p, q, r).
        """
        result = self.search_detailed(data)
        return result['signature']

    def search_detailed(self, data: torch.Tensor) -> Dict:
        """Returns signature and full diagnostics.

        Args:
            data (torch.Tensor): Input data [N, D].

        Returns:
            Dict: Diagnostics with 'signature', 'coherence', 'curvature',
                'energy_breakdown', 'per_probe_results'.
        """
        data = data.to(self.device)
        N, X = data.shape

        # Conformal lift
        mv_data, algebra = self._lift_data(data)

        # Bias schedule: first 3 are deterministic, rest random
        bias_types = ['euclidean', 'minkowski', 'projective']
        while len(bias_types) < self.num_probes:
            bias_types.append('random')
        bias_types = bias_types[:self.num_probes]

        # Train probes
        def _run_probe(bias_type):
            return self._train_probe(mv_data, algebra, bias_type)

        if self.num_probes <= 2:
            # Sequential for small probe counts
            probe_results = [_run_probe(bt) for bt in bias_types]
        else:
            # Parallel with threads (avoids pickle issues, PyTorch releases GIL)
            max_w = self.max_workers or min(self.num_probes, 4)
            with concurrent.futures.ThreadPoolExecutor(max_workers=max_w) as pool:
                futures = [pool.submit(_run_probe, bt) for bt in bias_types]
                probe_results = [f.result() for f in futures]

        # Select best probe by loss
        best_idx = min(range(len(probe_results)), key=lambda i: probe_results[i]['loss'])
        best = probe_results[best_idx]

        # Analyze bivector energy
        signature, energy_breakdown = self._analyze_bivector_energy(
            best['probe'], algebra, X
        )

        return {
            'signature': signature,
            'coherence': best['coherence'],
            'curvature': best['curvature'],
            'energy_breakdown': energy_breakdown,
            'per_probe_results': [
                {'loss': r['loss'], 'coherence': r['coherence'], 'curvature': r['curvature']}
                for r in probe_results
            ],
        }

__init__(device='cpu', num_probes=6, probe_epochs=80, probe_lr=0.005, probe_channels=4, k=8, energy_threshold=0.05, curvature_weight=0.3, sparsity_weight=0.01, max_workers=None, micro_batch_size=None, early_stop_patience=0)

Initialize Metric Search.

Parameters:

Name	Type	Description	Default
device	str	Computation device.	'cpu'
num_probes	int	Number of parallel probes to train.	6
probe_epochs	int	Epochs per probe training.	80
probe_lr	float	Learning rate for probes.	0.005
probe_channels	int	Channels in probe models.	4
k	int	Nearest neighbors for flow analysis.	8
energy_threshold	float	Energy cutoff for active base vectors.	0.05
curvature_weight	float	Curvature penalty in probe loss.	0.3
sparsity_weight	float	Sparsity penalty in probe loss.	0.01
max_workers	int	Max threads for parallel training.	None
micro_batch_size	int	If set, train probes on random micro-batches of this size each epoch instead of full data.	None
early_stop_patience	int	Stop probe training if best loss has not improved for this many epochs. 0 disables early stopping.	0

Source code in core/search.py

def __init__(
    self,
    device: str = 'cpu',
    num_probes: int = 6,
    probe_epochs: int = 80,
    probe_lr: float = 0.005,
    probe_channels: int = 4,
    k: int = 8,
    energy_threshold: float = 0.05,
    curvature_weight: float = 0.3,
    sparsity_weight: float = 0.01,
    max_workers: Optional[int] = None,
    micro_batch_size: Optional[int] = None,
    early_stop_patience: int = 0,
):
    """Initialize Metric Search.

    Args:
        device (str): Computation device.
        num_probes (int): Number of parallel probes to train.
        probe_epochs (int): Epochs per probe training.
        probe_lr (float): Learning rate for probes.
        probe_channels (int): Channels in probe models.
        k (int): Nearest neighbors for flow analysis.
        energy_threshold (float): Energy cutoff for active base vectors.
        curvature_weight (float): Curvature penalty in probe loss.
        sparsity_weight (float): Sparsity penalty in probe loss.
        max_workers (int, optional): Max threads for parallel training.
        micro_batch_size (int, optional): If set, train probes on random
            micro-batches of this size each epoch instead of full data.
        early_stop_patience (int): Stop probe training if best loss has
            not improved for this many epochs. 0 disables early stopping.
    """
    self.device = device
    self.num_probes = num_probes
    self.probe_epochs = probe_epochs
    self.probe_lr = probe_lr
    self.probe_channels = probe_channels
    self.k = k
    self.energy_threshold = energy_threshold
    self.curvature_weight = curvature_weight
    self.sparsity_weight = sparsity_weight
    self.max_workers = max_workers
    self.micro_batch_size = micro_batch_size
    self.early_stop_patience = early_stop_patience

search(data)

Returns optimal (p, q, r) signature for the data.

Parameters:

Name	Type	Description	Default
data	Tensor	Input data [N, D].	required

Returns:

Type	Description
Tuple[int, int, int]	Tuple[int, int, int]: Optimal signature (p, q, r).

Source code in core/search.py

def search(self, data: torch.Tensor) -> Tuple[int, int, int]:
    """Returns optimal (p, q, r) signature for the data.

    Args:
        data (torch.Tensor): Input data [N, D].

    Returns:
        Tuple[int, int, int]: Optimal signature (p, q, r).
    """
    result = self.search_detailed(data)
    return result['signature']

search_detailed(data)

Returns signature and full diagnostics.

Parameters:

Name	Type	Description	Default
data	Tensor	Input data [N, D].	required

Returns:

Name	Type	Description
Dict	Dict	Diagnostics with 'signature', 'coherence', 'curvature', 'energy_breakdown', 'per_probe_results'.

Source code in core/search.py

def search_detailed(self, data: torch.Tensor) -> Dict:
    """Returns signature and full diagnostics.

    Args:
        data (torch.Tensor): Input data [N, D].

    Returns:
        Dict: Diagnostics with 'signature', 'coherence', 'curvature',
            'energy_breakdown', 'per_probe_results'.
    """
    data = data.to(self.device)
    N, X = data.shape

    # Conformal lift
    mv_data, algebra = self._lift_data(data)

    # Bias schedule: first 3 are deterministic, rest random
    bias_types = ['euclidean', 'minkowski', 'projective']
    while len(bias_types) < self.num_probes:
        bias_types.append('random')
    bias_types = bias_types[:self.num_probes]

    # Train probes
    def _run_probe(bias_type):
        return self._train_probe(mv_data, algebra, bias_type)

    if self.num_probes <= 2:
        # Sequential for small probe counts
        probe_results = [_run_probe(bt) for bt in bias_types]
    else:
        # Parallel with threads (avoids pickle issues, PyTorch releases GIL)
        max_w = self.max_workers or min(self.num_probes, 4)
        with concurrent.futures.ThreadPoolExecutor(max_workers=max_w) as pool:
            futures = [pool.submit(_run_probe, bt) for bt in bias_types]
            probe_results = [f.result() for f in futures]

    # Select best probe by loss
    best_idx = min(range(len(probe_results)), key=lambda i: probe_results[i]['loss'])
    best = probe_results[best_idx]

    # Analyze bivector energy
    signature, energy_breakdown = self._analyze_bivector_energy(
        best['probe'], algebra, X
    )

    return {
        'signature': signature,
        'coherence': best['coherence'],
        'curvature': best['curvature'],
        'energy_breakdown': energy_breakdown,
        'per_probe_results': [
            {'loss': r['loss'], 'coherence': r['coherence'], 'curvature': r['curvature']}
            for r in probe_results
        ],
    }

GeodesicFlow

Geodesic flow analysis in Clifford algebra.

Interprets data points as grade-1 multivectors and computes the flow field - a bivector at each point that encodes the direction of shortest algebraic paths to its k-nearest neighbours.

The flow is computed as the mean of connection bivectors:

B_ij = <x_i . ~x_j>_2   (grade-2 part of the geometric product)

This bivector encodes the rotational "turn" needed to map x_i toward x_j, analogous to the parallel transport connection on a Lie group.

The coherence and curvature of this field reveal whether the data has causal (directional) structure:

• High coherence, low curvature -> the flow is smooth and aligned in one direction. Causality is visible.
• Low coherence, high curvature -> the flow is fragmented and collides with itself. The signal is dominated by noise.

Source code in core/search.py

class GeodesicFlow:
    """Geodesic flow analysis in Clifford algebra.

    Interprets data points as grade-1 multivectors and computes the *flow
    field* - a bivector at each point that encodes the direction of shortest
    algebraic paths to its k-nearest neighbours.

    The flow is computed as the mean of *connection bivectors*:

        B_ij = <x_i . ~x_j>_2   (grade-2 part of the geometric product)

    This bivector encodes the rotational "turn" needed to map x_i toward
    x_j, analogous to the parallel transport connection on a Lie group.

    The coherence and curvature of this field reveal whether the data has
    causal (directional) structure:

    - **High coherence, low curvature** -> the flow is smooth and aligned in
      one direction.  Causality is visible.
    - **Low coherence, high curvature** -> the flow is fragmented and
      collides with itself.  The signal is dominated by noise.
    """

    def __init__(self, algebra: CliffordAlgebra, k: int = 8):
        """Initialize Geodesic Flow.

        Args:
            algebra (CliffordAlgebra): The algebra instance.
            k (int): Number of nearest neighbours for the flow field.
        """
        self.algebra = algebra
        self.k = k

    def _embed(self, data: torch.Tensor) -> torch.Tensor:
        """Embeds raw vectors into the grade-1 multivector subspace.

        Args:
            data (torch.Tensor): ``[N, d]`` where d <= algebra.n.

        Returns:
            torch.Tensor: ``[N, algebra.dim]`` grade-1 multivectors.
        """
        N, d = data.shape
        n = self.algebra.n
        if d > n:
            raise ValueError(
                f"Data dimension {d} exceeds algebra dimension {n}. "
                f"Use DimensionLifter to lift data before flow analysis."
            )
        # Pad to algebra dimension with zeros if d < n
        if d < n:
            pad = torch.zeros(N, n - d, device=data.device, dtype=data.dtype)
            data = torch.cat([data, pad], dim=-1)
        return self.algebra.embed_vector(data)

    def _knn(self, mv: torch.Tensor) -> torch.Tensor:
        """Returns k-nearest neighbour indices in multivector coefficient space.

        Args:
            mv (torch.Tensor): ``[N, dim]`` multivectors.

        Returns:
            torch.Tensor: ``[N, k]`` neighbour indices.
        """
        N = mv.shape[0]
        k = min(self.k, N - 1)
        diff = mv.unsqueeze(1) - mv.unsqueeze(0)     # [N, N, dim]
        dists = diff.norm(dim=-1)                    # [N, N]
        dists.fill_diagonal_(float('inf'))
        _, idx = dists.topk(k, dim=-1, largest=False)
        return idx                                    # [N, k]

    def _connection_bivectors(self, mv: torch.Tensor) -> torch.Tensor:
        """Computes unit connection bivectors for all (point, neighbour) pairs.

        The connection bivector from x_i to x_j encodes the rotational "turn"
        in the algebra needed to map one vector toward the other:

            B_ij = unit( <x_i . ~x_j>_2 )

        Args:
            mv (torch.Tensor): ``[N, dim]`` grade-1 multivectors.

        Returns:
            torch.Tensor: ``[N, k, dim]`` unit connection bivectors.
        """
        N, D = mv.shape
        k = min(self.k, N - 1)
        nn_idx = self._knn(mv)

        neighbors = mv[nn_idx]                          # [N, k, dim]
        xi = mv.unsqueeze(1).expand(N, k, D).reshape(N * k, D)
        xj_rev = self.algebra.reverse(neighbors.reshape(N * k, D))

        prod = self.algebra.geometric_product(xi, xj_rev)        # [N*k, dim]
        bv_raw = self.algebra.grade_projection(prod, 2)           # [N*k, dim]
        bv_norm = bv_raw.norm(dim=-1, keepdim=True).clamp(min=1e-8)
        return (bv_raw / bv_norm).reshape(N, k, D)               # [N, k, dim]

    def flow_bivectors(self, mv: torch.Tensor) -> torch.Tensor:
        """Computes the mean flow bivector at each data point.

        For each point x_i, aggregates the unit connection bivectors to its
        k-nearest neighbours:

            B_i = mean_j { unit( <x_i . ~x_j>_2 ) }

        .. note::
            For perfectly symmetric data (e.g. a closed circle) the mean
            cancels to zero - which is geometrically correct since there is
            no preferred flow direction.  Use :meth:`coherence` to measure
            structure without this cancellation.

        Args:
            mv (torch.Tensor): ``[N, dim]`` grade-1 multivectors.

        Returns:
            torch.Tensor: ``[N, dim]`` mean flow bivectors.
        """
        bv = self._connection_bivectors(mv)   # [N, k, dim]
        return bv.mean(dim=1)                 # [N, dim]

    def _coherence_tensor(self, mv: torch.Tensor) -> torch.Tensor:
        """Differentiable coherence -- returns a scalar tensor with grad_fn.

        Args:
            mv (torch.Tensor): ``[N, dim]`` multivectors.

        Returns:
            torch.Tensor: Scalar coherence in [0, 1].
        """
        bv = self._connection_bivectors(mv)   # [N, k, dim]
        N, k, D = bv.shape

        bi = bv.unsqueeze(2)                  # [N, k, 1, dim]
        bj = bv.unsqueeze(1)                  # [N, 1, k, dim]
        abs_cos = (bi * bj).sum(dim=-1).abs() # [N, k, k]

        mask = ~torch.eye(k, dtype=torch.bool, device=mv.device)  # [k, k]
        off_diag = abs_cos[:, mask]           # [N, k*(k-1)]
        return off_diag.mean()

    def coherence(self, mv: torch.Tensor) -> float:
        """Measures concentration of connection bivectors within each neighbourhood.

        For each point, computes the mean **absolute** cosine similarity between
        all pairs of its k connection bivectors.  This captures how consistently
        the neighbourhood connections lie on the same rotation plane.

        - **1.0**: all connections at every point are parallel or anti-parallel
          (maximally structured).
        - **1/num_bivectors** (~= baseline): connections point in random directions.

        .. note::
            In Cl(2,0) the grade-2 space is 1-dimensional (only e_12), so
            coherence is trivially 1.0 for any data - use at least Cl(3,0)
            for meaningful discrimination.

        Args:
            mv (torch.Tensor): ``[N, dim]`` multivectors.

        Returns:
            float: Coherence score in [0, 1].
        """
        return self._coherence_tensor(mv).item()

    def _curvature_tensor(self, mv: torch.Tensor) -> torch.Tensor:
        """Differentiable curvature -- returns a scalar tensor with grad_fn.

        Args:
            mv (torch.Tensor): ``[N, dim]`` multivectors.

        Returns:
            torch.Tensor: Scalar curvature in [0, 1].
        """
        bv = self._connection_bivectors(mv)   # [N, k, dim]
        N, k, D = bv.shape
        nn_idx = self._knn(mv)                # [N, k_nn]

        bi = bv.unsqueeze(2)                             # [N, k, 1, dim]
        bj_all = bv[nn_idx]                              # [N, k_nn, k, dim]

        bj = bj_all[:, 0]                                # [N, k, dim]
        bj = bj.unsqueeze(1)                             # [N, 1, k, dim]

        cross_cos = (bi * bj).sum(dim=-1).abs()          # [N, k, k]
        alignment = cross_cos.mean(dim=(-1, -2))         # [N]

        return 1.0 - alignment.mean()

    def curvature(self, mv: torch.Tensor) -> float:
        """Measures how much connection structure changes across the manifold.

        Computes the mean **dissimilarity** of connection bivectors between
        neighbouring pairs of points:

            dissimilarity(i, j) = 1 - mean_abs_cos( {B_ia}, {B_jb} )

        where {B_ia} is the set of k unit connection bivectors at point i and
        {B_jb} at point j, and mean_abs_cos is the cross-set absolute cosine
        similarity.

        - **0.0**: all neighbouring points have the same connection structure
          (flat geodesics, smooth manifold).
        - **High**: the connection direction changes rapidly between neighbours
          (high curvature, fragmented flow).

        Args:
            mv (torch.Tensor): ``[N, dim]`` multivectors.

        Returns:
            float: Curvature score in [0, 1].
        """
        return self._curvature_tensor(mv).item()

    def interpolate(
        self,
        a: torch.Tensor,
        b: torch.Tensor,
        steps: int = 10,
    ) -> torch.Tensor:
        """Interpolates along the geodesic from ``a`` to ``b``.

        Uses the *Lie group exponential map* on the transition element:

            T = a_inv . b
            log(T) ~= <T - 1>_2     (grade-2 approximation for small angles)
            gamma(t) = a . exp(t . log(T))

        Exact when a and b are close; a first-order approximation otherwise.

        Args:
            a (torch.Tensor): Start multivector ``[dim]``.
            b (torch.Tensor): End multivector ``[dim]``.
            steps (int): Number of interpolation steps (including endpoints).

        Returns:
            torch.Tensor: ``[steps, dim]`` sequence of multivectors.
        """
        a = a.unsqueeze(0)  # [1, dim]
        b = b.unsqueeze(0)  # [1, dim]

        # a_inv = ~a / <a . ~a>_0
        a_rev = self.algebra.reverse(a)
        a_sq = (a * a_rev)[..., 0:1].clamp(min=1e-8)   # grade-0 scalar
        a_inv = a_rev / a_sq                            # [1, dim]

        # Transition element T = a_inv . b
        T = self.algebra.geometric_product(a_inv, b)    # [1, dim]

        # Log approximation: grade-2 part of (T - 1)
        T_shift = T.clone()
        T_shift[..., 0] -= 1.0
        log_T = self.algebra.grade_projection(T_shift, 2)  # [1, dim]

        # Sample t in [0, 1]
        ts = torch.linspace(0.0, 1.0, steps, device=a.device, dtype=a.dtype)

        frames = []
        for t in ts:
            exp_tlogT = self.algebra.exp(t * log_T)       # [1, dim]
            frame = self.algebra.geometric_product(a, exp_tlogT)  # [1, dim]
            frames.append(frame)

        return torch.cat(frames, dim=0)  # [steps, dim]

    def causal_report(self, data: torch.Tensor) -> Dict:
        """Full geodesic flow analysis with a causal interpretation.

        Embeds data, computes coherence and curvature, and returns a
        human-readable verdict:

        - **Causal**: coherence > 0.5 and curvature < 0.5
          -> flow is smooth and aligned in one direction.
        - **Noisy**: otherwise
          -> flow is fragmented and collides with itself.

        Args:
            data (torch.Tensor): ``[N, d]`` raw data.

        Returns:
            Dict: report with keys ``coherence``, ``curvature``, ``causal``, ``label``.
        """
        mv = self._embed(data)
        coh = self.coherence(mv)
        curv = self.curvature(mv)
        is_causal = (coh > 0.5) and (curv < 0.5)
        return {
            'coherence': coh,
            'curvature': curv,
            'causal': is_causal,
            'label': (
                'Causal - smooth, aligned flow (low curvature)'
                if is_causal else
                'Noisy - fragmented, colliding flow (high curvature)'
            ),
        }

__init__(algebra, k=8)

Initialize Geodesic Flow.

Parameters:

Name	Type	Description	Default
algebra	CliffordAlgebra	The algebra instance.	required
k	int	Number of nearest neighbours for the flow field.	8

Source code in core/search.py

def __init__(self, algebra: CliffordAlgebra, k: int = 8):
    """Initialize Geodesic Flow.

    Args:
        algebra (CliffordAlgebra): The algebra instance.
        k (int): Number of nearest neighbours for the flow field.
    """
    self.algebra = algebra
    self.k = k

flow_bivectors(mv)

Computes the mean flow bivector at each data point.

For each point x_i, aggregates the unit connection bivectors to its k-nearest neighbours:

B_i = mean_j { unit( <x_i . ~x_j>_2 ) }

.. note:: For perfectly symmetric data (e.g. a closed circle) the mean cancels to zero - which is geometrically correct since there is no preferred flow direction. Use :meth:coherence to measure structure without this cancellation.

Parameters:

Name	Type	Description	Default
mv	Tensor	[N, dim] grade-1 multivectors.	required

Returns:

Type	Description
Tensor	torch.Tensor: [N, dim] mean flow bivectors.

Source code in core/search.py

def flow_bivectors(self, mv: torch.Tensor) -> torch.Tensor:
    """Computes the mean flow bivector at each data point.

    For each point x_i, aggregates the unit connection bivectors to its
    k-nearest neighbours:

        B_i = mean_j { unit( <x_i . ~x_j>_2 ) }

    .. note::
        For perfectly symmetric data (e.g. a closed circle) the mean
        cancels to zero - which is geometrically correct since there is
        no preferred flow direction.  Use :meth:`coherence` to measure
        structure without this cancellation.

    Args:
        mv (torch.Tensor): ``[N, dim]`` grade-1 multivectors.

    Returns:
        torch.Tensor: ``[N, dim]`` mean flow bivectors.
    """
    bv = self._connection_bivectors(mv)   # [N, k, dim]
    return bv.mean(dim=1)                 # [N, dim]

coherence(mv)

Measures concentration of connection bivectors within each neighbourhood.

For each point, computes the mean absolute cosine similarity between all pairs of its k connection bivectors. This captures how consistently the neighbourhood connections lie on the same rotation plane.

• 1.0: all connections at every point are parallel or anti-parallel (maximally structured).
• 1/num_bivectors (~= baseline): connections point in random directions.

.. note:: In Cl(2,0) the grade-2 space is 1-dimensional (only e_12), so coherence is trivially 1.0 for any data - use at least Cl(3,0) for meaningful discrimination.

Parameters:

Name	Type	Description	Default
mv	Tensor	[N, dim] multivectors.	required

Returns:

Name	Type	Description
float	float	Coherence score in [0, 1].

Source code in core/search.py

def coherence(self, mv: torch.Tensor) -> float:
    """Measures concentration of connection bivectors within each neighbourhood.

    For each point, computes the mean **absolute** cosine similarity between
    all pairs of its k connection bivectors.  This captures how consistently
    the neighbourhood connections lie on the same rotation plane.

    - **1.0**: all connections at every point are parallel or anti-parallel
      (maximally structured).
    - **1/num_bivectors** (~= baseline): connections point in random directions.

    .. note::
        In Cl(2,0) the grade-2 space is 1-dimensional (only e_12), so
        coherence is trivially 1.0 for any data - use at least Cl(3,0)
        for meaningful discrimination.

    Args:
        mv (torch.Tensor): ``[N, dim]`` multivectors.

    Returns:
        float: Coherence score in [0, 1].
    """
    return self._coherence_tensor(mv).item()

curvature(mv)

Measures how much connection structure changes across the manifold.

Computes the mean dissimilarity of connection bivectors between neighbouring pairs of points:

dissimilarity(i, j) = 1 - mean_abs_cos( {B_ia}, {B_jb} )

where {B_ia} is the set of k unit connection bivectors at point i and {B_jb} at point j, and mean_abs_cos is the cross-set absolute cosine similarity.

• 0.0: all neighbouring points have the same connection structure (flat geodesics, smooth manifold).
• High: the connection direction changes rapidly between neighbours (high curvature, fragmented flow).

Parameters:

Name	Type	Description	Default
mv	Tensor	[N, dim] multivectors.	required

Returns:

Name	Type	Description
float	float	Curvature score in [0, 1].

Source code in core/search.py

def curvature(self, mv: torch.Tensor) -> float:
    """Measures how much connection structure changes across the manifold.

    Computes the mean **dissimilarity** of connection bivectors between
    neighbouring pairs of points:

        dissimilarity(i, j) = 1 - mean_abs_cos( {B_ia}, {B_jb} )

    where {B_ia} is the set of k unit connection bivectors at point i and
    {B_jb} at point j, and mean_abs_cos is the cross-set absolute cosine
    similarity.

    - **0.0**: all neighbouring points have the same connection structure
      (flat geodesics, smooth manifold).
    - **High**: the connection direction changes rapidly between neighbours
      (high curvature, fragmented flow).

    Args:
        mv (torch.Tensor): ``[N, dim]`` multivectors.

    Returns:
        float: Curvature score in [0, 1].
    """
    return self._curvature_tensor(mv).item()

interpolate(a, b, steps=10)

Interpolates along the geodesic from a to b.

Uses the Lie group exponential map on the transition element:

T = a_inv . b
log(T) ~= <T - 1>_2     (grade-2 approximation for small angles)
gamma(t) = a . exp(t . log(T))

Exact when a and b are close; a first-order approximation otherwise.

Parameters:

Name	Type	Description	Default
a	Tensor	Start multivector [dim].	required
b	Tensor	End multivector [dim].	required
steps	int	Number of interpolation steps (including endpoints).	10

Returns:

Type	Description
Tensor	torch.Tensor: [steps, dim] sequence of multivectors.

Source code in core/search.py

def interpolate(
    self,
    a: torch.Tensor,
    b: torch.Tensor,
    steps: int = 10,
) -> torch.Tensor:
    """Interpolates along the geodesic from ``a`` to ``b``.

    Uses the *Lie group exponential map* on the transition element:

        T = a_inv . b
        log(T) ~= <T - 1>_2     (grade-2 approximation for small angles)
        gamma(t) = a . exp(t . log(T))

    Exact when a and b are close; a first-order approximation otherwise.

    Args:
        a (torch.Tensor): Start multivector ``[dim]``.
        b (torch.Tensor): End multivector ``[dim]``.
        steps (int): Number of interpolation steps (including endpoints).

    Returns:
        torch.Tensor: ``[steps, dim]`` sequence of multivectors.
    """
    a = a.unsqueeze(0)  # [1, dim]
    b = b.unsqueeze(0)  # [1, dim]

    # a_inv = ~a / <a . ~a>_0
    a_rev = self.algebra.reverse(a)
    a_sq = (a * a_rev)[..., 0:1].clamp(min=1e-8)   # grade-0 scalar
    a_inv = a_rev / a_sq                            # [1, dim]

    # Transition element T = a_inv . b
    T = self.algebra.geometric_product(a_inv, b)    # [1, dim]

    # Log approximation: grade-2 part of (T - 1)
    T_shift = T.clone()
    T_shift[..., 0] -= 1.0
    log_T = self.algebra.grade_projection(T_shift, 2)  # [1, dim]

    # Sample t in [0, 1]
    ts = torch.linspace(0.0, 1.0, steps, device=a.device, dtype=a.dtype)

    frames = []
    for t in ts:
        exp_tlogT = self.algebra.exp(t * log_T)       # [1, dim]
        frame = self.algebra.geometric_product(a, exp_tlogT)  # [1, dim]
        frames.append(frame)

    return torch.cat(frames, dim=0)  # [steps, dim]

causal_report(data)

Full geodesic flow analysis with a causal interpretation.

Embeds data, computes coherence and curvature, and returns a human-readable verdict:

• Causal: coherence > 0.5 and curvature < 0.5 -> flow is smooth and aligned in one direction.
• Noisy: otherwise -> flow is fragmented and collides with itself.

Parameters:

Name	Type	Description	Default
data	Tensor	[N, d] raw data.	required

Returns:

Name	Type	Description
Dict	Dict	report with keys coherence, curvature, causal, label.

Source code in core/search.py

def causal_report(self, data: torch.Tensor) -> Dict:
    """Full geodesic flow analysis with a causal interpretation.

    Embeds data, computes coherence and curvature, and returns a
    human-readable verdict:

    - **Causal**: coherence > 0.5 and curvature < 0.5
      -> flow is smooth and aligned in one direction.
    - **Noisy**: otherwise
      -> flow is fragmented and collides with itself.

    Args:
        data (torch.Tensor): ``[N, d]`` raw data.

    Returns:
        Dict: report with keys ``coherence``, ``curvature``, ``causal``, ``label``.
    """
    mv = self._embed(data)
    coh = self.coherence(mv)
    curv = self.curvature(mv)
    is_causal = (coh > 0.5) and (curv < 0.5)
    return {
        'coherence': coh,
        'curvature': curv,
        'causal': is_causal,
        'label': (
            'Causal - smooth, aligned flow (low curvature)'
            if is_causal else
            'Noisy - fragmented, colliding flow (high curvature)'
        ),
    }

DimensionLifter

Tests whether lifting data to a higher-dimensional algebra reveals structure.

The hypothesis: data living on an n-dimensional manifold may possess latent structure that only becomes visible in Cl(n+1, q) or Cl(n, q+1).

Lifting appends extra coordinates to the grade-1 embedding:

• Positive lift Cl(p, q) -> Cl(p+1, q): adds a spacelike dimension. The extra coordinate is set to 1 (projective / homogeneous lift).
• Null lift Cl(p, q) -> Cl(p, q+1): adds a timelike dimension. The extra coordinate is set to 0 (null vector lift for conformal-like embeddings).

After lifting, geodesic flow coherence is re-measured. An improvement indicates that the extra dimension captures hidden geometric structure.

Source code in core/search.py

class DimensionLifter:
    """Tests whether lifting data to a higher-dimensional algebra reveals structure.

    The hypothesis: data living on an n-dimensional manifold may possess
    latent structure that only becomes visible in Cl(n+1, q) or Cl(n, q+1).

    Lifting appends extra coordinates to the grade-1 embedding:

    - **Positive lift** ``Cl(p, q) -> Cl(p+1, q)``: adds a spacelike dimension.
      The extra coordinate is set to 1 (projective / homogeneous lift).
    - **Null lift** ``Cl(p, q) -> Cl(p, q+1)``: adds a timelike dimension.
      The extra coordinate is set to 0 (null vector lift for conformal-like embeddings).

    After lifting, geodesic flow coherence is re-measured.  An improvement
    indicates that the extra dimension captures hidden geometric structure.
    """

    def __init__(self, device: str = 'cpu'):
        """Initialize Dimension Lifter.

        Args:
            device (str): Computation device.
        """
        self.device = device

    def lift(
        self,
        data: torch.Tensor,
        target_algebra: CliffordAlgebra,
        fill: float = 1.0,
    ) -> torch.Tensor:
        """Lifts data into the grade-1 subspace of a higher-dimensional algebra.

        Pads each data vector with ``fill`` values in the new dimensions,
        then embeds as a grade-1 multivector.

        Args:
            data (torch.Tensor): ``[N, d]`` source data.
            target_algebra (CliffordAlgebra): Target algebra with n >= d.
            fill (float): Coordinate value for the new dimensions.
                Use 1.0 for a projective (homogeneous) lift,
                0.0 for a null-vector (conformal) lift.

        Returns:
            torch.Tensor: ``[N, target_algebra.dim]`` grade-1 multivectors.
        """
        N, d = data.shape
        n = target_algebra.n
        if n < d:
            raise ValueError(f"Target algebra dim {n} < source data dim {d}.")
        if n == d:
            return target_algebra.embed_vector(data.to(self.device))

        pad = torch.full(
            (N, n - d), fill,
            device=self.device, dtype=data.dtype
        )
        lifted = torch.cat([data.to(self.device), pad], dim=-1)  # [N, n]
        return target_algebra.embed_vector(lifted)

    def test(
        self,
        data: torch.Tensor,
        p: int,
        q: int,
        k: int = 8,
    ) -> Dict:
        """Compares geodesic flow coherence before and after dimension lifting.

        Tests three algebras in parallel:

        1. **Original** Cl(p, q): baseline coherence and curvature.
        2. **Positive lift** Cl(p+1, q): spacelike extra dimension, fill=1.
        3. **Null lift** Cl(p, q+1): timelike extra dimension, fill=0.

        Args:
            data (torch.Tensor): ``[N, d]`` data where d = p + q.
            p (int): Original positive signature.
            q (int): Original negative signature.
            k (int): Nearest neighbours for geodesic flow.

        Returns:
            Dict: results with keys ``original``, ``lift_positive``, ``lift_null``,
                and ``best`` (name of the algebra with highest coherence).
        """
        data = data.to(self.device)
        N, d = data.shape
        results: Dict = {}

        def _measure(alg: CliffordAlgebra, mv: torch.Tensor) -> Dict:
            gf = GeodesicFlow(alg, k=k)
            coh = gf.coherence(mv)
            curv = gf.curvature(mv)
            return {
                'signature': (alg.p, alg.q),
                'coherence': coh,
                'curvature': curv,
                'causal': (coh > 0.5) and (curv < 0.5),
            }

        # 1. Original
        alg_orig = CliffordAlgebra(p, q, device=self.device)
        mv_orig = alg_orig.embed_vector(data[..., :alg_orig.n])
        results['original'] = _measure(alg_orig, mv_orig)

        # 2. Positive lift: Cl(p+1, q)
        alg_pos = CliffordAlgebra(p + 1, q, device=self.device)
        mv_pos = self.lift(data, alg_pos, fill=1.0)
        results['lift_positive'] = _measure(alg_pos, mv_pos)

        # 3. Null lift: Cl(p, q+1)
        alg_null = CliffordAlgebra(p, q + 1, device=self.device)
        mv_null = self.lift(data, alg_null, fill=0.0)
        results['lift_null'] = _measure(alg_null, mv_null)

        # Which lifting gave the highest coherence?
        best = max(
            ('original', 'lift_positive', 'lift_null'),
            key=lambda key: results[key]['coherence'],
        )
        results['best'] = best

        return results

    def format_report(self, results: Dict) -> str:
        """Renders a lifting test result as a human-readable string.

        Args:
            results (Dict): Output of :meth:`test`.

        Returns:
            str: Multi-line report.
        """
        lines = ['Dimension Lifting Report', '=' * 40]
        for key in ('original', 'lift_positive', 'lift_null'):
            r = results[key]
            p, q = r['signature']
            coh = r['coherence']
            curv = r['curvature']
            causal = 'O Causal' if r['causal'] else 'X Noisy'
            lines.append(
                f"  Cl({p},{q})  coherence={coh:+.3f}  curvature={curv:.3f}  {causal}"
            )
        lines.append(f"\n  Best algebra: {results['best']}")
        return '\n'.join(lines)

__init__(device='cpu')

Initialize Dimension Lifter.

Parameters:

Name	Type	Description	Default
device	str	Computation device.	'cpu'

Source code in core/search.py

def __init__(self, device: str = 'cpu'):
    """Initialize Dimension Lifter.

    Args:
        device (str): Computation device.
    """
    self.device = device

lift(data, target_algebra, fill=1.0)

Lifts data into the grade-1 subspace of a higher-dimensional algebra.

Pads each data vector with fill values in the new dimensions, then embeds as a grade-1 multivector.

Parameters:

Name	Type	Description	Default
data	Tensor	[N, d] source data.	required
target_algebra	CliffordAlgebra	Target algebra with n >= d.	required
fill	float	Coordinate value for the new dimensions. Use 1.0 for a projective (homogeneous) lift, 0.0 for a null-vector (conformal) lift.	1.0

Returns:

Type	Description
Tensor	torch.Tensor: [N, target_algebra.dim] grade-1 multivectors.

Source code in core/search.py

def lift(
    self,
    data: torch.Tensor,
    target_algebra: CliffordAlgebra,
    fill: float = 1.0,
) -> torch.Tensor:
    """Lifts data into the grade-1 subspace of a higher-dimensional algebra.

    Pads each data vector with ``fill`` values in the new dimensions,
    then embeds as a grade-1 multivector.

    Args:
        data (torch.Tensor): ``[N, d]`` source data.
        target_algebra (CliffordAlgebra): Target algebra with n >= d.
        fill (float): Coordinate value for the new dimensions.
            Use 1.0 for a projective (homogeneous) lift,
            0.0 for a null-vector (conformal) lift.

    Returns:
        torch.Tensor: ``[N, target_algebra.dim]`` grade-1 multivectors.
    """
    N, d = data.shape
    n = target_algebra.n
    if n < d:
        raise ValueError(f"Target algebra dim {n} < source data dim {d}.")
    if n == d:
        return target_algebra.embed_vector(data.to(self.device))

    pad = torch.full(
        (N, n - d), fill,
        device=self.device, dtype=data.dtype
    )
    lifted = torch.cat([data.to(self.device), pad], dim=-1)  # [N, n]
    return target_algebra.embed_vector(lifted)

test(data, p, q, k=8)

Compares geodesic flow coherence before and after dimension lifting.

Tests three algebras in parallel:

1. Original Cl(p, q): baseline coherence and curvature.
2. Positive lift Cl(p+1, q): spacelike extra dimension, fill=1.
3. Null lift Cl(p, q+1): timelike extra dimension, fill=0.

Parameters:

Name	Type	Description	Default
data	Tensor	[N, d] data where d = p + q.	required
p	int	Original positive signature.	required
q	int	Original negative signature.	required
k	int	Nearest neighbours for geodesic flow.	8

Returns:

Name	Type	Description
Dict	Dict	results with keys original, lift_positive, lift_null, and best (name of the algebra with highest coherence).

Source code in core/search.py

def test(
    self,
    data: torch.Tensor,
    p: int,
    q: int,
    k: int = 8,
) -> Dict:
    """Compares geodesic flow coherence before and after dimension lifting.

    Tests three algebras in parallel:

    1. **Original** Cl(p, q): baseline coherence and curvature.
    2. **Positive lift** Cl(p+1, q): spacelike extra dimension, fill=1.
    3. **Null lift** Cl(p, q+1): timelike extra dimension, fill=0.

    Args:
        data (torch.Tensor): ``[N, d]`` data where d = p + q.
        p (int): Original positive signature.
        q (int): Original negative signature.
        k (int): Nearest neighbours for geodesic flow.

    Returns:
        Dict: results with keys ``original``, ``lift_positive``, ``lift_null``,
            and ``best`` (name of the algebra with highest coherence).
    """
    data = data.to(self.device)
    N, d = data.shape
    results: Dict = {}

    def _measure(alg: CliffordAlgebra, mv: torch.Tensor) -> Dict:
        gf = GeodesicFlow(alg, k=k)
        coh = gf.coherence(mv)
        curv = gf.curvature(mv)
        return {
            'signature': (alg.p, alg.q),
            'coherence': coh,
            'curvature': curv,
            'causal': (coh > 0.5) and (curv < 0.5),
        }

    # 1. Original
    alg_orig = CliffordAlgebra(p, q, device=self.device)
    mv_orig = alg_orig.embed_vector(data[..., :alg_orig.n])
    results['original'] = _measure(alg_orig, mv_orig)

    # 2. Positive lift: Cl(p+1, q)
    alg_pos = CliffordAlgebra(p + 1, q, device=self.device)
    mv_pos = self.lift(data, alg_pos, fill=1.0)
    results['lift_positive'] = _measure(alg_pos, mv_pos)

    # 3. Null lift: Cl(p, q+1)
    alg_null = CliffordAlgebra(p, q + 1, device=self.device)
    mv_null = self.lift(data, alg_null, fill=0.0)
    results['lift_null'] = _measure(alg_null, mv_null)

    # Which lifting gave the highest coherence?
    best = max(
        ('original', 'lift_positive', 'lift_null'),
        key=lambda key: results[key]['coherence'],
    )
    results['best'] = best

    return results

format_report(results)

Renders a lifting test result as a human-readable string.

Parameters:

Name	Type	Description	Default
results	Dict	Output of :meth:test.	required

Returns:

Name	Type	Description
str	str	Multi-line report.

Source code in core/search.py

def format_report(self, results: Dict) -> str:
    """Renders a lifting test result as a human-readable string.

    Args:
        results (Dict): Output of :meth:`test`.

    Returns:
        str: Multi-line report.
    """
    lines = ['Dimension Lifting Report', '=' * 40]
    for key in ('original', 'lift_positive', 'lift_null'):
        r = results[key]
        p, q = r['signature']
        coh = r['coherence']
        curv = r['curvature']
        causal = 'O Causal' if r['causal'] else 'X Noisy'
        lines.append(
            f"  Cl({p},{q})  coherence={coh:+.3f}  curvature={curv:.3f}  {causal}"
        )
    lines.append(f"\n  Best algebra: {results['best']}")
    return '\n'.join(lines)