ANN: Clifford Algebras in Python
kern at caltech.edu
kern at caltech.edu
Fri Aug 25 20:43:28 EDT 2000
Hello all three of you who would care about this module!
At long last, I felt confident enough in my understanding of
Geometric (Clifford) Algebras to code up an implementation in
Python (gracefully assisted by NumPy) which I call pyCA.
(Note: this is numerical, not symbolic.)
The documentation strings are fairly complete (read: I was too
lazy to write more documentation without a good reason). pyCA
is a fairly general framework for doing work in geometric algebras.
It's not especially fast for large numbers of computations or
efficient for many multivectors. It also concedes to Python's
operator precedence. I plan to add modules over the framework
to implement a convenient shell and implement specific algebras
and subalgebras (3D Euclidean, quaternions, STA, PGA, etc.).
Well, have fun with the code below. I'll place it on the Starship
once I get around to getting my account re-activated.
Apologies for the line-breaks.
----------%<-------------------------------------------------------
"""\
pyCA 0.5 -- Geometric Algebra for Python
Robert Kern
kern at caltech.edu
This module implements geometric algebras (a.k.a. Clifford algebras).
For the uninitiated, a geometric algebra is an algebra of vectors of
given dimensions and signature. The familiar inner (dot) product and the
outer product, a generalized relative of the three-dimensional cross
product are unified in an invertible geometric product. Scalars, vectors,
and higher-grade entities can be mixed freely and consistently in the
form of mixed-grade multivectors.
For more information, see the Cambridge Clifford Algebras website:
http://www.mrao.cam.ac.uk/~clifford
and David Hestenes' website:
http://modelingnts.la.asu.edu/GC_R&D.html
This implementation is inspired by Christian Perwass' XGAlib, a C++
library for GA. Although there is practically no code shared with
that LGPLed library, I have still obtained permission from Mr. Perwass
to release this code into the public domain. In addition, I have
obtained permission from James Lehmann to use a snippet of his code
posted on comp.lang.python . Thus:
I hereby release this code into the PUBLIC DOMAIN AS IS. There is no
support, warranty, or guarantee. I will gladly accept comments, bug
reports, and patches, however.
Two classes, Layout and MultiVector, and several helper functions are
provided to implement the algebras.
Helper Functions
================
Cl(p, q=0, names=None, firstIdx=0) -- generates the geometric algebra
Cl_p,q and returns the appropriate Layout and dictionary of basis
element names to their MultiVector instances.
bases(layout) -- returns the dictionary of basis element names: instances
as above
randomMV(layout, grades=None) -- random MultiVector using layout. If
grades is a sequence of integrers, only those grades will be non-zero
pretty() -- set repr() to pretty-print MultiVectors
ugly() -- set repr() to return eval-able strings
eps(newEps) -- set _eps, the epsilon for comparisons and zero-testing
Issues
======
* Due to Python's order of operations, the bit operators & ^ << follow
the normal arithmetic operators + - * /, so
1^e0 + 2^e1 != (1^e0) + (2^e1)
as is probably intended. Additionally (comments on this pun will be
summarily ignored),
M = MultiVector(layout2D) # null multivector
M << 1^e0 << 2^e1 == 10.0^e1 + 1.0^e01
M == 1.0
e0 == 2 + 1^e0
as is definitely not intended. However,
M = MultiVector(layout2D)
M << (2^e0) << e1 << (3^e01) == M == 2^e0 + 1^e1 + 3^e01
e0 == 1^e0
e1 == 1^e1
e01 == 1^e01
A shell that reorders the operators and adds a few conveniences is in
the works.
* 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.
* If the LinearAlgebra module from the NumPy distribution is available,
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**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 = ~M/(M&~M) iff
M&~M == |M|**2] fails for those multivectors where M&~M is not a scalar.
It is only used if the inv method is manually set to point to normalInv
or the module LinearAlgebra from NumPy is not available.
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.
* 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 is in preparation
that will calculate these tables for you (and include methods for
generating reciprocal bases and the like). I would greatly appreciate
help on the algorithm to generate these tables from an arbitrary
bilinear form.
* No care is taken to preserve the typecode 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.
* The way I make variables of the bases is by .update()'ing the global
(or local for functions) namespace dictionary. I use a convenient
function for that in the global case. I have not included it here
as I don't want to canonize something that is, well, let's face it,
evil. But it's a venial kind of evil, so do it anyways.
* Integration with PyOpenGL for some kick-ass demonstrations are left
as an exercise for the reader.
Happy hacking!
Robert Kern
kern at caltech.edu
"""
import Numeric
import types
import operator
import math
from Numeric import pi, e
try:
import LinearAlgebra
LA = LinearAlgebra
except ImportError:
LA = None
_eps = 1e-15 # float epsilon for float comparisons
_pretty = 0 # pretty-print global
class Layout:
"""\
Layout stores information regarding the geometric algebra itself and
the internal representation of multivectors. It is constructed like this:
Layout(signature, bladeList, firstIdx=0, names=None)
The arguments to be passed are interpreted as follows:
signature -- the signature of the vector space. This should be
a list of positive and negative numbers where the sign determines
the sign of the inner product of the corresponding vector with itself.
The values are irrelevant except for sign. This list also determines
the dimensionality of the vectors. Signatures with zeroes are not
permitted at this time.
Examples:
signature = [+1, -1, -1, -1] # Hestenes', et al. Space-Time Algebra
signature = [+1, +1, +1] # 3-D Euclidean signature
bladeList -- list of tuples corresponding to the blades in the whole
algebra. This list determines the order of coefficients in the
internal representation of multivectors. The entry for the scalar
must be an empty tuple, and the entries for grade-1 vectors must be
singleton tuples. Remember, the length of the list will be 2**dims.
Example:
bladeList = [(), (0,), (1,), (0,1)] # 2-D
firstIdx -- the index of the first vector. That is, some systems number
the base vectors starting with 0, some with 1. Choose by passing
the correct number as firstIdx. 0 is the default.
names -- list of names of each blade. When pretty-printing multivectors,
use these symbols for the blades. names should be in the same order
as bladeList. You may use an empty string for scalars. By default,
the name for each non-scalar blade is 'e' plus the indices of the blade
as given in bladeList.
Example:
names = ['', 's0', 's1', 'i'] # 2-D
Layout's Members:
dims -- dimensionality of vectors (== len(signature))
sig -- normalized signature (i.e. all values are +1 or -1)
firstIdx -- starting point for vector indices
bladeList -- list of blades
gradeList -- corresponding list of the grades of each blade
gaDims -- 2**dims
names -- pretty-printing symbols for the blades
even -- dictionary of even permutations of blades to the canonical blades
odd -- dictionary of odd permutations of blades to the canonical blades
gmt -- multiplication table for geometric product [1]
imt -- multiplication table for inner product [1]
omt -- multiplication table for outer product [1]
[1] The multiplication tables are NumPy arrays of rank 3 with indices like
the tensor g_ijk discussed above.
"""
def __init__(self, sig, bladeList, firstIdx=0, names=None):
self.dims = len(sig)
self.sig = Numeric.divide(sig,
Numeric.absolute(sig)).astype(Numeric.Int)
self.firstIdx = firstIdx
self.bladeList = map(tuple, bladeList)
self._checkList()
self.gaDims = len(self.bladeList)
self.gradeList = map(len, self.bladeList)
if names is None:
e = 'e'
self.names = []
for i in range(self.gaDims):
if self.gradeList[i] > 1:
self.names.append(e + str(Numeric.add.reduce(map(str, self.bladeList[i]))))
elif self.gradeList[i] == 1:
self.names.append(e + str(self.bladeList[i][0]))
else:
self.names.append('')
elif len(names) == self.gaDims:
self.names = names
else:
raise ValueError, "names list of length %i needs to be of length %i" % \
(len(names), self.gaDims)
self._genEvenOdd()
self._genTables()
def __repr__(self):
s = "Layout(%s, %s, firstIdx=%s, names=%s)" % (list(self.sig),
self.bladeList,
self.firstIdx,
self.names)
return s
def _sign(self, seq, orig):
"""Determine {even,odd}-ness of permutation seq or orig.
Returns 1 if even; -1 if odd.
"""
sign = 1
seq = list(seq)
for i in range(len(seq)):
if seq[i] != orig[i]:
j = seq.index(orig[i])
sign = -sign
seq[i], seq[j] = seq[j], seq[i]
return sign
def _containsDups(self, list):
"Checks if list contains duplicates."
for k in list:
if list.count(k) != 1:
return 1
return 0
def _genEvenOdd(self):
"Create mappings of even and odd permutations to their canonical blades."
self.even = {}
self.odd = {}
class NoMorePermutations(StandardError):
pass
for i in range(self.gaDims):
blade = self.bladeList[i]
grade = self.gradeList[i]
if grade in (0, 1):
# handled easy cases
self.even[blade] = blade
continue
elif grade == 2:
# another easy case
self.even[blade] = blade
self.odd[(blade[1], blade[0])] = blade
continue
else:
# general case, lifted from Chooser.py released on
# comp.lang.python by James Lehmann.
idx = range(grade)
try:
for i in range(Numeric.multiply.reduce(range(1, grade+1))):
# grade! permutations
done = 0
j = grade - 1
while not done:
idx[j] = idx[j] + 1
while idx[j] == grade:
idx[j] = 0
j = j - 1
idx[j] = idx[j] + 1
if j == -1:
raise NoMorePermutations()
j = grade - 1
if not self._containsDups(idx):
done = 1
perm = []
for k in idx:
perm.append(blade[k])
perm = tuple(perm)
if self._sign(perm, blade) == 1:
self.even[perm] = blade
else:
self.odd[perm] = blade
except NoMorePermutations:
pass
self.even[blade] = blade
def _checkList(self):
"Ensure validity of arguments."
# check for uniqueness
for blade in self.bladeList:
if self.bladeList.count(blade) != 1:
raise ValueError, "blades not unique"
# check for right dimensionality
if len(self.bladeList) != 2**self.dims:
raise ValueError, "incorrect number of blades"
# check for valid ranges of indices
valid = range(self.firstIdx, self.firstIdx+self.dims)
try:
for blade in self.bladeList:
for idx in blade:
if (idx not in valid) or (list(blade).count(idx) != 1):
raise ValueError
except (ValueError, TypeError):
raise ValueError, "invalid bladeList; must be a list of tuples"
def _gmtElement(self, a, b):
"Returns the element of the geometric multiplication table given blades a, b."
mul = 1 # multiplier
newBlade = list(a) + list(b)
unique = 0
while unique == 0:
for i in range(len(newBlade)):
if newBlade.count(newBlade[i]) != 1:
delta = newBlade[i+1:].index(newBlade[i])
eo = 1 - 2*(delta % 2)
# eo == 1 if the elements are an even number of flips away
# eo == -1 otherwise
del newBlade[i+delta+1]
del newBlade[i]
mul = eo * mul * self.sig[a[i] - self.firstIdx]
break
else:
unique = 1
newBlade = tuple(newBlade)
if newBlade in self.bladeList:
idx = self.bladeList.index(newBlade)
# index of the product blade in the bladeList
elif newBlade in self.even.keys():
# even permutation
idx = self.bladeList.index(self.even[newBlade])
else:
# odd permutation
idx = self.bladeList.index(self.odd[newBlade])
mul = -mul
element = Numeric.zeros((self.gaDims,), Numeric.Int)
element[idx] = mul
return element, idx
def _genTables(self):
"Generate the multiplication tables."
# geometric multiplication table
gmt = Numeric.zeros((self.gaDims, self.gaDims, self.gaDims),
Numeric.Int)
# inner product table
imt = Numeric.zeros((self.gaDims, self.gaDims, self.gaDims),
Numeric.Int)
# outer product table
omt = Numeric.zeros((self.gaDims, self.gaDims, self.gaDims),
Numeric.Int)
for i in range(self.gaDims):
for j in range(self.gaDims):
gmt[i,:,j], idx = self._gmtElement(list(self.bladeList[i]),
list(self.bladeList[j]))
if self.gradeList[idx] == abs(self.gradeList[i] - self.gradeList[j]) and \
self.gradeList[i] != 0 and self.gradeList[j] != 0:
# A_r . B_s = <A_r B_s>_|r-s|
# if r,s != 0
imt[i,:,j] = gmt[i,:,j]
if self.gradeList[idx] == (self.gradeList[i] + self.gradeList[j]):
# A_r ^ B_s = <A_r B_s>_|r+s|
omt[i,:,j] = gmt[i,:,j]
self.gmt = gmt
self.imt = imt
self.omt = omt
def pseudoScalar(self):
"Returns a MultiVector that is the pseudoscalar of this space."
psIdx = self.gradeList.index(self.dims)
# pick out out blade with grade equal to dims
pScalar = MultiVector(self)
pScalar.value[psIdx] = 1
return pScalar
def invPS(self):
"Returns the inverse of the pseudoscalar of the algebra."
ps = self.pseudoScalar()
return ps.inv()
class MultiVector:
"""\
The elements of the algebras, the multivectors, are implemented in the
MultiVector class. It has the following constructor:
MultiVector(layout, value=None)
The meaning of the arguments is as follows:
layout -- instance of Layout
value -- a sequence, of length layout.gaDims, of coefficients of the base
blades
MultiVector's Members
=====================
layout -- instance of Layout
value -- a NumPy array, of length layout.gaDims, of coefficients of the base
blades
MultiVector's Public Methods
============================
__and__, __rand__ -- geometric product. M & N
__xor__, __rxor__ -- outer product. M ^ N
__mul__, __rmul__ -- inner product. M * N
__add__, __radd__ -- addition. M + N
__sub__, __rsub__ -- subtraction. M - N
__div__, __rdiv__ -- division. M / N <==> M & (N.inv())
__pow__ -- exponentiation of a multivector by an integer. M ** n
__rpow__ -- exponentiation of a scalar by a multivector. x ** M
__lshift__ -- in-place addition. M << N
Adds N to M, modifying M and returning M.
mag2() -- magnitude squared. (M & ~M)[()] == |M|**2
__abs__ -- magnitude. abs(M) == |M|
adjoint(), __invert__ -- adjoint/reversion. ~M
__int__, __long__, __float__ -- cast to scalar Python number iff only the
multivector's scalar coefficient is non-zero
__getitem__, __setitem__ -- get/set coefficient to/from a Python number.
Item can be an index into self.value or a
blade (or a permutation of a blade).
__delitem__ -- sets coefficient of item (same rules as above) to 0
__getslice__ -- returns a new MultiVector with only the selected slice
non-zero
__setslice__ -- sets the slice of self.value to a value or values
__delslice__ -- sets the slice of self.value to 0
__call__ -- projects multivector onto the specified grade. M(grade)
__cmp__ -- returns 0 if the multivectors' coefficients are each equal to
within epsilon. Otherwise, returns a well-defined, but
meaningless, ordering.
__str__ -- pretty-printed representation
__repr__ -- ugly, but eval-able by default, but pretty if global variable
_pretty is true
isScalar() -- returns 1 if multivector is a scalar
isBlade() -- returns 1 if multivector is a blade
grades() -- returns the grades contained in the multivector
normalInv() -- multivector inverse iff (M & ~M) == |M|**2. ~M/|M|**2
laInv() -- experimental, more general inverse using computational
linear algebra to solve for the inverse.
inv() -- one of normalInv or laInv. By default, inv is laInv if the
LinearAlgebra module is available, normalInv otherwise.
dual() -- returns the dual of the multivector
commutator(other) -- returns the commutator of two multivectors
join(other) -- returns the join of two blades if it exists
meet(other, subspace=None) -- returns the meet of two blades wrt some
subspace which defaults to the pseudoscalar
gradeInvol() -- returns the grade involution of the multivector
conjugate() -- returns the Clifford conjugate of the multivector
"""
def __init__(self, layout, value=None):
"""Constructor.
MultiVector(layout, value=None) --> MultiVector
"""
self.layout = layout
if value is None:
self.value = Numeric.zeros((self.layout.gaDims,), Numeric.Float)
else:
self.value = Numeric.array(value)
if self.value.shape != (self.layout.gaDims,):
raise ValueError, "value must be a sequence of length %s" % self.layout.gaDims
def _checkOther(self, other, coerce=1):
"""Ensure that the other argument has the same Layout or coerce value if
necessary/requested.
_checkOther(other, coerce=1) --> newOther, isMultiVector
"""
if type(other) in (types.IntType, types.FloatType, types.LongType):
if coerce:
# numeric scalar
newOther = MultiVector(self.layout)
newOther[()] = other
return newOther, 1
else:
return other, 0
elif isinstance(other, self.__class__) and \
other.layout is not self.layout:
raise ValueError, "cannot operate on MultiVectors with different Layouts"
return other, 1
## numeric special methods
# binary
def __and__(self, other):
"""Geometric product
M & N --> MN
__and__(other) --> MultiVector
"""
other, mv = self._checkOther(other, coerce=0)
if mv:
newValue = Numeric.dot(self.value, Numeric.dot(self.layout.gmt,
other.value))
else:
newValue = other * self.value
return MultiVector(self.layout, newValue)
def __rand__(self, other):
"""Right-hand geometric product
N & M --> NM
__rand__(other) --> MultiVector
"""
other, mv = self._checkOther(other, coerce=0)
if mv:
newValue = Numeric.dot(other.value, Numeric.dot(self.layout.gmt,
self.value))
else:
newValue = other * self.value
return MultiVector(self.layout, newValue)
def __xor__(self, other):
"""Outer product
M ^ N
__xor__(other) --> MultiVector
"""
other, mv = self._checkOther(other, coerce=0)
if mv:
newValue = Numeric.dot(self.value, Numeric.dot(self.layout.omt,
other.value))
else:
newValue = other * self.value
return MultiVector(self.layout, newValue)
def __rxor__(self, other):
"""Right-hand outer product
N ^ M
__rxor__(other) --> MultiVector
"""
other, mv = self._checkOther(other, coerce=0)
if mv:
newValue = Numeric.dot(other.value, Numeric.dot(self.layout.omt,
self.value))
else:
newValue = other * self.value
return MultiVector(self.layout, newValue)
def __mul__(self, other):
"""Inner product
M * N
__mul__(other) --> MultiVector
"""
other, mv = self._checkOther(other)
if mv:
newValue = Numeric.dot(self.value, Numeric.dot(self.layout.imt,
other.value))
else:
return MultiVector(self.layout) # l * M = M * l = 0 for scalar l
return MultiVector(self.layout, newValue)
__rmul__ = __mul__
def __add__(self, other):
"""Addition
M + N
__add__(other) --> MultiVector
"""
other, mv = self._checkOther(other)
newValue = self.value + other.value
return MultiVector(self.layout, newValue)
__radd__ = __add__
def __sub__(self, other):
"""Subtraction
M - N
__sub__(other) --> MultiVector
"""
other, mv = self._checkOther(other)
newValue = self.value - other.value
return MultiVector(self.layout, newValue)
def __rsub__(self, other):
"""Right-hand subtraction
N - M
__rsub__(other) --> MultiVector
"""
other, mv = self._checkOther(other)
newValue = other.value - self.value
return MultiVector(self.layout, newValue)
def __div__(self, other):
"""Division
-1
M / N --> M & N
__div__(other) --> MultiVector
"""
other, mv = self._checkOther(other, coerce=0)
if mv:
return self & other.inv()
else:
newValue = self.value / other
return MultiVector(self.layout, newValue)
def __rdiv__(self, other):
"""Right-hand division
-1
N / M --> N & M
__rdiv__(other) --> MultiVector
"""
other, mv = self._checkOther(other)
return other & self.inv()
def __pow__(self, other):
"""Exponentiation of a multivector by an integer
n
M ** n --> M
__pow__(other) --> MultiVector
"""
if type(other) not in (types.IntType, types.FloatType):
raise ValueError, "exponent must be a Python Int or Float"
if abs(round(other) - other) > _eps:
raise ValueError, "exponent must have no fractional part"
other = int(round(other))
newMV = MultiVector(self.layout, list(self.value))
for i in range(1, other):
newMV = newMV & self
return newMV
def __rpow__(self, other):
"""Exponentiation of a real by a multivector
M
r**M --> r
__rpow__(other) --> MultiVector
"""
# check that other is a Python number, not a MultiVector
if type(other) not in (types.IntType, types.FloatType, types.LongType):
raise ValueError, "base must be a Python number"
intMV = math.log(other) & self
# pow(x, y) == exp(y & log(x))
newMV = MultiVector(self.layout) # null
nextTerm = MultiVector(self.layout)
nextTerm[()] = 1 # order 0 term of exp(x) Taylor expansion
n = 1
while nextTerm != 0:
# iterate until the added term is within _eps of 0
newMV << nextTerm
nextTerm = nextTerm & intMV / n
n = n + 1
else:
# squeeze out that extra little bit of accuracy
newMV << nextTerm
return newMV
def __lshift__(self, other):
"""In-place addition
M << N --> M + N
__lshift__(other) --> MultiVector
"""
other, mv = self._checkOther(other)
self.value = self.value + other.value
return self
# unary
def __neg__(self):
"""Negation
-M
__neg__() --> MultiVector
"""
newValue = -self.value
return MultiVector(self.layout, newValue)
def __pos__(self):
"""Positive (just a copy)
+M
__pos__(self) --> MultiVector
"""
newValue = self.value + 0 # copy
return MultiVector(self.layout, newValue)
def mag2(self):
"""Magnitude (modulus) squared
2
|M|
mag2() --> PyFloat | PyInt
"""
return (~self & self)[()]
def __abs__(self):
"""Magnitude (modulus)
abs(M) --> |M|
__abs__() --> PyFloat
"""
return Numeric.sqrt(self.mag2())
def adjoint(self):
"""Adjoint / reversion
_
~M --> M (any one of several conflicting notations)
~(N & M) --> ~M & ~N
adjoint() --> MultiVector
"""
# The multivector created by reversing all multiplications
newMV = MultiVector(self.layout)
grades = Numeric.array(self.layout.gradeList)
signs = Numeric.power(-1, grades*(grades-1)/2)
newMV.value = signs * self.value
return newMV
__invert__ = adjoint
# builtin
def __int__(self):
"""Coerce to an integer iff scalar.
int(M)
__int__() --> PyInt
"""
return int(self.__float__())
def __long__(self):
"""Coerce to a long iff scalar.
long(M)
__long__() --> PyLong
"""
return long(self.__float__())
def __float__(self):
""""Coerce to a float iff scalar.
float(M)
__float__() --> PyFloat
"""
if self.isScalar():
return float(self[()])
else:
raise ValueError, "non-scalar coefficients are non-zero"
## sequence special methods
def __len__(self):
"""Returns length of value array.
len(M)
__len__() --> PyInt
"""
return self.layout.gaDims
def __getitem__(self, key):
"""If key is a blade tuple (e.g. (0,1) or (1,3)), then return
the (real) value of that blade's coefficient.
Otherwise, treat key as an index into the list of coefficients.
M[blade] --> PyFloat | PyInt
M[index] --> PyFloat | PyInt
__getitem__(key) --> PyFloat | PyInt
"""
if key in self.layout.bladeList:
return self.value[self.layout.bladeList.index(key)]
elif key in self.layout.even.keys():
return self.value[self.layout.bladeList.index(self.layout.even[key])]
elif key in self.layout.odd.keys():
return -self.value[self.layout.bladeList.index(self.layout.odd[key])]
else:
return self.value[key]
def __setitem__(self, key, value):
"""If key is a blade tuple (e.g. (0,1) or (1,3)), then set
the (real) value of that blade's coeficient.
Otherwise treat key as an index into the list of coefficients.
M[blade] = PyFloat | PyInt
M[index] = PyFloat | PyInt
__setitem__(key, value)
"""
if key in self.layout.bladeList:
self.value[self.layout.bladeList.index(key)] = value
elif key in self.layout.even.keys():
self.value[self.layout.bladeList.index(self.layout.even[key])] = value
elif key in self.layout.odd.keys():
self.value[self.layout.bladeList.index(self.layout.odd[key])] = -value
else:
self.value[key] = value
def __delitem__(self, key):
"""Set the selected coefficient to 0.
del M[blade]
del M[index]
__delitem__(key)
"""
if key in self.layout.bladeList:
self.value[self.layout.bladeList.index(key)] = 0
elif key in self.layout.even.keys():
self.value[self.layout.bladeList.index(self.layout.even[key])] = 0
elif key in self.layout.odd.keys():
self.value[self.layout.bladeList.index(self.layout.odd[key])] = 0
else:
self.value[key] = 0
def __getslice__(self, i, j):
"""Return a copy with only the slice non-zero.
M[i:j]
__getslice__(i, j) --> MultiVector
"""
newMV = MultiVector(self.layout)
newMV.value[i:j] = self.value[i:j]
return newMV
def __setslice__(self, i, j, sequence):
"""Paste a sequence into coefficients array.
M[i:j] = sequence
__setslice__(i, j, sequence)
"""
self.value[i:j] = sequence
def __delslice__(self, i, j):
"""Set slice to zeros.
del M[i:j]
__delslice__(i, j)
"""
self.value[i:j] = 0
## grade projection
def __call__(self, grade):
"""Return a new multi-vector projected onto the specified grade.
M(grade) --> <M>
grade
__call__(grade) --> MultiVector
"""
if grade not in self.layout.gradeList:
raise ValueError, "algebra does not have grade %s" % grade
if type(grade) is not types.IntType:
raise ValueError, "grade must be an integer"
mask = Numeric.equal(grade, self.layout.gradeList)
newValue = Numeric.multiply(mask, self.value)
return MultiVector(self.layout, newValue)
## fundamental special methods
def __str__(self):
"""Return pretty-printed representation.
str(M)
__str__() --> PyString
"""
s = ''
for i in range(self.layout.gaDims):
# if we have nothing yet, don't use + and - as operators but
# use - as an unary prefix if necessary
if s:
seps = (' + ', ' - ')
else:
seps = ('', '-')
if self.layout.gradeList[i] == 0:
# scalar
if abs(self.value[i]) >= _eps:
if self.value[i] > 0:
s = '%s%s%s' % (s, seps[0], self.value[i])
else:
s = '%s%s%s' % (s, seps[1], -self.value[i])
else:
if abs(self.value[i]) >= _eps:
# not a scalar
if self.value[i] > 0:
s = '%s%s%s^%s' % (s, seps[0], self.value[i],
self.layout.names[i])
else:
s = '%s%s%s^%s' % (s, seps[1], -self.value[i],
self.layout.names[i])
if s:
# non-zero
return s
else:
# return scalar 0
return '0'
def __repr__(self):
"""Return eval-able representation if global _pretty is false.
Otherwise, return str(self).
repr(M)
__repr__() --> PyString
"""
if _pretty:
return self.__str__()
s = "MultiVector(%s, value=%s)" % \
(repr(self.layout), list(self.value))
return s
def __nonzero__(self):
"""Instance is nonzero iff at least one of the coefficients
is nonzero.
__nonzero() --> Boolean
"""
nonzeroes = Numeric.greater(Numeric.absolute(self.value), _eps)
if nonzeroes:
return 1
else:
return 0
def __cmp__(self, other):
"""Compares two multivectors.
This is mostly defined for equality testing (within epsilon).
In the case that the two multivectors have different Layouts,
we will raise an error. In the case that they are not equal,
we will compare the tuple represenations of the coefficients
lists just so as to return something valid. Therefore,
inequalities are well-nigh meaningless (since they are
meaningless for multivectors while equality is meaningful).
Oh, how I wish for rich comparisons.
M == N
__cmp__(other) --> -1|0|1
"""
other, mv = self._checkOther(other)
if Numeric.alltrue(Numeric.less(Numeric.absolute( \
self.value - other.value), _eps)):
# equal within epsilon
return 0
else:
return cmp(tuple(self.value), tuple(other.value))
## Geometric Algebraic functions
def isScalar(self):
"""Returns true iff self is a scalar.
isScalar() --> Boolean
"""
indices = range(self.layout.gaDims)
indices.remove(self.layout.gradeList.index(0))
for i in indices:
if abs(self.value[i]) < _eps:
continue
else:
return 0
return 1
def isBlade(self):
"""Returns true if multivector is a blade.
isBlade() --> Boolean
"""
grade = None
for i in range(self.layout.gaDims):
if abs(self.value[i]) > _eps:
if grade is None:
grade = self.layout.gradeList[i]
elif self.layout.gradeList[i] != grade:
return 0
return 1
def grades(self):
"""Return the grades contained in the multivector.
grades() --> PyInt"""
grades = []
for i in range(self.layout.gaDims):
if abs(self.value[i]) > _eps and \
self.layout.gradeList[i] not in grades:
grades.append(self.layout.gradeList[i])
return grades
def laInv(self):
"""Return inverse using a computational linear algebra method
proposed by Christian Perwass.
-1
M
laInv() --> MultiVector
"""
identity = Numeric.zeros((self.layout.gaDims,))
identity[self.layout.gradeList.index(0)] = 1
intermed = Numeric.dot(self.value, self.layout.gmt)
intermed = Numeric.transpose(intermed)
if abs(LA.determinant(intermed)) < _eps:
raise ValueError, "multivector has no inverse"
sol = LA.solve_linear_equations(intermed, identity)
return MultiVector(self.layout, sol)
def normalInv(self):
"""Returns the inverse of itself if M&~M == |M|**2.
-1
M = ~M / (M & ~M)
normalInv() --> MultiVector
"""
Madjoint = ~self
MadjointM = (Madjoint & self)
if MadjointM.isScalar() and abs(MadjointM[()]) > _eps:
# inverse exists
return Madjoint / MadjointM[()]
else:
raise ValueError, "no inverse exists for this multivector"
if LA:
inv = laInv
else:
inv = normalInv
def dual(self):
"""Returns the dual of the multivector.
~ -1
M = M * I where I is the pseudoscalar.
dual() --> MultiVector
"""
return self * self.layout.invPS()
def commutator(self, other):
"""Returns the commutator product of two multivectors.
[M, N] = M X N = (M&N - N&M)/2
commutator(other) --> MultiVector
"""
return ((self & other) - (other & self)) / 2
def join(self, other):
"""Returns the join of two blades.
J = M ^ N iff M, N are blades and M^N != 0
join(other) --> MultiVector
"""
other, mv = self._checkOther(other)
if self.isBlade() and other.isBlade():
J = self ^ other
if J:
return J
else:
raise ValueError, "join does not exist, blades share common factors"
else:
raise ValueError, "not blades"
def meet(self, other, subspace=None):
"""Returns the meet of an r-blade and an s-blade iff r+s >= dims.
Computation is done with respect to a subspace that defaults to
the pseudo scalar if none is given.
-1
M \/i N = (Mi ) * N
meet(other, subspace=None) --> MultiVector
"""
other, mv = self._checkOther(other)
r = self.grades()
s = other.grades()
if len(r) > 1 or len(s) > 1:
raise ValueError, "not blades"
r = r[0]
s = s[0]
if subspace is None:
subspace = self.layout.pseudoScalar()
dims = subspace.grades()[0]
if r + s < dims:
raise ValueError, "r+s < dims; meet not defined"
elif r + s == dims:
# special case; M v N == 0
return MultiVector(self.layout)
else:
return (self & subspace.inv()) * other
def gradeInvol(self):
"""Returns the grade involution of the multivector.
* i
M = Sum[i, dims, (-1) <M> ]
i
gradeInvol() --> MultiVector
"""
signs = Numeric.power(-1, self.layout.gradeList)
newValue = signs * self.value
return MultiVector(self.layout, newValue)
def conjugate(self):
"""Returns the Clifford conjugate (reversion and grade involution).
*
M --> (~M).gradeInvol()
conjugate() --> MultiVector
"""
return (~self).gradeInvol()
def comb(n, k):
"""\
Returns /n\\
\\k/
comb(n, k) --> PyInt"""
def fact(n):
return Numeric.multiply.reduce(range(1, n+1))
return fact(n) / (fact(k) * fact(n-k))
def elements(dims, firstIdx=0):
"""Return a list of tuples representing all 2**dims of blades
in a dims-dimensional GA.
elements(dims, firstIdx=0) --> bladeList"""
indcs = range(firstIdx, firstIdx + dims)
blades = [()]
for k in range(1, dims+1):
# k = grade
if k == 1:
for i in indcs:
blades.append((i,))
continue
curBladeX = indcs[:k]
marker = -2
# curBladeX[marker:] will be replaced with a range each time
# the final element overflows; you'll see
for i in range(comb(dims, k)):
if curBladeX[-1] < firstIdx+dims-1:
# increment last index
blades.append(tuple(curBladeX))
curBladeX[-1] = curBladeX[-1] + 1
else:
if curBladeX[marker:] == range(firstIdx+dims+marker,
firstIdx+dims):
# decrement marker
marker = marker - 1
if marker < -k:
blades.append(tuple(curBladeX))
continue
# replace
blades.append(tuple(curBladeX))
curBladeX[marker:] = range(curBladeX[marker]+1,
curBladeX[marker]+1 - marker)
return blades
def Cl(p, q=0, names=None, firstIdx=0):
"""Returns a Layout and basis blades for the geometric algebra Cl_p,q.
The notation Cl_p,q means that the algebra is p+q dimensional, with
the first p vectors with positive signature and the final q vectors
negative.
Cl(p, q=0, names=None, firstIdx=0) --> Layout, {'name': basisElement, ...}"""
sig = [+1]*p + [-1]*q
bladeList = elements(p+q, firstIdx)
layout = Layout(sig, bladeList, firstIdx=firstIdx, names=names)
blades = bases(layout)
return layout, blades
def bases(layout):
"""Returns a dictionary mapping basis element names to their MultiVector
instances.
bases(layout) --> {'name': baseElement, ...}"""
dict = {}
for i in range(layout.gaDims):
if layout.gradeList[i] != 0:
v = Numeric.zeros((layout.gaDims,))
v[i] = 1
dict[layout.names[i]] = MultiVector(layout, v)
return dict
def randomMV(layout, min=-2.0, max=2.0, grades=None):
"""Random MultiVector with given layout.
Coefficients are between min and max, and if grades is a list of integers,
only those grades will be non-zero. Uses the RandomArray module.
randomMV(layout, min=-2.0, max=2.0, grades=None)"""
from RandomArray import uniform
if grades is None:
return MultiVector(layout, uniform(min, max, (layout.gaDims,)))
else:
newValue = Numeric.zeros((layout.gaDims,)).astype(Numeric.Float)
for i in range(layout.gaDims):
if layout.gradeList[i] in grades:
newValue[i] = uniform(min, max)
return MultiVector(layout, newValue)
def pretty():
"""Makes repr(M) default to pretty-print.
pretty()"""
global _pretty
_pretty = 1
def ugly():
"""Makes repr(M) default to eval-able representation.
ugly()"""
global _pretty
_pretty = 0
def eps(newEps):
"""Set the epsilon for float comparisons.
eps(newEps)"""
global _eps
_eps = newEps
------------%<------------------------------------------
--
Robert Kern
kern at caltech.edu
"In the fields of hell where the grass grows high
Are the graves of dreams allowed to die."
-- Richard Harter
More information about the Python-list
mailing list