# Source code for galgebra.lt

"""
Multivector Linear Transformation
"""

import inspect
import types
import itertools
import warnings
from copy import copy
from functools import reduce
from typing import Mapping

from sympy import (
expand, symbols, Matrix, Transpose, zeros, Symbol, Function, S, Add, Expr
)
from sympy.printing.latex import LatexPrinter as _LatexPrinter
from sympy.printing.str import StrPrinter as _StrPrinter

from ._utils import cached_property as _cached_property

from . import printer
from . import metric
from . import mv

# Add custom settings to the builtin latex printer
_LatexPrinter._default_settings.update({
'galgebra_mlt_lcnt': 1
})
_StrPrinter._default_settings.update({
'galgebra_mlt_lcnt': 1
})

def Symbolic_Matrix(root, coords=None, mode='g', f=False, sub=True):
if sub:
pos = '_'
else:
pos = '__'
if isinstance(coords, (list, tuple)):
n = len(coords)
n_range = range(n)
mat = zeros(n)
if mode == 'g':  # General symbolic matrix
for row in n_range:
row_index = str(coords[row])
for col in n_range:
col_index = str(coords[col])
element = root + pos + row_index + col_index
if not f:
mat[row, col] = Symbol(element, real=True)
else:
mat[row, col] = Function(element)(*coords)

elif mode == 's':  # Symmetric symbolic matrix
for row in n_range:
row_index = str(coords[row])
for col in n_range:
col_index = str(coords[col])
if row <= col:
element = root + pos + row_index + col_index
else:
element = root + pos + col_index + row_index
if not f:
mat[row, col] = Symbol(element, real=True)
else:
mat[row, col] = Function(element)(*coords)

elif mode == 'a':  # Asymmetric symbolic matrix
for row in n_range:
row_index = str(coords[row])
for col in n_range:
col_index = str(coords[col])
if row <= col:
sign = S.One
element = root + pos + row_index + col_index
else:
sign = -S.One
element = root + pos + col_index + row_index
if row == col:
sign = S.Zero
if not f:
mat[row, col] = sign * Symbol(element, real=True)
else:
mat[row, col] = sign * Function(element)(*coords)
else:
raise ValueError('In Symbolic_Matrix mode = ' + str(mode))
else:
raise ValueError('In Symbolic_Matrix coords = ' + str(coords))
return mat

[docs]def Matrix_to_dictionary(mat_rep, basis):
""" Convert matrix representation of linear transformation to dictionary """
n = len(basis)
if mat_rep.rows != n or mat_rep.cols != n:
raise ValueError('Matrix and Basis dimensions not equal for Matrix = ' + str(mat_rep))
n_range = list(range(n))
return {
basis[col]: sum(
(mat_rep[row, col]*basis[row] for row in n_range), S.Zero
)
for col in n_range
}

[docs]def Dictionary_to_Matrix(dict_rep, ga):
""" Convert dictionary representation of linear transformation to matrix """
lst_mat = []  # list representation of sympy matrix
for e_row in ga.basis:
lst_mat_row = len(ga.basis) * [S.Zero]

element = dict_rep.get(e_row, S.Zero)
if isinstance(element, mv.Mv):
element = element.obj
for coef, base in metric.linear_expand_terms(element):
index = ga.basis.index(base)
lst_mat_row[index] = coef

lst_mat.append(lst_mat_row)
# expand the transpose
return Transpose(Matrix(lst_mat)).doit()

[docs]class Lt(printer.GaPrintable):
r"""
A Linear Transformation

Except for the spinor representation the linear transformation
is stored as a dictionary with basis vector keys and vector
values self.lt_dict so that a is a vector :math:a = a^{i}e_{i} then

.. math::
\mathtt{self(}a\mathtt{)}
= a^{i} * \mathtt{self.lt\_dict[}e_{i}\mathtt{]}.

For the spinor representation the linear transformation is
stored as the even multivector self.R so that if a is a
vector::

self(a) = self.R * a * self.R.rev().

Attributes
----------
lt_dict : dict
the keys are the basis symbols, :math:e_i, and the dictionary
entries are the object vector images (linear combination of sympy
non-commutative basis symbols) of the keys so that if L is the
linear transformation then::

L(e_i) = self.Ga.mv(L.lt_dict[e_i])

"""

@staticmethod
def setup(ga):
# galgebra 0.5.0
warnings.warn(
"Lt.setup(ga) is deprecated, use ga.coords and ga.coord_vec "
"directly.", DeprecationWarning, stacklevel=2)
return ga.coords, ga.coord_vec

@property
def coords(self):
# galgebra 0.6.0
warnings.warn(
"lt.coords is deprecated, use lt.Ga.coords instead.",
DeprecationWarning, stacklevel=2)
return self.Ga.coords

@property
def X(self):
# galgebra 0.6.0
warnings.warn(
"lt.X is deprecated, use lt.Ga.coord_vec instead.",
DeprecationWarning, stacklevel=2)
return self.Ga.coord_vec

@property
def mode(self):
# galgebra 0.6.0
warnings.warn(
"lt.mode is deprecated, inspect lt.matrix() and its transpose to "
"determine symmetry",
DeprecationWarning, stacklevel=2)
m = self.matrix()
if m == m.T:
return 's'
elif m == -m.T:
return 'a'
else:
return 'g'

@property
def fct_flg(self):
# galgebra 0.6.0
warnings.warn(
"lt.fct_flg is deprecated, inspect lt.matrix().free_symbols to "
"determine coordinate-dependence",
DeprecationWarning, stacklevel=2)
if self.Ga.coords is None:
return False
return set(self.Ga.coords) <= self.matrix().free_symbols

def __init__(self, *args, ga, f=False, mode='g'):
"""
__init__(self, *args, ga, **kwargs)

Note this constructor is overloaded, based on the type of the
positional argument:

.. class:: Lt(lt_dict: Dict[Expr, Expr], /, *, ga)
:noindex:

Construct from a dictionary mapping source basis blade expressions
to multivectors.

.. class:: Lt(lt_matrix: Matrix, /, *, ga)
:noindex:

Construct from the operation of matrix pre-multiplication.

.. class:: Lt(spinor: mv.Mv, /, *, ga)
:noindex:

Construct from a spinor / rotor, which need not square to one.

.. class:: Lt(func: Callable[[mv.Mv], mv.Mv], /, *, ga)
:noindex:

Construct from a function, which is tested for linearity.

.. class:: Lt(s: str, /, *, ga, f=False, mode='g')
:noindex:

Construct an appropriate matrix from a string s.

Parameters
----------
ga : Ga
Geometric algebra which is the domain and codomain of this transform
f : bool
True if Lt if function of coordinates. Only supported in the string
constructor
mode : str
g:general, s:symmetric, a:antisymmetric transformation.
Only supported in the string constructor.
"""
mat_rep = args[0]
self.Ga = ga
self.spinor = False
self.rho_sq = None

self.lt_dict = {}
self.mat = None

if isinstance(mat_rep, dict):  # Dictionary input
for key in mat_rep:
self.lt_dict[key] = mat_rep[key]

elif isinstance(mat_rep, list):  # List of lists input
if not isinstance(mat_rep[0], list):
for lt_i, base in zip(mat_rep, self.Ga.basis):
self.lt_dict[base] = lt_i
else:
# mat_rep = map(list, zip(*mat_rep))  # Transpose list of lists
for row, base1 in zip(mat_rep, self.Ga.basis):
tmp = 0
for col, base2 in zip(row, self.Ga.basis):
tmp += col * base2
self.lt_dict[base1] = tmp

elif isinstance(mat_rep, Matrix):  # Matrix input
self.mat = mat_rep
mat_rep = self.mat * self.Ga.g_inv
self.lt_dict = Matrix_to_dictionary(mat_rep, self.Ga.basis)

elif isinstance(mat_rep, mv.Mv):  # Spinor input
self.spinor = True
self.R = mat_rep
self.Rrev = mat_rep.rev()
self.rho_sq = self.R * self.Rrev
if self.rho_sq.is_scalar():
self.rho_sq = self.rho_sq.scalar()
if self.rho_sq == S.One:
self.rho_sq = None
else:
raise ValueError('In Spinor input for Lt, S*S.rev() not a scalar!\n')

elif isinstance(mat_rep, str):  # String input
Amat = Symbolic_Matrix(mat_rep, coords=self.Ga.coords, mode=mode, f=f)
self.__init__(Amat, ga=self.Ga)

elif callable(mat_rep):  # Linear multivector function input
# F is a multivector function to be tested for linearity
F = mat_rep
a = mv.Mv('a', 'vector', ga=self.Ga)
b = mv.Mv('b', 'vector', ga=self.Ga)
if F(a + b) != F(a) + F(b):
raise ValueError('{} is not linear'.format(F))
self.lt_dict = {}
for base in self.Ga.basis:
out = F(mv.Mv(base, ga=self.Ga))
if not out.is_vector():
raise ValueError('{} must return vectors'.format(F))
self.lt_dict[base] = out.obj
else:
raise TypeError("Unsupported argument type {}".format(type(mat_rep)))

@_cached_property
def mv_dict(self) -> Mapping[Expr, Expr]:
# dict for linear transformation of multivector
if self.spinor:
# no lt_dict
return None

return {
self.Ga.wedge,
(self.Ga.basis[i].xreplace(self.lt_dict) for i in index),
S.One
)
}

[docs]    def __call__(self, v, obj=False):
r"""
Returns the image of the multivector :math:A under the linear transformation :math:L.

:math:{{L}\lp {A} \rp } is defined by the linearity of :math:L, the
vector values :math:{{L}\lp {{{\eb}}_{i}} \rp }, and the definition
:math:{{L}\lp {{{\eb}}_{i_{1}}{\wedge}\dots{\wedge}{{\eb}}_{i_{r}}} \rp } = {{L}\lp {{{\eb}}_{i_{1}}} \rp }{\wedge}\dots{\wedge}{{L}\lp {{{\eb}}_{i_{r}}} \rp }.
"""

if isinstance(v, mv.Mv) and self.Ga != v.Ga:
raise ValueError('In Lt call Lt and argument refer to different vector spaces')

if self.spinor:
if not isinstance(v, mv.Mv):
v = mv.Mv(v, ga=self.Ga)
if self.rho_sq is None:
R_v_Rrev = self.R * v * self.Rrev
else:
R_v_Rrev = self.rho_sq * self.R * v * self.Rrev
if obj:
return R_v_Rrev.obj
else:
return R_v_Rrev

if isinstance(v, mv.Mv):
if v.is_vector():
lt_v = v.obj.xreplace(self.lt_dict)
if obj:
return lt_v
else:
return mv.Mv(lt_v, ga=self.Ga)
else:
mv_obj = v.obj
else:
mv_obj = mv.Mv(v, ga=self.Ga).obj

lt_v = mv_obj.xreplace(self.mv_dict)
if obj:
return lt_v
else:
return mv.Mv(lt_v, ga=self.Ga)

if self.Ga != LT.Ga:
raise ValueError("Attempting addition of Lt's from different geometric algebras")

for key in list(LT.lt_dict.keys()):
if key in self_add_LT:
else:

def __sub__(self, LT):

if self.Ga != LT.Ga:
raise ValueError("Attempting subtraction of Lt's from different geometric algebras")

for key in list(LT.lt_dict.keys()):
if key in self_add_LT:
else:

def __mul__(self, LT):

if isinstance(LT, Lt):

if self.Ga != LT.Ga:
raise ValueError("Attempting multiplication of Lt's from different geometric algebras")
self_mul_LT = {}
for base in LT.lt_dict:
self_mul_LT[base] = self(LT(base, obj=True), obj=True)
for key in self_mul_LT:
self_mul_LT[key] = metric.collect(expand(self_mul_LT[key]), self.Ga.basis)
return Lt(self_mul_LT, ga=self.Ga)
else:
self_mul_LT = {}
for key in self.lt_dict:
self_mul_LT[key] = LT * self.lt_dict[key]
return Lt(self_mul_LT, ga=self.Ga)

def __rmul__(self, LT):

if not isinstance(LT, Lt):
self_mul_LT = {}
for key in self.lt_dict:
self_mul_LT[key] = LT * self.lt_dict[key]
return Lt(self_mul_LT, ga=self.Ga)
else:
raise TypeError('Cannot have LT as left argument in Lt __rmul__\n')

[docs]    def det(self) -> Expr:  # det(L) defined by L(I) = det(L)I
r"""
Returns the determinant (a scalar) of the linear transformation,
:math:L, defined by :math:{{\det}\lp {L} \rp }I = {{L}\lp {I} \rp }.
"""

lt_I = self(self.Ga.i, obj=True)
det_lt_I = lt_I.subs(self.Ga.i.obj, S.One)
return det_lt_I

[docs]    def tr(self) -> Expr:  # tr(L) defined by tr(L) = grad|L(x)
r"""
Returns the trace (a scalar) of the linear transformation,
:math:L, defined by :math:{{\operatorname{tr}}\lp {L} \rp }=\nabla_{a}\cdot{{L}\lp {a} \rp }
where :math:a is a vector in the tangent space.
"""

connect_flg = self.Ga.connect_flg
self.Ga.connect_flg = False

F_x = mv.Mv(self(self.Ga.coord_vec, obj=True), ga=self.Ga)
tr_F = (self.Ga.grad | F_x).scalar()
self.Ga.connect_flg = connect_flg
return tr_F

[docs]    def adj(self) -> 'Lt':
r"""
Returns the adjoint (a linear transformation) of the linear
transformation, :math:L, defined by :math:a\cdot{{L}\lp {b} \rp } = b\cdot{{\bar{L}}\lp {a} \rp }
where :math:a and :math:b are any two vectors in the tangent space
and :math:\bar{L} is the adjoint of :math:L.
"""

for e_j in self.Ga.basis:
s = S.Zero
for e_i, er_i in zip(self.Ga.basis, self.Ga.r_basis):
s += er_i * self.Ga.hestenes_dot(e_j, self(e_i, obj=True))
if self.Ga.is_ortho:
else:

def inv(self):
if self.spinor:
Lt_inv = Lt(self.Rrev, ga=self.Ga)
Lt_inv.rho_sq = S.One/(self.rho_sq**2)
else:
raise ValueError('Lt inverse currently implemented only for spinor!\n')
return Lt_inv

def _sympystr(self, print_obj):

if self.spinor:
return 'R = ' + print_obj._print(self.R)
else:
pre = 'Lt('
s = ''
for base in self.Ga.basis:
if base in self.lt_dict:
s += pre + print_obj._print(base) + ') = ' + print_obj._print(mv.Mv(self.lt_dict[base], ga=self.Ga)) + '\n'
else:
s += pre + print_obj._print(base) + ') = 0\n'
return s[:-1]

def _latex(self, print_obj):
parts = []
for base in self.Ga.basis:
if self.spinor:
val = self.R * mv.Mv(base, ga=self.Ga) * self.Rrev
else:
val = mv.Mv(self.lt_dict.get(base, S.Zero), ga=self.Ga)
parts.append(print_obj._print(base) + ' &\\mapsto ' + print_obj._print(val))
return '\\left\\{ \\begin{aligned} ' + ' \\\\ '.join(parts) + ' \\end{aligned} \\right\\}'

def Fmt(self, fmt=1, title=None) -> printer.GaPrintable:
return printer._FmtResult(self, title)

[docs]    def matrix(self) -> Matrix:
r"""
Returns the matrix representation of the linear transformation,
:math:L, defined by :math:{{L}\lp {{{\eb}}_{i}} \rp } = L_{ij}{{\eb}}_{j}
where :math:L_{ij} is the matrix representation.
"""
if self.mat is not None:
return self.mat
else:
if self.spinor:
self.lt_dict = {}
for base in self.Ga.basis:
self.lt_dict[base] = self(base).simplify()
self.spinor = False
mat = self.matrix()
self.spinor = True
return mat
else:
"""
mat_rep = []
for base in self.Ga.basis:
if base in self.lt_dict:
row = []
image = (self.lt_dict[base])
if isinstance(image, mv.Mv):
image = image.obj
coefs, bases = metric.linear_expand(image)
for base in self.Ga.basis:
try:
i = bases.index(base)
row.append(coefs[i])
except:
row.append(0)
mat_rep.append(row)
else:
mat_rep.append(self.Ga.n * [0])
return Matrix(mat_rep).transpose()
"""
self.mat = Dictionary_to_Matrix(self.lt_dict, self.Ga) * self.Ga.g
return self.mat

[docs]class Mlt(printer.GaPrintable):
r"""
A multilinear transformation (mlt) is a multilinear multivector function of
a list of vectors (*args) :math:F(v_1,...,v_r) where for any argument slot
:math:j we have (:math:a is a scalar and :math:u_j a vector)

.. math::
F(v_1,...,a*v_j,...,v_r) &= a*F(v_1,...,v_j,...,v_r) \\
F(v_1,...,v_j+u_j,...,v_r) &= F(v_1,...,v_j,...,v_r) + F(v_1,...,u_j,...,v_r).

If F and G are two :class:Mlt\ s with the same number of argument slots then the sum is

.. math:: (F+G)F(v_1,...,v_r) = F(v_1,...,v_r) + G(v_1,...,v_r).

If :math:F and :math:G are two :class:Mlt\ s with :math:r and :math:s
argument slots then their product is

.. math:: (F*G)(v_1,...,v_r,...,v_{r+s}) = F(v_1,...,v_r)*G(v_{r+1},...,v_{r+s}),

where :math:* is any of the multivector multiplicative operations.
The derivative of a :class:Mlt with is defined as the directional derivative with respect
to the coordinate vector (we assume :math:F is implicitely a function of the
coordinates)

.. math:: F(v_1,...,v_r;v_{r+1}) = (v_{r+1} \bullet \nabla)F(v_1,...,v_j,...,v_r).

The contraction of a :class:Mlt between slots :math:j and :math:k is defined as the
geometric derivative of :math:F with respect to slot :math:k and the inner geometric
derivative with respect to slot :math:j (this gives the standard tensor
definition of contraction for the case that :math:F is a scalar function)

.. math::

\operatorname{Contract}(i,j,F)
&= \nabla_i \bullet (\nabla_j F(v_1,...,v_i,...,v_j,...,v_r)) \\
&= \nabla_j \bullet (\nabla_i F(v_1,...,v_i,...,v_j,...,v_r)).

This returns a :class:Mlt\ with slot :math:i and :math:j removed.
"""

@staticmethod
def subs(Ga, anew):
#  Generate coefficient substitution list for new Mlt slot
#  vectors (arguments) where anew is a list of slot vectors
#  to be substituted for the old slot vectors.
#  This is used when one wishes to substitute specific vector
#  values into the Mlt such as the basis/reciprocal basis vectors.
sub_lst = []
for i, a in enumerate(anew):
acoefs = a.get_coefs(1)
sub_lst += list(zip(Ga._mlt_pdiffs[i], acoefs))
return sub_lst

@staticmethod
def increment_slots(nargs, Ga):
# Increment cache of available slots (vector variables) if needed for Mlt class
n_a = len(Ga._mlt_a)
if n_a < nargs:
for i in range(n_a, nargs):
#  New slot variable with coefficients a_{n_a}__k
a = Ga.mv('a_' + str(i + 1), 'vector')
#  Append new slot variable a_j
Ga._mlt_a.append(a)
#  Append slot variable coefficients a_j__k for purpose
#  of differentiation
coefs = a.get_coefs(1)
Ga._mlt_pdiffs.append(coefs)
Ga._mlt_acoefs += coefs

@staticmethod
def extact_basis_indexes(Ga):
# galgebra 0.5.0
warnings.warn(
"Mlt.extact_basis_indexes(ga) is deprecated, use ga.basis_super_scripts",
DeprecationWarning, stacklevel=2)
return Ga.basis_super_scripts

def _sympystr(self, print_obj):
return print_obj._print(self.fvalue)

def _latex(self, print_obj):
if self.nargs <= 1:
return print_obj._print(self.fvalue)
expr_lst = Mlt.expand_expr(self.fvalue, self.Ga)
latex_str = '\\begin{aligned} '
first = True
lcnt = print_obj._settings['galgebra_mlt_lcnt']
cnt = 1  # Component count on line
for term in expr_lst:
coef_str = str(term[0])
coef_latex = print_obj._print(term[0])
coef_latex = r'\left ( ' + coef_latex + r'\right ) '
if first:
first = False
else:
if coef_str[0].strip() != '-' or term_add_flg:
coef_latex = ' + ' + coef_latex
for aij in term[1]:
coef_latex += print_obj._print(aij) + ' '
if cnt == 1:
latex_str += ' & ' + coef_latex
else:
latex_str += coef_latex
if cnt % lcnt == 0:
latex_str += '\\\\ '
cnt = 1
else:
cnt += 1
if lcnt == len(expr_lst) or lcnt == 1:
latex_str = latex_str[:-3]
latex_str = latex_str + ' \\end{aligned} '
return latex_str

[docs]    def Fmt(self, lcnt=1, title=None) -> printer.GaPrintable:
"""
Set format for printing of Tensors

Parameters
----------
lcnt :
Number of components per line

Notes
-----
Usage for tensor T example is::

T.fmt('2', 'T')

output is::

print 'T = '+str(A)

with two components per line.  Works for both standard printing and
for latex.
"""
obj = printer._WithSettings(self, dict(galgebra_mlt_lcnt=lcnt))
return printer._FmtResult(obj, title)

@staticmethod
def expand_expr(expr, ga):
lst_expr = []
expr = expand(expr)
for term in expr.args:
coef = S.One
a_lst = []
for factor in term.args:
if factor in ga._mlt_acoefs:
a_lst.append(factor)
else:
coef *= factor
a_lst = tuple([x for x in a_lst if x in ga._mlt_acoefs])
b_lst = tuple([ga._mlt_acoefs.index(x) for x in a_lst])
lst_expr.append((coef, a_lst, b_lst))
lst_expr = sorted(lst_expr, key=lambda x: x[2])
new_lst_expr = []
previous = (-1,)
first = True
a = None
for term in lst_expr:
if previous == term[2]:
coef += term[0]
previous = term[2]
else:
if not first:
new_lst_expr.append((coef, a))
else:
first = False
coef = term[0]
previous = term[2]
a = term[1]
new_lst_expr.append((coef, a))
return new_lst_expr

def __init__(self, f, Ga, nargs=None, fct=False):
#  f is a function, a multivector, a string, or a component expression
#  self.f is a function or None such as T | a_1 where T and a_1 are vectors
#  self.fvalue is a component expression such as
#  T_x*a_1__x+T_y*a_1__y+T_z*a_1__z for a rank 1 tensor in 3 space and all
#  symbols are sympy real scalar symbols
self.Ga = Ga
if isinstance(f, mv.Mv):
if f.is_vector():  # f is vector T = f | a1
self.nargs = 1
Mlt.increment_slots(1, Ga)
self.fvalue = (f | Ga._mlt_a[0]).obj
self.f = None
else:  # To be inplemented for f a general pure grade mulitvector
self.nargs = nargs
self.fvalue = f
self.f = None
elif isinstance(f, Lt):  # f is linear transformation T = a1 | f(a2)
self.nargs = 2
Mlt.increment_slots(2, Ga)
self.fvalue = (Ga._mlt_a[0] | f(Ga._mlt_a[1])).obj
self.f = None
elif isinstance(f, str) and nargs is not None:
self.f = None
self.nargs = nargs
Mlt.increment_slots(nargs, Ga)
self.fvalue = S.Zero
for t_index, a_prod in zip(itertools.product(self.Ga.basis_super_scripts, repeat=self.nargs),
itertools.product(*self.Ga._mlt_pdiffs)):
name = '{}_{}'.format(f, ''.join(map(str, t_index)))
if fct:  # Tensor field
coef = Function(name, real=True)(*self.Ga.coords)
else:  # Constant Tensor
coef = symbols(name, real=True)
self.fvalue += reduce(lambda x, y: x*y, a_prod, coef)

else:
if isinstance(f, types.FunctionType):  # Tensor defined by general multi-linear function
args, _varargs, _kwargs, _defaults = inspect.getargspec(f)
self.nargs = len(args)
self.f = f
Mlt.increment_slots(self.nargs, Ga)
self.fvalue = f(*tuple(Ga._mlt_a[0:self.nargs]))
else:  # Tensor defined by component expression
self.f = None
self.nargs = len(args)
Mlt.increment_slots(self.nargs, Ga)
self.fvalue = f

[docs]    def __call__(self, *args):
"""
Evaluate the multilinear function for the given vector arguments.
Note that a sympy scalar is returned, *not* a multilinear function.
"""
if len(args) == 0:
return self.fvalue
if self.f is not None:
return self.f(*args)
else:
sub_lst = []
for x, ai in zip(args, self.Ga._mlt_pdiffs):
for r_base, aij in zip(self.Ga.r_basis_mv, ai):
sub_lst.append((aij, (r_base | x).scalar()))
return self.fvalue.subs(sub_lst, simultaneous=True)

if isinstance(X, Mlt):
if self.nargs == X.nargs:
return Mlt(self.fvalue + X.fvalue, self.Ga, self.nargs)
else:
raise ValueError('In Mlt add number of args not the same\n')
else:
raise TypeError('In Mlt add second argument not an Mkt\n')

def __sub__(self, X):
if isinstance(X, Mlt):
if self.nargs == X.nargs:
return Mlt(self.fvalue - X.fvalue, self.Ga, self.nargs)
else:
raise ValueError('In Mlt sub number of args not the same\n')
else:
raise TypeError('In Mlt sub second argument not an Mlt\n')

def __mul__(self, X):
if isinstance(X, Mlt):
nargs = self.nargs + X.nargs
Mlt.increment_slots(nargs, self.Ga)
self_args = self.Ga._mlt_a[:self.nargs]
X_args = X.Ga._mlt_a[self.nargs:nargs]
value = (self(*self_args) * X(*X_args)).expand()
return Mlt(value, self.Ga, nargs)
else:
return Mlt(X * self.fvalue, self.Ga, self.nargs)

def __xor__(self, X):
if isinstance(X, Mlt):
nargs = self.nargs + X.nargs
Mlt.increment_slots(nargs, self.Ga)
value = self(*self.Ga._mlt_a[:self.nargs]) ^ X(*X.Ga._mlt_a[self.nargs:nargs])
return Mlt(value, self.Ga, nargs)
else:
return Mlt(X * self.fvalue, self.Ga, self.nargs)

def __or__(self, X):
if isinstance(X, Mlt):
nargs = self.nargs + X.nargs
Mlt.increment_slots(nargs, self.Ga)
value = self(*self.Ga._mlt_a[:self.nargs]) | X(*X.Ga._mlt_a[self.nargs:nargs])
return Mlt(value, self.Ga, nargs)
else:
return Mlt(X * self.fvalue, self.Ga, self.nargs)

def dd(self):
Mlt.increment_slots(self.nargs + 1, self.Ga)
dd_fvalue = (self.Ga._mlt_a[self.nargs] | self.Ga.grad) * self.fvalue
return Mlt(dd_fvalue, self.Ga, self.nargs + 1)

[docs]    def pdiff(self, slot: int):
r"""
Returns gradient of tensor, T, with respect to slot vector.

For example if the tensor is :math:{{T}\lp {a_{1},a_{2}} \rp } then T.pdiff(2) is :math:\nabla_{a_{2}}T. Since T is a scalar function,
T.pdiff(2) is a vector function.
"""
# Take geometric derivative of mlt with respect to slot argument
self.Ga.dslot = slot - 1
return self.Ga.grad * self.Ga.mv(self.fvalue)

@staticmethod
def remove_slot(mv, slot, nargs, ga):
if slot == nargs:
return mv
for islot in range(slot, nargs):
mv = mv.subs(list(zip(ga._mlt_pdiffs[islot], ga._mlt_pdiffs[islot - 1])))
return mv

[docs]    def contract(self, slot1: int, slot2: int):
"""
Returns contraction of tensor between slot1 and slot2 where
slot1 is the index of the first vector argument and slot2 is the
index of the second vector argument of the tensor.

For example if we have a rank two tensor, T(a1, a2), then
T.contract(1, 2) is the contraction of T.
For this case since there are only two slots, there can only be one
contraction.
"""
min_slot = min(slot1, slot2)
max_slot = max(slot1, slot2)
cnargs = self.nargs - 2
self.Ga.dslot = min_slot - 1
self.Ga.dslot = max_slot - 2
div_grad_self = Mlt.remove_slot(div_grad_self.obj, max_slot - 1, self.nargs - 1, self.Ga)
return Mlt(div_grad_self, self.Ga, cnargs)

[docs]    def cderiv(self):
"""
Returns covariant derivative of tensor field.

If T is a tensor of rank :math:k then T.cderiv() is a tensor
of rank :math:k+1. The operation performed is defined in section
:ref:MLtrans.
"""
Mlt.increment_slots(self.nargs + 1, self.Ga)
CD = Mlt((agrad * self.Ga.mv(self.fvalue)).obj, self.Ga, self.nargs + 1)
if CD != 0:
CD = CD.fvalue
for i in range(self.nargs):
args = self.Ga._mlt_a[:self.nargs]
tmp = agrad * self.Ga._mlt_a[i]
if tmp.obj != 0:
args[i] = tmp
CD = CD - self(*args)
CD = Mlt(CD, self.Ga, self.nargs + 1)
return CD

def expand(self):
self.fvalue = expand(self.fvalue)
return self

def comps(self):
basis = self.Ga.mv()
rank = self.nargs
ndim = len(basis)
i_indexes = itertools.product(list(range(ndim)), repeat=rank)
indexes = itertools.product(basis, repeat=rank)
output = ''
for i, (e, i_index) in enumerate(zip(indexes, i_indexes)):
if i_index[-1] % ndim == 0:
print('')
output += str(i)+':'+str(i_index)+':'+str(self(*e)) + '\n'
return output