# Source code for galgebra.mv

"""
Multivector and Linear Multivector Differential Operator
"""

import itertools
import copy
import numbers
import operator
# from compiler.ast import flatten
# https://stackoverflow.com/questions/16176742/python-3-replacement-for-deprecated-compiler-ast-flatten-function
from .utils import flatten
from operator import itemgetter, mul, add
from itertools import combinations
from sympy import Symbol, Function, S, expand, Add, Mul, Pow, Basic, \
sin, cos, sinh, cosh, sqrt, trigsimp, expand, \
simplify, diff, Rational, Expr, Abs, collect, combsimp
from sympy import exp as sympy_exp
from sympy import N as Nsympy
from . import printer
from . import metric
from . import utils
from .printer import ZERO_STR
import sys
from functools import reduce, cmp_to_key

ONE = S(1)
ZERO = S(0)
HALF = Rational(1, 2)

half = Rational(1, 2)

modules = \
"""
from sympy import symbols, sin, Function
from mv import Mv
from ga import Ga, half
from printer import Eprint, xdvi
from lt import Lt
"""

########################### Multivector Class ##########################

[docs] @staticmethod def setup(ga): """ Set up constant mutilvectors reqired for multivector class for a given geometric algebra, 'ga'. """ Mv.fmt = 1 basis = [Mv(x, ga=ga) for x in ga.basis] I = Mv(ga.iobj, ga=ga) # default pseudoscalar x = Mv('XxXx', 'vector', ga=ga) # testing vectors # return default basis vectors and grad vector if coords defined return I, basis, x
@staticmethod def Format(mode=1): Mv.latex_flg = True Mv.fmt = mode return
[docs] @staticmethod def Mul(A, B, op): """ Function for all types of geometric multiplications called by overloaded operators for *, ^, |, <, and >. """ if not isinstance(A, Mv): A = B.Ga.mv(A) if not isinstance(B, Mv): B = A.Ga.mv(B) if op == '*': return A * B elif op == '^': return A ^ B elif op == '|': return A | B elif op == '<': return A < B elif op == '>': return A > B else: raise ValueError('Operation ' + op + 'not allowed in Mv.Mul!') return
[docs] def collect(self,deep=False): """ # group coeffients of blades of multivector # so there is only one coefficient per grade self.obj = expand(self.obj) if self.is_blade_rep or Mv.Ga.is_ortho: c = self.Ga.blades_lst else: c = self.Ga.bases_lst self.obj = self.obj.collect(c) return self """ coefs, bases = metric.linear_expand(self.obj) obj_dict = {} for (coef, base) in zip(coefs, bases): if base in list(obj_dict.keys()): obj_dict[base] += coef else: obj_dict[base] = coef obj = 0 for base in list(obj_dict.keys()): if deep: obj += collect(obj_dict[base])*base else: obj += obj_dict[base]*base self.obj = obj return(self)
[docs] def is_versor(self): # Test for versor (geometric product of vectors) """ This follows Leo Dorst's test for a versor. Leo Dorst, 'Geometric Algebra for Computer Science,' p.533 Sets self.versor_flg and returns value """ if self.versor_flg is not None: return self.versor_flg self.characterise_Mv() self.versor_flg = False self_rev = self.rev() # see if self*self.rev() is a scalar test = self*self_rev if not test.is_scalar(): return self.versor_flg # see if self*x*self.rev() returns a vector for x an arbitrary vector test = self * self.Ga.XOX * self.rev() self.versor_flg = test.is_vector() return self.versor_flg
def is_zero(self): if self.obj == 0: return True return False def scalar(self): # return scalar part of multivector # as sympy expression return self.Ga.scalar_part(self.obj) def get_grade(self, r): # return r-th grade of multivector as # a multivector return Mv(self.Ga.get_grade(self.obj, r), ga=self.Ga) def components(self): (coefs, bases) = metric.linear_expand(self.obj) bases_lst = self.Ga.blades_lst cb = list(zip(coefs, bases)) cb = sorted(cb, key=lambda x: self.Ga.blades_lst0.index(x[1])) terms = [] for (coef, base) in cb: terms.append(self.Ga.mv(coef * base)) return terms def get_coefs(self, grade): (coefs, bases) = metric.linear_expand(self.obj) bases_lst = self.Ga.blades_lst cb = list(zip(coefs, bases)) cb = sorted(cb, key=lambda x: self.Ga.blades[grade].index(x[1])) (coefs, bases) = list(zip(*cb)) return coefs
[docs] def proj(self, bases_lst): """ Project multivector onto a given list of bases. That is find the part of multivector with the same bases as in the bases_lst. """ bases_lst = [x.obj for x in bases_lst] (coefs, bases) = metric.linear_expand(self.obj) obj = 0 for (coef, base) in zip(coefs, bases): if base in bases_lst: obj += coef * base return Mv(obj, ga=self.Ga)
def dual(self): mode = self.Ga.dual_mode_value sign = S(1) if '-' in mode: sign = -sign if 'Iinv' in mode: I = self.Ga.i_inv else: I = self.Ga.i if mode[0] == '+' or mode[0] == '-': return sign * I * self else: return sign * self * I def even(self): # return even parts of multivector return Mv(self.Ga.even_odd(self.obj, True), ga=self.Ga) def odd(self): # return odd parts of multivector return Mv(self.Ga.even_odd(self.obj, False), ga=self.Ga) def rev(self): self = self.blade_rep() return Mv(self.Ga.reverse(self.obj), ga=self.Ga) __invert__ = rev # allow ~x to call x.rev() def diff(self, coord): Dself = Mv(ga=self.Ga) if self.Ga.coords is None: Dself.obj = diff(self.obj, coord) return Dself elif coord not in self.Ga.coords: if self.Ga.par_coords is None: Dself.obj = diff(self.obj, coord) elif coord not in self.Ga.par_coords: Dself.obj = diff(self.obj, coord) else: Dself.obj = diff(self.obj, coord) for x_coord in self.Ga.coords: f = self.Ga.par_coords[x_coord] if f != S(0): tmp1 = self.Ga.pDiff(self.obj, x_coord) tmp2 = diff(f, coord) Dself.obj += tmp1 * tmp2 Dself.characterise_Mv() return Dself else: Dself.obj = self.Ga.pDiff(self.obj, coord) Dself.characterise_Mv() return Dself def pdiff(self, var): return Mv(self.Ga.pDiff(self.obj, var), ga=self.Ga)
[docs] def Grad(self, coords, mode='*', left=True): """ Returns various derivatives (*,^,|,<,>) of multivector functions with respect to arbitrary coordinates, 'coords'. This would be used where you have a multivector function of both the basis coordinate set and and auxilliary coordinate set. Consider for example a linear transformation in which the matrix coefficients depend upon the manifold coordinates, but the vector being transformed does not and you wish to take the divergence of the linear transformation with respect to the linear argument. """ return Mv(self.Ga.Diff(self, mode, left, coords=coords), ga=self.Ga)
[docs] def exp(self, hint='-'): # Calculate exponential of multivector """ Only works if square of multivector is a scalar. If square is a number we can determine if square is > or < zero and hence if one should use trig or hyperbolic functions in expansion. If square is not a number use 'hint' to determine which type of functions to use in expansion """ self = self.blade_rep() self_sq = self * self if self_sq.is_scalar(): sq = simplify(self_sq.obj) # sympy expression for self**2 if sq == S(0): # sympy expression for self**2 = 0 return self + S(1) (coefs,bases) = metric.linear_expand(self.obj) if len(coefs) == 1: # Exponential of scalar * base base = bases[0] base_Mv = self.Ga.mv(base) base_sq = (base_Mv*base_Mv).scalar() if hint == '-': # base^2 < 0 base_n = sqrt(-base_sq) return self.Ga.mv(cos(base_n*coefs[0]) + sin(base_n*coefs[0])*(bases[0]/base_n)) else: # base^2 > 0 base_n = sqrt(base_sq) return self.Ga.mv(cosh(base_n*coefs[0]) + sinh(base_n*coefs[0])*(bases[0]/base_n)) if sq.is_number: # Square is number, can test for sign if sq > S(0): norm = sqrt(sq) value = self.obj / norm tmp = Mv(cosh(norm) + sinh(norm) * value, ga=self.Ga) tmp.is_blade_rep = True return tmp else: norm = sqrt(-sq) value = self.obj / norm tmp = Mv(cos(norm) + sin(norm) * value, ga=self.Ga) tmp.is_blade_rep = True return tmp else: if hint == '+': norm = simplify(sqrt(sq)) value = self.obj / norm tmp = Mv(cosh(norm) + sinh(norm) * value, ga=self.Ga) tmp.is_blade_rep = True return tmp else: norm = simplify(sqrt(-sq)) value = self.obj / norm obj = cos(norm) + sin(norm) * value tmp = Mv(cos(norm) + sin(norm) * value, ga=self.Ga) tmp.is_blade_rep = True return tmp else: raise ValueError('"' + str(self) + '**2" is not a scalar in exp.')
def set_coef(self, igrade, ibase, value): if self.blade_rep: base = self.Ga.blades[igrade][ibase] else: base = self.Ga.bases[igrade][ibase] (coefs, bases) = metric.linear_expand(self.obj) bases_lst = list(bases) # python 2.5 if base in bases: self.obj += (value - coefs[bases_lst.index(base)]) * base else: self.obj += value * base return
[docs] def Fmt(self, fmt=1, title=None): """ Set format for printing of multivectors - fmt = 1 - One multivector per line fmt = 2 - One grade per line fmt = 3 - one base per line Usage for multivector A example is - A.Fmt('2','A') output is 'A = '+str(A) with one grade per line. Works for both standard printing and for latex. """ if printer.GaLatexPrinter.latex_flg: printer.GaLatexPrinter.prev_fmt = printer.GaLatexPrinter.fmt printer.GaLatexPrinter.fmt = fmt else: printer.GaPrinter.prev_fmt = printer.GaPrinter.fmt printer.GaPrinter.fmt = fmt if title is not None: self.title = title if printer.isinteractive(): return self if Mv.latex_flg: latex_str = printer.GaLatexPrinter.latex(self) printer.GaLatexPrinter.fmt = printer.GaLatexPrinter.prev_fmt if title is not None: return title + ' = ' + latex_str else: return latex_str else: s = str(self) printer.GaPrinter.fmt = printer.GaPrinter.prev_fmt if title is not None: return title + ' = ' + s else: return s return
def _repr_latex_(self): latex_str = printer.GaLatexPrinter.latex(self) if r'\begin{align*}' not in latex_str: if self.title is None: latex_str = r'\begin{equation*} ' + latex_str + r' \end{equation*}' else: latex_str = r'\begin{equation*} ' + self.title + ' = ' + latex_str + r' \end{equation*}' else: if self.title is not None: latex_str = latex_str.replace('&',' ' + self.title + ' =&',1) return latex_str def norm2(self): reverse = self.rev() product = self * reverse if product.is_scalar(): return product.scalar() else: raise TypeError('"(' + str(product) + ')**2" is not a scalar in norm2.')
[docs] def norm(self, hint='+'): """ If A is a multivector and A*A.rev() is a scalar then A.norm() = sqrt(Abs(A*A.rev())) The problem in simplifing the norm is that if A is symbolic you don't know if A*A.rev() is positive or negative. The use of the hint argument is as follows: hint A.norm() '+' sqrt(A*A.rev()) '-' sqrt(-A*A.rev()) '0' sqrt(Abs(A*A.rev())) The default hint='+' is correct for vectors in a Euclidean vector space. For bivectors in a Euclidean vector space use hint='-'. In a mixed signature space all bets are off for the norms of symbolic expressions. """ reverse = self.rev() product = self * reverse if product.is_scalar(): product = product.scalar() if product.is_number: if product >= S(0): return sqrt(product) else: return sqrt(-product) else: if hint == '+': return metric.square_root_of_expr(product) elif hint == '-': return metric.square_root_of_expr(-product) else: return sqrt(Abs(product)) else: raise TypeError('"(' + str(product) + ')" is not a scalar in norm.')
__abs__=norm # allow abs(x) to call z.norm() def inv(self): if self.is_scalar(): # self is a scalar return self.Ga.mv(S(1)/self.obj) self_sq = self * self if self_sq.is_scalar(): # self*self is a scalar """ if self_sq.scalar() == S(0): raise ValueError('!!!!In multivector inverse, A*A is zero!!!!') """ return (S(1)/self_sq.obj)*self self_rev = self.rev() self_self_rev = self * self_rev if(self_self_rev.is_scalar()): # self*self.rev() is a scalar """ if self_self_rev.scalar() == S(0): raise ValueError('!!!!In multivector inverse A*A.rev() is zero!!!!') """ return (S(1)/self_self_rev.obj) * self_rev raise TypeError('In inv() for self =' + str(self) + 'self, or self*self or self*self.rev() is not a scalar') def func(self, fct): # Apply function, fct, to each coefficient of multivector (coefs, bases) = metric.linear_expand(self.obj) s = S(0) for (coef, base) in zip(coefs, bases): s += fct(coef) * base fct_self = Mv(s, ga=self.Ga) fct_self.characterise_Mv() return fct_self def trigsimp(self): return self.func(trigsimp) def simplify(self, modes=simplify): (coefs, bases) = metric.linear_expand(self.obj) obj = S(0) if isinstance(modes, list) or isinstance(modes, tuple): for (coef, base) in zip(coefs, bases): for mode in modes: coef = mode(coef) obj += coef * base else: for (coef, base) in zip(coefs, bases): obj += modes(coef) * base self.obj = obj return self def subs(self, d): # For each scalar coef of the multivector apply substitution argument d (coefs, bases) = metric.linear_expand(self.obj) obj = S(0) for (coef, base) in zip(coefs, bases): obj += coef.subs(d) * base #self.obj = obj #return self return self.Ga.mv(obj) def expand(self): coefs,bases = metric.linear_expand(self.obj) new_coefs = [] for coef in coefs: new_coefs.append(expand(coef)) obj = 0 for coef,base in zip(new_coefs,bases): obj += coef * base self.obj = obj return self def list(self): (coefs, bases) = metric.linear_expand(self.obj) indexes = [] key_coefs = [] for (coef, base) in zip(coefs, bases): if base in self.Ga.basis: index = self.Ga.basis.index(base) key_coefs.append((coef, index)) indexes.append(index) for index in self.Ga.n_range: if index not in indexes: key_coefs.append((S(0), index)) key_coefs = sorted(key_coefs, key=itemgetter(1)) coefs = [x[0] for x in key_coefs] return coefs def grade(self, r=0): return self.get_grade(r)
[docs]def compare(A,B): """ Determine is B = c*A where c is a scalar. If true return c otherwise return 0. """ if isinstance(A, Mv) and isinstance(B, Mv): Acoefs, Abases = metric.linear_expand(A.obj) Bcoefs, Bbases = metric.linear_expand(B.obj) if len(Acoefs) != len(Bcoefs): return 0 if Abases != Bbases: return 0 if Bcoefs[0] != 0 and Abases[0] == Bbases[0]: c = simplify(Acoefs[0]/Bcoefs[0]) print('c =',c) else: return 0 for acoef,abase,bcoef,bbase in zip(Acoefs[1:],Abases[1:],Bcoefs[1:],Bbases[1:]): print(acoef,'\n',abase,'\n',bcoef,'\n',bbase) if bcoef != 0 and abase == bbase: print('c-a/b =',simplify(c-(acoef/bcoef))) if simplify(acoef/bcoef) != c: return 0 else: pass else: return 0 return c else: raise TypeError('In compare both arguments are not multivectors\n')
################ Scalar Partial Differential Operator Class ############
[docs]class Sdop(object): """ Scalar differential operator is of the form (Einstein summation) D = c_{i}*D_{i} where the c_{i}'s are scalar coefficient (they could be functions) and the D_{i}'s are partial differential operators. """ init_slots = {'ga': (None, 'Associated geometric algebra')} ga = None str_mode = False @staticmethod def setGa(ga): Sdop.ga = ga Pdop.setGa(ga) return def TSimplify(self): new_terms = [] for (coef, pdiff) in self.terms: new_terms.append((metric.Simp.apply(coef), pdiff)) self.terms = new_terms return
[docs] @staticmethod def consolidate_coefs(sdop): """ Remove zero coefs and consolidate coefs with repeated pdiffs. """ if isinstance(sdop, Sdop): terms = sdop.terms else: terms = sdop new_coefs = [] new_pdiffs = [] for (coef, pd) in terms: if coef != S(0): if pd in new_pdiffs: index = new_pdiffs.index(pd) new_coefs[index] += coef else: new_coefs.append(coef) new_pdiffs.append(pd) new_terms = list(zip(new_coefs, new_pdiffs)) if isinstance(sdop, Sdop): return Sdop(new_terms, ga=sdop.Ga) else: return new_terms