"""
Metric Tensor and Derivatives of Basis Vectors.
"""
import sys
import copy
import itertools
from sympy import diff, trigsimp, Matrix, Rational, \
sqf_list, Symbol, sqrt, eye, zeros, S, expand, Mul, \
Add, simplify, together, ratsimp, Expr, latex, \
Function
from . import printer
from . import utils
half = Rational(1, 2)
def apply_function_list(f,x):
if isinstance(f,(tuple,list)):
fx = x
for fi in f:
fx = fi(fx)
return fx
else:
return f(x)
def str_to_lst(s):
if '[' in s:
s = s.replace('[', '')
if ']' in s:
s = s.replace(']', '')
s_lst = s.split(',')
v_lst = []
for x in s_lst:
try:
v_lst.append(int(s))
except ValueError:
v_lst.append(Symbol(s, real=True))
return v_lst
def linear_expand(expr, mode=True):
if isinstance(expr, Expr):
expr = expand(expr)
if expr == 0:
coefs = [expr]
bases = [S(1)]
return (coefs, bases)
if isinstance(expr, Add):
args = expr.args
else:
if expr.is_commutative:
return ([expr], [S(1)])
else:
args = [expr]
coefs = []
bases = []
for term in args:
if term.is_commutative:
if S(1) in bases:
coefs[bases.index(S(1))] += term
else:
bases.append(S(1))
coefs.append(term)
else:
c, nc = term.args_cnc()
base = nc[0]
coef = Mul._from_args(c)
if base in bases:
coefs[bases.index(base)] += coef
else:
bases.append(base)
coefs.append(coef)
if mode:
return (coefs, bases)
else:
return list(zip(coefs, bases))
[docs]def square_root_of_expr(expr):
"""
If expression is product of even powers then every power is divided
by two and the product is returned. If some terms in product are
not even powers the sqrt of the absolute value of the expression is
returned. If the expression is a number the sqrt of the absolute
value of the number is returned.
"""
if expr.is_number:
if expr > 0:
return(sqrt(expr))
else:
return(sqrt(-expr))
else:
expr = trigsimp(expr)
(coef, pow_lst) = sqf_list(expr)
if coef != S(1):
if coef.is_number:
coef = square_root_of_expr(coef)
else:
coef = sqrt(abs(coef)) # Product coefficient not a number
for p in pow_lst:
(f, n) = p
if n % 2 != 0:
return(sqrt(abs(expr))) # Product not all even powers
else:
coef *= f ** (n / 2) # Positive sqrt of the square of an expression
return coef
def symbols_list(s, indices=None, sub=True, commutative=False):
if isinstance(s, list): # s is already a list of symbols
return(s)
if sub is True: # subscripted list
pos = '_'
else: # superscripted list
pos = '__'
if indices is None: # symbol list completely generated by s
if '*' in s:
[base, index] = s.split('*')
if '|' in s:
index = index.split('|')
s_lst = [base + pos + i for i in index]
else: # symbol list indexed with integers 0 to n-1
try:
n = int(index)
except ValueError:
raise ValueError(index + 'is not an integer')
s_lst = [base + pos + str(i) for i in range(n)]
else:
if ',' in s:
s_lst = s.split(',')
else:
s_lst = s.split(' ')
if not sub:
s_lst = [x.replace('_', '__', 1) for x in s_lst]
else: # indices symbol list used for sub/superscripts of generated symbol list
s_lst = [s + pos + str(i) for i in indices]
return [Symbol(printer.Eprint.Base(s), commutative=commutative) for s in s_lst]
[docs]def test_init_slots(init_slots, **kwargs):
"""
Tests kwargs for allowed keyword arguments as defined by dictionary
init_slots. If keyword argument defined by init_slots is not present
set default value asdefined by init_slots. Allow for backward
compatible keyword arguments by equivalencing keywords by setting
default value of backward compatible keyword to new keyword and then
referencing new keywork (see init_slots for Metric class and equivalence
between keywords 'g' and 'metric')
"""
for slot in kwargs:
if slot not in init_slots:
print('Allowed keyed input arguments')
for key in init_slots:
print(key + ': ' + init_slots[key][1])
raise ValueError('"' + slot + ' = " not in allowed values.')
for slot in init_slots:
if slot in kwargs:
if init_slots[slot][0] in init_slots: # redirect for backward compatibility
kwargs[init_slots[slot][0]] = kwargs[slot]
else: # use default value
if init_slots[slot][0] in init_slots: # redirect for backward compatibility
kwargs[init_slots[slot][0]] = init_slots[init_slots[slot][0]][0]
kwargs[slot] = init_slots[slot][0]
return kwargs
class Simp:
modes = [simplify]
@staticmethod
def profile(s):
Simp.modes = s
return
@staticmethod
def apply(expr):
(coefs, bases) = linear_expand(expr)
obj = S(0)
if isinstance(Simp.modes, list) or isinstance(Simp.modes, tuple):
for (coef, base) in zip(coefs, bases):
for mode in Simp.modes:
coef = mode(coef)
obj += coef * base
else:
for (coef, base) in zip(coefs, bases):
obj += Simp.modes(coef) * base
return obj
@staticmethod
def applymv(mv):
mv.obj = Simp.apply(mv.obj)
return mv
[docs]class Metric(object):
"""
Data Variables -
g[,]: metric tensor (sympy matrix)
g_inv[,]: inverse of metric tensor (sympy matirx)
norm: normalized diagonal metric tensor (list of sympy numbers)
coords[]: coordinate variables (list of sympy symbols)
is_ortho: True if basis is orthogonal (bool)
connect_flg: True if connection is non-zero (bool)
basis[]: basis vector symbols (list of non-commutative sympy variables)
r_symbols[]: reciprocal basis vector symbols (list of non-commutative sympy variables)
n: dimension of vector space/manifold (integer)
n_range: list of basis indices
de[][]: derivatives of basis functions. Two dimensional list. First
entry is differentiating coordiate. Second entry is basis
vector. Quantities are linear combinations of basis vector
symbols.
sig: Signature of metric (p,q) where n = p+q. If metric tensor
is numerical and orthogonal it is calculated. Otherwise the
following inputs are used -
input signature Type
"e" (n,0) Euclidean
"m+" (n-1,1) Minkowski (One negative square)
"m-" (1,n-1) Minkowski (One positive square)
p (p,n-p) General (integer not string input)
"""
count = 1
init_slots = {'g': (None, 'metric tensor'),
'coords': (None, 'manifold/vector space coordinate list/tuple'),
'X': (None, 'vector manifold function'),
'norm': (False, 'True to normalize basis vectors'),
'debug': (False, 'True to print out debugging information'),
'gsym': (None, 'String s to use "det("+s+")" function in reciprocal basis'),
'sig': ('e', 'Signature of metric, default is (n,0) a Euclidean metric'),
'Isq': ('-', "Sign of square of pseudo-scalar, default is '-'"),
'wedge': (True, 'Use ^ symbol to print basis blades')}
[docs] @staticmethod
def dot_orthogonal(V1, V2, g=None):
"""
Returns the dot product of two vectors in an orthogonal coordinate
system. V1 and V2 are lists of sympy expressions. g is
a list of constants that gives the signature of the vector space to
allow for non-euclidian vector spaces.
This function is only used to form the dot product of vectors in the
embedding space of a vector manifold or in the case where the basis
vectors are explicitly defined by vector fields in the embedding
space.
A g of None is for a Euclidian embedding space.
"""
if g is None:
dot = 0
for (v1, v2) in zip(V1, V2):
dot += v1 * v2
return dot
else:
if len(g) == len(V1):
dot = 0
for (v1, v2, gii) in zip(V1, V2, g):
dot += v1 * v2 * gii
return dot
else:
raise ValueError('In dot_orthogonal dimension of metric ' +
'must equal dimension of vector')
[docs] def metric_symbols_list(self, s=None): # input metric tensor as string
"""
rows of metric tensor are separated by "," and elements
of each row separated by " ". If the input is a single
row it is assummed that the metric tensor is diagonal.
Output is a square matrix.
"""
if s is None:
s = self.n * '# '
s = self.n * (s[:-1] + ',')
s = s[:-1]
if utils.isstr(s):
rows = s.split(',')
n_rows = len(rows)
if n_rows == 1: # orthogonal metric
m_lst = s.split(' ')
m = []
for (s, base) in zip(m_lst, self.basis):
if s == '#':
s_symbol = Symbol('(' + str(base) + '.' + str(base) + ')', real=True)
else:
if '/' in s:
[num, dem] = s.split('/')
s_symbol = Rational(num, dem)
else:
s_symbol = Rational(s)
m.append(s_symbol)
if len(m) != self.n:
raise ValueError('Input metric "' + s + '" has' +
' different rank than bases "' + str(self.basis) + '"')
diagonal = eye(self.n)
for i in self.n_range:
diagonal[i, i] = m[i]
return diagonal
else: # non orthogonal metric
rows = s.split(',')
n_rows = len(rows)
m_lst = []
for row in rows:
cols = row.strip().split(' ')
n_cols = len(cols)
if n_rows != n_cols: # non square metric
raise ValueError("'" + s + "' does not represent square metric")
m_lst.append(cols)
m = []
n = len(m_lst)
if n != self.n:
raise ValueError('Input metric "' + s + '" has' +
' different rank than bases "' + str(self.basis) + '"')
n_range = list(range(n))
for (row, i1) in zip(m_lst, n_range):
row_symbols = []
for (s, i2) in zip(row, n_range):
if s == '#':
if i1 <= i2: # for default elment insure symmetry
row_symbols.append(Symbol('(' + str(self.basis[i1]) +
'.' + str(self.basis[i2]) + ')', real=True))
else:
row_symbols.append(Symbol('(' + str(self.basis[i2]) +
'.' + str(self.basis[i1]) + ')', real=True))
else:
if '/' in s: # element is fraction
[num, dem] = s.split('/')
row_symbols.append(Rational(num, dem))
else: # element is integer
row_symbols.append(Rational(s))
m.append(row_symbols)
m = Matrix(m)
return m
def derivatives_of_g(self):
# dg[i][j][k] = \partial_{x_{k}}g_{ij}
dg = [[[
diff(self.g[i, j], x_k)
for x_k in self.coords]
for j in self.n_range]
for i in self.n_range]
return dg
def init_connect_flg(self):
# See if metric is flat
self.connect_flg = False
for i in self.n_range:
for j in self.n_range:
for k in self.n_range:
if self.dg[i][j][k] != 0:
self.connect_flg = True
break
def derivatives_of_basis(self): # Derivatives of basis vectors from Christoffel symbols
n_range = self.n_range
self.dg = dg = self.derivatives_of_g()
self.init_connect_flg()
if not self.connect_flg:
self.de = None
return
de = [] # de[i][j] = \partial_{x_{i}}e^{x_{j}}
# Christoffel symbols of the first kind, \Gamma_{ijk}
# TODO handle None
dG = self.Christoffel_symbols(mode=1)
# \frac{\partial e_{j}}{\partial x^{i}} = \Gamma_{ijk} e^{k}
de = [[
sum([Gamma_ijk * e__k for (Gamma_ijk, e__k) in zip(dG[i][j], self.r_symbols)])
for j in n_range]
for i in n_range]
if self.debug:
printer.oprint('D_{i}e^{j}', de)
self.de = de
return
def inverse_metric(self):
if self.g_inv is not None:
return
if self.is_ortho: # Orthogonal metric
self.g_inv = eye(self.n)
for i in range(self.n):
self.g_inv[i,i] = S(1)/self.g(i,i)
else:
if self.gsym is None:
self.g_inv = simplify(self.g.inv())
else:
self.detg = Function('|' +self.gsym +'|',real=True)(*self.coords)
self.g_adj = simplify(self.g.adjugate())
self.g_inv = self.g_adj/self.detg
return
[docs] def Christoffel_symbols(self,mode=1):
"""
mode = 1 Christoffel symbols of the first kind
mode = 2 Christoffel symbols of the second kind
"""
# See if connection is zero
if not self.connect_flg:
return
n_range = self.n_range
# dg[i][j][k] = \partial_{x_{k}}g_{ij}
dg = self.dg
if mode == 1:
dG = [] # dG[i][j][k] = half * (dg[j][k][i] + dg[i][k][j] - dg[i][j][k])
# Christoffel symbols of the first kind, \Gamma_{ijk}
# \partial_{x^{i}}e_{j} = \Gamma_{ijk}e^{k}
def Gamma_ijk(i, j, k):
return half * (dg[j][k][i] + dg[i][k][j] - dg[i][j][k])
dG = [[[
Simp.apply(Gamma_ijk(i, j, k))
for k in n_range]
for j in n_range]
for i in n_range]
if self.debug:
printer.oprint('Gamma_{ijk}', dG)
return dG
elif mode == 2:
# TODO handle None
Gamma1 = self.Christoffel_symbols(mode=1)
self.inverse_metric()
# Christoffel symbols of the second kind, \Gamma_{ij}^{k} = \Gamma_{ijl}g^{lk}
# \partial_{x^{i}}e_{j} = \Gamma_{ij}^{k}e_{k}
def Gamma2_ijk(i, j, k):
return sum([Gamma_ijl * self.g_inv[l, k] for l, Gamma_ijl in enumerate(Gamma1[i][j])])
Gamma2 = [[[
Simp.apply(Gamma2_ijk(i, j, k))
for k in n_range]
for j in n_range]
for i in n_range]
return Gamma2
else:
raise ValueError('In Christoffle_symobols mode = ' + str(mode) +' is not allowed\n')
def normalize_metric(self):
if self.de is None:
return
renorm = []
# Generate mapping for renormalizing reciprocal basis vectors
for ib in self.n_range: # e^{ib} --> e^{ib}/|e_{ib}|
renorm.append((self.r_symbols[ib], self.r_symbols[ib] / self.e_norm[ib]))
# Normalize derivatives of basis vectors
for x_i in self.n_range:
for jb in self.n_range:
self.de[x_i][jb] = Simp.apply((((self.de[x_i][jb].subs(renorm)
- diff(self.e_norm[jb], self.coords[x_i]) *
self.basis[jb]) / self.e_norm[jb])))
if self.debug:
for x_i in self.n_range:
for jb in self.n_range:
print(r'\partial_{' + str(self.coords[x_i]) + r'}\hat{e}_{' + str(self.coords[jb]) + '} =', self.de[x_i][jb])
# Normalize metric tensor
for ib in self.n_range:
for jb in self.n_range:
self.g[ib, jb] = Simp.apply(self.g[ib, jb] / (self.e_norm[ib] * self.e_norm[jb]))
if self.debug:
printer.oprint('e^{i}->e^{i}/|e_{i}|', renorm)
printer.oprint('renorm(g)', self.g)
return
def signature(self):
if self.is_ortho:
p = 0
q = 0
for i in self.n_range:
g_ii = self.g[i,i]
if g_ii.is_number:
if g_ii > 0:
p += 1
else:
q += 1
else:
break
if p + q == self.n:
self.sig = (p,q)
return
if isinstance(self.sig,int): # General signature
if self.sig <= self.n:
self.sig = (self.sig,self.n - self.sig)
return
else:
raise ValueError('self.sig = ' + str(self.sig) + ' > self.n, not an allowed hint')
if utils.isstr(self.sig):
if self.sig == 'e': # Euclidean metric signature
self.sig = (self.n, 0)
elif self.sig == 'm+': # Minkowski metric signature (n-1,1)
self.sig = (self.n - 1, 1)
elif self.sig == 'm-': # Minkowski metric signature (1,n-1)
self.sig = (1, self.n - 1)
else:
raise ValueError('self.sig = ' + str(self.sig) + ' is not an allowed hint')
return
raise ValueError(str(self.sig) + ' is not allowed value for self.sig')
def __init__(self, basis, **kwargs):
kwargs = test_init_slots(Metric.init_slots, **kwargs)
self.name = 'GA' + str(Metric.count)
Metric.count += 1
if not utils.isstr(basis):
raise TypeError('"' + str(basis) + '" must be string')
X = kwargs['X'] # Vector manifold
g = kwargs['g'] # Explicit metric or base metric for vector manifold
debug = kwargs['debug']
coords = kwargs['coords'] # Manifold coordinates (sympy symbols)
norm = kwargs['norm'] # Normalize basis vectors
self.sig = kwargs['sig'] # Hint for metric signature
"""
String for symbolic metric determinant. If self.gsym = 'g'
then det(g) is sympy scalar function of coordinates with
name 'det(g)'. Useful for complex non-orthogonal coordinate
systems or for calculations with general metric.
"""
self.gsym = kwargs['gsym']
self.Isq = kwargs['Isq'] # Sign of I**2, only needed if I**2 not a number
self.debug = debug
self.is_ortho = False # Is basis othogonal
self.coords = coords # Manifold coordinates
if self.coords is None:
self.connect_flg = False
else:
self.connect_flg = True # Connection needed for postion dependent metric
self.norm = norm # True to normalize basis vectors
self.detg = None # Determinant of g
self.g_adj = None # Adjugate of g
self.g_inv = None # Inverse of g
# Generate list of basis vectors and reciprocal basis vectors
# as non-commutative symbols
if ' ' in basis or ',' in basis or '*' in basis: # bases defined by substrings separated by spaces or commas
self.basis = symbols_list(basis)
self.r_symbols = symbols_list(basis, sub=False)
else:
if coords is not None: # basis defined by root string with symbol list as indices
self.basis = symbols_list(basis, coords)
self.r_symbols = symbols_list(basis, coords, sub=False)
self.coords = coords
if self.debug:
printer.oprint('x^{i}', self.coords)
else:
raise ValueError('for basis "' + basis + '" coords must be entered')
if self.debug:
printer.oprint('e_{i}', self.basis, 'e^{i}', self.r_symbols)
self.n = len(self.basis)
self.n_range = list(range(self.n))
# Generate metric as list of lists of symbols, rationals, or functions of coordinates
if g is None:
if X is None: # default metric from dot product of basis as symbols
self.g = self.metric_symbols_list()
else: # Vector manifold
if coords is None:
raise ValueError('For metric derived from vector field ' +
' coordinates must be defined.')
else: # Vector manifold defined by vector field
dX = []
for coord in coords: # Get basis vectors by differentiating vector field
dX.append([diff(x, coord) for x in X])
g_tmp = []
for dx1 in dX:
g_row = []
for dx2 in dX:
dx1_dot_dx2 = trigsimp(Metric.dot_orthogonal(dx1, dx2, g))
g_row.append(dx1_dot_dx2)
g_tmp.append(g_row)
self.g = Matrix(g_tmp)
if self.debug:
printer.oprint('X_{i}', X, 'D_{i}X_{j}', dX)
else: # metric is symbolic or list of lists of functions of coordinates
if utils.isstr(g): # metric elements are symbols or constants
if g == 'g': # general symbolic metric tensor (g_ij functions of position)
g_lst = []
g_inv_lst = []
for coord in self.coords:
i1 = str(coord)
tmp = []
tmp_inv = []
for coord2 in self.coords:
i2 = str(coord2)
tmp.append(Function('g_'+i1+'_'+i2)(*self.coords))
tmp_inv.append(Function('g__'+i1+'__'+i2)(*self.coords))
g_lst.append(tmp)
g_inv_lst.append(tmp_inv)
self.g = Matrix(g_lst)
self.g_inv = Matrix(g_inv_lst)
else: # specific symbolic metric tensor (g_ij are symbolic or numerical constants)
self.g = self.metric_symbols_list(g) # construct symbolic metric from string and basis
else: # metric is given as list of function or list of lists of function or matrix of functions
if isinstance(g, Matrix):
self.g = g
else:
if isinstance(g[0], list):
self.g = Matrix(g)
else:
m = eye(len(g))
for i in range(len(g)):
m[i, i] = g[i]
self.g = m
self.g_raw = copy.deepcopy(self.g) # save original metric tensor for use with submanifolds
if self.debug:
printer.oprint('g', self.g)
# Determine if metric is orthogonal
self.is_ortho = True
for i in self.n_range:
for j in self.n_range:
if i < j:
if self.g[i, j] != 0:
self.is_ortho = False
break
self.g_is_numeric = True
for i in self.n_range:
for j in self.n_range:
if i < j:
if not self.g[i, j].is_number:
self.g_is_numeric = False
break
if self.coords is not None:
self.derivatives_of_basis() # calculate derivatives of basis
if self.norm: # normalize basis, metric, and derivatives of normalized basis
if not self.is_ortho:
raise ValueError('!!!!Basis normalization only implemented for orthogonal basis!!!!')
self.e_norm = []
for i in self.n_range:
self.e_norm.append(square_root_of_expr(self.g[i, i]))
if debug:
printer.oprint('|e_{i}|', self.e_norm)
else:
self.e_norm = None
if self.norm:
if self.is_ortho:
self.normalize_metric()
else:
raise ValueError('!!!!Basis normalization only implemented for orthogonal basis!!!!')
if not self.g_is_numeric:
self.signature()
# Sign of square of pseudo scalar
self.e_sq_sgn = '+'
if ((self.n*(self.n-1))//2+self.sig[1])%2 == 1:
self.e_sq_sgn = '-'
if self.debug:
print('signature =', self.sig)
if __name__ == "__main__":
pass