"""
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, simplify
)
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
})
# ## GSG code starts ###
[docs]
def Symbolic_Matrix(kernel, coords=None, f=False, mode='g'):
"""
Returns a square real matrix the entries of which are symbolic
constants or symbolic functions of the coordinates.
- `kernel` is a one-letter string. It specifies the kernel letter of
indexed symbols or functions used to specify the matrix's entries
- `coords` is a list or tuple. Its entries are used to label the
components of a vector.
- `f`, a boolean, specifies that matrix entries are symbolic functions
of the coordinates or are symbolic constants, according to whether
`f` is True or False.
- `mode` is a one-letter string. When`mode` is 'g', 's', or 'a' the
matrix will be general, symmetric, or antisymmetric.
"""
def general_matrix(kernel, coords=None, f=False):
"""Returns a general square matrix. The row index of each entry
appears as a superscript, while the column index appears as a
subscript."""
n = len(coords)
# Create matrix entries and store in appropriate locations in `G`:
G = zeros(n, n)
if f: # entries are symbolic functions
for i in range(n):
for j in range(n):
entry = '{' + kernel + '__' + str(coords[i]) + '}_' + str(coords[j])
G[i, j] = Function(entry)(*coords)
else: # entries are symbolic constants
for i in range(n):
for j in range(n):
entry = '{' + kernel + '__' + str(coords[i]) + '}_' + str(coords[j])
G[i, j] = Symbol(entry, real=True)
return G
def symmetric_matrix(kernel, coords=None, f=False):
"""Returns a symmetric matrix. Entries have a single index, which
appears as a subscript."""
n = len(coords)
# Create and temporarily store matrix entries in `parameters`
parameters = []
if f: # entries are symbolic functions
for i in range((n*(n+1)//2), 0, -1):
parameters.append(Function(kernel + '_' + str(i))(*coords))
else: # entries are symbolic constants
for i in range((n*(n+1)//2), 0, -1):
parameters.append(Symbol(kernel + '_' + str(i), real=True))
# Transfer entries to symmetric matrix `S`.
S = zeros(n, n)
for i in range(n):
for j in range(i, n):
S[i, j] = parameters.pop()
S[j, i] = S[i, j]
return S
def antisymmetric_matrix(kernel, coords=None, f=False):
"""Returns an antisymmetric matrix. Entries have a a single index,
which appears as a subscript."""
n = len(coords)
# Create and temporarily store matrix entries in `parameters`
parameters = []
if f: # entries are symbolic functions
for i in range((n*(n-1)//2), 0, -1):
parameters.append(Function(kernel + '_' + str(i))(*coords))
else: # entries are symbolic constants
for i in range((n*(n-1)//2), 0, -1): # each parameter is a symbol
parameters.append(Symbol(kernel + '_' + str(i), real=True))
# Transfer entries to antisymmetric matrix `A`.
A = zeros(n, n)
for i in range(n):
for j in range(i+1, n):
A[i, j] = parameters.pop()
A[j, i] = - A[i, j]
return A
# Check legitimacy of parameter values:
if not isinstance(coords, (list, tuple)):
raise ValueError('coords = ' + str(coords) + ' in Symbolic_Matrix')
if mode not in ['g', 's', 'a']:
raise ValueError('mode = ' + str(mode) + ' in Symbolic_Matrix')
if mode == 'g':
return general_matrix(kernel, coords, f)
if mode == 's':
return symmetric_matrix(kernel, coords, f)
if mode == 'a':
return antisymmetric_matrix(kernel, coords, f)
# ## GSG code ends ###
[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
}
# ## GSG code starts ###
[docs]
def Dictionary_to_Matrix(dict_rep, ga):
"""Returns the matrix representation of that linear transformation on
geometric algebra ga which has dictionary representation dict_rep."""
# columns[j] is a list of the entries in the matrix's jth column.
# columns[j][i] is the (i,j)th entry in the matrix.
# Matrix[columns] instantiates the transpose of the desired matrix.
columns = []
for b in ga.basis: # b is a basis symbol for ga.
column = ga.n * [S.Zero] # Initialize column for dict_rep value at b.
dict_value = dict_rep[b] # dict_rep's value at b
if isinstance(dict_value, mv.Mv):
dict_value = dict_value.obj
if dict_value is not S.Zero:
for coef, base in metric.linear_expand_terms(dict_value):
row_index = ga.basis.index(base)
column[row_index] = coef
columns.append(column)
return Transpose(Matrix(columns)).doit()
# ## GSG code ends ###
[docs]
class Lt(printer.GaPrintable):
r"""
A Linear Transformation
Except for the versor 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 versor representation, the linear transformation is
stored as a versor ``self.V`` so that if a is a
vector::
self(a) = self.V.g_invol() * a * self.V.inv()
where ``self.V.g_invol()`` is the grade involute of ``self.V``.
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.
# ## GSG code starts ###
.. class:: Lt(lt_list: list, /, *, ga)
:noindex:
Construct from a list of lists, the j_th list of which contains
the coefficients of j_th image vector's basis expansion.
# ## GSG code ends ###
.. class:: Lt(versor: mv.Mv, /, *, ga)
:noindex:
Construct from a not-necessarily-normalized versor.
.. 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 both domain and codomain of this transformation
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.lt_dict = {}
self.mat = None
self.versor = False
# self.V, self.Vrev, and self.Vqform are never actually used in the current
# implementation of orthogonal outermorphisms created by a versor input.
self.V = None
self.Vrev = None
self.Vqform = 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 input
if not isinstance(mat_rep[0], list):
# At this point mat_rep[i] is the desired image vector for the
# i_th basis image vectors.
for lt_i, base in zip(mat_rep, self.Ga.basis):
self.lt_dict[base] = sym(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
# ## GSG code starts ###
elif isinstance(mat_rep, Matrix): # Matrix input
self.lt_dict = Matrix_to_dictionary(mat_rep, self.Ga.basis)
# ## GSG code ends ###
# ## GSG code starts. ###
# This code segment uses versor `mat_rep` and a sandwich product only to
# create a linear vector-valued function of vector. That function is then
# used to create a dictionary-based outermorphism. Evaluation of the
# outermorphism on a multivector is by dictionary lookup, not by a
# sandwich product of the multivector with the versor.
elif isinstance(mat_rep, mv.Mv): # Versor input
if not mat_rep.is_versor:
raise ValueError(mat_rep, 'is not a versor in Versor input for Lt!\n')
V = mat_rep
Vg_invol = V.g_invol()
Vinv = V.inv()
outermorphism = ga.lt(lambda x: Vg_invol * x * Vinv)
self.lt_dict = simplify(outermorphism.lt_dict)
# ## GSG code ends ###
# ## GSG code starts ###
elif isinstance(mat_rep, str): # (One-letter) string input
Amat = Symbolic_Matrix(mat_rep, coords=self.Ga.coords, f=f, mode=mode)
if mode == 'g':
self.__init__(Amat, ga=self.Ga)
elif mode in ['s', 'a']:
self.__init__(self.Ga.g_inv * Amat, ga=self.Ga)
# ## GSG code ends ###
elif callable(mat_rep): # Linear multivector function input
# Function is tested for linearity before use.
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.versor:
# 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()
}
[docs]
def __call__(self, M, obj=False):
r"""
Returns the image of multivector :math:`M` under the linear transformation
:math:`L`. :math:`{{L}\lp{M}\rp}` is defined by
the linearity of :math:`L`,
the vector values :math:`{{L}\lp{{{\eb}}_{j}}\rp }`, and the definition
:math:`{{L}\lp{{{\eb}}_{j_{1}}{\wedge}\dots{\wedge}{{\eb}}_{j_{r}}}\rp}={{L}\lp{{{\eb}}_{j_{1}}}\rp}{\wedge}\dots{\wedge}{{L}\lp{{{\eb}}_{j_{r}}}\rp}`.
"""
if isinstance(M, mv.Mv) and self.Ga != M.Ga:
raise ValueError('In Lt call Lt and argument refer to different vector spaces')
# ## GSG code starts ###
# Given the current way an outermorphism is created from a versor input,
# self.versor will always be false; hence the following code fragment will
# never execute.
if self.versor:
# Sandwich M or M's grade involute depending on whether versor self.V
# is even or odd.
if self.V == self.V.odd():
V_M_Vrev = self.V * M.g_invol() * self.Vrev
elif self.V == self.V.even():
V_M_Vrev = self.V * M * self.Vrev
else:
raise ValueError('self.V is not a versor in __call__')
# Divide by normalization factor self.Vqform to convert sandwiching
# between self.V and its reverse to sandwiching between self.V and
# its inverse.
V_M_Vinv = 1/(self.Vqform) * V_M_Vrev
if obj:
return V_M_Vinv.obj
else:
return V_M_Vinv
# ## GSG code ends ###
if isinstance(M, mv.Mv):
if M.is_vector():
lt_M = M.obj.xreplace(self.lt_dict)
if obj:
return lt_M
else:
return mv.Mv(lt_M, ga=self.Ga)
else:
mv_obj = M.obj
else:
mv_obj = mv.Mv(M, ga=self.Ga).obj
lt_M = mv_obj.xreplace(self.mv_dict)
if obj:
return lt_M
else:
return mv.Mv(lt_M, 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')
# ## GSG code starts ###
[docs]
def det(self) -> Expr: # det(L) defined by L(E) = det(L)E
r"""
- Returns the determinant of the linear transformation :math:`L`,
defined by :math:`\det(L) = L(E) E^{-1}`, where :math:`E` is the
basis blade for the pseudoscalar grade space.
- Expression returned is a real SymPy scalar, not a GAlgebra 0-vector.
"""
return (self(self.Ga.e) * self.Ga.e.inv()).scalar()
# ## GSG code ends ###
[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
r'''
def adj(self) -> 'Lt':
r"""
Returns the adjoint :math:`{\bar{L}}`(a linear transformation) of 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.
"""
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)
'''
# ## GSG code starts ###
[docs]
def adj(self) -> 'Lt':
r"""
Returns the adjoint transformation :math:`{\bar{L}}` of 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.
"""
matrix_of_adjoint = self.Ga.g_inv * self.matrix().T * self.Ga.g
return self.Ga.lt(matrix_of_adjoint)
# ## GSG code ends ###
# ## GSG code starts ###
[docs]
def is_singular(self):
"""Returns `True` if and only if linear transformation `self` is singular."""
E = self.Ga.E()
return simplify((self(E) < E.inv()).scalar()) == S.Zero
# ## GSG code ends
# ## GSG code starts ###
[docs]
def inv(self):
"""Returns compositional inverse of linear transformation`self`.
Assumes transformation is nonsingular. If `self` is a versor based
transformation, its inverse will also be versor based."""
if self.versor:
return self.Ga.lt(self.V.rev())
if not self.is_singular():
return self.Ga.lt(Matrix(self.matrix().inv()))
else:
raise ValueError('transformation in inv() is non-invertible')
# ## GSG code ends ###
def _sympystr(self, print_obj):
if self.versor: # ## GSG: changed `self.spinor` to `self.versor` ###
return 'R = ' + print_obj._print(self.V)
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]
# ## GSG code starts ###
def _latex(self, print_obj):
parts = []
for base in self.Ga.basis: # base is a basis symbol
if self.versor:
b = mv.Mv(base, ga=self.Ga) # b is the corresponding basis vector
if self.V == self.V.odd():
unnormalized_image = self.V * (b.g_invol()) * self.Vrev
elif self.V == self.V.even():
unnormalized_image = self.V * b * self.Vrev
else:
raise ValueError('self.V is not a versor in _latex')
image = 1/(self.Vqform) * unnormalized_image
else:
image = mv.Mv(self.lt_dict.get(base, S.Zero), ga=self.Ga)
parts.append(print_obj._print(base) + ' &\\mapsto ' + print_obj._print(image))
return '\\left\\{ \\begin{aligned} ' + ' \\\\ '.join(parts) + ' \\end{aligned} \\right\\}'
# ## GSG code ends ###
def Fmt(self, fmt=1, title=None) -> printer.GaPrintable:
return printer._FmtResult(self, title)
# ## GSG code starts ###
[docs]
def matrix(self) -> Matrix:
r"""
Returns the matrix :math:`[{L__i}_j]` defined for linear transformation
:math:`L` by :math:`L({\eb}_j)=\sum_i {L__i}_j \eb}_i`.
"""
if self.mat is not None:
return self.mat
elif self.versor:
self.lt_dict = {}
for base in self.Ga.basis:
self.lt_dict[base] = self(base).simplify()
self.versor = False # temporary change of self.versor
mat = self.matrix()
self.versor = True # reverse change to self.versor
return mat
else:
self.mat = Dictionary_to_Matrix(self.lt_dict, self.Ga)
return self.mat.doit()
# ## GSG code ends ###
[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])
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
[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 = 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
[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)
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)
[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
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)
[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)
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
# ## GSG code starts ###
[docs]
def det(L: Lt) -> Expr: # det(L) defined by L(E) = det(L)E
r"""
- Returns the determinant of the linear transformation :math:`L`,
defined by :math:`\det(L) = L(E) E^{-1}`, where :math:`E` is the
basis blade for the pseudoscalar grade space.
- Expression returned is a real SymPy scalar, not a GAlgebra 0-vector.
"""
return L.det()
# ## GSG code ends ###
# ## GSG code starts ###
[docs]
def sym(v):
"""
Returns that linear combination of basis vector symbols which corresponds
to vector v, itself a linear combination of basis vectors.
"""
# Obtain the coefficients in basis vector expansion of `v`.
# Then construct and return corresponding basis vector symbol expansion.
coefs = v.blade_coefs(v.Ga.mv())
return sum(coefs[j]*v.Ga.basis[j] for j in range(v.Ga.n))
# ## GSG code ends ###