sympy/galgebra/GAsympy.py

   1 #!/usr/bin/python
   2
   3 """
   4 The module symbolicGA implements symbolic Geometric Algebra in python.
   5 The relevant references for this module are:
   6
   7     1. "Geometric Algebra for Physicists" by C. Doran and A. Lazenby,
   8        Cambridge University Press, 2003.
   9
  10     2. GiNaC 1.4.1, An open framwork for symbolic computation with the
  11        C++ programming language, 9/5/2007, http://www.ginac.de
  12
  13     3. Symbolic Computation with GiNaC and Python 1.3.5, 2007-02-01,
  14        http://swiginac.berlios.de/ginac-tutorial.py.html
  15
  16     4. Sympy Tutorial, http://code.google.com/p/sympy/wiki/Tutorial
  17
  18     5. "Design of a Python Module for Symbolic Geometric Algebra
  19        Calculations" by Alan Bromborsky, included as symbolGA.pdf
  20 """
  21
  22 import numpy
  23 import os,sys,string,types,re,copy
  24 import sympy
  25 from latex_out import *
  26
  27 NUMPAT = re.compile( '([\-0-9])|([\-0-9]/[0-9])')
  28 """Re pattern for rational number"""
  29
  30 ZERO = sympy.Rational(0)
  31 ONE  = sympy.Rational(1)
  32 TWO  = sympy.Rational(2)
  33 HALF = sympy.Rational(1,2)
  34
  35 global MAIN_PROGRAM
  36
  37 MAIN_PROGRAM = ''
  38
  39 def set_main(main_program):
  40     global MAIN_PROGRAM
  41     MAIN_PROGRAM = main_program
  42     return
  43
  44 def plist(lst):
  45     if type(lst) == types.ListType:
  46         for x in lst:
  47             plist(x)
  48     else:
  49         print lst
  50     return
  51
  52 def numeric(num_str):
  53     """
  54     Returns rational numbers compatible with symbols.
  55     Input is a string representing a fraction or integer
  56     or a simple integer.
  57     """
  58     if type(num_str) == types.IntType:
  59         a = num_str
  60         b = 1
  61     else:
  62         tmp = num_str.split('/')
  63         if len(tmp) == 1:
  64             a = int(tmp[0])
  65             b = 1
  66         else:
  67             a = int(tmp[0])
  68             b = int(tmp[1])
  69     return(sympy.Rational(a,b))
  70
  71
  72 def symbol(sym_str):
  73     """
  74     Symbol converts a string to a sympy/sympy symbol.
  75     """
  76     sym = sympy.Symbol(sym_str)
  77     return(sym)
  78
  79 def expand(expr):
  80     return(sympy.expand(expr))
  81
  82
  83 def collect(expr,lst):
  84     """
  85     Wrapper for sympy.collect and sympy.collect.
  86     See references 2, 3, and 4.
  87     """
  88     lst = MV.scalar_to_symbol(lst)
  89     return(sympy.collect(expr,lst))
  90
  91 def sqrfree(expr,lst):
  92     """
  93     Wrapper for sympy.sqrfree and sympy.factor.
  94     See references 2, 3, and 4.
  95     """
  96     return(sympy.sqrfree(expr,lst))
  97
  98 def collect_common_factors(expr):
  99     """
 100     Wrapper for sympy.collect_common_factors and sympy.factor.
 101     See references 2, 3, and 4.
 102     """
 103     return(sympy.collect_common_factors(expr))
 104
 105 def sqrt(expr):
 106     return(sympy.sqrt(expr))
 107
 108 def isint(a):
 109     """
 110     Test for integer.
 111     """
 112     return(type(a) == types.IntType)
 113
 114 def make_null_array(n):
 115     """
 116     Return list of n empty lists.
 117     """
 118     a = []
 119     for i in range(n):
 120         a.append([])
 121     return(a)
 122
 123 def test_int_flgs(lst):
 124     """
 125     Test if all elements in list are 0.
 126     """
 127     for i in lst:
 128         if i:
 129             return(1)
 130     return(0)
 131
 132 def comb(N,P):
 133     """
 134     Calculates the combinations of the integers [0,N-1] taken P at a time.
 135     The output is a list of lists of integers where the inner lists are
 136     the different combinations. Each combination is sorted in ascending
 137     order.
 138     """
 139     def rloop(n,p,combs,comb):
 140         if p:
 141             for i in range(n):
 142                 newcomb = comb+[i]
 143                 np = p-1
 144                 rloop(i,np,combs,newcomb)
 145         else:
 146             combs.append(comb)
 147     combs = []
 148     rloop(N,P,combs,[])
 149     for comb in combs:
 150         comb.sort()
 151     return(combs)
 152
 153 def diagpq(p,q=0):
 154     """
 155     Return string equivalent metric tensor for signature (p,q).
 156     """
 157     n = p+q
 158     D = ''
 159     rn = range(n)
 160     for i in rn:
 161         for j in rn:
 162             if i ==j:
 163                 if i < p:
 164                     D += '1 '
 165                 else:
 166                     D += '-1 '
 167             else:
 168                 D += '0 '
 169         D = D[:-1]+','
 170     return(D)
 171
 172 def make_scalars(symnamelst):
 173     """
 174     make_symbols takes a string of symbol names separated by
 175     blanks and converts them to MV scalars separately
 176     accessible by the main program and addition returns a list
 177     of the symbols.
 178     """
 179     global MAIN_PROGRAM
 180     if type(symnamelst) == types.StringType:
 181         symnamelst = symnamelst.split()
 182     scalar_lst = []
 183     isym = 0
 184     for name in symnamelst:
 185         tmp = symbol(name)
 186         tmp = MV(tmp,'scalar')
 187         scalar_lst.append(tmp)
 188         setattr(MAIN_PROGRAM,name,tmp)
 189         isym += 1
 190     return(scalar_lst)
 191
 192 def make_symbols(symnamelst):
 193     """
 194     make_symbols takes a string of symbol names separated by
 195     blanks and converts them to MV scalars separately
 196     accessible by the main program and addition returns a list
 197     of the symbols.
 198     """
 199     global MAIN_PROGRAM
 200     if type(symnamelst) == types.StringType:
 201         symnamelst = symnamelst.split()
 202     sym_lst = []
 203     isym = 0
 204     for name in symnamelst:
 205         tmp = symbol(name)
 206         sym_lst.append(tmp)
 207         setattr(MAIN_PROGRAM,name,tmp)
 208         isym += 1
 209     return(sym_lst)
 210
 211 def israt(numstr):
 212     """
 213     Test if string represents a rational number.
 214     """
 215     global NUMPAT
 216     if NUMPAT.match(numstr):
 217         return(1)
 218     return(0)
 219
 220 def dualsort(lst1, lst2):
 221     """
 222     Inplace dual sort of lst1 and lst2 keyed on sorted lst1.
 223     """
 224     _indices = range(len(lst1))
 225     _indices.sort(key=lst2.__getitem__)
 226     lst1[:] = map(lst1.__getitem__, _indices)
 227     lst2[:] = map(lst2.__getitem__, _indices)
 228     return
 229
 230 def cp(A,B):
 231     """
 232     Calculates the comutator product (A*B-B*A)/2 for
 233     the objects A and B.
 234     """
 235     return(HALF*(A*B-B*A))
 236
 237 class MV:
 238     def pad_zeros(value,n):
 239         """
 240         Pad list with zeros to lenght n. If length is > n
 241         truncate list.  Return padded list.
 242         """
 243         nvalue = len(value)
 244         if nvalue < n:
 245             value = value+(n-nvalue)*[ZERO]
 246         if nvalue > n:
 247             value = value[:n]
 248         return(value)
 249     pad_zeros = staticmethod(pad_zeros)
 250
 251     def define_basis(basis):
 252         """
 253         Calculates all the MV static variables needed for
 254         basis operations.  See reference 5 section 2.
 255         """
 256         MV.vbasis     = basis
 257         MV.vsyms      = make_symbols(MV.vbasis)
 258         MV.n          = len(MV.vbasis)
 259         MV.nrg        = range(MV.n)
 260         MV.n1         = MV.n+1
 261         MV.n1rg       = range(MV.n1)
 262         MV.npow       = 2**MV.n
 263         MV.index      = range(MV.n)
 264         MV.gabasis    = [[]]
 265         MV.basis      = (MV.n+1)*[0]
 266         MV.basislabel = (MV.n+1)*[0]
 267         MV.basis[0]   = []
 268         MV.basislabel[0] = '1'
 269         MV.nbasis = numpy.array((MV.n+1)*[1],dtype=numpy.object)
 270         for igrade in range(1,MV.n+1):
 271             tmp = comb(MV.n,igrade)
 272             MV.gabasis += [tmp]
 273             ntmp = len(tmp)
 274             MV.nbasis[igrade] = ntmp
 275             MV.basis[igrade] = tmp
 276             gradelabels = []
 277             for i in range(ntmp):
 278                 bstr = ''
 279                 for j in tmp[i]:
 280                     bstr += MV.vbasis[j]
 281                 gradelabels.append(bstr)
 282             MV.basislabel[igrade] = gradelabels
 283         if MV.debug:
 284             print 'basis strings =',MV.vbasis
 285             print 'basis symbols =',MV.vsyms
 286             print 'basis labels  =',MV.basislabel
 287             print 'basis         =',MV.basis
 288             print 'grades        =',MV.nbasis
 289             print 'index         =',MV.index
 290         return
 291     define_basis = staticmethod(define_basis)
 292
 293     def define_metric(metric):
 294         """
 295         Calculates all the MV static variables needed for
 296         metric operations.  See reference 5 section 2.
 297         """
 298         name_flg = False
 299         MV.g = []
 300         MV.metric = numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object)
 301         if not metric:
 302             metric = numpy.array(MV.n*[MV.n*['#']],dtype=numpy.object)
 303         for i in MV.index:
 304             for j in MV.index:
 305                 gij = metric[i][j]
 306                 if israt(gij):
 307                     MV.metric[i][j] = numeric(gij)
 308                 else:
 309                     if gij == '#':
 310                         name_flg = True
 311                         if i == j:
 312                             gij = '('+MV.vbasis[j]+'**2)'
 313                             name = MV.vbasis[j]+'sq'
 314                         else:
 315                             gij = '('+MV.vbasis[min(i,j)]+'.'+MV.vbasis[max(i,j)]+')'
 316                             name = MV.vbasis[min(i,j)]+'dot'+MV.vbasis[max(i,j)]
 317                     tmp = symbol(gij)
 318                     MV.metric[i][j] = tmp
 319                     if i <= j:
 320                         MV.g.append(tmp)
 321                         if name_flg:
 322                             setattr(MAIN_PROGRAM,name,tmp)
 323                             name_flg = False
 324         if MV.debug:
 325             print 'metric =',MV.metric
 326         return
 327     define_metric = staticmethod(define_metric)
 328
 329     def define_reciprocal_frame():
 330         """
 331         Calculates unscaled reciprocal vectors (MV.brecp) and scale
 332         factor (MV.Esq). The ith scaled reciprocal vector is
 333         (1/MV.Esq)*MV.brecp[i].  The pseudoscalar for the set of
 334         basis vectors is MV.E.
 335         """
 336         if MV.tables_flg:
 337             if MV.rframe_flg:
 338                 return
 339             MV.rframe_flg = True
 340             MV.E = MV.bvec[0]
 341             MV.brecp = []
 342             for i in range(1,MV.n):
 343                 MV.E = MV.E^MV.bvec[i]
 344             for i in range(MV.n):
 345                 tmp = ONE
 346                 if i%2 != 0:
 347                     tmp = -ONE
 348                 for j in range(MV.n):
 349                     if i != j:
 350                         tmp = tmp^MV.bvec[j]
 351                 tmp = tmp*MV.E
 352                 MV.brecp.append(tmp)
 353             MV.Esq = (MV.E*MV.E)()
 354             MV.Esq_inv = ONE/MV.Esq
 355             for i in range(MV.n):
 356                 MV.brecp[i] = MV.brecp[i]*MV.Esq_inv
 357         return
 358     define_reciprocal_frame = staticmethod(define_reciprocal_frame)
 359
 360     def reduce_basis_loop(blst):
 361         """
 362         Makes one pass through basis product representation for
 363         reduction of representation to normal form.
 364         See reference 5 section 3.
 365         """
 366         nblst = len(blst)
 367         if nblst <= 1:
 368             return(1)
 369         jstep = 1
 370         while jstep <nblst:
 371             istep = jstep-1
 372             if blst[istep] == blst[jstep]:
 373                 i = blst[istep]
 374                 if len(blst) >2:
 375                     blst = blst[:istep]+blst[jstep+1:]
 376                 else:
 377                     blst = []
 378                 if len(blst) <= 1 or jstep == nblst-1:
 379                     blst_flg = 0
 380                 else:
 381                     blst_flg = 1
 382                 return(MV.metric[i][i],blst,blst_flg)
 383             if blst[istep] > blst[jstep]:
 384                 blst1 = blst[:istep]+blst[jstep+1:]
 385                 a1 = TWO*MV.metric[blst[jstep]][blst[istep]]
 386                 blst = blst[:istep]+[blst[jstep]]+[blst[istep]]+blst[jstep+1:]
 387                 if len(blst1) <= 1:
 388                     blst1_flg = 0
 389                 else:
 390                     blst1_flg = 1
 391                 return(a1,blst1,blst1_flg,blst)
 392             jstep +=1
 393         return(1)
 394     reduce_basis_loop = staticmethod(reduce_basis_loop)
 395
 396     def reduce_basis(blst):
 397         """
 398         Repetively applies reduce_basis_loop to basis
 399         product representation untill normal form is
 400         realized.  See reference 5 section 3.
 401         """
 402         if blst == []:
 403             blst_coef   = []
 404             blst_expand = []
 405             for i in MV.n1rg:
 406                 blst_coef.append([])
 407                 blst_expand.append([])
 408             blst_expand[0].append([])
 409             blst_coef[0].append(ONE)
 410             return(blst_coef,blst_expand)
 411         blst_expand = [blst]
 412         blst_coef      = [ONE]
 413         blst_flg         = [1]
 414         while test_int_flgs(blst_flg):
 415             for i in range(len(blst_flg)):
 416                 if blst_flg[i]:
 417                     tmp = MV.reduce_basis_loop(blst_expand[i])
 418                     if tmp ==1:
 419                         blst_flg[i] = 0
 420                     else:
 421                         if len(tmp) == 3:
 422                             blst_coef[i] = tmp[0]*blst_coef[i]
 423                             blst_expand[i] = tmp[1]
 424                             blst_flg[i] = tmp[2]
 425                         else:
 426                             blst_coef[i]      = -blst_coef[i]
 427                             blst_flg[i]        = 1
 428                             blst_expand[i] = tmp[3]
 429                             blst_coef.append(-blst_coef[i]*tmp[0])
 430                             blst_expand.append(tmp[1])
 431                             blst_flg.append(tmp[2])
 432         (blst_coef,blst_expand) = MV.combine_common_factors(blst_coef,blst_expand)
 433         return(blst_coef,blst_expand)
 434     reduce_basis = staticmethod(reduce_basis)
 435
 436     def combine_common_factors(blst_coef,blst_expand):
 437         new_blst_coef = []
 438         new_blst_expand = []
 439         for i in range(MV.n1):
 440             new_blst_coef.append([])
 441             new_blst_expand.append([])
 442         nfac = len(blst_coef)
 443         for ifac in range(nfac):
 444             blen = len(blst_expand[ifac])
 445             new_blst_coef[blen].append(blst_coef[ifac])
 446             new_blst_expand[blen].append(blst_expand[ifac])
 447         for i in range(MV.n1):
 448             if len(new_blst_coef[i]) > 1:
 449                 MV.contract(new_blst_coef[i],new_blst_expand[i])
 450         return(new_blst_coef,new_blst_expand)
 451     combine_common_factors = staticmethod(combine_common_factors)
 452
 453     def contract(coefs,bases):
 454         dualsort(coefs,bases)
 455         n = len(bases)-1
 456         i = 0
 457         while i < n:
 458             j = i+1
 459             if bases[i] == bases[j]:
 460                 coefs[i] += coefs[j]
 461                 bases.pop(j)
 462                 coefs.pop(j)
 463                 n -= 1
 464             else:
 465                 i += 1
 466         n = len(coefs)
 467         i = 0
 468         while i < n:
 469             if coefs[i] == ZERO:
 470                 coefs.pop(i)
 471                 bases.pop(i)
 472                 n -= 1
 473             else:
 474                 i +=1
 475         return
 476     contract = staticmethod(contract)
 477
 478     def convert(coefs,bases):
 479         mv = MV()
 480         mv.bladeflg = 0
 481         for igrade in MV.n1rg:
 482             coef = coefs[igrade]
 483             base = bases[igrade]
 484             if len(coef) > 0:
 485                 nbases = MV.nbasis[igrade]
 486                 mv.mv[igrade] = numpy.array(nbases*[ZERO],dtype=numpy.object)
 487                 nbaserg = range(len(base))
 488                 for ibase in nbaserg:
 489                     if igrade > 0:
 490                         k = MV.basis[igrade].index(base[ibase])
 491                         mv.mv[igrade][k] = coef[ibase]
 492                     else:
 493                         mv.mv[igrade] = numpy.array([coef[0]],dtype=numpy.object)
 494         return(mv)
 495     convert = staticmethod(convert)
 496
 497     def set_str_format(str_mode=0):
 498         MV.str_mode = str_mode
 499         return
 500     set_str_format = staticmethod(set_str_format)
 501
 502     def str_rep(mv,lst_mode=0):
 503         """
 504         Converts internal representation of a multivector to a string
 505         for outputing.  If lst_mode = 1, str_rep outputs a list of
 506         strings where each string contains one multivector coefficient
 507         concatenated with the corressponding base or blade symbol.
 508         """
 509         if lst_mode:
 510             outlst = []
 511         if MV.bladeprint:
 512             mv.convert_to_blades()
 513             labels = MV.bladelabel
 514         else:
 515             if not mv.bladeflg:
 516                 labels = MV.basislabel
 517             else:
 518                 labels = MV.bladelabel
 519         value = ''
 520         for igrade in MV.n1rg:
 521             tmp = []
 522             if isinstance(mv.mv[igrade],numpy.ndarray):
 523                 j = 0
 524                 for x in mv.mv[igrade]:
 525                     if x != ZERO:
 526                         xstr = x.__str__()
 527                         if xstr == '+1' or xstr == '1' or xstr == '-1':
 528                             if xstr == '+1' or xstr == '1':
 529                                 xstr = '+'
 530                             else:
 531                                 xstr = '-'
 532                         else:
 533                             if xstr[0] != '+':
 534                                 xstr = '+{'+xstr+'}'
 535                             else:
 536                                 xstr = '+{'+xstr[1:]+'}'
 537                         value += xstr+labels[igrade][j]
 538                         if MV.str_mode and not lst_mode:
 539                             value += '\n'
 540                         if lst_mode:
 541                             tmp.append(value)
 542                     j += 1
 543             if lst_mode:
 544                 if len(tmp) > 0:
 545                     outlst.append(tmp)
 546                 value = ''
 547         if not lst_mode:
 548             if len(value) > 1 and value[0] == '+':
 549                 value = value[1:]
 550             if len(value) == 0:
 551                 value = '0'
 552         else:
 553             value = outlst
 554         if MV.latexflg:
 555             value = LaTeXstr(value)
 556         return(value)
 557     str_rep = staticmethod(str_rep)
 558
 559     def LaTeX():
 560         MV.latexflg = not MV.latexflg
 561         return
 562     LaTeX = staticmethod(LaTeX)
 563
 564     def setup(basis,metric='',debug=0):
 565         """
 566         MV.setup initializes the MV class by calculating the static
 567         multivector tables required for geometric algebra operations
 568         on multivectors.  See reference 5 section 2 for details on
 569         basis and metric arguments.
 570         """
 571         global MAIN_PROGRAM
 572         MV.latexflg = False
 573         MV.debug = debug
 574         MV.bladeprint = 0
 575         MV.tables_flg = 0
 576         MV.str_mode  = 0
 577         MV.rframe_flg = False
 578         if type(basis) == types.StringType:
 579             basislst = basis.split()
 580             MV.define_basis(basislst)
 581         if len(metric) > 0 and type(metric) == types.StringType:
 582             tmps = metric.split(',')
 583             metric = []
 584             for tmp in tmps:
 585                 xlst = tmp.split()
 586                 xtmp = []
 587                 for x in xlst:
 588                     xtmp.append(x)
 589                 metric.append(xtmp)
 590         MV.define_metric(metric)
 591         MV.multiplication_table()
 592         MV.blade_table()
 593         MV.inverse_blade_table()
 594         MV.tables_flg = 1
 595         isym = 0
 596         MV.bvec = []
 597         for name in MV.vbasis:
 598             bvar = MV(value=isym,mvtype='basisvector',mvname=name)
 599             bvar.bladeflg = 1
 600             MV.bvec.append(bvar)
 601             setattr(MAIN_PROGRAM,name,bvar)
 602             isym += 1
 603         return('Setup of '+basis+' complete!')
 604     setup = staticmethod(setup)
 605
 606     def print_blades():
 607         """
 608         Set multivector output to blade representation.
 609         """
 610         MV.bladeprint = 1
 611         return
 612     print_blades=staticmethod(print_blades)
 613
 614     def print_bases():
 615         """
 616         Set multivector output to base representation.
 617         """
 618         MV.bladeprint = 0
 619         return
 620     print_bases=staticmethod(print_bases)
 621
 622     def multiplication_table():
 623         """
 624         Calculate geometric product base multiplication table.
 625         See reference 5 section 3 for details.
 626         """
 627         MV.mtable = []
 628         for igrade in MV.n1rg:
 629             MV.mtable.append([])
 630             for ibase in range(MV.nbasis[igrade]):
 631                 MV.mtable[igrade].append([])
 632                 if igrade == 0:
 633                     base1 = []
 634                 else:
 635                     base1 = MV.basis[igrade][ibase]
 636                 for jgrade in MV.n1rg:
 637                     MV.mtable[igrade][ibase].append([])
 638                     for jbase in range(MV.nbasis[jgrade]):
 639                         if jgrade == 0:
 640                             base2 = []
 641                         else:
 642                             base2 = MV.basis[jgrade][jbase]
 643                         base = base1+base2
 644                         (coefs,bases) = MV.reduce_basis(base)
 645                         product = MV.convert(coefs,bases)
 646                         product.name = '('+MV.basislabel[igrade][ibase]+')('+MV.basislabel[jgrade][jbase]+')'
 647                         MV.mtable[igrade][ibase][jgrade].append(product)
 648         if MV.debug:
 649             print 'Multiplication Table:'
 650             for level1 in MV.mtable:
 651                 for level2 in level1:
 652                     for level3 in level2:
 653                         for mv in level3:
 654                             mv.printmv()
 655         return
 656     multiplication_table = staticmethod(multiplication_table)
 657
 658     def geometric_product(mv1,mv2):
 659         """
 660         MV.geometric_product(mv1,mv2) calculates the geometric
 661         product the multivectors mv1 and mv2 (mv1*mv2). See
 662         reference 5 section 3.
 663         """
 664         product = MV()
 665         if isinstance(mv1,MV) and isinstance(mv2,MV):
 666             bladeflg1 = mv1.bladeflg
 667             bladeflg2 = mv2.bladeflg
 668             if bladeflg1:
 669                 mv1.convert_from_blades()
 670             if bladeflg2:
 671                 mv2.convert_from_blades()
 672             mul = MV()
 673             for igrade in MV.n1rg:
 674                 gradei = mv1.mv[igrade]
 675                 if isinstance(gradei,numpy.ndarray):
 676                     for ibase in range(MV.nbasis[igrade]):
 677                         xi = gradei[ibase]
 678                         if xi != ZERO:
 679                             for jgrade in MV.n1rg:
 680                                 gradej = mv2.mv[jgrade]
 681                                 if isinstance(gradej,numpy.ndarray):
 682                                     for jbase in range(MV.nbasis[jgrade]):
 683                                         xj = gradej[jbase]
 684                                         if xj != ZERO:
 685                                             xixj = MV.mtable[igrade][ibase][jgrade][jbase].scalar_mul(xi*xj)
 686                                             product.add_in_place(xixj)
 687             product.bladeflg = 0
 688             if bladeflg1:
 689                 mv1.convert_to_blades()
 690             if bladeflg2:
 691                 mv1.convert_to_blades()
 692             if bladeflg1 and bladeflg2:
 693                 product.convert_to_blades()
 694         else:
 695             if isinstance(mv1,MV):
 696                 product = mv1.scalar_mul(mv2)
 697             else:
 698                 product = mv2.scalar_mul(mv1)
 699         return(product)
 700     geometric_product = staticmethod(geometric_product)
 701
 702     def wedge(igrade1,blade1,vector2,name=''):
 703         """
 704         Calculate the outer product of a multivector blade
 705         and a vector.  See reference 5 section 5 for details.
 706         """
 707         w12 = blade1*vector2
 708         w21 = vector2*blade1
 709         if igrade1%2 != 0:
 710             w = w12-w21
 711         else:
 712             w = w12+w21
 713         w.name = name
 714         return(w*HALF)
 715     wedge = staticmethod(wedge)
 716
 717     def blade_table():
 718         """
 719         Calculate basis blades in terms of bases. See reference 5
 720         section 5 for details. Used to convert from blade to base
 721         representation.
 722         """
 723         MV.btable = []
 724         MV.bladelabel = []
 725         basis_str = MV.basislabel[1]
 726         for igrade in MV.n1rg:
 727             MV.bladelabel.append([])
 728             if igrade == 0:
 729                 MV.bladelabel[0].append('1')
 730                 tmp = [MV(value=ONE,mvtype='scalar',mvname='1')]
 731             if igrade == 1:
 732                 tmp = []
 733                 for ibase in range(MV.nbasis[1]):
 734                     MV.bladelabel[1].append(basis_str[ibase])
 735                     tmp.append(MV(value=ibase,mvtype='basisvector',mvname=basis_str[ibase]))
 736             if igrade >= 2:
 737                 tmp = []
 738                 basis = MV.basis[igrade]
 739                 nblades = MV.nbasis[igrade]
 740                 iblade = 0
 741                 for blade in basis:
 742                     name = ''
 743                     for i in blade:
 744                         name += basis_str[i]+'^'
 745                     name = name[:-1]
 746                     MV.bladelabel[igrade].append(name)
 747                     lblade = MV.basis[igrade-1].index(blade[:-1])
 748                     rblade = blade[-1]
 749                     igrade1 = igrade-1
 750                     blade1  = MV.btable[igrade1][lblade]
 751                     vector2 = MV.btable[1][rblade]
 752                     b1Wv2 = MV.wedge(igrade1,blade1,vector2,name)
 753                     tmp.append(b1Wv2)
 754             MV.btable.append(tmp)
 755         if MV.debug:
 756             print 'Blade Tabel:'
 757             for grade in MV.btable:
 758                 for mv in grade:
 759                     print mv
 760             print 'Blade Labels:'
 761             print MV.bladelabel
 762         return
 763     blade_table = staticmethod(blade_table)
 764
 765     def inverse_blade_table():
 766         """
 767         Calculate bases in terms of basis blades. See reference 5
 768         section 5 for details. Used to convert from base to blade
 769         representation.
 770         """
 771         MV.ibtable = []
 772         for igrade in MV.n1rg:
 773             if igrade == 0:
 774                 tmp = [MV(value=ONE,mvtype='scalar',mvname='1')]
 775             if igrade == 1:
 776                 tmp = []
 777                 for ibase in range(MV.nbasis[1]):
 778                     tmp.append(MV(value=ibase,mvtype='basisvector'))
 779             if igrade >= 2:
 780                 tmp = []
 781                 iblade = 0
 782                 for blade in MV.btable[igrade]:
 783                     invblade = MV()
 784                     for i in range(igrade-1):
 785                         invblade.mv[i] = -blade.mv[i]
 786                         invblade.mv[igrade] = +blade.mv[igrade]
 787                         invblade.bladeflg = 1
 788                     if igrade >= 4:
 789                         jgrade = igrade-2
 790                         while jgrade > 1:
 791                             for ibase in range(MV.nbasis[jgrade]):
 792                                 invblade.substitute_base(jgrade,ibase,MV.ibtable[jgrade][ibase])
 793                             jgrade -= 2
 794                     invblade.name = MV.basislabel[igrade][iblade]
 795                     iblade += 1
 796                     tmp.append(invblade)
 797             MV.ibtable.append(tmp)
 798         if MV.debug:
 799             print 'Inverse Blade Tabel:'
 800             for grade in MV.ibtable:
 801                 for mv in grade:
 802                     mv.printmv()
 803         return
 804     inverse_blade_table = staticmethod(inverse_blade_table)
 805
 806     def outer_product(mv1,mv2):
 807         """
 808         MV.outer_product(mv1,mv2) calculates the outer (exterior,wedge)
 809         product of the multivectors mv1 and mv2 (mv1^mv2). See
 810         reference 5 section 6.
 811         """
 812         if isinstance(mv1,MV) and isinstance(mv2,MV):
 813             product = MV()
 814             product.bladeflg = 1
 815             mv1.convert_to_blades()
 816             mv2.convert_to_blades()
 817             for igrade1 in MV.n1rg:
 818                 if isinstance(mv1.mv[igrade1],numpy.ndarray):
 819                     pg1 = mv1.project(igrade1)
 820                     for igrade2 in MV.n1rg:
 821                         igrade = igrade1+igrade2
 822                         if igrade <= MV.n:
 823                             if isinstance(mv2.mv[igrade2],numpy.ndarray):
 824                                 pg2 = mv2.project(igrade2)
 825                                 pg1pg2 = pg1*pg2
 826                                 product.add_in_place(pg1pg2.project(igrade))
 827         else:
 828             if isinstance(mv1,MV):
 829                 product = mv1.scalar_mul(mv2)
 830             else:
 831                 product = mv2.scalar_mul(mv1)
 832         return(product)
 833     outer_product = staticmethod(outer_product)
 834
 835     def inner_product(mv1,mv2):
 836         """
 837         MV.inner_product(mv1,mv2) calculates the inner (scalar,dot)
 838         product of the multivectors mv1 and mv2 (mv1|mv2).  See
 839         reference 5 section 6.
 840         """
 841         if isinstance(mv1,MV) and isinstance(mv2,MV):
 842             product = MV()
 843             product.bladeflg = 1
 844             mv1.convert_to_blades()
 845             mv2.convert_to_blades()
 846             for igrade1 in range(1,MV.n1):
 847                 if isinstance(mv1.mv[igrade1],numpy.ndarray):
 848                     pg1 = mv1.project(igrade1)
 849                     for igrade2 in range(1,MV.n1):
 850                         igrade = (igrade1-igrade2).__abs__()
 851                         if isinstance(mv2.mv[igrade2],numpy.ndarray):
 852                             pg2 = mv2.project(igrade2)
 853                             pg1pg2 = pg1*pg2
 854                             product.add_in_place(pg1pg2.project(igrade))
 855             return(product)
 856         return(MV())
 857     inner_product = staticmethod(inner_product)
 858
 859     def addition(mv1,mv2):
 860         """
 861         MV.addition(mv1,mv2) calculates the sum
 862         of the multivectors mv1 and mv2 (mv1+mv2).
 863         """
 864         sum = MV()
 865         if isinstance(mv1,MV) and isinstance(mv2,MV):
 866             if mv1.bladeflg or mv2.bladeflg:
 867                 mv1.convert_to_blades()
 868                 mv2.convert_to_blades()
 869                 sum.bladeflg = 1
 870             for i in MV.n1rg:
 871                 if isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray):
 872                     sum.mv[i] = mv1.mv[i]+mv2.mv[i]
 873                 else:
 874                     if isinstance(mv1.mv[i],numpy.ndarray) and not isinstance(mv2.mv[i],numpy.ndarray):
 875                         sum.mv[i] = +mv1.mv[i]
 876                     else:
 877                         if isinstance(mv2.mv[i],numpy.ndarray) and not isinstance(mv1.mv[i],numpy.ndarray):
 878                             sum.mv[i] = +mv2.mv[i]
 879         else:
 880             if isinstance(mv1,MV):
 881                 sum = mv1.copy()
 882                 if isinstance(sum.mv[0],numpy,ndarray):
 883                     sum.mv[0] += numpy.array([mv2],dtype=numpy.object)
 884                 else:
 885                     sum.mv[0] = numpy.array([mv2],dtype=numpy.object)
 886             else:
 887                 sum = mv2.copy()
 888                 if isinstance(sum.mv[0],numpy.ndarray):
 889                     sum.mv[0] += numpy.array([mv1],dtype=numpy.object)
 890                 else:
 891                     sum.mv[0] = numpy.array([mv1],dtype=numpy.object)
 892         return(sum)
 893     addition = staticmethod(addition)
 894
 895     def subtraction(mv1,mv2):
 896         """
 897         MV.subtraction(mv1,mv2) calculates the difference
 898         of the multivectors mv1 and mv2 (mv1-mv2).
 899         """
 900         diff = MV()
 901         if isinstance(mv1,MV) and isinstance(mv2,MV):
 902             if mv1.bladeflg or mv2.bladeflg:
 903                 mv1.convert_to_blades()
 904                 mv2.convert_to_blades()
 905                 diff.bladeflg = 1
 906             for i in MV.n1rg:
 907                 if isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray):
 908                     diff.mv[i] = mv1.mv[i]-mv2.mv[i]
 909                 else:
 910                     if isinstance(mv1.mv[i],numpy.ndarray) and not isinstance(mv2.mv[i],numpy.ndarray):
 911                         diff.mv[i] = +mv1.mv[i]
 912                     else:
 913                         if not isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray):
 914                             diff.mv[i] = -mv2.mv[i]
 915         else:
 916             if isinstance(mv1,MV):
 917                 diff = mv1.copy()
 918                 if isinstandce(diff.mv[0],numpy.ndarray):
 919                     diff.mv[0] -= numpy.array([mv2],dtype=numpy.object)
 920                 else:
 921                     diff.mv[0] = -numpy.array([mv2],dtype=numpy.object)
 922             else:
 923                 diff = mv2.copy(1)
 924                 if isinstance(diff.mv[0],numpy.ndarray):
 925                     diff.mv[0] += numpy.array([mv1],dtype=numpy.object)
 926                 else:
 927                     diff.mv[0] = +numpy.array([mv1],dtype=numpy.object)
 928         return(diff)
 929     subtraction = staticmethod(subtraction)
 930
 931     def scalar_to_symbol(scalar):
 932         if isinstance(scalar,MV):
 933             return(scalar.mv[0][0])
 934         if type(scalar) == types.ListType:
 935             sym = []
 936             for x in scalar:
 937                 sym.append(MV.scalar_to_symbol(x))
 938             return(sym)
 939         return(scalar)
 940     scalar_to_symbol = staticmethod(scalar_to_symbol)
 941
 942     def __init__(self,value='',mvtype='',mvname=''):
 943         """
 944         Initialization of multivector X. Inputs are as follows
 945
 946         mvtype           value                       result
 947
 948         default         default                  Zero multivector
 949         'basisvector'   int i                    ith basis vector
 950         'basisbivector' int i                    ith basis bivector
 951         'scalar'        symbol x                 scalar of value x
 952         'grade'         [int i, symbol array A]  X.grade(i) = A
 953         'vector'        symbol array A           X.grade(1) = A
 954         'grade2'        symbol array A           X.grade(2) = A
 955
 956         mvname is name of multivector.
 957         """
 958         self.name      = mvname
 959         self.mv        = MV.n1*[0]
 960         self.bladeflg  = 0  #1 for blade expansion
 961         self.puregrade = 1
 962         if mvtype == 'basisvector':
 963             self.mv[1] = numpy.array(MV.nbasis[1]*[ZERO],dtype=numpy.object)
 964             self.mv[1][value] = ONE
 965         if mvtype == 'basisbivector':
 966             self.mv[2] = numpy.array(MV.nbasis[2]*[ZERO],dtype=numpy.object)
 967             self.mv[2][value] = ONE
 968         if mvtype == 'scalar':
 969             self.mv[0] = numpy.array([value],dtype=numpy.object)
 970         if mvtype == 'grade':
 971             igrade          = value[0]
 972             coefs           = value[1]
 973             coefs = MV.pad_zeros(coefs,MV.nbasis[igrade])
 974             self.mv[igrade] = numpy.array(coefs,dtype=numpy.object)
 975         if mvtype == 'vector':
 976             value = MV.pad_zeros(value,MV.nbasis[1])
 977             self.mv[1]    = numpy.array(value,dtype=numpy.object)
 978         if mvtype == 'grade2':
 979             value = MV.pad_zeros(value,MV.nbasis[2])
 980             self.mv[2]    = numpy.array(value,dtype=numpy.object)
 981
 982     def max_grade(self):
 983         """
 984         X.max_grade() is maximum grade of non-zero grades of X.
 985         """
 986         for i in range(MV.n,-1,-1):
 987             if isinstance(self.mv[i],numpy.ndarray):
 988                 return(i)
 989         return(-1)
 990
 991     def coord(xname,offset=0):
 992         xi_str = ''
 993         for i in MV.nrg:
 994             xi_str += xname+str(i+offset)+' '
 995         xi = make_symbols(xi_str)
 996         x = MV(xi,'vector')
 997         return(x)
 998     coord = staticmethod(coord)
 999
1000     def x(self,i):
1001         if isint(self.mv[1]):
1002             return(ZERO)
1003         return(self.mv[1][i])
1004
1005     def named(mvname,value='',mvtype=''):
1006         name = mvname
1007         tmp = MV(value=value,mvtype=mvtype,mvname=name)
1008         setattr(sys.modules[__name__],name,tmp)
1009         return
1010     named=staticmethod(named)
1011
1012     def printnm(tpl):
1013         for a in tpl:
1014             print a.name,' =',a.mv
1015         return
1016     printnm = staticmethod(printnm)
1017
1018     def __str__(self):
1019         """See MV.str_rep(self)"""
1020         return(MV.str_rep(self))
1021
1022     def printmv(self,name=''):
1023         title = ''
1024         if name:
1025             title += name+' = '
1026         else:
1027             if self.name:
1028                 title += self.name+' = '
1029         print title+MV.str_rep(self)
1030         return
1031
1032     def add_in_place(self,mv):
1033         """
1034         X.add_in_place(mv) increments multivector X by multivector
1035         mv.
1036         """
1037         if self.bladeflg or mv.bladeflg:
1038             self.convert_to_blades()
1039             mv.convert_to_blades()
1040         for i in MV.n1rg:
1041             if isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
1042                 self.mv[i] += mv.mv[i]
1043             else:
1044                 if not isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
1045                     self.mv[i] = +mv.mv[i]
1046         return
1047
1048     def sub_in_place(self,mv):
1049         """
1050         X.sub_in_place(mv) decrements multivector X by multivector
1051         mv.
1052         """
1053         if self.bladeflg or mv.bladeflg:
1054             self.convert_to_blades()
1055             mv.convert_to_blades()
1056         for i in MV.n1rg:
1057             if isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
1058                 self.mv[i] -= mv.mv[i]
1059             else:
1060                 if not isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
1061                     self.mv[i] = +mv.mv[i]
1062         return
1063
1064     def __pos__(self):
1065         p = MV()
1066         p.bladeflg = self.bladeflg
1067         for i in MV.n1rg:
1068             if isinstance(self.mv[i],numpy.ndarray):
1069                 p.mv[i] = +self.mv[i]
1070         return(p)
1071
1072     def __neg__(self):
1073         n = MV()
1074         n.bladeflg = self.bladeflg
1075         for i in MV.n1rg:
1076             if isinstance(self.mv[i],numpy.ndarray):
1077                 n.mv[i] = -self.mv[i]
1078         return(n)
1079
1080     def __add_ab__(self,mv):
1081         self.add_in_place(mv)
1082         return
1083
1084     def __sub_ab__(self,mv):
1085         self.sub_in_place(mv)
1086         return
1087
1088     def __add__(self,mv):
1089         """See MV.addition(self,mv)"""
1090         return(MV.addition(self,mv))
1091
1092     def __radd__(self,mv):
1093         """See MV.addition(mv,self)"""
1094         return(MV.addition(mv,self))
1095
1096     def __sub__(self,mv):
1097         """See MV.subtraction(self,mv)"""
1098         return(MV.subtraction(self,mv))
1099
1100     def __rsub__(self,mv):
1101         """See MV.subtraction(mv,self)"""
1102         return(MV.subtraction(mv,self))
1103
1104     def __xor__(self,mv):
1105         """See MV.outer_product(self,mv)"""
1106         return(MV.outer_product(self,mv))
1107
1108     def __pow__(self,mv):
1109         """See MV.outer_product(self,mv)"""
1110         return(MV.outer_product(self,mv))
1111
1112     def __rxor__(self,mv):
1113         """See MV.outer_product(mv,self)"""
1114         return(MV.outer_product(mv,self))
1115
1116     def __or__(self,mv):
1117         """See MV.inner_product(self,mv)"""
1118         return(MV.inner_product(self,mv))
1119
1120     def __ror__(self,mv):
1121         """See MV.inner_product(mv,self)"""
1122         return(MV.inner_product(mv,self))
1123
1124     def scalar_mul(self,c):
1125         """
1126         Y = X.scalar_mul(c), multiply multivector X by scalar c and return
1127         result.
1128         """
1129         mv = MV()
1130         mv.bladeflg = self.bladeflg
1131         mv.puregrade = self.puregrade
1132         for i in MV.n1rg:
1133             if isinstance(self.mv[i],numpy.ndarray):
1134                 mv.mv[i] = self.mv[i]*c
1135         return(mv)
1136
1137     def scalar_mul_inplace(self,c):
1138         """
1139         X.scalar_mul_inplace(c), multiply multivector X by scalar c and save
1140         result in X.
1141         """
1142         for i in MV.n1rg:
1143             if isinstance(self.mv[i],numpy.ndarray):
1144                 self.mv[i] = self.mv[i]*c
1145         return(mv)
1146
1147     def __mul__(self,mv):
1148         """See MV.geometric_product(self,mv)"""
1149         return(MV.geometric_product(self,mv))
1150
1151     def __rmul__(self,mv):
1152         """See MV.geometric_product(mv,self)"""
1153         return(MV.geometric_product(mv,self))
1154
1155     def __div__(self,scalar):
1156         div = MV()
1157         div.bladeflg = self.bladeflg
1158         for i in MV.n1rg:
1159             if isinstance(self.mv[i],numpy.ndarray):
1160                 div.mv[i] = self.mv[i]/scalar
1161         return(div)
1162
1163     def __div_ab__(self,scalar):
1164         for i in MV.n1rg:
1165             if isinstance(self.mv[i],numpy.ndarray):
1166                 self.mv[i] /= scalar
1167         return
1168
1169     def __call__(self,igrade=0,ibase=0):
1170         """
1171         X(i,j) returns symbol in ith grade and jth base or blade of
1172         multivector X.
1173         """
1174         if not isinstance(self.mv[igrade],numpy.ndarray):
1175             return(ZERO)
1176         return(self.mv[igrade][ibase])
1177
1178     def __eq__(self,mv):
1179         if not isinstance(mv,MV):
1180             return(False)
1181         for (mvi,mvj) in zip(self.mv,mv.mv):
1182             if isint(mvi) ^ isint(mvj):
1183                 return(False)
1184             if isinstance(mvi,numpy.ndarray) and isinstance(mvj,numpy.ndarray):
1185                 for (x,y) in zip(mvi,mvj):
1186                     if x != y:
1187                         return(False)
1188         return(True)
1189
1190     def copy(self,sub=0):
1191         """
1192         Y = X.copy(), make a deep copy of multivector X in multivector
1193         Y so that Y can be modified without affecting X.
1194         """
1195         cpy = MV()
1196         cpy.name      = self.name
1197         cpy.bladeflg  = self.bladeflg
1198         cpy.puregrade = self.puregrade
1199         for i in MV.n1rg:
1200             if sub:
1201                 if isinstance(self.mv[i],numpy.ndarray):
1202                     cpy.mv[i] = -self.mv[i]
1203             else:
1204                 if isinstance(self.mv[i],numpy.ndarray):
1205                     cpy.mv[i] = +self.mv[i]
1206         return(cpy)
1207
1208     def substitute_base(self,igrade,base,mv):
1209         if not isinstance(self.mv[igrade],numpy.ndarray):
1210             return
1211         if isinstance(base,numpy.ndarray):
1212             ibase = MV.basis[igrade].index(base)
1213         else:
1214             ibase = base
1215         coef = self.mv[igrade][ibase]
1216         if coef == ZERO:
1217             return
1218         self.mv[igrade][ibase] = ZERO
1219         self.add_in_place(mv*coef)
1220         return
1221
1222     def convert_to_blades(self):
1223         """
1224         X.convert_to_blades(), inplace convert base representation
1225         to blade representation.  See reference 5 section 5.
1226         """
1227         if self.bladeflg:
1228             return
1229         self.bladeflg = 1
1230         for igrade in range(2,MV.n1):
1231             if isinstance(self.mv[igrade],numpy.ndarray):
1232                 for ibase in range(MV.nbasis[igrade]):
1233                     coef = self.mv[igrade][ibase]
1234                     if (not coef == ZERO):
1235                         self.mv[igrade][ibase] = ZERO
1236                         self.add_in_place(MV.ibtable[igrade][ibase]*coef)
1237         return
1238
1239     def convert_from_blades(self):
1240         """
1241         X.convert_from_blades(), inplace convert blade representation
1242         to base representation.  See reference 5 section 5.
1243         """
1244         if not self.bladeflg:
1245             return
1246         self.bladeflg = 0
1247         for igrade in range(2,MV.n1):
1248             if isinstance(self.mv[igrade],numpy.ndarray):
1249                 for ibase in range(MV.nbasis[igrade]):
1250                     coef = self.mv[igrade][ibase]
1251                     if (not coef == ZERO):
1252                         self.mv[igrade][ibase] = ZERO
1253                         self.add_in_place(MV.btable[igrade][ibase]*coef)
1254         return
1255
1256     def project(self,r):
1257         """
1258         Grade projection operator. For multivector X, X.project(r)
1259         returns multivector of grade r components of X.
1260         """
1261         grade_r = MV()
1262         if r > MV.n:
1263             return(grade_r)
1264         self.convert_to_blades()
1265         if not isinstance(self.mv[r],numpy.ndarray):
1266             return(grade_r)
1267         grade_r.bladeflg = 1
1268         grade_r.puregrade = 1
1269         grade_r.mv[r] = +self.mv[r]
1270         return(grade_r)
1271
1272     def even(self):
1273         """
1274         Even grade projection operator. For multivector X, X.even()
1275         returns multivector of even grade components of X.
1276         """
1277         egrades = MV()
1278         self.convert_to_blades()
1279         egrades.bladeflg  = self.bladeflg
1280         egrades.puregrade = self.puregrade
1281         for igrade in range(0,MV.n1,2):
1282             egrades.mv[igrade] = +self.mv[igrade]
1283         return(egrades)
1284
1285     def odd(self):
1286         """
1287         Odd grade projection operator. For multivector X, X.odd()
1288         returns multivector of odd grade components of X.
1289         """
1290         ogrades = MV()
1291         self.convert_to_blades()
1292         ogrades.bladeflg  = self.bladeflg
1293         ogrades.puregrade = self.puregrade
1294         for igrade in range(1,MV.n1,2):
1295             ogrades.mv[igrade] = +self.mv[igrade]
1296         return(ogrades)
1297
1298     def rev(self):
1299         """
1300         Revisioin operator. For multivector X, X.rev()
1301         returns reversed multivector of X.
1302         """
1303         revmv = MV()
1304         self.convert_to_blades()
1305         revmv.bladeflg  = self.bladeflg
1306         revmv.puregrade = self.puregrade
1307         for igrade in MV.n1rg:
1308             if isinstance(self.mv[igrade],numpy.ndarray):
1309                 if igrade < 2 or (not (((igrade*(igrade-1))/2)%2)):
1310                     revmv.mv[igrade] = +self.mv[igrade]
1311                 else:
1312                     revmv.mv[igrade] = -self.mv[igrade]
1313         return(revmv)
1314
1315     def collect(self,lst):
1316         """
1317         Applies sympy/sympy collect function to each component
1318         of multivector.
1319         """
1320         for igrade in MV.n1rg:
1321             if isinstance(self.mv[igrade],numpy.ndarray):
1322                 for ibase in range(MV.nbasis[igrade]):
1323                     if self.mv[igrade][ibase] != ZERO:
1324                         self.mv[igrade][ibase] = \
1325                         collect(self.mv[igrade][ibase],lst)
1326                         self.mv[igrade][ibase] = \
1327                         collect(self.mv[igrade][ibase],lst)
1328         return
1329
1330     def sqrfree(self,lst):
1331         """
1332         Applies sympy/sympy sqrfree function to each component
1333         of multivector.
1334         """
1335         for igrade in MV.n1rg:
1336             if isinstance(self.mv[igrade],numpy.ndarray):
1337                 for ibase in range(MV.nbasis[igrade]):
1338                     if self.mv[igrade][ibase] != ZERO:
1339                         self.mv[igrade][ibase] = \
1340                         sqrfree(self.mv[igrade][ibase],lst)
1341         return
1342
1343     def collect(self,faclst):
1344         """
1345         Applies sympy collect_common_factors function
1346         to each component of multivector.
1347         """
1348         for igrade in MV.n1rg:
1349             if isinstance(self.mv[igrade],numpy.ndarray):
1350                 for ibase in range(MV.nbasis[igrade]):
1351                     if self.mv[igrade][ibase] != ZERO:
1352                         self.mv[igrade][ibase] = \
1353                         sympy.collect(self.mv[igrade][ibase],faclst)
1354         return
1355
1356     def subs(self,var,substitute):
1357         """
1358         Applies sympy subs function
1359         to each component of multivector.
1360         """
1361         for igrade in MV.n1rg:
1362             if isinstance(self.mv[igrade],numpy.ndarray):
1363                 for ibase in range(MV.nbasis[igrade]):
1364                     if self.mv[igrade][ibase] != ZERO:
1365                         self.mv[igrade][ibase] = \
1366                         self.mv[igrade][ibase].subs(var,substitute)
1367         return
1368
1369     def simplify(self):
1370         """
1371         Applies sympy subs function
1372         to each component of multivector.
1373         """
1374         for igrade in MV.n1rg:
1375             if isinstance(self.mv[igrade],numpy.ndarray):
1376                 for ibase in range(MV.nbasis[igrade]):
1377                     if self.mv[igrade][ibase] != ZERO:
1378                         self.mv[igrade][ibase] = \
1379                         sympy.simplify(self.mv[igrade][ibase])
1380         return
1381
1382     def cancel(self):
1383         """
1384         Applies sympy cancle function
1385         to each component of multivector.
1386         """
1387         for igrade in MV.n1rg:
1388             if isinstance(self.mv[igrade],numpy.ndarray):
1389                 for ibase in range(MV.nbasis[igrade]):
1390                     if self.mv[igrade][ibase] != ZERO:
1391                         self.mv[igrade][ibase] = \
1392                         sympy.cancel(self.mv[igrade][ibase])
1393         return
1394
1395     def trim(self):
1396         """
1397         Applies sympy cancle function
1398         to each component of multivector.
1399         """
1400         for igrade in MV.n1rg:
1401             if isinstance(self.mv[igrade],numpy.ndarray):
1402                 for ibase in range(MV.nbasis[igrade]):
1403                     if self.mv[igrade][ibase] != ZERO:
1404                         self.mv[igrade][ibase] = \
1405                         sympy.trim(self.mv[igrade][ibase])
1406         return
1407
1408     def expand(self):
1409         """
1410         Applies sympy/sympy expand function to each component
1411         of multivector.
1412         """
1413         for igrade in MV.n1rg:
1414             if isinstance(self.mv[igrade],numpy.ndarray):
1415                 for ibase in range(MV.nbasis[igrade]):
1416                     if self.mv[igrade][ibase] != ZERO:
1417                         self.mv[igrade][ibase] = \
1418                         self.mv[igrade][ibase].expand()
1419         return
1420
1421     def is_pure(self):
1422         igrade = -1
1423         ngrade = 0
1424         for i in MV.n1rg:
1425             if isinstance(self.mv[i],numpy.ndarray):
1426                 igrade = i
1427                 ngrade += 1
1428                 if ngrade > 1:
1429                     return(-1)
1430         return(igrade)
1431
1432     def compact(self):
1433         """
1434         Convert zero numpy arrays to single interge zero place holder
1435         in grade list for instanciated multivector. For example if
1436         numpy array of grade one components is a zero array then replace
1437         with single integer equal to zero.
1438         """
1439         for i in MV.n1rg:
1440             if not isint(self.mv[i]):
1441                 zero_flg = True
1442                 for x in self.mv[i]:
1443                     if x != 0:
1444                         zero_flg = False
1445                         break
1446                 if zero_flg:
1447                     self.mv[i] = 0
1448         return
1449
1450 ##class LT:
1451 ##    def __init__(self,fct_lst):
1452 ##        """
1453 ##        Initialize linear transformation.  fct_lst is a list of vectors A_{i}
1454 ##        such that A_{i} = LT(a_{i}) where the a_{i} are the basis vectos of
1455 ##        the MV class.
1456 ##        """
1457 ##        self.fct_lst = fct_lst
1458 ##        self.norm = ONE
1459 ##        self.ltblades = [[],self.fct_lst]
1460 ##        for igrade in range(2,MV.n1):
1461 ##            tmp = []
1462 ##            for ibasis in range(MV.nbasis[igrade]):
1463 ##                ilst = MV.basis[igrade][ibasis][:-1]
1464 ##                jlst = MV.basis[igrade-1]
1465 ##                iblade = jlst.index(ilst)
1466 ##                jblade = MV.basis[igrade][ibasis][-1]
1467 ##                tmp.append(self.ltblades[igrade-1][iblade]^self.fct_lst[jblade])
1468 ##            self.ltblades.append(tmp)
1469 ##        return
1470 ##
1471 ##    def Fvector(self,mv):
1472 ##        """
1473 ##        Return the linear transformation of an arbitrary vector.
1474 ##        """
1475 ##        mv.convert_to_blades()
1476 ##        if mv.is_pure() == 1:
1477 ##            fofa = MV()
1478 ##            fofa.bladeflg = 1
1479 ##            for i in range(MV.nbasis[1]):
1480 ##                fofa += mv.mv[1][i]*self.fct_lst[i]
1481 ##            fofa.collect_common_factors()
1482 ##            return(fofa/self.norm)
1483 ##        print 'Error in LT.F, mv =',mv,' not a vector!'
1484 ##        return
1485 ##
1486 ##    def __call__(self,mv):
1487 ##        """
1488 ##        Return the linear transformation of an arbitrary multivector.
1489 ##        """
1490 ##        if mv.is_pure() == 1:
1491 ##            return(self.Fvector(mv))
1492 ##        fofx = MV()
1493 ##        mv.convert_to_blades()
1494 ##        fofx.bladeflg = 1
1495 ##        for igrade in range(1,MV.n1):
1496 ##            if isinstance(mv.mv[igrade],numpy.ndarray):
1497 ##                for ibase in range(MV.nbasis[igrade]):
1498 ##                    if mv.mv[igrade][ibase] != ZERO:
1499 ##                        fofx.add_in_place(mv.mv[igrade][ibase]*self.ltblades[igrade][ibase])
1500 ##        fofx.collect_common_factors()
1501 ##        return(fofx/self.norm)
1502 ##
1503 ##    def __str__(self):
1504 ##        LTstr = ''
1505 ##        for ibasis in range(MV.n):
1506 ##            if self.norm != ONE:
1507 ##                LTstr += 'F('+MV.basislabel[1][ibasis]+') = ('
1508 ##            else:
1509 ##                LTstr += 'F('+MV.basislabel[1][ibasis]+') = '
1510 ##            LTstr += MV.str_rep(self.fct_lst[ibasis])
1511 ##            if self.norm != ONE:
1512 ##                LTstr += ')/('+str(self.norm)+')\n'
1513 ##            else:
1514 ##                LTstr += '\n'
1515 ##        return(LTstr[:-1])
1516 ##
1517 ##    def adj(self):
1518 ##        MV.define_reciprocal()
1519 ##        adj = []
1520 ##        norm = MV.Esq
1521 ##        norm = norm*self.norm
1522 ##        for ibasis in range(MV.n):
1523 ##            tmp = MV()
1524 ##            tmp.bladeflg = 1
1525 ##            for jbasis in range(MV.n):
1526 ##                tmp.add_in_place(MV.brecp[jbasis]*(MV.bvec[ibasis]|self.fct_lst[jbasis]))
1527 ##            tmp.expand()
1528 ##            tmp.sqrfree([])
1529 ##            adj.append(tmp)
1530 ##        adj = LT(adj)
1531 ##        adj.norm = norm
1532 ##        return(adj)