From: Conrad Schneiker Date: 2000-08-27T10:10:02+09:00 Subject: [ruby-talk:4595] Re: ANN: Clifford Algebras in Python FYI -- in case someone is interested in doing a Ruby variant.... Even if not, the URLs cited below will be interesting for many of you with an interest in math and physics. In my non-expert opinion, Geometric calculus:multidimensional math =~= Ruby:OOP kern@caltech.edu wrote: > > 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@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@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 = _|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 = _|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) --> > 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) ] > 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@caltech.edu > > "In the fields of hell where the grass grows high > Are the graves of dreams allowed to die." > -- Richard Harter -- Conrad Schneiker (This note is unofficial and subject to improvement without notice.)