Source code for sympy.galgebra.GA

#!/usr/bin/python

"""
The module symbolicGA implements symbolic Geometric Algebra in python.
The relevant references for this module are:

    1. "Geometric Algebra for Physicists" by C. Doran and A. Lazenby,
       Cambridge University Press, 2003.

    2. "Geometric Algebra for Computer Science" by Leo Dorst,
       Daniel Fontijne, and Stephen Mann, Morgan Kaufmann Publishers, 2007.

    3. SymPy Tutorial, http://docs.sympy.org/
"""
import sys
import numpy, sympy
import re as regrep
import sympy.galgebra.latex_ex
from sympy.core.decorators import deprecated

NUMPAT = regrep.compile( '([\-0-9])|([\-0-9]/[0-9])')
"""Re pattern for rational number"""

ZERO = sympy.Rational(0)
ONE  = sympy.Rational(1)
TWO  = sympy.Rational(2)
HALF = sympy.Rational(1,2)

from sympy.core import Pow as pow_type
from sympy import Abs as abs_type
from sympy.core import Mul as mul_type
from sympy.core import Add as add_type

@sympy.vectorize(0)
def substitute_array(array,*args):
    return(array.subs(*args))

[docs]def is_quasi_unit_numpy_array(array): """ Determine if a array is square and diagonal with entries of +1 or -1. """ shape = numpy.shape(array) if len(shape) == 2 and (shape[0] == shape[1]): n = shape[0] ix = 0 while ix < n: iy = 0 while iy < ix: if array[ix][iy] != ZERO: return(False) iy += 1 if sympy.Abs(array[ix][ix]) != ONE: return(False) ix += 1 return(True) else: return(False)
@deprecated(value="This function is no longer needed.", issue=3379, deprecated_since_version="0.7.2") def set_main(main_program): pass def plist(lst): if type(lst) == list: for x in lst: plist(x) else: sys.stderr.write(lst+'\n') return
[docs]def numeric(num_str): """ Returns rational numbers compatible with symbols. Input is a string representing a fraction or integer or a simple integer. """ if type(num_str) == int: a = num_str b = 1 else: tmp = num_str.split('/') if len(tmp) == 1: a = int(tmp[0]) b = 1 else: a = int(tmp[0]) b = int(tmp[1]) return(sympy.Rational(a,b))
[docs]def collect(expr,lst): """ Wrapper for sympy.collect. """ lst = MV.scalar_to_symbol(lst) return(sympy.collect(expr,lst))
def sqrt(expr): return(sympy.sqrt(expr))
[docs]def isint(a): """ Test for integer. """ return(type(a) == int)
[docs]def make_null_array(n): """ Return list of n empty lists. """ a = [] for i in range(n): a.append([]) return(a)
[docs]def test_int_flgs(lst): """ Test if all elements in list are 0. """ for i in lst: if i: return(1) return(0)
[docs]def comb(N,P): """ Calculates the combinations of the integers [0,N-1] taken P at a time. The output is a list of lists of integers where the inner lists are the different combinations. Each combination is sorted in ascending order. """ def rloop(n,p,combs,comb): if p: for i in range(n): newcomb = comb+[i] np = p-1 rloop(i,np,combs,newcomb) else: combs.append(comb) combs = [] rloop(N,P,combs,[]) for comb in combs: comb.sort() return(combs)
[docs]def diagpq(p, q=0): """ Return string equivalent metric tensor for signature (p, q). """ n = p + q D = [] rn = list(range(n)) for i in rn: for j in rn: if i ==j: if i < p: D.append('1 ') else: D.append('-1 ') else: D.append('0 ') return ','.join(D)
[docs]def make_scalars(symnamelst): """ make_scalars takes a string of symbol names separated by blanks and converts them to MV scalars and returns a list of the symbols. """ symlst = sympy.symbols(symnamelst) scalar_lst = [] for s in symlst: tmp = MV(s,'scalar') scalar_lst.append(tmp) return(scalar_lst)
@deprecated(useinstead="sympy.symbols()", issue=3379, deprecated_since_version="0.7.2") def make_symbols(symnamelst): return sympy.symbols(symnamelst)
[docs]def israt(numstr): """ Test if string represents a rational number. """ global NUMPAT if NUMPAT.match(numstr): return(1) return(0)
[docs]def dualsort(lst1, lst2): """ Inplace dual sort of lst1 and lst2 keyed on sorted lst1. """ _indices = list(range(len(lst1))) _indices.sort(key=lst2.__getitem__) lst1[:] = list(map(lst1.__getitem__, _indices)) lst2[:] = list(map(lst2.__getitem__, _indices)) return
[docs]def cp(A,B): """ Calculates the commutator product (A*B-B*A)/2 for the objects A and B. """ return(HALF*(A*B-B*A))
[docs]def reduce_base(k,base): """ If base is a list of sorted integers [i_1,...,i_R] then reduce_base sorts the list [k,i_1,...,i_R] and calculates whether an odd or even number of permutations is required to sort the list. The sorted list is returned and +1 for even permutations or -1 for odd permutations. """ if k in base: return(0,base) grade = len(base) if grade == 1: if k < base[0]: return(1,[k,base[0]]) else: return(-1,[base[0],k]) ilo = 0 ihi = grade-1 if k < base[0]: return(1,[k]+base) if k > base[ihi]: if grade%2 == 0: return(1,base+[k]) else: return(-1,base+[k]) imid = ihi+ilo if grade == 2: return(-1,[base[0],k,base[1]]) while True: if ihi-ilo == 1: break if base[imid] > k: ihi = imid else: ilo = imid imid = (ilo+ihi)/2 if ilo%2 == 1: return(1,base[:ihi]+[k]+base[ihi:]) else: return(-1,base[:ihi]+[k]+base[ihi:])
[docs]def sub_base(k,base): """ If base is a list of sorted integers [i_1,...,i_R] then sub_base returns a list with the k^th element removed. Note that k=0 removes the first element. There is no test to see if k is in the range of the list. """ n = len(base) if n == 1: return([]) if n == 2: if k == base[0]: return([base[1]]) else: return([base[0]]) return(base[:k]+base[k+1:])
[docs]def magnitude(vector): """ Calculate magnitude of vector containing trig expressions and simplify. This is a hack because of way he sign of magsq is determined and because of the way that absolute values are removed. """ magsq = sympy.expand((vector|vector)()) magsq = sympy.trigsimp(magsq,deep=True,recursive=True) #print magsq magsq_str = sympy.galgebra.latex_ex.LatexPrinter()._print(magsq) if magsq_str[0] == '-': magsq = -magsq mag = unabs(sqrt(magsq)) #print mag return(mag)
[docs]def LaTeX_lst(lst,title=''): """ Output a list in LaTeX format. """ if title != '': sympy.galgebra.latex_ex.LaTeX(title) for x in lst: sympy.galgebra.latex_ex.LaTeX(x) return
[docs]def unabs(x): """ Remove absolute values from expressions so a = sqrt(a**2). This is a hack. """ if type(x) == mul_type: y = unabs(x.args[0]) for yi in x.args[1:]: y *= unabs(yi) return(y) if type(x) == pow_type: if x.args[1] == HALF and type(x.args[0]) == add_type: return(x) y = 1/unabs(x.args[0]) return(y) if len(x.args) == 0: return(x) if type(x) == abs_type: return(x.args[0]) return(x)
def function_lst(fstr,xtuple): sys.stderr.write(fstr+'\n') fct_lst = [] for xstr in fstr.split(): f = sympy.Function(xstr)(*xtuple) fct_lst.append(f) return(fct_lst)
[docs]def vector_fct(Fstr,x): """ Create a list of functions of arguments x. One function is created for each variable in x. Fstr is a string that is the base name of each function while each function in the list is given the name Fstr+'__'+str(x[ix]) so that if Fstr = 'f' and str(x[1]) = 'theta' then the LaTeX output of the second element in the output list would be 'f^{\\theta}'. """ nx = len(x) Fvec = [] for ix in range(nx): ftmp = sympy.Function(Fstr+'__'+sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[ix]))(*tuple(x)) Fvec.append(ftmp) return(Fvec)
def print_lst(lst): for x in lst: print(x) return
[docs]def normalize(elst,nname_lst): """ Normalize a list of vectors and rename the normalized vectors. 'elist' is the list (or array) of vectors to be normalized and nname_lst is a list of the names for the normalized vectors. The function returns the numpy arrays enlst and mags containing the normalized vectors (enlst) and the magnitudes of the original vectors (mags). """ i = 0 mags = numpy.array(MV.n*[ZERO],dtype=numpy.object) enlst= numpy.array(MV.n*[ZERO],dtype=numpy.object) for (e,nname) in zip(elst,nname_lst): emag = magnitude(e) emaginv = 1/emag mags[i] = emag enorm = emaginv*e enorm.name = nname enlst[i] = enorm i += 1 return(enlst,mags)
def build_base(base_index,base_vectors,reverse=False): base = base_vectors[base_index[0]] if len(base_index) > 1: for i in base_index[1:]: base = base^base_vectors[i] if reverse: base = base.rev() return(base) class MV(object): is_setup = False basislabel_lst = 0 curvilinear_flg = False coords = None @staticmethod def pad_zeros(value,n): """ Pad list with zeros to length n. If length is > n truncate list. Return padded list. """ nvalue = len(value) if nvalue < n: value = value+(n-nvalue)*[ZERO] if nvalue > n: value = value[:n] return(value) @staticmethod def define_basis(basis): """ Calculates all the MV static variables needed for basis operations. See reference 5 section 2. """ MV.vbasis = basis MV.vsyms = sympy.symbols(MV.vbasis) MV.n = len(MV.vbasis) MV.nrg = list(range(MV.n)) MV.n1 = MV.n+1 MV.n1rg = list(range(MV.n1)) MV.npow = 2**MV.n MV.index = list(range(MV.n)) MV.gabasis = [[]] MV.basis = (MV.n+1)*[0] MV.basislabel = (MV.n+1)*[0] MV.basis[0] = [] MV.basislabel[0] = '1' MV.basislabel_lst = [['1']] MV.nbasis = numpy.array((MV.n+1)*[1],dtype=numpy.object) for igrade in range(1,MV.n+1): tmp = comb(MV.n,igrade) MV.gabasis += [tmp] ntmp = len(tmp) MV.nbasis[igrade] = ntmp MV.basis[igrade] = tmp gradelabels = [] gradelabel_lst = [] for i in range(ntmp): tmp_lst = [] bstr = '' for j in tmp[i]: bstr += MV.vbasis[j] tmp_lst.append(MV.vbasis[j]) gradelabel_lst.append(tmp_lst) gradelabels.append(bstr) MV.basislabel_lst.append(gradelabel_lst) MV.basislabel[igrade] = gradelabels MV.basis_map = [{'':0}] igrade = 1 while igrade <= MV.n: tmpdict = {} bases = MV.gabasis[igrade] ibases = 0 for base in bases: tmpdict[str(base)] = ibases ibases += 1 MV.basis_map.append(tmpdict) igrade += 1 if MV.debug: print('basis strings =',MV.vbasis) print('basis symbols =',MV.vsyms) print('basis labels =',MV.basislabel) print('basis =',MV.basis) print('grades =',MV.nbasis) print('index =',MV.index) return @staticmethod def define_metric(metric): """ Calculates all the MV static variables needed for metric operations. See reference 5 section 2. """ if MV.metric_str: MV.g = [] MV.metric = numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object) if metric == '': metric = numpy.array(MV.n*[MV.n*['#']],dtype=numpy.object) for i in MV.index: for j in MV.index: gij = metric[i][j] if israt(gij): MV.metric[i][j] = numeric(gij) else: if gij == '#': if i == j: gij = '('+MV.vbasis[j]+'**2)' else: gij = '('+MV.vbasis[min(i,j)]+'.'+MV.vbasis[max(i,j)]+')' tmp = sympy.Symbol(gij) MV.metric[i][j] = tmp if i <= j: MV.g.append(tmp) else: MV.metric = metric MV.g = [] for row in metric: g_row = [] for col in metric: g_row.append(col) MV.g.append(g_row) if MV.debug: print('metric =',MV.metric) return @staticmethod def define_reciprocal_frame(): """ Calculates unscaled reciprocal vectors (MV.brecp) and scale factor (MV.Esq). The ith scaled reciprocal vector is (1/MV.Esq)*MV.brecp[i]. The pseudoscalar for the set of basis vectors is MV.E. """ if MV.tables_flg: MV.E = MV.bvec[0] MV.brecp = [] for i in range(1,MV.n): MV.E = MV.E^MV.bvec[i] for i in range(MV.n): tmp = ONE if i%2 != 0: tmp = -ONE for j in range(MV.n): if i != j: tmp = tmp^MV.bvec[j] tmp = tmp*MV.E MV.brecp.append(tmp) MV.Esq = (MV.E*MV.E)() MV.Esq_inv = ONE/MV.Esq for i in range(MV.n): MV.brecp[i] = MV.brecp[i]*MV.Esq_inv return @staticmethod def reduce_basis_loop(blst): """ Makes one pass through basis product representation for reduction of representation to normal form. See reference 5 section 3. """ nblst = len(blst) if nblst <= 1: return(1) jstep = 1 while jstep <nblst: istep = jstep-1 if blst[istep] == blst[jstep]: i = blst[istep] if len(blst) >2: blst = blst[:istep]+blst[jstep+1:] else: blst = [] if len(blst) <= 1 or jstep == nblst-1: blst_flg = 0 else: blst_flg = 1 return(MV.metric[i][i],blst,blst_flg) if blst[istep] > blst[jstep]: blst1 = blst[:istep]+blst[jstep+1:] a1 = TWO*MV.metric[blst[jstep]][blst[istep]] blst = blst[:istep]+[blst[jstep]]+[blst[istep]]+blst[jstep+1:] if len(blst1) <= 1: blst1_flg = 0 else: blst1_flg = 1 return(a1,blst1,blst1_flg,blst) jstep +=1 return(1) @staticmethod def reduce_basis(blst): """ Repetitively applies reduce_basis_loop to basis product representation until normal form is realized. See reference 5 section 3. """ if blst == []: blst_coef = [] blst_expand = [] for i in MV.n1rg: blst_coef.append([]) blst_expand.append([]) blst_expand[0].append([]) blst_coef[0].append(ONE) return(blst_coef,blst_expand) blst_expand = [blst] blst_coef = [ONE] blst_flg = [1] while test_int_flgs(blst_flg): for i in range(len(blst_flg)): if blst_flg[i]: tmp = MV.reduce_basis_loop(blst_expand[i]) if tmp ==1: blst_flg[i] = 0 else: if len(tmp) == 3: blst_coef[i] = tmp[0]*blst_coef[i] blst_expand[i] = tmp[1] blst_flg[i] = tmp[2] else: blst_coef[i] = -blst_coef[i] blst_flg[i] = 1 blst_expand[i] = tmp[3] blst_coef.append(-blst_coef[i]*tmp[0]) blst_expand.append(tmp[1]) blst_flg.append(tmp[2]) (blst_coef,blst_expand) = MV.combine_common_factors(blst_coef,blst_expand) return(blst_coef,blst_expand) @staticmethod def combine_common_factors(blst_coef,blst_expand): new_blst_coef = [] new_blst_expand = [] for i in range(MV.n1): new_blst_coef.append([]) new_blst_expand.append([]) nfac = len(blst_coef) for ifac in range(nfac): blen = len(blst_expand[ifac]) new_blst_coef[blen].append(blst_coef[ifac]) new_blst_expand[blen].append(blst_expand[ifac]) for i in range(MV.n1): if len(new_blst_coef[i]) > 1: MV.contract(new_blst_coef[i],new_blst_expand[i]) return(new_blst_coef,new_blst_expand) @staticmethod def contract(coefs,bases): dualsort(coefs,bases) n = len(bases)-1 i = 0 while i < n: j = i+1 if bases[i] == bases[j]: coefs[i] += coefs[j] bases.pop(j) coefs.pop(j) n -= 1 else: i += 1 n = len(coefs) i = 0 while i < n: if coefs[i] == ZERO: coefs.pop(i) bases.pop(i) n -= 1 else: i +=1 return @staticmethod def convert(coefs,bases): mv = MV() mv.bladeflg = 0 for igrade in MV.n1rg: coef = coefs[igrade] base = bases[igrade] if len(coef) > 0: nbases = MV.nbasis[igrade] mv.mv[igrade] = numpy.array(nbases*[ZERO],dtype=numpy.object) nbaserg = list(range(len(base))) for ibase in nbaserg: if igrade > 0: k = MV.basis[igrade].index(base[ibase]) mv.mv[igrade][k] = coef[ibase] else: mv.mv[igrade] = numpy.array([coef[0]],dtype=numpy.object) return(mv) @staticmethod def set_str_format(str_mode=0): MV.str_mode = str_mode return @staticmethod def str_rep(mv): """ Converts internal representation of a multivector to a string for outputting. If lst_mode = 1, str_rep outputs a list of strings where each string contains one multivector coefficient concatenated with the corresponding base or blade symbol. MV.str_mode Effect 0 Print entire multivector on one line (default) 1 Print each grade on a single line 2 Print each base on a single line """ if MV.bladeprint: mv.convert_to_blades() labels = MV.bladelabel else: if not mv.bladeflg: labels = MV.basislabel else: labels = MV.bladelabel mv.compact() if isinstance(mv.mv[0],int): value = '' else: value = (mv.mv[0][0]).__str__() value = value.replace(' ','') dummy = sympy.Dummy('dummy') for igrade in MV.n1rg[1:]: if isinstance(mv.mv[igrade],numpy.ndarray): j = 0 for x in mv.mv[igrade]: if x != ZERO: xstr = (x*dummy).__str__() xstr = xstr.replace(' ','') if xstr[0] != '-' and len(value) > 0: xstr = '+'+xstr if xstr.find('_dummy') < 2 and xstr[-5:] != '_dummy': xstr = xstr.replace('_dummy*','')+'*'+labels[igrade][j] else: xstr = xstr.replace('_dummy',labels[igrade][j]) if MV.str_mode == 2: xstr += '\n' value += xstr j += 1 if MV.str_mode == 1: value += '\n' if value == '': value = '0' #value = value.replace(' ','') value = value.replace('dummy','1') return(value) @staticmethod def xstr_rep(mv,lst_mode=0): """ Converts internal representation of a multivector to a string for outputing. If lst_mode = 1, str_rep outputs a list of strings where each string contains one multivector coefficient concatenated with the corresponding base or blade symbol. """ if lst_mode: outlst = [] if MV.bladeprint: mv.convert_to_blades() labels = MV.bladelabel else: if not mv.bladeflg: labels = MV.basislabel else: labels = MV.bladelabel value = '' for igrade in MV.n1rg: tmp = [] if isinstance(mv.mv[igrade],numpy.ndarray): j = 0 for x in mv.mv[igrade]: if x != ZERO: xstr = x.__str__() if xstr == '+1' or xstr == '1' or xstr == '-1': if xstr == '+1' or xstr == '1': xstr = '+' else: xstr = '-' else: if xstr[0] != '+': xstr = '+('+xstr+')' else: xstr = '+('+xstr[1:]+')' value += xstr+labels[igrade][j] if MV.str_mode and not lst_mode: value += value+'\n' if lst_mode: tmp.append(value) j += 1 if lst_mode: if len(tmp) > 0: outlst.append(tmp) value = '' if not lst_mode: if len(value) > 1 and value[0] == '+': value = value[1:] if len(value) == 0: value = '0' else: value = outlst return(value) @staticmethod def setup(basis,metric='',rframe=False,coords=None,debug=False,offset=0): """ MV.setup initializes the MV class by calculating the static multivector tables required for geometric algebra operations on multivectors. See reference 5 section 2 for details on basis and metric arguments. """ MV.is_setup = True MV.metric_str = False MV.debug = debug MV.bladeprint = 0 MV.tables_flg = 0 MV.str_mode = 0 MV.basisroot = '' MV.index_offset = offset if coords == None: MV.coords = None else: MV.coords = tuple(coords) rframe= True if type(basis) == str: basislst = basis.split() if len(basislst) == 1: MV.basisroot = basislst[0] basislst = [] for coord in coords: basislst.append(MV.basisroot+'_'+str(coord)) MV.define_basis(basislst) if type(metric) == str: MV.metric_str = True if len(metric) > 0: if metric[0] == '[' and metric[-1] == ']': tmps = metric[1:-1].split(',') N = len(tmps) metric = [] itmp = 0 for tmp in tmps: xtmp = N*['0'] xtmp[itmp] = tmp itmp += 1 metric.append(xtmp) else: tmps = metric.split(',') metric = [] for tmp in tmps: xlst = tmp.split() xtmp = [] for x in xlst: xtmp.append(x) metric.append(xtmp) MV.define_metric(metric) MV.multiplication_table() MV.blade_table() MV.inverse_blade_table() MV.tables_flg = 1 isym = 0 MV.bvec = [] for name in MV.vbasis: bvar = MV(value=isym,mvtype='basisvector',mvname=name) bvar.bladeflg = 1 MV.bvec.append(bvar) isym += 1 if rframe: MV.define_reciprocal_frame() MV.I = MV(ONE,'pseudo','I') MV.ZERO = MV() Isq = (MV.I*MV.I)() MV.Iinv = (1/Isq)*MV.I return MV.bvec @staticmethod def set_coords(coords): MV.coords = coords return @staticmethod def scalar_fct(fct_name): """ Create multivector scalar function with name fct_name (string) and independent variables coords (list of variable). Default variables are those associated with each dimension of vector space. """ phi = sympy.Function(fct_name)(*MV.coords) Phi = MV(phi,'scalar') Phi.name = fct_name return(Phi) @staticmethod def vector_fct(fct_name,vars=''): """ Create multivector vector function with name fct_name (string) and independent variables coords (list of variable). Default variables are those associated with each dimension of vector space. """ if isinstance(vars,str): Acoefs = vector_fct(fct_name,MV.coords) else: Acoefs =numpy.array(MV.n*[ZERO],dtype=numpy.object) x = MV.coords if isinstance(vars,sympy.core.symbol.Symbol): for icoef in MV.nrg: Acoefs[icoef] = sympy.Function(fct_name+'__'+\ sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[icoef]))(vars) else: for icoef in MV.nrg: Acoefs[icoef] = sympy.Function(fct_name+'__'+\ sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[icoef]))(*tuple(vars)) A = MV(Acoefs,'vector',fct_name) return(A) @staticmethod def rebase(x,coords,base_name='',debug=False,debug_level=0): """ Define curvilinear coordinates for previously defined vector (multivector) space (MV.setup has been run) with position vector, x, that is a vector function of the independent coordinates, coords (list of sympy variables equal in length to dimension of vector space), and calculate: 1. Frame (basis) vectors 2. Normalized frame (basis) vectors. 3. Metric tensor 4. Reciprocal frame vectors 5. Reciprocal metric tensor 6. Connection multivectors The basis vectors are named with the base_name (string) and a subscript derived from the name of each coordinate. So that if the base name is 'e' and the coordinated are [r,theta,z] the variable names of the frame vectors would be e_r, e_theta, and e_z. For LaTeX output the names of the frame vectors would be e_{r}, e_{\theta}, and e_{z}. Everything needed to compute the geometric, outer, and inner derivatives of multivector functions in curvilinear coordinates is calculated. If debug is True all the quantities in the above list are output in LaTeX format. Currently rebase works with cylindrical and spherical coordinates in any dimension. The limitation is the ability to automatically simplify complex sympy expressions generated while calculating the quantities in the above list. This is why the debug option is included. The debug_level can equal 0,1,2, or 3 and determines how far in the list to calculate (input 0 to do the entire list) while debugging. """ #Form root names for basis, reciprocal basis, normalized basis, and normalized reciprocal basis if base_name == '': base_name = MV.basisroot+'prm' LaTeX_base = sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(base_name) bm = '\\bm{'+LaTeX_base+'}' bmhat = '\\hat{'+bm+'}' bstr = bmhat+'_{[i_{1},\dots, i_{R}]}' base_name += 'bm' base_name_hat = base_name+'hat' base_name_lst = [] nbase_name_lst = [] rbase_name_lst = [] rnbase_name_lst = [] coords_lst = [] for coord in coords: coord_str = sympy.galgebra.latex_ex.LatexPrinter.str_basic(coord) coords_lst.append(coord_str) base_name_lst.append(base_name+'_'+coord_str) rbase_name_lst.append(base_name+'__'+coord_str) nbase_name_lst.append(base_name_hat+'_'+coord_str) rnbase_name_lst.append(base_name_hat+'__'+coord_str) if not (MV.n == len(coords) == len(base_name_lst)): print('rebaseMV inputs not congruent:') print('MV.n =',MV.n) print('coords =',coords) print('bases =',base_name) sys.exit(1) if isinstance(x,MV): #Calculate basis vectors from derivatives of position vector x bases = numpy.array(MV.n*[ZERO],dtype=numpy.object) i = 0 for coord in coords: ei = x.diff(coords[i]) ei.set_name(base_name_lst[i]) bases[i] = ei i += 1 #Calculate normalizee basis vectors and basis vector magnitudes if debug: print('Coordinate Generating Vector') print(x) print('Basis Vectors') for base in bases: print(base) else: #Input basis vectors as N vector fields bases = x for (base,name) in zip(bases,base_name_lst): base.set_name(name) if debug: print('Basis Vectors') for base in bases: print(base) if debug_level == 1: return if debug_level == 1: return #Calculate normalized basis vectors and magnitudes of #unormalized basis vectors (nbases,mags) = normalize(bases,nbase_name_lst) if debug: print('Magnitudes') print('\\abs{'+LaTeX_base+'_{i}} = ',mags) print('Normalized Basis Vectors') for nbase in nbases: print(nbase) g = numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object) for irow in MV.nrg: for icol in MV.nrg: magsq = sympy.expand((nbases[irow]|nbases[icol])()) g[irow][icol] = sympy.simplify(sympy.trigsimp(magsq,deep=True,recursive=True)) if debug: print('Metric $\\hat{g}_{ij} = \\hat{'+LaTeX_base+'}_{i}\\cdot \\hat{'+\ LaTeX_base+'}_{j}$') print(r'\hat{g}_{ij} =',sympy.galgebra.latex_ex.LaTeX(g)) if debug_level == 2: return #Calculate reciprocal normalized basis vectors rnbases = [] if is_quasi_unit_numpy_array(g): ibasis = 0 while ibasis < MV.n: base = g[ibasis][ibasis]*nbases[ibasis] base.set_name(rnbase_name_lst[ibasis]) base.simplify() base.trigsimp() rnbases.append(base) ibasis += 1 else: rnbases = reciprocal_frame(nbases,rnbase_name_lst) ibase = 0 for base in rnbases: base.simplify() base.trigsimp() rnbases[ibase] = base ibase += 1 if debug: if debug_level != 0: sympy.galgebra.latex_ex.MV_format(1) print('Reciprocal Normalized Basis Vectors') for rnbase in rnbases: print(rnbase) if debug_level == 3: return #Calculate components of inverse vectors Acoef = [] for ibasis in MV.nrg: evec = numpy.array(MV.n*[ZERO],dtype=numpy.object) for jbasis in MV.nrg: evec[jbasis] = (MV.bvec[ibasis]|rnbases[jbasis])() Acoef.append(evec) #Calculat metric tensors gr = numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object) for irow in MV.nrg: for icol in MV.nrg: magsq = sympy.expand((rnbases[irow]|rnbases[icol])()) gr[irow][icol] = sympy.simplify(sympy.trigsimp(magsq,deep=True,recursive=True)) if debug: print('Metric $\\hat{g}^{ij} = \\hat{'+LaTeX_base+'}^{i}\\cdot \\hat{'+\ LaTeX_base+'}^{j}$') print(r'\hat{g}^{ij} =',sympy.galgebra.latex_ex.LaTeX(gr)) if debug_level == 4: return #Calculate bases and reciprocal bases for curvilinear mulitvectors MV_bases = [[ONE]] MV_rbases = [[ONE]] igrade = 1 while igrade <= MV.n: base_index = MV.gabasis[igrade] grade_bases = [] rgrade_bases = [] for index in base_index: base = build_base(index,nbases) base.simplify() base.trigsimp() rbase = build_base(index,rnbases,True) rbase.simplify() rbase.trigsimp() grade_bases.append(base) rgrade_bases.append(rbase) igrade += 1 MV_bases.append(grade_bases) MV_rbases.append(rgrade_bases) #Calculate connection multivectors for geometric derivative MV_connect = [[ZERO]] igrade = 1 while igrade <= MV.n: grade_connect = [] ibase = 0 for base in MV_bases[igrade]: sum = MV() itheta = 0 for (theta,etheta) in zip(coords,rnbases): psum = (1/mags[itheta])*etheta*base.diff(theta) psum.trigsimp() sum += psum itheta += 1 sum.simplify() sum.trigsimp() grade_connect.append(sum) ibase += 1 MV_connect.append(grade_connect) igrade += 1 if debug: print('Curvilinear Bases: $'+bstr+' = '+bmhat+'_{i_{1}}\\W\\dots\\W'+bmhat+'_{i_{R}}$') igrade = 1 for grade in MV_bases[1:]: ibase = 0 for base in grade: index = MV.gabasis[igrade][ibase] sub_str = '' for i in index: sub_str += sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[i]) base_str = bmhat+'_{['+sub_str+']} = ' print(base_str,base) ibase += 1 igrade += 1 if debug_level == 5: return #Calculate representation of connection multivectors in curvilinear system MV_Connect = [[ZERO]] igrade = 1 while igrade <= MV.n: grade_connect = [] ibase = 0 nbase = len(MV_bases[igrade]) if igrade < MV.n: ibase = 0 p1base = len(MV_bases[igrade+1]) m1base = len(MV_bases[igrade-1]) while ibase < nbase: Cm1 = numpy.array(m1base*[ZERO],dtype=numpy.object) Cp1 = numpy.array(p1base*[ZERO],dtype=numpy.object) C = MV_connect[igrade][ibase] if igrade == 1: X = C(0) else: X = MV() jbase = 0 while jbase < m1base: Cm1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade-1][jbase],C))(),deep=True,recursive=True) jbase += 1 jbase = 0 while jbase < p1base: Cp1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade+1][jbase],C))(),deep=True,recursive=True) jbase += 1 X += MV((igrade-1,Cm1),'grade')+MV((igrade+1,Cp1),'grade') X.simplify() X.trigsimp() grade_connect.append(X) ibase += 1 else: ibase = 0 m1base = len(MV_bases[igrade-1]) while ibase < nbase: Cm1 = numpy.array(m1base*[ZERO],dtype=numpy.object) C = MV_connect[igrade][ibase] jbase = 0 while jbase < m1base: Cm1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade-1][jbase],C))(),deep=True,recursive=True) jbase += 1 X = MV() X.mv[MV.n-1] = Cm1 X.simplify() X.trigsimp() grade_connect.append(X) ibase += 1 MV_Connect.append(grade_connect) igrade += 1 base_str = '' for coord in coords: base_str += base_name+'_'+sympy.galgebra.latex_ex.LatexPrinter.str_basic(coord)+' ' base_str = base_str[:-1] old_names = MV.vbasis MV.setup(base_str,g,True,coords) MV.curvilinear_flg = True MV.Connect = MV_Connect sympy.galgebra.latex_ex.LatexPrinter.latex_bases() MV.Rframe = numpy.array(MV.n*[ZERO],dtype=numpy.object) ibasis = 0 while ibasis < MV.n: base = MV() jbasis = 0 while jbasis < MV.n: base.add_in_place(gr[ibasis][jbasis]*MV.bvec[jbasis]) jbasis += 1 base.scalar_mul_inplace(1/mags[ibasis]) MV.Rframe[ibasis] = base ibasis += 1 MV.org_basis = [] for ibasis in MV.nrg: evec = MV(Acoef[ibasis],'vector',old_names[ibasis]) MV.org_basis.append(evec) if MV.coords[0] == sympy.Symbol('t'): MV.dedt = [] # I can't fix this no name error because I have no clue what the # original writer was trying to do. for coef in dedt_coef: MV.dedt.append(MV(coef,'vector')) else: MV.dedt = None if debug: print('Representation of Original Basis Vectors') for evec in MV.org_basis: print(evec) print('Renormalized Reciprocal Vectors '+\ '$\\bfrac{'+bmhat+'^{k}}{\\abs{\\bm{'+LaTeX_base+'}_{k}}}$') ibasis = 0 while ibasis < MV.n: c_str = sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[ibasis]) print('\\bfrac{\\bm{\\hat{'+LaTeX_base+\ '}}^{'+c_str+'}}{\\abs{\\bm{'+LaTeX_base+\ '}_{'+c_str+'}}} =',MV.Rframe[ibasis]) ibasis += 1 title_str = 'Connection Multivectors: $C\\lbrc'+bstr+'\\rbrc = '+\ '\\bfrac{'+bmhat+'^{k}}{\\abs{'+bmhat+\ '_{k}}}\\pdiff{'+bstr+'}{\\theta^{k}}$' print(title_str) igrade = 1 for grade in MV.Connect[1:]: ibase = 0 for base in grade: index = MV.gabasis[igrade][ibase] sub_str = '' for i in index: sub_str += sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[i]) base_str = 'C\\lbrc\\hat{'+LaTeX_base+'}_{['+sub_str+']}\\rbrc = ' print(base_str,base) ibase += 1 igrade += 1 return @staticmethod def print_blades(): """ Set multivector output to blade representation. """ MV.bladeprint = 1 return @staticmethod def print_bases(): """ Set multivector output to base representation. """ MV.bladeprint = 0 return @staticmethod def multiplication_table(): """ Calculate geometric product base multiplication table. See reference 5 section 3 for details. """ MV.mtable = [] for igrade in MV.n1rg: MV.mtable.append([]) for ibase in range(MV.nbasis[igrade]): MV.mtable[igrade].append([]) if igrade == 0: base1 = [] else: base1 = MV.basis[igrade][ibase] for jgrade in MV.n1rg: MV.mtable[igrade][ibase].append([]) for jbase in range(MV.nbasis[jgrade]): if jgrade == 0: base2 = [] else: base2 = MV.basis[jgrade][jbase] base = base1+base2 (coefs,bases) = MV.reduce_basis(base) product = MV.convert(coefs,bases) product.name = '('+MV.basislabel[igrade][ibase]+')('+\ MV.basislabel[jgrade][jbase]+')' MV.mtable[igrade][ibase][jgrade].append(product) if MV.debug: print('Multiplication Table:') for level1 in MV.mtable: for level2 in level1: for level3 in level2: for mv in level3: mv.printmv() return @staticmethod def geometric_product(mv1,mv2): """ MV.geometric_product(mv1,mv2) calculates the geometric product the multivectors mv1 and mv2 (mv1*mv2). See reference 5 section 3. """ product = MV() if isinstance(mv1,MV) and isinstance(mv2,MV): bladeflg1 = mv1.bladeflg bladeflg2 = mv2.bladeflg if bladeflg1: mv1.convert_from_blades() if bladeflg2: mv2.convert_from_blades() for igrade in MV.n1rg: gradei = mv1.mv[igrade] if isinstance(gradei,numpy.ndarray): for ibase in range(MV.nbasis[igrade]): xi = gradei[ibase] if xi != ZERO: for jgrade in MV.n1rg: gradej = mv2.mv[jgrade] if isinstance(gradej,numpy.ndarray): for jbase in range(MV.nbasis[jgrade]): xj = gradej[jbase] if xj != ZERO: xixj = MV.mtable[igrade][ibase][jgrade][jbase].scalar_mul(xi*xj) product.add_in_place(xixj) product.bladeflg = 0 if bladeflg1: mv1.convert_to_blades() if bladeflg2: mv1.convert_to_blades() if bladeflg1 and bladeflg2: product.convert_to_blades() else: if isinstance(mv1,MV): product = mv1.scalar_mul(mv2) else: product = mv2.scalar_mul(mv1) return(product) @staticmethod def wedge(igrade1,blade1,vector2,name=''): """ Calculate the outer product of a multivector blade and a vector. See reference 5 section 5 for details. """ w12 = blade1*vector2 w21 = vector2*blade1 if igrade1%2 != 0: w = w12-w21 else: w = w12+w21 w.name = name return(w*HALF) @staticmethod def blade_table(): """ Calculate basis blades in terms of bases. See reference 5 section 5 for details. Used to convert from blade to base representation. """ MV.btable = [] MV.bladelabel = [] basis_str = MV.basislabel[1] for igrade in MV.n1rg: MV.bladelabel.append([]) if igrade == 0: MV.bladelabel[0].append('1') tmp = [MV(value=ONE,mvtype='scalar',mvname='1')] if igrade == 1: tmp = [] for ibase in range(MV.nbasis[1]): MV.bladelabel[1].append(basis_str[ibase]) tmp.append(MV(value=ibase,mvtype='basisvector',mvname=basis_str[ibase])) if igrade >= 2: tmp = [] basis = MV.basis[igrade] for blade in basis: name = '' for i in blade: name += basis_str[i]+'^' name = name[:-1] MV.bladelabel[igrade].append(name) lblade = MV.basis[igrade-1].index(blade[:-1]) rblade = blade[-1] igrade1 = igrade-1 blade1 = MV.btable[igrade1][lblade] vector2 = MV.btable[1][rblade] b1Wv2 = MV.wedge(igrade1,blade1,vector2,name) tmp.append(b1Wv2) MV.btable.append(tmp) if MV.debug: print('Blade Tabel:') for grade in MV.btable: for mv in grade: print(mv) print('Blade Labels:') print(MV.bladelabel) return @staticmethod def inverse_blade_table(): """ Calculate bases in terms of basis blades. See reference 5 section 5 for details. Used to convert from base to blade representation. """ MV.ibtable = [] for igrade in MV.n1rg: if igrade == 0: tmp = [MV(value=ONE,mvtype='scalar',mvname='1')] if igrade == 1: tmp = [] for ibase in range(MV.nbasis[1]): tmp.append(MV(value=ibase,mvtype='basisvector')) if igrade >= 2: tmp = [] iblade = 0 for blade in MV.btable[igrade]: invblade = MV() for i in range(igrade-1): invblade.mv[i] = -blade.mv[i] invblade.mv[igrade] = +blade.mv[igrade] invblade.bladeflg = 1 if igrade >= 4: jgrade = igrade-2 while jgrade > 1: for ibase in range(MV.nbasis[jgrade]): invblade.substitute_base(jgrade,ibase,MV.ibtable[jgrade][ibase]) jgrade -= 2 invblade.name = MV.basislabel[igrade][iblade] iblade += 1 tmp.append(invblade) MV.ibtable.append(tmp) if MV.debug: print('Inverse Blade Tabel:') for grade in MV.ibtable: for mv in grade: mv.printmv() return @staticmethod def outer_product(mv1,mv2): """ MV.outer_product(mv1,mv2) calculates the outer (exterior,wedge) product of the multivectors mv1 and mv2 (mv1^mv2). See reference 5 section 6. """ if type(mv1) == type(MV) and type(mv2) == type(MV): if mv1.is_scalar() and mv2.is_scalar(): return(mv1*mv2) if isinstance(mv1,MV) and isinstance(mv2,MV): product = MV() product.bladeflg = 1 mv1.convert_to_blades() mv2.convert_to_blades() for igrade1 in MV.n1rg: if isinstance(mv1.mv[igrade1],numpy.ndarray): pg1 = mv1.project(igrade1) for igrade2 in MV.n1rg: igrade = igrade1+igrade2 if igrade <= MV.n: if isinstance(mv2.mv[igrade2],numpy.ndarray): pg2 = mv2.project(igrade2) pg1pg2 = pg1*pg2 product.add_in_place(pg1pg2.project(igrade)) else: if isinstance(mv1,MV): product = mv1.scalar_mul(mv2) if isinstance(mv2,MV): product = mv2.scalar_mul(mv1) return(product) @staticmethod def inner_product(mv1,mv2,mode='s'): """ MV.inner_product(mv1,mv2) calculates the inner mode = 's' - symmetric (Doran & Lasenby) mode = 'l' - left contraction (Dorst) mode = 'r' - right contraction (Dorst) """ if type(mv1) == type(MV) and type(mv2) == type(MV): if mv1.is_scalar() and mv2.is_scalar(): return(mv1*mv2) if isinstance(mv1,MV) and isinstance(mv2,MV): product = MV() product.bladeflg = 1 mv1.convert_to_blades() mv2.convert_to_blades() for igrade1 in range(MV.n1): if isinstance(mv1.mv[igrade1],numpy.ndarray): pg1 = mv1.project(igrade1) for igrade2 in range(MV.n1): igrade = igrade1-igrade2 if mode == 's': igrade = igrade.__abs__() else: if mode == 'l': igrade = -igrade if igrade >= 0: if isinstance(mv2.mv[igrade2],numpy.ndarray): pg2 = mv2.project(igrade2) pg1pg2 = pg1*pg2 product.add_in_place(pg1pg2.project(igrade)) return(product) else: if mode == 's': if isinstance(mv1,MV): product = mv1.scalar_mul(mv2) if isinstance(mv2,MV): product = mv2.scalar_mul(mv1) else: product = None return(product) @staticmethod def addition(mv1,mv2): """ MV.addition(mv1,mv2) calculates the sum of the multivectors mv1 and mv2 (mv1+mv2). """ sum = MV() if isinstance(mv1,MV) and isinstance(mv2,MV): if mv1.bladeflg or mv2.bladeflg: mv1.convert_to_blades() mv2.convert_to_blades() sum.bladeflg = 1 for i in MV.n1rg: if isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray): sum.mv[i] = mv1.mv[i]+mv2.mv[i] else: if isinstance(mv1.mv[i],numpy.ndarray) and not isinstance(mv2.mv[i],numpy.ndarray): sum.mv[i] = +mv1.mv[i] else: if isinstance(mv2.mv[i],numpy.ndarray) and not isinstance(mv1.mv[i],numpy.ndarray): sum.mv[i] = +mv2.mv[i] return(sum) else: if isinstance(mv1,MV): return(mv1+MV(mv2,'scalar')) else: return(MV(mv1,'scalar')+mv2) @staticmethod def subtraction(mv1,mv2): """ MV.subtraction(mv1,mv2) calculates the difference of the multivectors mv1 and mv2 (mv1-mv2). """ diff = MV() if isinstance(mv1,MV) and isinstance(mv2,MV): if mv1.bladeflg or mv2.bladeflg: mv1.convert_to_blades() mv2.convert_to_blades() diff.bladeflg = 1 for i in MV.n1rg: if isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray): diff.mv[i] = mv1.mv[i]-mv2.mv[i] else: if isinstance(mv1.mv[i],numpy.ndarray) and not isinstance(mv2.mv[i],numpy.ndarray): diff.mv[i] = +mv1.mv[i] else: if not isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray): diff.mv[i] = -mv2.mv[i] return(diff) else: if isinstance(mv1,MV): return(mv1-MV(mv2,'scalar')) else: return(MV(mv1,'scalar')-mv2) @staticmethod def vdiff(vec,x): dvec = numpy.array(len(vec)*[ZERO]) ivec = 0 for veci in vec: dvec[ivec] = sympy.diff(veci,x) ivec += 1 return(dvec) @staticmethod def scalar_to_symbol(scalar): if isinstance(scalar,MV): return(scalar.mv[0][0]) if type(scalar) == list: sym = [] for x in scalar: sym.append(MV.scalar_to_symbol(x)) return(sym) return(scalar) def __init__(self,value='',mvtype='',mvname='',fct=False,vars=None): """ Initialization of multivector X. Inputs are as follows mvtype value result default default Zero multivector 'basisvector' int i ith basis vector 'basisbivector' int i ith basis bivector 'scalar' symbol x scalar of value x string s 'grade' [int i, symbol array A] X.grade(i) = A [int i, string s] 'vector' symbol array A X.grade(1) = A string s 'grade2' symbol array A X.grade(2) = A string s 'pseudo' symbol x X.grade(n) = x string s 'spinor' string s spinor with coefficients s__indices and name sbm mvname is name of multivector. If fct is 'True' and MV.coords is defined in MV.setup then a multivector field of MV.coords is instantiated. """ self.name = mvname self.mv = MV.n1*[0] self.bladeflg = 0 #1 for blade expansion self.puregrade = 1 if mvtype == 'basisvector': self.mv[1] = numpy.array(MV.nbasis[1]*[ZERO],dtype=numpy.object) self.mv[1][value] = ONE if mvtype == 'basisbivector': self.mv[2] = numpy.array(MV.nbasis[2]*[ZERO],dtype=numpy.object) self.mv[2][value] = ONE if mvtype == 'scalar': if isinstance(value,str): value = sympy.Symbol(value) if isinstance(value,int): value = sympy.Rational(value) self.mv[0] = numpy.array([value],dtype=numpy.object) if mvtype == 'pseudo': if isinstance(value,str): value = sympy.Symbol(value) self.mv[MV.n] = numpy.array([value],dtype=numpy.object) if mvtype == 'vector': if isinstance(value,str): #Most general vector symbol_str = '' for ibase in MV.nrg: if MV.coords == None: symbol = value+'__'+str(ibase+MV.index_offset) symbol_str += symbol+' ' else: symbol = value+'__'+(MV.coords[ibase]).name symbol_str += symbol+' ' symbol_lst = sympy.symbols(symbol_str) self.mv[1] = numpy.array(symbol_lst,dtype=numpy.object) self.name = value else: value = MV.pad_zeros(value,MV.nbasis[1]) self.mv[1] = numpy.array(value,dtype=numpy.object) if mvtype == 'grade2': if isinstance(value,str): #Most general grade-2 multivector if value != '': symbol_str = '' for base in MV.basis[2]: symbol = value+MV.construct_index(base) symbol_str += symbol+' ' symbol_lst = sympy.symbols(symbol_str) self.mv[2] = numpy.array(symbol_lst,dtype=numpy.object) else: value = MV.pad_zeros(value,MV.nbasis[2]) self.mv[2] = numpy.array(value,dtype=numpy.object) if mvtype == 'grade': igrade = value[0] coefs = value[1] if isinstance(coefs,str): #Most general pure grade multivector base_symbol = coefs coefs = [] bases = MV.basis[igrade] if igrade == 0: self.mv[0] = numpy.array([sympy.Symbol(base_symbol)],dtype=numpy.object) else: for base in bases: coef = base_symbol+MV.construct_index(base) coef = sympy.Symbol(coef) coefs.append(coef) self.mv[igrade] = numpy.array(coefs,dtype=numpy.object) else: self.mv[igrade] = coefs if mvtype == 'base': self.mv[value[0]] = numpy.array(MV.nbasis[value[0]]*[ZERO],dtype=numpy.object) self.mv[value[0]][value[1]] = ONE if mvtype == 'spinor': if isinstance(value,str): #Most general spinor for grade in MV.n1rg: if grade%2 == 0: symbol_str = '' if grade != 0: symbol_lst = [] for base in MV.basis[grade]: symbol = sympy.Symbol(value+MV.construct_index(base)) symbol_lst.append(symbol) self.mv[grade] = numpy.array(symbol_lst,dtype=numpy.object) else: self.mv[0] = numpy.array([sympy.Symbol(value)],dtype=numpy.object) self.name = value+'bm' if isinstance(value,str) and mvtype == '': #Most general multivector if value != '': for grade in MV.n1rg: symbol_str = '' if grade != 0: symbol_lst = [] for base in MV.basis[grade]: symbol = sympy.Symbol(value+MV.construct_index(base)) symbol_lst.append(symbol) self.mv[grade] = numpy.array(symbol_lst,dtype=numpy.object) else: self.mv[0] = numpy.array([sympy.Symbol(value)],dtype=numpy.object) self.name = value+'bm' if fct: if vars is not None: vars = tuple(vars) for grade in MV.n1rg: if not isinstance(self.mv[grade],int): if grade == 0: coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(self.mv[0][0]) if vars == None and MV.coords is not None: self.mv[0]= numpy.array([sympy.Function(coef)(*MV.coords)],dtype=numpy.object) else: self.mv[0]= numpy.array([sympy.Function(coef)(*vars)],dtype=numpy.object) else: for base in range(MV.nbasis[grade]): coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(self.mv[grade][base]) if vars == None and MV.coords is not None: self.mv[grade][base] = sympy.Function(coef)(*MV.coords) else: self.mv[grade][base] = sympy.Function(coef)(*vars) @staticmethod def construct_index(base): index_str = '' if len(base) == 0: return('') if MV.coords == None: for ix in base: index_str += str(ix+MV.index_offset) else: for ix in base: index_str += (MV.coords[ix]).name return('__'+index_str) def set_name(self,namestr): self.name = namestr return def max_grade(self): """ X.max_grade() is maximum grade of non-zero grades of X. """ for i in range(MV.n,-1,-1): if isinstance(self.mv[i],numpy.ndarray): return(i) return(-1) @staticmethod def coord(xname,offset=0): xi_str = '' for i in MV.nrg: xi_str += xname+str(i+offset)+' ' xi = sympy.symbols(xi_str) x = MV(xi,'vector') return(x) def x(self,i): if isint(self.mv[1]): return(ZERO) return(self.mv[1][i]) def set_coef(self,grade,base,value): if isinstance(self.mv[grade],int): self.mv[grade] = numpy.array(MV.nbasis[grade]*[ZERO],dtype=numpy.object) self.mv[grade][base] = value return @staticmethod def named(mvname,value='',mvtype=''): name = mvname tmp = MV(value=value,mvtype=mvtype,mvname=name) setattr(sys.modules[__name__],name,tmp) return @staticmethod def printnm(tpl): for a in tpl: print(a.name,' =',a.mv) return def __str__(self): return(MV.str_rep(self)) def printmv(self,name=''): title = '' if name: title += name+' = ' else: if self.name: title += self.name+' = ' print(title+MV.str_rep(self)) return def set_value(self,igrade,ibase,value): if isinstance(self.mv[igrade],numpy.ndarray): self.mv[igrade][ibase] = value else: self.mv[igrade] = numpy.array(MV.nbasis[igrade]*[ZERO],dtype=numpy.object) self.mv[igrade][ibase] = value return def add_in_place(self,mv): """ X.add_in_place(mv) increments multivector X by multivector mv. """ if self.bladeflg or mv.bladeflg: self.convert_to_blades() mv.convert_to_blades() for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray): self.mv[i] += mv.mv[i] else: if not isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray): self.mv[i] = +mv.mv[i] return def sub_in_place(self,mv): """ X.sub_in_place(mv) decrements multivector X by multivector mv. """ if self.bladeflg or mv.bladeflg: self.convert_to_blades() mv.convert_to_blades() for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray): self.mv[i] -= mv.mv[i] else: if not isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray): self.mv[i] = +mv.mv[i] return def __pos__(self): p = MV() p.bladeflg = self.bladeflg for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): p.mv[i] = +self.mv[i] return(p) def __neg__(self): n = MV() n.bladeflg = self.bladeflg for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): n.mv[i] = -self.mv[i] return(n) def __add_ab__(self,mv): self.add_in_place(mv) return def __sub_ab__(self,mv): self.sub_in_place(mv) return def __add__(self,mv): """See MV.addition(self,mv)""" return(MV.addition(self,mv)) def __radd__(self,mv): """See MV.addition(mv,self)""" return(MV.addition(mv,self)) def __sub__(self,mv): """See MV.subtraction(self,mv)""" return(MV.subtraction(self,mv)) def __rsub__(self,mv): """See MV.subtraction(mv,self)""" return(MV.subtraction(mv,self)) def __xor__(self,mv): """See MV.outer_product(self,mv)""" return(MV.outer_product(self,mv)) def __pow__(self,mv): """See MV.outer_product(self,mv)""" return(MV.outer_product(self,mv)) def __rxor__(self,mv): """See MV.outer_product(mv,self)""" return(MV.outer_product(mv,self)) def __or__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(self,mv)) def __ror__(self,mv): """See MV.inner_product(mv,self)""" return(MV.inner_product(mv,self)) def __lt__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(self,mv,'l')) def __lshift__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(self,mv,'l')) def __rlshift__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(mv,self,'l')) def lc(self,mv): return(MV.inner_product(self,mv,'l')) def __gt__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(self,mv,'r')) def __rshift__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(self,mv,'r')) def __rrshift__(self,mv): """See MV.inner_product(self,mv)""" return(MV.inner_product(mv,self,'r')) def rc(self,mv): return(MV.inner_product(self,mv,'r')) def scalar_mul(self,c): """ Y = X.scalar_mul(c), multiply multivector X by scalar c and return result. """ mv = MV() mv.bladeflg = self.bladeflg mv.puregrade = self.puregrade for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): #print self.mv[i] #print c,type(c) mv.mv[i] = self.mv[i]*c return(mv) def scalar_mul_inplace(self,c): """ X.scalar_mul_inplace(c), multiply multivector X by scalar c and save result in X. """ for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): self.mv[i] = self.mv[i]*c return def __mul__(self,mv): """See MV.geometric_product(self,mv)""" return(MV.geometric_product(self,mv)) def __rmul__(self,mv): """See MV.geometric_product(mv,self)""" return(MV.geometric_product(mv,self)) def __div__(self,scalar): div = MV() div.bladeflg = self.bladeflg for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): div.mv[i] = self.mv[i]/scalar return(div) def __div_ab__(self,scalar): for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): self.mv[i] /= scalar return def __call__(self,igrade=0,ibase=0): """ X(i,j) returns symbol in ith grade and jth base or blade of multivector X. """ if not isinstance(self.mv[igrade],numpy.ndarray): return(ZERO) return(self.mv[igrade][ibase]) @staticmethod def equal(mv1,mv2): mv1.compact() if isinstance(mv2,MV): mv2.compact() pure_grade = mv1.is_pure() if not isinstance(mv2,MV) and pure_grade != 0: return(False) if not isinstance(mv2,MV) and pure_grade == 0: if isinstance(mv1.mv[0],int): return(mv2 == 0) else: return(mv1.mv[0][0] == mv2) for (mvi,mvj) in zip(mv1.mv,mv2.mv): if isint(mvi) ^ isint(mvj): return(False) if isinstance(mvi,numpy.ndarray) and isinstance(mvj,numpy.ndarray): for (x,y) in zip(mvi,mvj): if x != y: return(False) return(True) def __eq__(self,mv): return(MV.equal(self,mv)) def copy(self,sub=0): """ Y = X.copy(), make a deep copy of multivector X in multivector Y so that Y can be modified without affecting X. """ cpy = MV() cpy.name = self.name cpy.bladeflg = self.bladeflg cpy.puregrade = self.puregrade for i in MV.n1rg: if sub: if isinstance(self.mv[i],numpy.ndarray): cpy.mv[i] = -self.mv[i] else: if isinstance(self.mv[i],numpy.ndarray): cpy.mv[i] = +self.mv[i] return(cpy) def substitute_base(self,igrade,base,mv): if not isinstance(self.mv[igrade],numpy.ndarray): return if isinstance(base,numpy.ndarray): ibase = MV.basis[igrade].index(base) else: ibase = base coef = self.mv[igrade][ibase] if coef == ZERO: return self.mv[igrade][ibase] = ZERO self.add_in_place(mv*coef) return def convert_to_blades(self): """ X.convert_to_blades(), inplace convert base representation to blade representation. See reference 5 section 5. """ if self.bladeflg: return self.bladeflg = 1 for igrade in range(2,MV.n1): if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): coef = self.mv[igrade][ibase] if (not coef == ZERO): self.mv[igrade][ibase] = ZERO self.add_in_place(MV.ibtable[igrade][ibase]*coef) return def convert_from_blades(self): """ X.convert_from_blades(), inplace convert blade representation to base representation. See reference 5 section 5. """ if not self.bladeflg: return self.bladeflg = 0 for igrade in range(2,MV.n1): if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): coef = self.mv[igrade][ibase] if (not coef == ZERO): self.mv[igrade][ibase] = ZERO self.add_in_place(MV.btable[igrade][ibase]*coef) return def project(self,r): """ Grade projection operator. For multivector X, X.project(r) returns multivector of grade r components of X if r is an integer. If r is a multivector X.project(r) returns a multivector consisting of the grade of X for which r has non- zero grades. For example if X is a general multivector and r is a general spinor then X.project(r) will return the even grades of X. """ if isinstance(r,int): grade_r = MV() if r > MV.n: return(grade_r) self.convert_to_blades() if not isinstance(self.mv[r],numpy.ndarray): return(grade_r) grade_r.bladeflg = 1 grade_r.puregrade = 1 grade_r.mv[r] = +self.mv[r] return(grade_r) if isinstance(r,MV): self.convert_to_blades() r.convert_to_blades() proj = MV() for i in MV.n1rg: if not isinstance(r.mv[i],int): proj.mv[i] = self.mv[i] proj.bladeflg = self.bladeflg return(proj) return(None) def even(self): """ Even grade projection operator. For multivector X, X.even() returns multivector of even grade components of X. """ egrades = MV() self.convert_to_blades() egrades.bladeflg = self.bladeflg egrades.puregrade = self.puregrade for igrade in range(0,MV.n1,2): egrades.mv[igrade] = +self.mv[igrade] return(egrades) def odd(self): """ Odd grade projection operator. For multivector X, X.odd() returns multivector of odd grade components of X. """ ogrades = MV() self.convert_to_blades() ogrades.bladeflg = self.bladeflg ogrades.puregrade = self.puregrade for igrade in range(1,MV.n1,2): ogrades.mv[igrade] = +self.mv[igrade] return(ogrades) def rev(self): """ Revision operator. For multivector X, X.rev() returns reversed multivector of X. """ revmv = MV() self.convert_to_blades() revmv.bladeflg = self.bladeflg revmv.puregrade = self.puregrade for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): if igrade < 2 or (not (((igrade*(igrade-1))/2)%2)): revmv.mv[igrade] = +self.mv[igrade] else: revmv.mv[igrade] = -self.mv[igrade] return(revmv) def cse(self,grade): cse_lst = [] if isinstance(self.mv[grade],numpy.ndarray): for ibase in range(MV.nbasis[grade]): if self.mv[grade][ibase] != ZERO: cse_lst.append(sympy.cse(self.mv[grade][ibase])) return(cse_lst) def div(self,grade,divisor): div_lst = [] if isinstance(self.mv[grade],numpy.ndarray): for ibase in range(MV.nbasis[grade]): if self.mv[grade][ibase] != ZERO: div_lst.append(self.mv[grade][ibase].as_coefficient(divisor)) return(div_lst) # don't know which one is correctly named def div(self): return self.grad_int() def collect(self,lst): """ Applies sympy collect function to each component of multivector. """ for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): if self.mv[igrade][ibase] != ZERO: self.mv[igrade][ibase] = \ sympy.collect(self.mv[igrade][ibase],lst) return def sqrfree(self,lst): """ Applies sympy sqrfree function to each component of multivector. """ for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): if self.mv[igrade][ibase] != ZERO: self.mv[igrade][ibase] = \ sympy.sqrfree(self.mv[igrade][ibase],lst) return def flatten(self): flst = [] for igrade in MV.n1rg: if isinstance(self.mv[igrade],int): flst += MV.nbasis[igrade]*[ZERO] else: for coef in self.mv[igrade]: flst.append(coef) return(flst) def subs(self,*args): X = MV() X.bladeflg = self.bladeflg X.puregrade = self.puregrade for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): X.mv[igrade] = numpy.array(substitute_array(self.mv[igrade],*args)) return(X) def sub_mv(self,mv1,mv2): mv1_flat = mv1.flatten() mv2_flat = mv2.flatten() self.sub_scalar(mv1_flat,mv2_flat) return def sub_scalar(self,expr1,expr2): if (isinstance(expr1,list) and isinstance(expr2,list)) or \ (isinstance(expr1,tuple) and isinstance(expr2,tuple)): for (var1,var2) in zip(expr1,expr2): self.sub_scalar(var1,var2) else: for igrade in MV.n1rg: if not isinstance(self.mv[igrade],int): for ibase in range(MV.nbasis[igrade]): if expr1 != ZERO: self.mv[igrade][ibase] = self.mv[igrade][ibase].subs(expr1,expr2) return def simplify(self): """ Applies sympy simplify function to each component of multivector. """ for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): if self.mv[igrade][ibase] != ZERO: self.mv[igrade][ibase] = \ sympy.simplify(self.mv[igrade][ibase]) return def trigsimp(self): """ Applies sympy trigsimp function to each component of multivector using the recursive option. """ for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): if self.mv[igrade][ibase] != ZERO: self.mv[igrade][ibase] = \ sympy.trigsimp(self.mv[igrade][ibase],deep=True, recursive=True) return def cancel(self): """ Applies sympy cancel function to each component of multivector. """ for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): if self.mv[igrade][ibase] != ZERO: self.mv[igrade][ibase] = \ sympy.cancel(self.mv[igrade][ibase]) return def expand(self): """ Applies sympy expand function to each component of multivector. """ for igrade in MV.n1rg: if isinstance(self.mv[igrade],numpy.ndarray): for ibase in range(MV.nbasis[igrade]): if self.mv[igrade][ibase] != ZERO: self.mv[igrade][ibase] = \ self.mv[igrade][ibase].expand() return def is_pure(self): igrade = -1 ngrade = 0 self.compact() for i in MV.n1rg: if isinstance(self.mv[i],numpy.ndarray): for base in self.mv[i]: if base != 0: igrade = i ngrade += 1 break if ngrade > 1: return(-1) if igrade == -1: return(0) return(igrade) def is_scalar(self): if self.is_pure() == 0: return(True) return(False) def compact(self): """ Convert zero numpy arrays to single integer zero place holder in grade list for instantiated multivector. For example if numpy array of grade one components is a zero array then replace with single integer equal to zero. """ for i in MV.n1rg: if not isint(self.mv[i]): zero_flg = True for x in self.mv[i]: if x != 0: zero_flg = False break if zero_flg: self.mv[i] = 0 return def diff(self,x): """ Calculate partial derivative of multivector with respect to argument x. """ D = MV() igrade = 0 for grade in self.mv: if not isinstance(grade,int): D.mv[igrade] = MV.vdiff(grade,x) igrade += 1 return(D) def ddt(self): if MV.coords[0] != sympy.Symbol('t'): return(MV()) dxdt = self.diff(MV.coords[0]) for ibase in MV.nrg: dxdt += self.mv[1][ibase]*MV.dedt[ibase] #dxdt.simplify() #dxdt.trigsimp() return(dxdt) def grad(self): """ Calculate geometric (grad) derivative of multivector function """ D = [] dD = MV() for theta in MV.coords: D.append(self.diff(theta)) if MV.curvilinear_flg: recp = MV.Rframe else: recp = MV.brecp for (rbase,iD) in zip(recp,D): dD.add_in_place(rbase*iD) if MV.curvilinear_flg: #Add Connection igrade = 1 while igrade <= MV.n: coefs = self.mv[igrade] if type(coefs) != int: for (coef,connect) in zip(coefs,MV.Connect[igrade]): dD.add_in_place(coef*connect) igrade += 1 return(dD) def grad_ext(self): """ Calculate outer (exterior,curl) derivative of multivector function. """ D = [] dD = MV() for ix in MV.coords: D.append(self.diff(ix)) if MV.curvilinear_flg: recp = MV.Rframe else: recp = MV.brecp for (irbase,iD) in zip(recp,D): dD.add_in_place(irbase^iD) if MV.curvilinear_flg: #Add Connection igrade = 1 while igrade <= MV.n: coefs = self.mv[igrade] if type(coefs) != int: for (coef,connect) in zip(coefs,MV.Connect[igrade]): if igrade < MV.n: dD.add_in_place(coef*connect.project(igrade+1)) igrade += 1 return(dD) def curl(self): return(self.grad_ext()) def grad_int(self): """ Calculate inner (interior,div) derivative of multivector function. """ D = [] dD = MV() for ix in MV.coords: D.append(self.diff(ix)) if MV.curvilinear_flg: recp = MV.Rframe else: recp = MV.brecp for (irbase,iD) in zip(recp,D): dD.add_in_place(irbase|iD) if MV.curvilinear_flg: #Add Connection igrade = 1 while igrade <= MV.n: coefs = self.mv[igrade] if type(coefs) != int: for (coef,connect) in zip(coefs,MV.Connect[igrade]): dD.add_in_place(coef*connect.project(igrade-1)) igrade += 1 return(dD) def mag2(self): """ Calculate scalar component of square of multivector. """ return((self|self)()) def Name(self): """ Get LaTeX name of multivector. """ return(sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(self.name))
[docs]def set_names(var_lst,var_str): """ Set the names of a list of multivectors (var_lst) for a space delimited string (var_str) containing the names. """ var_str_lst = var_str.split() if len(var_lst) == len(var_str_lst): for (var,var_name) in zip(var_lst,var_str_lst): var.set_name(var_name) return sys.stderr.write('Error in set_names. Lists incongruent!\n') sys.exit() return
[docs]def reciprocal_frame(vlst,names=''): """ Calculate reciprocal frame of list (vlst) of vectors. If desired name each vector in list of reciprocal vectors with names in space delimited string (names). """ E = vlst[0] recp = [] if type(names) != str: name_lst = names else: if names != '': name_lst = names.split() for i in range(1,MV.n): E = E^vlst[i] for i in range(MV.n): tmp = ONE if i%2 != 0: tmp = -ONE for j in range(MV.n): if i != j: tmp = tmp^vlst[j] tmp = tmp*E recp.append(tmp) Esq = sympy.trigsimp(E.mag2(),deep=True,recursive=True) print(Esq) print(sympy.simplify(Esq)) Esq_inv = ONE/Esq i = 0 for i in range(MV.n): recp[i].trigsimp() recp[i] = recp[i]*Esq_inv if names != '': recp[i].set_name(name_lst[i]) i += 1 return(recp)
def S(value): return(MV(value,'scalar'))