# Source code for galgebra.mv

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

import copy
import numbers
import operator
from functools import reduce
import sys
import warnings

from sympy import (
sin, cos, sinh, cosh, sqrt, trigsimp, expand,
simplify, diff, Rational, Expr, Abs, collect,
)
from sympy import exp as sympy_exp
from sympy import N as Nsympy

from . import printer
from . import metric
from . import utils
from .printer import ZERO_STR
from .utils import _KwargParser

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

half = Rational(1, 2)

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

[docs]class Mv(object):
"""
Wrapper class for multivector objects (self.obj) so that it is easy
to overload operators (*, ^, |, <, >)  for the various
multivector products and for printing.

Also provides a constructor to easily instantiate multivector objects.

of the multivector derivative have been added via the special vector
grad so that one can take the geometric derivative of a multivector
function A by applying grad from the left, grad*A, or the
right A*grad for both the left and right derivatives.  The operator
between the grad and the 'A' can be any of the multivector product
operators.

If f is a scalar function grad*f is the usual gradient of a function.
If A is a vector function grad|f is the divergence of A and
-I*(grad^A) is the curl of A (I is the pseudo scalar for the geometric
algebra)

Attributes
----------
obj : sympy.core.Expr
The underlying sympy expression for this multivector
"""

################### Multivector initialization #####################

fmt = 1
latex_flg = False
restore = False
dual_mode_lst = ['+I','I+','+Iinv','Iinv+','-I','I-','-Iinv','Iinv-']

[docs]    @staticmethod
def setup(ga):
"""
Set up constant multivectors required for multivector class for
a given geometric algebra, ga.
"""
Mv.fmt = 1
# copy basis in case the caller wanted to change it
return ga.mv_I, list(ga.mv_basis), ga.mv_x

@staticmethod
def Format(mode=1):
Mv.latex_flg = True
Mv.fmt = mode
return

[docs]    @staticmethod
def Mul(A, B, op):
"""
Function for all types of geometric multiplications called by
overloaded operators for *, ^, |, <, and >.
"""
if not isinstance(A, Mv):
A = B.Ga.mv(A)
if not isinstance(B, Mv):
B = A.Ga.mv(B)

if op == '*':
return A * B
elif op == '^':
return A ^ B
elif op == '|':
return A | B
elif op == '<':
return A < B
elif op == '>':
return A > B
else:
raise ValueError('Operation ' + op + 'not allowed in Mv.Mul!')

def characterise_Mv(self):
if self.char_Mv:
return
obj = expand(self.obj)
if isinstance(obj, numbers.Number):
return
if  obj.is_commutative:
return
args = obj.args
else:
self.char_Mv = True
return
else:
args = [obj]

#print 'args =', args
for term in args:
if term.is_commutative:
else:
c, nc = term.args_cnc(split_1=False)
else:
self.char_Mv = True
return
else:
self.char_Mv = True
return

# helper methods called by __init__. Note that these names must not change,
# as the part of the name after _make_ is public API via the string
# argument passed to __init__.
#
# The double underscores in argument names are to force the passing
# positionally. When python 3.8 is the lowest supported version, we can
# switch to using the / syntax from PEP570

@staticmethod
""" Make a pure grade multivector. """
if not s:
return root
return '{}__{}'.format(root, s)
if utils.isstr(__name_or_coeffs):
name = __name_or_coeffs
f = kw.pop('f', False)
kw.reject_remaining()
if isinstance(f, bool):
if f:  # Is a multivector function of all coordinates
return sum([Function(add_superscript(name, super_script), real=True)(*ga.coords) * base
else:  # Is a constant multivector function
return sum([Symbol(add_superscript(name, super_script), real=True) * base
else:  # Is a multivector function of tuple f variables
return sum([Function(add_superscript(name, super_script), real=True)(*f) * base
elif isinstance(__name_or_coeffs, (list, tuple)):
coeffs = __name_or_coeffs
kw.reject_remaining()
return sum([coef * base
else:
raise ValueError("Too many coefficients")
else:
raise TypeError("Expected a string, list, or tuple")

@staticmethod
def _make_scalar(ga, __name_or_value, **kwargs):
""" Make a scalar multivector """
if utils.isstr(__name_or_value):
name = __name_or_value
else:
value = __name_or_value
return value

@staticmethod
def _make_vector(ga, __name_or_coeffs, **kwargs):
""" Make a vector multivector """

@staticmethod
def _make_bivector(ga, __name_or_coeffs, **kwargs):
""" Make a bivector multivector """

@staticmethod
def _make_pseudo(ga, __name_or_coeffs, **kwargs):
""" Make a pseudo scalar multivector """

@staticmethod
def _make_mv(ga, __name, **kwargs):
""" Make a general (2**n components) multivector """
if not isinstance(__name, str):
raise TypeError("Must be a string")
for grade in range(ga.n + 1)
))

@staticmethod
def _make_spinor(ga, __name, **kwargs):
""" Make a general even (spinor) multivector """
if not isinstance(__name, str):
raise TypeError("Must be a string")
for grade in range(0, ga.n + 1, 2)
))

@staticmethod
def _make_odd(ga, __name, **kwargs):
""" Make a general odd multivector """
if not isinstance(__name, str):
raise TypeError("Must be a string")
for grade in range(1, ga.n + 1, 2)
), S(0))  # base case needed in case n == 0

# aliases
_make_even = _make_spinor

def __init__(self, *args, **kwargs):
"""
__init__(self, *args, ga, recp=None, **kwargs)

Note this constructor is overloaded, based on the type and number of
positional arguments:

.. class:: Mv(*, ga, recp=None)

Create a zero multivector
.. class:: Mv(expr, /, *, ga, recp=None)

Create a multivector from an existing vector or sympy expression
.. class:: Mv(coeffs, grade, /, ga, recp=None)

Create a multivector constant with a given grade
.. class:: Mv(name, category, /, *cat_args, ga, recp=None, f=False)

Create a multivector constant with a given category
.. class:: Mv(name, grade, /, ga, recp=None, f=False)

Create a multivector variable or function of a given grade
.. class:: Mv(coeffs, category, /, *cat_args, ga, recp=None)

Create a multivector variable or function of a given category

* and / in the signatures above are python
3.8 syntax, and respectively indicate the boundaries between
positional-only, normal, and keyword-only arguments.

Parameters
----------
ga : ~galgebra.ga.Ga
Geometric algebra to be used with multivectors
recp : object, optional
Normalization for reciprocal vector. Unused.
name : str
Name of this multivector, if it is a variable or function
coeffs : sequence
Sequence of coefficients for the given category.
This is only meaningful
category : str
One of:

* "grade" - this takes an additional argument, the grade to
create, in cat_args
* "scalar"
* "vector"
* "bivector" / "grade2"
* "pseudo"
* "mv"
* "even" / "spinor"
* "odd"

f : bool, tuple
True if function of coordinates, or a tuple of those coordinates.
Only valid if a name is passed

coords :
This argument is always accepted but ignored.

It is incorrectly described internally as the coordinates to be
used with multivector functions.
"""
kw = _KwargParser('__init__', kwargs)
self.Ga = kw.pop('ga')
self.recp = kw.pop('recp', None)  # not used
kw.pop('coords', None)  # ignored

self.char_Mv = False
self.versor_flg = None  # if is_versor is called flag is set
self.coords = self.Ga.coords
self.title = None

if len(args) == 0:  # default constructor 0
self.obj = S(0)
kw.reject_remaining()
elif len(args) == 1 and not utils.isstr(args[0]):  # copy constructor
x = args[0]
if isinstance(x, Mv):
self.obj = x.obj
else:
if isinstance(x, Expr):  #copy constructor for obj expression
self.obj = x
else:  #copy constructor for scalar obj expression
self.obj = S(x)
self.characterise_Mv()
kw.reject_remaining()
else:
if utils.isstr(args[1]):
make_args = list(args)
mode = make_args.pop(1)
make_func = getattr(Mv, '_make_{}'.format(mode), None)
if make_func is None:
raise ValueError('{!r} is not an allowed multivector type.'.format(mode))
self.obj = make_func(self.Ga, *make_args, **kwargs)
elif isinstance(args[1], int):  # args[1] = r (integer) Construct grade r multivector
if args[1] == 0:
# _make_scalar interprets its coefficient argument differently
make_args = list(args)
make_args.pop(1)
self.obj = Mv._make_scalar(self.Ga, *make_args, **kwargs)
else:
else:
raise TypeError("Expected string or int")

if utils.isstr(args[0]):
self.title = args[0]
self.characterise_Mv()

################# Multivector member functions #####################

# See Mv class functions documentation
self.characterise_Mv()
reflect = Mv(0,'scalar',ga=self.Ga)
else:
return reflect
else:

# See Mv class functions documentation
else:

def rotate_multivector(self,itheta,hint='-'):
Rm = (-itheta/S(2)).exp(hint)
Rp = (itheta/S(2)).exp(hint)
return Rm * self * Rp

[docs]    def base_rep(self):
""" Express as a linear combination of geometric products """
return self

b = copy.copy(self)
return b

""" Express as a linear combination of blades """
return self

b = copy.copy(self)
return b

def __hash__(self):
if self.is_scalar():
# ensure we match equality
return hash(self.obj)
else:
return hash((self.Ga, self.obj))

def __eq__(self, A):
if isinstance(A, Mv):
diff = (self - A).expand().simplify()
#diff = (self - A).expand()
if diff.obj == S(0):
return True
else:
return False
else:
if self.is_scalar() and self.obj == A:
return True
else:
return False

if sys.version_info.major < 3:
def __ne__(self, other):
return not (self == other)

"""
def __eq__(self, A):
if not isinstance(A, Mv):
if not self.is_scalar():
return False
if expand(self.obj) == expand(A):
return True
else:
return False
coefs, bases = metric.linear_expand(self.obj)
Acoefs, Abases = metric.linear_expand(A.obj)
if len(bases) != len(Abases):
return False
if set(bases) != set(Abases):
return False
for base in bases:
index = bases.index(base)
indexA = Abases.index(base)
if expand(coefs[index]) != expand(Acoefs[index]):
return False
return True
"""

def __neg__(self):
return Mv(-self.obj, ga=self.Ga)

if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):
return Mv(self.obj + A, ga=self.Ga)

if self.Ga != A.Ga:
raise ValueError('In + operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):

return Mv(self.obj + A.obj, ga=self.Ga)
else:
else:
return Mv(self.obj + A.obj, ga=self.Ga)

return(self + A)

def __sub__(self, A):

if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):
return Mv(self.obj - A, ga=self.Ga)

if self.Ga != A.Ga:
raise ValueError('In - operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):

return Mv(self.obj - A.obj, ga=self.Ga)
else:
else:
return Mv(self.obj - A.obj, ga=self.Ga)

def __rsub__(self, A):
return -self + A

def __mul__(self, A):

if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):
return Mv(expand(A * self.obj), ga=self.Ga)

if self.Ga != A.Ga:
raise ValueError('In * operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):
return A.Mul(self, A, op='*')

if self.is_scalar():
return Mv(self.obj * A, ga=self.Ga)

self = self.base_rep()
A = A.base_rep()

selfxA = Mv(self.Ga.mul(self.obj, A.obj), ga=self.Ga)

self = self.base_rep()

selfxA = Mv(self.Ga.mul(self.obj, A.obj), ga=self.Ga)

A = A.base_rep()

selfxA = Mv(self.Ga.mul(self.obj, A.obj), ga=self.Ga)
else:
return Mv(self.Ga.mul(self.obj, A.obj), ga=self.Ga)

def __rmul__(self, A):
return Mv(expand(A * self.obj), ga=self.Ga)

def __truediv__(self, A):
if isinstance(A,Mv):
return self * A.inv()
else:
return self * (S(1)/A)

if sys.version_info.major < 3:
__div__ = __truediv__

def __str__(self):
if printer.GaLatexPrinter.latex_flg:
Printer = printer.GaLatexPrinter
else:
Printer = printer.GaPrinter
return Printer().doprint(self)

def __repr__(self):
return str(self)

def __getitem__(self,key):
'''
get a specified grade of a multivector
'''

def Mv_str(self):
global print_replace_old, print_replace_new
return str(self.obj)

# note: this just replaces self for the rest of this function
obj = expand(self.obj)
obj = metric.Simp.apply(obj)
self = Mv(obj, ga=self.Ga)

else:
base_keys = self.Ga._all_bases_lst
if isinstance(self.obj, Add):  # collect coefficients of bases
if self.obj.is_commutative:
return self.obj
args = self.obj.args
terms = {}  # dictionary with base indexes as keys
for arg in args:
c, nc = arg.args_cnc()
c = reduce(operator.mul, c, S(1))
if len(nc) > 0:
base = nc[0]
if base in base_keys:
index = base_keys.index(base)
if index in terms:
(c_tmp, base, g_keys) = terms[index]
terms[index] = (c_tmp + c, base, g_keys)
else:
else:
terms = list(terms.items())
sorted_terms = sorted(terms, key=operator.itemgetter(0))  # sort via base indexes

s = str(sorted_terms[0][1][0] * sorted_terms[0][1][1])
if printer.GaPrinter.fmt == 3:
s = ' ' + s + '\n'
if printer.GaPrinter.fmt == 2:
s = ' ' + s
for (key, (c, base, grade)) in sorted_terms[1:]:
term = str(c * base)
s += '\n'
if term[0] == '-':
term = ' - ' + term[1:]
else:
term = ' + ' + term
if printer.GaPrinter.fmt == 3:  # one base per line
s += term + '\n'
else:  # one multivector per line
s += term
if s[-1] == '\n':
s = s[:-1]
if printer.print_replace_old is not None:
s = s.replace(printer.print_replace_old,printer.print_replace_new)
return s
else:
return str(self.obj)

def Mv_latex_str(self):

if self.obj == 0:
return ZERO_STR

# todo: use the nonlocal keyword here instead once we drop python 2
class first_line:
value = True

def append_plus(c_str):
if first_line.value:
first_line.value = False
return c_str
else:
c_str = c_str.strip()
if c_str[0] == '-':
return ' ' + c_str
else:
return ' + ' + c_str

# str representation of multivector
# note: this just replaces self for the rest of this function
obj = expand(self.obj)
obj = metric.Simp.apply(obj)
self = Mv(obj, ga=self.Ga)

if self.obj == S(0):
return ZERO_STR

else:
base_keys = self.Ga._all_bases_lst
args = self.obj.args
else:
args = [self.obj]
terms = {}  # dictionary with base indexes as keys
for arg in args:
c, nc = arg.args_cnc(split_1=False)
c = reduce(operator.mul, c, S(1))
if len(nc) > 0:
base = nc[0]
if base in base_keys:
index = base_keys.index(base)
if index in terms:
(c_tmp, base, g_keys) = terms[index]
terms[index] = (c_tmp + c, base, g_keys)
else:
else:
terms = list(terms.items())

sorted_terms = sorted(terms, key=operator.itemgetter(0))  # sort via base indexes

if len(sorted_terms) == 1 and sorted_terms[0][1][2] == 0:  # scalar
return printer.latex(printer.coef_simplify(sorted_terms[0][1][0]))

lines = []
s = ''
for (index, (coef, base, grade)) in sorted_terms:
coef = printer.coef_simplify(coef)
#coef = simplify(coef)
l_coef = printer.latex(coef)
if l_coef == '1' and base != S(1):
l_coef = ''
if l_coef == '-1' and base != S(1):
l_coef = '-'
if base == S(1):
l_base = ''
else:
l_base = printer.latex(base)
cb_str = '\\left ( ' + l_coef + '\\right ) ' + l_base
else:
cb_str = l_coef + ' ' + l_base
if printer.GaLatexPrinter.fmt == 3:  # One base per line
lines.append(append_plus(cb_str))
elif printer.GaLatexPrinter.fmt == 2:  # One grade per line
if not first_line.value:
lines.append(s)
s = append_plus(cb_str)
else:
s += append_plus(cb_str)
else:  # One multivector per line
s += append_plus(cb_str)
if printer.GaLatexPrinter.fmt == 2:
lines.append(s)
if printer.GaLatexPrinter.fmt >= 2:
if len(lines) == 1:
return lines[0]
s = ' \\begin{align*} '
for line in lines:
s += ' & ' + line + ' \\\\ '
s = s[:-3] + ' \\end{align*} \n'
return s

def __xor__(self, A):  # wedge (^) product

if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):
return Mv(A * self.obj, ga=self.Ga)

if self.Ga != A.Ga:
raise ValueError('In ^ operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):
return A.Mul(self, A, op='^')

if self.is_scalar():
return self * A

self_W_A = self.Ga.wedge(self.obj, A.obj)
self_W_A = Mv(self_W_A, ga=self.Ga)
return self_W_A

def __rxor__(self, A):  # wedge (^) product
if not isinstance(A, Mv):
return Mv(A * self.obj, ga=self.Ga)
else:
return A * self

def __or__(self, A):  # dot (|) product
if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):
return Mv(ga=self.Ga)

if self.Ga != A.Ga:
raise ValueError('In | operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):
return A.Mul(self, A, op='|')

return Mv(self.Ga.hestenes_dot(self.obj, A.obj), ga=self.Ga)

def __ror__(self, A):  # dot (|) product
if not isinstance(A, Mv):
return Mv(ga=self.Ga)
else:
return A | self

def __pow__(self,n):  # Integer power operator
if not isinstance(n,int):
raise ValueError('!!!!Multivector power can only be to integer power!!!!')

result = S(1)
for x in range(n):
result *= self
return result

def __lshift__(self, A): # anti-comutator (<<)
return half * (self * A + A * self)

def __rshift__(self, A): # comutator (>>)
return half * (self * A - A * self)

def __rlshift__(self, A): # anti-comutator (<<)
return half * (A * self + self * A)

def __rrshift__(self, A): # comutator (>>)
return half * (A * self - self * A)

def __lt__(self, A):  # left contraction (<)

if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):  # sympy scalar
return Mv(A * self.obj, ga=self.Ga)

if self.Ga != A.Ga:
raise ValueError('In < operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):
return A.Mul(self, A, op='<')

"""
if A.is_scalar():
if self.is_scalar():
return self.obj * A.obj
else:
return S(0)
"""

self_lc_A = Mv(self.Ga.left_contract(self.obj, A.obj), ga=self.Ga)
return self_lc_A

def __gt__(self, A):  # right contraction (>)

if (not isinstance(A, Mv)) and (not isinstance(A, Dop)):  # sympy scalar
return self.Ga.mv(A * self.scalar())

if self.Ga != A.Ga:
raise ValueError('In > operation Mv arguments are not from same geometric algebra')

if isinstance(A, Dop):
return A.Mul(self, A, op='>')

"""
if self.is_scalar():
if A.is_scalar():
return self.obj * A.obj
else:
return S(0)
"""

self_rc_A = Mv(self.Ga.right_contract(self.obj, A.obj), ga=self.Ga)
return self_rc_A

[docs]    def collect(self,deep=False):
"""
group coeffients of blades of multivector
so there is only one coefficient per grade
"""
self.obj = expand(self.obj)
else:
c = self.Ga.bases_lst
self.obj = self.obj.collect(c)
return self
"""
obj_dict = {}
for coef, base in metric.linear_expand_terms(self.obj):
if base in list(obj_dict.keys()):
obj_dict[base] += coef
else:
obj_dict[base] = coef
obj = S(0)
for base in list(obj_dict.keys()):
if deep:
obj += collect(obj_dict[base])*base
else:
obj += obj_dict[base]*base
return Mv(obj, ga=self.Ga)

def is_scalar(self):
return True
else:
return False

def is_vector(self):
return True
else:
return False

"""
True is self is blade, otherwise False
"""
else:
if self.is_versor():
else:
else:

def is_base(self):
(coefs, _bases) = metric.linear_expand(self.obj)
if len(coefs) > 1:
return False
else:
return coefs[0] == ONE

[docs]    def is_versor(self):
"""
Test for versor (geometric product of vectors)

This follows Leo Dorst's test for a versor.
Leo Dorst, 'Geometric Algebra for Computer Science,' p.533
Sets self.versor_flg and returns value
"""

if self.versor_flg is not None:
return self.versor_flg
self.characterise_Mv()
self.versor_flg = False
self_rev = self.rev()
# see if self*self.rev() is a scalar
test = self*self_rev
if not test.is_scalar():
return self.versor_flg
# see if self*x*self.rev() returns a vector for x an arbitrary vector
test = self * self.Ga._XOX * self.rev()
self.versor_flg = test.is_vector()
return self.versor_flg

def is_zero(self):
if self.obj == 0:
return True
return False

[docs]    def scalar(self):
""" return scalar part of multivector as sympy expression """
return self.Ga.scalar_part(self.obj)

""" return r-th grade of multivector as a multivector """

def components(self):
cb = metric.linear_expand_terms(self.obj)
cb = sorted(cb, key=lambda x: self.Ga._all_blades_lst.index(x[1]))
return [self.Ga.mv(coef * base) for (coef, base) in cb]

cb = metric.linear_expand_terms(self.obj)
(coefs, bases) = list(zip(*cb))
return coefs

"""
For a multivector, A, and a list of basis blades, blade_lst return
a list (sympy expressions) of the coefficients of each basis blade
"""

(coefs, bases) = metric.linear_expand(self.obj)
coef_lst = []
else:
coef_lst.append(ZERO)
return coef_lst

[docs]    def proj(self, bases_lst):
"""
Project multivector onto a given list of bases.  That is find the
part of multivector with the same bases as in the bases_lst.
"""
bases_lst = [x.obj for x in bases_lst]
obj = 0
for coef, base in metric.linear_expand_terms(self.obj):
if base in bases_lst:
obj += coef * base
return Mv(obj, ga=self.Ga)

def dual(self):
mode = self.Ga.dual_mode_value
sign = S(1)
if '-' in mode:
sign = -sign
if 'Iinv' in mode:
I = self.Ga.i_inv
else:
I = self.Ga.i
if mode[0] == '+' or mode[0] == '-':
return sign * I * self
else:
return sign * self * I

[docs]    def even(self):
""" return even parts of multivector """
return Mv(self.Ga.even_odd(self.obj, True), ga=self.Ga)

[docs]    def odd(self):
""" return odd parts of multivector """
return Mv(self.Ga.even_odd(self.obj, False), ga=self.Ga)

def rev(self):
return Mv(self.Ga.reverse(self.obj), ga=self.Ga)

__invert__ = rev # allow ~x to call x.rev()

def diff(self, coord):
if self.Ga.coords is None:
obj = diff(self.obj, coord)
elif coord not in self.Ga.coords:
if self.Ga.par_coords is None:
obj = diff(self.obj, coord)
elif coord not in self.Ga.par_coords:
obj = diff(self.obj, coord)
else:
obj = diff(self.obj, coord)
for x_coord in self.Ga.coords:
f = self.Ga.par_coords[x_coord]
if f != S(0):
tmp1 = self.Ga.pDiff(self.obj, x_coord)
tmp2 = diff(f, coord)
obj += tmp1 * tmp2
else:
obj = self.Ga.pDiff(self.obj, coord)
return Mv(obj, ga=self.Ga)

def pdiff(self, var):
return Mv(self.Ga.pDiff(self.obj, var), ga=self.Ga)

[docs]    def Grad(self, coords, mode='*', left=True):
"""
Returns various derivatives (*,^,|,<,>) of multivector functions
with respect to arbitrary coordinates, 'coords'.  This would be
used where you have a multivector function of both the basis
coordinate set and and auxiliary coordinate set.  Consider for
example a linear transformation in which the matrix coefficients
depend upon the manifold coordinates, but the vector being
transformed does not and you wish to take the divergence of the
linear transformation with respect to the linear argument.
"""
return Mv(self.Ga.Diff(self, mode, left, coords=coords), ga=self.Ga)

[docs]    def exp(self, hint='-'):  # Calculate exponential of multivector
"""
Only works if square of multivector is a scalar.  If square is a
number we can determine if square is > or < zero and hence if
one should use trig or hyperbolic functions in expansion.  If
square is not a number use 'hint' to determine which type of
functions to use in expansion
"""
self_sq = self * self
if self_sq.is_scalar():
sq = simplify(self_sq.obj)  # sympy expression for self**2
if sq == S(0):  # sympy expression for self**2 = 0
return self + S(1)
(coefs,bases) = metric.linear_expand(self.obj)
if len(coefs) == 1:  # Exponential of scalar * base
base = bases[0]
base_Mv = self.Ga.mv(base)
base_sq = (base_Mv*base_Mv).scalar()
if hint == '-': # base^2 < 0
base_n = sqrt(-base_sq)
return self.Ga.mv(cos(base_n*coefs[0]) + sin(base_n*coefs[0])*(bases[0]/base_n))
else:  # base^2 > 0
base_n = sqrt(base_sq)
return self.Ga.mv(cosh(base_n*coefs[0]) + sinh(base_n*coefs[0])*(bases[0]/base_n))
if sq.is_number:  # Square is number, can test for sign
if sq > S(0):
norm = sqrt(sq)
value = self.obj / norm
tmp = Mv(cosh(norm) + sinh(norm) * value, ga=self.Ga)
return tmp
else:
norm = sqrt(-sq)
value = self.obj / norm
tmp = Mv(cos(norm) + sin(norm) * value, ga=self.Ga)
return tmp
else:
if hint == '+':
norm = simplify(sqrt(sq))
value = self.obj / norm
tmp = Mv(cosh(norm) + sinh(norm) * value, ga=self.Ga)
return tmp
else:
norm = simplify(sqrt(-sq))
value = self.obj / norm
obj = cos(norm) + sin(norm) * value
tmp = Mv(cos(norm) + sin(norm) * value, ga=self.Ga)
return tmp
else:
raise ValueError('"' + str(self) + '**2" is not a scalar in exp.')

else:
(coefs, bases) = metric.linear_expand(self.obj)
bases_lst = list(bases)  # python 2.5
if base in bases:
self.obj += (value - coefs[bases_lst.index(base)]) * base
else:
self.obj += value * base
return

[docs]    def Fmt(self, fmt=1, title=None):
"""
Set format for printing of multivectors

* fmt=1 - One multivector per line
* fmt=2 - One grade per line
* fmt=3 - one base per line

Usage for multivector A example is::

A.Fmt('2','A')

output is::

'A = '+str(A)

with one grade per line.  Works for both standard printing and
for latex.
"""
if printer.GaLatexPrinter.latex_flg:
printer.GaLatexPrinter.prev_fmt = printer.GaLatexPrinter.fmt
printer.GaLatexPrinter.fmt = fmt
else:
printer.GaPrinter.prev_fmt = printer.GaPrinter.fmt
printer.GaPrinter.fmt = fmt

if title is not None:
self.title = title

if printer.isinteractive():
return self

if Mv.latex_flg:
latex_str = printer.GaLatexPrinter.latex(self)
printer.GaLatexPrinter.fmt = printer.GaLatexPrinter.prev_fmt

if title is not None:
return title + ' = ' + latex_str
else:
return latex_str
else:
s = str(self)
printer.GaPrinter.fmt = printer.GaPrinter.prev_fmt
if title is not None:
return title + ' = ' + s
else:
return s

def _repr_latex_(self):
latex_str = printer.GaLatexPrinter.latex(self)
if r'\begin{align*}' not in latex_str:
if self.title is None:
latex_str = r'\begin{equation*} ' + latex_str + r' \end{equation*}'
else:
latex_str = r'\begin{equation*} ' + self.title + ' = ' + latex_str + r' \end{equation*}'
else:
if self.title is not None:
latex_str = latex_str.replace('&',' ' + self.title + ' =&',1)
return latex_str

def norm2(self):
reverse = self.rev()
product = self * reverse
if product.is_scalar():
return product.scalar()
else:
raise TypeError('"(' + str(product) + ')**2" is not a scalar in norm2.')

[docs]    def norm(self, hint='+'):
"""
If A is a multivector and A*A.rev() is a scalar then::

A.norm() == sqrt(Abs(A*A.rev()))

The problem in simplifying the norm is that if A is symbolic
you don't know if A*A.rev() is positive or negative. The use
of the hint argument is as follows:

=======  ========================
hint     A.norm()
=======  ========================
'+'  sqrt(A*A.rev())
'-'  sqrt(-A*A.rev())
'0'  sqrt(Abs(A*A.rev()))
=======  ========================

The default hint='+' is correct for vectors in a Euclidean vector
space.  For bivectors in a Euclidean vector space use hint='-'. In
a mixed signature space all bets are off for the norms of symbolic
expressions.
"""
reverse = self.rev()
product = self * reverse

if product.is_scalar():
product = product.scalar()
if product.is_number:
if product >= S(0):
return sqrt(product)
else:
return sqrt(-product)
else:
if hint == '+':
return metric.square_root_of_expr(product)
elif hint == '-':
return metric.square_root_of_expr(-product)
else:
return sqrt(Abs(product))
else:
raise TypeError('"(' + str(product) + ')" is not a scalar in norm.')

__abs__ = norm # allow abs(x) to call z.norm()

def inv(self):
if self.is_scalar():  # self is a scalar
return self.Ga.mv(S(1)/self.obj)
self_sq = self * self
if self_sq.is_scalar():  # self*self is a scalar
"""
if self_sq.scalar() == S(0):
raise ValueError('!!!!In multivector inverse, A*A is zero!!!!')
"""
return (S(1)/self_sq.obj)*self
self_rev = self.rev()
self_self_rev = self * self_rev
if(self_self_rev.is_scalar()): # self*self.rev() is a scalar
"""
if self_self_rev.scalar() == S(0):
raise ValueError('!!!!In multivector inverse A*A.rev() is zero!!!!')
"""
return (S(1)/self_self_rev.obj) * self_rev
raise TypeError('In inv() for self =' + str(self) + 'self, or self*self or self*self.rev() is not a scalar')

def func(self, fct):  # Apply function, fct, to each coefficient of multivector
s = S(0)
for coef, base in metric.linear_expand_terms(self.obj):
s += fct(coef) * base
fct_self = Mv(s, ga=self.Ga)
fct_self.characterise_Mv()
return fct_self

def trigsimp(self):
return self.func(trigsimp)

def simplify(self, modes=simplify):
if not isinstance(modes, (list, tuple)):
modes = [modes]

obj = S(0)
for coef, base in metric.linear_expand_terms(self.obj):
for mode in modes:
coef = mode(coef)
obj += coef * base
return Mv(obj, ga=self.Ga)

def subs(self, d):
# For each scalar coef of the multivector apply substitution argument d
obj = sum((
coef.subs(d) * base for coef, base in metric.linear_expand_terms(self.obj)
), S(0))
return Mv(obj, ga=self.Ga)

def expand(self):
obj = sum((
expand(coef) * base for coef, base in metric.linear_expand_terms(self.obj)
), S(0))
return Mv(obj, ga=self.Ga)

def list(self):
indexes = []
key_coefs = []
for coef, base in metric.linear_expand_terms(self.obj):
if base in self.Ga.basis:
index = self.Ga.basis.index(base)
key_coefs.append((coef, index))
indexes.append(index)

for index in self.Ga.n_range:
if index not in indexes:
key_coefs.append((S(0), index))

key_coefs = sorted(key_coefs, key=operator.itemgetter(1))
coefs = [x[0] for x in key_coefs]
return coefs

"""
"""
self.characterise_Mv()

[docs]def compare(A,B):
"""
Determine if B = c*A where c is a scalar.  If true return c
otherwise return 0.
"""
if isinstance(A, Mv) and isinstance(B, Mv):
Acoefs, Abases = metric.linear_expand(A.obj)
Bcoefs, Bbases = metric.linear_expand(B.obj)
if len(Acoefs) != len(Bcoefs):
return 0
if Abases != Bbases:
return 0
if Bcoefs[0] != 0 and Abases[0] == Bbases[0]:
c = simplify(Acoefs[0]/Bcoefs[0])
print('c =',c)
else:
return 0
for acoef,abase,bcoef,bbase in zip(Acoefs[1:],Abases[1:],Bcoefs[1:],Bbases[1:]):
print(acoef,'\n',abase,'\n',bcoef,'\n',bbase)
if bcoef != 0 and abase == bbase:
print('c-a/b =',simplify(c-(acoef/bcoef)))
if simplify(acoef/bcoef) != c:
return 0
else:
pass
else:
return 0
return c
else:
raise TypeError('In compare both arguments are not multivectors\n')

def _consolidate_terms(terms):
"""
Remove zero coefs and consolidate coefs with repeated pdiffs.
"""
new_coefs = []
new_pdiffs = []
for (coef, pd) in terms:
if coef != S(0):
if pd in new_pdiffs:
index = new_pdiffs.index(pd)
new_coefs[index] += coef
else:
new_coefs.append(coef)
new_pdiffs.append(pd)
return tuple(zip(new_coefs, new_pdiffs))

def _merge_terms(terms1, terms2):
""" Concatenate and consolidate two sets of already-consolidated terms """
pdiffs1 = [pdiff for _, pdiff in terms1]
pdiffs2 = [pdiff for _, pdiff in terms2]

pdiffs = pdiffs1 + [x for x in pdiffs2 if x not in pdiffs1]
coefs = len(pdiffs) * [S(0)]

for coef, pdiff in terms1:
index = pdiffs.index(pdiff)
coefs[index] += coef

for coef, pdiff in terms2:
index = pdiffs.index(pdiff)
coefs[index] += coef

# remove zeros
return [(coef, pdiff) for coef, pdiff in zip(coefs, pdiffs) if coef != S(0)]

################ Scalar Partial Differential Operator Class ############

[docs]class Sdop(object):
"""
Scalar differential operator is of the form (Einstein summation)

.. math:: D = c_{i}*D_{i}

where the :math:c_{i}'s are scalar coefficient (they could be functions)
and the :math:D_{i}'s are partial differential operators (:class:Pdop).

Attributes
----------
terms : tuple of tuple
the structure :math:((c_{1},D_{1}),(c_{2},D_{2}), ...)
"""

init_slots = {'ga': (None, 'Associated geometric algebra')}

str_mode = False

def TSimplify(self):
return Sdop([
(metric.Simp.apply(coef), pdiff) for (coef, pdiff) in self.terms
], ga=self.Ga)

[docs]    @staticmethod
def consolidate_coefs(sdop):
"""
Remove zero coefs and consolidate coefs with repeated pdiffs.
"""
if isinstance(sdop, Sdop):
return Sdop(_consolidate_terms(sdop.terms), ga=sdop.Ga)
else:
return _consolidate_terms(sdop)

def simplify(self, modes=simplify):
return Sdop([
(metric.apply_function_list(modes, coef), pdiff)
for coef, pdiff in self.terms
], ga=self.Ga)

def _with_sorted_terms(self):
new_terms = sorted(self.terms, key=lambda term: Pdop.sort_key(term[1]))
return Sdop(new_terms, ga=self.Ga)

def Sdop_str(self):
if len(self.terms) == 0:
return ZERO_STR

self = self._with_sorted_terms()
s = ''
for (coef, pdop) in self.terms:
coef_str = printer.latex(coef)
pd_str = printer.latex(pdop)

if coef == S(1):
s += pd_str
elif coef == S(-1):
s += '-' + pd_str
else:
s += '(' + coef_str + ')*' + pd_str
else:
s += coef_str + '*' + pd_str
s += ' + '

s = s.replace('+ -','- ')
s = s[:-3]
if Sdop.str_mode:
if len(self.terms) > 1 or isinstance(self.terms[0][0], Add):
s = '(' + s + ')'
return s

def Sdop_latex_str(self):
if len(self.terms) == 0:
return ZERO_STR

self = self._with_sorted_terms()

s = ''
for (coef, pdop) in self.terms:
coef_str = printer.latex(coef)
pd_str = printer.latex(pdop)
if coef == S(1):
if pd_str == '':
s += '1'
else:
s += pd_str
elif coef == S(-1):
if pd_str == '':
s += '-1'
else:
s += '-' + pd_str
else:
s += r'\left ( ' + coef_str + r'\right ) ' + pd_str
else:
s += coef_str + ' ' + pd_str
s += ' + '

s = s.replace('+ -','- ')
return s[:-3]

def _repr_latex_(self):
latex_str = printer.GaLatexPrinter.latex(self)
return ' ' + latex_str + ' '

def __str__(self):
if printer.GaLatexPrinter.latex_flg:
Printer = printer.GaLatexPrinter
else:
Printer = printer.GaPrinter

return Printer().doprint(self)

def __repr__(self):
return str(self)

def __init__(self, *args, **kwargs):
kwargs = metric.test_init_slots(Sdop.init_slots, **kwargs)

self.Ga = kwargs['ga']  # Associated geometric algebra (coords)

if self.Ga is None:
raise ValueError('In Sdop.__init__ self.Ga must be defined.')

if len(args) == 1 and isinstance(args[0],Symbol):  # Simple Pdop of order 1
self.terms = ((S(1), Pdop(args[0], ga=self.Ga)),)
else:
if len(args) == 2 and isinstance(args[0],list) and isinstance(args[1],list):
if len(args[0]) != len(args[1]):
raise ValueError('In Sdop.__init__ coefficent list and Pdop list must be same length.')
self.terms = tuple(zip(args[0], args[1]))
elif len(args) == 1 and isinstance(args[0], (list, tuple)):
self.terms = tuple(args[0])
else:
raise ValueError('In Sdop.__init__ length of args must be 1 or 2 args = '+str(args))

def __call__(self, arg):
if isinstance(arg, Sdop):
if self.Ga != arg.Ga:
raise ValueError(
'When chaining scalar differential operators, self.Ga != arg.Ga.')
terms = []
for (coef, pdiff) in self.terms:
new_terms = pdiff(arg.terms)
new_terms = [(coef * c, p) for c, p in new_terms]
terms += new_terms
product = Sdop(terms, ga=self.Ga)
return Sdop.consolidate_coefs(product)
else:
return sum([coef * pdiff(arg) for coef, pdiff in self.terms], S(0))

def __neg__(self):
return Sdop([(-coef, pdiff) for coef, pdiff in self.terms], ga=self.Ga)

@staticmethod
if isinstance(sdop1, Sdop) and isinstance(sdop2, Sdop):
if sdop1.Ga != sdop2.Ga:
raise ValueError('In Sdop.Add sdop1.Ga != sdop2.Ga.')
return Sdop(_merge_terms(sdop1.terms, sdop2.terms), ga=sdop1.Ga)
else:
# convert values to multiplicative operators
if isinstance(sdop1, Sdop):
sdop2 = Sdop([(sdop2, Pdop({}, ga=sdop1.Ga))], ga=sdop1.Ga)
elif isinstance(sdop2, Sdop):
sdop1 = Sdop([(sdop1, Pdop({}, ga=sdop2.Ga))], ga=sdop2.Ga)
else:
raise TypeError("Neither argument is a Dop instance")

def __eq__(self, other):
if isinstance(other, Sdop):
if self.Ga != other.Ga:
return NotImplemented

diff = self - other
return len(diff.terms) == 0
else:
return NotImplemented

if sys.version_info.major < 3:
def __ne__(self, other):
return not (self == other)

def __sub__(self, sdop):

def __rsub__(self, sdop):

def __mul__(self, sdopr):
# alias for applying the operator
return self.__call__(sdopr)

def __rmul__(self, sdop):
terms = [(sdop * coef, pdiff) for coef, pdiff in self.terms]
return Sdop(terms, ga=self.Ga)

#################### Partial Derivative Operator Class #################

[docs]class Pdop(object):
r"""
Partial derivative class for multivectors.  The partial derivatives
are of the form

.. math::
\partial_{i_{1}...i_{n}} =
\frac{\partial^{i_{1}+...+i_{n}}}{\partial{x_{1}^{i_{1}}}...\partial{x_{n}^{i_{n}}}}.

If :math:i_{j} = 0 then the partial derivative does not contain the
:math:x^{i_{j}} coordinate.

Attributes
----------
pdiffs : dict
A dictionary where coordinates are keys and key value are the number of
times one differentiates with respect to the key.
order : int
Total number of differentiations.
When this is zero (i.e. when :attr:pdiffs is {}) then this object
is the identity operator, and returns its operand unchanged.
"""

init_slots = {'ga': (None, 'Associated geometric algebra')}

def sort_key(self, order=None):
return (
# lower order derivatives first
self.order,
# sorted by symbol after that, after expansion
sorted([
x.sort_key(order)
if x not in self.Ga.coords else
self.Ga.coords.index(x)
for x, k in self.pdiffs.items()
for i in range(k)
])
)

def __eq__(self,A):
if isinstance(A, Pdop) and self.Ga == A.Ga and self.pdiffs == A.pdiffs:
return True
else:
if len(self.pdiffs) == 0 and A == S(1):
return True
return False

if sys.version_info.major < 3:
def __ne__(self, other):
return not (self == other)

def __init__(self, __arg, **kwargs):
"""
The partial differential operator is a partial derivative with
respect to a set of real symbols (variables).  The allowed
variables are in two lists. self.Ga.coords is a list of the
coordinates associated with the geometric algebra.  self.Ga.auxvars
is a list of auxiallary symbols that have be added to the geometric
algebra using the member function Ga.AddVars(self,auxvars).

"""

kwargs = metric.test_init_slots(Pdop.init_slots, **kwargs)

self.Ga = kwargs['ga']  # Associated geometric algebra

if self.Ga is None:
raise ValueError('In Pdop.__init__ self.Ga must be defined.')

# galgebra 0.4.5
if __arg is None:
warnings.warn(
"Pdop(None) is deprecated, use Pdop({}) instead",
DeprecationWarning, stacklevel=2)
__arg = {}

if isinstance(__arg, dict):  # Pdop defined by dictionary
self.pdiffs = __arg
elif isinstance(__arg, Symbol):  # First order derivative with respect to symbol
self.pdiffs = {__arg: 1}
else:
raise TypeError('A dictionary or symbol is required, got {!r}'.format(__arg))

self.order = sum(self.pdiffs.values())

[docs]    def factor(self):
"""
If partial derivative operator self.order > 1 factor out first
order differential operator.  Needed for application of partial
derivative operator to product of sympy expression and partial
differential operator.  For example if D = Pdop({x:3}) then::

(Pdop({x:2}), Pdop({x:1})) = D.factor()
"""
if self.order == 1:
return S(0), self
else:
new_pdiffs = self.pdiffs.copy()
x, n = next(iter(new_pdiffs.items()))
if n == 1:
del new_pdiffs[x]
else:
new_pdiffs[x] -= 1
return Pdop(new_pdiffs, ga=self.Ga), Pdop(x, ga=self.Ga)

def __call__(self, arg):
"""
Calculate nth order partial derivative (order defined by
self) of :class:Mv, :class:Dop, :class:Sdop or sympy expression
"""
if self.pdiffs == {}:
return arg  # result is Pdop identity (1)

if isinstance(arg, Pdop):  # arg is Pdop
if self.Ga != arg.Ga:
raise ValueError('In Pdop.__call__ arguments do not belong to same geometric algebra.')
elif arg.pdiffs == {}:  # arg is one
return self
#return S(0)  # derivative is zero
else:  # arg is partial derivative
pdiffs = copy.copy(arg.pdiffs)
for key in self.pdiffs:
if key in pdiffs:
pdiffs[key] += self.pdiffs[key]
else:
pdiffs[key] = self.pdiffs[key]
return Pdop(pdiffs,ga=self.Ga)  # result is Pdop

elif isinstance(arg, Mv):  # arg is multivector
for x in self.pdiffs:
for i in range(self.pdiffs[x]):
arg = self.Ga.pDiff(arg, x)
return arg  # result is multivector

elif isinstance(arg, (Expr, Symbol, numbers.Number)):  # arg is sympy expression
for x in self.pdiffs:
arg = diff(arg,x,self.pdiffs[x])
return arg  # derivative is sympy expression

elif isinstance(arg, (list, tuple)):  # arg is list of tuples (coef, partial derivative)
terms = list(arg)
D = self
while True:
D, D0 = D.factor()
for k, term in enumerate(terms):
dc = D0(term[0])
pd = D0(term[1])
#print 'D0, term, dc, pd =', D0, term, dc, pd
tmp = []
if dc != 0:
tmp.append((dc,term[1]))
if pd != 0 :
tmp.append((term[0],pd))
terms[k] = tmp
terms = [i for o in terms for i in o]  # flatten list one level
if D == 0:
break
terms = Sdop.consolidate_coefs(terms)
return terms  # result is list of tuples (coef, partial derivative)
elif isinstance(arg, Sdop):  # arg is scalar differential operator
if self.Ga != arg.Ga:
raise ValueError('In Pdop.__call__ self.Ga != arg.Ga.')
return self(arg.terms)  # result is list of tuples (coef, partial derivative)
else:
raise ValueError('In Pdop.__call__ type(arg) = ' + str(type(arg)) + ' not allowed.')

def __mul__(self, pdop):  # functional product of self and arg (self*arg)
return self(pdop)

def __rmul__(self, pdop):  # functional product of arg and self (arg*self)
if isinstance(pdop, Pdop):
return pdop(self)
return Sdop([(pdop, self)], ga=self.Ga)

def Pdop_str(self):
if self.order == 0:
return 'D{}'
s = 'D'
for x in self.pdiffs:
s += '{' + str(x) + '}'
n = self.pdiffs[x]
if n > 1:
s += '^' + str(n)
return s

def Pdop_latex_str(self):
if self.order == 0:
return ''
s = r'\frac{\partial'
if self.order > 1:
s += '^{' + printer.latex(self.order) + '}'
s += '}{'
keys = list(self.pdiffs.keys())
keys.sort(key=(self.Ga.coords + keys).index)
for key in keys:
i = self.pdiffs[key]
s += r'\partial ' + printer.latex(key)
if i > 1:
s += '^{' + printer.latex(i) + '}'
s += '}'
return s

def _repr_latex_(self):
latex_str = printer.GaLatexPrinter.latex(self)
return ' ' + latex_str + ' '

def __str__(self):
if printer.GaLatexPrinter.latex_flg:
Printer = printer.GaLatexPrinter
else:
Printer = printer.GaPrinter
return Printer().doprint(self)

def __repr__(self):
return str(self)

################# Multivector Differential Operator Class ##############

[docs]class Dop(object):
r"""
Differential operator class for multivectors.  The operators are of
the form

.. math:: D = D^{i_{1}...i_{n}}\partial_{i_{1}...i_{n}}

where the :math:D^{i_{1}...i_{n}} are multivector functions of the coordinates
:math:x_{1},...,x_{n} and :math:\partial_{i_{1}...i_{n}} are partial derivative
operators

.. math:: \partial_{i_{1}...i_{n}} =
\frac{\partial^{i_{1}+...+i_{n}}}{\partial{x_{1}^{i_{1}}}...\partial{x_{n}^{i_{n}}}}.

If :math:* is any multivector multiplicative operation then the operator D
operates on the multivector function :math:F by the following definitions

.. math:: D*F = D^{i_{1}...i_{n}}*\partial_{i_{1}...i_{n}}F

returns a multivector and

.. math:: F*D = F*D^{i_{1}...i_{n}}\partial_{i_{1}...i_{n}}

returns a differential operator.  If the :attr:cmpflg in the operator is
set to True the operation returns

.. math:: F*D = (\partial_{i_{1}...i_{n}}F)*D^{i_{1}...i_{n}}

a multivector function.  For example the representation of the grad
operator in 3d would be:

.. math::
D^{i_{1}...i_{n}} &= [e_x,e_y,e_z] \\
\partial_{i_{1}...i_{n}} &= [(1,0,0),(0,1,0),(0,0,1)].

See LaTeX documentation for definitions of operator algebraic
operations +, -, *, ^, |, <, and >.

Attributes
----------
ga : ~galgebra.ga.Ga
Associated geometric algebra
cmpflg : bool
Complement flag
terms : list of tuples
"""

init_slots = {'ga': (None, 'Associated geometric algebra'),
'cmpflg': (False, 'Complement flag for Dop'),
'debug': (False, 'True to print out debugging information'),
'fmt_dop': (1, '1 for normal dop partial derivative formating')}

def __init__(self, *args, **kwargs):

kwargs = metric.test_init_slots(Dop.init_slots, **kwargs)

self.cmpflg = kwargs['cmpflg']  # Complement flag (default False)
self.Ga = kwargs['ga']

if self.Ga is None:
raise ValueError('In Dop.__init__ self.Ga must be defined.')

self.dop_fmt = kwargs['fmt_dop']  # Partial derivative output format (default 1)
self.title = None

if len(args) == 2:
coefs, pdiffs = args
if len(coefs) != len(pdiffs):
raise ValueError('In Dop.__init__ coefficent list and Pdop list must be same length.')
self.terms = tuple(zip(coefs, pdiffs))
elif len(args) == 1:
arg, = args
if len(arg) == 0:
self.terms = ()
elif isinstance(arg[0][0], Mv):  # Mv expansion [(Mv, Pdop)]
self.terms = tuple(arg)
elif isinstance(arg[0][0], Sdop):  # Sdop expansion [(Sdop, Mv)]
self.terms = _consolidate_terms(
(coef * mv, pdiff)
for (sdop, mv) in arg
for (coef, pdiff) in sdop.terms
)
else:
raise ValueError('In Dop.__init__ args[0] form not allowed. args = ' + str(args))
else:
raise ValueError('In Dop.__init__ length of args must be 1 or 2.')

[docs]    def simplify(self, modes=simplify):
"""
Simplify each multivector coefficient of a partial derivative
"""
return Dop(
[(coef.simplify(modes=modes), pd) for coef, pd in self.terms],
ga=self.Ga, cmpflg=self.cmpflg
)

[docs]    def consolidate_coefs(self):
"""
Remove zero coefs and consolidate coefs with repeated pdiffs.
"""
return Dop(_consolidate_terms(self.terms), ga=self.Ga, cmpflg=self.cmpflg)

@staticmethod

if isinstance(dop1, Dop) and isinstance(dop2, Dop):
if dop1.Ga != dop2.Ga:
raise ValueError('In Dop.Add Dop arguments are not from same geometric algebra')

if dop1.cmpflg != dop2.cmpflg:
raise ValueError('In Dop.Add complement flags have different values: %s vs. %s' % (dop1.cmpflg, dop2.cmpflg))

return Dop(_merge_terms(dop1.terms, dop2.terms), cmpflg=dop1.cmpflg, ga=dop1.Ga)
else:
# convert values to multiplicative operators
if isinstance(dop1, Dop):
if not isinstance(dop2, Mv):
dop2 = dop1.Ga.mv(dop2)
dop2 = Dop([(dop2, Pdop({}, ga=dop1.Ga))], cmpflg=dop1.cmpflg, ga=dop1.Ga)
elif isinstance(dop2, Dop):
if not isinstance(dop1, Mv):
dop1 = dop2.Ga.mv(dop1)
dop1 = Dop([(dop1, Pdop({}, ga=dop2.Ga))], cmpflg=dop2.cmpflg, ga=dop2.Ga)
else:
raise TypeError("Neither argument is a Dop instance")

def __neg__(self):
return Dop(
[(-coef, pdiff) for coef, pdiff in self.terms],
ga=self.Ga, cmpflg=self.cmpflg
)

def __sub__(self, dop):

def __rsub__(self, dop):

@staticmethod
def Mul(dopl, dopr, op='*'):  # General multiplication of Dop's
# cmpflg is True if the Dop operates on the left argument and
# False if the Dop operates on the right argument

if isinstance(dopl, Dop) and isinstance(dopr, Dop):
if dopl.Ga != dopr.Ga:
raise ValueError('In Dop.Mul Dop arguments are not from same geometric algebra')
ga = dopl.Ga
if dopl.cmpflg != dopr.cmpflg:
raise ValueError('In Dop.Mul Dop arguments do not have same cmplfg')
if not dopl.cmpflg:  # dopl and dopr operate on right argument
terms = []
for (coef, pdiff) in dopl.terms:  #Apply each dopl term to dopr
Ddopl = pdiff(dopr.terms)  # list of terms
Ddopl = [(Mv.Mul(coef, c, op=op), p) for c, p in Ddopl]
terms += Ddopl
product = Dop(terms, ga=ga)
else:  # dopl and dopr operate on left argument
terms = []
for (coef, pdiff) in dopr.terms:
Ddopr = pdiff(dopl.terms)  # list of terms
Ddopr = [(Mv.Mul(c, coef, op=op), p) for c, p in Ddopr]
terms += Ddopr
product = Dop(terms, ga=ga, cmpflg=True)
else:
if not isinstance(dopl, Dop):  # dopl is a scalar or Mv and dopr is Dop
if isinstance(dopl, Mv) and dopl.Ga != dopr.Ga:
raise ValueError('In Dop.Mul Dop arguments are not from same geometric algebra')
else:
dopl = dopr.Ga.mv(dopl)
ga = dopl.Ga

if not dopr.cmpflg:  # dopr operates on right argument
terms = [(Mv.Mul(dopl, coef, op=op), pdiff) for coef, pdiff in dopr.terms]
return Dop(terms, ga=ga)  # returns Dop
else:
product = sum([Mv.Mul(pdiff(dopl), coef, op=op) for coef, pdiff in dopr.terms], Mv(0, ga=ga))  # returns multivector
else:  # dopr is a scalar or a multivector

if isinstance(dopr, Mv) and dopl.Ga != dopr.Ga:
raise ValueError('In Dop.Mul Dop arguments are not from same geometric algebra')
ga = dopl.Ga

if not dopl.cmpflg:  # dopl operates on right argument
return sum([Mv.Mul(coef, pdiff(dopr), op=op) for coef, pdiff in dopl.terms], Mv(0, ga=ga))  # returns multivector
else:
terms = [(Mv.Mul(coef, dopr, op=op), pdiff) for coef, pdiff in dopl.terms]
product = Dop(terms, ga=dopl.Ga, cmpflg=True)  # returns Dop complement
if isinstance(product, Dop):
product = product.consolidate_coefs()
return product

def TSimplify(self):
return Dop([
(metric.Simp.apply(coef), pdiff) for (coef, pdiff) in self.terms
], ga=self.Ga)

def __truediv__(self, dopr):
if isinstance(dopr, (Dop, Mv)):
raise TypeError('In Dop.__truediv__ dopr must be a sympy scalar.')
return Dop([
(coef / dopr, pdiff) for (coef, pdiff) in self.terms
], ga=self.Ga, cmpflg=self.cmpflg)

if sys.version_info.major < 3:
__div__ = __truediv__

def __mul__(self, dopr):  # * geometric product
return Dop.Mul(self, dopr, op='*')

def __rmul__(self, dopl):  # * geometric product
return Dop.Mul(dopl, self, op='*')

def __xor__(self, dopr):  # ^ outer product
return Dop.Mul(self, dopr, op='^')

def __rxor__(self, dopl):  # ^ outer product
return Dop.Mul(dopl, self, op='^')

def __or__(self, dopr):  # | inner product
return Dop.Mul(self, dopr, op='|')

def __ror__(self, dopl):  # | inner product
return Dop.Mul(dopl, self, op='|')

def __lt__(self, dopr):  # < left contraction
return Dop.Mul(self, dopr, op='<')

def __gt__(self, dopr):  # > right contraction
return Dop.Mul(self, dopr, op='>')

def __eq__(self, other):
if isinstance(other, Dop):
if self.Ga != other.Ga:
return NotImplemented

diff = self - other
return len(diff.terms) == 0
else:
return NotImplemented

if sys.version_info.major < 3:
def __ne__(self, other):
return not (self == other)

def __str__(self):
if printer.GaLatexPrinter.latex_flg:
Printer = printer.GaLatexPrinter
else:
Printer = printer.GaPrinter

return Printer().doprint(self)

def __repr__(self):
return str(self)

def _repr_latex_(self):
latex_str = printer.GaLatexPrinter.latex(self)
if r'\begin{align*}' not in latex_str:
if self.title is None:
latex_str = r'\begin{equation*} ' + latex_str + r' \end{equation*}'
else:
latex_str = r'\begin{equation*} ' + self.title + ' = ' + latex_str + r' \end{equation*}'
else:
if self.title is not None:
latex_str = latex_str.replace('&',' ' + self.title + ' =&',1)
return latex_str

def is_scalar(self):
for coef, pdiff in self.terms:
if isinstance(coef, Mv) and not coef.is_scalar():
return False
return True

def components(self):
return tuple(
Dop(_consolidate_terms(
(Mv(coef * base, ga=self.Ga), pdiff)
for (coef, pdiff) in sdop.terms
), ga=self.Ga)
for (sdop, base) in self.Dop_mv_expand()
)

def Dop_mv_expand(self, modes=None):
coefs = []
bases = []
self.consolidate_coefs()

for (coef, pdiff) in self.terms:
if isinstance(coef, Mv) and not coef.is_scalar():
for mv_coef, mv_base in metric.linear_expand_terms(coef.obj):
if mv_base in bases:
index = bases.index(mv_base)
coefs[index] += Sdop([(mv_coef, pdiff)], ga=self.Ga)
else:
bases.append(mv_base)
coefs.append(Sdop([(mv_coef, pdiff)], ga=self.Ga))
else:
if isinstance(coef, Mv):
mv_coef = coef.obj
else:
mv_coef = coef
if S(1) in bases:
index = bases.index(S(1))
coefs[index] += Sdop([(mv_coef, pdiff)], ga=self.Ga)
else:
bases.append(S(1))
coefs.append(Sdop([(mv_coef, pdiff)], ga=self.Ga))
if modes is not None:
for i in range(len(coefs)):
coefs[i] = coefs[i].simplify(modes)
terms = list(zip(coefs, bases))

def Dop_str(self):
if len(self.terms) == 0:
return ZERO_STR

mv_terms = self.Dop_mv_expand(modes=simplify)
s = ''

for (sdop, base) in mv_terms:
str_base = printer.latex(base)
str_sdop = printer.latex(sdop)
if base == S(1):
s += str_sdop
else:
if len(sdop.terms) > 1:
if self.cmpflg:
s += '(' + str_sdop + ')*' + str_base
else:
s += str_base + '*(' + str_sdop + ')'
else:
if str_sdop[0] == '-' and not isinstance(sdop.terms[0][0], Add):
if self.cmpflg:
s += str_sdop + '*' + str_base
else:
s += '-' + str_base + '*' + str_sdop[1:]
else:
if self.cmpflg:
s += str_sdop + '*' + str_base
else:
s += str_base + '*' + str_sdop
s += ' + '

s = s.replace('+ -','-')
return s[:-3]

def Dop_latex_str(self):
if len(self.terms) == 0:
return ZERO_STR

self.consolidate_coefs()

mv_terms = self.Dop_mv_expand(modes=simplify)
s = ''

for (sdop, base) in mv_terms:
str_base = printer.latex(base)
str_sdop = printer.latex(sdop)
if base == S(1):
s += str_sdop
else:
if str_sdop == '1':
s += str_base
if str_sdop == '-1':
s += '-' + str_base
if str_sdop[1:] != '1':
s += ' ' + str_sdop[1:]
else:
if len(sdop.terms) > 1:
if self.cmpflg:
s += r'\left ( ' + str_sdop + r'\right ) ' + str_base
else:
s += str_base + ' ' + r'\left ( ' + str_sdop + r'\right ) '
else:
if str_sdop[0] == '-' and not isinstance(sdop.terms[0][0], Add):
if self.cmpflg:
s += str_sdop + str_base
else:
s += '-' + str_base + ' ' + str_sdop[1:]
else:
if self.cmpflg:
s += str_sdop + ' ' + str_base
else:
s += str_base + ' ' + str_sdop
s += ' + '

s = s.replace('+ -','-')
Sdop.str_mode = False
return s[:-3]

def Fmt(self, fmt=1, title=None, dop_fmt=None):
if printer.GaLatexPrinter.latex_flg:
printer.GaLatexPrinter.prev_fmt = printer.GaLatexPrinter.fmt
printer.GaLatexPrinter.prev_dop_fmt = printer.GaLatexPrinter.dop_fmt
else:
printer.GaPrinter.prev_fmt = printer.GaPrinter.fmt
printer.GaPrinter.prev_dop_fmt = printer.GaPrinter.dop_fmt

if title is not None:
self.title = title

if printer.isinteractive():
return self

if Mv.latex_flg:
latex_str = printer.GaLatexPrinter.latex(self)
printer.GaLatexPrinter.fmt = printer.GaLatexPrinter.prev_fmt
printer.GaLatexPrinter.dop_fmt = printer.GaLatexPrinter.prev_dop_fmt

if title is not None:
return title + ' = ' + latex_str
else:
return latex_str
else:
s = str(self)
printer.GaPrinter.fmt = printer.GaPrinter.prev_fmt
printer.GaPrinter.dop_fmt = printer.GaPrinter.prev_dop_fmt

if title is not None:
return title + ' = ' + s
else:
return s

def Nga(x, prec=5):
if isinstance(x, Mv):
return Mv(Nsympy(x.obj, prec), ga=x.Ga)
else:
return Nsympy(x, prec)

def printeigen(M):    # Print eigenvalues, multiplicities, eigenvectors of M.
evects = M.eigenvects()
for i in range(len(evects)):                   # i iterates over eigenvalues
print(('Eigenvalue =', evects[i][0], '  Multiplicity =', evects[i][1], ' Eigenvectors:'))
for j in range(len(evects[i][2])):         # j iterates over eigenvectors of a given eigenvalue
result = '['
for k in range(len(evects[i][2][j])):  # k iterates over coordinates of an eigenvector
result += str(trigsimp(evects[i][2][j][k]).evalf(3))
if k != len(evects[i][2][j]) - 1:
result += ', '
result += '] '
print(result)

def printGS(M, norm=False):  # Print Gram-Schmidt output.
from sympy import GramSchmidt
global N
N = GramSchmidt(M, norm)
result = '[ '
for i in range(len(N)):
result += '['
for j in range(len(N[0])):
result += str(trigsimp(N[i][j]).evalf(3))
if j != len(N[0]) - 1:
result += ', '
result += '] '
if j != len(N[0]) - 1:
result += ' '
result += ']'
print(result)

def printrref(matrix, vars="xyzuvwrs"):   # Print rref of matrix with variables.
rrefmatrix = matrix.rref()[0]
rows, cols = rrefmatrix.shape
if len(vars) < cols - 1:
print('Not enough variables.')
return
for i in range(rows):
result = ''
for j in range(cols - 1):
result += str(rrefmatrix[i, j]) + vars[j]
if j != cols - 2:
result += ' + '
result += ' = ' + str(rrefmatrix[i, cols - 1])
print(result)

def com(A, B):
raise ImportError(

def correlation(u, v, dec=3):  # Compute the correlation coefficient of vectors u and v.
rows, cols = u.shape
uave = 0
vave = 0
for i in range(rows):
uave += u[i]
vave += v[i]
uave = uave / rows
vave = vave / rows
ulocal = u[:, :]  # Matrix copy
vlocal = v[:, :]
for i in range(rows):
ulocal[i] -= uave
vlocal[i] -= vave
return ulocal.dot(vlocal) / (ulocal.norm() * vlocal.norm()). evalf(dec)

def cross(v1, v2):
if v1.is_vector() and v2.is_vector() and v1.Ga == v2.Ga and v1.Ga.n == 3:
return -v1.Ga.I() * (v1 ^ v2)
else:
raise ValueError(str(v1) + ' and ' + str(v2) + ' not compatible for cross product.')

def dual(A):
if isinstance(A, Mv):
return A.dual()
else:
raise ValueError('A not a multivector in dual(A)')

def even(A):
if not isinstance(A,Mv):
raise ValueError('A = ' + str(A) + ' not a multivector in even(A).')
return A.even()

def odd(A):
if not isinstance(A,Mv):
raise ValueError('A = ' + str(A) + ' not a multivector in even(A).')
return A.odd()

def exp(A,hint='-'):
if isinstance(A,Mv):
return A.exp(hint)
else:
return sympy_exp(A)

if isinstance(A, Mv):
else:
raise ValueError('A not a multivector in grade(A,r)')

def inv(A):
if not isinstance(A,Mv):
raise ValueError('A = ' + str(A) + ' not a multivector in inv(A).')
return A.inv()

def norm(A, hint='+'):
if isinstance(A, Mv):
return A.norm(hint=hint)
else:
raise ValueError('A not a multivector in norm(A)')

def norm2(A):
if isinstance(A, Mv):
return A.norm2()
else:
raise ValueError('A not a multivector in norm(A)')

def proj(B, A):  # Project on the blade B the multivector A
if isinstance(A,Mv):
else:
raise ValueError('A not a multivector in proj(B,A)')

def rot(itheta, A, hint='-'):  # Rotate by the 2-blade itheta the multivector A
if isinstance(A,Mv):
return A.rotate_multivector(itheta, hint)
else:
raise ValueError('A not a multivector in rotate(A,itheta)')

def refl(B, A):  #  Project on the blade B the multivector A
if isinstance(A,Mv):
else:
raise ValueError('A not a multivector in reflect(B,A)')

def rev(A):
if isinstance(A, Mv):
return A.rev()
else:
raise ValueError('A not a multivector in rev(A)')

def scalar(A):
if not isinstance(A,Mv):
raise ValueError('A = ' + str(A) + ' not a multivector in inv(A).')
return A.scalar()