diff --git a/doc/src/modules/tensor/index.rst b/doc/src/modules/tensor/index.rst index 853f490aef7f..f64d3e80dba5 100644 --- a/doc/src/modules/tensor/index.rst +++ b/doc/src/modules/tensor/index.rst @@ -14,3 +14,5 @@ Contents indexed.rst index_methods.rst + tensor.rst + dgamma_matr.rst diff --git a/sympy/combinatorics/tensor_can.py b/sympy/combinatorics/tensor_can.py index 873956ea3424..792142eaaa48 100644 --- a/sympy/combinatorics/tensor_can.py +++ b/sympy/combinatorics/tensor_can.py @@ -1,9 +1,9 @@ from sympy.combinatorics.permutations import Permutation, _af_rmul, _af_rmuln,\ - _af_invert, _af_new + _af_invert, _af_new from sympy.combinatorics.perm_groups import PermutationGroup, _orbit, \ - _orbit_transversal + _orbit_transversal from sympy.combinatorics.util import _distribute_gens_by_base, \ - _orbits_transversals_from_bsgs + _orbits_transversals_from_bsgs """ References for tensor canonicalization: @@ -28,40 +28,48 @@ def dummy_sgs(dummies, sym, n): """ Return the strong generators for dummy indices - dummies list of dummy indices, - `dummies[2k], dummies[2k+1]` are paired indices - sym symmetry under interchange of contracted dummies - sym = None no symmetry - 0 symmetric - 1 antisymmetric - n number of indices + Parameters + ========== + + dummies : list of dummy indices + `dummies[2k], dummies[2k+1]` are paired indices + sym : symmetry under interchange of contracted dummies:: + * None no symmetry + * 0 commuting + * 1 anticommuting + + n : number of indices in base form the dummy indices are always in consecutive positions Examples ======== + >>> from sympy.combinatorics.tensor_can import dummy_sgs - >>> dummy_sgs([2,3,4,5,6,7], 0, 8) - [[0, 1, 3, 2, 4, 5, 6, 7, 8, 9], [0, 1, 2, 3, 5, 4, 6, 7, 8, 9], [0, 1, 2, 3, 4, 5, 7, 6, 8, 9], [0, 1, 4, 5, 2, 3, 6, 7, 8, 9], [0, 1, 2, 3, 6, 7, 4, 5, 8, 9]] + >>> dummy_sgs(range(2, 8), 0, 8) + [[0, 1, 3, 2, 4, 5, 6, 7, 8, 9], [0, 1, 2, 3, 5, 4, 6, 7, 8, 9], + [0, 1, 2, 3, 4, 5, 7, 6, 8, 9], [0, 1, 4, 5, 2, 3, 6, 7, 8, 9], + [0, 1, 2, 3, 6, 7, 4, 5, 8, 9]] """ assert len(dummies) <= n res = [] # exchange of contravariant and covariant indices - if sym != None: + if sym is not None: for j in dummies[::2]: - a = range(n+2) + a = range(n + 2) if sym == 1: - a[n] = n+1 - a[n+1] = n - a[j], a[j+1] = a[j+1], a[j] + a[n] = n + 1 + a[n + 1] = n + a[j], a[j + 1] = a[j + 1], a[j] res.append(a) # rename dummy indices for j in dummies[:-3:2]: - a = range(n+2) - a[j:j+4] = a[j+2],a[j+3],a[j],a[j+1] + a = range(n + 2) + a[j:j + 4] = a[j + 2], a[j + 3], a[j], a[j + 1] res.append(a) return res + def _min_dummies(dummies, sym, indices): """ Return list of minima of the orbits of indices in group of dummies @@ -70,8 +78,9 @@ def _min_dummies(dummies, sym, indices): Examples ======== + >>> from sympy.combinatorics.tensor_can import _min_dummies - >>> _min_dummies([[2,3,4,5,6,7]], [0], range(10)) + >>> _min_dummies([range(2, 8)], [0], range(10)) [0, 1, 2, 2, 2, 2, 2, 2, 8, 9] """ num_types = len(sym) @@ -102,6 +111,7 @@ def _trace_S(s, j, b, S_cosets): return h return None + def _trace_D(gj, p_i, Dxtrav): """ Return the representative h satisfying h[gj] == p_i @@ -113,6 +123,7 @@ def _trace_D(gj, p_i, Dxtrav): return h return None + def _dumx_remove(dumx, dumx_flat, p0): """ remove p0 from dumx @@ -124,15 +135,16 @@ def _dumx_remove(dumx, dumx_flat, p0): continue k = dx.index(p0) if k % 2 == 0: - p0_paired = dx[k+1] + p0_paired = dx[k + 1] else: - p0_paired = dx[k-1] + p0_paired = dx[k - 1] dx.remove(p0) dx.remove(p0_paired) dumx_flat.remove(p0) dumx_flat.remove(p0_paired) res.append(dx) + def transversal2coset(size, base, transversal): a = [] j = 0 @@ -142,10 +154,10 @@ def transversal2coset(size, base, transversal): j += 1 else: a.append([range(size)]) - j = len(a)-1 + j = len(a) - 1 while a[j] == [range(size)]: j -= 1 - return a[:j+1] + return a[:j + 1] def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): @@ -156,7 +168,7 @@ def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): list of lists of dummy indices, one list for each type of index; the dummy indices are put in order contravariant, covariant - [d0, -d0, d1,-d1,...]. + [d0, -d0, d1, -d1, ...]. sym list of the symmetries of the index metric for each type. @@ -379,7 +391,7 @@ def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): g = g.array_form num_dummies = size - 2 indices = range(num_dummies) - all_metrics_with_sym = all([_ != None for _ in sym]) + all_metrics_with_sym = all([_ is not None for _ in sym]) num_types = len(sym) dumx = dummies[:] dumx_flat = [] @@ -392,7 +404,7 @@ def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): # strong generating set for D dsgsx = [] for i in range(num_types): - dsgsx.extend(dummy_sgs(dumx[i], sym[i], num_dummies)) + dsgsx.extend(dummy_sgs(dumx[i], sym[i], num_dummies)) ginv = _af_invert(g) idn = range(size) # TAB = list of entries (s, d, h) where h = _af_rmuln(d,g,s) @@ -410,24 +422,25 @@ def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): if all_metrics_with_sym: md = _min_dummies(dumx, sym, indices) else: - md = [min(_orbit(size, [_af_new(ddx) for ddx in dsgsx], ii)) for ii in range(size-2)] + md = [min(_orbit(size, [_af_new( + ddx) for ddx in dsgsx], ii)) for ii in range(size - 2)] - p_i = min([min([md[h[x]] for x in deltab]) for s,d,h in TAB]) + p_i = min([min([md[h[x]] for x in deltab]) for s, d, h in TAB]) dsgsx1 = [_af_new(_) for _ in dsgsx] Dxtrav = _orbit_transversal(size, dsgsx1, p_i, False, af=True) \ - if dsgsx else None + if dsgsx else None if Dxtrav: Dxtrav = [_af_invert(x) for x in Dxtrav] # compute the orbit of p_i for ii in range(num_types): if p_i in dumx[ii]: # the orbit is made by all the indices in dum[ii] - if sym[ii] != None: + if sym[ii] is not None: deltap = dumx[ii] else: # the orbit is made by all the even indices if p_i # is even, by all the odd indices if p_i is odd - p_i_index = dumx[ii].index(p_i)%2 + p_i_index = dumx[ii].index(p_i) % 2 deltap = dumx[ii][p_i_index::2] break else: @@ -492,7 +505,7 @@ def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): # if TAB contains equal permutations up to the sign, return 0 TAB1.sort(key=lambda x: x[-1]) nTAB1 = len(TAB1) - prev = [0]*size + prev = [0] * size while TAB1: s, d, h = TAB1.pop() if h[:-2] == prev[:-2]: @@ -517,17 +530,19 @@ def canonical_free(base, gens, g, num_free): """ canonicalization of a tensor with respect to free indices choosing the minimum with respect to lexicographical ordering + in the free indices - base, gens BSGS for slot permutation group - g permutation representing the tensor - num_free number of free indices + ``base``, ``gens`` BSGS for slot permutation group + ``g`` permutation representing the tensor + ``num_free`` number of free indices The indices must be ordered with first the free indices see explanation in double_coset_can_rep - The algorithm is given in [2] + The algorithm is a variation of the one given in [2]. Examples ======== + >>> from sympy.combinatorics import Permutation >>> from sympy.combinatorics.tensor_can import canonical_free >>> gens = [[1,0,2,3,5,4], [2,3,0,1,4,5],[0,1,3,2,5,4]] @@ -538,20 +553,23 @@ def canonical_free(base, gens, g, num_free): [0, 3, 1, 2, 5, 4] Consider the product of Riemann tensors - `T = R^{a}_{d0}^{d1,d2}*R_{d2,d1}^{d0,b}` - The order of the indices is [a,b,d0,-d0,d1,-d1,d2,-d2] + ``T = R^{a}_{d0}^{d1,d2}*R_{d2,d1}^{d0,b}`` + The order of the indices is ``[a,b,d0,-d0,d1,-d1,d2,-d2]`` The permutation corresponding to the tensor is - g = [0,3,4,6,7,5,2,1,8,9] - Use the slot symmetries to get `T` is the form which is the minimum - in lexicographic order - `R^{a}_{d0}^{d1,d2}*R^{b,d0}_{d1,d2}` corresponding to - `[0, 3, 4, 6, 1, 2, 5, 7, 8, 9]` + ``g = [0,3,4,6,7,5,2,1,8,9]`` + + In particular ``a`` is position ``0``, ``b`` is in position ``9``. + Use the slot symmetries to get `T` is a form which is the minimal + in lexicographic order in the free indices ``a`` and ``b``, e.g. + ``-R^{a}_{d0}^{d1,d2}*R^{b,d0}_{d2,d1}`` corresponding to + ``[0, 3, 4, 6, 1, 2, 7, 5, 9, 8]`` + >>> from sympy.combinatorics.tensor_can import riemann_bsgs, tensor_gens >>> base, gens = riemann_bsgs >>> size, sbase, sgens = tensor_gens(base, gens, [[],[]], 0) >>> g = Permutation([0,3,4,6,7,5,2,1,8,9]) >>> canonical_free(sbase, [Permutation(h) for h in sgens], g, 2) - [0, 3, 4, 6, 1, 2, 5, 7, 8, 9] + [0, 3, 4, 6, 1, 2, 7, 5, 9, 8] """ g = g.array_form size = len(g) @@ -559,44 +577,44 @@ def canonical_free(base, gens, g, num_free): return g[:] transversals = get_transversals(base, gens) - cosets = transversal2coset(size, base, transversals) - K = [h._array_form for h in gens] m = len(base) for x in sorted(g[:-2]): if x not in base: base.append(x) - lambd = g - K1 = [] - for i in range(m): + h = g + for i, transv in enumerate(transversals): b = base[i] - delta = [x[b] for x in cosets[b]] - delta1 = [lambd[x] for x in delta] - delta1min = min(delta1) - k = delta1.index(delta1min) - p = delta[k] - for omega in cosets[b]: - if omega[b] == p: - break - lambd = _af_rmul(lambd, omega) - K = [px for px in K if px[b] == b] - return lambd + h_i = [size]*num_free + # find the element s in transversals[i] such that + # _af_rmul(h, s) has its free elements with the lowest position in h + s = None + for sk in transv.values(): + h1 = _af_rmul(h, sk) + hi = [h1.index(ix) for ix in range(num_free)] + if hi < h_i: + h_i = hi + s = sk + if s: + h = _af_rmul(h, s) + return h def _get_map_slots(size, fixed_slots): res = range(size) pos = 0 for i in range(size): - if i in fixed_slots: - continue - res[i] = pos - pos += 1 + if i in fixed_slots: + continue + res[i] = pos + pos += 1 return res + def _lift_sgens(size, fixed_slots, free, s): a = [] j = k = 0 fd = zip(fixed_slots, free) - fd = [y for x,y in sorted(fd)] + fd = [y for x, y in sorted(fd)] num_free = len(free) for i in range(size): if i in fixed_slots: @@ -607,44 +625,49 @@ def _lift_sgens(size, fixed_slots, free, s): j += 1 return a + def canonicalize(g, dummies, msym, *v): """ canonicalize tensor formed by tensors - g permutation representing the tensor - - dummies list representing the dummy indices - it can be a list of dummy indices of the same type - or a list of lists of dummy indices, one list for each - type of index; - the dummy indices must come after the free indices, - and put in order contravariant, covariant - [d0, -d0, d1,-d1,...] + Parameters + ========== - msym symmetry of the metric(s) - it can be an integer or a list; - in the first case it is the symmetry of the dummy index metric; - in the second case it is the list of the symmetries of the - index metric for each type + g : permutation representing the tensor - v is a list of (base_i, gens_i, n_i, sym_i) for tensors of type `i` - base_i, gens_i BSGS for tensors of this type. + dummies : list representing the dummy indices + it can be a list of dummy indices of the same type + or a list of lists of dummy indices, one list for each + type of index; + the dummy indices must come after the free indices, + and put in order contravariant, covariant + [d0, -d0, d1,-d1,...] + msym : symmetry of the metric(s) + it can be an integer or a list; + in the first case it is the symmetry of the dummy index metric; + in the second case it is the list of the symmetries of the + index metric for each type + v : list, (base_i, gens_i, n_i, sym_i) for tensors of type `i` - The BSGS should have minimal base under lexicographic ordering; - if not, an attempt is made do get the minimal BSGS; - in case of failure, - canonicalize_naive is used, which is much slower. + base_i, gens_i : BSGS for tensors of this type. + The BSGS should have minimal base under lexicographic ordering; + if not, an attempt is made do get the minimal BSGS; + in case of failure, + canonicalize_naive is used, which is much slower. - n_i number ot tensors of type `i`. + n_i : number of tensors of type `i`. - sym_i symmetry under exchange of component tensors of type `i`. + sym_i : symmetry under exchange of component tensors of type `i`. - Both for msym and sym_i the cases are + Both for msym and sym_i the cases are * None no symmetry * 0 commuting * 1 anticommuting - Return 0 if the tensor is zero, else return the array form of + Returns + ======= + + 0 if the tensor is zero, else return the array form of the permutation representing the canonical form of the tensor. Algorithm @@ -677,7 +700,7 @@ def canonicalize(g, dummies, msym, *v): g = [1,3,0,5,4,2,6,7] - T_c = 0 + `T_c = 0` >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, canonicalize, bsgs_direct_product >>> from sympy.combinatorics import Permutation @@ -692,7 +715,7 @@ def canonicalize(g, dummies, msym, *v): `T_c = -A^{d0 d1} * B_{d0}{}^{d2} * B_{d1 d2}` - `can = [0,2,1,4,3,5,7,6]` + can = [0,2,1,4,3,5,7,6] >>> t1 = (base2a, gens2a, 2, 1) >>> canonicalize(g, range(6), 0, t0, t1) @@ -733,10 +756,11 @@ def canonicalize(g, dummies, msym, *v): num_types = 1 else: num_types = len(msym) - if not all(msymx in [0,1,None] for msymx in msym): + if not all(msymx in [0, 1, None] for msymx in msym): raise ValueError('msym entries must be 0, 1 or None') if len(dummies) != num_types: - raise ValueError('dummies and msym must have the same number of elements') + raise ValueError( + 'dummies and msym must have the same number of elements') size = g.size num_tensors = 0 v1 = [] @@ -751,10 +775,10 @@ def canonicalize(g, dummies, msym, *v): can = canonicalize_naive(g, dummies, msym, *v) return can base_i, gens_i = mbsgs - v1.append((base_i, gens_i, [[]]*n_i, sym_i)) + v1.append((base_i, gens_i, [[]] * n_i, sym_i)) num_tensors += n_i - if num_types == 1: + if num_types == 1 and not isinstance(msym, list): dummies = [dummies] msym = [msym] flat_dummies = [] @@ -767,8 +791,9 @@ def canonicalize(g, dummies, msym, *v): # slot symmetry of the tensor size1, sbase, sgens = gens_products(*v1) if size != size1: - raise ValueError('g has size %d, generators have size %d' %(size, size1)) - free = [i for i in range(size-2) if i not in flat_dummies] + raise ValueError( + 'g has size %d, generators have size %d' % (size, size1)) + free = [i for i in range(size - 2) if i not in flat_dummies] num_free = len(free) # g1 minimal tensor under slot symmetry @@ -776,7 +801,7 @@ def canonicalize(g, dummies, msym, *v): if not flat_dummies: return g1 # save the sign of g1 - sign = 0 if g1[-1] == size-1 else 1 + sign = 0 if g1[-1] == size - 1 else 1 # the free indices are kept fixed. # Determine free_i, the list of slots of tensors which are fixed @@ -817,7 +842,8 @@ def canonicalize(g, dummies, msym, *v): dummies_red = [[x - num_free for x in y] for y in dummies] transv_red = get_transversals(sbase_red, sgens_red) g1_red = _af_new(g1_red) - g2 = double_coset_can_rep(dummies_red, msym, sbase_red, sgens_red, transv_red, g1_red) + g2 = double_coset_can_rep( + dummies_red, msym, sbase_red, sgens_red, transv_red, g1_red) if g2 == 0: return 0 # lift to the case with the free indices @@ -831,6 +857,7 @@ def perm_af_direct_product(gens1, gens2, signed=True): Examples ======== + >>> from sympy.combinatorics.tensor_can import perm_af_direct_product >>> gens1 = [[1,0,2,3], [0,1,3,2]] >>> gens2 = [[1,0]] @@ -849,7 +876,8 @@ def perm_af_direct_product(gens1, gens2, signed=True): start = list(range(n1)) end = list(range(n1, n1 + n2)) if signed: - gens1 = [gen[:-2] + end + [gen[-2]+n2, gen[-1]+n2] for gen in gens1] + gens1 = [gen[:-2] + end + [gen[-2] + n2, gen[-1] + n2] + for gen in gens1] gens2 = [start + [x + n1 for x in gen] for gen in gens2] else: gens1 = [gen + end for gen in gens1] @@ -859,6 +887,7 @@ def perm_af_direct_product(gens1, gens2, signed=True): return res + def bsgs_direct_product(base1, gens1, base2, gens2, signed=True): """ direct product of two BSGS @@ -873,6 +902,7 @@ def bsgs_direct_product(base1, gens1, base2, gens2, signed=True): Examples ======== + >>> from sympy.combinatorics import Permutation >>> from sympy.combinatorics.tensor_can import (get_symmetric_group_sgs, bsgs_direct_product) >>> Permutation.print_cyclic = True @@ -895,13 +925,19 @@ def bsgs_direct_product(base1, gens1, base2, gens2, signed=True): gens = [id_af] return base, [_af_new(h) for h in gens] -def get_symmetric_group_sgs(n, sym=0): + +def get_symmetric_group_sgs(n, antisym=False): """ - Return base, gens of the minimal BSGS for (anti)symmetric group - with n elements + Return base, gens of the minimal BSGS for (anti)symmetric tensor + + ``n`` rank of the tensor + + ``antisym = False`` symmetric tensor + ``antisym = True`` antisymmetric tensor Examples ======== + >>> from sympy.combinatorics import Permutation >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs >>> Permutation.print_cyclic = True @@ -910,16 +946,17 @@ def get_symmetric_group_sgs(n, sym=0): """ if n == 1: return [], [_af_new(range(3))] - gens = [Permutation(n-1)(i, i+1)._array_form for i in range(n-1)] - if sym == 0: - gens = [x + [n, n+1] for x in gens] + gens = [Permutation(n - 1)(i, i + 1)._array_form for i in range(n - 1)] + if antisym == 0: + gens = [x + [n, n + 1] for x in gens] else: - gens = [x + [n+1, n] for x in gens] - base = range(n-1) + gens = [x + [n + 1, n] for x in gens] + base = range(n - 1) return base, [_af_new(h) for h in gens] -riemann_bsgs = [0, 2], [Permutation(0,1)(4,5), Permutation(2,3)(4,5), \ - Permutation(5)(0,2)(1,3)] +riemann_bsgs = [0, 2], [Permutation(0, 1)(4, 5), Permutation(2, 3)(4, 5), + Permutation(5)(0, 2)(1, 3)] + def get_transversals(base, gens): """ @@ -927,10 +964,10 @@ def get_transversals(base, gens): """ if not base: return [] - stabs = _distribute_gens_by_base(base, gens) + stabs = _distribute_gens_by_base(base, gens) orbits, transversals = _orbits_transversals_from_bsgs(base, stabs) - transversals = [dict((x, h._array_form) for x, h in y.items()) for y in \ - transversals] + transversals = [dict((x, h._array_form) for x, h in y.items()) for y in + transversals] return transversals @@ -942,6 +979,7 @@ def _is_minimal_bsgs(base, gens): Examples ======== + >>> from sympy.combinatorics import Permutation >>> from sympy.combinatorics.tensor_can import riemann_bsgs, _is_minimal_bsgs >>> _is_minimal_bsgs(*riemann_bsgs) @@ -959,6 +997,7 @@ def _is_minimal_bsgs(base, gens): sgs1 = [h for h in sgs1 if h._array_form[i] == i] return base1 == base + def get_minimal_bsgs(base, gens): """ Compute a minimal GSGS @@ -987,6 +1026,7 @@ def get_minimal_bsgs(base, gens): return None return base, gens + def tensor_gens(base, gens, list_free_indices, sym=0): """ Returns size, res_base, res_gens BSGS for n tensors of the same type @@ -1036,7 +1076,7 @@ def _get_bsgs(G, base, gens, free_indices): if not base and list_free_indices.count([]) < 2: n = len(list_free_indices) size = gens[0].size - size = n*(gens[0].size - 2) + 2 + size = n * (gens[0].size - 2) + 2 return size, [], [_af_new(range(size))] # if any(list_free_indices) one needs to compute the pointwise @@ -1058,14 +1098,14 @@ def _get_bsgs(G, base, gens, free_indices): for i in range(1, len(list_free_indices)): base1, gens1 = _get_bsgs(G, base, gens, list_free_indices[i]) res_base, res_gens = bsgs_direct_product(res_base, res_gens, - base1, gens1, 1) + base1, gens1, 1) if not list_free_indices[i]: - no_free.append(range(size-2, size - 2 + num_indices)) + no_free.append(range(size - 2, size - 2 + num_indices)) size += num_indices nr = size - 2 res_gens = [h for h in res_gens if h._array_form != id_af] # if sym there are no commuting tensors stop here - if sym == None or not no_free: + if sym is None or not no_free: if not res_gens: res_gens = [_af_new(id_af)] return size, res_base, res_gens @@ -1085,7 +1125,7 @@ def _get_bsgs(G, base, gens, free_indices): a.extend(ind2) a.extend(ind1) base_comm.append(ind1[0]) - a.extend(range(ind2[-1]+1, nr)) + a.extend(range(ind2[-1] + 1, nr)) if sym == 0: a.extend([nr, nr + 1]) else: @@ -1130,7 +1170,7 @@ def gens_products(*v): for i in range(1, len(v)): size, base, gens = tensor_gens(*v[i]) res_base, res_gens = bsgs_direct_product(res_base, res_gens, base, - gens, 1) + gens, 1) res_size = res_gens[0].size id_af = range(res_size) res_gens = [h for h in res_gens if h != id_af] diff --git a/sympy/combinatorics/tests/test_tensor_can.py b/sympy/combinatorics/tests/test_tensor_can.py index 9f7912773bfd..8c2a4d5fc557 100644 --- a/sympy/combinatorics/tests/test_tensor_can.py +++ b/sympy/combinatorics/tests/test_tensor_can.py @@ -241,6 +241,18 @@ def test_no_metric_symmetry(): can = canonicalize(g, range(16), None, [[], [Permutation(range(4))], 8, 0]) assert can == [0,3,2,5,4,7,6,1,8,11,10,13,12,15,14,9,16,17] +def test_canonical_free(): + # t = A^{d0 a1}*A_d0^a0 + # ord = [a0,a1,d0,-d0]; g = [2,1,3,0,4,5]; dummies = [[2,3]] + # t_c = A_d0^a0*A^{d0 a1} + # can = [3,0, 2,1, 4,5] + base = [0] + gens = [Permutation(5)(0,2)(1,3)] + g = Permutation([2,1,3,0,4,5]) + num_free = 2 + dummies = [[2,3]] + can = canonicalize(g, dummies, [None], ([], [Permutation(3)], 2, 0)) + assert can == [3,0, 2,1, 4,5] def test_canonicalize1(): base1, gens1 = get_symmetric_group_sgs(1) @@ -305,12 +317,11 @@ def test_canonicalize1(): # g = [10,4,8, 0,7,9, 6,11,1, 2,3,5, 12,13] # T_c = A^{a0 d0 d1}*A^a1_d0^d2*A^{a2 a3 d3}*A_{d1 d2 d3} # can = [0,4,6, 1,5,8, 2,3,10, 7,9,11, 12,13] - base3, gens3 = get_symmetric_group_sgs(3) g = Permutation([10,4,8, 0,7,9, 6,11,1, 2,3,5, 12,13]) can = canonicalize(g, range(4,12), 0, (base3, gens3, 4, 0)) assert can == [0,4,6, 1,5,8, 2,3,10, 7,9,11, 12,13] - # A commuting symmetric, B anticommuting + # A commuting symmetric, B antisymmetric # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 # ord = [d0,-d0,d1,-d1,d2,-d2,d3,-d3] # g = [0,2,4,5,7,3,1,6,8,9] @@ -327,14 +338,15 @@ def test_canonicalize1(): # can = [0,2,4, 1,3,6, 5,7, 8,9] can = canonicalize(g, range(8), 0, (base3, gens3,2,1), (base2a,gens2a,1,0)) assert can == [0,2,4, 1,3,6, 5,7, 8,9] - # A anticommuting symmetric, B anticommuting, antisymmetric metric + # A anticommuting symmetric, B antisymmetric commuting, antisymmetric metric # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 # T_c = -A^{d0 d1 d2} * A_{d0 d1}^d3 * B_{d2 d3} # can = [0,2,4, 1,3,6, 5,7, 9,8] can = canonicalize(g, range(8), 1, (base3, gens3,2,1), (base2a,gens2a,1,0)) assert can == [0,2,4, 1,3,6, 5,7, 9,8] - # A anticommuting symmetric, B anticommuting, no metric symmetry + # A anticommuting symmetric, B anticommuting anticommuting, + # no metric symmetry # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 # T_c = A^{d0 d1 d2} * A_{d0 d1 d3} * B_d2^d3 # can = [0,2,4, 1,3,7, 5,6, 8,9] @@ -403,6 +415,8 @@ def test_riemann_invariants(): # R_d11^d1_d0^d5 * R^{d6 d4 d0}_d5 * R_{d7 d2 d8 d9} * # R_{d10 d3 d6 d4} * R^{d2 d7 d11}_d1 * R^{d8 d9 d3 d10} # ord: contravariant d_k ->2*k, covariant d_k -> 2*k+1 + # T_c = R^{d0 d1 d2 d3} * R_{d0 d1}^{d4 d5} * R_{d2 d3}^{d6 d7} * + # R_{d4 d5}^{d8 d9} * R_{d6 d7}^{d10 d11} * R_{d8 d9 d10 d11} g = Permutation([23,2,1,10,12,8,0,11,15,5,17,19,21,7,13,9,4,14,22,3,16,18,6,20,24,25]) can = canonicalize(g, range(24), 0, (baser, gensr, 6, 0)) assert can == [0,2,4,6,1,3,8,10,5,7,12,14,9,11,16,18,13,15,20,22,17,19,21,23,24,25] @@ -411,8 +425,6 @@ def test_riemann_invariants(): can = canonicalize(g, range(24), 0, ([2, 0], [Permutation([1,0,2,3,5,4]), Permutation([2,3,0,1,4,5])], 6, 0)) assert can == [0,2,4,6,1,3,8,10,5,7,12,14,9,11,16,18,13,15,20,22,17,19,21,23,24,25] - #R^{d0 d1 d2 d3} * R_{d0 d1}^{d4 d5} * R_{d2 d3}^{d6 d7} * - # R_{d4 d5}^{d8 d9} * R_{d6 d7}^{d10 d11} * R_{d8 d9 d10 d11} g = Permutation([0,2,5,7,4,6,9,11,8,10,13,15,12,14,17,19,16,18,21,23,20,22,25,27,24,26,29,31,28,30,33,35,32,34,37,39,36,38,1,3,40,41]) can = canonicalize(g, range(40), 0, (baser, gensr, 10, 0)) assert can == [0,2,4,6,1,3,8,10,5,7,12,14,9,11,16,18,13,15,20,22,17,19,24,26,21,23,28,30,25,27,32,34,29,31,36,38,33,35,37,39,40,41] @@ -470,6 +482,8 @@ def test_riemann_products(): # R^{d2 a0 a2 d0} * R^d1_d2^{a1 a3} * R^{a4 a5}_{d0 d1} # ord = [a0,a1,a2,a3,a4,a5,d0,-d0,d1,-d1,d2,-d2] # 0 1 2 3 4 5 6 7 8 9 10 11 + # can = [0, 6, 2, 8, 1, 3, 7, 10, 4, 5, 9, 11, 12, 13] + # T_c = R^{a0 d0 a2 d1}*R^{a1 a3}_d0^d2*R^{a4 a5}_{d1 d2} g = Permutation([10,0,2,6,8,11,1,3,4,5,7,9,12,13]) can = canonicalize(g, range(6,12), 0, (baser, gensr, 3, 0)) assert can == [0, 6, 2, 8, 1, 3, 7, 10, 4, 5, 9, 11, 12, 13] diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py index df8584991607..154b3e446f4c 100644 --- a/sympy/core/tests/test_args.py +++ b/sympy/core/tests/test_args.py @@ -2447,6 +2447,67 @@ def test_sympy__tensor__indexed__IndexedBase(): assert _test_args(IndexedBase('A', 1)) assert _test_args(IndexedBase('A')[0, 1]) +@XFAIL +def test_sympy__tensor__tensor__TensorIndexType(): + from sympy.tensor.tensor import TensorIndexType + from sympy import Symbol + assert _test_args(TensorIndexType('Lorentz', metric=False)) + +@XFAIL +def test_sympy__tensor__tensor__TensorSymmetry(): + from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorSymmetry, TensorType, get_symmetric_group_sgs + assert _test_args(TensorSymmetry(get_symmetric_group_sgs(2))) + + +@XFAIL +def test_sympy__tensor__tensor__TensorType(): + from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, get_symmetric_group_sgs, TensorType + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + sym = TensorSymmetry(get_symmetric_group_sgs(1)) + assert _test_args(TensorType([Lorentz], sym)) + +@XFAIL +def test_sympy__tensor__tensor__TensorHead(): + from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, TensorType, get_symmetric_group_sgs, TensorHead + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + sym = TensorSymmetry(get_symmetric_group_sgs(1)) + S1 = TensorType([Lorentz], sym) + assert _test_args(TensorHead('p', S1, 0)) + +@XFAIL +def test_sympy__tensor__tensor__TensorIndex(): + from sympy.tensor.tensor import TensorIndexType, TensorIndex, TensorSymmetry, TensorType, get_symmetric_group_sgs + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + assert _test_args(TensorIndex('i', Lorentz)) + +@SKIP("abstract class") +def test_sympy__tensor__tensor__TensExpr(): + pass + +def test_sympy__tensor__tensor__TensAdd(): + from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, TensorType, get_symmetric_group_sgs, tensor_indices, TensAdd + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + a, b = tensor_indices('a,b', Lorentz) + sym = TensorSymmetry(get_symmetric_group_sgs(1)) + S1 = TensorType([Lorentz], sym) + p, q = S1('p,q') + t1 = p(a) + t2 = q(a) + assert _test_args(TensAdd(t1, t2)) + + +def test_sympy__tensor__tensor__TensMul(): + from sympy.core import S + from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, TensorType, get_symmetric_group_sgs, tensor_indices, TensMul + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + a, b = tensor_indices('a,b', Lorentz) + sym = TensorSymmetry(get_symmetric_group_sgs(1)) + S1 = TensorType([Lorentz], sym) + p = S1('p') + free, dum = TensMul.from_indices(a) + assert _test_args(TensMul(S.One, [p], free, dum)) + + @XFAIL def test_as_coeff_add(): diff --git a/sympy/printing/str.py b/sympy/printing/str.py index eea89d07479a..f52a95ff6e3c 100644 --- a/sympy/printing/str.py +++ b/sympy/printing/str.py @@ -335,6 +335,18 @@ def _print_Permutation(self, expr): use = trim return 'Permutation(%s)' % use + def _print_TensorIndex(self, expr): + return expr._pretty() + + def _print_TensorHead(self, expr): + return expr._pretty() + + def _print_TensMul(self, expr): + return expr._pretty() + + def _print_TensAdd(self, expr): + return expr._pretty() + def _print_PermutationGroup(self, expr): p = [' %s' % str(a) for a in expr.args] return 'PermutationGroup([\n%s])' % ',\n'.join(p) diff --git a/sympy/tensor/dgamma_matr.py b/sympy/tensor/dgamma_matr.py new file mode 100644 index 000000000000..52a6b1ccfea4 --- /dev/null +++ b/sympy/tensor/dgamma_matr.py @@ -0,0 +1,534 @@ +from sympy import Symbol, S, I +from sympy.combinatorics import Permutation +from sympy.tensor.tensor import (TensExpr, TensorIndexType, tensor_indices, + TensorSymmetry, get_symmetric_group_sgs, TensorType, tensor_mul, TensMul, + TensAdd, TensorHead, tensorlist_contract_metric) + + +class GammaMatrices(object): + """ + Gamma matrices in dimensional regularization + with G5 naively anticommuting with all gamma matrices (NDR) + + NDR is inconsistent but widely used. + Notice in particular that the cyclic property of the + trace in presence of G5 does not hold. + + The original dimensional regularization scheme by 't Hooft and Veltman + has `G5` anticommuting only with `G(m)` for `m` in 0,1,2,3, + and commuting for `m > 3`. It is consistent, and it is used to + compute the chiral anomaly; however it breaks gauge invariance + involving G5 even when there are no anomalies, which complicates + a lot maintaining the Ward (or Slavnov-Taylor) identities. + + TODO: introduce the option of using the 't Hooft and Veltman scheme. + """ + + def __init__(self, typ, gctr=4, g5c=I): + """ + typ TensorIndexType (Lorentz or Eulidean) + + gct3 `tr(G(m)*G(n)) = gctr*g(m,n)` + in D dimensions, when D is an integer, the gamma + matrices are represented with matrices or rank `2**(D//2)`, + so `gctr = 2**(D//2)`; + in dimensional regularization one uses usually + `gctr = 4*g(m,n)`, so we choose `gctr=4` as default. + + g5c `tr(G(m0)*G(m1)*G(m2)*G(m3)) = gctr*g5c*epsilon(m0,m1,m2,m3)` + with Lorentz signature `g5c = I`, which is the default; + in Euclidean space `g5c = 1` + + The attribute `Gsymbol` can be used to set commutation relations + with other tensors using ``TensorManager``. + """ + self.g = typ.metric + if not typ.dim: + raise ValueError('Dimension not assigned') + self.D = typ.dim + self.typ = typ + sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + S1 = TensorType([typ], sym1) + self.Gsymbol = Symbol('Gsymbol') + self.G = S1('G', self.Gsymbol) + sym0 = TensorSymmetry(([], [Permutation(1)])) + S0 = TensorType([], sym0) + self.Gamma5 = S0('G5', self.Gsymbol) + self.G5 = TensMul(S.One, [self.Gamma5], [], []) + self.g5c = g5c + self.epsilon = typ.epsilon + self.eps_dim = typ.eps_dim + self.gctr = gctr + + def G5_to_right(self, t): + """ + move G5 to the right + """ + if not isinstance(t, TensExpr): + return t + if isinstance(t, TensAdd): + a = [self.G5_to_right(x) for x in t.args] + return TensAdd(*a) + components = t._components + if not components: + return t + ncomps = len(components) + G = self.G + Gamma5 = self.Gamma5 + G5 = self.G5 + # [g5,g,g,g5,g,g] + numG = 0 + vposG5 = [] + for i in range(ncomps): + if components[i] == Gamma5: + vposG5.append(i) + for j in range(i + 1, ncomps): + if components[j] == G: + numG += 1 + ct = t._coeff if 0 in vposG5 else S.One + a = t.split() + a1 = [] + for i in range(ncomps): + if not i in vposG5: + a1.append(a[i]) + for i in range(ncomps - 1, -1, -1): + if components[i] == G: + break + if len(vposG5) % 2: + a1 = a1[:i + 1] + [G5] + a1[i + 1:ncomps] + else: + a1 = a1[:i + 1] + a1[i + 1:ncomps] + t1 = ((-1)**numG*ct)*tensor_mul(*a1) + return t1 + + def match1_gamma(self, t, n): + #t = t.sorted_components() + components = t._components + ncomps = len(components) + G = self.G + for i in range(ncomps): + if not components[i] == G: + continue + for j in range(i + 1, ncomps): + if components[j] == G and abs(j - i) == n + 1 and \ + (0, 0, i, j) in t._dum: + return i, j + + + def rule1_gamma(self, t, n, r=None): + """ + simplify products of gamma matrices `G(m)*G(m_1)...*G(m_n)*G(-m)` + + t Gamma matrix monomial + n apply rule for `G(m)*G(m_1)...*G(m_n)*G(-m)` + + Examples + ======== + + >>> from sympy import Symbol, S + >>> from sympy.tensor.tensor import (TensorIndexType, tensor_indices,\ + TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul) + >>> from sympy.tensor.dgamma_matr import GammaMatrices + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + >>> S1 = TensorType([Lorentz], sym1) + >>> GM = GammaMatrices(Lorentz) + >>> G = GM.G + >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz) + >>> GM.gamma_trace(G(m0)*G(m1)*G(-m0)*G(m3)) + (-4*D + 8)*metric(m1, m3) + >>> t = G(m1)*G(m0)*G(-m1)*G(-m0) + >>> GM.rule1_gamma(t, 1) + (-D + 2)*G(L_0)*G(-L_0) + """ + if not r: + r = self.match1_gamma(t, n) + if not r: + return t + i, j = r + a = t.split() + za = zip(a, range(len(a))) + tc = t._coeff if i == 0 else S.One + a1 = [x for x, y in za if y not in r] + D = self.D + if n == 0: + # G(m)*G(-m) = d + t1 = tensor_mul(*a1) + t2 = TensMul(D*t1._coeff*tc, t1._components, t1._free, t1._dum) + return t2 + if n == 1: + # G(m)*G(n)*G(-m) = (2 - d)*G(n) + t1 = tensor_mul(*a1) + t2 = TensMul((2-D)*t1._coeff*tc, t1._components, t1._free, t1._dum) + return t2 + if n == 2: + # G(m)*G(n1)*G(n2)*G(-m) = (d-4)*G(n1)*G(n2) + 4*delta(n1, n2) + t1 = ((D-4)*tc)*tensor_mul(*a1) + a2 = a[:i] + a[j + 1:] + ind1 = a[r[0] + 1]._free[0][0] + ind2 = a[r[1] - 1]._free[0][0] + if a2: + a2 = tensorlist_contract_metric(a2, self.g(ind1, ind2)) + t2 = (4*tc)*tensor_mul(*a2) + else: + t2 = (4*tc)*self.g(ind1, ind2) + return t1 + t2 + if n == 3: + # G(m)*G(n1)*G(n2)*G(n3)*G(-m) = (4-d)*G(n1)*G(n2)*G(n3) - \ + # 2*G(n3)*G(n2)*G(n1) + t1 = ((4-D)*tc)*tensor_mul(*a1) + a2a = a[i + 1:j] + a2a.reverse() + a2 = a[:i] + a2a + a[j + 1:] + t2 = (-2*tc)*tensor_mul(*a2[:]) + return t1 + t2 + if n > 3: + # G(m0)*G(m_1)*...*G(m_k)*G(-m0) = + # 2*G(m_k)*G(m_1)*...*G(m_(k-1)) - + # G(m0)*G(m_1)*...*G(m_(k-1))G(-m0)*G(m_k)) + a1 = a[:i] + [a[j-1]] + a[i+1:j-1] + a[j + 1:] + t1 = (2*tc)*tensor_mul(*a1) + a2 = a[:j-1] + [a[j], a[j-1]] + a[j + 1:] + t2 = -tensor_mul(*a2) + return t1 + t2 + + + def do_rule1_gamma(self, t, nmax=4, doall=False): + """ + simplify products of gamma matrices `G(m)*G(m_1)...*G(m_n)*G(-m)` + + `t` Gamma matrix expression + + nmax maximum number for which `rule1\_gamma(t, n)` is applied + + doall if true apply `rule1\_gamma` till the expression does not change + + Examples + ======== + + >>> from sympy import Symbol, S + >>> from sympy.tensor.tensor import (TensorIndexType, tensor_indices,\ + TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul) + >>> from sympy.tensor.dgamma_matr import GammaMatrices + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + >>> S1 = TensorType([Lorentz], sym1) + >>> GM = GammaMatrices(Lorentz) + >>> G = GM.G + >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz) + >>> t = G(m1)*G(m0)*G(-m1)*G(-m0) + >>> GM.do_rule1_gamma(t, doall=True) + D*(-D + 2) + """ + if not isinstance(t, TensExpr): + return t + if isinstance(t, TensMul): + for n in range(nmax + 1): + r = self.match1_gamma(t, n) + if not r: + continue + t1 = self.rule1_gamma(t, n, r) + t1 = t1.contract_metric(self.g, contract_all=True) + if doall: + if t1 != t: + t1 = self.do_rule1_gamma(t1, nmax, doall) + return t1 + return t + if isinstance(t, TensAdd): + a = [self.do_rule1_gamma(tx, nmax, doall) for tx in t.args] + t1 = TensAdd(*a) + return t1.contract_metric(self.g, contract_all=True) + + def match2_gamma(self, t, n): + components = t._components + ncomps = len(components) + G = self.G + # G(a)*G(b)*...*p(-a)*p(-b) + for i in range(ncomps - 1 - n): + if not components[i] == G: + continue + if not components[i + n + 1] == G: + continue + p_pos1 = p_pos2 = 0 + for dx in t._dum: + if dx[2] == i: + p_pos1 = dx[3] + if dx[2] == i + n + 1: + p_pos2 = dx[3] + if p_pos1 > 0 and p_pos2 > 0: + if components[p_pos1] == components[p_pos2] \ + and components[p_pos1].comm == 0: + return i, p_pos1, p_pos2 + return None + + def rule2_gamma(self, t, n): + """ + simplify `G(m0)*p(-m0)*G(i_1)...*G(i_n)*G(m1)*p(-m1)` + + `t` Gamma matrix monomial + `n` as in `G(m0)*p(-m0)*G(i_1)...*G(i_n)*G(m1)*p(-m1)` + + Examples + ======== + + >>> from sympy import Symbol, S + >>> from sympy.tensor.tensor import (TensorIndexType, tensor_indices,\ + TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul) + >>> from sympy.tensor.dgamma_matr import GammaMatrices + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + >>> S1 = TensorType([Lorentz], sym1) + >>> GM = GammaMatrices(Lorentz) + >>> G = GM.G + >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz) + >>> p = S1('p') + >>> ps = G(m0)*p(-m0) + >>> t = G(m0)*ps*ps*G(-m0) + >>> GM.rule2_gamma(t, 0) + G(L_0)*G(-L_0)*p(L_1)*p(-L_1) + """ + t = t.canon_bp() + r = self.match2_gamma(t, n) + if not r: + return t + i0, i1, i2 = r + if n == 0: + a = t.split() + ct = t._coeff if i0 == 0 else S.One + p = a[i1] + ind_p = p._free[0][0] + # G(m0)*G(m1)*p(-m0)*p(-m1) = p(m0)*p(-m0) + a1 = a[:i0] + a[i0 + 2:i1] + a[i1 + 1:i2] + a[i2 + 1:] + a2 = [p.substitute_indices((ind_p, -ind_p)), p] + a3 = a1 + a2 + t = tensor_mul(*a3)*ct + return t + elif n == 1: + # G(m0)*G(ind1)*G(m1)*p(-m0)*p(-m1) = + # 2*G(m0)*p(-m0)*p(ind1) - G(ind1)*p(m0)*p(-m0) + a = t.split() + ind1 = a[i0 + 1]._free[0][0] + ct = t._coeff if i0 == 0 else S.One + p = a[i1] + ind_p = p._free[0][0] + a1 = a[:i0] + [a[i0 + 1]] + a[i0 + 3:i1] + a[i1 + 1:i2] + a[i2 + 1:] + a2 = [p.substitute_indices((ind_p, -ind_p)), p] + a3 = a1 + a2 + t1 = tensor_mul(*a3)*(-ct) + args = [t1] + a1 = a[:i0 + 1] + a[i0 + 3:i2] + a[i2 + 1:] + p = a[i2] + ind_p = p._free[0][0] + a2 = [p.substitute_indices((ind_p, ind1))] + a3 = a1 + a2 + t2 = tensor_mul(*a3)*2 + t3 = t1 + t2 + return t3 + else: + raise NotImplementedError + + def do_rule2_gamma(self, t, nmax=1): + """ + simplify `G(m0)*p(-m0)*G(m1)*p(-m1)` to `p(m)*p(-m)` + + `t` Gamma matrix expression + + Examples + ======== + + >>> from sympy import Symbol, S + >>> from sympy.tensor.tensor import (TensorIndexType, tensor_indices,\ + TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul) + >>> from sympy.tensor.dgamma_matr import GammaMatrices + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + >>> S1 = TensorType([Lorentz], sym1) + >>> GM = GammaMatrices(Lorentz) + >>> G = GM.G + >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz) + >>> p = S1('p') + >>> ps = G(m0)*p(-m0) + >>> M = Symbol('M') + >>> t = G(m0)*(ps + M)*(ps - M)*G(-m0) + >>> GM.do_rule2_gamma(t) + (-M**2)*G(L_0)*G(-L_0) + G(L_0)*G(-L_0)*p(L_1)*p(-L_1) + """ + if not isinstance(t, TensExpr): + return t + if isinstance(t, TensMul): + for n in range(nmax + 1): + t = self.rule2_gamma(t, n) + if isinstance(t, TensAdd): + a = [self.do_rule2_gamma(tx) for tx in t.args] + t = TensAdd(*a) + return t + + + def gamma_trace(self, t): + """ + Trace of gamma matrices + + Examples + ======== + + >>> from sympy import Symbol, S + >>> from sympy.tensor.tensor import (TensorIndexType, tensor_indices,\ + TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul) + >>> from sympy.tensor.dgamma_matr import GammaMatrices + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + >>> S1 = TensorType([Lorentz], sym1) + >>> GM = GammaMatrices(Lorentz) + >>> gamma_trace = GM.gamma_trace + >>> G = GM.G + >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz) + >>> gamma_trace(G(m0)*G(m1)*G(-m0)*G(m3)) + (-4*D + 8)*metric(m1, m3) + + >>> p, q = S1('p,q') + >>> t = G(m0)*G(m1)*G(m2)*G(m3)*p(-m1)*q(-m3) + >>> gamma_trace(t) + -4*metric(m0, m2)*p(L_0)*q(-L_0) + 4*p(m0)*q(m2) + 4*p(m2)*q(m0) + """ + if not isinstance(t, TensExpr): + return self.gctr*t + t = self.G5_to_right(t) + if isinstance(t, TensAdd): + a = [self.gamma_trace(x) for x in t.args] + return TensAdd(*a) + components = t._components + ncomps = len(components) + G = self.G + G5 = self.G5 + Gamma5 = self.Gamma5 + g = self.g + if all(x != G for x in components): + if any(x == Gamma5 for x in components): + return TensMul(S.Zero, [], [], []) + return self.gctr*t + withG5 = False + t = t.canon_bp() + components = t._components + for i in range(ncomps): + if components[i] == G: + break + for j in range(i + 1, ncomps): + if not components[j] == G: + break + else: + j = ncomps + numG = j - i + if any(x == Gamma5 for x in components): + withG5 = True + if withG5: + if numG < 4 or numG % 2 != self.eps_dim % 2: + return TensMul(S.Zero, [], [], []) + if numG == 4: + a = t.split() + indices = [x._free[0][0] for x in a[i:j]] + ct = t._coeff if i == 0 else S.One + t1 = (self.gctr*ct*self.g5c)*self.epsilon(*indices) + a2 = a[:i] + a[j + 1:] + t2 = tensor_mul(*a2) + res = t1*t2 + res = res.canon_bp() + return res + if numG > 4: + prev = TensMul(S.Zero, [], [], []) + t2 = t + while True: + t2 = self.do_rule1_gamma(t2) + t2 = t2.contract_metric(g, contract_all=True) + t2 = self.do_rule2_gamma(t2) + if t2 == prev: + break + prev = t2 + if t == t2: + a = t.split() + args = [] + for k1 in range(i, j): + ind1 = a[k1]._free[0][0] + for k2 in range(k1 + 1, j): + ind2 = a[k2]._free[0][0] + sign = 1 if (k1 - k2) % 2 == 1 else -1 + a1 = a[:k1] + a[k1 + 1:k2] + a[k2 + 1:] + a1 = tensorlist_contract_metric(a1, g(ind1, ind2)) + ct = t._coeff if k1 == 0 else S.One + t1 = (sign*ct)*tensor_mul(*a1) + args.append(t1) + t2 = TensAdd(*args) + return self.gamma_trace(t2) + else: + return self.gamma_trace(t2) + else: + # without G5 + if numG % 2 == 1: + return TensMul(S.Zero, [], [], []) + if numG > 4: + prev = TensMul(S.Zero, [], [], []) + t2 = t + while True: + t2 = self.do_rule1_gamma(t2) + t2 = t2.contract_metric(g, contract_all=True) + t2 = self.do_rule2_gamma(t2) + if t2 == prev: + break + prev = t2 + if t == t2: + a = t.split() + # G are from i to j1 included; anticommute G(j1) through + ind1 = a[i]._free[0][0] + ind2 = a[i + 1]._free[0][0] + aa = a[:i] + a[i + 2:] + aa = tensorlist_contract_metric(aa, g(ind1, ind2)) + ct = t._coeff if i == 0 else S.One + t1 = ct*tensor_mul(*aa) + args = [t1] + sign = 1 + for k in range(i + 2, j): + sign = -sign + ind2 = a[k]._free[0][0] + aa = a[:i] + a[i + 1:k] + a[k + 1:] + aa = tensorlist_contract_metric(aa, g(ind1, ind2)) + t2 = sign*tensor_mul(*aa) + args.append(t2) + t3 = TensAdd(*args) + return self.gamma_trace(t3) + return self.gamma_trace(t2) + + a = t.split() + t1 = self.gamma_trace1(*a[i:j]) + a2 = a[:i] + a[j:] + t2 = tensor_mul(*a2) + ct = t._coeff if i == 0 else S.One + t3 = ct*t1*t2 + if not t3: + return t3 + t3 = t3.contract_metric(g, contract_all=True) + return t3 + + + def gamma_trace1(self, *a): + if not a: + return self.gctr + n = len(a) + if n%2 == 1: + return TensMul(S.Zero, [], [], []) + if n == 2: + ind0 = a[0]._free[0][0] + ind1 = a[1]._free[0][0] + return self.gctr*self.g(ind0, ind1) + if n == 4: + ind0 = a[0]._free[0][0] + ind1 = a[1]._free[0][0] + ind2 = a[2]._free[0][0] + ind3 = a[3]._free[0][0] + return self.gctr*(self.g(ind0, ind1)*self.g(ind2, ind3) - \ + self.g(ind0, ind2)*self.g(ind1, ind3) + self.g(ind0, ind3)*self.g(ind1, ind2)) + else: + raise NotImplementedError diff --git a/sympy/tensor/group_factors.py b/sympy/tensor/group_factors.py new file mode 100644 index 000000000000..79c24441848a --- /dev/null +++ b/sympy/tensor/group_factors.py @@ -0,0 +1,316 @@ +from sympy import I, S +from sympy.tensor.tensor import (TensorIndexType, tensor_indices, \ +TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul, tensorlist_contract_metric, tensor_mul, TensAdd, TensExpr, tensorsymmetry) + + + +def match_CijCij(t, th): + """ + match for two tensors components ``th`` with two or more indices contracted + """ + components = t._components + rdum = [tuple(sorted([cpos1, cpos2])) for ipos1, ipos2, cpos1, cpos2 \ + in t._dum if components[cpos1] == th and components[cpos2] == th] + rdum.sort() + vcount = [rdum.count(x) for x in rdum] + repeated = {} + for i in range(len(rdum)): + if vcount[i] > 1 and rdum[i] not in repeated: + repeated[rdum[i]] = vcount[i] + if not repeated: + return None + repeated = repeated.items() + res = max(repeated, key=lambda x: x[1]) + return res + +def match_C_triangle(t, th): + """ + match for three tensor components ``th`` with triangle of contractions + """ + components = t._components + rdum = [sorted([cpos1, cpos2]) for ipos1, ipos2, cpos1, cpos2 in \ + t._dum if components[cpos1] == th and components[cpos2] == th] + rdum.sort() + for cpos1, cpos2 in rdum: + for i in range(cpos1 + 1, len(rdum)): + if [cpos1, i] in rdum: + if [cpos2, i] in rdum: + return cpos1, cpos2, i + elif [i, cpos2] in rdum: + return cpos1, i, cpos2 + +def match_C_square(t, th): + """ + match for four tensor components ``th`` with square of contractions + + It is assumed that there are not smaller cycles + """ + components = t._components + rdum = [sorted([cpos1, cpos2]) for ipos1, ipos2, cpos1, cpos2 in \ + t._dum if components[cpos1] == th and components[cpos2] == th] + rdum.sort() + for i in range(len(rdum)): + cpos1, cpos2 = rdum[i] + for j in range(i + 1, len(rdum)): + cpos3, cpos4 = rdum[j] + if cpos3 in rdum[i] or cpos4 in rdum[i]: + continue + if [cpos1, cpos3] in rdum: + if [cpos2, cpos4] in rdum or [cpos4, cpos2] in rdum: + return cpos1, cpos2, cpos4, cpos3 + if [cpos1, cpos4] in rdum: + if [cpos2, cpos3] in rdum or [cpos3, cpos2] in rdum: + return cpos1, cpos2, cpos3, cpos4 + + +class SuNGroupFactors(object): + """ + SU(N) generators in any representation satisfy the + Lie algebra ``[T(i), T(j)] = I*C(i,j, -k)*T(k)`` + + ``tr(rep, T(i)*T(j)) = a_rep*g(i, j)`` where + ``a_rep`` depends on the representation. + + ``g(i, j)`` is the metric in the adjoint representation space. + + For ``N > 2`` there is no metric in the defining representation. + + Let ``T(i, -a, b)`` be the generator in the defining representation ``N`` + normalized with ``a_rep = 1`` + + Define ``theta(i_1, ..., i_n) = tr(T(i_1)*...*T^(i_n)) + + Multiplying the Lie algebra relation by ``T(-k)`` and tracing one gets + ``C_(i j k) = -I*theta(i, j, k) + I*theta(j, i, k)`` + + For ``SU(N)`` one has the relation + ``T(i,-d,a)*T(-i,-b,c) =`` + ``delta(-b,a)*delta(-d,c) - 1/N*delta(-d,a)*delta(-b,c)`` + + Using this relation one obtains easily relations among contracted + ``theta`` tensors. + """ + def __init__(self, N): + self.N = N + self.DSuN = TensorIndexType('DSuN', metric=None, dim=N, dummy_fmt='D') + self.ASuN = TensorIndexType('ASuN', metric=0, dim = N**2-1, dummy_fmt='A') + self.delta = self.DSuN.delta + self.g = self.ASuN.metric + symT = tensorsymmetry(*[[1]]*3) + ST = TensorType([self.ASuN, self.DSuN, self.DSuN], symT) + self.T = ST('T') + symA = tensorsymmetry([3]) + SC = TensorType([self.ASuN]*3, symA) + self.C = SC('C') + self._theta = [None]*3 + [TensorType([self.ASuN]*i, \ + tensorsymmetry(*[[1]]*i))('theta') for i in range(3, 6)] + i0, i1, i2 = tensor_indices('i0,i1,i2', self.ASuN) + self.tC = self.C(i0,i1,i2) + self.tCs = -I*self.theta(i0, i1, i2) + I*self.theta(i1, i0, i2) + + def theta(self, *indices): + n = len(indices) + if n == 1: + return S.Zero + if n == 2: + return self.g(*indices) + try: + return self._theta[n](*indices) + except IndexError: + for i in range(len(self._theta), n + 1): + self._theta.append(TensorType([self.ASuN]*i, tensorsymmetry(*[[1]]*i))('theta')) + return self._theta[n](*indices) + + def match_thetaii(self, t): + """ + match a ``theta`` with two indices contracted + + Returns the slot position of the contracted indices and the position + of the ``theta``. + """ + components = t._components + dum = t._dum + for ipos1, ipos2, cpos1, cpos2 in dum: + if cpos1 == cpos2 and \ + components[cpos1].name == 'theta' and \ + components[cpos1].types == [self.ASuN]: + return ipos1, ipos2, cpos1 + + def match_thetaithetai(self, t): + """ + match for two component tensors ``theta`` with + an index contracted with each other + + Return the positions of the contracted indices and the positions + of the two tensors. + """ + components = t._components + dum = t._dum + for ipos1, ipos2, cpos1, cpos2 in dum: + if cpos1 != cpos2 and \ + components[cpos1].name == 'theta' and \ + components[cpos2].name == 'theta' and \ + components[cpos1].types == components[cpos1].types == [self.ASuN]: + return ipos1, ipos2, cpos1, cpos2 + + + def rule_C_triangle(self, t): + """ + evaluate a triangle of C's + """ + if isinstance(t, TensAdd): + args = t.args + args = [self.rule_C_triangle(x) for x in args] + return TensAdd(*args) + r = match_C_triangle(t, self.C) + if not r: + return t + i0, i1, i2 = sorted(r) + a = t.split() + a1 = a[:i0] + a[i0 + 1: i1] + a[1 + 1:i2] + a[i2 + 1:] + t1 = tensor_mul(*a1) + t2 = tensor_mul(a[i0], a[i1], a[i2]) + t2 = self.rule_C2theta_all(t2) + t = t1*t2 + prev = S.Zero + while t != prev: + prev = t + t = self.rule_thetaithetai(t) + t = t.contract_metric(self.g, True) + t = self.rule_thetaii(t) + return t + + def rule_C_square(self, t): + """ + evaluate a square of C's + """ + if isinstance(t, TensAdd): + args = t.args + args = [self.rule_C_square(x) for x in args] + return TensAdd(*args) + r = match_C_square(t, self.C) + if not r: + return t + i0, i1, i2, i3 = sorted(r) + a = t.split() + a1 = a[:i0] + a[i0 + 1: i1] + a[1 + 1:i2] + a[i2 + 1:i3] + a[i3 + 1:] + t1 = tensor_mul(*a1) + t2 = tensor_mul(a[i0], a[i1], a[i2], a[i3]) + t2 = self.rule_C2theta_all(t2) + t = t1*t2 + prev = S.Zero + while t != prev: + prev = t + t = self.rule_thetaithetai(t) + t = t.contract_metric(self.g, True) + t = self.rule_thetaii(t) + return t + + def rule_thetaii(self, t): + """ + theta(k,i_0,..,i_m,-k,j_0,..,j_n) = + theta(i_0,..,i_m)*theta(j_0,..,j_n) - 1/N*theta(i_0,..,i_m,j_0,..,j_n) + + """ + if isinstance(t, TensAdd): + args = t.args + args = [self.rule_thetaii(x) for x in args] + return TensAdd(*args) + r = self.match_thetaii(t) + if not r: + return t + ipos1, ipos2, cpos1 = r + a = t.split() + ct = t._coeff if cpos1 == 0 else S.One + a1 = a[:cpos1] + a[cpos1 + 1:] + t1 = a[cpos1] + indices = t1.get_indices() + if ipos1 < ipos2: + indices = indices[ipos1:] + indices[:ipos1] + i = ipos2 - ipos1 + else: + indices = indices[ipos2:] + indices[:ipos2] + i = ipos1 - ipos2 + indices1 = indices[1:i] + indices2 = indices[i + 1:] + theta = self.theta + N = self.N + if indices1 == []: + t2 = (N - 1/N)*theta(*indices2) + elif indices2 == []: + t2 = (N - 1/N)*theta(*indices1) + else: + t2 = theta(*indices1)*theta(*indices2) - 1/N*theta(*(indices1 + indices2)) + t3 = tensor_mul(*a1) + return t2*t3*ct + + def rule_thetaithetai(self, t): + """ + theta(k,i_0,..,i_m)*theta(k,j_0,...,j_n) = + theta(i_0,...,i_m,j_0,...,j_n) - 1/N*theta(i_0,..,i_m)*theta(j_0,...,j_n) + """ + if isinstance(t, TensAdd): + args = t.args + args = [self.rule_thetaithetai(x) for x in args] + return TensAdd(*args) + + r = self.match_thetaithetai(t) + if not r: + return t + ipos1, ipos2, cpos1, cpos2 = r + a = t.split() + if cpos1 < cpos2: + a1 = a[:cpos1] + a[cpos1 + 1:cpos2] + a[cpos2 + 1:] + else: + a1 = a[:cpos2] + a[cpos2 + 1:cpos1] + a[cpos1 + 1:] + if cpos1 == 0: + ct = t._coeff + elif cpos2 == 0: + ct = t._coeff + else: + ct = S.One + t1 = a[cpos1] + t2 = a[cpos2] + indices1 = t1.get_indices() + indices2 = t2.get_indices() + indices1 = indices1[ipos1 + 1:] + indices1[:ipos1] + indices2 = indices2[ipos2 + 1:] + indices2[:ipos2] + theta = self.theta + N = self.N + t3 = theta(*(indices1 + indices2)) + t4 = theta(*indices1)*theta(*indices2) + t5 = theta(*(indices1 + indices2)) - 1/N*theta(*indices1)*theta(*indices2) + res = tensor_mul(*a1)*t5*ct + return res + + def rule_C2theta_all(self, t): + """ + convert all ``C`` to ``theta`` and apply simplifying rules. + + Examples + ======== + + >>> from sympy.tensor.group_factors import SuNGroupFactors + >>> from sympy import Symbol + >>> from sympy.tensor.tensor import tensor_indices + >>> N = Symbol('N') + >>> SN = SuNGroupFactors(N) + >>> C = SN.C + >>> i0,i1,i2,i3,i4,i5,i6 = tensor_indices('i0:7', SN.ASuN) + >>> t = C(i0,i1,i2)*C(-i0,-i2,i3) + >>> t = SN.rule_C2theta_all(t) + >>> SN.rule_C2theta_all(t) + (-2*N)*metric(i1, i3) + """ + + tC = self.tC + tCs = self.tCs + t = t.substitute_tensor(tC, tCs, True) + prev = S.Zero + while t != prev: + prev = t + t = t.substitute_tensor(tC, tCs, True) + t = self.rule_thetaithetai(t) + t = t.contract_metric(self.g, True) + t = self.rule_thetaii(t) + return t diff --git a/sympy/tensor/tensor.py b/sympy/tensor/tensor.py new file mode 100644 index 000000000000..0f54d3f97ddd --- /dev/null +++ b/sympy/tensor/tensor.py @@ -0,0 +1,2331 @@ +""" +This module defines tensors with abstract index notation. + +The abstract index notation has been first formalized by Penrose. + +Tensor indices are formal objects, with a tensor type; there is no +notion of index range, it is only possible to assign the dimension, +used to trace the Kronecker delta; the dimension can be a Symbol. + +The Einstein summation convention is used. +The covariant indices are indicated with a minus sign in front of the index. + +For instance the tensor ``t = p(a)*A(b,c)*q(-c)`` has the index ``c`` +contracted. + +A tensor expression ``t`` can be called; called with its +indices in sorted order it is equal to itself: +in the above example ``t(a, b) == t``; +one can call ``t`` with different indices; ``t(c, d) == p(c)*A(d,a)*q(-a)``. + +The contracted indices are dummy indices, internally they have no name, +the indices being represented by a graph-like structure. + +Tensors are put in canonical form using ``canon_bp``, which uses +the Butler-Portugal algorithm for canonicalization. + +If there is a (anti)symmetric metric, the indices can be raised and +lowered when the tensor is put in canonical form. +""" + + +from collections import defaultdict +from sympy.core import Basic, sympify, Add, Mul, S +from sympy.core.symbol import Symbol, symbols +from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, bsgs_direct_product, canonicalize, riemann_bsgs + +class _TensorManager(object): + """ + Class to manage tensor properties. + + Notes + ===== + + Tensors belong to tensor commutation groups; each group has a label + ``comm``; there are predefined labels: + + ``0`` tensors commuting with any other tensor + + ``1`` tensors anticommuting among themselves + + ``2`` tensors not commuting, apart with those with ``comm=0`` + + Other groups can be defined using ``set_comm``; tensors in those + groups commute with those with ``comm=0``; by default they + do not commute with any other group. + """ + def __init__(self): + self._comm_init() + + def _comm_init(self): + self._comm = [{} for i in range(3)] + for i in range(3): + self._comm[0][i] = 0 + self._comm[i][0] = 0 + self._comm[1][1] = 1 + self._comm[2][1] = None + self._comm[1][2] = None + self._comm_symbols2i = {0:0, 1:1, 2:2} + self._comm_i2symbol = {0:0, 1:1, 2:2} + + @property + def comm(self): + return self._comm + + def comm_symbols2i(self, i): + """ + get the commutation group number corresponding to ``i`` + + ``i`` can be a symbol or a number or a string + + If ``i`` is not already defined its commutation group number + is set. + """ + if i not in self._comm_symbols2i: + n = len(self._comm) + self._comm.append({}) + self._comm[n][0] = 0 + self._comm[0][n] = 0 + self._comm_symbols2i[i] = n + self._comm_i2symbol[n] = i + return n + return self._comm_symbols2i[i] + + def comm_i2symbol(self, i): + """ + Returns the symbol corresponding to the commutation group number. + """ + return self._comm_i2symbol[i] + + def set_comm(self, i, j, c): + """ + set the commutation parameter ``c`` for commutation groups ``i, j`` + + Parameters + ========== + + i, j : symbols representing commutation groups + + c : group commutation number + + Notes + ===== + + ``i, j`` can be symbols, strings or numbers, + apart from ``0, 1`` and ``2`` which are reserved respectively + for commuting, anticommuting tensors and tensors not commuting + with any other group apart with the commuting tensors. + For the remaining cases, use this method to set the commutation rules; + by default ``c=None``. + + The group commutation number ``c`` is assigned in corrispondence + to the group commutation symbols; it can be + + 0 commuting + + 1 anticommuting + + None no commutation property + + Examples + ======== + + ``G`` and ``GH`` do not commute with themselves and commute with + each other; A is commuting. + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead, TensorManager + >>> Lorentz = TensorIndexType('Lorentz') + >>> i0,i1,i2,i3,i4 = tensor_indices('i0:5', Lorentz) + >>> A = tensorhead('A', [Lorentz], [[1]]) + >>> G = tensorhead('G', [Lorentz], [[1]], 'Gcomm') + >>> GH = tensorhead('GH', [Lorentz], [[1]], 'GHcomm') + >>> TensorManager.set_comm('Gcomm', 'GHcomm', 0) + >>> (GH(i1)*G(i0)).canon_bp() + G(i0)*GH(i1) + >>> (G(i1)*G(i0)).canon_bp() + G(i1)*G(i0) + >>> (G(i1)*A(i0)).canon_bp() + A(i0)*G(i1) + """ + if c not in (0, 1, None): + raise ValueError('`c` can assume only the values 0, 1 or None') + + if i not in self._comm_symbols2i: + n = len(self._comm) + self._comm.append({}) + self._comm[n][0] = 0 + self._comm[0][n] = 0 + self._comm_symbols2i[i] = n + self._comm_i2symbol[n] = i + if j not in self._comm_symbols2i: + n = len(self._comm) + self._comm.append({}) + self._comm[0][n] = 0 + self._comm[n][0] = 0 + self._comm_symbols2i[j] = n + self._comm_i2symbol[n] = j + ni = self._comm_symbols2i[i] + nj = self._comm_symbols2i[j] + self._comm[ni][nj] = c + self._comm[nj][ni] = c + + def set_comms(self, *args): + """ + set the commutation group numbers ``c`` for symbols ``i, j`` + + Parameters + ========== + + args : sequence of ``(i, j, c)`` + """ + for i, j, c in args: + self.set_comm(i, j, c) + + def get_comm(self, i, j): + """ + Return the commutation parameter for commutation group numbers ``i, j`` + + see ``_TensorManager.set_comm`` + """ + return self._comm[i].get(j, 0 if i == 0 or j == 0 else None) + + def clear(self): + """ + Clear the TensorManager. + """ + self._comm_init() + + +TensorManager = _TensorManager() + +class TensorIndexType(Basic): + """ + A TensorIndexType is characterized by its name and its metric. + + Parameters + ========== + + name : name of the tensor type + + metric : metric symmetry or metric object or ``None`` + + + dim : dimension, it can be a symbol or an integer or ``None`` + + eps_dim : dimension of the epsilon tensor + + dummy_fmt : name of the head of dummy indices + + Attributes + ========== + + ``name`` + ``metric_name`` : it is 'metric' or metric.name + ``metric_antisym`` + ``metric`` : the metric tensor + ``delta`` : ``Kronecker delta`` + ``epsilon`` : the ``Levi-Civita epsilon`` tensor + ``dim`` + ``dim_eps`` + ``dummy_fmt`` + + Notes + ===== + + The ``metric`` parameter can be: + ``metric = False`` symmetric metric (in Riemannian geometry) + + ``metric = True`` antisymmetric metric (for spinor calculus) + + ``metric = None`` there is no metric + + ``metric`` can be an object having ``name`` and ``antisym`` attributes. + + + If there is a metric the metric is used to raise and lower indices. + + In the case of antisymmetric metric, the following raising and + lowering conventions will be adopted: + + ``psi(a) = g(a, b)*psi(-b); chi(-a) = chi(b)*g(-b, -a)`` + + ``g(-a, b) = delta(-a, b); g(b, -a) = -delta(a, -b)`` + + where ``delta(-a, b) = delta(b, -a)`` is the ``Kronecker delta`` + (see ``TensorIndex`` for the conventions on indices). + + If there is no metric it is not possible to raise or lower indices; + e.g. the index of the defining representation of ``SU(N)`` + is 'covariant' and the conjugate representation is + 'contravariant'; for ``N > 2`` they are linearly independent. + + ``eps_dim`` is by default equal to ``dim``, if the latter is an integer; + else it can be assigned (for use in naive dimensional regularization); + if ``eps_dim`` is not an integer ``epsilon`` is ``None``. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> Lorentz.metric + metric(Lorentz,Lorentz) + """ + + def __new__(cls, name, metric=False, dim=None, eps_dim = None, + dummy_fmt=None): + obj = Basic.__new__(cls, name, metric) + obj._name = name + if not dummy_fmt: + obj._dummy_fmt = '%s_%%d' % obj.name + else: + obj._dummy_fmt = '%s_%%d' % dummy_fmt + if metric is None: + obj.metric_antisym = None + obj.metric = None + else: + if metric in (True, False, 0, 1): + metric_name = 'metric' + obj.metric_antisym = metric + else: + metric_name = metric.name + obj.metric_antisym = metric.antisym + sym2 = TensorSymmetry(get_symmetric_group_sgs(2, obj.metric_antisym)) + S2 = TensorType([obj]*2, sym2) + obj.metric = S2(metric_name) + + obj._dim = dim + obj._delta = obj.get_kronecker_delta() + obj._eps_dim = eps_dim if eps_dim else dim + obj._epsilon = obj.get_epsilon() + return obj + + @property + def name(self): + return self._name + + @property + def dim(self): + return self._dim + + @property + def delta(self): + return self._delta + + @property + def eps_dim(self): + return self._eps_dim + + @property + def epsilon(self): + return self._epsilon + + @property + def dummy_fmt(self): + return self._dummy_fmt + + def get_kronecker_delta(self): + sym2 = TensorSymmetry(get_symmetric_group_sgs(2)) + S2 = TensorType([self]*2, sym2) + delta = S2('KD') + return delta + + def get_epsilon(self): + if not isinstance(self._eps_dim, int): + return None + sym = TensorSymmetry(get_symmetric_group_sgs(self._eps_dim, 1)) + Sdim = TensorType([self]*self._eps_dim, sym) + epsilon = Sdim('Eps') + return epsilon + + def __lt__(self, other): + return self.name < other.name + + def __str__(self): + return self.name + + __repr__ = __str__ + +class TensorIndex(Basic): + """ + Represents an abstract tensor index. + + Parameters + ========== + + name : name of the index + tensortype : ``TensorIndexType`` of the index + is_up : flag for contravariant index + + Attributes + ========== + + ``name`` + ``tensortype`` + ``is_up`` + + Notes + ===== + + Tensor indices are contracted with the Einstein summation convention. + + An index can be in contravariant or in covariant form; in the latter + case it is represented prepending a ``-`` to the index name. + + Dummy indices have a name with head given by ``tensortype.dummy_fmt`` + + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, TensorIndex, TensorSymmetry, TensorType, get_symmetric_group_sgs + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> i = TensorIndex('i', Lorentz); i + i + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) + >>> S1 = TensorType([Lorentz], sym1) + >>> A, B = S1('A,B') + >>> A(i)*B(-i) + A(L_0)*B(-L_0) + """ + def __new__(cls, name, tensortype, is_up=True): + + obj = Basic.__new__(cls, name, tensortype, is_up) + obj._name = name + obj._tensortype = tensortype + obj._is_up = is_up + return obj + + @property + def name(self): + return self._name + + @property + def tensortype(self): + return self._tensortype + + @property + def is_up(self): + return self._is_up + + def _pretty(self): + s = self._name + if not self._is_up: + s = '-%s' % s + return s + + def __lt__(self, other): + return (self._tensortype, self._name) < (other._tensortype, other._name) + + def __neg__(self): + t1 = TensorIndex(self._name, self._tensortype, + (not self._is_up)) + return t1 + +def tensor_indices(s, typ): + """ + Returns list of tensor indices given their names and their types + + Parameters + ========== + + s : string of comma separated names of indices + + typ : list of ``TensorIndexType`` of the indices + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz) + """ + if isinstance(s, str): + a = [x.name for x in symbols(s, seq=True)] + else: + raise ValueError('expecting a string') + + return [TensorIndex(i, typ) for i in a] + + +class TensorSymmetry(Basic): + """ + Symmetry of a tensor + + Parameters + ========== + + bsgs : tuple ``(base, sgs)`` BSGS of the symmetry of the tensor + + Attributes + ========== + + ``base`` : base of the BSGS + ``generators`` : generators of the BSGS + ``rank`` : rank of the tensor + + Notes + ===== + + A tensor can have an arbitrary symmetry provided by its BSGS. + + See Also + ======== + + sympy.combinatorics.tensor_can.get_symmetric_group_sgs + + Examples + ======== + + Define a symmetric tensor + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorSymmetry, TensorType, get_symmetric_group_sgs + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> sym2 = TensorSymmetry(get_symmetric_group_sgs(2)) + >>> S2 = TensorType([Lorentz]*2, sym2) + >>> V = S2('V') + """ + def __new__(cls, bsgs, **kw_args): + base, generators = bsgs + obj = Basic.__new__(cls, base, generators, **kw_args) + return obj + + @property + def base(self): + return self.args[0] + + @property + def generators(self): + return self.args[1] + + @property + def rank(self): + return self.args[1][0].size - 2 + + def _hashable_content(self): + r = (tuple(self.base), tuple(self.generators)) + return r + + +def tensorsymmetry(*args): + """ + Return a ``TensorSymmetry`` object. + + One can represent a tensor with any slot symmetry group using a BSGS + ``args`` can be a BSGS + ``args[0]`` base + ``args[1]`` sgs + + Usually tensors are in (direct products of) irreducible representations + of the symmetric group; + ``args`` can be a list of lists representing Young tableaux + ``[[1]]`` vector + ``[[1]*n]`` symmetric tensor of rank ``n`` + ``[[n]`` antisymmetric tensor of rank ``n`` + ``[[2, 2]]`` slot symmetry of the Riemann tensor + ``[[1],[1]]`` vector*vector + + + Examples + ======== + + Symmetric tensor using a Young tableau + + >>> from sympy.tensor.tensor import TensorIndexType, TensorType, tensorsymmetry + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> sym2 = tensorsymmetry([1, 1]) + >>> S2 = TensorType([Lorentz]*2, sym2) + >>> V = S2('V') + + Symmetric tensor using a BSGS + >>> from sympy.tensor.tensor import TensorSymmetry, get_symmetric_group_sgs + >>> sym2 = tensorsymmetry(*get_symmetric_group_sgs(2)) + >>> S2 = TensorType([Lorentz]*2, sym2) + >>> V = S2('V') + """ + from sympy.combinatorics import Permutation + def tableau2bsgs(a): + if len(a) == 1: + # antisymmetric vector + n = a[0] + bsgs = get_symmetric_group_sgs(n, 1) + else: + if all(x == 1 for x in a): + # symmetric vector + n = len(a) + bsgs = get_symmetric_group_sgs(n) + elif a == [2, 2]: + bsgs = riemann_bsgs + else: + raise NotImplementedError + return bsgs + + if not args: + return TensorSymmetry([[], [Permutation(1)]]) + if len(args) == 2 and isinstance(args[1][0], Permutation): + return TensorSymmetry(args) + base, sgs = tableau2bsgs(args[0]) + for a in args[1:]: + basex, sgsx = tableau2bsgs(a) + base, sgs = bsgs_direct_product(base, sgs, basex, sgsx) + return TensorSymmetry((base, sgs)) + + +class TensorType(Basic): + """ + Class of tensor types. + + Parameters + ========== + + index_types : list of ``TensorIndexType`` of the tensor indices + symmetry : ``TensorSymmetry`` of the tensor + + Attributes + ========== + + ``index_types`` + ``symmetry`` + ``types`` : list of ``TensorIndexType`` without repetitions + + Examples + ======== + + Define a symmetric tensor + + >>> from sympy.tensor.tensor import TensorIndexType, tensorsymmetry, TensorType + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> sym2 = tensorsymmetry([1, 1]) + >>> S2 = TensorType([Lorentz]*2, sym2) + >>> V = S2('V') + """ + is_commutative = False + + def __new__(cls, index_types, symmetry, **kw_args): + assert symmetry.rank == len(index_types) + obj = Basic.__new__(cls, index_types, symmetry, **kw_args) + return obj + + @property + def index_types(self): + return self.args[0] + + @property + def symmetry(self): + return self.args[1] + + @property + def types(self): + return sorted(set(self.index_types), key=lambda x: x.name) + + def __str__(self): + return 'TensorType(%s)' %([str(x) for x in self.index_types]) + + def __call__(self, s, comm=0): + """ + Return a TensorHead object or a list of TensorHead objects. + + ``s`` name or string of names + + ``comm``: commutation group number + see ``_TensorManager.set_comm`` + + Examples + ======== + + Define symmetric tensors ``V``, ``W`` and ``G``, respectively + commuting, anticommuting and with no commutation symmetry + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorsymmetry, TensorType, canon_bp + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a, b = tensor_indices('a,b', Lorentz) + >>> sym2 = tensorsymmetry([1]*2) + >>> S2 = TensorType([Lorentz]*2, sym2) + >>> V = S2('V') + >>> W = S2('W', 1) + >>> G = S2('G', 2) + >>> canon_bp(V(a, b)*V(-b, -a)) + V(L_0, L_1)*V(-L_0, -L_1) + >>> canon_bp(W(a, b)*W(-b, -a)) + 0 + """ + if isinstance(s, str): + names = [x.name for x in symbols(s, seq=True)] + else: + raise ValueError('expecting a string') + if len(names) == 1: + return TensorHead(names[0], self, comm) + else: + return [TensorHead(name, self, comm) for name in names] + +def tensorhead(name, typ, sym, comm=0): + """ + Function generating tensorhead(s). + + Parameters + ========== + + name : name or sequence of names (as in ``symbol``) + + typ : index types + + sym : same as ``*args`` in ``tensorsymmetry`` + + comm : commutation group number + see ``_TensorManager.set_comm`` + + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a, b = tensor_indices('a,b', Lorentz) + >>> A = tensorhead('A', [Lorentz]*2, [[1]*2]) + >>> A(a, -b) + A(a, -b) + + """ + sym = tensorsymmetry(*sym) + S = TensorType(typ, sym) + return S(name, comm) + + +class TensorHead(Basic): + """ + Tensor head of the tensor + + Parameters + ========== + + name : name of the tensor + + typ : list of TensorIndexType + + comm : commutation group number + + Attributes + ========== + + ``name`` + ``index_types`` + ``rank`` + ``types`` : equal to ``typ.types`` + ``symmetry`` : equal to ``typ.symmetry`` + ``comm`` : commutation group + + Notes + ===== + + A ``TensorHead`` belongs to a commutation group, defined by a + symbol on number ``comm`` (see ``_TensorManager.set_comm``); + tensors in a commutation group have the same commutation properties; + by default ``comm`` is ``0``, the group of the commuting tensors. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensorsymmetry, TensorType + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> sym2 = tensorsymmetry([1]*2) + >>> S2 = TensorType([Lorentz]*2, sym2) + >>> A = S2('A') + """ + is_commutative = False + + def __new__(cls, name, typ, comm, **kw_args): + assert isinstance(name, basestring) + + obj = Basic.__new__(cls, name, typ, **kw_args) + obj._name = obj.args[0] + obj._rank = len(obj.index_types) + obj._types = typ.types + obj._symmetry = typ.symmetry + obj._comm = TensorManager.comm_symbols2i(comm) + return obj + + @property + def name(self): + return self._name + + @property + def rank(self): + return self._rank + + @property + def types(self): + return self._types[:] + + @property + def symmetry(self): + return self._symmetry + + @property + def typ(self): + return self.args[1] + + @property + def comm(self): + return self._comm + + @property + def index_types(self): + return self.args[1].index_types[:] + + def __lt__(self, other): + return (self.name, self.index_types) < (other.name, other.index_types) + + def _hashable_content(self): + r = (self._name, tuple(self._types), self._symmetry, self._comm) + return r + + def commutes_with(self, other): + """ + Returns 0 (1) if self and other (anti)commute. + + Returns None if self and other do not (anti)commute. + """ + r = TensorManager.get_comm(self._comm, other._comm) + return r + + + def _pretty(self): + return '%s(%s)' %(self.name, ','.join([str(x) for x in self.index_types])) + + def __call__(self, *indices): + """ + Returns a tensor with indices. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a, b = tensor_indices('a,b', Lorentz) + >>> A = tensorhead('A', [Lorentz]*2, [[1]*2]) + >>> t = A(a, -b) + """ + if not [indices[i]._tensortype for i in range(len(indices))] == self.index_types: + raise ValueError('wrong index type') + components = [self] + free, dum = TensMul.from_indices(*indices) + free.sort(key=lambda x: x[0].name) + dum.sort() + return TensMul(S.One, components, free, dum) + +class HoldTensorHead(TensorHead): + def __new__(cls, t, **kw_args): + """ + + ``t`` tensor expression + + Keyword arguments: + ``symmetry`` symmetry of the tensor expression; + if not assigned no symmetry is assumed + + ``comm`` commutation group; + if not assigned it is assumed to be not commuting + + Examples + ==== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead, HoldTensorHead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a,b,c,d = tensor_indices('a,b,c,d', Lorentz) + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> h1 = HoldTensorHead(p(a)*q(b) - p(b)*q(a), symm=[[2]], comm=0) + >>> h2 = HoldTensorHead(p(a)*q(b) - p(b)*q(a), symm=[[1]*2], comm=0) + >>> (h1(a,b)*h2(-a,-b)).canon_bp() + 0 + >>> t = h1(a,b)*h1(-a,-b) + >>> t.expand(True) - 2*p(a)*p(-a)*q(b)*q(-b) + 2*p(a)*q(-a)*q(b)*p(-b) + 0 + """ + types = [x.tensortype for x in t.free_args] + rank = len(t.free_args) + symmetry = kw_args.get('symm', [[1]]*rank) + if isinstance(symmetry, list): + symmetry = tensorsymmetry(*symmetry) + comm = kw_args.get('comm') + if not comm: + comm = 2 + else: + comm = TensorManager.comm_symbols2i(comm) + name = str(t) + typ = TensorType(types, symmetry) + obj = Basic.__new__(cls, name, typ, t) + obj._name = name + obj.free_args = t.free_args + obj._comm = comm + obj._symmetry = symmetry + obj._types = types + obj._rank = rank + return obj + + name = property(lambda self: self._name) + t = property(lambda self: self.args[2]) + + + def call(self, *indices): + t = self.t.substitute_indices(*zip(self.free_args, indices)) + return t + + +class TensExpr(Basic): + """ + Abstract base class for tensor expressions + + Notes + ===== + + A tensor expression is an expression formed by tensors; + currently the sums of tensors are distributed. + + A ``TensExpr`` can be a ``TensAdd`` or a ``TensMul``. + + ``TensAdd`` objects are put in canonical form using the Butler-Portugal + algorithm for canonicalization under monoterm symmetries. + + ``TensMul`` objects are formed by products of component tensors, + and include a coefficient, which is a SymPy expression. + + + In the internal representation contracted indices are represented + by ``(ipos1, ipos2, icomp1, icomp2)``, where ``icomp1`` is the position + of the component tensor with contravariant index, ``ipos1`` is the + slot which the index occupies in that component tensor. + + Contracted indices are therefore nameless in the internal representation. + """ + + _op_priority = 11.0 + is_commutative = False + + def __neg__(self): + return self*S.NegativeOne + + def __abs__(self): + raise NotImplementedError + + def __add__(self, other): + raise NotImplementedError + + def __radd__(self, other): + raise NotImplementedError + + def __sub__(self, other): + raise NotImplementedError + + def __rsub__(self, other): + raise NotImplementedError + + def __mul__(self, other): + raise NotImplementedError + + def __rmul__(self, other): + return self*other + + def __pow__(self, other): + raise NotImplementedError + + def __rpow__(self, other): + raise NotImplementedError + + def __div__(self, other): + raise NotImplementedError + + def __rdiv__(self, other): + raise NotImplementedError() + + __truediv__ = __div__ + __rtruediv__ = __rdiv__ + + +def _tensAdd_collect_terms(args): + """ + collect TensMul terms differing at most by their coefficient + """ + a = [] + pprev = None + prev = args[0] + prev_coeff = prev._coeff + changed = False + + for x in args[1:]: + # if x and prev have the same tensor, update the coeff of prev + if x._components == prev._components \ + and x._free == prev._free and x._dum == prev._dum: + prev_coeff = prev_coeff + x._coeff + changed = True + op = 0 + else: + # x and prev are different; if not changed, prev has not + # been updated; store it + if not changed: + a.append(prev) + else: + # get a tensor from prev with coeff=prev_coeff and store it + if prev_coeff: + t = TensMul(prev_coeff, prev._components, + prev._free, prev._dum) + a.append(t) + # move x to prev + op = 1 + pprev, prev = prev, x + pprev_coeff, prev_coeff = prev_coeff, x._coeff + changed = False + # if the case op=0 prev was not stored; store it now + # in the case op=1 x was not stored; store it now (as prev) + if op == 0 and prev_coeff: + prev = TensMul(prev_coeff, prev._components, prev._free, prev._dum) + a.append(prev) + elif op == 1: + a.append(prev) + return a + +def _tensAdd_flatten(args): + """ + flatten TensAdd, coerce terms which are not tensors to tensors + """ + if not all(isinstance(x, TensExpr) for x in args): + args1 = [] + for x in args: + if isinstance(x, TensExpr): + if isinstance(x, TensAdd): + args1.extend(list(x.args)) + else: + args1.append(x) + args1 = [x for x in args1 if isinstance(x, TensExpr) and x._coeff] + args2 = [x for x in args if not isinstance(x, TensExpr)] + t1 = TensMul(Add(*args2), [], [], []) + args = [t1] + args1 + a = [] + for x in args: + if isinstance(x, TensAdd): + a.extend(list(x.args)) + else: + a.append(x) + args = [x for x in a if x._coeff] + return args + +def _tensAdd_check(args): + """ + check that all addends have the same free indices + """ + indices0 = sorted([x[0] for x in args[0]._free], key=lambda x: x.name) + list_indices = [sorted([y[0] for y in x._free], key=lambda x: x.name) for x in args[1:]] + if not all(x == indices0 for x in list_indices): + raise ValueError('all tensors must have the same indices') + + +class TensAdd(TensExpr): + """ + Sum of tensors + + Parameters + ========== + + free_args : list of the free indices + + Attributes + ========== + + ``args`` : tuple of addends + ``rank`` : rank of the tensor + ``free_args`` : list of the free indices in sorted order + + Notes + ===== + + Sum of more than one tensor are put automatically in canonical form. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensorhead, tensor_indices + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a, b = tensor_indices('a,b', Lorentz) + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> t = p(a) + q(a); t + p(a) + q(a) + >>> t(b) + p(b) + q(b) + """ + + def __new__(cls, *args, **kw_args): + args = [sympify(x) for x in args if x] + args = _tensAdd_flatten(args) + if not args: + return S.Zero + + _tensAdd_check(args) + obj = Basic.__new__(cls, **kw_args) + if len(args) == 1 and isinstance(args[0], TensMul): + obj._args = tuple(args) + return obj + args = [x.canon_bp() for x in args if x] + args = [x for x in args if x] + if not args: + return S.Zero + + # collect canonicalized terms + args.sort(key=lambda x: (x._components, x._free, x._dum)) + a = _tensAdd_collect_terms(args) + if not a: + return S.Zero + # it there is only a component tensor return it + if len(a) == 1: + return a[0] + obj._args = tuple(a) + return obj + + @property + def free_args(self): + return self.args[0].free_args + + @property + def rank(self): + return self.args[0].rank + + def __call__(self, *indices): + """Returns tensor with ordered free indices replaced by ``indices`` + + Parameters + ========== + + indices + + Examples + ======== + + >>> from sympy import Symbol + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> i0,i1,i2,i3,i4 = tensor_indices('i0:5', Lorentz) + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> g = Lorentz.metric + >>> t = p(i0)*p(i1) + g(i0,i1)*q(i2)*q(-i2) + >>> t(i0,i2) + metric(i0, i2)*q(L_0)*q(-L_0) + p(i0)*p(i2) + >>> t(i0,i1) - t(i1,i0) + 0 + """ + free_args = self.free_args + indices = list(indices) + if [x._tensortype for x in indices] != [x._tensortype for x in free_args]: + raise ValueError('incompatible types') + if indices == free_args: + return self + index_tuples = zip(free_args, indices) + a = [x.fun_eval(*index_tuples) for x in self.args] + res = TensAdd(*a) + + return res + + def canon_bp(self): + """ + canonicalize using the Butler-Portugal algorithm for canonicalization + under monoterm symmetries. + """ + args = [x.canon_bp() for x in self.args] + res = TensAdd(*args) + return res + + def __eq__(self, other): + other = sympify(other) + if isinstance(other, TensMul) and other._coeff == 0: + return self == 0 + t = self - other + if not isinstance(t, TensExpr): + return t == 0 + else: + if isinstance(t, TensMul): + return t._coeff == 0 + else: + return all(x._coeff == 0 for x in t.args) + + def __add__(self, other): + return TensAdd(self, other) + + def __radd__(self, other): + return TensAdd(other, self) + + def __sub__(self, other): + return TensAdd(self, -other) + + def __rsub__(self, other): + return TensAdd(other, -self) + + def __mul__(self, other): + return TensAdd(*[x*other for x in self.args]) + + def __div__(self, other): + other = sympify(other) + if isinstance(other, TensExpr): + raise ValueError('cannot divide by a tensor') + return TensAdd(*[x/other for x in self.args]) + + def __rdiv__(self, other): + raise ValueError('cannot divide by a tensor') + + __truediv__ = __div__ + __truerdiv__ = __rdiv__ + + def _hashable_content(self): + return tuple(self.args) + + def __hash__(self): + return super(TensAdd, self).__hash__() + + def __ne__(self, other): + return not (self == other) + + def contract_delta(self, delta): + args = [x.contract_delta(delta) for x in self.args] + t = TensAdd(*args) + return canon_bp(t) + + def contract_metric(self, g, contract_all=False): + """ + Raise or lower indices with the metric ``g`` + + Parameters + ========== + + g : metric + + contract_all : if True, eliminate all ``g`` which are contracted + + Notes + ===== + + see the ``TensorIndexType`` docstring for the contraction conventions + """ + + args = [x.contract_metric(g, contract_all) for x in self.args] + t = TensAdd(*args) + return canon_bp(t) + + def substitute_tensor(self, t1, t2, subst_all=False): + args = self.args + args1 = [] + for x in args: + y = x.substitute_tensor(t1, t2, subst_all) + args1.append(y) + return TensAdd(*args1) + + def fun_eval(self, *index_tuples): + """ + Return a tensor with free indices substituted according to ``index_tuples`` + + Parameters + ========== + + index_types : list of tuples ``(old_index, new_index)`` + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz) + >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2]) + >>> t = A(i, k)*B(-k, -j) + A(i, -j) + >>> t.fun_eval((i, k),(-j, l)) + A(k, L_0)*B(l, -L_0) + A(k, l) + """ + args = self.args + args1 = [] + for x in args: + y = x.fun_eval(*index_tuples) + args1.append(y) + return TensAdd(*args1) + + def expand_coeff(self): + args = [_tensor_expand_coeff(x) for x in self.args] + return TensAdd(*args) + + def expand(self, expand_all=False): + """ + expand HoldTensorHead tensors + """ + a = [x.expand(expand_all) for x in self.args] + return TensAdd(*a) + + def substitute_indices(self, *index_tuples): + """ + Return a tensor with free indices substituted according to ``index_tuples`` + + Parameters + ========== + + index_types : list of tuples ``(old_index, new_index)`` + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz) + >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2]) + >>> t = A(i, k)*B(-k, -j); t + A(i, L_0)*B(-L_0, -j) + >>> t.substitute_indices((i,j), (j, k)) + A(j, L_0)*B(-L_0, -k) + """ + args = self.args + args1 = [] + for x in args: + y = x.substitute_indices(*index_tuples) + args1.append(y) + return TensAdd(*args1) + + def _pretty(self): + a = [] + args = self.args + for x in args: + a.append(str(x)) + a.sort() + s = ' + '.join(a) + s = s.replace('+ -', '- ') + return s + +class TensMul(TensExpr): + """ + Product of tensors + + Parameters + ========== + + coeff : SymPy coefficient of the tensor + args + + Attributes + ========== + + ``_components`` : list of ``TensorHead`` of the component tensors + ``types`` : list of nonrepeated ``TensorIndexType`` + ``free`` : list of ``(ind, ipos, icomp)``, see Notes + ``dum`` : list of ``(ipos1, ipos2, icomp1, icomp2)``, see Notes + ``ext_rank`` : rank of the tensor counting the dummy indices + ``rank`` : rank of the tensor + ``coeff`` : SymPy coefficient of the tensor + ``free_args`` : list of the free indices in sorted order + ``is_canon_bp`` : ``True`` if the tensor in in canonical form + + Notes + ===== + + ``args[0]`` list of ``TensorHead`` of the component tensors. + + ``args[1]`` list of ``(ind, ipos, icomp)`` + where ``ind`` is a free index, ``ipos`` is the slot position + of ``ind`` in the ``icomp``-th component tensor. + + ``args[2]`` list of tuples representing dummy indices. + ``(ipos1, ipos2, icomp1, icomp2)`` indicates that the contravariant + dummy index is the ``ipos1``-th slot position in the ``icomp1``-th + component tensor; the corresponding covariant index is + in the ``ipos2`` slot position in the ``icomp2``-th component tensor. + """ + + def __new__(cls, coeff, *args, **kw_args): + obj = Basic.__new__(cls) + obj._components = args[0] + obj._types = [] + for t in obj._components: + obj._types.extend(t._types) + obj._free = args[1] + obj._dum = args[2] + obj._ext_rank = len(obj._free) + 2*len(obj._dum) + obj._coeff = coeff + obj._is_canon_bp = kw_args.get('is_canon_bp', False) + + return obj + + @property + def free_args(self): + return sorted([x[0] for x in self._free]) + + @property + def components(self): + return self._components[:] + + @property + def free(self): + return self._free[:] + + @property + def coeff(self): + return self._coeff + + @property + def dum(self): + return self._dum[:] + + @property + def rank(self): + return len(self._free) + + @property + def types(self): + return self._types[:] + + def __eq__(self, other): + if other == 0: + return self._coeff == 0 + other = sympify(other) + if not isinstance(other, TensExpr): + assert not self._components + return self._coeff == other + res = self - other + return res == 0 + + def _hashable_content(self): + t = self.canon_bp() + r = (t._coeff, tuple(t._components), \ + tuple(sorted(t._free)), tuple(sorted(t._dum))) + return r + + def __hash__(self): + return super(TensMul, self).__hash__() + + def __ne__(self, other): + return not self == other + + def expand(self, expand_all=False): + """ + expand HoldTensorHead tensors + """ + if len(self._components) == 1: + if isinstance(self._components[0], HoldTensorHead): + return self._components[0].call(*self.free_args) + else: + return self + a = self.split() + if expand_all: + a1 = [x.expand(expand_all) for x in a] + else: + a1 = [] + expanded = False + for i, t in enumerate(a): + t1 = t.expand() + if t1 is not t: + expanded = True + a1.append(t1) + if expanded is True: + break + a1 += a[i + 1:] + return tensor_mul(*a1) + + @staticmethod + def from_indices(*indices): + """ + Convert ``indices`` into ``free``, ``dum`` for single component tensor + + ``free`` list of tuples ``(index, pos, 0)``, + where ``pos`` is the position of index in + the list of indices formed by the component tensors + + ``dum`` list of tuples ``(pos_contr, pos_cov, 0, 0)`` + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensMul + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz) + >>> TensMul.from_indices(m0, m1, -m1, m3) + ([(m0, 0, 0), (m3, 3, 0)], [(1, 2, 0, 0)]) + """ + n = len(indices) + if n == 1: + return [(indices[0], 0, 0)], [] + + # find the positions of the free indices and of the dummy indices + free = [True]*len(indices) + index_dict = {} + dum = [] + for i, index in enumerate(indices): + name = index._name + typ = index._tensortype + contr = index._is_up + if (name, typ) in index_dict: + # found a pair of dummy indices + is_contr, pos = index_dict[(name, typ)] + # check consistency and update free + if is_contr: + if contr: + raise ValueError('two equal contravariant indices in slots %d and %d' %(pos, i)) + else: + free[pos] = False + free[i] = False + else: + if contr: + free[pos] = False + free[i] = False + else: + raise ValueError('two equal covariant indices in slots %d and %d' %(pos, i)) + if contr: + dum.append((i, pos, 0, 0)) + else: + dum.append((pos, i, 0, 0)) + else: + index_dict[(name, typ)] = index._is_up, i + + free = [(index, i, 0) for i, index in enumerate(indices) if free[i]] + free.sort() + return free, dum + + def get_indices(self): + """ + Returns the list of indices of the tensor + + The indices are listed in the order in which they appear in the + component tensors. + The dummy indices are given a name which does not collide with + the names of the free indices. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz) + >>> g = Lorentz.metric + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> t = p(m1)*g(m0,m2) + >>> t.get_indices() + [m1, m0, m2] + """ + indices = [None]*self._ext_rank + start = 0 + pos = 0 + vpos = [] + components = self._components + for t in components: + vpos.append(pos) + pos += t._rank + cdt = defaultdict(int) + # if the free indices have names with dummy_fmt, start with an + # index higher than those for the dummy indices + # to avoid name collisions + for indx, ipos, cpos in self._free: + if indx._name.split('_')[0] == indx._tensortype._dummy_fmt[:-3]: + cdt[indx._tensortype] = max(cdt[indx._tensortype], int(indx._name.split('_')[1]) + 1) + start = vpos[cpos] + indices[start + ipos] = indx + for ipos1, ipos2, cpos1, cpos2 in self._dum: + start1 = vpos[cpos1] + start2 = vpos[cpos2] + typ1 = components[cpos1].index_types[ipos1] + assert typ1 == components[cpos2].index_types[ipos2] + fmt = typ1.dummy_fmt + nd = cdt[typ1] + indices[start1 + ipos1] = TensorIndex(fmt % nd, typ1) + indices[start2 + ipos2] = TensorIndex(fmt % nd, typ1, False) + cdt[typ1] += 1 + return indices + + def split(self): + """ + Returns a list of tensors, whose product is ``self`` + + Dummy indices contracted among different tensor components + become free indices with the same name as the one used to + represent the dummy indices. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz) + >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2]) + >>> t = A(a,b)*B(-b,c) + >>> t + A(a, L_0)*B(-L_0, c) + >>> t.split() + [A(a, L_0), B(-L_0, c)] + """ + indices = self.get_indices() + pos = 0 + components = self._components + if not components: + return [TensMul(self._coeff, [], [], [])] + res = [] + for t in components: + t1 = t(*indices[pos:pos + t._rank]) + pos += t._rank + res.append(t1) + res[0] = TensMul(self._coeff, res[0]._components, res[0]._free, res[0]._dum, is_canon_bp=res[0]._is_canon_bp) + return res + + def canon_args(self): + """ + Returns ``(g, dummies, msym, v)``, the entries of ``canonicalize`` + + see ``canonicalize`` in ``tensor_can.py`` + """ + # to be called after sorted_components + from sympy.combinatorics.permutations import _af_new + types = list(set(self._types)) + types.sort(key = lambda x: x._name) + n = self._ext_rank + g = [None]*n + [n, n+1] + pos = 0 + vpos = [] + components = self._components + for t in components: + vpos.append(pos) + pos += t._rank + # ordered indices: first the free indices, ordered by types + # then the dummy indices, ordered by types and contravariant before + # covariant + # g[position in tensor] = position in ordered indices + for i, (indx, ipos, cpos) in enumerate(self._free): + pos = vpos[cpos] + ipos + g[pos] = i + pos = len(self._free) + j = len(self._free) + dummies = [] + prev = None + a = [] + msym = [] + for ipos1, ipos2, cpos1, cpos2 in self._dum: + pos1 = vpos[cpos1] + ipos1 + pos2 = vpos[cpos2] + ipos2 + g[pos1] = j + g[pos2] = j + 1 + j += 2 + typ = components[cpos1].index_types[ipos1] + if typ != prev: + if a: + dummies.append(a) + a = [pos, pos + 1] + prev = typ + msym.append(typ.metric_antisym) + else: + a.extend([pos, pos + 1]) + pos += 2 + if a: + dummies.append(a) + numtyp = [] + prev = None + for t in components: + if t == prev: + numtyp[-1][1] += 1 + else: + prev = t + numtyp.append([prev, 1]) + v = [] + for h, n in numtyp: + if h._comm == 0 or h._comm == 1: + comm = h._comm + else: + comm = TensorManager.get_comm(h._comm, h._comm) + v.append((h._symmetry.base, h._symmetry.generators, n, comm)) + return _af_new(g), dummies, msym, v + + def __add__(self, other): + return TensAdd(self, other) + + def __radd__(self, other): + return TensAdd(other, self) + + def __sub__(self, other): + return TensAdd(self, -other) + + def __rsub__(self, other): + return TensAdd(other, -self) + + def __mul__(self, other): + """ + Multiply two tensors using Einstein summation convention. + + If the two tensors have an index in common, one contravariant + and the other covariant, in their product the indices are summed + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz) + >>> g = Lorentz.metric + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> t1 = p(m0) + >>> t2 = q(-m0) + >>> t1*t2 + p(L_0)*q(-L_0) + """ + other = sympify(other) + if not isinstance(other, TensExpr): + coeff = self._coeff*other + return TensMul(coeff, self._components, self._free, self._dum, is_canon_bp=self._is_canon_bp) + if isinstance(other, TensAdd): + return TensAdd(*[self*x for x in other.args]) + + components = self._components + other._components + # find out which free indices of self and other are contracted + free_dict1 = dict([(i.name, (pos, cpos, i)) for i, pos, cpos in self._free]) + free_dict2 = dict([(i.name, (pos, cpos, i)) for i, pos, cpos in other._free]) + + free_names = set(free_dict1.keys()) & set(free_dict2.keys()) + # find the new `free` and `dum` + nc1 = len(self._components) + dum2 = [(i1, i2, c1 + nc1, c2 + nc1) for i1, i2, c1, c2 in other._dum] + free1 = [(ind, i, c) for ind, i, c in self._free if ind.name not in free_names] + free2 = [(ind, i, c + nc1) for ind, i, c in other._free if ind.name not in free_names] + free = free1 + free2 + dum = self._dum + dum2 + for name in free_names: + ipos1, cpos1, ind1 = free_dict1[name] + ipos2, cpos2, ind2 = free_dict2[name] + cpos2 += nc1 + if ind1._is_up == ind2._is_up: + raise ValueError('wrong index contruction %s' % ind1) + if ind1._is_up: + new_dummy = (ipos1, ipos2, cpos1, cpos2) + else: + new_dummy = (ipos2, ipos1, cpos2, cpos1) + dum.append(new_dummy) + coeff = self._coeff*other._coeff + return TensMul(coeff, components, free, dum) + + def __rmul__(self, other): + other = sympify(other) + coeff = other*self._coeff + return TensMul(coeff, self._components, self._free, self._dum) + + def __div__(self, other): + other = sympify(other) + if isinstance(other, TensExpr): + raise ValueError('cannot divide by a tensor') + coeff = self._coeff/other + return TensMul(coeff, self._components, self._free, self._dum, is_canon_bp=self._is_canon_bp) + + def __rdiv__(self, other): + raise ValueError('cannot divide by a tensor') + + __truediv__ = __div__ + __truerdiv__ = __rdiv__ + + def sorted_components(self): + """ + Returns a tensor with sorted components + + The sorting is done taking into account the commutation group + of the component tensors. + """ + from sympy.combinatorics.permutations import _af_invert + cv = zip(self._components, range(len(self._components))) + sign = 1 + n = len(cv) - 1 + for i in range(n): + for j in range(n, i, -1): + c = cv[j-1][0].commutes_with(cv[j][0]) + if c not in [0, 1]: + continue + if (cv[j-1][0]._types, cv[j-1][0]._name) > \ + (cv[j][0]._types, cv[j][0]._name): + cv[j-1], cv[j] = cv[j], cv[j-1] + if c: + sign = -sign + # perm_inv[new_pos] = old_pos + components = [x[0] for x in cv] + perm_inv = [x[1] for x in cv] + perm = _af_invert(perm_inv) + free = [(ind, i, perm[c]) for ind, i, c in self._free] + free.sort() + dum = [(i1, i2, perm[c1], perm[c2]) for i1, i2, c1, c2 in self._dum] + dum.sort(key = lambda x: components[x[2]].index_types[x[0]]) + coeff = -self._coeff if sign == -1 else self._coeff + t = TensMul(coeff, components, free, dum) + return t + + def perm2tensor(self, g, canon_bp=False): + """ + Returns the tensor corresponding to the permutation ``g`` + + ``g`` permutation corrisponding to the tensor in the representation + used in canonicalization + + ``canon_bp`` if True, then ``g`` is the permutation + corresponding to the canonical form of the tensor + """ + from bisect import bisect_right + vpos = [] + components = self._components + pos = 0 + for t in components: + vpos.append(pos) + pos += t._rank + sorted_free = [x[0] for x in self._free] + sorted_free.sort() + nfree = len(sorted_free) + rank = self._ext_rank + indices = [None]*rank + dum = [[None]*4 for i in range((rank - nfree)//2)] + free = [] + icomp = -1 + for i in range(rank): + if i in vpos: + icomp += vpos.count(i) + pos0 = i + ipos = i - pos0 + gi = g[i] + if gi < nfree: + ind = sorted_free[gi] + free.append((ind, ipos, icomp)) + else: + j = gi - nfree + idum, cov = divmod(j, 2) + if cov: + dum[idum][1] = ipos + dum[idum][3] = icomp + else: + dum[idum][0] = ipos + dum[idum][2] = icomp + dum = [tuple(x) for x in dum] + coeff = self._coeff + if g[-1] != len(g) - 1: + coeff = -coeff + res = TensMul(coeff, components, free, dum, is_canon_bp=canon_bp) + return res + + def canon_bp(self): + """ + canonicalize using the Butler-Portugal algorithm for canonicalization + under monoterm symmetries. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz) + >>> A = tensorhead('A', [Lorentz]*2, [[2]]) + >>> t = A(m0,-m1)*A(m1,-m0) + >>> t.canon_bp() + -A(L_0, L_1)*A(-L_0, -L_1) + >>> t = A(m0,-m1)*A(m1,-m2)*A(m2,-m0) + >>> t.canon_bp() + 0 + """ + from sympy.combinatorics.tensor_can import canonicalize + if self._is_canon_bp: + return self + if not self._components: + return self + t = self.sorted_components() + g, dummies, msym, v = t.canon_args() + can = canonicalize(g, dummies, msym, *v) + if can == 0: + return S.Zero + return t.perm2tensor(can, True) + + def _contract(self, g, antisym, contract_all=False): + """ + helper method for ``contract_metric`` and ``contract_delta`` + + ``g`` metric to be contracted + + ``antisym``: + False symmetric metric + True antisymmetric metric + None delta + """ + if not self._components: + return self + free_indices = [x[0] for x in self._free] + a = self.split() + typ = g.index_types[0] + for i, tg in enumerate(a): + if tg._components[0] == g: + tg_free = [x[0] for x in tg._free] + if len(tg_free) == 0: + t = _contract_g_with_itself(a, i, tg, tg_free, g, antisym) + if contract_all == True and g in t._components: + return t._contract(g, antisym, True) + return t + if all(indx in free_indices for indx in tg_free): + continue + else: + break + else: + # all metric tensors have only free indices, + # there is no contraction + return self + + # tg has one or two indices contracted with other tensors + # i position of tg in a + coeff = S.One + tg_free = tg._free + if antisym: + # order by slot position + tg_free = sorted(tg_free, key=lambda x: x[1]) + + if tg_free[0][0] in free_indices or tg_free[1][0] in free_indices: + # tg has one free index + res = _contract_g_with_free_index(a, free_indices, i, tg, tg_free, g, antisym) + else: + # tg has two indices contracted with other tensors + res = _contract_g_without_free_index(a, free_indices, i, tg, tg_free, g, typ, antisym) + if contract_all == True and g in res._components: + return res._contract(g, antisym, True) + return res + + def contract_delta(self, delta): + typ = delta._types[0] + t = self._contract(delta, None, True) + return t + + def contract_metric(self, g, contract_all=False): + """ + Raise or lower indices with the metric ``g`` + + ``g`` metric + + ``contract_all`` if True, eliminate all ``g`` which are contracted + + Notes + ===== + + see the ``TensorIndexType`` docstring for the contraction conventions + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz) + >>> g = Lorentz.metric + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> t = p(m0)*q(m1)*g(-m0, -m1) + >>> t.canon_bp() + metric(L_0, L_1)*p(-L_0)*q(-L_1) + >>> t.contract_metric(g).canon_bp() + p(L_0)*q(-L_0) + """ + return self._contract(g, g.index_types[0].metric_antisym, contract_all) + + def substitute_indices(self, *index_tuples): + """ + Return a tensor with free indices substituted according to ``index_tuples`` + + ``index_types`` list of tuples ``(old_index, new_index)`` + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz) + >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2]) + >>> t = A(i, k)*B(-k, -j); t + A(i, L_0)*B(-L_0, -j) + >>> t.substitute_indices((i,j), (j, k)) + A(j, L_0)*B(-L_0, -k) + """ + free = self._free + free1 = [] + for j, ipos, cpos in free: + for i, v in index_tuples: + if i._name == j._name and i._tensortype == j._tensortype: + if i._is_up == j._is_up: + free1.append((v, ipos, cpos)) + else: + free1.append((-v, ipos, cpos)) + break + else: + free1.append((j, ipos, cpos)) + + return TensMul(self._coeff, self._components, free1, self._dum) + + def substitute_tensor(self, t1, t2, subst_all=False): + """ + Return tensor obtained substituting ``t1`` with ``t2`` in ``self``. + """ + # FIXME fix the check + #assert sorted(t1._free) == sorted(t2._free) + components = self._components + assert len(t1._components) == 1 and t1._coeff == 1 + component1 = t1._components[0] + for i in range(len(components)): + if component1 == components[i]: + break + else: + return self + a = self.split() + if a[i]._dum != t1._dum: + return self + + ct = self._coeff if i == 0 else S.One + # sort free indices according to the position + free0 = sorted(a[i]._free, key=lambda x: x[1]) + free1 = sorted(t1._free, key=lambda x: x[1]) + index_tuples = [] + for j in range(len(free0)): + index_tuples.append((free1[j][0], free0[j][0])) + t2s = t2.substitute_indices(*index_tuples) + a1 = a[:i] + [t2s] + a[i + 1:] + t = tensor_mul(*a1)*ct + if subst_all: + return t.substitute_tensor(t1, t2, subst_all) + return t + + def substitute_tensors(self, i, j, t, a=None): + """ + substitute tensors in positions ``i, j`` with tensor ``t`` + """ + if not a: + a = self.split() + assert i < j + ct = self._coeff if i == 0 else S.One + #a = self.split() + a1 = a[:i] + a[i + 1:j] + a[j + 1:] + [t] + t = tensor_mul(*a1)*ct + return t + + def fun_eval(self, *index_tuples): + """ + Return a tensor with free indices substituted according to ``index_tuples`` + + ``index_types`` list of tuples ``(old_index, new_index)`` + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz) + >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2]) + >>> t = A(i, k)*B(-k, -j); t + A(i, L_0)*B(-L_0, -j) + >>> t.fun_eval((i, k),(-j, l)) + A(k, L_0)*B(-L_0, l) + """ + free = self._free + free1 = [] + for j, ipos, cpos in free: + # search j in index_tuples + for i, v in index_tuples: + if i == j: + free1.append((v, ipos, cpos)) + break + else: + free1.append((j, ipos, cpos)) + return TensMul(self._coeff, self._components, free1, self._dum) + + + def __call__(self, *indices): + """Returns tensor with ordered free indices replaced by ``indices`` + + Examples + ======== + + >>> from sympy import Symbol + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead + >>> D = Symbol('D') + >>> Lorentz = TensorIndexType('Lorentz', dim=D, dummy_fmt='L') + >>> i0,i1,i2,i3,i4 = tensor_indices('i0:5', Lorentz) + >>> g = Lorentz.metric + >>> p, q = tensorhead('p,q', [Lorentz], [[1]]) + >>> t = p(i0)*q(i1)*q(-i1) + >>> t(i1) + p(i1)*q(L_0)*q(-L_0) + """ + free_args = self.free_args + indices = list(indices) + if [x._tensortype for x in indices] != [x._tensortype for x in free_args]: + raise ValueError('incompatible types') + if indices == free_args: + return self + t = self.fun_eval(*zip(free_args, indices)) + return t + + def _as_coeff_mul(self): + a = self.split() + coeff = a[0]._coeff() + t0 = TensMul(S.One, a[0]._components, a[0]._free, a[0]._dum) + a[0] = t0 + t1 = tensor_mul(*a) + return coeff, t1 + + def expand_coeff(self): + coeff, t = self._as_coeff_mul() + return coeff.expand()*t + + def _pretty(self): + if self._components == []: + return str(self._coeff) + indices = self.get_indices() + sindices = [str(ind) for ind in indices] + pos = 0 + a = [] + for t in self._components: + if isinstance(t, HoldTensorHead): + rt = '(%s)' % (t.call(*indices[pos:pos + t.rank])) + a.append(rt) + else: + if t.rank > 0: + a.append('%s(%s)' % (t.name, ', '.join(sindices[pos:pos + t.rank]))) + else: + a.append('%s' % t.name) + pos += t._rank + res = '*'. join(a) + if self._coeff == S.One: + return res + elif self._coeff == -S.One: + return '-%s' % res + if self._coeff.is_Atom: + return '%s*%s' % (self._coeff, res) + else: + return '(%s)*%s' %(self._coeff, res) + + + +def canon_bp(p): + """ + Butler-Portugal canonicalization + """ + if isinstance(p, TensExpr): + return p.canon_bp() + return p + +def tensor_mul(*a): + """ + product of tensors + """ + if not a: + return TensMul(S.One, [], [], []) + t = a[0] + for tx in a[1:]: + t = t*tx + return t + +def _tensor_expand_coeff(p): + if not isinstance(p, TensExpr): + return p + return p.expand_coeff() + + + +def riemann_cyclic_replace(t_r): + """ + replace Riemann tensor with an equivalent expression + + ``R(m,n,p,q) -> 2/3*R(m,n,p,q) - 1/3*R(m,q,n,p) + 1/3*R(m,p,n,q)`` + + """ + free = sorted(t_r._free, key=lambda x: x[1]) + m, n, p, q = [x[0] for x in free] + t0 = S(2)/3*t_r + t1 = - S(1)/3*t_r.substitute_indices((m,m),(n,q),(p,n),(q,p)) + t2 = S(1)/3*t_r.substitute_indices((m,m),(n,p),(p,n),(q,q)) + t3 = t0 + t1 + t2 + return t3 + +def riemann_cyclic(t2): + """ + replace each Riemann tensor with an equivalent expression + satisfying the cyclic identity. + + This trick is discussed in the reference guide to Cadabra. + + Examples + ======== + + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead, riemann_cyclic + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz) + >>> R = tensorhead('R', [Lorentz]*4, [[2, 2]]) + >>> t = R(i,j,k,l)*(R(-i,-j,-k,-l) - 2*R(-i,-k,-j,-l)) + >>> riemann_cyclic(t) + 0 + """ + if isinstance(t2, TensMul): + args = [t2] + else: + args = t2.args + a1 = [x.split() for x in args] + a2 = [[riemann_cyclic_replace(tx) for tx in y] for y in a1] + a3 = [tensor_mul(*v) for v in a2] + t3 = TensAdd(*a3) + if not t3: + return t3 + else: + return canon_bp(t3) + + + +def tensorlist_contract_metric(a, tg): + """ + contract `tg` with a tensor in the list `a = t.split()` + Only for symmetric metric. + """ + ind1, ind2 = [x[0] for x in tg._free] + mind1 = -ind1 + mind2 = -ind2 + for i in range(len(a)): + t1 = a[i] + for j in range(len(t1._free)): + indx, ipos, _ = t1._free[j] + if indx == mind1 or indx == mind2: + ind3 = ind2 if indx == mind1 else ind1 + free1 = t1._free[:] + free1[j] = (ind3, ipos, 0) + t2 = TensMul(t1._coeff, t1._components, free1, t1._dum) + a[i] = t2 + return a + a.append(tg) + return a + +def _contract_g_with_itself(a, i, tg, tg_free, g, antisym): + """ + helper function for _contract + """ + typ = g.index_types[0] + a1 = a[:i] + a[i + 1:] + t11 = tensor_mul(*a1) + if typ._dim is None: + raise ValueError('dimension not assigned') + coeff = typ._dim*a[i]._coeff + if antisym and tg._dum[0][0] == 0: + # g(i, -i) = -D + coeff = -coeff + t = tensor_mul(*a1)*coeff + return t + +def _contract_g_with_free_index(a, free_indices, i, tg, tg_free, g, antisym): + """ + helper function for _contract + """ + if tg_free[0][0] in free_indices: + ind_free = tg_free[0][0] + ind, ipos1, _ = tg_free[1] + else: + ind_free = tg_free[1][0] + ind, ipos1, _ = tg_free[0] + + ind1 = -ind + # search ind1 in the other component tensors + for j, tx in enumerate(a): + if ind1 in [x[0] for x in tx._free]: + break + # replace ind1 with ind_free + free1 = [] + for indx, iposx, _ in tx._free: + if indx == ind1: + free1.append((ind_free, iposx, 0)) + else: + free1.append((indx, iposx, 0)) + coeff = tx._coeff + if antisym: + if ind._is_up and ind == tg_free[0][0] or \ + (not ind._is_up) and ind == tg_free[1][0]: + # g(i1, i0)*psi(-i1) = -psi(i0) + # g(-i0, -i1)*psi(i1) = -psi(-i0) + coeff = -coeff + t1 = TensMul(coeff, tx._components, free1, tx._dum) + a[j] = t1 + a = a[:i] + a[i + 1:] + coeff = tg._coeff + res = tensor_mul(*a) + return coeff*res + +def _contract_g_without_free_index(a, free_indices, i, tg, tg_free, g, typ, antisym): + """ + helper function for _contract + """ + coeff = S.One + ind1 = tg_free[0][0] + ind2 = tg_free[1][0] + ind1m = -ind1 + ind2m = -ind2 + for k, ty in enumerate(a): + if ind2m in [x[0] for x in ty._free]: + break + # ty has the index ind2m + ty_free = ty._free[:] + if ty._components == [g]: + ty_indices = [x[0] for x in ty_free] + if all(x in [ind1m, ind2m] for x in ty_indices): + # the two `g` are completely contracted + # i < k always + a = a[:i] + a[i+1:k] + a[k+1:] + coeff = coeff*typ._dim*tg._coeff*ty._coeff + if antisym: + ty_free = sorted(ty_free, key=lambda x: x[1]) + if ind1._is_up == ind2._is_up: + # g(i,j)*g(-i,-j) = g(-i,-j)*g(i,j) = dim + # g(i,j)*g(-j,-i) = g(-i,-j)*g(j,i) = -dim + if ind1m == ty_free[1][0]: + coeff = -coeff + else: + # g(-i,j)*g(i,-j) = g(i,-j)^g(-i,j) = -dim + # g(-i,j)*g(-j,i) = g(i,-j)*g(j,i) = dim + if ind1m == ty_free[0][0]: + coeff = -coeff + + if a: + res = tensor_mul(*a) + res = coeff*res + else: + res = TensMul(coeff, [],[],[], is_canon_bp=True) + return res + + free2 = [] + ty_freeindices = [x[0] for x in ty_free] + if ind1m in ty_freeindices: + # tg has both indices contracted with ty + free2 = [(indx, iposx, cposx) for indx, iposx, cposx in ty._free if indx != ind1m and indx != ind2m] + dum2 = ty._dum[:] + for indx, iposx, _ in ty_free: + if indx == ind1m: + iposx1 = iposx + if indx == ind2m: + iposx2 = iposx + if antisym: + if ind1._is_up == ind2._is_up: + if iposx1 < iposx2: + coeff = -coeff + dum2.append((iposx1, iposx2, 0, 0)) + else: + dum2.append((iposx2, iposx1, 0, 0)) + else: + if iposx1 > iposx2: + coeff = -coeff + dum2.append((iposx2, iposx1, 0, 0)) + else: + dum2.append((iposx1, iposx2, 0, 0)) + else: + dum2.append((iposx1, iposx2, 0, 0)) + else: + # replace ind2m with ind1 in the free indices of ty + + free2 = [] + if not antisym: + for indx, iposx, _ in ty_free: + if indx == ind2m: + free2.append((ind1, iposx, 0)) + else: + free2.append((indx, iposx, 0)) + else: + for indx, iposx, _ in ty_free: + if indx == ind2m: + free2.append((ind1, iposx, 0)) + if indx._is_up: + coeff = -coeff + else: + free2.append((indx, iposx, 0)) + dum2 = ty._dum + + t2 = TensMul(ty._coeff, ty._components, free2, dum2) + a[k] = t2 + a = a[:i] + a[i + 1:] + coeff = coeff*tg._coeff + res = tensor_mul(*a) + return coeff*res diff --git a/sympy/tensor/tests/test_dgamma_matr.py b/sympy/tensor/tests/test_dgamma_matr.py new file mode 100644 index 000000000000..689f8b35a0a7 --- /dev/null +++ b/sympy/tensor/tests/test_dgamma_matr.py @@ -0,0 +1,278 @@ +from sympy import Symbol, S, I +from sympy.tensor.tensor import (TensorIndexType, tensor_indices, \ +TensorSymmetry, get_symmetric_group_sgs, TensorType, TensMul) +from sympy.tensor.dgamma_matr import GammaMatrices +from sympy.utilities.pytest import skip, XFAIL + + +D = Symbol('D') +Lorentz = TensorIndexType('Lorentz', dim=D, eps_dim=4, dummy_fmt='L') +m0, m1, m2, m3, m4, m5 = tensor_indices('m0,m1,m2,m3,m4,m5', Lorentz) +n0, n1, n2, n3, n4, n5 = tensor_indices('n0,n1,n2,n3,n4,n5', Lorentz) +sym1 = TensorSymmetry(get_symmetric_group_sgs(1)) +S1 = TensorType([Lorentz], sym1) + +GM = GammaMatrices( Lorentz) +match1_gamma = GM.match1_gamma +G = GM.G +g = GM.g +epsilon = GM.epsilon +rule1_gamma = GM.rule1_gamma +do_rule1_gamma = GM.do_rule1_gamma +match2_gamma = GM.match2_gamma +rule2_gamma = GM.rule2_gamma +do_rule2_gamma = GM.do_rule2_gamma +gamma_trace = GM.gamma_trace +gamma_trace1 = GM.gamma_trace1 +gctr = GM.gctr + + +def test_rule1_gamma(): + t = 2*G(m2)*G(m0)*G(m1)*G(-m0)*G(-m1) + r = match1_gamma(t, 1) + t = rule1_gamma(t, 1) + r = match1_gamma(t, 0) + t = rule1_gamma(t, 0) + assert t == (D*(-2*D + 4))*G(m2) + + t = 3*G(m2)*G(m0)*G(m1)*G(-m0) + t = rule1_gamma(t, 1) + + t = G(m2)*G(m0)*G(m1)*G(-m0)*G(-m2) + r = match1_gamma(t, 1) + assert r == (1, 3) + t = rule1_gamma(t, 1) + t = rule1_gamma(t, 1) + assert t == ((-D + 2)**2)*G(m1) + + t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1) + r = match1_gamma(t, 2) + assert r == (1, 4) + t = rule1_gamma(t, 2) + assert t == (D - 4)*G(m0)*G(m2)*G(m3) + 4*g(m2, m3)*G(m0) + + t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1)*G(-m0) + r = match1_gamma(t, 2) + t = rule1_gamma(t, 2) + + t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1)*G(-m0) + t = do_rule1_gamma(t) + t = do_rule1_gamma(t) + assert t == ((D - 4)**2)*G(m2)*G(m3) + (8*D - 16)*g(m2, m3) + + t = G(m2)*G(m0)*G(m1)*G(-m2)*G(-m0) + t = do_rule1_gamma(t) + t = do_rule1_gamma(t) + assert t == ((-D + 2)*(D - 4) + 4)*G(m1) + + t = G(m3)*G(m1)*G(m0)*G(m2)*G(-m3)*G(-m0)*G(-m2) + t = do_rule1_gamma(t) + t = do_rule1_gamma(t) + t = do_rule1_gamma(t) + assert t == (-4*D + (-D + 2)**2*(D - 4) + 8)*G(m1) + + M = Symbol('M') + + p = S1('p') + ps = p(m0)*G(-m0) + t0 = ps + M + t = G(m0)*(ps + M)*G(-m0) + t = do_rule1_gamma(t) + assert t == (2 - D)*ps + D*M + + t = G(m0)*(ps + M)*G(-m0)*(ps + M) + t = do_rule1_gamma(t) + + t = 2*G(m0)*G(m1)*G(m2)*G(m3)*G(-m0) + t = do_rule1_gamma(t) + assert t == (-2*D + 8)*G(m1)*G(m2)*G(m3) - 4*G(m3)*G(m2)*G(m1) + + t = G(m5)*G(m0)*G(m1)*G(m4)*G(m2)*G(-m4)*G(m3)*G(-m0) + t = do_rule1_gamma(t) + t = do_rule1_gamma(t) + assert t == ((-D + 2)*(-D + 4))*G(m5)*G(m1)*G(m2)*G(m3) + (2*D - 4)*G(m5)*G(m3)*G(m2)*G(m1) + + t = -G(m0)*G(m1)*G(m2)*G(m3)*G(-m0)*G(m4) + t = do_rule1_gamma(t) + assert t == (D - 4)*G(m1)*G(m2)*G(m3)*G(m4) + 2*G(m3)*G(m2)*G(m1)*G(m4) + + t = G(-m5)*G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)*G(m5) + t = do_rule1_gamma(t, doall=True) + assert t == ((-D + 4)**2 + 4)*G(m1)*G(m2)*G(m3)*G(m4) + (4*D - 16)*G(m3)*G(m2)*G(m1)*G(m4) + (4*D - 16)*G(m4)*G(m1)*G(m2)*G(m3) + 4*G(m2)*G(m1)*G(m4)*G(m3) + 4*G(m3)*G(m4)*G(m1)*G(m2) + 4*G(m4)*G(m3)*G(m2)*G(m1) + + +def test_rule2_gamma(): + p = S1('p') + M = Symbol('M') + ps = p(m0)*G(-m0) + t = ps*ps + t = t.canon_bp() + t = rule2_gamma(t, 0) + assert t == p(m0)*p(-m0) + + t = 2*G(m1)*ps*ps*G(m2) + t = rule2_gamma(t, 0) + assert t == 2*G(m1)*G(m2)*p(m0)*p(-m0) + + t0 = ps + M + t = G(m0)*(ps + M)*G(-m0) + t = do_rule1_gamma(t) + assert t == (2 - D)*ps + D*M + + t = G(m0)*(ps + M)*G(-m0)*(ps + M) + t = do_rule1_gamma(t) + t = do_rule2_gamma(t) + assert t == (-D + 2)*p(m0)*p(-m0) + (D*M + M*(-D + 2))*G(m0)*p(-m0) + D*M**2 + + +def test_gamma_trace(): + a = [G(m0), G(m1)] + t1 = gamma_trace1(*a) + a = [G(m0), G(m1), G(m2), G(m3)] + t1 = gamma_trace1(*a) + + t = G(m0)*G(m1) + t1 = gamma_trace(t) + assert t1 == 4*g(m0, m1) + t = G(m0)*G(m1)*G(m2)*G(m3) + t1 = gamma_trace(t) + t2 = -4*g(m0, m2)*g(m1, m3) + 4*g(m0, m1)*g(m2, m3) + 4*g(m0, m3)*g(m1, m2) + st2 = str(t2) + assert t1 == t2 + + t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0) + t1 = gamma_trace(t) + assert t1 == (-4*D)*g(m1, m3)*g(m2, m4) + (4*D)*g(m1, m2)*g(m3, m4) + \ + (4*D)*g(m1, m4)*g(m2, m3) + + t = G(-m5)*G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)*G(m5) + t1 = gamma_trace(t) + assert t1 == (32*D + 4*(-D + 4)**2 - 64)*(g(m1, m2)*g(m3, m4) - \ + g(m1, m3)*g(m2, m4) + g(m1, m4)*g(m2, m3)) + + t = G(m0)*G(m1)*G(-m0)*G(m3) + t1 = gamma_trace(t) + assert t1 == (-4*D + 8)*g(m1, m3) + + t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(m5) + t1 = gamma_trace(t) + t2 = t1*g(-m0, -m5) + t2 = t2.contract_metric(g) + assert t2 == D*gamma_trace(G(m1)*G(m2)*G(m3)*G(m4)) + + p, q = S1('p,q') + ps = p(m0)*G(-m0) + qs = q(m0)*G(-m0) + t = ps*qs*ps*qs + t1 = gamma_trace(t) + assert t1 == 8*p(m0)*q(-m0)*p(m1)*q(-m1) - 4*p(m0)*p(-m0)*q(m1)*q(-m1) + +def test_gamma5(): + G5 = GM.G5 + epsilon = GM.epsilon + t1 = GM.G5_to_right(2*G5) + assert t1 == 2*G5 + assert gamma_trace(t1) == 0 + + t = G5*G(m0) + t1 = GM.G5_to_right(t) + assert t1 == - G(m0)*G5 + + t = G5*G(m0)*G5*G(m1) + t1 = GM.G5_to_right(t) + assert t1 == -G(m0)*G(m1) + t = G5*G(m0)*G5*G(m1)*G5 + t1 = GM.G5_to_right(t) + assert t1 == -G(m0)*G(m1)*G5 + + t = (1 + G5)*(1 - G5) + t1 = GM.G5_to_right(t) + assert t1 == 0 + + t = G5*G(m0) + t1 = gamma_trace(t) + assert t1 == 0 + + t = G(m0)*G5*G(m1) + t1 = gamma_trace(t) + assert t1 == 0 + + t= G5*G(m0)*G(m1)*G(m2)*G(m3) + t1 = gamma_trace(t) + assert t1 == 4*I*epsilon(m0, m1, m2, m3) + + t= G5*G(m0)*G(-m0)*G(m1)*G(m2)*G(m3) + t1 = gamma_trace(t) + assert t1 == 0 + + t= G5*G(m0)*G(-m0)*G(m1)*G(m2)*G(m3)*G(m4) + t1 = gamma_trace(t) + assert t1 == 4*D*I*epsilon(m1, m2, m3, m4) + + t= G5*G(m0)*G(m1)*G(-m0)*G(m2)*G(m3)*G(m4) + t1 = gamma_trace(t) + assert t1 == (I*(-4*D + 8))*epsilon(m1, m2, m3, m4) + + t = G(m0)*(1 - G5) + t1 = gamma_trace(t) + assert t1 == 0 + + t = G(m0)*(1 - G5)*G(m1)*G(m2)*G(m3) + t1 = gamma_trace(t) + assert t1 == 4*I*epsilon(m0, m1, m2, m3) + 4*(g(m0,m1)*g(m2,m3) - \ + g(m0,m2)*g(m1,m3) + g(m0,m3)*g(m1,m2)) + + t = G(m0)*(1 - G5)*G(m1)*(1 + G5)*G(m2)*(1 - G5)*G(m3) + t1 = gamma_trace(t) + assert t1 == 16*I*epsilon(m0, m1, m2, m3) + 16*(g(m0,m1)*g(m2,m3) - \ + g(m0,m2)*g(m1,m3) + g(m0,m3)*g(m1,m2)) + + M = Symbol('M') + p, q = S1('p,q') + ps = p(m0)*G(-m0) + qs = q(m0)*G(-m0) + p2 = p(m0)*p(-m0) + q2 = q(m0)*q(-m0) + pq = p(m0)*q(-m0) + t = (ps + M)*G5*(qs + M)*G5 + t1 = gamma_trace(t) + assert t1 == -4*p(m0)*q(-m0) + 4*M**2 + + t = G(m0)*(ps + M)*G(m2)*(ps + qs + M)*G5 + t1 = gamma_trace(t) + assert t1 == -4*I*epsilon(m0,m2,m1,m3)*p(-m1)*q(-m3) + + t = G(m0)*(ps + M)*G(m1)*(ps + qs)*G(m2)*(1 - G5) + t1 = gamma_trace(t) + assert t1 == 4*M*(-g(m1, m2)*p(m0) -g(m1, m2)*q(m0) + \ + I*epsilon(m0, m1, m2, m3)*p(-m3) + I*epsilon(m0, m1, m2, m3)*q(-m3) + \ + g(m0, m1)*p(m2) + g(m0, m1)*q(m2) + g(m0, m2)*p(m1) + g(m0, m2)*q(m1)) + + t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(m5)*G5 + t1 = gamma_trace(t) + t2 = t1*g(-m0, -m2) + t2 = t2.contract_metric(g, contract_all=True) + assert t2 == (-4*I*D + 8*I)*epsilon(m1, m3, m4, m5) + +def test_trace1(): + M = Symbol('M') + p, q = S1('p,q') + ps = p(m0)*G(-m0) + qs = q(m0)*G(-m0) + p2 = p(m0)*p(-m0) + q2 = q(m0)*q(-m0) + pq = p(m0)*q(-m0) + + t = ps*qs*ps*qs*ps*qs + t1 = gamma_trace(t) + assert t1 == -12*p2*pq*q2 + 16*pq*pq*pq + + t = ps*qs*ps*qs*ps*qs*ps*qs + t1 = gamma_trace(t) + assert t1 == -32*pq*pq*p2*q2 + 32*pq*pq*pq*pq + 4*p2*p2*q2*q2 + +def test_epsilon(): + t = G(m0)*G(m1)*G(m2)*G(m3)*G(n0)*G(n1)*G(n2)*G(n3)*epsilon(-n0,-n1,-n2,-n3) + + t1 = gamma_trace(t) + assert t1 == 96*epsilon(m0, m1, m2, m3) diff --git a/sympy/tensor/tests/test_group_factors.py b/sympy/tensor/tests/test_group_factors.py new file mode 100644 index 000000000000..25ff790e64fe --- /dev/null +++ b/sympy/tensor/tests/test_group_factors.py @@ -0,0 +1,67 @@ +from time import time +from sympy import Symbol, S +from sympy.tensor.tensor import (tensor_indices) +from sympy.tensor.group_factors import SuNGroupFactors, match_CijCij + +""" +References + +[1] P. Cvitanovic "Group Theory" version 9.0.1 +""" + +def test_SuN(): + N = Symbol('N') + SN = SuNGroupFactors(N) + C = SN.C + g = SN.ASuN.metric + i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14 = \ + tensor_indices('i0:15', SN.ASuN) + theta = SN.theta + + t = C(i0,i1,i2)*C(-i0,-i2,i3) + t = SN.rule_C2theta_all(t) + assert t == -2*N*g(i1, i3) + + t = C(i0,i1,i2)*C(-i1,i3,i4)*C(i5,-i0,-i2)*C(-i5,-i4,-i3) + t = SN.rule_C2theta_all(t) + t = t.expand_coeff() + assert t == 4*N**4 - 4*N**2 + + # see eq.(1.1) and page 11 in Ref.[1] see also page 72 + # We apply rule_C_square to two parts of the tensor; + # applying rule_C_square once to the whole tensor it is 2x slower. + # Without applying rule_C_square it is roughly 50x slower. + t1 = C(i0,i1,i2)*C(-i1,i5,i3)*C(-i2,i4,i7)*C(-i3,-i4,i6) + t2 = C(-i6,i10,i11)*C(-i5,-i11,i8)*C(-i7,-i10,i9)*C(-i8,-i9,i12) + t1 = SN.rule_C_square(t1) + t2 = SN.rule_C_square(t2) + t = t1*t2 + t = SN.rule_C2theta_all(t) + t = t.expand_coeff() + assert t == (2*N**4 + 24*N**2)*g(i0, i12) + + # without using rule_C_triangle it is 2.5x slower + t = C(i0,i1,i2)*C(-i2,i3,i5)*C(-i1,-i3,i4)*C(-i4,-i5,i6) + t = SN.rule_C_triangle(t) + t = SN.rule_C2theta_all(t) + t = t.expand_coeff() + assert t == 2*N**2*g(i0, i6) + + # first graph in page 72 in Ref. [1] + t = theta(i0,i1,i2,i3)*theta(-i3,-i2,-i1,-i0) + t = SN.rule_C2theta_all(t) + t = t.expand_coeff() + assert (t._coeff - (N**4 - 3*N**2 + 3)*(N**2 - 1)/N**2).expand() == 0 + + + # Identically vanishing tensors + # contracted according to the Kuratowski graph eq.(6.59) in Ref. [1] + t = C(i0,i1,i2)*C(-i1,i3,i4)*C(-i3,i7,i5)*C(-i2,-i5,i6)*C(-i4,-i6,i8) + t1 = t.canon_bp() + assert t1 == 0 + + # contracted according to the Peterson graph eq.(6.60) in Ref. [1] + t = C(i0,i1,i2)*C(-i1,i3,i4)*C(-i2,i5,i6)*C(-i3,i7,i8)*C(-i6,-i7,i9)*\ + C(-i8,i10,i13)*C(-i5,-i10,i11)*C(-i4,-i11,i12)*C(-i9,-i12,i14) + t1 = t.canon_bp() + assert t1 == 0 diff --git a/sympy/tensor/tests/test_tensor.py b/sympy/tensor/tests/test_tensor.py new file mode 100644 index 000000000000..ab047803173b --- /dev/null +++ b/sympy/tensor/tests/test_tensor.py @@ -0,0 +1,1108 @@ +from sympy.core import S, Rational, Symbol, Basic +from sympy.combinatorics import Permutation +from sympy.combinatorics.tensor_can import (bsgs_direct_product, riemann_bsgs) +from sympy.tensor.tensor import (TensorIndexType, tensor_indices, + TensorSymmetry, get_symmetric_group_sgs, TensorType, TensorIndex, + tensor_mul, canon_bp, TensAdd, riemann_cyclic_replace, riemann_cyclic, + tensorlist_contract_metric, TensMul, tensorsymmetry, tensorhead, + TensorManager, HoldTensorHead, TensExpr) +from sympy.utilities.pytest import raises + +#################### Tests from tensor_can.py ####################### + +def test_canonicalize_no_slot_sym(): + # A_d0 * B^d0; T_c = A^d0*B_d0 + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + a, b, d0, d1 = tensor_indices('a,b,d0,d1', Lorentz) + sym1 = tensorsymmetry([1]) + S1 = TensorType([Lorentz], sym1) + A, B = S1('A,B') + t = A(-d0)*B(d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0)*B(-L_0)' + + # A^a * B^b; T_c = T + t = A(a)*B(b) + tc = t.canon_bp() + assert tc == t + # B^b * A^a + t1 = B(b)*A(a) + tc = t1.canon_bp() + assert str(tc) == 'A(a)*B(b)' + + # A symmetric + # A^{b}_{d0}*A^{d0, a}; T_c = A^{a d0}*A{b}_{d0} + sym2 = tensorsymmetry([1]*2) + S2 = TensorType([Lorentz]*2, sym2) + A = S2('A') + t = A(b, -d0)*A(d0, a) + tc = t.canon_bp() + assert str(tc) == 'A(a, L_0)*A(b, -L_0)' + + # A^{d1}_{d0}*B^d0*C_d1 + # T_c = A^{d0 d1}*B_d0*C_d1 + B, C = S1('B,C') + t = A(d1, -d0)*B(d0)*C(-d1) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-L_0)*C(-L_1)' + + # A without symmetry + # A^{d1}_{d0}*B^d0*C_d1 ord=[d0,-d0,d1,-d1]; g = [2,1,0,3,4,5] + # T_c = A^{d0 d1}*B_d1*C_d0; can = [0,2,3,1,4,5] + nsym2 = tensorsymmetry([1],[1]) + NS2 = TensorType([Lorentz]*2, nsym2) + A = NS2('A') + B, C = S1('B, C') + t = A(d1, -d0)*B(d0)*C(-d1) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-L_1)*C(-L_0)' + + # A, B without symmetry + # A^{d1}_{d0}*B_{d1}^{d0} + # T_c = A^{d0 d1}*B_{d0 d1} + B = NS2('B') + t = A(d1, -d0)*B(-d1, d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-L_0, -L_1)' + # A_{d0}^{d1}*B_{d1}^{d0} + # T_c = A^{d0 d1}*B_{d1 d0} + t = A(-d0, d1)*B(-d1, d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-L_1, -L_0)' + + # A, B, C without symmetry + # A^{d1 d0}*B_{a d0}*C_{d1 b} + # T_c=A^{d0 d1}*B_{a d1}*C_{d0 b} + C = NS2('C') + t = A(d1, d0)*B(-a, -d0)*C(-d1, -b) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-a, -L_1)*C(-L_0, -b)' + + # A symmetric, B and C without symmetry + # A^{d1 d0}*B_{a d0}*C_{d1 b} + # T_c = A^{d0 d1}*B_{a d0}*C_{d1 b} + A = S2('A') + t = A(d1, d0)*B(-a, -d0)*C(-d1, -b) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-a, -L_0)*C(-L_1, -b)' + + # A and C symmetric, B without symmetry + # A^{d1 d0}*B_{a d0}*C_{d1 b} ord=[a,b,d0,-d0,d1,-d1] + # T_c = A^{d0 d1}*B_{a d0}*C_{b d1} + C = S2('C') + t = A(d1, d0)*B(-a, -d0)*C(-d1, -b) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1)*B(-a, -L_0)*C(-b, -L_1)' + +def test_canonicalize_no_dummies(): + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + a, b, c, d = tensor_indices('a, b, c, d', Lorentz) + sym1 = tensorsymmetry([1]) + sym2 = tensorsymmetry([1]*2) + sym2a = tensorsymmetry([2]) + + # A commuting + # A^c A^b A^a + # T_c = A^a A^b A^c + S1 = TensorType([Lorentz], sym1) + A = S1('A') + t = A(c)*A(b)*A(a) + tc = t.canon_bp() + assert str(tc) == 'A(a)*A(b)*A(c)' + + # A anticommuting + # A^c A^b A^a + # T_c = -A^a A^b A^c + A = S1('A', 1) + t = A(c)*A(b)*A(a) + tc = t.canon_bp() + assert str(tc) == '-A(a)*A(b)*A(c)' + + # A commuting and symmetric + # A^{b,d}*A^{c,a} + # T_c = A^{a c}*A^{b d} + S2 = TensorType([Lorentz]*2, sym2) + A = S2('A') + t = A(b, d)*A(c, a) + tc = t.canon_bp() + assert str(tc) == 'A(a, c)*A(b, d)' + + # A anticommuting and symmetric + # A^{b,d}*A^{c,a} + # T_c = -A^{a c}*A^{b d} + A = S2('A', 1) + t = A(b, d)*A(c, a) + tc = t.canon_bp() + assert str(tc) == '-A(a, c)*A(b, d)' + + # A^{c,a}*A^{b,d} + # T_c = A^{a c}*A^{b d} + t = A(c, a)*A(b, d) + tc = t.canon_bp() + assert str(tc) == 'A(a, c)*A(b, d)' + +def test_no_metric_symmetry(): + # no metric symmetry; A no symmetry + # A^d1_d0 * A^d0_d1 + # T_c = A^d0_d1 * A^d1_d0 + Lorentz = TensorIndexType('Lorentz', metric=None, dummy_fmt='L') + d0, d1, d2, d3 = tensor_indices('d:4', Lorentz) + A = tensorhead('A', [Lorentz]*2, [[1], [1]]) + t = A(d1, -d0)*A(d0, -d1) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, -L_1)*A(L_1, -L_0)' + + # A^d1_d2 * A^d0_d3 * A^d2_d1 * A^d3_d0 + # T_c = A^d0_d1 * A^d1_d0 * A^d2_d3 * A^d3_d2 + t = A(d1, -d2)*A(d0, -d3)*A(d2,-d1)*A(d3,-d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, -L_1)*A(L_1, -L_0)*A(L_2, -L_3)*A(L_3, -L_2)' + + # A^d0_d2 * A^d1_d3 * A^d3_d0 * A^d2_d1 + # T_c = A^d0_d1 * A^d1_d2 * A^d2_d3 * A^d3_d0 + t = A(d0, -d1)*A(d1, -d2)*A(d2, -d3)*A(d3,-d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, -L_1)*A(L_1, -L_2)*A(L_2, -L_3)*A(L_3, -L_0)' + +def test_canonicalize1(): + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + a, a0, a1, a2, a3, b, d0, d1, d2, d3 = \ + tensor_indices('a,a0,a1,a2,a3,b,d0,d1,d2,d3', Lorentz) + sym1 = tensorsymmetry([1]) + base3, gens3 = get_symmetric_group_sgs(3) + sym2 = tensorsymmetry([1]*2) + sym2a = tensorsymmetry([2]) + sym3 = tensorsymmetry([1]*3) + sym3a = tensorsymmetry([3]) + + # A_d0*A^d0; ord = [d0,-d0] + # T_c = A^d0*A_d0 + S1 = TensorType([Lorentz], sym1) + A = S1('A') + t = A(-d0)*A(d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0)*A(-L_0)' + + # A commuting + # A_d0*A_d1*A_d2*A^d2*A^d1*A^d0 + # T_c = A^d0*A_d0*A^d1*A_d1*A^d2*A_d2 + t = A(-d0)*A(-d1)*A(-d2)*A(d2)*A(d1)*A(d0) + tc = t.canon_bp() + assert str(tc) == 'A(L_0)*A(-L_0)*A(L_1)*A(-L_1)*A(L_2)*A(-L_2)' + + # A anticommuting + # A_d0*A_d1*A_d2*A^d2*A^d1*A^d0 + # T_c 0 + A = S1('A', 1) + t = A(-d0)*A(-d1)*A(-d2)*A(d2)*A(d1)*A(d0) + tc = t.canon_bp() + assert tc == 0 + + # A commuting symmetric + # A^{d0 b}*A^a_d1*A^d1_d0 + # T_c = A^{a d0}*A^{b d1}*A_{d0 d1} + S2 = TensorType([Lorentz]*2, sym2) + A = S2('A') + t = A(d0, b)*A(a, -d1)*A(d1, -d0) + tc = t.canon_bp() + assert str(tc) == 'A(a, L_0)*A(b, L_1)*A(-L_0, -L_1)' + + # A, B commuting symmetric + # A^{d0 b}*A^d1_d0*B^a_d1 + # T_c = A^{b d0}*A_d0^d1*B^a_d1 + B = S2('B') + t = A(d0, b)*A(d1, -d0)*B(a, -d1) + tc = t.canon_bp() + assert str(tc) == 'A(b, L_0)*A(-L_0, L_1)*B(a, -L_1)' + + # A commuting symmetric + # A^{d1 d0 b}*A^{a}_{d1 d0}; ord=[a,b, d0,-d0,d1,-d1] + # T_c = A^{a d0 d1}*A^{b}_{d0 d1} + S3 = TensorType([Lorentz]*3, sym3) + A = S3('A') + t = A(d1, d0, b)*A(a, -d1, -d0) + tc = t.canon_bp() + assert str(tc) == 'A(a, L_0, L_1)*A(b, -L_0, -L_1)' + + # A^{d3 d0 d2}*A^a0_{d1 d2}*A^d1_d3^a1*A^{a2 a3}_d0 + # T_c = A^{a0 d0 d1}*A^a1_d0^d2*A^{a2 a3 d3}*A_{d1 d2 d3} + t = A(d3, d0, d2)*A(a0, -d1, -d2)*A(d1, -d3, a1)*A(a2, a3, -d0) + tc = t.canon_bp() + assert str(tc) == 'A(a0, L_0, L_1)*A(a1, -L_0, L_2)*A(a2, a3, L_3)*A(-L_1, -L_2, -L_3)' + + # A commuting symmetric, B antisymmetric + # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 + # in this esxample and in the next three, + # renaming dummy indices and using symmetry of A, + # T = A^{d0 d1 d2} * A_{d0 d1 d3} * B_d2^d3 + # can = 0 + S2a = TensorType([Lorentz]*2, sym2a) + A = S3('A') + B = S2a('B') + t = A(d0, d1, d2)*A(-d2, -d3, -d1)*B(-d0, d3) + tc = t.canon_bp() + assert tc == 0 + + # A anticommuting symmetric, B anticommuting + # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 + # T_c = A^{d0 d1 d2} * A_{d0 d1}^d3 * B_{d2 d3} + A = S3('A', 1) + B = S2a('B') + t = A(d0, d1, d2)*A(-d2, -d3, -d1)*B(-d0, d3) + tc = t.canon_bp() + assert str(tc) == 'A(L_0, L_1, L_2)*A(-L_0, -L_1, L_3)*B(-L_2, -L_3)' + + # A anticommuting symmetric, B antisymmetric commuting, antisymmetric metric + # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 + # T_c = -A^{d0 d1 d2} * A_{d0 d1}^d3 * B_{d2 d3} + Spinor = TensorIndexType('Spinor', metric=1, dummy_fmt='S') + a, a0, a1, a2, a3, b, d0, d1, d2, d3 = \ + tensor_indices('a,a0,a1,a2,a3,b,d0,d1,d2,d3', Spinor) + S3 = TensorType([Spinor]*3, sym3) + S2a = TensorType([Spinor]*2, sym2a) + A = S3('A', 1) + B = S2a('B') + t = A(d0, d1, d2)*A(-d2, -d3, -d1)*B(-d0, d3) + tc = t.canon_bp() + assert str(tc) == '-A(S_0, S_1, S_2)*A(-S_0, -S_1, S_3)*B(-S_2, -S_3)' + + # A anticommuting symmetric, B antisymmetric anticommuting, + # no metric symmetry + # A^{d0 d1 d2} * A_{d2 d3 d1} * B_d0^d3 + # T_c = A^{d0 d1 d2} * A_{d0 d1 d3} * B_d2^d3 + Mat = TensorIndexType('Mat', metric=None, dummy_fmt='M') + a, a0, a1, a2, a3, b, d0, d1, d2, d3 = \ + tensor_indices('a,a0,a1,a2,a3,b,d0,d1,d2,d3', Mat) + S3 = TensorType([Mat]*3, sym3) + S2a = TensorType([Mat]*2, sym2a) + A = S3('A', 1) + B = S2a('B') + t = A(d0, d1, d2)*A(-d2, -d3, -d1)*B(-d0, d3) + tc = t.canon_bp() + assert str(tc) == 'A(M_0, M_1, M_2)*A(-M_0, -M_1, -M_3)*B(-M_2, M_3)' + + # Gamma anticommuting + # Gamma_{mu nu} * gamma^rho * Gamma^{nu mu alpha} + # T_c = -Gamma^{mu nu} * gamma^rho * Gamma_{alpha mu nu} + S1 = TensorType([Lorentz], sym1) + S2a = TensorType([Lorentz]*2, sym2a) + S3a = TensorType([Lorentz]*3, sym3a) + alpha, beta, gamma, mu, nu, rho = \ + tensor_indices('alpha,beta,gamma,mu,nu,rho', Lorentz) + Gamma = S1('Gamma', 2) + Gamma2 = S2a('Gamma', 2) + Gamma3 = S3a('Gamma', 2) + t = Gamma2(-mu,-nu)*Gamma(rho)*Gamma3(nu, mu, alpha) + tc = t.canon_bp() + assert str(tc) == '-Gamma(L_0, L_1)*Gamma(rho)*Gamma(alpha, -L_0, -L_1)' + + # Gamma_{mu nu} * Gamma^{gamma beta} * gamma_rho * Gamma^{nu mu alpha} + # T_c = Gamma^{mu nu} * Gamma^{beta gamma} * gamma_rho * Gamma^alpha_{mu nu} + t = Gamma2(mu, nu)*Gamma2(beta, gamma)*Gamma(-rho)*Gamma3(alpha, -mu, -nu) + tc = t.canon_bp() + assert str(tc) == 'Gamma(L_0, L_1)*Gamma(beta, gamma)*Gamma(-rho)*Gamma(alpha, -L_0, -L_1)' + + # f^a_{b,c} antisymmetric in b,c; A_mu^a no symmetry + # f^c_{d a} * f_{c e b} * A_mu^d * A_nu^a * A^{nu e} * A^{mu b} + # g = [8,11,5, 9,13,7, 1,10, 3,4, 2,12, 0,6, 14,15] + # T_c = -f^{a b c} * f_a^{d e} * A^mu_b * A_{mu d} * A^nu_c * A_{nu e} + Flavor = TensorIndexType('Flavor', dummy_fmt='F') + a, b, c, d, e, ff = tensor_indices('a,b,c,d,e,f', Flavor) + mu, nu = tensor_indices('mu,nu', Lorentz) + sym_f = tensorsymmetry([1], [2]) + S_f = TensorType([Flavor]*3, sym_f) + sym_A = tensorsymmetry([1], [1]) + S_A = TensorType([Lorentz, Flavor], sym_A) + f = S_f('f') + A = S_A('A') + t = f(c, -d, -a)*f(-c, -e, -b)*A(-mu, d)*A(-nu, a)*A(nu, e)*A(mu, b) + tc = t.canon_bp() + assert str(tc) == '-f(F_0, F_1, F_2)*f(-F_0, F_3, F_4)*A(L_0, -F_1)*A(-L_0, -F_3)*A(L_1, -F_2)*A(-L_1, -F_4)' + + +def test_riemann_invariants(): + Lorentz = TensorIndexType('Lorentz', dummy_fmt='L') + d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11 = \ + tensor_indices('d0:12', Lorentz) + # R^{d0 d1}_{d1 d0}; ord = [d0,-d0,d1,-d1] + # T_c = -R^{d0 d1}_{d0 d1} + R = tensorhead('R', [Lorentz]*4, [[2, 2]]) + t = R(d0, d1, -d1, -d0) + tc = t.canon_bp() + assert str(tc) == '-R(L_0, L_1, -L_0, -L_1)' + + # R_d11^d1_d0^d5 * R^{d6 d4 d0}_d5 * R_{d7 d2 d8 d9} * + # R_{d10 d3 d6 d4} * R^{d2 d7 d11}_d1 * R^{d8 d9 d3 d10} + # can = [0,2,4,6, 1,3,8,10, 5,7,12,14, 9,11,16,18, 13,15,20,22, + # 17,19,21