galgebra/lt.py
"""
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
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
}
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()
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 {
blade: reduce(
self.Ga.wedge,
(self.Ga.basis[i].xreplace(self.lt_dict) for i in index),
S.One
)
for index, blade in self.Ga.indexes_to_blades_dict.items()
}
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)
def __add__(self, LT):
if self.Ga != LT.Ga:
raise ValueError("Attempting addition of Lt's from different geometric algebras")
self_add_LT = copy(self.lt_dict)
for key in list(LT.lt_dict.keys()):
if key in self_add_LT:
self_add_LT[key] = metric.collect(self_add_LT[key] + LT.lt_dict[key], self.Ga.basis)
else:
self_add_LT[key] = LT.lt_dict[key]
return Lt(self_add_LT, ga=self.Ga)
def __sub__(self, LT):
if self.Ga != LT.Ga:
raise ValueError("Attempting subtraction of Lt's from different geometric algebras")
self_add_LT = copy(self.lt_dict)
for key in list(LT.lt_dict.keys()):
if key in self_add_LT:
self_add_LT[key] = metric.collect(self_add_LT[key] - LT.lt_dict[key], self.Ga.basis)
else:
self_add_LT[key] = -LT.lt_dict[key]
return Lt(self_add_LT, ga=self.Ga)
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')
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
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
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`.
"""
self_adj = []
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:
self_adj.append(expand(s))
else:
self_adj.append(expand(s) / self.Ga.e_sq)
return Lt(self_adj, ga=self.Ga)
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)
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
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])
term_add_flg = isinstance(term[0], Add)
if term_add_flg:
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
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 = inspect.getfullargspec(f)[0]
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
raise NotImplementedError
# self.f = None
# self.nargs = len(args) # args isn't defined, which is why we raise NotImplementedError
# Mlt.increment_slots(self.nargs, Ga)
# self.fvalue = f
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)
def __add__(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 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)
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
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
grad_self = self.Ga.grad * self.Ga.mv(self.fvalue)
grad_self = Mlt.remove_slot(grad_self.obj, min_slot, self.nargs, self.Ga)
self.Ga.dslot = max_slot - 2
div_grad_self = self.Ga.grad | self.Ga.mv(grad_self)
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)
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)
agrad = self.Ga._mlt_a[self.nargs] | self.Ga.grad
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