clifford: Geometric Algebra for Python

In [1]: from clifford.g3 import *  # import GA for 3D space

In [2]: import math

In [3]: a = e1 + 2*e2 + 3*e3  # vector

In [4]: R = math.e**(math.pi/4*e12)  # rotor

In [5]: R*a*~R    # rotate the vector
Out[5]: (2.0^e1) - (1.0^e2) + (3.0^e3)

This module implements Geometric Algebras (a.k.a. Clifford algebras). Geometric Algebra (GA) is a universal algebra which subsumes complex algebra, quaternions, linear algebra and several other independent mathematical systems. Scalars, vectors, and higher-grade entities can be mixed freely and consistently in the form of mixed-grade multivectors.

Installation

Clifford can be installed via the standard mechanisms for Python packages, and is available through both PyPI and Anaconda.

conda

If you are using the Anaconda python distribution, then you should install using the conda command.

Once you have anaconda installed you can install clifford from the conda-forge channel with the command:

conda install -c conda-forge clifford

PyPI

If you are not using Anaconda, you can install clifford from PyPI using the pip command that comes with Python:

pip3 install clifford

Note

The 1.3 release and older are not compatible with numba 0.50.0, and emit a large number of warnings on 0.49.0. To avoid these warnings, install with:

pip3 install clifford "numba~=0.48.0" "sparse~=0.9.0"

If you want the development version of clifford you can pip install from the latest master with:

pip3 install git+https://github.com/pygae/clifford.git@master

If you want to modify clifford itself, then you should use an editable (-e) installation:

git clone https://github.com/pygae/clifford.git
pip3 install -e ./clifford

If you are not running as sudo, and your python installation is system-wide, you will need to pass --user to the invocations above.

API

clifford (clifford)

The top-level module. Provides two core classes, Layout and MultiVector, along with several helper functions to implement the algebras.

Constructing algebras

Note that typically the Predefined Algebras are sufficient, and there is no need to build an algebra from scratch.

Cl([p, q, r, sig, names, firstIdx, mvClass])	Returns a Layout and basis blade MultiVectors for the geometric algebra \(Cl_{p,q,r}\).
conformalize(layout[, added_sig, mvClass])	Conformalize a Geometric Algebra

Whether you construct your algebras from scratch, or use the predefined ones, you’ll end up working with the following types:

MultiVector(layout[, value, string])	An element of the algebra
Layout(sig, *[, ids, order, names])	Layout stores information regarding the geometric algebra itself and the internal representation of multivectors.
ConformalLayout(*args[, layout])	A layout for a conformal algebra, which adds extra constants and helpers.

Advanced algebra configuration

It is unlikely you will need these features, but they remain as a better spelling for features which have always been in clifford.

BasisBladeOrder(bitmaps)	Represents the storage order in memory of basis blade coefficients.
BasisVectorIds(blade_ids)	Stores ids for the ordered set of basis vectors, typically integers.

Global configuration functions

These functions are used to change the global behavior of clifford.

clifford.eps(newEps=None)

Get/Set the epsilon for float comparisons.

clifford.pretty(precision=None)

Makes repr(MultiVector) default to pretty-print.

precision arg can be used to set the printed precision.

Parameters

precision (int) – number of sig figs to print past decimal

Examples

>>> pretty(5)

clifford.ugly()

Makes repr(MultiVector) default to eval-able representation.

clifford.print_precision(newVal)

Set the epsilon for float comparisons.

Parameters

newVal (int) – number of sig figs to print (see builtin round)

Examples

>>> print_precision(5)

Miscellaneous classes

MVArray(input_array)	MultiVector Array
Frame(input_array)	A frame of vectors
BladeMap(blades_map[, map_scalars])	A Map Relating Blades in two different algebras

Miscellaneous functions

grade_obj(objin[, threshold])	Returns the modal grade of a multivector
randomMV(layout[, min, max, grades, …])	n Random MultiVectors with given layout.

cga (clifford.cga)

Object Oriented Conformal Geometric Algebra.

Examples

>>> from clifford import Cl
>>> from clifford.cga import CGA
>>> g3, blades = Cl(3)
>>> locals().update(blades)
>>> g3c = CGA(g3)
>>> C = g3c.round(3)             # create random sphere
>>> T = g3c.translation(e1+e2)   # create translation
>>> C_ = T(C)                    # translate the sphere
>>> C_.center                    # compute center of sphere
-(1.0^e4) - (1.0^e5)

The CGA

CGA(layout_orig)

Conformal Geometric Algebra

Objects

Flat(cga, *args)	A line, plane, or hyperplane.
Round(cga, *args)	A point pair, circle, sphere or hyper-sphere.

Operators

Rotation(cga, *args)	A Rotation
Dilation(cga, *args)	A global dilation
Translation(cga, *args)	A Translation
Transversion(cga, *args)	A Transversion

Meta-Class

CGAThing(cga)

base class for cga objects and operators.

tools (clifford.tools)

Algorithms and tools of various kinds.

Tools for specific ga’s

g3	Tools for 3DGA (g3)
g3c	Tools for 3DCGA (g3c)

Classifying conformal GAs

classify

Tools for interpreting conformal blades

Determining Rotors From Frame Pairs or Orthogonal Matrices

Given two frames that are related by a orthogonal transform, we seek a rotor which enacts the transform. Details of the mathematics and pseudo-code used to create the algorithms below can be found at Allan Cortzen’s website, [Cor17].

There are also some helper functions which can be used to translate matrices into GA frames, so an orthogonal (or complex unitary) matrix can be directly translated into a Versor.

orthoFrames2Versor(B[, A, delta, eps, det, …])	Determines versor for two frames related by an orthogonal transform
orthoMat2Versor(A[, eps, layout, is_complex])	Translates an orthogonal (or unitary) matrix to a Versor
mat2Frame(A[, layout, is_complex])	Translates a (possibly complex) matrix into a real vector frame
omoh(A, B)	Determines homogenization scaling for two Frames related by a Rotor

operator functions (clifford.operator)

This module exists to enable functional programming via functools.reduce(). It can be thought of as equivalent to the builtin operator module, but for the operators from geometric algebra.

>>> import functools
>>> import clifford.operator
>>> from clifford.g3 import *
>>> Ms = [e1, e1 + e2, e2 + e3]  # list of multivectors
>>> assert functools.reduce(clifford.operator.op, Ms) == Ms[0] ^ Ms[1] ^ Ms[2]

clifford.operator.gp(M, N)

Geometric product function \(MN\), equivalent to M * N.

M and N must be from the same layout

clifford.operator.op(M, N)

Outer product function \(M \wedge N\), equivalent to M ^ N.

M and N must be from the same layout

clifford.operator.ip(M, N)

Hestenes inner product function \(M \bullet N\), equivalent to M | N.

M and N must be from the same layout

Changed in version 1.3.0: These functions used to be in clifford, but have been moved to this submodule.

transformations (clifford.transformations)

New in version 1.3.0.

This module may in future become the home for optimized rotor transformations, or non-linear transformations.

See the Linear transformations tutorial for an introduction to how to use parts of this module.

Base classes

clifford.transformations.Transformation

A callable mapping one MultiVector to another.

alias of Callable[[clifford._multivector.MultiVector], clifford._multivector.MultiVector]

class clifford.transformations.Linear(*args, **kwds)

A transformation which is linear, such that for scalar \(a_i\), \(f(a_1 x_1 + a_2 x_2) = a_1 f(x_1) + a_2 f(x_2)\).

class clifford.transformations.FixedLayout(layout_src: clifford._layout.Layout, layout_dst: Optional[clifford._layout.Layout] = None)

A transformation with a fixed source and destination layout.

Parameters

• layout_src (Layout of S dimensions) – The layout from which this transformation takes multivectors as input

• layout_dst (Layout of D dimensions) – The layout in which this transformation produces multivectors as output. Defaults to the same as the input.

Matrix-backed implementations

class clifford.transformations.LinearMatrix(matrix, layout_src: clifford._layout.Layout, layout_dst: Optional[clifford._layout.Layout] = None)

Linear transformation implemented by a matrix

Transformations need not be grade preserving.

Parameters

• matrix ((2**D, 2**S) array_like) – A matrix that transforms multivectors from layout_src with \(2^S\) elements to multivectors in layout_dst with \(2^D\) elements, by left-multiplication.

• layout_src (Layout of S dimensions) – Passed on to FixedLayout.

• layout_dst (Layout of D dimensions) – Passed on to FixedLayout.

See also

clifford.BladeMap

A faster but less general approach that works on basis blades

property adjoint: clifford.transformations.LinearMatrix

The adjoint transformation

classmethod from_function(func: Callable[[clifford._multivector.MultiVector], clifford._multivector.MultiVector], layout_src: clifford._layout.Layout, layout_dst: Optional[clifford._layout.Layout] = None) → clifford.transformations.LinearMatrix

Build a linear transformation from the result of a function applied to each basis blade.

Parameters

• func – A function taking basis blades from layout_src that produces multivectors in layout_dst.

• layout_src (Layout of S dimensions) – The layout to pass into the generating function

• layout_dst (Layout of D dimensions) – The layout the generating function is expected to produce. If not passed, this is inferred.

Example

>>> from clifford import transformations, Layout
>>> l = Layout([1, 1])
>>> e1, e2 = l.basis_vectors_lst
>>> rot_90 = transformations.LinearMatrix.from_function(lambda x: (1 + e1*e2)*x*(1 - e1*e2)/2, l)
>>> rot_90(e1)
(1.0^e2)
>>> rot_90(e2)
-(1.0^e1)
>>> rot_90(e1*e2)
(1.0^e12)

See also

LinearMatrix.from_rotor

a shorter way to spell the above example

clifford.linear_operator_as_matrix

a lower-level function for working with a subset of basis blades

classmethod from_rotor(rotor: clifford._multivector.MultiVector) → clifford.transformations.LinearMatrix

Build a linear transformation from the result of applying a rotor sandwich.

The resulting transformation operates within the algebra of the provided rotor.

Parameters

rotor – The rotor to apply

Example

>>> from clifford import transformations, Layout
>>> l = Layout([1, 1])
>>> e1, e2 = l.basis_vectors_lst
>>> rot_90 = transformations.LinearMatrix.from_rotor(1 + e1*e2)
>>> rot_90(e1)
(1.0^e2)
>>> rot_90(e2)
-(1.0^e1)
>>> rot_90(e1*e2)
(1.0^e12)

class clifford.transformations.OutermorphismMatrix(matrix, layout_src: clifford._layout.Layout, layout_dst: Optional[clifford._layout.Layout] = None)

A generalization of a linear transformation to vectors via the outer product.

Namely, given a linear transformation \(F(u) \to v\), this generalizes to the blades by outermorphism, \(F(u_1 \wedge u_2) \to F(u_1) \wedge F(u_2)\), and to the multivectors by distributivity.

Such a transformation is grade preserving.

See [DFM09] Chapter 4 for more information

Parameters

• matrix ((D, S) array_like) – A matrix that transforms vectors from layout_src of size S to vectors in layout_dst of size D by left-multiplication.

• layout_src (Layout of S dimensions) – Passed on to FixedLayout.

• layout_dst (Layout of D dimensions) – Passed on to FixedLayout.

Example

We can construct a simple transformation that permutes and non-uniformly scales the basis vectors:

>>> from clifford import transformations, Layout
>>> layout = Layout([1, 1, 1])
>>> e1, e2, e3 = layout.basis_vectors_lst
>>> layout_new = Layout([1, 1, 1], names='f')
>>> m = np.array([[0, 1, 0],
...               [0, 0, 2],
...               [3, 0, 0]])
>>> lt = transformations.OutermorphismMatrix(m, layout, layout_new)

Applying it to some multivectors:

>>> # the transformation we specified
>>> lt(e1), lt(e2), lt(e3)
((3^f3), (1^f1), (2^f2))

>>> # the one deduced by outermorphism
>>> lt(e1^e2), lt(e2^e3), lt(e1^e3)
(-(3^f13), (2^f12), -(6^f23))

>>> # and by distributivity
>>> lt(1 + (e1^e2^e3))
1 + (6^f123)

Helper functions

clifford.transformations.between_basis_vectors(layout_src: clifford._layout.Layout, layout_dst: clifford._layout.Layout, mapping: Optional[Dict[Any, Any]] = None) → clifford.transformations.OutermorphismMatrix

Construct an outermorphism that maps basis vectors from one layout to basis vectors in another.

Parameters

• layout_src (Layout) – Passed on to FixedLayout.

• layout_dst (Layout) – Passed on to FixedLayout.

• mapping – If provided, a dictionary mapping the ids of source basis vectors to the ids of destination basis vectors. For example, {1: 2, 2: 3, 3: 1} would permute the basis vectors of clifford.g3.

Example

See the tutorial on Working with custom algebras for a motivating example.

taylor_expansions (clifford.taylor_expansions)

New in version 1.4.0.

This file implements various Taylor expansions for useful functions of multivectors. For some algebra signatures there may exist closed forms of these functions which would likely be faster and more accurate. Nonetheless, having pre-written taylor expansions for the general case is useful.

Note

Many of these functions are also exposed as MultiVector methods, such as clifford.MultiVector.sin(). This means that mv.sin() or even np.sin(mv) can be used as a convenient interface to functions in this module, without having to import it directly.

For example:

>>> from clifford.g3 import *
>>> import numpy as np
>>> np.sin(np.pi*e12/4)
(0.86867^e12)

Implemented functions

clifford.taylor_expansions.exp(x, max_order=15)

This implements the series expansion of \(\exp x\) where \(x\) is a multivector The parameter max_order is the maximum order of the taylor series to use

clifford.taylor_expansions.sin(X, max_order=30)

A taylor series expansion for sin The parameter max_order is the maximum order of the taylor series to use

clifford.taylor_expansions.cos(X, max_order=30)

A taylor series expansion for cos The parameter max_order is the maximum order of the taylor series to use

clifford.taylor_expansions.tan(X, max_order=30)

The tan function as the ratio of sin and cos The parameter max_order is the maximum order of the taylor series to use

Note

It would probably be better to implement this as its own taylor series. This function is not JITed as currently we do not overload the truediv operator for multivectors.

clifford.taylor_expansions.sinh(X, max_order=30)

A taylor series expansion for sinh The parameter max_order is the maximum order of the taylor series to use

clifford.taylor_expansions.cosh(X, max_order=30)

A taylor series expansion for cosh The parameter max_order is the maximum order of the taylor series to use

clifford.taylor_expansions.tanh(X, max_order=30)

The tanh function as the ratio of sinh and cosh The parameter max_order is the maximum order of the taylor series to use

Note

It would probably be better to implement this as its own taylor series. This function is not JITed as currently we do not overload the truediv operator for multivectors.

numba extension support (clifford.numba)

New in version 1.4.0.

This module provides numba extension types MultiVectorType and LayoutType corresponding to MultiVector and Layout.

You do not need to import this module to take advantage of these types; they are needed only directly when writing numba overloads via numba.extending.overload() and similar.

As a simple example, the following code defines a vectorized up() function for CGA

from clifford.g3c import *

@numba.njit
def jit_up(x):
    return eo + x + 0.5*abs(x)**2*einf

assert up(e1) == jit_up(e1)

Note that a rough equivalent to this particular function is provided elsewhere as clifford.tools.g3c.fast_up().

Supported operations

The following list of operations are supported in a jitted context:

• MultiVector: A limited version of the constructor supporting only MultiVector(layout, value) and MultiVector(layout, dtype=dtype).

• layout.MultiVector(), with the same caveats as above.

• layout.dims

• layout.gaDims

• layout.sig

• MultiVector.value

• MultiVector.layout

• Arithmetic:

  • MultiVector.__add__()

  • MultiVector.__sub__()

  • MultiVector.__mul__()

  • MultiVector.__xor__()

  • MultiVector.__or__()

  • MultiVector.__pow__()

  • MultiVector.__truediv__()

  • MultiVector.__invert__()

  • MultiVector.__pos__()

  • MultiVector.__neg__()

• MultiVector.__call__()

• MultiVector.mag2()

• MultiVector.__abs__()

• MultiVector.normal()

• MultiVector.leftLaInv()

• MultiVector.shirokov_inverse()

• MultiVector.hitzer_inverse()

• MultiVector.gradeInvol()

• MultiVector.conjugate()

• MultiVector.commutator()

• MultiVector.anticommutator()

• MultiVector.even

• MultiVector.odd

Performance considerations

While the resulted jitted code is typically faster, there are two main performance issues to consider. The first is the startup time of @jiting. This can be quite substantial, although can be somewhat alleviated by using the cache=True argument to numba.jit(). The second is the time taken for numba to find the appropriate dispatch loop given the Python types of its arguments, which adds overhead to every call.

clifford tries as hard as possible to reduce this second overhead, by using the undocumented _numba_type_ attribute and keeping our own optimized cache instead of going through the recommended @numba.extending.typeof_impl.register(LayoutType) mechanism. However, this overhead can still be slow compared to elementary operations. The following code is significantly impacted by this overhead:

from clifford.g3c import *
import numba

@numba.njit
def mul(a, b):
    return a * b

# 286 ms, ignoring jitting time
x = e1
for i in range(10000):
    x = mul(x, x + e2)

as each iteration of the loop pays it again. The overhead can be avoided by jitting the entire loop:

from clifford.g3c import *
import numba

@numba.njit
def mul(a, b):
    return a * b

@numba.njit
def the_loop(x):
    for i in range(1000):
        x = mul(x, x + e1)
    return x

# 2.4 ms, ignoring jitting time
x = the_loop(eq)

Predefined Algebras

The easiest way to get started with clifford is to use one of several predefined algebras:

• g2: 2D Euclidean, Cl(2). See Quick Start (G2) for some examples.

• g3: 3D Euclidean, Cl(3). See The Algebra Of Space (G3) for some examples.

• g4: 4D Euclidean, Cl(4).

• g2c: Conformal space for G2, Cl(3, 1). See Conformal Geometric Algebra for some examples.

• g3c: Conformal space for G3, Cl(4, 1).

• pga: Projective space for G3 Cl(3, 0, 1).

• gac: Geometric Algebra for Conics, Cl(5, 3).

• dpga: Double PGA also referred to as the Mother Algebra, Cl(4, 4).

• dg3c: Double Conformal Geometric Algebra, effectively two g3c algebras glued together Cl(8, 2).

By using the pre-defined algebras in place of calling Cl directly, you will often find that your program starts up faster.

clifford.<predefined>.e<ijk>

All of these modules expose the basis blades as attributes, and can be used like so

In [1]: from clifford import g2

In [2]: g2.e1 * g2.e2
Out[2]: (1^e12)

Additionally, they define the following attributes, which contain the return values of clifford.Cl():

clifford.<predefined>.layout

The associated clifford.Layout

In [3]: g2.layout
Out[3]: 
Layout([1, 1],
       ids=BasisVectorIds.ordered_integers(2),
       order=BasisBladeOrder.shortlex(2),
       names=['', 'e1', 'e2', 'e12'])

clifford.<predefined>.blades

A shorthand for Layout.blades()

In [4]: g2.blades
Out[4]: {'': 1, 'e1': (1^e1), 'e2': (1^e2), 'e12': (1^e12)}

For interactive use, it’s very handy to use import *

In [5]: from clifford.g2 import *

In [6]: e1, e2, e12
Out[6]: ((1^e1), (1^e2), (1^e12))

For the conformal layouts g2c and g3c, the full contents of the stuff result of clifford.conformalize() is also exposed as members of the module.

Changelog

Changes in 1.4.x

• Python 3.9 is officially supported.

• MultiVector is now supported directly within @numba.njited code. See the documentation for this feature in the clifford.numba module for more details.

• New algorithms for the multivector inverse MultiVector.inv():

  • Inverting a non-blade multivector in algebras where \(p = q \le 5\) now falls back on the approach described in [HS17] instead of using a linear algebra approach. This algorithm can be used directly via MultiVector.hitzer_inverse().

  • An additional method is available, MultiVector.shirokov_inverse(), which is the arbitrary signature algorithm described in [Shi20].

• A new clifford.taylor_expansions module for Taylor series of various multivector functions, starting with common trigonometric functions. These functions are additionally exposed via methods on MultiVector like MultiVector.cos().

• Random functions now accept an rng keyword argument that accepts the object returned by numpy.random.default_rng(), for deterministic randomness.

• Projection using MultiVector.__call__() as mv(grade) no longer raises ValueError for grades not present in the algebra, and instead just returns zero.

• Some JIT-ed code is now cached when clifford is imported for the very first time, speeding up subsequent imports.

• Improved citations in the documentation.

Bugs fixed

• Where possible, MultiVectors preserve their data type in the dual, and the right and left complements.

• MVArray no longer errantly promotes 0-dimensional arrays to 1-dimensional arrays.

• MultiVectors with complex coefficients are now printed correctly.

• cga.Round.radius() is no longer always 1.

Compatibility notes

• This will be the last release to support Python 3.5 and 3.6, as the former reached its end-of-life during our last release cycle, and the latter is expected to before our next release.

• clifford.general_exp() is deprecated in favor of clifford.taylor_expansions.exp(), although typically clifford.MultiVector.exp() is a better choice

• Transparently treating a MultiVector as a flat array of coefficients is deprecated, and so mv[i] and len(mv) both emit DeprecationWarnings. If the underlying storage order is of interest, use mv.value[i] and len(mv) respectively. To obtain the scalar part of a MultiVector, use mv[()] instead of mv[0].

• Layout.gradeList has been removed. If still needed, it can be recovered as an ndarray isntead of a list via the private attribute layout._basis_blade_order.grades (BasisBladeOrder.grades).

Changes in 1.3.x

• Python 3.8 is officially supported. 1.2.0 was pinned to a bad numba version that was incompatible with 3.8.

• A new clifford.operator module to contain the previously undocumented gp(), op(), and ip() helpers.

• A new clifford.transformations module for linear transformations.

• Two new Predefined Algebras, clifford.dpga and clifford.dg3c.

• Improvements throughout the documentation:

  • Better overall structure, visible in the docs sidebar.

  • New tutorials for Conformal Geometric Algebra on visualization and applications.

  • New tutorial on Working with custom algebras.

  • New tutorial on Linear transformations.

  • New links at the top of each notebook tutorial, to run examples from the browser.

• Faster algebra construction. Cl(3) is now 100× faster, and Cl(6) is 20× faster. This is achieved by deferring product JIT compilation until the product is used for the first time.

• Additional testing and assorted improvements for clifford.tools.g3c:

  • closest_point_on_circle_from_line() has now been implemented roughly following the procedure described in Appendix A of Andreas Aristidou’s PhD thesis.

  • closest_point_on_line_from_circle() has now also been added, projecting the result of closest_point_on_circle_from_line() to the line.

• clifford.ugly() now results in less ugly output for Predefined Algebras.

Bugs fixed

• MultiVector.meet() no longer produces zero erroneously.

• mv[e1 + e12] now raises ValueError, rather than being interpreted as mv[e1].

• ip() (the inner product) no longer performs the outer product.

• Layout.parse_multivector() now throws SyntaxError on invalid input, rather than silenltly producing nonsense.

• Layout.parse_multivector() supports basis vector names which do not start with e.

• In clifford.tools.g3c:

  • val_midpoint_between_lines() now handles the case that the two lines are touching.

  • val_fit_circle() now correctly selects the first and second eigenvalue regardless of order.

  • sphere_beyond_plane() now tested and correct.

  • sphere_behind_plane() now tested and correct.

  • val_unsign_sphere() is now jitted, as it should have been from the start.

  • get_nearest_plane_point() correctly returns the conformal point rather than the 3D point.

Compatibility notes

• clifford.grades_present is deprecated in favor of MultiVector.grades(), the latter of which now takes an eps argument.

• del mv[i] is no longer legal, the equivalent mv[i] = 0 should be used instead.

• Layout.dict_to_multivector has been removed. It was accidentally broken in 1.0.5, so there is little point deprecating it.

• Layout.basis_names() now returns a list of str, rather than a numpy array of bytes. The result now matches the construction order, rather than being sorted alphabetically. The order of Layout.metric() has been adjusted for consistency.

• The imt_prod_mask, omt_prod_mask, and lcmt_prod_mask attributes of Layout objects have been removed, as these were an unnecessary intermediate computation that had no need to be public.

• Some functions in clifford.tools.g3c have been renamed:

  • closest_points_on_circles has been renamed to iterative_closest_points_on_circles().

  • closest_points_circle_line has been renamed to iterative_closest_points_circle_line().

  • furthest_points_on_circles has been renamed to iterative_furthest_points_on_circles().

• While this release is compatible with numba version 0.49.0, it is recommended to use 0.48.0 which does not emit as many warnings. See the Installation instructions for how to follow this guidance.

Patch releases

• 1.3.1: Added compatibility with numba version 0.50.0.

Changes in 1.2.x

• layout.isconformal, layout.einf, and layout.eo, which were added in 1.0.4, have been removed. The first can now be spelt isinstance(layout, clifford.ConformalLayout), and the other properties now exist only on ConformalLayout objects.

• MultiVector.left_complement() has been added for consistency with MultiVector.right_complement().

• A new clifford.tools.classify module has been added for classifying blades.

• Layout objects print slightly more cleanly in Jupyter notebooks.

• Layout.scalar is now integral rather than floating point

Bugs fixed

• pow(mv, 0) gives the right result

• nan is now printed correctly when it appears in multivectors. Previously it was hidden

• MultiVector.right_complement() no longer performs the left complement.

• MultiVector.vee() has been corrected to have the same sign as MultiVector.meet()

Compatibility notes

• Layout.scalar is now integral rather than floating point, to match Layout.pseudoScalar.

Changes in 1.1.x

• Restores layout.gmt, Layout.omt, Layout.imt, and Layout.lcmt. A few releases ago, these existed but were dense. For memory reasons, they were then removed entirely. They have now been reinstated as sparse.COO matrix objects, which behave much the same as the original dense arrays.

• MultiVectors preserve their data type in addition, subtraction, and products. This means that integers remain integers until combined with floats. Note that this means in principle integer overflow is possible, so working with floats is still recommended. This also adds support for floating point types of other precision, such as np.float32.

• setup.py is now configured such that pip2 install clifford will not attempt to download this version, since it does not work at all on python 2.

• Documentation now includes examples of pyganja visualizations.

Compatibility notes

• Layout.blades() now includes the scalar 1, as do other similar functions.

• MultiVector.grades() now returns a set not a list. This means code like mv.grades() == [0] will need to change to mv.grades() == {0}, or to work both before and after this change, set(mv.grades()) == {0}.

Bugs fixed

• mv[(i, j)] would sometimes fail if the indices were not in canonical order.

• mv == None and layout == None would crash rather than return False.

• blade.isVersor() would return False.

• layout.blades_of_grade(0) would not return the list it claimed to return.

Internal changes

• Switch to pytest for testing.

• Enable code coverage.

• Split into smaller files.

• Remove python 2 compatibility code, which already no longer worked.

Changes 0.6-0.7

• Added a real license.

• Convert to NumPy instead of Numeric.

Changes 0.5-0.6

• join() and meet() actually work now, but have numerical accuracy problems

• added clean() to MultiVector

• added leftInv() and rightInv() to MultiVector

• moved pseudoScalar() and invPS() to MultiVector (so we can derive new classes from MultiVector)

• changed all of the instances of creating a new MultiVector to create an instance of self.__class__ for proper inheritance

• fixed bug in laInv()

• fixed the massive confusion about how dot() works

• added left-contraction

• fixed embarrassing bug in gmt generation

• added normal() and anticommutator() methods

• fixed dumb bug in elements() that limited it to 4 dimensions

Acknowledgements

Konrad Hinsen fixed a few bugs in the conversion to numpy and adding some unit tests.

Issues

Warning

This document is kept for historic reasons, but may no longer reflect the current state of the latest release of clifford. For the most up to date source of issues, look at the GitHub issue tracker.

• Currently, algebras over 6 dimensions are very slow. this is because this module was written for pedagogical purposes. However, because the syntax for this module is so attractive, we plan to fix the performance problems, in the future…

• Due to Python’s order of operations, the bit operators ^ << | are evaluated after the normal arithmetic operators + - * /, which do not follow the precedence expected in GA

  # written        meaning            possibly intended
  1^e1 + 2^e2   == 1^(e1+2)^e2     != (1^e0) + (2^e1)
  e2 + e1|e2    == (e2 + e1)|e2    != e1 + (e1|e2)
  

  This can also cause confusion within the bitwise operators:

  # written        meaning            possibly intended
  e1 << e2 ^ e1 == (e1 << e2) ^ e1 != e1 << (e2 ^ e1)
  e1 ^ e2 | e1  == (e1 << e2) ^ e1 != e1 << (e2 ^ e1)
  

• Since | is the inner product and the inner product with a scalar vanishes by definition, an expression like:

  (1|e0) + (2|e1)
  

  is null. Use the outer product or full geometric product, to multiply scalars with MultiVectors. This can cause problems if one has code that mixes Python numbers and MultiVectors. If the code multiplies two values that can each be either type without checking, one can run into problems as 1 | 2 has a very different result from the same multiplication with scalar MultiVectors.

• Taking the inverse of a MultiVector will use a method proposed by Christian Perwass that involves the solution of a matrix equation. A description of that method follows:

  Representing multivectors as \(2^\text{dims}\)-vectors (in the matrix sense), we can carry out the geometric product with a multiplication table. In pseudo-tensorish language (using summation notation)

  \[m_i g_{ijk} n_k = v_j\]

  Suppose \(m_i\) are known (M is the vector we are taking the inverse of), the \(g_{ijk}\) have been computed for this algebra, and \(v_j = 1\) if the \(j\)’th element is the scalar element and 0 otherwise, we can compute the dot product \(m_i g_{ijk}\). This yields a rank-2 matrix. We can then use well-established computational linear algebra techniques to solve this matrix equation for \(n_k\). The laInv method does precisely that.

  The usual, analytic, method for computing inverses (\(M^{-1} = \tilde M/(M \tilde M)\) iff \(M\tilde M = {|M|}^2\)) fails for those multivectors where M*~M is not a scalar. It is onl)y used if the inv method is manually set to point to normalInv.

  My testing suggests that laInv works. In the cases where normalInv works, laInv returns the same result (within _eps). In all cases, M * M.laInv() == 1.0 (within _eps). Use whichever you feel comfortable with.

  Of course, a new issue arises with this method. The inverses found are sometimes dependant on the order of multiplication. That is:

  M.laInv() * M == 1.0
  M * M.laInv() != 1.0
  

  XXX Thus, there are two other methods defined, leftInv and rightInv which point to leftLaInv and rightLaInv. The method inv points to rightInv. Should the user choose, leftInv and rightInv will both point to normalInv, which yields a left- and right-inverse that are the same should either exist (the proof is fairly simple).

• The basis vectors of any algebra will be orthonormal unless you supply your own multiplication tables (which you are free to do after the Layout constructor is called). A derived class could be made to calculate these tables for you (and include methods for generating reciprocal bases and the like).

• No care is taken to preserve the dtype of the arrays. The purpose of this module is pedagogical. If your application requires so many multivectors that storage becomes important, the class structure here is unsuitable for you anyways. Instead, use the algorithms from this module and implement application-specific data structures.

• Conversely, explicit typecasting is rare. MultiVectors will have integer coefficients if you instantiate them that way. Dividing them by Python integers will have the same consequences as normal integer division. Public outcry will convince me to add the explicit casts if this becomes a problem.

---

Happy hacking!

Robert Kern

robert.kern@gmail.com

This page was generated from docs/tutorials/g2-quick-start.ipynb. Interactive online version: .

Quick Start (G2)

This notebook gives a terse introduction to using the clifford module, using a two-dimensional geometric algebra as the context.

Setup

First, import clifford and instantiate a two-dimensional algebra (G2),

import clifford as cf

layout, blades = cf.Cl(2) # creates a 2-dimensional clifford algebra

Inspect blades.

blades

{'': 1, 'e1': (1^e1), 'e2': (1^e2), 'e12': (1^e12)}

Assign blades to variables

e1 = blades['e1']
e2 = blades['e2']
e12 = blades['e12']

Basics

e1*e2  # geometric product

(1^e12)

e1|e2  # inner product

0

e1^e2  # outer product

(1^e12)

Reflection

a = e1+e2     # the vector
n = e1        # the reflector
-n*a*n.inv()  # reflect `a` in hyperplane normal to `n`

-(1.0^e1) + (1.0^e2)

Rotation

from math import e, pi

R = e**(pi/4*e12)  # enacts rotation by pi/2
R

0.70711 + (0.70711^e12)

R*e1*~R    # rotate e1 by pi/2 in the e12-plane

-(1.0^e2)

This page was generated from docs/tutorials/g3-algebra-of-space.ipynb. Interactive online version: .

The Algebra Of Space (G3)

In this notebook, we give a more detailed look at how to use clifford, using the algebra of three dimensional space as a context.

Setup

First, we import clifford as cf, and instantiate a three dimensional geometric algebra using Cl() (docs).

import clifford as cf

layout, blades = cf.Cl(3)  # creates a 3-dimensional clifford algebra

Given a three dimensional GA with the orthonormal basis,

\[e_{i}\cdot e_{j}=\delta_{ij}\]

The basis consists of scalars, three vectors, three bivectors, and a trivector.

\[\{\hspace{0.5em} \underbrace{\hspace{0.5em}\alpha,\hspace{0.5em}}_{\mbox{scalar}}\hspace{0.5em} \underbrace{\hspace{0.5em}e_{1},\hspace{1.5em}e_{2},\hspace{1.5em}e_{3},\hspace{0.5em}}_{\mbox{vectors}}\hspace{0.5em} \underbrace{\hspace{0.5em}e_{12},\hspace{1.5em}e_{23},\hspace{1.5em}e_{13},\hspace{0.5em}}_{\mbox{bivectors}}\hspace{0.5em} \underbrace{\hspace{0.5em}e_{123}\hspace{0.5em}}_{\text{trivector}} \hspace{0.5em} \}\]

Cl() creates the algebra and returns a layout and blades. The layout holds information and functions related this instance of G3, and the blades is a dictionary which contains the basis blades, indexed by their string representations,

blades

{'': 1,
 'e1': (1^e1),
 'e2': (1^e2),
 'e3': (1^e3),
 'e12': (1^e12),
 'e13': (1^e13),
 'e23': (1^e23),
 'e123': (1^e123)}

You may wish to explicitly assign the blades to variables like so,

e1 = blades['e1']
e2 = blades['e2']
# etc ...

Or, if you’re lazy and just working in an interactive session you can use locals() to update your namespace with all of the blades at once.

locals().update(blades)

Now, all the blades have been defined in the local namespace

e3, e123

((1^e3), (1^e123))

Basics

Products

The basic products are available

e1*e2  # geometric product

(1^e12)

e1|e2  # inner product

0

e1^e2  # outer product

(1^e12)

e1^e2^e3  # even more outer products

(1^e123)

Defects in Precedence

Python’s operator precedence makes the outer product evaluate after addition. This requires the use of parentheses when using outer products. For example

e1^e2 + e2^e3  # fail, evaluates as

(2^e123)

(e1^e2) + (e2^e3)  # correct

(1^e12) + (1^e23)

Also the inner product of a scalar and a Multivector is 0,

4|e1

0

So for scalars, use the outer product or geometric product instead

4*e1

(4^e1)

Multivectors

Multivectors can be defined in terms of the basis blades. For example you can construct a rotor as a sum of a scalar and bivector, like so

import math

theta = math.pi/4
R = math.cos(theta) - math.sin(theta)*e23
R

0.70711 - (0.70711^e23)

You can also mix grades without any reason

A = 1 + 2*e1 + 3*e12 + 4*e123
A

1 + (2^e1) + (3^e12) + (4^e123)

Reversion

The reversion operator is accomplished with the tilde ~ in front of the Multivector on which it acts

~A

1 + (2^e1) - (3^e12) - (4^e123)

Grade Projection

Taking a projection onto a specific grade \(n\) of a Multivector is usually written

\[\langle A \rangle _n\]

can be done by using soft brackets, like so

A(0)  # get grade-0 elements of R

1

A(1)  # get grade-1 elements of R

(2^e1)

A(2)  # you get it

(3^e12)

Magnitude

Using the reversion and grade projection operators, we can define the magnitude of \(A\)

\[|A|^2 = \langle A\tilde{A}\rangle\]

(A*~A)(0)

30

This is done in the abs() operator

abs(A)**2

30.0

Inverse

The inverse of a Multivector is defined as \(A^{-1}A=1\)

A.inv()*A

1.0

A.inv()

0.13415 + (0.12195^e1) - (0.14634^e3) + (0.18293^e12) + (0.09756^e23) - (0.29268^e123)

Dual

The dual of a multivector \(A\) can be defined as

\[AI^{-1}\]

Where, \(I\) is the pseudoscalar for the GA. In \(G_3\), the dual of a vector is a bivector,

a = 1*e1 + 2*e2 + 3*e3
a.dual()

-(3^e12) + (2^e13) - (1^e23)

Pretty, Ugly, and Display Precision

You can toggle pretty printing with with pretty() or ugly(). ugly returns an eval-able string.

cf.ugly()
A.inv()

MultiVector(Layout([1, 1, 1],
                   ids=BasisVectorIds.ordered_integers(3),
                   order=BasisBladeOrder.shortlex(3),
                   names=['', 'e1', 'e2', 'e3', 'e12', 'e13', 'e23', 'e123']),
            [0.13414634146341464, 0.12195121951219512, 0.0, -0.14634146341463414, 0.18292682926829268, 0.0, 0.0975609756097561, -0.2926829268292683])

You can also change the displayed precision

cf.pretty(precision=2)

A.inv()

0.13 + (0.12^e1) - (0.15^e3) + (0.18^e12) + (0.1^e23) - (0.29^e123)

This does not effect the internal precision used for computations.

Applications

Reflections

Reflecting a vector \(c\) about a normalized vector \(n\) is pretty simple,

\[c \rightarrow ncn\]

c = e1+e2+e3    # a vector
n = e1          # the reflector
n*c*n           # reflect `a` in hyperplane normal to `n`

(1^e1) - (1^e2) - (1^e3)

Because we have the inv() available, we can equally well reflect in un-normalized vectors using,

\[a \rightarrow nan^{-1}\]

a = e1+e2+e3    # the vector
n = 3*e1        # the reflector
n*a*n.inv()

(1.0^e1) - (1.0^e2) - (1.0^e3)

Reflections can also be made with respect to the a ‘hyperplane normal to the vector \(n\)’, in this case the formula is negated

\[c \rightarrow -ncn^{-1}\]

Rotations

A vector can be rotated using the formula

\[a \rightarrow Ra\tilde{R}\]

Where \(R\) is a rotor. A rotor can be defined by multiple reflections,

\[R=mn\]

or by a plane and an angle,

\[R = e^{-\frac{\theta}{2}\hat{B}}\]

For example

R = math.e**(-math.pi/4*e12)  # enacts rotation by pi/2
R

0.71 - (0.71^e12)

R*e1*~R  # rotate e1 by pi/2 in the e12-plane

(1.0^e2)

Some Ways to use Functions

Maybe we want to define a function which can return rotor of some angle \(\theta\) in the \(e_{12}\)-plane,

\[R_{12} = e^{-\frac{\theta}{2}e_{12}}\]

R12 = lambda theta: math.e**(-theta/2*e12)
R12(math.pi/2)

0.71 - (0.71^e12)

And use it like this

a = e1+e2+e3
R = R12(math.pi/2)
R*a*~R

-(1.0^e1) + (1.0^e2) + (1.0^e3)

You might as well make the angle argument a bivector, so that you can control the plane of rotation as well as the angle

\[R_B = e^{-\frac{B}{2}}\]

R_B = lambda B: math.e**(-B/2)

Then you could do

R12 = R_B(math.pi/4*e12)
R23 = R_B(math.pi/5*e23)

or

R_B(math.pi/6*(e23+e12))  # rotor enacting a pi/6-rotation in the e23+e12-plane

0.93 - (0.26^e12) - (0.26^e23)

Maybe you want to define a function which returns a function that enacts a specified rotation,

\[f(B) \rightarrow \underline{R_B}(a) = R_Ba\tilde{R_B}\]

This just saves you having to write out the sandwich product, which is nice if you are cascading a bunch of rotors, like so

\[\underline{R_C}( \underline{R_B}( \underline{R_A}(a)))\]

def R_factory(B):
    def apply_rotation(a):
        R = math.e**(-B/2)
        return R*a*~R
    return apply_rotation

R = R_factory(math.pi/6*(e23+e12))  # this returns a function
R(a)  # which acts on a vector

(0.52^e1) + (0.74^e2) + (1.48^e3)

Then you can do things like

R12 = R_factory(math.pi/3*e12)
R23 = R_factory(math.pi/3*e23)
R13 = R_factory(math.pi/3*e13)

R12(R23(R13(a)))

(0.41^e1) - (0.66^e2) + (1.55^e3)

To make cascading a sequence of rotations as concise as possible, we could define a function which takes a list of bivectors \(A,B,C,..\) , and enacts the sequence of rotations which they represent on a some vector \(x\).

\[f(A,B,C,x) = \underline{R_A} (\underline{R_B} (\underline{R_C}(x)))\]

from functools import reduce

# a sequence of rotations
def R_seq(*args):
    *Bs, x = args
    R_lst = [math.e**(-B/2) for B in Bs]  # create list of Rotors from list of Bivectors
    R = reduce(cf.gp, R_lst)          # apply the geometric product to list of Rotors
    return R*x*~R

# rotation sequence by  pi/2-in-e12 THEN pi/2-in-e23
R_seq(math.pi/2*e23, math.pi/2*e12, e1)

(1.0^e3)

Changing Basis Names

If you want to use different names for your basis as opposed to e’s with numbers, supply the Cl() with a list of names. For example for a two dimensional GA,

layout, blades = cf.Cl(2, names=['','x','y','i'])

blades

{'': 1, 'x': (1^x), 'y': (1^y), 'i': (1^i)}

locals().update(blades)

1*x + 2*y

(1^x) + (2^y)

1 + 4*i

1 + (4^i)

This page was generated from docs/tutorials/euler-angles.ipynb. Interactive online version: .

Rotations in Space: Euler Angles, Matrices, and Quaternions

This notebook demonstrates how to use clifford to implement rotations in three dimensions using euler angles, rotation matices and quaternions. All of these forms are derived from the more general rotor form, which is provided by GA. Conversion from the rotor form to a matrix representation is shown, and takes about three lines of code. We start with euler angles.

Euler Angles with Rotors

A common way to parameterize rotations in three dimensions is through Euler Angles.

Any orientation can be achieved by composing three elemental rotations. The elemental rotations can either occur about the axes of the fixed coordinate system (extrinsic rotations) or about the axes of a rotating coordinate system, which is initially aligned with the fixed one, and modifies its orientation after each elemental rotation (intrinsic rotations). — wikipedia

The animation below shows an intrinsic rotation model as each elemental rotation is applied.Label the left, right, and vertical blue-axes as \(e_1, e_2,\) and \(e_3\), respectively. The series of rotations can be described:

1. rotate about \(e_3\)-axis

2. rotate about the rotated \(e_1\)-axis, call it \(e_1^{'}\)

3. rotate about the twice rotated axis of \(e_3\)-axis, call it \(e_3^{''}\)

So the elemental rotations are about \(e_3, e_{1}^{'}, e_3^{''}\)-axes, respectively.

(taken from wikipedia)

Following Sec. 2.7.5 from [DL03], we first rotate an angle \(\phi\) about the \(e_3\)-axis, which is equivalent to rotating in the \(e_{12}\)-plane. This is done with the rotor

\[R_{\phi} = e^{-\frac{\phi}{2} e_{12}}\]

Next we rotate about the rotated \(e_1\)-axis, which we label \(e_1^{'}\). To find where this is, we can rotate the axis,

\[e_1^{'} =R_{\phi} e_1 \tilde{R_{\phi}}\]

The plane corresponding to this axis is found by taking the dual of \(e_1^{'}\)

\[I R_{\phi} e_1 \tilde{R_{\phi}} = R_{\phi} e_{23} \tilde{R_{\phi}}\]

Where we have made use of the fact that the pseudo-scalar commutes in G3. Using this result, the second rotation by angle \(\theta\) about the \(e_1^{'}\)-axis is then ,

\[R_{\theta} = e^{\frac{\theta}{2} R_{\phi} e_{23} \tilde{R_{\phi}}}\]

However, noting that

\[e^{R_{\phi} e_{23} \tilde{R_{\phi}}} =R_{\phi} e^{e_{23}} \tilde{R_{\phi}}\]

Allows us to write the second rotation by angle \(\theta\) about the \(e_1^{'}\)-axis as

\[R_{\theta} = R_{\phi} e^{\frac{\theta}{2}e_{23}} \tilde{R_{\phi}}\]

So, the combination of the first two elemental rotations equals,

\begin{align*} R_{\theta} R_{\phi} &= R_{\phi} e^{\frac{\theta}{2}e_{23}} \tilde{R_{\phi}} R_{\phi} \\ &= e^{-\frac{\phi}{2} e_{12}}e^{-\frac{\theta}{2} e_{23}} \end{align*}

This pattern can be extended to the third elemental rotation of angle \(\psi\) in the twice-rotated \(e_1\)-axis, creating the total rotor

\[R = e^{-\frac{\phi}{2} e_{12}} e^{-\frac{\theta}{2} e_{23}} e^{-\frac{\psi}{2} e_{12}}\]

Implementation of Euler Angles

First, we initialize the algebra and assign the variables

from numpy import e,pi
from clifford import Cl

layout, blades = Cl(3)   # create a 3-dimensional clifford algebra

locals().update(blades) # lazy way to put entire basis in the namespace

Next we define a function to produce a rotor given euler angles

def R_euler(phi, theta,psi):
    Rphi = e**(-phi/2.*e12)
    Rtheta = e**(-theta/2.*e23)
    Rpsi = e**(-psi/2.*e12)

    return Rphi*Rtheta*Rpsi

For example, using this to create a rotation similar to that shown in the animation above,

R = R_euler(pi/4, pi/4, pi/4)
R

0.65328 - (0.65328^e12) - (0.38268^e23)

Convert to Quaternions

A Rotor in 3D space is a unit quaternion, and so we have essentially created a function that converts Euler angles to quaternions. All you need to do is interpret the bivectors as \(i,j,\) and \(k\)’s. See Interfacing Other Mathematical Systems, for more on quaternions.

Convert to Rotation Matrix

The matrix representation for a rotation can defined as the result of rotating an ortho-normal frame. Rotating an ortho-normal frame can be done easily,

A = [e1, e2, e3]         # initial ortho-normal frame
B = [R*a*~R for a in A]  # resultant frame after rotation

B

[(0.14645^e1) + (0.85355^e2) + (0.5^e3),
 -(0.85355^e1) - (0.14645^e2) + (0.5^e3),
 (0.5^e1) - (0.5^e2) + (0.70711^e3)]

The components of this frame are the rotation matrix (although note that each vector in the frame corresponds to a row of the rotation matrix), so we just enter the frame components as columns into a matrix.

import numpy as np

# b|a is a multivector, `[()]` selects the scalar part
M = np.array([
    [(b|a)[()] for b in B]
    for a in A
])
M  # rows correspond to A, columns to B

array([[ 0.14644661, -0.85355339,  0.5       ],
       [ 0.85355339, -0.14644661, -0.5       ],
       [ 0.5       ,  0.5       ,  0.70710678]])

That’s a rotation matrix. We can compare transforming the vector \(e_1 + 2e_2\) via the rotor and via the matrix:

X = e1+2*e2
R*X*~R

-(1.56066^e1) + (0.56066^e2) + (1.5^e3)

Xv = np.array([1, 2, 0])
M @ Xv

array([-1.56066017,  0.56066017,  1.5       ])

Convert a Rotation Matrix to a Rotor

In 3 Dimensions, there is a simple formula which can be used to directly transform a rotations matrix into a rotor. For arbitrary dimensions you have to use a different algorithm (see clifford.tools.orthoMat2Versor() (docs)). Anyway, in 3 dimensions there is a closed form solution, as described in Sec. 4.3.3 of [DL03]. Given a rotor \(R\) which transforms an orthonormal frame \(A={a_k}\) into \(B={b_k}\) as such,

\[b_k = Ra_k\tilde{R}\]

\(R\) is given by

\[R= \frac{1+a_kb_k}{|1+a_kb_k|}\]

So, if you want to convert from a rotation matrix into a rotor, start by converting the matrix M into a frame \(B\).(You could do this with loop if you want.)

B = [M[0,0]*e1 + M[0,1]*e2 + M[0,2]*e3,
     M[1,0]*e1 + M[1,1]*e2 + M[1,2]*e3,
     M[2,0]*e1 + M[2,1]*e2 + M[2,2]*e3]
B

[(0.14645^e1) - (0.85355^e2) + (0.5^e3),
 (0.85355^e1) - (0.14645^e2) - (0.5^e3),
 (0.5^e1) + (0.5^e2) + (0.70711^e3)]

Then implement the formula

A = [e1,e2,e3]
R = 1+sum([A[k]*B[k] for k in range(3)])
R = R/abs(R)

R

0.65328 - (0.65328^e12) - (0.38268^e23)

which matches the original rotor

R_euler(pi/4, pi/4, pi/4)

0.65328 - (0.65328^e12) - (0.38268^e23)

This page was generated from docs/tutorials/space-time-algebra.ipynb. Interactive online version: .

Space Time Algebra

Intro

This notebook demonstrates how to use clifford to work with Space Time Algebra. The Pauli algebra of space \(\mathbb{P}\), and Dirac algebra of space-time \(\mathbb{D}\), are related using the spacetime split. The split is implemented by using a BladeMap (docs), which maps a subset of blades in \(\mathbb{D}\) to the blades in \(\mathbb{P}\). This split allows a spacetime bivector \(F\) to be broken up into relative electric and magnetic fields in space. Lorentz transformations are implemented as rotations in \(\mathbb{D}\), and the effects on the relative fields are computed with the split.

Setup

First we import clifford, instantiate the two algebras, and populate the namespace with the blades of each algebra. The elements of \(\mathbb{D}\) are prefixed with \(d\), while the elements of \(\mathbb{P}\) are prefixed with \(p\). Although unconventional, it is easier to read and to translate into code.

from clifford import Cl, pretty

pretty(precision=1)

# Dirac Algebra  `D`
D, D_blades = Cl(1,3, firstIdx=0, names='d')

# Pauli Algebra  `P`
P, P_blades = Cl(3, names='p')

# put elements of each in namespace
locals().update(D_blades)
locals().update(P_blades)

The Space Time Split

To two algebras can be related by the spacetime-split. First, we create a BladeMap which relates the bivectors in \(\mathbb{D}\) to the vectors/bivectors in \(\mathbb{P}\). The scalars and pseudo-scalars in each algebra are equated.

from clifford import BladeMap

bm = BladeMap([(d01,p1),
               (d02,p2),
               (d03,p3),
               (d12,p12),
               (d23,p23),
               (d13,p13),
               (d0123, p123)])

Splitting a space-time vector (an event)

A vector in \(\mathbb{D}\), represents a unique place in space and time, i.e. an event. To illustrate the split, create a random event \(X\).

X = D.randomV()*10
X

-(7.0^d0) - (9.4^d1) - (0.6^d2) + (8.8^d3)

This can be split into time and space components by multiplying with the time-vector \(d_0\),

X*d0

-7.0 + (9.4^d01) + (0.6^d02) - (8.8^d03)

and applying the BladeMap, which results in a scalar+vector in \(\mathbb{P}\)

bm(X*d0)

-7.0 + (9.4^p1) + (0.6^p2) - (8.8^p3)

The space and time components can be separated by grade projection,

x = bm(X*d0)
x(0) # the time component

-7.0

x(1) # the space component

(9.4^p1) + (0.6^p2) - (8.8^p3)

We therefor define a split() function, which has a simple condition allowing it to act on a vector or a multivector in \(\mathbb{D}\). Splitting a spacetime bivector will be treated in the next section.

def split(X):
    return bm(X.odd*d0+X.even)

split(X)

-7.0 + (9.4^p1) + (0.6^p2) - (8.8^p3)

The split can be inverted by applying the BladeMap again, and multiplying by \(d_0\)

x = split(X)
bm(x)*d0

-(7.0^d0) - (9.4^d1) - (0.6^d2) + (8.8^d3)

Splitting a Bivector

Given a random bivector \(F\) in \(\mathbb{D}\),

F = D.randomMV()(2)
F

-(0.9^d01) + (0.2^d02) - (0.4^d03) - (1.7^d12) - (0.0^d13) - (1.1^d23)

\(F\) splits into a vector/bivector in \(\mathbb{P}\)

split(F)

-(0.9^p1) + (0.2^p2) - (0.4^p3) - (1.7^p12) - (0.0^p13) - (1.1^p23)

If \(F\) is interpreted as the electromagnetic bivector, the Electric and Magnetic fields can be separated by grade

E = split(F)(1)
iB = split(F)(2)

E

-(0.9^p1) + (0.2^p2) - (0.4^p3)

iB

-(1.7^p12) - (0.0^p13) - (1.1^p23)

Lorentz Transformations

Lorentz Transformations are rotations in \(\mathbb{D}\), which are implemented with Rotors. A rotor in G4 will, in general, have scalar, bivector, and quadvector components.

R = D.randomRotor()
R

25.4 + (17.5^d01) + (18.9^d02) - (0.7^d03) + (1.5^d12) + (8.4^d13) - (2.1^d23) - (7.8^d0123)

In this way, the effect of a lorentz transformation on the electric and magnetic fields can be computed by rotating the bivector with \(F \rightarrow RF\tilde{R}\)

F_ = R*F*~R
F_

-(2554.1^d01) + (2172.2^d02) - (1599.7^d03) - (3328.5^d12) + (688.3^d13) + (1499.4^d23)

Then splitting into \(E\) and \(B\) fields

E_ = split(F_)(1)
E_

-(2554.1^p1) + (2172.2^p2) - (1599.7^p3)

iB_ = split(F_)(2)
iB_

-(3328.5^p12) + (688.3^p13) + (1499.4^p23)

Lorentz Invariants

Since lorentz rotations in \(\mathbb{D}\), the magnitude of elements of \(\mathbb{D}\) are invariants of the lorentz transformation. For example, the magnitude of electromagnetic bivector \(F\) is invariant, and it can be related to \(E\) and \(B\) fields in \(\mathbb{P}\) through the split,

i = p123
E = split(F)(1)
B = -i*split(F)(2)

F**2

-2.9 + (3.2^d0123)

split(F**2) == E**2 - B**2 + (2*E|B)*i

True

This page was generated from docs/tutorials/InterfacingOtherMathSystems.ipynb. Interactive online version: .

Interfacing Other Mathematical Systems

Geometric Algebra is known as a universal algebra because it subsumes several other mathematical systems. Two algebras commonly used by engineers and scientists are complex numbers and quaternions. These algebras can be subsumed as the even sub-algebras of 2 and 3 dimensional geometric algebras, respectively. This notebook demonstrates how clifford can be used to incorporate data created with these systems into geometric algebra.

Complex Numbers

Given a two dimensional GA with the orthonormal basis,

\[e_{i}\cdot e_{j}=\delta_{ij}\]

The geometric algebra consists of scalars, two vectors, and a bivector,

\[\{\quad\underbrace{\alpha,}_{\mbox{scalar}}\qquad\underbrace{e_{1},\qquad e_{2},}_{\mbox{vector}}\qquad\underbrace{e_{12}}_{\mbox{bivector}}\quad\}\]

A complex number can be directly associated with a 2D spinor in the \(e_{12}\)-plane,

\[\underbrace{\mathbf{z}=\alpha+\beta i}_{\mbox{complex number}}\quad\Longrightarrow\quad\underbrace{Z=\alpha+\beta e_{12}}_{\mbox{2D spinor}}\]

The even subalgebra of a two dimensional geometric algebra is isomorphic to the complex numbers. We can setup translating functions which converts a 2D spinor into a complex number and vice-versa. In two dimensions the spinor can be also be mapped into vectors if desired.

Below is an illustration of the three different planes, the later two being contained within the geometric algebra of two dimensions, \(G_2\). Both spinors and vectors in \(G_2\) can be modeled as points on a plane, but they have distinct algebraic properties.

import clifford as cf
layout, blades = cf.Cl(2) # instantiate a 2D- GA
locals().update(blades)   # put all blades into local namespace


def c2s(z):
    '''convert a complex number to a spinor'''
    return z.real + z.imag*e12

def s2c(S):
    '''convert a spinor to a complex number'''
    S0 = float(S(0))
    S2 = float(-S|e12)
    return S0 + S2*1j

Convert a complex number to a spinor

c2s(1+2j)

1.0 + (2.0^e12)

Convert a spinor to a complex number

s2c(1+2*e12)

(1+2j)

Make sure we get what we started with when we make a round trip

s2c(c2s(1+2j)) == 1+2j

True

The spinor is then mapped to a vector by choosing a reference direction. This may be done by left multiplying with \(e_{1}\) .

\[Z\Longrightarrow e_{1}Z=e_{1}\alpha+\beta e_{1}e_{12}=\underbrace{\alpha e_{1}+\beta e_{2}}_{\mbox{vector}}\]

s = 1+2*e12
e1*s

(1^e1) + (2^e2)

Geometrically, this is interpreted as having the spinor rotate a specific vector, in this case \(e_1\). Building off of the previously defined functions

def c2v(c):
    '''convert a complex number to a vector'''
    return e1*c2s(c)

def v2c(v):
    '''convert a vector to a complex number'''
    return s2c(e1*v)

c2v(1+2j)

(1.0^e1) + (2.0^e2)

v2c(1*e1+2*e2)

(1+2j)

Depending on your applications, you may wish to have the bivector be an argument to the c2s and s2c functions. This allows you to map input data given in the form of complex number onto the planes of your choice. For example, in three dimensional space there are three bivector-planes; \(e_{12}, e_{23}\) and \(e_{13}\), so there are many bivectors which could be interpreted as the unit imaginary.

Complex numbers mapped in this way can be used to enact rotations within the specified planes.

import clifford as cf
layout, blades = cf.Cl(3)
locals().update(blades)

def c2s(z,B):
    '''convert a complex number to a spinor'''
    return z.real + z.imag*B

def s2c(S,B):
    '''convert a spinor to a complex number'''
    S0 = float(S(0))
    S2 = float(-S|B)
    return S0 + S2*1j

c2s(1+2j, e23)

1.0 + (2.0^e23)

c2s(3+4j, e13)

3.0 + (4.0^e13)

This brings us to the subject of quaternions, which are used to handle rotations in three dimensions much like complex numbers do in two dimensions. With geometric algebra, they are just spinors acting in a different geometry.

Quaternions

Note:

There is support for quaternions in numpy through the package quaternion.

For some reason people think quaternions (wiki page) are mystical or something. They are just spinors in a three dimensional geometric algebra.

In either case, we can pass the names parameters to Cl() to explicitly label the bivectors i,j, and k.

import clifford as cf

# the vector/bivector order is reversed because Hamilton defined quaternions using a
# left-handed frame. wtf.
names = ['', 'z', 'y', 'x', 'k', 'j', 'i', 'I']

layout, blades = cf.Cl(3, names=names)
locals().update(blades)

This leads to the commutations relations familiar to quaternion users

for m in [i, j, k]:
    for n in [i, j, k]:
        print ('{}*{}={}'.format(m, n, m*n))

(1^i)*(1^i)=-1
(1^i)*(1^j)=(1^k)
(1^i)*(1^k)=-(1^j)
(1^j)*(1^i)=-(1^k)
(1^j)*(1^j)=-1
(1^j)*(1^k)=(1^i)
(1^k)*(1^i)=(1^j)
(1^k)*(1^j)=-(1^i)
(1^k)*(1^k)=-1

Quaternion data could be stored in a variety of ways. Assuming you have the scalar components for the quaternion, all you will need to do is setup a map each component to the correct bivector.

def q2S(*args):
    '''convert tuple of quaternion coefficients to a spinor'''
    q = args
    return q[0] + q[1]*i + q[2]*j + q[3]*k

Then all the quaternion computations can be done using GA

q1 = q2S(1,2,3,4)
q1

1 + (4^k) + (3^j) + (2^i)

This prints \(i,j\) and \(k\) in reverse order but whatever,

# 'scalar' part
q1(0)

1

# 'vector' part (more like bivector part!)
q1(2)

(4^k) + (3^j) + (2^i)

quaternion conjugation is implemented with reversion

~q1

1 - (4^k) - (3^j) - (2^i)

The norm

abs(q1)

5.477225575051661

Taking the dual() of the “vector” part actually returns a vector,

q1(2).dual()

(2^z) - (3^y) + (4^x)

q1 = q2S(1, 2, 3, 4)
q2 = q2S(5, 6, 7, 8)

# quaternion product
q1*q2

-60 + (24^k) + (30^j) + (12^i)

If you want to keep using a left-handed frame and names like \(i,j\) and \(k\) to label the planes in 3D space, ok. If you think it makes more sense to use the consistent and transparent approach provided by GA, we think you have good taste. If we make this switch, the basis and q2S() functions will be changed to

import clifford as cf
layout, blades = cf.Cl(3)
locals().update(blades)

blades

{'': 1,
 'e1': (1^e1),
 'e2': (1^e2),
 'e3': (1^e3),
 'e12': (1^e12),
 'e13': (1^e13),
 'e23': (1^e23),
 'e123': (1^e123)}

def q2S(*args):
    '''
    convert tuple of quaternion coefficients to a spinor'''
    q = args
    return q[0] + q[1]*e13 +q[2]*e23 + q[3]*e12

q1 = q2S(1,2,3,4)
q1

1 + (4^e12) + (2^e13) + (3^e23)

This page was generated from docs/tutorials/PerformanceCliffordTutorial.ipynb. Interactive online version: .

Writing high(ish) performance code with Clifford and Numba via Numpy

This document describes how to take algorithms developed in the clifford package with notation that is close to the maths and convert it into numerically efficient and fast running code. To do this we will expose the underlying representation of multivector as a numpy array of canonical basis vector coefficients and operate directly on these arrays in a manner that is conducive to JIT compilation with numba.

First import the Clifford library as well as numpy and numba

import clifford as cf
import numpy as np
import numba

Choose a specific space

For this document we will use 3d euclidean space embedded in the conformal framework giving a Cl(4,1) algebra. We will also rename some of the variables to match the notation that used by Lasenby et al. in “A Covariant Approach to Geometry using Geometric Algebra”

from clifford import g3c
# Get the layout in our local namespace etc etc
layout = g3c.layout
locals().update(g3c.blades)

ep, en, up, down, homo, E0, ninf, no = (g3c.stuff["ep"], g3c.stuff["en"],
                                        g3c.stuff["up"], g3c.stuff["down"], g3c.stuff["homo"],
                                        g3c.stuff["E0"], g3c.stuff["einf"], -g3c.stuff["eo"])
# Define a few useful terms
E = ninf^(no)
I5 = e12345
I3 = e123

Performance of mathematically idiomatic Clifford algorithms

By default the Clifford library sacrifices performance for syntactic convenience.

Consider a function that applies a rotor to a multivector:

def apply_rotor(R,mv):
    return R*mv*~R

We will define a rotor that takes one line to another:

line_one = (up(0)^up(e1)^ninf).normal()
line_two = (up(0)^up(e2)^ninf).normal()
R = 1 + line_two*line_one

Check that this works

print(line_two)
print(apply_rotor(R,line_one).normal())

(1.0^e245)
(1.0^e245)

We would like to improve the speed of our algorithm, first we will profile it and see where it spends its time

#%%prun -s cumtime
#for i in range(1000000):
#    apply_rotor(R,line_one)

An example profile output from running this notebook on the author’s laptop is as follows:

 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      1    0.000    0.000   66.290   66.290 {built-in method builtins.exec}
      1    0.757    0.757   66.290   66.290 <string>:2(<module>)
1000000    3.818    0.000   65.534    0.000 <ipython-input-13-70a01003bf51>:1(apply_rotor)
2000000    9.269    0.000   55.641    0.000 __init__.py:751(__mul__)
2000000    3.167    0.000   29.900    0.000 __init__.py:717(_checkOther)
2000000    1.371    0.000   19.906    0.000 __init__.py:420(__ne__)
2000000    6.000    0.000   18.535    0.000 numeric.py:2565(array_equal)
2000000   10.505    0.000   10.505    0.000 __init__.py:260(mv_mult)

We can see that the function spends almost all of its time in __mul__ and within __mul__ it spends most of its time in _checkOther. In fact it only spends a small fraction of its time in mv_mult which does the numerical multivector multiplication. To write more performant code we need to strip away the high level abstractions and deal with the underlying representations of the blade component data.

Canonical blade coefficient representation in Clifford

In Clifford a multivector is internally represented as a numpy array of the coefficients of the canonical basis vectors, they are arranged in order of grade. So for our 4,1 algebra the first element is the scalar part, the next 5 are the vector coefficients, the next 10 are the bivectors, the next 10 the triectors, the next 5 the quadvectors and the final value is the pseudoscalar coefficient.

(5.0*e1 - e2 + e12 + e135 + np.pi*e1234).value

array([ 0.        ,  5.        , -1.        ,  0.        ,  0.        ,
        0.        ,  1.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  3.14159265,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ])

Exploiting blade representation to write a fast function

We can rewrite our rotor application function using the functions that the layout exposes for operations on the numpy arrays themselves.

def apply_rotor_faster(R,mv):
    return layout.MultiVector(layout.gmt_func(R.value,layout.gmt_func(mv.value,layout.adjoint_func(R.value))) )

#%%prun -s cumtime
#for i in range(1000000):
#    apply_rotor_faster(R,line_one)

This gives a much faster function

 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      1    0.000    0.000   19.567   19.567 {built-in method builtins.exec}
      1    0.631    0.631   19.567   19.567 <string>:2(<module>)
1000000    7.373    0.000   18.936    0.000 <ipython-input-35-6a5344d83bdb>:1(apply_rotor_faster)
2000000    9.125    0.000    9.125    0.000 __init__.py:260(mv_mult)
1000000    1.021    0.000    1.619    0.000 __init__.py:677(__init__)
1000000    0.819    0.000    0.819    0.000 __init__.py:244(adjoint_func)

We have successfully skipped past the higher level checks on the multivectors while maintaining exactly the same function signature. It is important to check that we still have the correct answer:

print(line_two)
print(apply_rotor_faster(R,line_one).normal())

(1.0^e245)
(1.0^e245)

The performance improvements gained by rewriting our function are significant but it comes at the cost of readability.

By loading the layouts gmt_func and adjoint_func into the global namespace before the function is defined and separating the value operations from the multivector wrapper we can make our code more concise.

gmt_func = layout.gmt_func
adjoint_func = layout.adjoint_func

def apply_rotor_val(R_val,mv_val):
    return gmt_func(R_val,gmt_func(mv_val,adjoint_func(R_val)))

def apply_rotor_wrapped(R,mv):
    return cf.MultiVector(layout,apply_rotor_val(R.value,mv.value))

#%%prun -s cumtime
#for i in range(1000000):
#    apply_rotor_wrapped(R,line_one)

The time cost is essentially the same, there is probably some minor overhead from the function call itself

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      1    0.000    0.000   19.621   19.621 {built-in method builtins.exec}
      1    0.557    0.557   19.621   19.621 <string>:2(<module>)
1000000    1.421    0.000   19.064    0.000 <ipython-input-38-a1e0b5c53cdc>:7(apply_rotor_wrapped)
1000000    6.079    0.000   16.033    0.000 <ipython-input-38-a1e0b5c53cdc>:4(apply_rotor_val)
2000000    9.154    0.000    9.154    0.000 __init__.py:260(mv_mult)
1000000    1.017    0.000    1.610    0.000 __init__.py:677(__init__)
1000000    0.800    0.000    0.800    0.000 __init__.py:244(adjoint_func)

print(line_two)
print(apply_rotor_wrapped(R,line_one).normal())

(1.0^e245)
(1.0^e245)

The additional advantage of splitting the function like this is that the numba JIT compiler can reason about the memory layout of numpy arrays in no python mode as long as no pure python objects are operated upon within the function. This means we can JIT our function that operates on the value directly.

@numba.njit
def apply_rotor_val_numba(R_val,mv_val):
    return gmt_func(R_val,gmt_func(mv_val,adjoint_func(R_val)))

def apply_rotor_wrapped_numba(R,mv):
    return cf.MultiVector(layout,apply_rotor_val_numba(R.value,mv.value))

#%%prun -s cumtime
#for i in range(1000000):
#    apply_rotor_wrapped_numba(R,line_one)

This gives a small improvement in performance but more importantly it allows us to write larger functions that also use the jitted apply_rotor_val_numba and are themselves jitted.

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       1    0.000    0.000   16.033   16.033 {built-in method builtins.exec}
       1    0.605    0.605   16.033   16.033 <string>:2(<module>)
 1000000    2.585    0.000   15.428    0.000 <ipython-input-42-1142126d93ca>:5(apply_rotor_wrapped_numba)
 1000000    8.606    0.000    8.606    0.000 <ipython-input-42-1142126d93ca>:1(apply_rotor_val_numba)
       1    0.000    0.000    2.716    2.716 dispatcher.py:294(_compile_for_args)
     7/1    0.000    0.000    2.716    2.716 dispatcher.py:554(compile)

Composing larger functions

By chaining together functions that operate on the value arrays of multivectors it is easy to construct fast and readable code

I5_val = I5.value
omt_func = layout.omt_func

def dual_mv(mv):
    return -I5*mv

def meet_unwrapped(mv_a,mv_b):
    return -dual_mv(dual_mv(mv_a)^dual_mv(mv_b))

@numba.njit
def dual_val(mv_val):
    return -gmt_func(I5_val,mv_val)

@numba.njit
def meet_val(mv_a_val,mv_b_val):
    return -dual_val( omt_func( dual_val(mv_a_val) , dual_val(mv_b_val)) )

def meet_wrapped(mv_a,mv_b):
    return cf.layout.MultiVector(meet_val(mv_a.value, mv_b.value))

sphere = (up(0)^up(e1)^up(e2)^up(e3)).normal()
print(sphere.meet(line_one).normal().normal())
print(meet_unwrapped(sphere,line_one).normal())
print(meet_wrapped(line_one,sphere).normal())

(1.0^e14) - (1.0^e15) - (1.0^e45)
(1.0^e14) - (1.0^e15) - (1.0^e45)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [16], in <module>
     22 print(sphere.meet(line_one).normal().normal())
     23 print(meet_unwrapped(sphere,line_one).normal())
---> 24 print(meet_wrapped(line_one,sphere).normal())

Input In [16], in meet_wrapped(mv_a, mv_b)
     18 def meet_wrapped(mv_a,mv_b):
---> 19     return cf.layout.MultiVector(meet_val(mv_a.value, mv_b.value))

AttributeError: module 'clifford' has no attribute 'layout'

#%%prun -s cumtime
#for i in range(100000):
#    meet_unwrapped(sphere,line_one)

ncalls  tottime  percall  cumtime  percall filename:lineno(function)         1    0.000    0.000   13.216   13.216 {built-in method builtins.exec}         1    0.085    0.085   13.216   13.216 <string>:2(<module>)    100000    0.418    0.000   13.131    0.000 <ipython-input-98-f91457c8741a>:7(meet_unwrapped)    300000    0.681    0.000    9.893    0.000 <ipython-input-98-f91457c8741a>:4(dual_mv)    300000    1.383    0.000    8.127    0.000 __init__.py:751(__mul__)    400000    0.626    0.000    5.762    0.000 __init__.py:717(_checkOther)    400000    0.270    0.000    3.815    0.000 __init__.py:420(__ne__)    400000    1.106    0.000    3.544    0.000 numeric.py:2565(array_equal)    100000    0.460    0.000    2.439    0.000 __init__.py:783(__xor__)    800000    0.662    0.000    2.053    0.000 __init__.py:740(_newMV)    400000    1.815    0.000    1.815    0.000 __init__.py:260(mv_mult)

#%%prun -s cumtime
#for i in range(100000):
#    meet_wrapped(sphere,line_one)

ncalls  tottime  percall  cumtime  percall filename:lineno(function)         1    0.000    0.000    1.951    1.951 {built-in method builtins.exec}         1    0.063    0.063    1.951    1.951 <string>:2(<module>)    100000    0.274    0.000    1.888    0.000 <ipython-input-98-f91457c8741a>:18(meet_wrapped)    100000    1.448    0.000    1.448    0.000 <ipython-input-98-f91457c8741a>:14(meet_val)    100000    0.096    0.000    0.166    0.000 __init__.py:677(__init__)

Algorithms exploiting known sparseness of MultiVector value array

The standard multiplication generator function for two general multivectors is as follows:

def get_mult_function(mult_table,n_dims):
    '''
    Returns a function that implements the mult_table on two input multivectors
    '''
    non_zero_indices = mult_table.nonzero()
    k_list = non_zero_indices[0]
    l_list = non_zero_indices[1]
    m_list = non_zero_indices[2]
    mult_table_vals = np.array([mult_table[k,l,m] for k,l,m in np.transpose(non_zero_indices)],dtype=int)

    @numba.njit
    def mv_mult(value,other_value):
        output = np.zeros(n_dims)
        for ind,k in enumerate(k_list):
            l = l_list[ind]
            m = m_list[ind]
            output[l] += value[k]*mult_table_vals[ind]*other_value[m]
        return output
    return mv_mult

There are however instances in which we might be able to use the known sparseness of the input data value representation to speed up the operations. For example, in Cl(4,1) rotors only contain even grade blades and we can therefore remove all the operations accessing odd grade objects.

def get_grade_from_index(index_in):
    if index_in == 0:
        return 0
    elif index_in < 6:
        return 1
    elif index_in < 16:
        return 2
    elif index_in < 26:
        return 3
    elif index_in < 31:
        return 4
    elif index_in == 31:
        return 5
    else:
        raise ValueError('Index is out of multivector bounds')

def get_sparse_mult_function(mult_table,n_dims,grades_a,grades_b):
    '''
    Returns a function that implements the mult_table on two input multivectors
    '''
    non_zero_indices = mult_table.nonzero()
    k_list = non_zero_indices[0]
    l_list = non_zero_indices[1]
    m_list = non_zero_indices[2]
    mult_table_vals = np.array([mult_table[k,l,m] for k,l,m in np.transpose(non_zero_indices)],dtype=int)

    # Now filter out the sparseness
    filter_mask = np.zeros(len(k_list), dtype=bool)
    for i in range(len(filter_mask)):
        if get_grade_from_index(k_list[i]) in grades_a:
            if get_grade_from_index(m_list[i]) in grades_b:
                filter_mask[i] = 1

    k_list = k_list[filter_mask]
    l_list = l_list[filter_mask]
    m_list = m_list[filter_mask]
    mult_table_vals = mult_table_vals[filter_mask]

    @numba.njit
    def mv_mult(value,other_value):
        output = np.zeros(n_dims)
        for ind,k in enumerate(k_list):
            l = l_list[ind]
            m = m_list[ind]
            output[l] += value[k]*mult_table_vals[ind]*other_value[m]
        return output
    return mv_mult

left_rotor_mult = get_sparse_mult_function(layout.gmt,layout.gaDims,[0,2,4],[0,1,2,3,4,5])
right_rotor_mult = get_sparse_mult_function(layout.gmt,layout.gaDims,[0,1,2,3,4,5],[0,2,4])

@numba.njit
def sparse_apply_rotor_val(R_val,mv_val):
    return left_rotor_mult(R_val,right_rotor_mult(mv_val,adjoint_func(R_val)))

def sparse_apply_rotor(R,mv):
    return cf.MultiVector(layout,sparse_apply_rotor_val(R.value,mv.value))

#%%prun -s cumtime
#for i in range(1000000):
#    sparse_apply_rotor(R,line_one)

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       1    0.000    0.000    9.490    9.490 {built-in method builtins.exec}
       1    0.624    0.624    9.489    9.489 <string>:2(<module>)
 1000000    2.684    0.000    8.865    0.000 <ipython-input-146-f75aae3ce595>:8(sparse_apply_rotor)
 1000000    4.651    0.000    4.651    0.000 <ipython-input-146-f75aae3ce595>:4(sparse_apply_rotor_val)
 1000000    0.934    0.000    1.530    0.000 __init__.py:677(__init__)
 1000000    0.596    0.000    0.596    0.000 {built-in method numpy.core.multiarray.array}

print(line_two)
print(sparse_apply_rotor(R,line_one).normal())

(1.0^e245)
(1.0^e245)

We can do the same with the meet operation that we defined earlier if we know what grade objects we are meeting

left_pseudo_mult = get_sparse_mult_function(layout.gmt,layout.gaDims,[5],[0,1,2,3,4,5])
sparse_omt_2_1 = get_sparse_mult_function(layout.omt,layout.gaDims,[2],[1])

@numba.njit
def dual_sparse_val(mv_val):
    return -left_pseudo_mult(I5_val,mv_val)

@numba.njit
def meet_sparse_3_4_val(mv_a_val,mv_b_val):
    return -dual_sparse_val( sparse_omt_2_1( dual_sparse_val(mv_a_val) , dual_sparse_val(mv_b_val)) )

def meet_sparse_3_4(mv_a,mv_b):
    return cf.layout.MultiVector(meet_sparse_3_4_val(mv_a.value, mv_b.value))

print(sphere.meet(line_one).normal().normal())
print(meet_sparse_3_4(line_one,sphere).normal())

(1.0^e14) - (1.0^e15) - (1.0^e45)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [24], in <module>
     13     return cf.layout.MultiVector(meet_sparse_3_4_val(mv_a.value, mv_b.value))
     15 print(sphere.meet(line_one).normal().normal())
---> 16 print(meet_sparse_3_4(line_one,sphere).normal())

Input In [24], in meet_sparse_3_4(mv_a, mv_b)
     12 def meet_sparse_3_4(mv_a,mv_b):
---> 13     return cf.layout.MultiVector(meet_sparse_3_4_val(mv_a.value, mv_b.value))

AttributeError: module 'clifford' has no attribute 'layout'

#%%prun -s cumtime
#for i in range(100000):
#    meet_sparse_3_4(line_one,sphere)

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     1    0.000    0.000    0.725    0.725 {built-in method builtins.exec}
     1    0.058    0.058    0.725    0.725 <string>:2(<module>)
100000    0.252    0.000    0.667    0.000 <ipython-input-156-f346d0563682>:12(meet_sparse_3_4)
100000    0.267    0.000    0.267    0.000 <ipython-input-156-f346d0563682>:8(meet_sparse_3_4_val)
100000    0.088    0.000    0.148    0.000 __init__.py:677(__init__)

Future work on performance

Investigate efficient operations on containers of large numbers of multivectors.

Possibly investigate http://numba.pydata.org/numba-doc/0.13/CUDAJit.html for larger algebras/other areas in which GPU memory latency will not be such a large factor, ie, lots of bulk parallel numerical operations.

This page was generated from docs/tutorials/cga/index.ipynb. Interactive online version: .

Conformal Geometric Algebra

Conformal Geometric Algebra (CGA) is a projective geometry tool which allows conformal transformations to be implemented with rotations. To do this, the original geometric algebra is extended by two dimensions, one of positive signature \(e_+\) and one of negative signature \(e_-\). Thus, if we started with \(G_p\), the conformal algebra is \(G_{p+1,1}\).

It is convenient to define a null basis given by

\[\begin{split}e_{o} = \frac{1}{2}(e_{-} -e_{+})\\e_{\infty} = e_{-}+e_{+}\end{split}\]

A vector in the original space \(x\) is up-projected into a conformal vector \(X\) by

\[X = x + \frac{1}{2} x^2 e_{\infty} +e_o\]

To map a conformal vector back into a vector from the original space, the vector is first normalized, then rejected from the minkowski plane \(E_0\),

\[X = \frac{X}{X \cdot e_{\infty}}\]

then

\[x = X \wedge E_0\, E_0^{-1}\]

To implement this in clifford we could create a CGA by instantiating the it directly, like Cl(3,1) for example, and then making the definitions and maps described above relating the various subspaces. Or, we you can use the helper function conformalize().

Using conformalize()

The purpose of conformalize() is to remove the redundancy associated with creating a conformal geometric algebras. conformalize() takes an existing geometric algebra layout and conformalizes it by adding two dimensions, as described above. Additionally, this function returns a new layout for the CGA, a dict of blades for the CGA, and dictionary containing the added basis vectors and up/down projection functions.

To demonstrate we will conformalize \(G_2\), producing a CGA of \(G_{3,1}\).

from numpy import pi,e
from clifford import Cl, conformalize

G2, blades_g2 = Cl(2)

blades_g2 # inspect the G2 blades

{'': 1, 'e1': (1^e1), 'e2': (1^e2), 'e12': (1^e12)}

Now, conformalize it

G2c, blades_g2c, stuff = conformalize(G2)

blades_g2c   #inspect the CGA blades

{'': 1,
 'e1': (1^e1),
 'e2': (1^e2),
 'e3': (1^e3),
 'e4': (1^e4),
 'e12': (1^e12),
 'e13': (1^e13),
 'e14': (1^e14),
 'e23': (1^e23),
 'e24': (1^e24),
 'e34': (1^e34),
 'e123': (1^e123),
 'e124': (1^e124),
 'e134': (1^e134),
 'e234': (1^e234),
 'e1234': (1^e1234)}

Additionally lets inspect stuff

stuff

{'ep': (1^e3),
 'en': (1^e4),
 'eo': -(0.5^e3) + (0.5^e4),
 'einf': (1^e3) + (1^e4),
 'E0': (1.0^e34),
 'up': <bound method ConformalLayout.up of ConformalLayout([1, 1, 1, -1], ids=BasisVectorIds.ordered_integers(4), order=BasisBladeOrder.shortlex(4), names=['', 'e1', 'e2', 'e3', 'e4', 'e12', 'e13', 'e14', 'e23', 'e24', 'e34', 'e123', 'e124', 'e134', 'e234', 'e1234'])>,
 'down': <bound method ConformalLayout.down of ConformalLayout([1, 1, 1, -1], ids=BasisVectorIds.ordered_integers(4), order=BasisBladeOrder.shortlex(4), names=['', 'e1', 'e2', 'e3', 'e4', 'e12', 'e13', 'e14', 'e23', 'e24', 'e34', 'e123', 'e124', 'e134', 'e234', 'e1234'])>,
 'homo': <bound method ConformalLayout.homo of ConformalLayout([1, 1, 1, -1], ids=BasisVectorIds.ordered_integers(4), order=BasisBladeOrder.shortlex(4), names=['', 'e1', 'e2', 'e3', 'e4', 'e12', 'e13', 'e14', 'e23', 'e24', 'e34', 'e123', 'e124', 'e134', 'e234', 'e1234'])>,
 'I_base': (1.0^e12)}

It contains the following:

• ep - positive basis vector added

• en - negative basis vector added

• eo - zero vector of null basis (=.5*(en-ep))

• einf - infinity vector of null basis (=en+ep)

• E0 - minkowski bivector (=einf^eo)

• up() - function to up-project a vector from GA to CGA

• down() - function to down-project a vector from CGA to GA

• homo() - function to homogenize a CGA vector

We can put the blades and the stuff into the local namespace,

locals().update(blades_g2c)
locals().update(stuff)

Now we can use the up() and down() functions to go in and out of CGA

x = e1+e2
X = up(x)
X

(1.0^e1) + (1.0^e2) + (0.5^e3) + (1.5^e4)

down(X)

(1.0^e1) + (1.0^e2)

Operations

Conformal transformations in \(G_n\) are achieved through versors in the conformal space \(G_{n+1,1}\). These versors can be categorized by their relation to the added minkowski plane, \(E_0\). There are three categories,

• versor purely in \(E_0\)

• versor partly in \(E_0\)

• versor out of \(E_0\)

A three dimensional projection for conformal space with the relevant subspaces labeled is shown below.

Versors purely in \(E_0\)

First we generate some vectors in G2, which we can operate on

a= 1*e1 + 2*e2
b= 3*e1 + 4*e2

Inversions

\[e_{+} X e_{+}\]

Inversion is a reflection in \(e_+\), this swaps \(e_o\) and \(e_{\infty}\), as can be seen from the model above.

assert(down(ep*up(a)*ep)  == a.inv())

Involutions

\[E_0 X E_0\]

assert(down(E0*up(a)*E0) == -a)

Dilations

\[D_{\alpha} = e^{-\frac{\ln{\alpha}}{2} \,E_0}\]

\[D_{\alpha} \, X \, \tilde{D_{\alpha}}\]

from scipy import rand,log

D = lambda alpha: e**((-log(alpha)/2.)*(E0))
alpha = rand()
assert(down( D(alpha)*up(a)*~D(alpha)) == (alpha*a))

/tmp/ipykernel_695/858833368.py:4: DeprecationWarning: scipy.rand is deprecated and will be removed in SciPy 2.0.0, use numpy.random.rand instead
  alpha = rand()
/tmp/ipykernel_695/858833368.py:3: DeprecationWarning: scipy.log is deprecated and will be removed in SciPy 2.0.0, use numpy.lib.scimath.log instead
  D = lambda alpha: e**((-log(alpha)/2.)*(E0))
/tmp/ipykernel_695/858833368.py:3: DeprecationWarning: scipy.log is deprecated and will be removed in SciPy 2.0.0, use numpy.lib.scimath.log instead
  D = lambda alpha: e**((-log(alpha)/2.)*(E0))

Versors partly in \(E_0\)

Translations

\[V = e ^{\frac{1}{2} e_{\infty} a } = 1 + e_{\infty}a\]

T = lambda x: e**(1/2.*(einf*x))
assert(down( T(a)*up(b)*~T(a)) == b+a)

Transversions

A transversion is an inversion, followed by a translation, followed by a inversion. The versor is

\[V= e_+ T_a e_+\]

which is recognised as the translation bivector reflected in the \(e_+\) vector. From the diagram, it is seen that this is equivalent to the bivector in \(x\wedge e_o\),

\[e_+ (1+e_{\infty}a)e_+\]

\[e_+^2 + e_+e_{\infty}a e_+\]

\[2 +2e_o a\]

the factor of 2 may be dropped, because the conformal vectors are null

V = ep * T(a) * ep
assert ( V == 1+(eo*a))

K = lambda x: 1+(eo*a)

B= up(b)
assert( down(K(a)*B*~K(a)) == 1/(a+1/b) )

Versors Out of \(E_0\)

Versors that are out of \(E_0\) are made up of the versors within the original space. These include reflections and rotations, and their conformal representation is identical to their form in \(G^n\), except the minus sign is dropped for reflections,

Reflections

\[-mam^{-1} \rightarrow MA\tilde{M}\]

m = 5*e1 + 6*e2
n = 7*e1 + 8*e2


assert(down(m*up(a)*m) == -m*a*m.inv())

Rotations

\[mnanm = Ra\tilde{R} \rightarrow RA\tilde{R}\]

R = lambda theta: e**((-.5*theta)*(e12))
theta = pi/2
assert(down( R(theta)*up(a)*~R(theta)) == R(theta)*a*~R(theta))

Combinations of Operations

simple example

As a simple example consider the combination operations of translation,scaling, and inversion.

\[b=-2a+e_0 \quad \rightarrow \quad B= (T_{e_0}E_0 D_2) A \tilde{ (D_2 E_0 T_{e_0})}\]

A = up(a)
V = T(e1)*E0*D(2)
B = V*A*~V
assert(down(B) == (-2*a)+e1 )

/tmp/ipykernel_695/858833368.py:3: DeprecationWarning: scipy.log is deprecated and will be removed in SciPy 2.0.0, use numpy.lib.scimath.log instead
  D = lambda alpha: e**((-log(alpha)/2.)*(E0))

Transversion

A transversion may be built from a inversion, translation, and inversion.

\[c = (a^{-1}+b)^{-1}\]

In conformal GA, this is accomplished by

\[C = VA\tilde{V}\]

\[V= e_+ T_b e_+\]

A = up(a)
V = ep*T(b)*ep
C = V*A*~V
assert(down(C) ==1/(1/a +b))

Rotation about a point

Rotation about a point \(a\) can be achieved by translating the origin to \(a\), then rotating, then translating back. Just like the transversion can be thought of as translating the involution operator, rotation about a point can also be thought of as translating the Rotor itself. Covariance.

More examples

This page was generated from docs/tutorials/cga/visualization-tools.ipynb. Interactive online version: .

Visualization tools

In this example we will look at some external tools that can be used with clifford to help visualize geometric objects.

The two tools available are:

• pyganja (github)

• mpl_toolkits.clifford (github)

Both of these can be installed with pip install followed by the package name above.

G2C

Let’s start by creating some objects in 2d Conformal Geometric Algebra to visualize:

from clifford.g2c import *

point = up(2*e1+e2)
line = up(3*e1 + 2*e2) ^ up(3*e1 - 2*e2) ^ einf
circle = up(e1) ^ up(-e1 + 2*e2) ^ up(-e1 - 2*e2)

We’ll create copies of the point and line reflected in the circle, using \(X = C\hat X\tilde C\), where \(\hat X\) is the grade involution.

point_refl = circle * point.gradeInvol() * ~circle
line_refl = circle * line.gradeInvol() * ~circle

pyganja

pyganja is a python interface to the ganja.js (github) library. To use it, typically we need to import two names from the library:

from pyganja import GanjaScene, draw
import pyganja; pyganja.__version__

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

'0.0.15'

GanjaScene lets us build scenes out of geometric objects, with attached labels and RGB colors:

sc = GanjaScene()
sc.add_object(point, color=(255, 0, 0), label='point')
sc.add_object(line, color=(0, 255, 0), label='line')
sc.add_object(circle, color=(0, 0, 255), label='circle')

sc_refl = GanjaScene()
sc_refl.add_object(point_refl, color=(128, 0, 0), label='point_refl')
sc_refl.add_object(line_refl, color=(0, 128, 0), label='line_refl')

Once we’ve built our scene, we can draw it, specifying a scale (which here we use to zoom out), and the signature of our algebra (which defaults to conformal 3D):

draw(sc, sig=layout.sig, scale=0.5)

A cool feature of GanjaScene is the ability to use + to draw both scenes together:

draw(sc + sc_refl, sig=layout.sig, scale=0.5)

mpl_toolkits.clifford

While ganja.js produces great diagrams, it’s hard to combine them with other plotting tools. mpl_toolkits.clifford works within matplotlib.

from matplotlib import pyplot as plt
plt.ioff()  # we'll ask for plotting when we want it

# if you're editing this locally, you'll get an interactive UI if you uncomment the following
#
#    %matplotlib notebook

from mpl_toolkits.clifford import plot
import mpl_toolkits.clifford; mpl_toolkits.clifford.__version__

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/traitlets/traitlets.py:3044: FutureWarning: --rc={'figure.dpi': 96, 'savefig.transparent': True} for dict-traits is deprecated in traitlets 5.0. You can pass --rc <key=value> ... multiple times to add items to a dict.
  warn(

'0.0.3'

Assembling the plot is a lot more work, but we also get much more control:

# standard matplotlib stuff - construct empty plots side-by-side, and set the scaling
fig, (ax_before, ax_both) = plt.subplots(1, 2, sharex=True, sharey=True)
ax_before.set(xlim=[-4, 4], ylim=[-4, 4], aspect='equal')
ax_both.set(xlim=[-4, 4], ylim=[-4, 4], aspect='equal')

# plot the objects before reflection on both plots
for ax in (ax_before, ax_both):
    plot(ax, [point], color='tab:blue', label='point', marker='x', linestyle=' ')
    plot(ax, [line], color='tab:green', label='line')
    plot(ax, [circle], color='tab:red', label='circle')

# plot the objects after reflection, with thicker lines
plot(ax_both, [point_refl], color='tab:blue', label='point_refl',  marker='x', linestyle=' ', markeredgewidth=2)
plot(ax_both, [line_refl], color='tab:green', label='line_refl', linewidth=2)

fig.tight_layout()
ax_both.legend()

# show the figure
fig

G3C

Let’s repeat the above, but with 3D Conformal Geometric Algebra. Note that if you’re viewing these docs in a jupyter notebook, the lines below will replace all your 2d variables with 3d ones

from clifford.g3c import *

point = up(2*e1+e2)
line = up(3*e1 + 2*e2) ^ up(3*e1 - 2*e2) ^ einf
circle = up(e1) ^ up(-e1 + 1.6*e2 + 1.2*e3) ^ up(-e1 - 1.6*e2 - 1.2*e3)
sphere = up(3*e1) ^ up(e1) ^ up(2*e1 + e2) ^ up(2*e1 + e3)

# note that due to floating point rounding, we need to truncate back to a single grade here, with ``(grade)``
point_refl = homo((circle * point.gradeInvol() * ~circle)(1))
line_refl = (circle * line.gradeInvol() * ~circle)(3)
sphere_refl = (circle * sphere.gradeInvol() * ~circle)(4)

pyganja

Once again, we can create a pair of scenes exactly as before

sc = GanjaScene()
sc.add_object(point, color=(255, 0, 0), label='point')
sc.add_object(line, color=(0, 255, 0), label='line')
sc.add_object(circle, color=(0, 0, 255), label='circle')
sc.add_object(sphere, color=(0, 255, 255), label='sphere')

sc_refl = GanjaScene()
sc_refl.add_object(point_refl, color=(128, 0, 0), label='point_refl')
sc_refl.add_object(line_refl.normal(), color=(0, 128, 0), label='line_refl')
sc_refl.add_object(sphere_refl.normal(), color=(0, 128, 128), label='sphere_refl')

But this time, when we draw them we don’t need to pass sig. Better yet, we can rotate the 3D world around using left click, pan with right click, and zoom with the scroll wheel.

draw(sc + sc_refl, scale=0.5)

Some more example of using pyganja to visualize 3D CGA can be found in the interpolation and clustering notebooks.

mpl_toolkits.clifford

The 3D approach for matplotlib is much the same. Note that due to poor handling of rounding errors in clifford.tools.classify, a call to .normal() is needed. Along with explicit grade selection, this is a useful trick to try and get something to render which otherwise would not.

# standard matplotlib stuff - construct empty plots side-by-side, and set the scaling
fig, (ax_before, ax_both) = plt.subplots(1, 2, subplot_kw=dict(projection='3d'), figsize=(8, 4))
ax_before.set(xlim=[-4, 4], ylim=[-4, 4], zlim=[-4, 4])
ax_both.set(xlim=[-4, 4], ylim=[-4, 4], zlim=[-4, 4])

# plot the objects before reflection on both plots
for ax in (ax_before, ax_both):
    plot(ax, [point], color='tab:red', label='point', marker='x', linestyle=' ')
    plot(ax, [line], color='tab:green', label='line')
    plot(ax, [circle], color='tab:blue', label='circle')
    plot(ax, [sphere], color='tab:cyan')  # labels do not work for spheres: pygae/mpl_toolkits.clifford#5

# plot the objects after reflection
plot(ax_both, [point_refl], color='tab:red', label='point_refl', marker='x', linestyle=' ', markeredgewidth=2)
plot(ax_both, [line_refl.normal()], color='tab:green', label='line_refl', linewidth=2)
plot(ax_both, [sphere_refl], color='tab:cyan')

fig.tight_layout()
ax_both.legend()

# show the figure
fig

Some more example of using mpl_toolkits.clifford to visualize 3D CGA can be found in the examples folder of the mpl_toolkits.clifford repositiory, here.

This page was generated from docs/tutorials/cga/object-oriented.ipynb. Interactive online version: .

Object Oriented CGA

This is a shelled out demo for a object-oriented approach to CGA with clifford. The CGA object holds the original layout for an arbitrary geometric algebra , and the conformalized version. It provides up/down projections, as well as easy ways to generate objects and operators.

Quick Use Demo

from clifford.cga import CGA, Round, Translation
from clifford import Cl

g3,blades = Cl(3)

cga = CGA(g3)  # make cga from existing ga
# or
cga = CGA(3)   # generate cga from dimension of 'base space'

locals().update(cga.blades)        # put ga's blades in local namespace

C = cga.round(e1,e2,e3,-e2)      # generate unit sphere from points
C

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

Sphere

## Objects
cga.round()              # from None
cga.round(3)             # from dim of space
cga.round(e1,e2,e3,-e2)  # from points
cga.round(e1,e2,e3)      # from points
cga.round(e1,e2)         # from points
cga.round((e1,3))        # from center, radius
cga.round(cga.round().mv)# from existing multivector

cga.flat()               # from None
cga.flat(2)              # from dim of space
cga.flat(e1,e2)          # from points
cga.flat(cga.flat().mv)  # from existing multivector


## Operations
cga.dilation()          # from from None
cga.dilation(.4)        # from int

cga.translation()       # from None
cga.translation(e1+e2)  # from vector
cga.translation(cga.down(cga.null_vector()))

cga.rotation()          # from None
cga.rotation(e12+e23)   # from bivector

cga.transversion(e1+e2).mv

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/clifford/_multivector.py:281: RuntimeWarning: divide by zero encountered in true_divide
  newValue = self.value / other
/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/clifford/_multivector.py:281: RuntimeWarning: invalid value encountered in true_divide
  newValue = self.value / other

1.0 + (0.5^e14) - (0.5^e15) + (0.5^e24) - (0.5^e25)

cga.round().inverted()

(0.74226^e1234) - (1.21395^e1235) + (0.19913^e1245) + (0.17632^e1345) - (0.08081^e2345)

D = cga.dilation(5)
cga.down(D(e1))

(5.0^e1)

C.mv # any CGA object/operator has a multivector

(1.0^e1235)

C.center_down,C.radius # some properties of spheres

(0, 1.0)

T = cga.translation(e1+e2) # make a translation
C_ = T(C)                  # translate the sphere
cga.down(C_.center)        # compute center again

(1.0^e1) + (1.0^e2)

cga.round()       #  no args == random sphere
cga.translation() #             random translation

Translation

if 1 in map(int, [1,2]):
    print(3)

3

Objects

Vectors

a = cga.base_vector()  # random vector with components in base space only
a

-(0.74713^e1) + (0.19482^e2) - (0.63152^e3)

cga.up(a)

-(0.74713^e1) + (0.19482^e2) - (0.63152^e3) - (0.00251^e4) + (0.99749^e5)

cga.null_vector()  # create null vector directly

(0.02501^e1) + (1.41201^e2) + (1.47494^e3) + (1.58494^e4) + (2.58494^e5)

Sphere (point pair, circles)

C = cga.round(e1, e2, -e1, e3) # generates sphere from points
C = cga.round(e1, e2, -e1)     # generates circle from points
C = cga.round(e1, e2)          # generates point-pair from points
#or
C2 = cga.round(2)            # random 2-sphere  (sphere)
C1 = cga.round(1)            # random 1-sphere, (circle)
C0 = cga.round(0)            # random 0-sphere, (point pair)


C1.mv                        # access the multivector

(0.28986^e123) + (0.83473^e124) + (1.26772^e125) - (0.22725^e134) - (0.37397^e135) - (0.08308^e145) + (0.17663^e234) + (0.17591^e235) - (0.26592^e245) + (0.08997^e345)

C = cga.round(e1, e2, -e1, e3)
C.center,C.radius        # spheres have properties

(-(1.0^e4) + (1.0^e5), 1.0)

cga.down(C.center) == C.center_down

True

C_ = cga.round().from_center_radius(C.center,C.radius)
C_.center,C_.radius

(-(2.0^e4) + (2.0^e5), 0.7071067811865476)

Operators

T = cga.translation(e1) # generate translation
T.mv

1.0 - (0.5^e14) - (0.5^e15)

C = cga.round(e1, e2, -e1)
T.mv*C.mv*~T.mv         # translate a sphere

-(0.5^e124) + (0.5^e125) - (1.0^e245)

T(C)                # shorthand call, same as above. returns type of arg

Circle

T(C).center

(2.0^e1) + (2.0^e5)

This page was generated from docs/tutorials/cga/interpolation.ipynb. Interactive online version: .

Example 1 Interpolating Conformal Objects

In this example we will look at a few of the tools provided by the clifford package for (4,1) conformal geometric algebra (CGA) and see how we can use them in a practical setting to interpolate geometric primitives.

The first step in using the package for CGA is to generate and import the algebra:

from clifford.g3c import *
blades

{'': 1,
 'e1': (1^e1),
 'e2': (1^e2),
 'e3': (1^e3),
 'e4': (1^e4),
 'e5': (1^e5),
 'e12': (1^e12),
 'e13': (1^e13),
 'e14': (1^e14),
 'e15': (1^e15),
 'e23': (1^e23),
 'e24': (1^e24),
 'e25': (1^e25),
 'e34': (1^e34),
 'e35': (1^e35),
 'e45': (1^e45),
 'e123': (1^e123),
 'e124': (1^e124),
 'e125': (1^e125),
 'e134': (1^e134),
 'e135': (1^e135),
 'e145': (1^e145),
 'e234': (1^e234),
 'e235': (1^e235),
 'e245': (1^e245),
 'e345': (1^e345),
 'e1234': (1^e1234),
 'e1235': (1^e1235),
 'e1245': (1^e1245),
 'e1345': (1^e1345),
 'e2345': (1^e2345),
 'e12345': (1^e12345)}

This creates an algebra with the required signature and imports the basis blades into the current workspace. We can check our metric by squaring our grade 1 multivectors.

print('e1*e1 ', e1*e1)
print('e2*e2 ', e2*e2)
print('e3*e3 ', e3*e3)
print('e4*e4 ', e4*e4)
print('e5*e5 ', e5*e5)

e1*e1  1
e2*e2  1
e3*e3  1
e4*e4  1
e5*e5  -1

As expected this gives us 4 basis vectors that square to 1.0 and one that squares to -1.0, therefore confirming our metric is (4,1).

The up() function implements the mapping of vectors from standard 3D space to conformal space. We can use this to construct conformal objects to play around with.

For example a line at the origin:

line_a = ( up(0)^up(e1+e2)^einf ).normal()
print(line_a)

(0.70711^e145) + (0.70711^e245)

The tools submodule of the clifford package contains a wide array of algorithms and tools that can be useful for manipulating objects in CGA. We will use these tools to generate rotors that rotate and translate objects:

from clifford.tools.g3 import *
from clifford.tools.g3c import *
from numpy import pi

rotation_radians = pi/4
euc_vector_in_plane_m = e1
euc_vector_in_plane_n = e3

euc_translation = -5.2*e1 + 3*e2 - pi*e3

rotor_rotation = generate_rotation_rotor(rotation_radians, euc_vector_in_plane_m, euc_vector_in_plane_n)
rotor_translation = generate_translation_rotor(euc_translation)
print(rotor_rotation)
print(rotor_translation)

combined_rotor = (rotor_translation*rotor_rotation).normal()

line_b = (combined_rotor*line_a*~combined_rotor).normal()
print(line_b)

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

0.92388 - (0.38268^e13)
1.0 + (2.6^e14) + (2.6^e15) - (1.5^e24) - (1.5^e25) + (1.5708^e34) + (1.5708^e35)
-(5.17696^e124) - (5.17696^e125) - (1.0292^e134) - (1.0292^e135) + (0.5^e145) + (3.72144^e234) + (3.72144^e235) + (0.70711^e245) + (0.5^e345)

In the above snippet of code we have generated rotors for translation and rotation, then combined these, then applied the combined rotor to the original line that we made.

Visualizations

The clifford package can be used alongside pyganja to render CGA objects which can be rotated interactively in a jupyter notebook:

from pyganja import GanjaScene, draw
sc = GanjaScene()
sc.add_object(line_a,color=0xFF0000, label='a')
sc.add_object(line_b,color=0x00FF00, label='b')
draw(sc, scale=0.1)

We can also interpolate the objects using the tools in clifford and can visualise the result

def interpolate_objects_linearly(L1, L2, n_steps=10, color_1=np.array([255,0,0]), color_2=np.array([0,255,0])):
    alpha_list = np.linspace(0, 1, num=n_steps)
    intermediary_list = []
    sc = GanjaScene()
    for alpha in alpha_list:
        intermediate_color = (alpha*color_1 + (1-alpha)*color_2).astype(np.uint8)
        intermediate_object = interp_objects_root(L1, L2, alpha)
        intermediary_list.append(intermediate_object)
        color_string = int(
            (intermediate_color[0] << 16) | (intermediate_color[1] << 8) | intermediate_color[2]
        )
        sc.add_object(intermediate_object, color_string)
    return intermediary_list, sc

intermediary_list, finished_scene = interpolate_objects_linearly(line_a, line_b)
draw(finished_scene, scale=0.1)

We can do the same for all the other geometric primitives as well

Circles

circle_a = (up(0)^up(e1)^up(e2)).normal()
circle_b = (combined_rotor*circle_a*~combined_rotor).normal()
intermediary_list, finished_scene = interpolate_objects_linearly(circle_a, circle_b)
draw(finished_scene, scale=0.1)

Point pairs

point_pair_a = (up(e3)^up(e1+e2)).normal()
point_pair_b = (combined_rotor*point_pair_a*~combined_rotor).normal()
intermediary_list, finished_scene = interpolate_objects_linearly(point_pair_a, point_pair_b)
draw(finished_scene, scale=0.1)

Planes

plane_a = (up(0)^up(e1)^up(e2)^einf).normal()
plane_b = (combined_rotor*plane_a*~combined_rotor).normal()
intermediary_list, finished_scene = interpolate_objects_linearly(plane_a, plane_b)
draw(finished_scene)

Spheres

sphere_a = (up(0)^up(e1)^up(e2)^up(e3)).normal()
sphere_b = (combined_rotor*sphere_a*~combined_rotor).normal()
intermediary_list, finished_scene = interpolate_objects_linearly(sphere_a, sphere_b)
draw(finished_scene, scale=0.1)

This page was generated from docs/tutorials/cga/clustering.ipynb. Interactive online version: .

Example 2 Clustering Geometric Objects

In this example we will look at a few of the tools provided by the clifford package for (4,1) conformal geometric algebra (CGA) and see how we can use them in a practical setting to cluster geometric objects via the simple K-means clustering algorithm provided in clifford.tools

As before the first step in using the package for CGA is to generate and import the algebra:

from clifford.g3c import *
print('e1*e1 ', e1*e1)
print('e2*e2 ', e2*e2)
print('e3*e3 ', e3*e3)
print('e4*e4 ', e4*e4)
print('e5*e5 ', e5*e5)

e1*e1  1
e2*e2  1
e3*e3  1
e4*e4  1
e5*e5  -1

The tools submodule of the clifford package contains a wide array of algorithms and tools that can be useful for manipulating objects in CGA. In this case we will be generating a large number of objects and then segmenting them into clusters.

We first need an algorithm for generating a cluster of objects in space. We will construct this cluster by generating a random object and then repeatedly disturbing this object by some small fixed amount and storing the result:

from clifford.tools.g3c import *
import numpy as np

def generate_random_object_cluster(n_objects, object_generator, max_cluster_trans=1.0, max_cluster_rot=np.pi/8):
    """ Creates a cluster of random objects """
    ref_obj = object_generator()
    cluster_objects = []
    for i in range(n_objects):
        r = random_rotation_translation_rotor(maximum_translation=max_cluster_trans, maximum_angle=max_cluster_rot)
        new_obj = apply_rotor(ref_obj, r)
        cluster_objects.append(new_obj)
    return cluster_objects

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

We can use this function to create a cluster and then we can visualise this cluster with pyganja.

from pyganja import *
clustered_circles = generate_random_object_cluster(10, random_circle)
sc = GanjaScene()
for c in clustered_circles:
    sc.add_object(c, rgb2hex([255,0,0]))
draw(sc, scale=0.05)

This cluster generation function appears in clifford tools by default and it can be imported as follows:

from clifford.tools.g3c import generate_random_object_cluster

Now that we can generate individual clusters we would like to generate many:

def generate_n_clusters( object_generator, n_clusters, n_objects_per_cluster ):
    object_clusters = []
    for i in range(n_clusters):
        cluster_objects = generate_random_object_cluster(n_objects_per_cluster, object_generator,
                                                         max_cluster_trans=0.5, max_cluster_rot=np.pi / 16)
        object_clusters.append(cluster_objects)
    all_objects = [item for sublist in object_clusters for item in sublist]
    return all_objects, object_clusters

Again this function appears by default in clifford tools and we can easily visualise the result:

from clifford.tools.g3c import generate_n_clusters

all_objects, object_clusters = generate_n_clusters(random_circle, 2, 5)
sc = GanjaScene()
for c in all_objects:
    sc.add_object(c, rgb2hex([255,0,0]))
draw(sc, scale=0.05)

Given that we can now generate multiple clusters of objects we can test algorithms for segmenting them.

The function run_n_clusters below generates a lot of objects distributed into n clusters and then attempts to segment the objects to recover the clusters.

from clifford.tools.g3c.object_clustering import n_clusters_objects
import time

def run_n_clusters( object_generator, n_clusters, n_objects_per_cluster, n_shotgunning):
    all_objects, object_clusters = generate_n_clusters( object_generator, n_clusters, n_objects_per_cluster )
    [new_labels, centroids, start_labels, start_centroids] = n_clusters_objects(n_clusters, all_objects,
                                                                                initial_centroids=None,
                                                                                n_shotgunning=n_shotgunning,
                                                                                averaging_method='unweighted')
    return all_objects, new_labels, centroids

Lets try it!

def visualise_n_clusters(all_objects, centroids, labels,
                         color_1=np.array([255, 0, 0]),
                         color_2=np.array([0, 255, 0])):
    """
    Utility method for visualising several clusters and their respective centroids
    using pyganja
    """
    alpha_list = np.linspace(0, 1, num=len(centroids))
    sc = GanjaScene()
    for ind, this_obj in enumerate(all_objects):
        alpha = alpha_list[labels[ind]]
        cluster_color = (alpha * color_1 + (1 - alpha) * color_2).astype(np.int)
        sc.add_object(this_obj, rgb2hex(cluster_color))

    for c in centroids:
        sc.add_object(c, Color.BLACK)

    return sc



object_generator = random_circle

n_clusters = 3
n_objects_per_cluster = 10
n_shotgunning = 60
all_objects, labels, centroids = run_n_clusters(object_generator, n_clusters,
                                                     n_objects_per_cluster, n_shotgunning)

sc = visualise_n_clusters(all_objects, centroids, labels,
                          color_1=np.array([255, 0, 0]),
                          color_2=np.array([0, 255, 0]))
draw(sc, scale=0.05)

/tmp/ipykernel_651/1063371436.py:12: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  cluster_color = (alpha * color_1 + (1 - alpha) * color_2).astype(np.int)

This page was generated from docs/tutorials/cga/robotic-manipulators.ipynb. Interactive online version: .

Application to Robotic Manipulators

This notebook is intended to expand upon the ideas in part of the presentation Robots, Ganja & Screw Theory

Serial manipulator

(slides)

Let’s consider a 2-link 3 DOF arm. We’ll model the links within the robot with rotors, which transform to the coordinate frame of the end of each link. This is very similar to the approach that would classically be taken with 4×4 matrices.

We’re going to define our class piecewise as we go along here. To aid that, we’ll write a simple base class to let us do just that. In your own code, there’s no need to do this.

class AddMethodsAsWeGo:
    @classmethod
    def _add_method(cls, m):
        if isinstance(m, property):
            name = (m.fget or m.fset).__name__
        else:
            name = m.__name__
        setattr(cls, name, m)

Let’s start by defining some names for the links, and a place to store our parameters:

from enum import Enum

class Links(Enum):
    BASE = 'b'
    SHOULDER = 's'
    UPPER = 'u'
    ELBOW = 'e'
    FOREARM = 'f'
    ENDPOINT = 'n'


class SerialRobot(AddMethodsAsWeGo):
    def __init__(self, rho, l):
        self.l = l
        self.rho = rho
        self._thetas = (0, 0, 0)

    @property
    def thetas(self):
        return self._thetas

Forward kinematics

(slides)

As a reminder, we can construct rotation and translation motors as:

\[\begin{split}\begin{align} T(a) &= \exp \left(\frac{1}{2} n_{\infty} \wedge a \right) \\ &= 1 + \frac{1}{2}n_{\infty} \wedge a \\ R(\theta, \hat B) &= \exp (\frac{1}{2} \theta \hat B) \\ &= \cos \frac{\theta}{2} + \sin \frac{\theta}{2} \hat B \end{align}\end{split}\]

Applying these to our geometry, we get

\[\begin{split}\begin{align} R_{\text{base} \gets \text{shoulder}} &= R(\theta_0, e_1 \wedge e_3) \\ R_{\text{shoulder} \gets \text{upper arm}} &= R(\theta_1, e_1 \wedge e_2) \\ R_{\text{upper arm} \gets \text{elbow}} &= T(\rho e_1) \\ R_{\text{elbow} \gets \text{forearm}} &= R(\theta_2, e_1 \wedge e_2) \\ R_{\text{forearm} \gets \text{endpoint}} &= T(-l e_1)\\ \end{align}\end{split}\]

From which we can get the overall rotor to the frame of the endpoint, and the positions \(X\) and \(Y\):

\[\begin{split}\begin{align} R_{\text{base} \gets \text{elbow}} &= R_{\text{base} \gets \text{shoulder}} R_{\text{shoulder} \gets \text{upper arm}} R_{\text{upper arm} \gets \text{elbow}} \\ X &= R_{\text{base} \gets \text{elbow}} n_0 \tilde R_{\text{base} \gets \text{elbow}} \\ R_{\text{base} \gets \text{endpoint}} &= R_{\text{base} \gets \text{shoulder}} R_{\text{shoulder} \gets \text{upper arm}} R_{\text{upper arm} \gets \text{elbow}} R_{\text{elbow} \gets \text{forearm}} R_{\text{forearm} \gets \text{endpoint}} \\ Y &= R_{\text{base} \gets \text{endpoint}} n_0 \tilde R_{\text{base} \gets \text{endpoint}} \\ \end{align}\end{split}\]

We can write this as:

from clifford.g3c import *
from clifford.tools.g3c import generate_translation_rotor, apply_rotor
from clifford.tools.g3 import generate_rotation_rotor

def _update_chain(rotors, a, b, c):
    rotors[a, c] = rotors[a, b] * rotors[b, c]

@SerialRobot._add_method
@SerialRobot.thetas.setter
def thetas(self, value):
    theta0, theta1, theta2 = self._thetas = value
    # shorthands for brevity
    R = generate_rotation_rotor
    T = generate_translation_rotor

    rotors = {}
    rotors[Links.BASE, Links.SHOULDER] = R(theta0, e1, e3)
    rotors[Links.SHOULDER, Links.UPPER] = R(theta1, e1, e2)
    rotors[Links.UPPER, Links.ELBOW] = T(self.rho * e1)
    rotors[Links.ELBOW, Links.FOREARM] = R(theta2, e1, e2)
    rotors[Links.FOREARM, Links.ENDPOINT] = T(-self.l * e1)

    _update_chain(rotors, Links.BASE, Links.SHOULDER, Links.UPPER)
    _update_chain(rotors, Links.BASE, Links.UPPER, Links.ELBOW)
    _update_chain(rotors, Links.BASE, Links.ELBOW, Links.FOREARM)
    _update_chain(rotors, Links.BASE, Links.FOREARM, Links.ENDPOINT)
    self.rotors = rotors

@SerialRobot._add_method
@property
def y_pos(self):
    return apply_rotor(eo, self.rotors[Links.BASE, Links.ENDPOINT])


@SerialRobot._add_method
@property
def x_pos(self):
    return apply_rotor(eo, self.rotors[Links.BASE, Links.ELBOW])

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

Let’s write a renderer so we can check this all works

from pyganja import GanjaScene

def add_rotor(sc: GanjaScene, r, *, label=None, color=None, scale=0.1):
    """ show how a rotor transforms the axes at the origin """
    y = apply_rotor(eo, r)
    y_frame = [
        apply_rotor(d, r)
        for d in [up(scale*e1), up(scale*e2), up(scale*e3)]
    ]
    sc.add_object(y, label=label, color=color)
    sc.add_facet([y, y_frame[0]], color=(255, 0, 0))
    sc.add_facet([y, y_frame[1]], color=(0, 255, 0))
    sc.add_facet([y, y_frame[2]], color=(0, 0, 255))


@SerialRobot._add_method
def to_scene(self):
    sc = GanjaScene()
    axis_scale = 0.1
    link_scale = 0.05
    arm_color = (192, 192, 192)

    base_obj = (up(0.2*e1)^up(0.2*e3)^up(-0.2*e1)).normal()
    sc.add_object(base_obj, color=0)

    shoulder_axis = [
        apply_rotor(p, self.rotors[Links.BASE, Links.UPPER])
        for p in [up(axis_scale*e3), up(-axis_scale*e3)]
    ]
    sc.add_facet(shoulder_axis, color=(0, 0, 128))
    shoulder_angle = [
        apply_rotor(eo, self.rotors[Links.BASE, Links.SHOULDER]),
        apply_rotor(up(axis_scale*e1), self.rotors[Links.BASE, Links.SHOULDER]),
        apply_rotor(up(axis_scale*e1), self.rotors[Links.BASE, Links.UPPER]),
    ]
    sc.add_facet(shoulder_angle, color=(0, 0, 128))

    upper_arm_points = [
        apply_rotor(up(link_scale*e3), self.rotors[Links.BASE, Links.UPPER]),
        apply_rotor(up(-link_scale*e3), self.rotors[Links.BASE, Links.UPPER]),
        apply_rotor(up(link_scale*e3), self.rotors[Links.BASE, Links.ELBOW]),
        apply_rotor(up(-link_scale*e3), self.rotors[Links.BASE, Links.ELBOW])
    ]
    sc.add_facet(upper_arm_points[:3], color=arm_color)
    sc.add_facet(upper_arm_points[1:], color=arm_color)

    elbow_axis = [
        apply_rotor(p, self.rotors[Links.BASE, Links.ELBOW])
        for p in [up(axis_scale*e3), up(-axis_scale*e3)]
    ]
    sc.add_facet(elbow_axis, color=(0, 0, 128))

    forearm_points = [
        apply_rotor(up(link_scale*e3), self.rotors[Links.BASE, Links.FOREARM]),
        apply_rotor(up(-link_scale*e3), self.rotors[Links.BASE, Links.FOREARM]),
        apply_rotor(up(link_scale*e3), self.rotors[Links.BASE, Links.ENDPOINT]),
        apply_rotor(up(-link_scale*e3), self.rotors[Links.BASE, Links.ENDPOINT])
    ]
    sc.add_facet(forearm_points[:3], color=arm_color)
    sc.add_facet(forearm_points[1:], color=arm_color)

    add_rotor(sc, self.rotors[Links.BASE, Links.ELBOW], label='x', color=(128, 128, 128))
    add_rotor(sc, self.rotors[Links.BASE, Links.ENDPOINT], label='y', color=(128, 128, 128))

    return sc

We can now instantiate our robot

serial_robot = SerialRobot(rho=1, l=0.5)

Choose a trajectory

import math
theta_traj = [
    (math.pi/6 + i*math.pi/12, math.pi/3 - math.pi/12*i, 3*math.pi/4)
    for i in range(3)
]

And plot the robot in each state, using ipywidgets (docs) to let us plot ganja side-by-side. Unfortunately, pyganja provides no mechanism to animate these plots from python.

This will not render side-by-side in the online clifford documentation, but will in a local notebook.

import ipywidgets
from IPython.display import Latex, display
from pyganja import draw

outputs = [
    ipywidgets.Output(layout=ipywidgets.Layout(flex='1'))
    for i in range(len(theta_traj))
]
for output, thetas in zip(outputs, theta_traj):
    with output:
        # interesting part here - run the forward kinematics, print the angles we used
        serial_robot.thetas = thetas
        display(Latex(r"$\theta_i = {:.2f}, {:.2f}, {:.2f}$".format(*thetas)))
        draw(serial_robot.to_scene(), scale=1.5)
ipywidgets.HBox(outputs)

Inverse kinematics

(slides)

For the forward kinematics, we didn’t actually need conformal geometric algebra at all—PGA would have done just fine, as all we needed were rotations and translations. The inverse kinematics of a serial manipulator is where CGA provide some nice tricks.

There are three facts we know about the position \(X\), each of which describes a constraint surface

• \(X\) must lie on a sphere with radius \(l\) centered at \(Y\), which can be written

  \[S^* = Y - \frac{1}{2}l^2n_\infty\]

• \(X\) must lie on a sphere with radius \(\rho\) centered at \(n_o\), which can be written

  \[S_\text{base}^* = n_0 - \frac{1}{2}\rho^2n_\infty\]

• \(X\) must lie on a plane through \(n_o\), \(e_3\), and \(Y\), which can be written

  \[\Pi = n_0\wedge \operatorname{up}(e_3)\wedge Y\wedge n_\infty\]

Note that \(\Pi = 0\) is possible iff \(Y = \operatorname{up}(ke_3)\).

For \(X\) to satisfy all three constraints. we have

\begin{align} S \wedge X = S_\text{base} \wedge X = \Pi \wedge X &= 0 \\ X \wedge (\underbrace{S \vee S_\text{base} \vee \Pi}_P) &= 0 \quad\text{If $\Pi \ne 0$} \\ X \wedge (\underbrace{S \vee S_\text{base}}_C) &= 0 \quad\text{otherwise} \\ \end{align}

By looking at the grade of the term labelled \(P\), we conclude it must be a point-pair—which tells us \(X\) must lie in one of two locations. Similarly, \(C\) must be a circle.

@SerialRobot._add_method
def _get_x_constraints_for(self, Y):
    """ Get the space containing all possible elbow positions """
    # strictly should be undual, but we don't have that in clifford
    S = (Y - 0.5*self.l**2*einf).dual()
    S_base = (eo - 0.5*self.rho**2*einf).dual()
    Pi = eo ^ up(e2) ^ Y ^ einf
    return S, S_base, Pi

@SerialRobot._add_method
def _get_x_positions_for(self, Y):
    """ Get the space containing all possible elbow positions """
    S, S_base, Pi = self._get_x_constraints_for(Y)
    if Pi == 0:
        # any solution on the circle is OK
        return S & S_base
    else:
        # there are just two solutions
        return S & S_base & Pi

From the pointpair \(P\) we can extract the two possible \(X\) locations with:

\[X = \left[1 \pm \frac{P}{\sqrt{P\tilde{P}}}\right](P\cdot n_\infty)\]

To be considered a full solution to the inverse kinematics problem, we need to produce the angles \(\theta_0, \theta_1, \theta_2\). We can do this as follows

@SerialRobot._add_method
@SerialRobot.y_pos.setter
def y_pos(self, Y):
    R = generate_rotation_rotor
    T = generate_translation_rotor

    rotors = {}
    rotors[Links.UPPER, Links.ELBOW] = T(self.rho * e1)
    rotors[Links.FOREARM, Links.ENDPOINT] = T(-self.l * e1)

    x_options = self._get_x_positions_for(Y)
    if x_options.grades == {3}:
        # no need to adjust the base angle
        theta_0 = self.thetas[0]
        rotors[Links.BASE, Links.SHOULDER] = self.rotors[Links.BASE, Links.SHOULDER]
        # remove the rotation from x, intersect it with the plane of the links
        x_options = x_options & (eo ^ up(e3) ^ up(e1) ^ einf)
    else:
        y_down = down(Y)
        theta0 = math.atan2(y_down[(3,)], y_down[(1,)])
        rotors[Links.BASE, Links.SHOULDER] = R(theta0, e1, e3)

        # remove the first rotor from x
        x_options = apply_rotor(x_options, ~rotors[Links.BASE, Links.SHOULDER])

    # project out one end of the point-pair
    x = (1 - x_options.normal()) * (x_options | einf)

    x_down = down(x)
    theta1 = math.atan2(x_down[(2,)], x_down[(1,)])
    rotors[Links.SHOULDER, Links.UPPER] = R(theta1, e1, e2)

    _update_chain(rotors, Links.BASE, Links.SHOULDER, Links.UPPER)
    _update_chain(rotors, Links.BASE, Links.UPPER, Links.ELBOW)

    # remove the second rotor
    Y = apply_rotor(Y, ~rotors[Links.BASE, Links.ELBOW])
    y_down = down(Y)

    theta2 = math.atan2(-y_down[(2,)], -y_down[(1,)])
    rotors[Links.ELBOW, Links.FOREARM] = R(theta2, e1, e2)
    _update_chain(rotors, Links.BASE, Links.ELBOW, Links.FOREARM)
    _update_chain(rotors, Links.BASE, Links.FOREARM, Links.ENDPOINT)

    self._thetas = (theta0, theta1, theta2)
    self.rotors = rotors

Define a trajectory again, this time with a scene to render it:

y_traj = [
    up(0.3*e3 + 0.8*e2 - 0.25*e1),
    up(0.6*e3 + 0.8*e2),
    up(0.9*e3 + 0.8*e2 + 0.25*e1)
]

expected_scene = GanjaScene()
expected_scene.add_facet(y_traj[0:2], color=(255, 128, 128))
expected_scene.add_facet(y_traj[1:3], color=(255, 128, 128))

And we can run the inverse kinematics by setting serial_robot.y_pos:

outputs = [
    ipywidgets.Output(layout=ipywidgets.Layout(flex='1'))
    for i in range(len(y_traj))
]
first = True
for output, y in zip(outputs, y_traj):
    with output:
        # interesting part here - run the reverse kinematics, print the angles we used
        serial_robot.y_pos = y
        display(Latex(r"$\theta_i = {:.2f}, {:.2f}, {:.2f}$".format(*serial_robot.thetas)))
        sc = serial_robot.to_scene()

        # Show the spheres we used to construct the solution
        sc += expected_scene
        if first:
            extra_scene = GanjaScene()
            S, S_base, Pi = serial_robot._get_x_constraints_for(y)
            extra_scene.add_object(S_base, label='S_base', color=(255, 255, 128))
            extra_scene.add_object(S, label='S', color=(255, 128, 128))
            extra_scene.add_object(Pi, label='Pi', color=(128, 255, 192, 128))
            sc += extra_scene
        draw(sc, scale=1.5)
    first = False
ipywidgets.HBox(outputs)

Parallel manipulators

For now, refer to the presentation

(slides)

Inverse kinematics

(slides)

For now, refer to the presentation

Forward kinematics

(slides)

For now, refer to the presentation

This page was generated from docs/tutorials/linear-transformations.ipynb. Interactive online version: .

Linear transformations

When working in regular vector spaces, a common tool is a linear transformation, typically in the form of a matrix.

While geometric algebra already provides the rotors as a means of describing transformations (see the CGA tutorial section), there are types of linear transformation that are not suitable for this representation.

This tutorial leans heavily on the explanation of linear transformations in [DFM09], chapter 4. It explores the clifford.transformations submodule.

Vector transformations in linear algebra

As a brief reminder, we can represent transforms in \(\mathbb{R}^3\) using the matrices in \(\mathbb{R}^{3 \times 3}\):

import numpy as np

rot_and_scale_x = np.array([
    [1, 0, 0],
    [0, 1, -1],
    [0, 1, 1],
])

We can read this as a table, where each column corresponds to a component of the input vector, and each row a component of the output:

def show_table(data, cols, rows):
    # trick to get a nice looking table in a notebook
    import pandas as pd; return pd.DataFrame(data, columns=cols, index=rows)

show_table(rot_and_scale_x, ["$\mathit{in}_%s$" % c for c in "xyz"], ["$\mathit{out}_%s$" % c for c in "xyz"])

	$\mathit{in}_x$	$\mathit{in}_y$	$\mathit{in}_z$
$\mathit{out}_x$	1	0	0
$\mathit{out}_y$	0	1	-1
$\mathit{out}_z$	0	1	1

We can apply it to some vectors using the @ matrix multiply operator:

v1 = np.array([1, 0, 0])
v2 = np.array([0, 1, 0])
v3 = np.array([0, 0, 1])

(
    rot_and_scale_x @ v1,
    rot_and_scale_x @ v2,
    rot_and_scale_x @ v3,
)

(array([1, 0, 0]), array([0, 1, 1]), array([ 0, -1,  1]))

We say this transformation is linear because \(f(a + b) = f(a) + f(b)\):

assert np.array_equal(
    rot_and_scale_x @ (2*v1 + 3*v2),
    2 * (rot_and_scale_x @ v1) + 3 * (rot_and_scale_x @ v2)
)

Multivector transformations in geometric algebra

How would we go about applying rot_and_scale_x in a geometric algebra? Clearly we can apply it to vectors in the same way as before, which we can do by unpacking coefficients and repacking them:

from clifford.g3 import *

v = 2*e1 + 3*e2
v_trans = layout.MultiVector()
v_trans[1,], v_trans[2,], v_trans[3,] = rot_and_scale_x @ [v[1,], v[2,], v[3,]]
v_trans

(2.0^e1) + (3.0^e2) + (3.0^e3)

However, in geometric algebra we don’t only care about the vectors, we want to transform the the higher-order blades too. This can be done via an outermorphism, which extends \(f(a)\) to \(f(a \wedge b) = f(a) \wedge f(b)\). This is where the clifford.transformations submodule comes in handy:

from clifford import transformations

rot_and_scale_x_ga = transformations.OutermorphismMatrix(rot_and_scale_x, layout)

To apply these transformations, we use the () operator, rather than @:

rot_and_scale_x_ga(e12)

(1^e12) + (1^e13)

# check it's an outermorphism
rot_and_scale_x_ga(e1) ^ rot_and_scale_x_ga(e2)

(1^e12) + (1^e13)

It shouldn’t come as a surprise that applying the transformation to the psuedoscalar will tell us the determinant of our original matrix - the determinant tells us how a transformation scales volumes, and layout.I is a representation of the unit volume element!

np.linalg.det(rot_and_scale_x), rot_and_scale_x_ga(layout.I)

(2.0, (2^e123))

Matrix representation

Under the hood, clifford implements this using a matrix too - it’s just now a matrix operating over all of the basis blades, not just over the vectors. We can see this by looking at the private _matrix attribute:

show_table(rot_and_scale_x_ga._matrix, ["$\mathit{in}_{%s}$" % c for c in layout.names], ["$\mathit{out}_{%s}$" % c for c in layout.names])

	$\mathit{in}_{}$	$\mathit{in}_{e1}$	$\mathit{in}_{e2}$	$\mathit{in}_{e3}$	$\mathit{in}_{e12}$	$\mathit{in}_{e13}$	$\mathit{in}_{e23}$	$\mathit{in}_{e123}$
$\mathit{out}_{}$	1	0	0	0	0	0	0	0
$\mathit{out}_{e1}$	0	1	0	0	0	0	0	0
$\mathit{out}_{e2}$	0	0	1	-1	0	0	0	0
$\mathit{out}_{e3}$	0	0	1	1	0	0	0	0
$\mathit{out}_{e12}$	0	0	0	0	1	-1	0	0
$\mathit{out}_{e13}$	0	0	0	0	1	1	0	0
$\mathit{out}_{e23}$	0	0	0	0	0	0	2	0
$\mathit{out}_{e123}$	0	0	0	0	0	0	0	2

This page was generated from docs/tutorials/apollonius-cga-augmented.ipynb. Interactive online version: .

Working with custom algebras

This notebook explores the algebra defined in The Lie Model for Euclidean Geometry (Hongbo Li), and its application to solving Apollonius’ Problem. It also shows

The algebra is constructed with basis elements \(e_{-2}, e_{-1}, e_1, \cdots, e_n, e_{n+1}\), where \(e_{-2}^2 = -e_{-1}^2 = -e_{n+1}^2 = 1\). This is an extension of a standard conformal algebra, with an extra \(e_{n+1}\) basis vector.

Note that we permuted the order in the source code below to make ConformalLayout happy.

from clifford import ConformalLayout, BasisVectorIds, MultiVector, transformations

class OurCustomLayout(ConformalLayout):
    def __init__(self, ndims):
        self.ndims = ndims

        euclidean_vectors = [str(i + 1) for i in range(ndims)]
        conformal_vectors = ['m2', 'm1']

        # Construct our custom algebra. Note that ConformalLayout requires the e- and e+ basis vectors to be last.
        ConformalLayout.__init__(
            self,
            [1]*ndims + [-1] + [1, -1],
            ids=BasisVectorIds(euclidean_vectors + ['np1'] + conformal_vectors)
        )
        self.enp1 = self.basis_vectors_lst[ndims]

        # Construct a base algebra without the extra `enp1`, which would not be understood by pyganja.
        self.conformal_base = ConformalLayout(
            [1]*ndims + [1, -1],
            ids=BasisVectorIds(euclidean_vectors + conformal_vectors)
        )

        # this lets us convert between the two layouts
        self.to_conformal = transformations.between_basis_vectors(self, self.conformal_base)

The code above also defines a stardard conformal \(\mathbb{R}^{N+1,1}\) layout without this new basis vector. This is primarily to support rendering with pyganja, which doesn’t support the presence of this extra vector. BasisVectorMap defaults to preserving vectors by name between one algebra and another, while throwing away blades containing vectors missing from the destination algebra.

We define an ups function which maps conformal dual-spheres into this algebra, as \(s^\prime = s + \left|s\right|e_{n+1}\), and a downs that applies the correct sign. The s suffix here is chosen to mean sphere.

def ups(self, s):
    return s + self.enp1*abs(s)

OurCustomLayout.ups = ups; del ups

def downs(self, mv):
    if (mv | self.enp1)[()] > 0:
        mv = -mv
    return mv

OurCustomLayout.downs = downs; del downs

Before we start looking at specified dimensions of euclidean space, we build a helper to construct conformal dual circles and spheres, with the word round being a general term intended to cover both circles and spheres.

def dual_round(at, r):
    l = at.layout
    return l.up(at) - 0.5*l.einf*r*r

In order to render with pyganja, we’ll need a helper to convert from our custom \(\mathbb{R}^{N+1,2}\) layout into a standard conformal \(\mathbb{R}^{N+1,1}\) layout. clifford maps indices in .value to basis blades via layout._basis_blade_order.index_to_bitmap, which we can use to convert the indices in one layout to the indices in another.

Visualization

Finally, we’ll define a plotting function, which plots the problem and solution circles in suitable colors via pyganja. Note that this all works because of our definition of the to_conformal BasisVectorMap.

import itertools
from pyganja import GanjaScene, draw

def plot_rounds(in_rounds, out_rounds, scale=1):
    colors = itertools.cycle([
        (255, 0, 0),
        (0, 255, 0),
        (0, 0, 255),
        (0, 255, 255),
    ])
    # note: .dual() neede here because we're passing in dual rounds, but ganja expects direct rounds
    s = GanjaScene()
    for r, color in zip(in_rounds, colors):
        s.add_object(r.layout.to_conformal(r).dual(), color=color)
    for r in out_rounds:
        s.add_object(r.layout.to_conformal(r).dual(), color=(64, 64, 64))
    draw(s, sig=r.layout.conformal_base.sig, scale=scale)

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

Apollonius’ problem in \(\mathbb{R}^2\) with circles

l2 = OurCustomLayout(ndims=2)
e1, e2 = l2.basis_vectors_lst[:2]

This gives us the Layout l2 with the desired metric,

import pandas as pd  # convenient but somewhat slow trick for showing tables
pd.DataFrame(l2.metric, index=l2.basis_names, columns=l2.basis_names)

	e1	e2	enp1	em2	em1
e1	1.0	0.0	0.0	0.0	0.0
e2	0.0	1.0	0.0	0.0	0.0
enp1	0.0	0.0	-1.0	0.0	0.0
em2	0.0	0.0	0.0	1.0	0.0
em1	0.0	0.0	0.0	0.0	-1.0

Now we can build some dual circles:

# add minus signs before `dual_round` to flip circle directions
c1 = dual_round(-e1-e2, 1)
c2 = dual_round(e1-e2, 0.75)
c3 = dual_round(e2, 0.5)

Compute the space orthogonal to all of them, which is an object of grade 2:

pp = (l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual()
pp.grades()

{2}

We hypothesize that this object is of the form l2.ups(c4) ^ l2.ups(c5). Taking a step not mentioned in the original paper, we decide to treat this as a regular conformal point pair, which allows us to project out the two factors with the approach taken in [LLW04]. Here, we normalize with \(e_{n+1}\) instead of the usual \(n_\infty\):

def pp_ends(pp):
    P = (1 + pp.normal()) / 2
    return P * (pp | pp.layout.enp1), ~P * (pp | pp.layout.enp1)

c4u, c5u = pp_ends(pp)

And finally, plot our circles:

plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)], scale=0.75)

This works for colinear circles too:

c1 = dual_round(-1.5*e1, 0.5)
c2 = dual_round(e1*0, 0.5)
c3 = dual_round(1.5*e1, 0.5)
c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())

plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])

c1 = dual_round(-3*e1, 1.5)
c2 = dual_round(-2*e1, 1)
c3 = -dual_round(2*e1, 1)
c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())

plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])

Apollonius’ problem in \(\mathbb{R}^3\) with spheres

l3 = OurCustomLayout(ndims=3)
e1, e2, e3 = l3.basis_vectors_lst[:3]

Again, we can check the metric:

pd.DataFrame(l3.metric, index=l3.basis_names, columns=l3.basis_names)

	e1	e2	e3	enp1	em2	em1
e1	1.0	0.0	0.0	0.0	0.0	0.0
e2	0.0	1.0	0.0	0.0	0.0	0.0
e3	0.0	0.0	1.0	0.0	0.0	0.0
enp1	0.0	0.0	0.0	-1.0	0.0	0.0
em2	0.0	0.0	0.0	0.0	1.0	0.0
em1	0.0	0.0	0.0	0.0	0.0	-1.0

And apply the solution to some spheres, noting that we now need 4 in order to constrain our solution

c1 = dual_round(e1+e2+e3, 1)
c2 = dual_round(-e1+e2-e3, 0.25)
c3 = dual_round(e1-e2-e3, 0.5)
c4 = dual_round(-e1-e2+e3, 1)
c5u, c6u = pp_ends((l3.ups(c1) ^ l3.ups(c2) ^ l3.ups(c3) ^ l3.ups(c4)).dual())

plot_rounds([c1, c2, c3, c4], [l3.downs(c6u), l3.downs(c5u)], scale=0.25)

Note that the figure above can be rotated!

Other resources for clifford

Slide Decks

• Installation

• Conformal Geometric Algebra with Python, Part 1

• Conformal Geometric Algebra with Python, Part 2

Videos

• Intro to Clifford

Other resources for Geometric Algebra

If you think Geometric Algebra looks interesting and want to learn more, check out these websites and textbooks

Links

• galgebra, a symbolic geometric algebra module for Python

• The Cambridge University Geometric Algebra Research Group home page

• David Hestenes’ Geometric Calculus R & D Home Page

Introductory textbooks

• Geometric Algebra for Physicists, by Doran and Lasenby

• Geometric Algebra for Computer Science, by Dorst, Fontijne and Mann

• New Foundations for Classical Mechanics, by David Hestenes

Bibliography

Cor17

Allan Cortzen. Geometric algebra algorithm representing an orthogonal isomorphism in versor form v2.03. Available online, 10 2017. URL: https://ctz.dk/geometric-algebra/frames-to-versor-algorithm/.

DL03

Chris Doran and Anthony Lasenby. Geometric Algebra for Physicists. Cambridge University Press, 2003. URL: http://www.mrao.cam.ac.uk/~clifford.

DFM09

Leo Dorst, Daniel Fontijne, and Stephen Mann. Geometric Algebra for Computer Science: An Object-Oriented Approach to Geometry. Morgan Kaufmann Publishers Inc., 2009. ISBN 9780080553108.

DV11

Leo Dorst and Robert Valkenburg. Square Root and Logarithm of Rotors in 3D Conformal Geometric Algebra Using Polar Decomposition, pages 81–104. Springer London, 01 2011. doi:10.1007/978-0-85729-811-9_5.

FD11

Daniel Fontijne and Leo Dorst. Reconstructing rotations and rigid body motions from exact point correspondences through reflections. Journal of Theoretical Biology - J THEOR BIOL, pages 63–78, 01 2011. doi:10.1007/978-0-85729-811-9_4.

HL19

Hugo Hadfield and Joan Lasenby. Direct linear interpolation of geometric objects in conformal geometric algebra. Advances in Applied Clifford Algebras, 29:, 09 2019. doi:10.1007/s00006-019-1003-y.

Hes15

David Hestenes. Space-time algebra. Springer Berlin Heidelberg, New York, second edition edition, 2015. ISBN 978-3-319-18412-8 978-3-319-18413-5.

HS17

Eckhard Hitzer and Stephen Sangwine. Multivector and multivector matrix inverses in real clifford algebras. Applied Mathematics and Computation, 311:375–389, Oct 2017. doi:10.1016/j.amc.2017.05.027.

LLW04

Anthony Lasenby, Joan Lasenby, and Rich Wareham. A Covariant Approach to Geometry using Geometric Algebra. Technical Report F-INFENG/TR-483, Department of Engineering, University of Cambridge, 2004. URL: https://pdfs.semanticscholar.org/baba/976fd7f6577eeaa1d3ef488c1db13ec24652.pdf.

LHL19

Joan Lasenby, Hugo Hadfield, and Anthony Lasenby. Calculating the rotor between conformal objects. Advances in Applied Clifford Algebras, 29:102, 10 2019. doi:10.1007/s00006-019-1014-8.

Shi20

D. S. Shirokov. On determinant, other characteristic polynomial coefficients, and inverses in clifford algebras of arbitrary dimension. 2020. arXiv:2005.04015.

WCL05

Rich Wareham, Jonathan Cameron, and Joan Lasenby. Applications of conformal geometric algebra in computer vision and graphics. In Hongbo Li, Peter J. Olver, and Gerald Sommer, editors, Computer Algebra and Geometric Algebra with Applications, 329–349. Berlin, Heidelberg, 2005. Springer Berlin Heidelberg.

WL08

Rich Wareham and Joan Lasenby. Mesh vertex pose and position interpolation using geometric algebra. In Articulated Motion and Deformable Objects, 122–131. Berlin, Heidelberg, 2008. Springer Berlin Heidelberg.