Geometric Algebra¶
See Wikipedia for an idea of what this is.
New 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.
 Parameters
basis – A sequence of names of basis vectors, or an integer (the number of dimensions) to use the default names
e0
througheN
.metric_matrix – See
metric_matrix
. If None, the Euclidean metric is assumed.

dimensions
¶

is_orthogonal
¶

is_euclidean
¶

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 scalarish 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 ObjectOriented 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
sympy.galgebra.GA.MV
object insympy
, especially with respect to the supported operators:Operation
Result
A+B
Sum of multivectors
AB
Difference of multivectors
A*B
Geometric product \(AB\)
A^B
Outer product \(A\wedge B\) of multivectors
AB
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 oddgrade 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].

norm_squared
()¶

__abs__
()¶
Comparisons
MultiVector
objects have a truth value corresponding to whether they have any blades with nonzero 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
, returnself.project(0).as_scalar()
. Ifr == 1
, returnself.project(1).as_vector(dtype)
. Otherwise, returnself.project(r)
.

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

odd
()¶ Extract the oddgrade blades.

even
()¶ Extract the evengrade blades.

project_min_grade
()¶ New in version 2014.2.

project_max_grade
()¶ New in version 2014.2.

as_scalar
()¶

as_vector
(dtype=None)¶ Return a
numpy
vector corresponding to the grade1MultiVector
self.If self is not grade1,
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.
 Parameters
data –
This may be one of the following:
a
numpy.ndarray
, which will be turned into a grade1 multivector,a mapping from tuples of basis indices (together indicating a blade, order matters and will be mapped to ‘normalized’ blades) to coefficients,
an array as described in
data
,a scalar–where everything that doesn’t fall into the above cases is viewed as a scalar.
space – A
Space
instance. If None or an integer,get_euclidean_space()
is called to obtain a default space with the right number of dimensions for data. Note: dimension guessing only works when anumpy.ndarray
is being passed for data.

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(1.85*e0 + 2.988*e1 + 5.2016*e2)
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(1*e0)
>>> print one.I ** 2
MV(1.0)
>>> print (3+5j)*(2+3j)/(3j)
(6.33333333333+3j)
>>> print (3+5*one.I)*(2+3*one.I)/(3*one.I)
MV(6.33333333333 + 3.0*e0)
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) + 1e13
# contractions
# (3.18) in [DFM]
assert abs(b.scalar_product(a ^ c)  (b >> a).scalar_product(c)) < 1e13
# 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))