Geometric Algebra

See Wikipedia for an idea of what this is.

Added in version 2013.2.

Also see Example usage.

Spaces

class pymbolic.geometric_algebra.Space(basis=None, metric_matrix=None)

basis_names

A sequence of names of basis vectors.

metric_matrix

A (dims,dims)-shaped matrix, whose (i,j)-th entry represents the inner product of basis vector i and basis vector j.

dimensions

is_orthogonal

is_euclidean

pymbolic.geometric_algebra.get_euclidean_space(n)

Return the canonical n-dimensional Euclidean Space.

Multivectors

class pymbolic.geometric_algebra.MultiVector(data, space=None)

An immutable multivector type. Its implementation follows [DFM]. It is pickleable, and not picky about what data is used as coefficients. It supports pymbolic.primitives.Expression objects of course, but it can take just about any other scalar-ish coefficients.

data

A mapping from a basis vector bitmap (indicating blades) to coefficients. (see [DFM], Chapter 19 for idea and rationale)

space

See the following literature:

[DFM] L. Dorst, D. Fontijne, and S. Mann, Geometric Algebra for Computer Science: An Object-Oriented Approach to Geometry. Morgan Kaufmann, 2010.

[HS] D. Hestenes and G. Sobczyk, Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics. Springer, 1987.

The object behaves much like the corresponding galgebra.mv.Mv object in galgebra, especially with respect to the supported operators:

Operation	Result
A+B	Sum of multivectors
A-B	Difference of multivectors
A*B	Geometric product \(AB\)
A^B	Outer product \(A\wedge B\) of multivectors
A|B	Inner product \(A\cdot B\) of multivectors
A<<B	Left contraction \(A\lrcorner B\) (_|) of multivectors, also read as ‘\(A\) removed from \(B\)’.
A>>B	Right contraction \(A\llcorner B\) (|_)of multivectors, also read as ‘\(A\) without \(B\)’.

Warning

Many of the multiplicative operators bind more weakly than even addition. Python’s operator precedence further does not match geometric algebra, which customarily evaluates outer, inner, and then geometric.

In other words: Use parentheses everywhere.

mapper_method = 'map_multivector'

More products

scalar_product(other)

Return the scalar product, as a scalar, not a MultiVector.

Often written \(A*B\).

x(other)

Return the commutator product.

See (1.1.55) in [HS].

Often written \(A\times B\).

__pow__(other)

Return self to the integer power other.

Unary operators

inv()

Return the multiplicative inverse of the blade self.

Often written \(A^{-1}\).

rev()

Return the reverse of self, i.e. the multivector obtained by reversing the order of all component blades.

Often written \(A^\dagger\).

invol()

Return the grade involution (see Section 2.9.5 of [DFM]), i.e. all odd-grade blades have their signs flipped.

Often written \(\widehat A\).

dual()

Return the dual of self, see (1.2.26) in [HS].

Written \(\widetilde A\) by [HS] and \(A^\ast\) by [DFW].

__inv__()

Return the dual of self, see dual().

norm_squared()

__abs__()

I

Return the pseudoscalar associated with this object’s Space.

Comparisons

MultiVector objects have a truth value corresponding to whether they have any blades with non-zero coefficients. They support testing for (exact) equality.

zap_near_zeros(tol=None)

Remove blades whose coefficient is close to zero relative to the norm of self.

close_to(other, tol=None)

Grade manipulation

gen_blades(grade=None)

Generate all blades in self, optionally only those of a specific grade.

project(r)

Return a new multivector containing only the blades of grade r.

Often written \(\langle A\rangle_r\).

xproject(r, dtype=None)

If r == 0, return self.project(0).as_scalar(). If r == 1, return self.project(1).as_vector(dtype). Otherwise, return self.project(r).

all_grades()

Return a set of grades occurring in self.

get_pure_grade()

If self only has components of a single grade, return that as an integer. Otherwise, return None.

odd()

Extract the odd-grade blades.

even()

Extract the even-grade blades.

project_min_grade()

Added in version 2014.2.

project_max_grade()

Added in version 2014.2.

as_scalar()

as_vector(dtype=None)

Return a numpy vector corresponding to the grade-1 MultiVector self.

If self is not grade-1, ValueError is raised.

Helper functions

map(f)

Return a new MultiVector with coefficients mapped by function f, which takes a single coefficient as input and returns the new coefficient.

Example usage

This first example demonstrates how to compute a cross product using MultiVector:

>>> import numpy as np
>>> import pymbolic.geometric_algebra as ga
>>> MV = ga.MultiVector

>>> a = np.array([3.344, 1.2, -0.5])
>>> b = np.array([7.4, 1.1, -2.0])
>>> np.cross(a, b)
array([-1.85  ,  2.988 , -5.2016])

>>> mv_a = MV(a)
>>> mv_b = MV(b)
>>> print(-mv_a.I*(mv_a^mv_b))
MV(
    e0 * -1.8499999999999999
    + e1 * 2.9879999999999995
    + e2 * -5.201600000000001)

This simple example demonstrates how a complex number is simply a special case of a MultiVector:

>>> import numpy as np
>>> import pymbolic.geometric_algebra as ga
>>> MV = ga.MultiVector
>>>
>>> sp = ga.Space(metric_matrix=-np.eye(1))
>>> sp
Space(['e0'], array([[-1.]]))

>>> one = MV(1, sp)
>>> one
MultiVector({0: 1}, Space(['e0'], array([[-1.]])))
>>> print(one)
MV(1)
>>> print(one.I)
MV(e0 * 1)
>>> print(one.I ** 2)
MV(-1.0)

>>> print((3+5j)*(2+3j)/(3j))
(6.333333333333333+3j)
>>> print((3+5*one.I)*(2+3*one.I)/(3*one.I))
MV(
    6.333333333333333
    + e0 * 3.0)

The following test demonstrates the use of the object and shows many useful properties:

def test_geometric_algebra(dims):
    pytest.importorskip("numpy")

    import numpy as np
    from pymbolic.geometric_algebra import MultiVector as MV  # noqa

    vec1 = MV(np.random.randn(dims))
    vec2 = MV(np.random.randn(dims))
    vec3 = MV(np.random.randn(dims))
    vec4 = MV(np.random.randn(dims))
    vec5 = MV(np.random.randn(dims))

    # Fundamental identity
    assert ((vec1 ^ vec2) + (vec1 | vec2)).close_to(vec1*vec2)

    # Antisymmetry
    assert (vec1 ^ vec2 ^ vec3).close_to(- vec2 ^ vec1 ^ vec3)

    vecs = [vec1, vec2, vec3, vec4, vec5]

    if len(vecs) > dims:
        from operator import xor as outer
        assert reduce(outer, vecs).close_to(0)

    assert (vec1.inv()*vec1).close_to(1)
    assert (vec1*vec1.inv()).close_to(1)
    assert ((1/vec1)*vec1).close_to(1)
    assert (vec1/vec1).close_to(1)

    for a, b, c in [
            (vec1, vec2, vec3),
            (vec1*vec2, vec3, vec4),
            (vec1, vec2*vec3, vec4),
            (vec1, vec2, vec3*vec4),
            (vec1, vec2, vec3*vec4*vec5),
            (vec1, vec2*vec1, vec3*vec4*vec5),
            ]:

        # Associativity
        assert ((a*b)*c).close_to(a*(b*c))
        assert ((a ^ b) ^ c).close_to(a ^ (b ^ c))
        # The inner product is not associative.

        # scalar product
        assert ((c*b).project(0)) .close_to(b.scalar_product(c))
        assert ((c.rev()*b).project(0)) .close_to(b.rev().scalar_product(c))
        assert ((b.rev()*b).project(0)) .close_to(b.norm_squared())

        assert b.norm_squared() >= 0
        assert c.norm_squared() >= 0

        # Cauchy's inequality
        assert b.scalar_product(c) <= abs(b)*abs(c) + 1e-13

        # contractions

        # (3.18) in [DFM]
        assert abs(b.scalar_product(a ^ c) - (b >> a).scalar_product(c)) < 1e-13

        # duality, (3.20) in [DFM]
        assert ((a ^ b) << c) .close_to(a << (b << c))

        # two definitions of the dual agree: (1.2.26) in [HS]
        # and (sec 3.5.3) in [DFW]
        assert (c << c.I.rev()).close_to(c | c.I.rev())

        # inverse
        for div in list(b.gen_blades()) + [vec1, vec1.I]:
            assert (div.inv()*div).close_to(1)
            assert (div*div.inv()).close_to(1)
            assert ((1/div)*div).close_to(1)
            assert (div/div).close_to(1)
            assert ((c/div)*div).close_to(c)
            assert ((c*div)/div).close_to(c)

        # reverse properties (Sec 2.9.5 [DFM])
        assert c.rev().rev() == c
        assert (b ^ c).rev() .close_to(c.rev() ^ b.rev())

        # dual properties
        # (1.2.26) in [HS]
        assert c.dual() .close_to(c | c.I.rev())
        assert c.dual() .close_to(c*c.I.rev())

        # involution properties (Sec 2.9.5 DFW)
        assert c.invol().invol() == c
        assert (b ^ c).invol() .close_to(b.invol() ^ c.invol())

        # commutator properties

        # Jacobi identity (1.1.56c) in [HS] or (8.2) in [DFW]
        assert (a.x(b.x(c)) + b.x(c.x(a)) + c.x(a.x(b))).close_to(0)

        # (1.57) in [HS]
        assert a.x(b*c) .close_to(a.x(b)*c + b*a.x(c))