# clifford: Geometric Algebra for Python¶

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

In : import math

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

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

In : R*a*~R    # rotate the vector
Out: (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. ## 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(*args, **kw) 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.

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(*args, **kwds) 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)[source]

Get/Set the epsilon for float comparisons.

clifford.pretty(precision=None)[source]

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()[source]

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

clifford.print_precision(newVal)[source]

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
>>> 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 psuedo-code used the create the algorithms below can be found at Allan Cortzen’s website.

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

### 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 ^ Ms ^ Ms

clifford.operator.gp(M, N)[source]

Geometric product function $$MN$$, equivalent to M * N.

M and N must be from the same layout

clifford.operator.op(M, N)[source]

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

M and N must be from the same layout

clifford.operator.ip(M, N)[source]

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 = typing.Callable[[clifford._multivector.MultiVector], clifford._multivector.MultiVector]

A callable mapping one MultiVector to another.

class clifford.transformations.Linear(*args, **kwds)[source]

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(*args, **kwds)[source]

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(*args, **kwds)[source]

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.

clifford.BladeMap

A faster but less general approach that works on basis blades

property adjoint

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

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)


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[source]

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(*args, **kwds)[source]

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.

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: Dict[Any, Any] = None)clifford.transformations.OutermorphismMatrix[source]

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

Parameters

Example

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

## 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 : from clifford import g2

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


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

clifford.<predefined>.layout

The associated clifford.Layout

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

clifford.<predefined>.blades

A shorthand for Layout.blades()

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


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

In : from clifford.g2 import *

In : e1, e2, e12
Out: ((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¶

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

### Changes in 1.3.x¶

#### Patch releases¶

• 1.3.1: Added compatibility with numba version 0.50.0.

### 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() ==  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¶

• 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

• 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),

:

from numpy import e,pi
import clifford as cf

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


:

blades

:

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


:

e1 = blades['e1']


### 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 numpy  import 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']
# 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

:

from scipy import cos, sin

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-14-7dddb582d1e8> in <module>
1 from scipy import cos, sin
2
----> 3 theta = pi/4
4 R = cos(theta) - sin(theta)*e23
5 R

NameError: name 'pi' is not defined


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)


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.0^e12) + (2.0^e13) - (1.0^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),
names=['', 'e1', 'e2', 'e3', 'e12', 'e13', 'e23', 'e123']),
[0.13414634146341464, 0.12195121951219513, -0.0, -0.14634146341463417, 0.18292682926829273, -2.0816681711721682e-17, 0.0975609756097561, -0.29268292682926833])


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

:

import math

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: e**(-theta/2*e12)
R12(pi/2)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-31-517400d2f6b9> in <module>
1 R12 = lambda theta: e**(-theta/2*e12)
----> 2 R12(pi/2)

NameError: name 'pi' is not defined


And use it like this

:

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


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-32-62fda315a2cc> in <module>
1 a = e1+e2+e3
----> 2 R = R12(math.pi/2)
3 R*a*~R

<ipython-input-31-517400d2f6b9> in <lambda>(theta)
----> 1 R12 = lambda theta: e**(-theta/2*e12)
2 R12(pi/2)

NameError: name 'e' is not defined


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(pi/6*(e23+e12))  # this returns a function
R(a)  # which acts on a vector

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-36-402a96f6039a> in <module>
5     return apply_rotation
6
----> 7 R = R_factory(pi/6*(e23+e12))  # this returns a function
8 R(a)  # which acts on a vector

NameError: name 'pi' is not defined


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(pi/2*e23, pi/2*e12, e1)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-38-21fc9c28e3a1> in <module>
9
10 # rotation sequence by  pi/2-in-e12 THEN pi/2-in-e23
---> 11 R_seq(pi/2*e23, pi/2*e12, e1)

NameError: name 'pi' is not defined


### 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'])


:

{'': 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 “Geometric Algebra for Physicists”, 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, so we just enter the frame components into a matrix.

:

from numpy import array

M = [float(b|a) for b in B for a in A] # you need float() due to bug in clifford
M = array(M).reshape(3,3)

M

:

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


Thats a rotation matrix.

### Convert a Rotation Matrix to a Rotor¶

In 3 Dimenions, 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 “Geometric Algebra for Physicists”. 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[1,0]*e2 + M[2,0]*e3,
M[0,1]*e1 + M[1,1]*e2 + M[2,1]*e3,
M[0,2]*e1 + M[1,2]*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)


blam.

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

# put elements of each in namespace


### 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

(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

:

-(13.2^d0) + (7.2^d1) - (8.1^d2) + (0.8^d3)


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

:

X*d0

:

-13.2 - (7.2^d01) + (8.1^d02) - (0.8^d03)


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

:

bm(X*d0)

:

-13.2 - (7.2^p1) + (8.1^p2) - (0.8^p3)


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

:

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

:

-13.2

:

x(1) # the space component

:

-(7.2^p1) + (8.1^p2) - (0.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)

:

-13.2 - (7.2^p1) + (8.1^p2) - (0.8^p3)


The split can be inverted by applying the BladeMap again, and multiplying by $$d_0$$

:

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

:

-(13.2^d0) + (7.2^d1) - (8.1^d2) + (0.8^d3)


### Splitting a Bivector¶

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

:

F = D.randomMV()(2)
F

:

(0.5^d01) + (1.5^d02) + (1.0^d03) - (0.9^d12) - (2.0^d13) - (1.5^d23)


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

:

split(F)

:

(0.5^p1) + (1.5^p2) + (1.0^p3) - (0.9^p12) - (2.0^p13) - (1.5^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.5^p1) + (1.5^p2) + (1.0^p3)

:

iB

:

-(0.9^p12) - (2.0^p13) - (1.5^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

:

-0.1 + (0.2^d01) - (0.3^d02) - (1.3^d03) - (0.1^d12) - (0.0^d13) - (0.9^d23) - (0.3^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_

:

(0.8^d01) - (0.7^d02) - (7.3^d03) - (1.0^d12) - (0.1^d13) - (7.5^d23)


Then splitting into $$E$$ and $$B$$ fields

:

E_ = split(F_)(1)
E_

:

(0.8^p1) - (0.7^p2) - (7.3^p3)

:

iB_ = split(F_)(2)
iB_

:

-(1.0^p12) - (0.1^p13) - (7.5^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

:

-3.4 + (2.6^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

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

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']



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 + q*i + q*j + q*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.0^z) - (3.0^y) + (4.0^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


:

{'': 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 + q*e13 +q*e23 + q*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

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):

:

#%%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

def apply_rotor_val(R_val,mv_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):

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)
<ipython-input-16-059d06919085> 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())

<ipython-input-16-059d06919085> in meet_wrapped(mv_a, mv_b)
17
18 def meet_wrapped(mv_a,mv_b):
---> 19     return cf.layout.MultiVector(meet_val(mv_a.value, mv_b.value))
20
21 sphere = (up(0)^up(e1)^up(e2)^up(e3)).normal()

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
l_list = non_zero_indices
m_list = non_zero_indices
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')

'''
Returns a function that implements the mult_table on two input multivectors
'''
non_zero_indices = mult_table.nonzero()
k_list = non_zero_indices
l_list = non_zero_indices
m_list = non_zero_indices
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

@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])

@numba.njit
def sparse_apply_rotor_val(R_val,mv_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,,[0,1,2,3,4,5])

@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)
<ipython-input-24-b25f7cce90d8> in <module>
14
15 print(sphere.meet(line_one).normal().normal())
---> 16 print(meet_sparse_3_4(line_one,sphere).normal())

<ipython-input-24-b25f7cce90d8> in meet_sparse_3_4(mv_a, mv_b)
11
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))
14
15 print(sphere.meet(line_one).normal().normal())

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


:

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


Now, conformalize it

:

G2c, blades_g2c, stuff = conformalize(G2)


:

{'': 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))


<ipython-input-10-0e2e8e068b72>:4: DeprecationWarning: scipy.rand is deprecated and will be removed in SciPy 2.0.0, use numpy.random.rand instead
alpha = rand()
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 )

<ipython-input-10-0e2e8e068b72>: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 $$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.12'


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

:

sc = GanjaScene()

:

sc_refl = GanjaScene()


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__

:

'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_refl = GanjaScene()


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

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

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(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:264: 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:264: 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.70621^e1234) + (1.15174^e1235) + (0.0035^e1245) - (0.18868^e1345) - (0.36962^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

:

-(1.41653^e1) + (0.4833^e2) - (0.15024^e3)

:

cga.up(a)

:

-(1.41653^e1) + (0.4833^e2) - (0.15024^e3) + (0.63135^e4) + (1.63135^e5)

:

cga.null_vector()  # create null vector directly

:

-(1.3863^e1) + (0.22667^e2) + (0.76895^e3) + (0.78224^e4) + (1.78224^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.07715^e123) - (0.50708^e124) - (1.0196^e125) + (0.19953^e134) + (0.36228^e135) - (0.25575^e145) - (0.13719^e234) - (0.25387^e235) + (0.14448^e245) + (0.01234^e345)

:

C = cga.round(e1, e2, -e1, e3)

:

(-(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)

:

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

###### 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 *

:

{'': 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

euc_vector_in_plane_m = e1
euc_vector_in_plane_n = e3

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

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()
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 << 16) | (intermediate_color << 8) | intermediate_color
)
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)