Coverage for /usr/lib/python3/dist-packages/sympy/tensor/tensor.py: 20%
2495 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2This module defines tensors with abstract index notation.
4The abstract index notation has been first formalized by Penrose.
6Tensor indices are formal objects, with a tensor type; there is no
7notion of index range, it is only possible to assign the dimension,
8used to trace the Kronecker delta; the dimension can be a Symbol.
10The Einstein summation convention is used.
11The covariant indices are indicated with a minus sign in front of the index.
13For instance the tensor ``t = p(a)*A(b,c)*q(-c)`` has the index ``c``
14contracted.
16A tensor expression ``t`` can be called; called with its
17indices in sorted order it is equal to itself:
18in the above example ``t(a, b) == t``;
19one can call ``t`` with different indices; ``t(c, d) == p(c)*A(d,a)*q(-a)``.
21The contracted indices are dummy indices, internally they have no name,
22the indices being represented by a graph-like structure.
24Tensors are put in canonical form using ``canon_bp``, which uses
25the Butler-Portugal algorithm for canonicalization using the monoterm
26symmetries of the tensors.
28If there is a (anti)symmetric metric, the indices can be raised and
29lowered when the tensor is put in canonical form.
30"""
32from __future__ import annotations
33from typing import Any
34from functools import reduce
35from math import prod
37from abc import abstractmethod, ABC
38from collections import defaultdict
39import operator
40import itertools
41from sympy.core.numbers import (Integer, Rational)
42from sympy.combinatorics import Permutation
43from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, \
44 bsgs_direct_product, canonicalize, riemann_bsgs
45from sympy.core import Basic, Expr, sympify, Add, Mul, S
46from sympy.core.containers import Tuple, Dict
47from sympy.core.sorting import default_sort_key
48from sympy.core.symbol import Symbol, symbols
49from sympy.core.sympify import CantSympify, _sympify
50from sympy.core.operations import AssocOp
51from sympy.external.gmpy import SYMPY_INTS
52from sympy.matrices import eye
53from sympy.utilities.exceptions import (sympy_deprecation_warning,
54 SymPyDeprecationWarning,
55 ignore_warnings)
56from sympy.utilities.decorator import memoize_property, deprecated
57from sympy.utilities.iterables import sift
60def deprecate_data():
61 sympy_deprecation_warning(
62 """
63 The data attribute of TensorIndexType is deprecated. Use The
64 replace_with_arrays() method instead.
65 """,
66 deprecated_since_version="1.4",
67 active_deprecations_target="deprecated-tensorindextype-attrs",
68 stacklevel=4,
69 )
71def deprecate_fun_eval():
72 sympy_deprecation_warning(
73 """
74 The Tensor.fun_eval() method is deprecated. Use
75 Tensor.substitute_indices() instead.
76 """,
77 deprecated_since_version="1.5",
78 active_deprecations_target="deprecated-tensor-fun-eval",
79 stacklevel=4,
80 )
83def deprecate_call():
84 sympy_deprecation_warning(
85 """
86 Calling a tensor like Tensor(*indices) is deprecated. Use
87 Tensor.substitute_indices() instead.
88 """,
89 deprecated_since_version="1.5",
90 active_deprecations_target="deprecated-tensor-fun-eval",
91 stacklevel=4,
92 )
95class _IndexStructure(CantSympify):
96 """
97 This class handles the indices (free and dummy ones). It contains the
98 algorithms to manage the dummy indices replacements and contractions of
99 free indices under multiplications of tensor expressions, as well as stuff
100 related to canonicalization sorting, getting the permutation of the
101 expression and so on. It also includes tools to get the ``TensorIndex``
102 objects corresponding to the given index structure.
103 """
105 def __init__(self, free, dum, index_types, indices, canon_bp=False):
106 self.free = free
107 self.dum = dum
108 self.index_types = index_types
109 self.indices = indices
110 self._ext_rank = len(self.free) + 2*len(self.dum)
111 self.dum.sort(key=lambda x: x[0])
113 @staticmethod
114 def from_indices(*indices):
115 """
116 Create a new ``_IndexStructure`` object from a list of ``indices``.
118 Explanation
119 ===========
121 ``indices`` ``TensorIndex`` objects, the indices. Contractions are
122 detected upon construction.
124 Examples
125 ========
127 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, _IndexStructure
128 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
129 >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz)
130 >>> _IndexStructure.from_indices(m0, m1, -m1, m3)
131 _IndexStructure([(m0, 0), (m3, 3)], [(1, 2)], [Lorentz, Lorentz, Lorentz, Lorentz])
132 """
134 free, dum = _IndexStructure._free_dum_from_indices(*indices)
135 index_types = [i.tensor_index_type for i in indices]
136 indices = _IndexStructure._replace_dummy_names(indices, free, dum)
137 return _IndexStructure(free, dum, index_types, indices)
139 @staticmethod
140 def from_components_free_dum(components, free, dum):
141 index_types = []
142 for component in components:
143 index_types.extend(component.index_types)
144 indices = _IndexStructure.generate_indices_from_free_dum_index_types(free, dum, index_types)
145 return _IndexStructure(free, dum, index_types, indices)
147 @staticmethod
148 def _free_dum_from_indices(*indices):
149 """
150 Convert ``indices`` into ``free``, ``dum`` for single component tensor.
152 Explanation
153 ===========
155 ``free`` list of tuples ``(index, pos, 0)``,
156 where ``pos`` is the position of index in
157 the list of indices formed by the component tensors
159 ``dum`` list of tuples ``(pos_contr, pos_cov, 0, 0)``
161 Examples
162 ========
164 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, \
165 _IndexStructure
166 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
167 >>> m0, m1, m2, m3 = tensor_indices('m0,m1,m2,m3', Lorentz)
168 >>> _IndexStructure._free_dum_from_indices(m0, m1, -m1, m3)
169 ([(m0, 0), (m3, 3)], [(1, 2)])
170 """
171 n = len(indices)
172 if n == 1:
173 return [(indices[0], 0)], []
175 # find the positions of the free indices and of the dummy indices
176 free = [True]*len(indices)
177 index_dict = {}
178 dum = []
179 for i, index in enumerate(indices):
180 name = index.name
181 typ = index.tensor_index_type
182 contr = index.is_up
183 if (name, typ) in index_dict:
184 # found a pair of dummy indices
185 is_contr, pos = index_dict[(name, typ)]
186 # check consistency and update free
187 if is_contr:
188 if contr:
189 raise ValueError('two equal contravariant indices in slots %d and %d' %(pos, i))
190 else:
191 free[pos] = False
192 free[i] = False
193 else:
194 if contr:
195 free[pos] = False
196 free[i] = False
197 else:
198 raise ValueError('two equal covariant indices in slots %d and %d' %(pos, i))
199 if contr:
200 dum.append((i, pos))
201 else:
202 dum.append((pos, i))
203 else:
204 index_dict[(name, typ)] = index.is_up, i
206 free = [(index, i) for i, index in enumerate(indices) if free[i]]
207 free.sort()
208 return free, dum
210 def get_indices(self):
211 """
212 Get a list of indices, creating new tensor indices to complete dummy indices.
213 """
214 return self.indices[:]
216 @staticmethod
217 def generate_indices_from_free_dum_index_types(free, dum, index_types):
218 indices = [None]*(len(free)+2*len(dum))
219 for idx, pos in free:
220 indices[pos] = idx
222 generate_dummy_name = _IndexStructure._get_generator_for_dummy_indices(free)
223 for pos1, pos2 in dum:
224 typ1 = index_types[pos1]
225 indname = generate_dummy_name(typ1)
226 indices[pos1] = TensorIndex(indname, typ1, True)
227 indices[pos2] = TensorIndex(indname, typ1, False)
229 return _IndexStructure._replace_dummy_names(indices, free, dum)
231 @staticmethod
232 def _get_generator_for_dummy_indices(free):
233 cdt = defaultdict(int)
234 # if the free indices have names with dummy_name, start with an
235 # index higher than those for the dummy indices
236 # to avoid name collisions
237 for indx, ipos in free:
238 if indx.name.split('_')[0] == indx.tensor_index_type.dummy_name:
239 cdt[indx.tensor_index_type] = max(cdt[indx.tensor_index_type], int(indx.name.split('_')[1]) + 1)
241 def dummy_name_gen(tensor_index_type):
242 nd = str(cdt[tensor_index_type])
243 cdt[tensor_index_type] += 1
244 return tensor_index_type.dummy_name + '_' + nd
246 return dummy_name_gen
248 @staticmethod
249 def _replace_dummy_names(indices, free, dum):
250 dum.sort(key=lambda x: x[0])
251 new_indices = list(indices)
252 assert len(indices) == len(free) + 2*len(dum)
253 generate_dummy_name = _IndexStructure._get_generator_for_dummy_indices(free)
254 for ipos1, ipos2 in dum:
255 typ1 = new_indices[ipos1].tensor_index_type
256 indname = generate_dummy_name(typ1)
257 new_indices[ipos1] = TensorIndex(indname, typ1, True)
258 new_indices[ipos2] = TensorIndex(indname, typ1, False)
259 return new_indices
261 def get_free_indices(self) -> list[TensorIndex]:
262 """
263 Get a list of free indices.
264 """
265 # get sorted indices according to their position:
266 free = sorted(self.free, key=lambda x: x[1])
267 return [i[0] for i in free]
269 def __str__(self):
270 return "_IndexStructure({}, {}, {})".format(self.free, self.dum, self.index_types)
272 def __repr__(self):
273 return self.__str__()
275 def _get_sorted_free_indices_for_canon(self):
276 sorted_free = self.free[:]
277 sorted_free.sort(key=lambda x: x[0])
278 return sorted_free
280 def _get_sorted_dum_indices_for_canon(self):
281 return sorted(self.dum, key=lambda x: x[0])
283 def _get_lexicographically_sorted_index_types(self):
284 permutation = self.indices_canon_args()[0]
285 index_types = [None]*self._ext_rank
286 for i, it in enumerate(self.index_types):
287 index_types[permutation(i)] = it
288 return index_types
290 def _get_lexicographically_sorted_indices(self):
291 permutation = self.indices_canon_args()[0]
292 indices = [None]*self._ext_rank
293 for i, it in enumerate(self.indices):
294 indices[permutation(i)] = it
295 return indices
297 def perm2tensor(self, g, is_canon_bp=False):
298 """
299 Returns a ``_IndexStructure`` instance corresponding to the permutation ``g``.
301 Explanation
302 ===========
304 ``g`` permutation corresponding to the tensor in the representation
305 used in canonicalization
307 ``is_canon_bp`` if True, then ``g`` is the permutation
308 corresponding to the canonical form of the tensor
309 """
310 sorted_free = [i[0] for i in self._get_sorted_free_indices_for_canon()]
311 lex_index_types = self._get_lexicographically_sorted_index_types()
312 lex_indices = self._get_lexicographically_sorted_indices()
313 nfree = len(sorted_free)
314 rank = self._ext_rank
315 dum = [[None]*2 for i in range((rank - nfree)//2)]
316 free = []
318 index_types = [None]*rank
319 indices = [None]*rank
320 for i in range(rank):
321 gi = g[i]
322 index_types[i] = lex_index_types[gi]
323 indices[i] = lex_indices[gi]
324 if gi < nfree:
325 ind = sorted_free[gi]
326 assert index_types[i] == sorted_free[gi].tensor_index_type
327 free.append((ind, i))
328 else:
329 j = gi - nfree
330 idum, cov = divmod(j, 2)
331 if cov:
332 dum[idum][1] = i
333 else:
334 dum[idum][0] = i
335 dum = [tuple(x) for x in dum]
337 return _IndexStructure(free, dum, index_types, indices)
339 def indices_canon_args(self):
340 """
341 Returns ``(g, dummies, msym, v)``, the entries of ``canonicalize``
343 See ``canonicalize`` in ``tensor_can.py`` in combinatorics module.
344 """
345 # to be called after sorted_components
346 from sympy.combinatorics.permutations import _af_new
347 n = self._ext_rank
348 g = [None]*n + [n, n+1]
350 # Converts the symmetry of the metric into msym from .canonicalize()
351 # method in the combinatorics module
352 def metric_symmetry_to_msym(metric):
353 if metric is None:
354 return None
355 sym = metric.symmetry
356 if sym == TensorSymmetry.fully_symmetric(2):
357 return 0
358 if sym == TensorSymmetry.fully_symmetric(-2):
359 return 1
360 return None
362 # ordered indices: first the free indices, ordered by types
363 # then the dummy indices, ordered by types and contravariant before
364 # covariant
365 # g[position in tensor] = position in ordered indices
366 for i, (indx, ipos) in enumerate(self._get_sorted_free_indices_for_canon()):
367 g[ipos] = i
368 pos = len(self.free)
369 j = len(self.free)
370 dummies = []
371 prev = None
372 a = []
373 msym = []
374 for ipos1, ipos2 in self._get_sorted_dum_indices_for_canon():
375 g[ipos1] = j
376 g[ipos2] = j + 1
377 j += 2
378 typ = self.index_types[ipos1]
379 if typ != prev:
380 if a:
381 dummies.append(a)
382 a = [pos, pos + 1]
383 prev = typ
384 msym.append(metric_symmetry_to_msym(typ.metric))
385 else:
386 a.extend([pos, pos + 1])
387 pos += 2
388 if a:
389 dummies.append(a)
391 return _af_new(g), dummies, msym
394def components_canon_args(components):
395 numtyp = []
396 prev = None
397 for t in components:
398 if t == prev:
399 numtyp[-1][1] += 1
400 else:
401 prev = t
402 numtyp.append([prev, 1])
403 v = []
404 for h, n in numtyp:
405 if h.comm in (0, 1):
406 comm = h.comm
407 else:
408 comm = TensorManager.get_comm(h.comm, h.comm)
409 v.append((h.symmetry.base, h.symmetry.generators, n, comm))
410 return v
413class _TensorDataLazyEvaluator(CantSympify):
414 """
415 EXPERIMENTAL: do not rely on this class, it may change without deprecation
416 warnings in future versions of SymPy.
418 Explanation
419 ===========
421 This object contains the logic to associate components data to a tensor
422 expression. Components data are set via the ``.data`` property of tensor
423 expressions, is stored inside this class as a mapping between the tensor
424 expression and the ``ndarray``.
426 Computations are executed lazily: whereas the tensor expressions can have
427 contractions, tensor products, and additions, components data are not
428 computed until they are accessed by reading the ``.data`` property
429 associated to the tensor expression.
430 """
431 _substitutions_dict: dict[Any, Any] = {}
432 _substitutions_dict_tensmul: dict[Any, Any] = {}
434 def __getitem__(self, key):
435 dat = self._get(key)
436 if dat is None:
437 return None
439 from .array import NDimArray
440 if not isinstance(dat, NDimArray):
441 return dat
443 if dat.rank() == 0:
444 return dat[()]
445 elif dat.rank() == 1 and len(dat) == 1:
446 return dat[0]
447 return dat
449 def _get(self, key):
450 """
451 Retrieve ``data`` associated with ``key``.
453 Explanation
454 ===========
456 This algorithm looks into ``self._substitutions_dict`` for all
457 ``TensorHead`` in the ``TensExpr`` (or just ``TensorHead`` if key is a
458 TensorHead instance). It reconstructs the components data that the
459 tensor expression should have by performing on components data the
460 operations that correspond to the abstract tensor operations applied.
462 Metric tensor is handled in a different manner: it is pre-computed in
463 ``self._substitutions_dict_tensmul``.
464 """
465 if key in self._substitutions_dict:
466 return self._substitutions_dict[key]
468 if isinstance(key, TensorHead):
469 return None
471 if isinstance(key, Tensor):
472 # special case to handle metrics. Metric tensors cannot be
473 # constructed through contraction by the metric, their
474 # components show if they are a matrix or its inverse.
475 signature = tuple([i.is_up for i in key.get_indices()])
476 srch = (key.component,) + signature
477 if srch in self._substitutions_dict_tensmul:
478 return self._substitutions_dict_tensmul[srch]
479 array_list = [self.data_from_tensor(key)]
480 return self.data_contract_dum(array_list, key.dum, key.ext_rank)
482 if isinstance(key, TensMul):
483 tensmul_args = key.args
484 if len(tensmul_args) == 1 and len(tensmul_args[0].components) == 1:
485 # special case to handle metrics. Metric tensors cannot be
486 # constructed through contraction by the metric, their
487 # components show if they are a matrix or its inverse.
488 signature = tuple([i.is_up for i in tensmul_args[0].get_indices()])
489 srch = (tensmul_args[0].components[0],) + signature
490 if srch in self._substitutions_dict_tensmul:
491 return self._substitutions_dict_tensmul[srch]
492 #data_list = [self.data_from_tensor(i) for i in tensmul_args if isinstance(i, TensExpr)]
493 data_list = [self.data_from_tensor(i) if isinstance(i, Tensor) else i.data for i in tensmul_args if isinstance(i, TensExpr)]
494 coeff = prod([i for i in tensmul_args if not isinstance(i, TensExpr)])
495 if all(i is None for i in data_list):
496 return None
497 if any(i is None for i in data_list):
498 raise ValueError("Mixing tensors with associated components "\
499 "data with tensors without components data")
500 data_result = self.data_contract_dum(data_list, key.dum, key.ext_rank)
501 return coeff*data_result
503 if isinstance(key, TensAdd):
504 data_list = []
505 free_args_list = []
506 for arg in key.args:
507 if isinstance(arg, TensExpr):
508 data_list.append(arg.data)
509 free_args_list.append([x[0] for x in arg.free])
510 else:
511 data_list.append(arg)
512 free_args_list.append([])
513 if all(i is None for i in data_list):
514 return None
515 if any(i is None for i in data_list):
516 raise ValueError("Mixing tensors with associated components "\
517 "data with tensors without components data")
519 sum_list = []
520 from .array import permutedims
521 for data, free_args in zip(data_list, free_args_list):
522 if len(free_args) < 2:
523 sum_list.append(data)
524 else:
525 free_args_pos = {y: x for x, y in enumerate(free_args)}
526 axes = [free_args_pos[arg] for arg in key.free_args]
527 sum_list.append(permutedims(data, axes))
528 return reduce(lambda x, y: x+y, sum_list)
530 return None
532 @staticmethod
533 def data_contract_dum(ndarray_list, dum, ext_rank):
534 from .array import tensorproduct, tensorcontraction, MutableDenseNDimArray
535 arrays = list(map(MutableDenseNDimArray, ndarray_list))
536 prodarr = tensorproduct(*arrays)
537 return tensorcontraction(prodarr, *dum)
539 def data_tensorhead_from_tensmul(self, data, tensmul, tensorhead):
540 """
541 This method is used when assigning components data to a ``TensMul``
542 object, it converts components data to a fully contravariant ndarray,
543 which is then stored according to the ``TensorHead`` key.
544 """
545 if data is None:
546 return None
548 return self._correct_signature_from_indices(
549 data,
550 tensmul.get_indices(),
551 tensmul.free,
552 tensmul.dum,
553 True)
555 def data_from_tensor(self, tensor):
556 """
557 This method corrects the components data to the right signature
558 (covariant/contravariant) using the metric associated with each
559 ``TensorIndexType``.
560 """
561 tensorhead = tensor.component
563 if tensorhead.data is None:
564 return None
566 return self._correct_signature_from_indices(
567 tensorhead.data,
568 tensor.get_indices(),
569 tensor.free,
570 tensor.dum)
572 def _assign_data_to_tensor_expr(self, key, data):
573 if isinstance(key, TensAdd):
574 raise ValueError('cannot assign data to TensAdd')
575 # here it is assumed that `key` is a `TensMul` instance.
576 if len(key.components) != 1:
577 raise ValueError('cannot assign data to TensMul with multiple components')
578 tensorhead = key.components[0]
579 newdata = self.data_tensorhead_from_tensmul(data, key, tensorhead)
580 return tensorhead, newdata
582 def _check_permutations_on_data(self, tens, data):
583 from .array import permutedims
584 from .array.arrayop import Flatten
586 if isinstance(tens, TensorHead):
587 rank = tens.rank
588 generators = tens.symmetry.generators
589 elif isinstance(tens, Tensor):
590 rank = tens.rank
591 generators = tens.components[0].symmetry.generators
592 elif isinstance(tens, TensorIndexType):
593 rank = tens.metric.rank
594 generators = tens.metric.symmetry.generators
596 # Every generator is a permutation, check that by permuting the array
597 # by that permutation, the array will be the same, except for a
598 # possible sign change if the permutation admits it.
599 for gener in generators:
600 sign_change = +1 if (gener(rank) == rank) else -1
601 data_swapped = data
602 last_data = data
603 permute_axes = list(map(gener, range(rank)))
604 # the order of a permutation is the number of times to get the
605 # identity by applying that permutation.
606 for i in range(gener.order()-1):
607 data_swapped = permutedims(data_swapped, permute_axes)
608 # if any value in the difference array is non-zero, raise an error:
609 if any(Flatten(last_data - sign_change*data_swapped)):
610 raise ValueError("Component data symmetry structure error")
611 last_data = data_swapped
613 def __setitem__(self, key, value):
614 """
615 Set the components data of a tensor object/expression.
617 Explanation
618 ===========
620 Components data are transformed to the all-contravariant form and stored
621 with the corresponding ``TensorHead`` object. If a ``TensorHead`` object
622 cannot be uniquely identified, it will raise an error.
623 """
624 data = _TensorDataLazyEvaluator.parse_data(value)
625 self._check_permutations_on_data(key, data)
627 # TensorHead and TensorIndexType can be assigned data directly, while
628 # TensMul must first convert data to a fully contravariant form, and
629 # assign it to its corresponding TensorHead single component.
630 if not isinstance(key, (TensorHead, TensorIndexType)):
631 key, data = self._assign_data_to_tensor_expr(key, data)
633 if isinstance(key, TensorHead):
634 for dim, indextype in zip(data.shape, key.index_types):
635 if indextype.data is None:
636 raise ValueError("index type {} has no components data"\
637 " associated (needed to raise/lower index)".format(indextype))
638 if not indextype.dim.is_number:
639 continue
640 if dim != indextype.dim:
641 raise ValueError("wrong dimension of ndarray")
642 self._substitutions_dict[key] = data
644 def __delitem__(self, key):
645 del self._substitutions_dict[key]
647 def __contains__(self, key):
648 return key in self._substitutions_dict
650 def add_metric_data(self, metric, data):
651 """
652 Assign data to the ``metric`` tensor. The metric tensor behaves in an
653 anomalous way when raising and lowering indices.
655 Explanation
656 ===========
658 A fully covariant metric is the inverse transpose of the fully
659 contravariant metric (it is meant matrix inverse). If the metric is
660 symmetric, the transpose is not necessary and mixed
661 covariant/contravariant metrics are Kronecker deltas.
662 """
663 # hard assignment, data should not be added to `TensorHead` for metric:
664 # the problem with `TensorHead` is that the metric is anomalous, i.e.
665 # raising and lowering the index means considering the metric or its
666 # inverse, this is not the case for other tensors.
667 self._substitutions_dict_tensmul[metric, True, True] = data
668 inverse_transpose = self.inverse_transpose_matrix(data)
669 # in symmetric spaces, the transpose is the same as the original matrix,
670 # the full covariant metric tensor is the inverse transpose, so this
671 # code will be able to handle non-symmetric metrics.
672 self._substitutions_dict_tensmul[metric, False, False] = inverse_transpose
673 # now mixed cases, these are identical to the unit matrix if the metric
674 # is symmetric.
675 m = data.tomatrix()
676 invt = inverse_transpose.tomatrix()
677 self._substitutions_dict_tensmul[metric, True, False] = m * invt
678 self._substitutions_dict_tensmul[metric, False, True] = invt * m
680 @staticmethod
681 def _flip_index_by_metric(data, metric, pos):
682 from .array import tensorproduct, tensorcontraction
684 mdim = metric.rank()
685 ddim = data.rank()
687 if pos == 0:
688 data = tensorcontraction(
689 tensorproduct(
690 metric,
691 data
692 ),
693 (1, mdim+pos)
694 )
695 else:
696 data = tensorcontraction(
697 tensorproduct(
698 data,
699 metric
700 ),
701 (pos, ddim)
702 )
703 return data
705 @staticmethod
706 def inverse_matrix(ndarray):
707 m = ndarray.tomatrix().inv()
708 return _TensorDataLazyEvaluator.parse_data(m)
710 @staticmethod
711 def inverse_transpose_matrix(ndarray):
712 m = ndarray.tomatrix().inv().T
713 return _TensorDataLazyEvaluator.parse_data(m)
715 @staticmethod
716 def _correct_signature_from_indices(data, indices, free, dum, inverse=False):
717 """
718 Utility function to correct the values inside the components data
719 ndarray according to whether indices are covariant or contravariant.
721 It uses the metric matrix to lower values of covariant indices.
722 """
723 # change the ndarray values according covariantness/contravariantness of the indices
724 # use the metric
725 for i, indx in enumerate(indices):
726 if not indx.is_up and not inverse:
727 data = _TensorDataLazyEvaluator._flip_index_by_metric(data, indx.tensor_index_type.data, i)
728 elif not indx.is_up and inverse:
729 data = _TensorDataLazyEvaluator._flip_index_by_metric(
730 data,
731 _TensorDataLazyEvaluator.inverse_matrix(indx.tensor_index_type.data),
732 i
733 )
734 return data
736 @staticmethod
737 def _sort_data_axes(old, new):
738 from .array import permutedims
740 new_data = old.data.copy()
742 old_free = [i[0] for i in old.free]
743 new_free = [i[0] for i in new.free]
745 for i in range(len(new_free)):
746 for j in range(i, len(old_free)):
747 if old_free[j] == new_free[i]:
748 old_free[i], old_free[j] = old_free[j], old_free[i]
749 new_data = permutedims(new_data, (i, j))
750 break
751 return new_data
753 @staticmethod
754 def add_rearrange_tensmul_parts(new_tensmul, old_tensmul):
755 def sorted_compo():
756 return _TensorDataLazyEvaluator._sort_data_axes(old_tensmul, new_tensmul)
758 _TensorDataLazyEvaluator._substitutions_dict[new_tensmul] = sorted_compo()
760 @staticmethod
761 def parse_data(data):
762 """
763 Transform ``data`` to array. The parameter ``data`` may
764 contain data in various formats, e.g. nested lists, SymPy ``Matrix``,
765 and so on.
767 Examples
768 ========
770 >>> from sympy.tensor.tensor import _TensorDataLazyEvaluator
771 >>> _TensorDataLazyEvaluator.parse_data([1, 3, -6, 12])
772 [1, 3, -6, 12]
774 >>> _TensorDataLazyEvaluator.parse_data([[1, 2], [4, 7]])
775 [[1, 2], [4, 7]]
776 """
777 from .array import MutableDenseNDimArray
779 if not isinstance(data, MutableDenseNDimArray):
780 if len(data) == 2 and hasattr(data[0], '__call__'):
781 data = MutableDenseNDimArray(data[0], data[1])
782 else:
783 data = MutableDenseNDimArray(data)
784 return data
786_tensor_data_substitution_dict = _TensorDataLazyEvaluator()
789class _TensorManager:
790 """
791 Class to manage tensor properties.
793 Notes
794 =====
796 Tensors belong to tensor commutation groups; each group has a label
797 ``comm``; there are predefined labels:
799 ``0`` tensors commuting with any other tensor
801 ``1`` tensors anticommuting among themselves
803 ``2`` tensors not commuting, apart with those with ``comm=0``
805 Other groups can be defined using ``set_comm``; tensors in those
806 groups commute with those with ``comm=0``; by default they
807 do not commute with any other group.
808 """
809 def __init__(self):
810 self._comm_init()
812 def _comm_init(self):
813 self._comm = [{} for i in range(3)]
814 for i in range(3):
815 self._comm[0][i] = 0
816 self._comm[i][0] = 0
817 self._comm[1][1] = 1
818 self._comm[2][1] = None
819 self._comm[1][2] = None
820 self._comm_symbols2i = {0:0, 1:1, 2:2}
821 self._comm_i2symbol = {0:0, 1:1, 2:2}
823 @property
824 def comm(self):
825 return self._comm
827 def comm_symbols2i(self, i):
828 """
829 Get the commutation group number corresponding to ``i``.
831 ``i`` can be a symbol or a number or a string.
833 If ``i`` is not already defined its commutation group number
834 is set.
835 """
836 if i not in self._comm_symbols2i:
837 n = len(self._comm)
838 self._comm.append({})
839 self._comm[n][0] = 0
840 self._comm[0][n] = 0
841 self._comm_symbols2i[i] = n
842 self._comm_i2symbol[n] = i
843 return n
844 return self._comm_symbols2i[i]
846 def comm_i2symbol(self, i):
847 """
848 Returns the symbol corresponding to the commutation group number.
849 """
850 return self._comm_i2symbol[i]
852 def set_comm(self, i, j, c):
853 """
854 Set the commutation parameter ``c`` for commutation groups ``i, j``.
856 Parameters
857 ==========
859 i, j : symbols representing commutation groups
861 c : group commutation number
863 Notes
864 =====
866 ``i, j`` can be symbols, strings or numbers,
867 apart from ``0, 1`` and ``2`` which are reserved respectively
868 for commuting, anticommuting tensors and tensors not commuting
869 with any other group apart with the commuting tensors.
870 For the remaining cases, use this method to set the commutation rules;
871 by default ``c=None``.
873 The group commutation number ``c`` is assigned in correspondence
874 to the group commutation symbols; it can be
876 0 commuting
878 1 anticommuting
880 None no commutation property
882 Examples
883 ========
885 ``G`` and ``GH`` do not commute with themselves and commute with
886 each other; A is commuting.
888 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorHead, TensorManager, TensorSymmetry
889 >>> Lorentz = TensorIndexType('Lorentz')
890 >>> i0,i1,i2,i3,i4 = tensor_indices('i0:5', Lorentz)
891 >>> A = TensorHead('A', [Lorentz])
892 >>> G = TensorHead('G', [Lorentz], TensorSymmetry.no_symmetry(1), 'Gcomm')
893 >>> GH = TensorHead('GH', [Lorentz], TensorSymmetry.no_symmetry(1), 'GHcomm')
894 >>> TensorManager.set_comm('Gcomm', 'GHcomm', 0)
895 >>> (GH(i1)*G(i0)).canon_bp()
896 G(i0)*GH(i1)
897 >>> (G(i1)*G(i0)).canon_bp()
898 G(i1)*G(i0)
899 >>> (G(i1)*A(i0)).canon_bp()
900 A(i0)*G(i1)
901 """
902 if c not in (0, 1, None):
903 raise ValueError('`c` can assume only the values 0, 1 or None')
905 if i not in self._comm_symbols2i:
906 n = len(self._comm)
907 self._comm.append({})
908 self._comm[n][0] = 0
909 self._comm[0][n] = 0
910 self._comm_symbols2i[i] = n
911 self._comm_i2symbol[n] = i
912 if j not in self._comm_symbols2i:
913 n = len(self._comm)
914 self._comm.append({})
915 self._comm[0][n] = 0
916 self._comm[n][0] = 0
917 self._comm_symbols2i[j] = n
918 self._comm_i2symbol[n] = j
919 ni = self._comm_symbols2i[i]
920 nj = self._comm_symbols2i[j]
921 self._comm[ni][nj] = c
922 self._comm[nj][ni] = c
924 def set_comms(self, *args):
925 """
926 Set the commutation group numbers ``c`` for symbols ``i, j``.
928 Parameters
929 ==========
931 args : sequence of ``(i, j, c)``
932 """
933 for i, j, c in args:
934 self.set_comm(i, j, c)
936 def get_comm(self, i, j):
937 """
938 Return the commutation parameter for commutation group numbers ``i, j``
940 see ``_TensorManager.set_comm``
941 """
942 return self._comm[i].get(j, 0 if i == 0 or j == 0 else None)
944 def clear(self):
945 """
946 Clear the TensorManager.
947 """
948 self._comm_init()
951TensorManager = _TensorManager()
954class TensorIndexType(Basic):
955 """
956 A TensorIndexType is characterized by its name and its metric.
958 Parameters
959 ==========
961 name : name of the tensor type
962 dummy_name : name of the head of dummy indices
963 dim : dimension, it can be a symbol or an integer or ``None``
964 eps_dim : dimension of the epsilon tensor
965 metric_symmetry : integer that denotes metric symmetry or ``None`` for no metric
966 metric_name : string with the name of the metric tensor
968 Attributes
969 ==========
971 ``metric`` : the metric tensor
972 ``delta`` : ``Kronecker delta``
973 ``epsilon`` : the ``Levi-Civita epsilon`` tensor
974 ``data`` : (deprecated) a property to add ``ndarray`` values, to work in a specified basis.
976 Notes
977 =====
979 The possible values of the ``metric_symmetry`` parameter are:
981 ``1`` : metric tensor is fully symmetric
982 ``0`` : metric tensor possesses no index symmetry
983 ``-1`` : metric tensor is fully antisymmetric
984 ``None``: there is no metric tensor (metric equals to ``None``)
986 The metric is assumed to be symmetric by default. It can also be set
987 to a custom tensor by the ``.set_metric()`` method.
989 If there is a metric the metric is used to raise and lower indices.
991 In the case of non-symmetric metric, the following raising and
992 lowering conventions will be adopted:
994 ``psi(a) = g(a, b)*psi(-b); chi(-a) = chi(b)*g(-b, -a)``
996 From these it is easy to find:
998 ``g(-a, b) = delta(-a, b)``
1000 where ``delta(-a, b) = delta(b, -a)`` is the ``Kronecker delta``
1001 (see ``TensorIndex`` for the conventions on indices).
1002 For antisymmetric metrics there is also the following equality:
1004 ``g(a, -b) = -delta(a, -b)``
1006 If there is no metric it is not possible to raise or lower indices;
1007 e.g. the index of the defining representation of ``SU(N)``
1008 is 'covariant' and the conjugate representation is
1009 'contravariant'; for ``N > 2`` they are linearly independent.
1011 ``eps_dim`` is by default equal to ``dim``, if the latter is an integer;
1012 else it can be assigned (for use in naive dimensional regularization);
1013 if ``eps_dim`` is not an integer ``epsilon`` is ``None``.
1015 Examples
1016 ========
1018 >>> from sympy.tensor.tensor import TensorIndexType
1019 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1020 >>> Lorentz.metric
1021 metric(Lorentz,Lorentz)
1022 """
1024 def __new__(cls, name, dummy_name=None, dim=None, eps_dim=None,
1025 metric_symmetry=1, metric_name='metric', **kwargs):
1026 if 'dummy_fmt' in kwargs:
1027 dummy_fmt = kwargs['dummy_fmt']
1028 sympy_deprecation_warning(
1029 f"""
1030 The dummy_fmt keyword to TensorIndexType is deprecated. Use
1031 dummy_name={dummy_fmt} instead.
1032 """,
1033 deprecated_since_version="1.5",
1034 active_deprecations_target="deprecated-tensorindextype-dummy-fmt",
1035 )
1036 dummy_name = dummy_fmt
1038 if isinstance(name, str):
1039 name = Symbol(name)
1041 if dummy_name is None:
1042 dummy_name = str(name)[0]
1043 if isinstance(dummy_name, str):
1044 dummy_name = Symbol(dummy_name)
1046 if dim is None:
1047 dim = Symbol("dim_" + dummy_name.name)
1048 else:
1049 dim = sympify(dim)
1051 if eps_dim is None:
1052 eps_dim = dim
1053 else:
1054 eps_dim = sympify(eps_dim)
1056 metric_symmetry = sympify(metric_symmetry)
1058 if isinstance(metric_name, str):
1059 metric_name = Symbol(metric_name)
1061 if 'metric' in kwargs:
1062 SymPyDeprecationWarning(
1063 """
1064 The 'metric' keyword argument to TensorIndexType is
1065 deprecated. Use the 'metric_symmetry' keyword argument or the
1066 TensorIndexType.set_metric() method instead.
1067 """,
1068 deprecated_since_version="1.5",
1069 active_deprecations_target="deprecated-tensorindextype-metric",
1070 )
1071 metric = kwargs.get('metric')
1072 if metric is not None:
1073 if metric in (True, False, 0, 1):
1074 metric_name = 'metric'
1075 #metric_antisym = metric
1076 else:
1077 metric_name = metric.name
1078 #metric_antisym = metric.antisym
1080 if metric:
1081 metric_symmetry = -1
1082 else:
1083 metric_symmetry = 1
1085 obj = Basic.__new__(cls, name, dummy_name, dim, eps_dim,
1086 metric_symmetry, metric_name)
1088 obj._autogenerated = []
1089 return obj
1091 @property
1092 def name(self):
1093 return self.args[0].name
1095 @property
1096 def dummy_name(self):
1097 return self.args[1].name
1099 @property
1100 def dim(self):
1101 return self.args[2]
1103 @property
1104 def eps_dim(self):
1105 return self.args[3]
1107 @memoize_property
1108 def metric(self):
1109 metric_symmetry = self.args[4]
1110 metric_name = self.args[5]
1111 if metric_symmetry is None:
1112 return None
1114 if metric_symmetry == 0:
1115 symmetry = TensorSymmetry.no_symmetry(2)
1116 elif metric_symmetry == 1:
1117 symmetry = TensorSymmetry.fully_symmetric(2)
1118 elif metric_symmetry == -1:
1119 symmetry = TensorSymmetry.fully_symmetric(-2)
1121 return TensorHead(metric_name, [self]*2, symmetry)
1123 @memoize_property
1124 def delta(self):
1125 return TensorHead('KD', [self]*2, TensorSymmetry.fully_symmetric(2))
1127 @memoize_property
1128 def epsilon(self):
1129 if not isinstance(self.eps_dim, (SYMPY_INTS, Integer)):
1130 return None
1131 symmetry = TensorSymmetry.fully_symmetric(-self.eps_dim)
1132 return TensorHead('Eps', [self]*self.eps_dim, symmetry)
1134 def set_metric(self, tensor):
1135 self._metric = tensor
1137 def __lt__(self, other):
1138 return self.name < other.name
1140 def __str__(self):
1141 return self.name
1143 __repr__ = __str__
1145 # Everything below this line is deprecated
1147 @property
1148 def data(self):
1149 deprecate_data()
1150 with ignore_warnings(SymPyDeprecationWarning):
1151 return _tensor_data_substitution_dict[self]
1153 @data.setter
1154 def data(self, data):
1155 deprecate_data()
1156 # This assignment is a bit controversial, should metric components be assigned
1157 # to the metric only or also to the TensorIndexType object? The advantage here
1158 # is the ability to assign a 1D array and transform it to a 2D diagonal array.
1159 from .array import MutableDenseNDimArray
1161 data = _TensorDataLazyEvaluator.parse_data(data)
1162 if data.rank() > 2:
1163 raise ValueError("data have to be of rank 1 (diagonal metric) or 2.")
1164 if data.rank() == 1:
1165 if self.dim.is_number:
1166 nda_dim = data.shape[0]
1167 if nda_dim != self.dim:
1168 raise ValueError("Dimension mismatch")
1170 dim = data.shape[0]
1171 newndarray = MutableDenseNDimArray.zeros(dim, dim)
1172 for i, val in enumerate(data):
1173 newndarray[i, i] = val
1174 data = newndarray
1175 dim1, dim2 = data.shape
1176 if dim1 != dim2:
1177 raise ValueError("Non-square matrix tensor.")
1178 if self.dim.is_number:
1179 if self.dim != dim1:
1180 raise ValueError("Dimension mismatch")
1181 _tensor_data_substitution_dict[self] = data
1182 _tensor_data_substitution_dict.add_metric_data(self.metric, data)
1183 with ignore_warnings(SymPyDeprecationWarning):
1184 delta = self.get_kronecker_delta()
1185 i1 = TensorIndex('i1', self)
1186 i2 = TensorIndex('i2', self)
1187 with ignore_warnings(SymPyDeprecationWarning):
1188 delta(i1, -i2).data = _TensorDataLazyEvaluator.parse_data(eye(dim1))
1190 @data.deleter
1191 def data(self):
1192 deprecate_data()
1193 with ignore_warnings(SymPyDeprecationWarning):
1194 if self in _tensor_data_substitution_dict:
1195 del _tensor_data_substitution_dict[self]
1196 if self.metric in _tensor_data_substitution_dict:
1197 del _tensor_data_substitution_dict[self.metric]
1199 @deprecated(
1200 """
1201 The TensorIndexType.get_kronecker_delta() method is deprecated. Use
1202 the TensorIndexType.delta attribute instead.
1203 """,
1204 deprecated_since_version="1.5",
1205 active_deprecations_target="deprecated-tensorindextype-methods",
1206 )
1207 def get_kronecker_delta(self):
1208 sym2 = TensorSymmetry(get_symmetric_group_sgs(2))
1209 delta = TensorHead('KD', [self]*2, sym2)
1210 return delta
1212 @deprecated(
1213 """
1214 The TensorIndexType.get_epsilon() method is deprecated. Use
1215 the TensorIndexType.epsilon attribute instead.
1216 """,
1217 deprecated_since_version="1.5",
1218 active_deprecations_target="deprecated-tensorindextype-methods",
1219 )
1220 def get_epsilon(self):
1221 if not isinstance(self._eps_dim, (SYMPY_INTS, Integer)):
1222 return None
1223 sym = TensorSymmetry(get_symmetric_group_sgs(self._eps_dim, 1))
1224 epsilon = TensorHead('Eps', [self]*self._eps_dim, sym)
1225 return epsilon
1227 def _components_data_full_destroy(self):
1228 """
1229 EXPERIMENTAL: do not rely on this API method.
1231 This destroys components data associated to the ``TensorIndexType``, if
1232 any, specifically:
1234 * metric tensor data
1235 * Kronecker tensor data
1236 """
1237 if self in _tensor_data_substitution_dict:
1238 del _tensor_data_substitution_dict[self]
1240 def delete_tensmul_data(key):
1241 if key in _tensor_data_substitution_dict._substitutions_dict_tensmul:
1242 del _tensor_data_substitution_dict._substitutions_dict_tensmul[key]
1244 # delete metric data:
1245 delete_tensmul_data((self.metric, True, True))
1246 delete_tensmul_data((self.metric, True, False))
1247 delete_tensmul_data((self.metric, False, True))
1248 delete_tensmul_data((self.metric, False, False))
1250 # delete delta tensor data:
1251 delta = self.get_kronecker_delta()
1252 if delta in _tensor_data_substitution_dict:
1253 del _tensor_data_substitution_dict[delta]
1256class TensorIndex(Basic):
1257 """
1258 Represents a tensor index
1260 Parameters
1261 ==========
1263 name : name of the index, or ``True`` if you want it to be automatically assigned
1264 tensor_index_type : ``TensorIndexType`` of the index
1265 is_up : flag for contravariant index (is_up=True by default)
1267 Attributes
1268 ==========
1270 ``name``
1271 ``tensor_index_type``
1272 ``is_up``
1274 Notes
1275 =====
1277 Tensor indices are contracted with the Einstein summation convention.
1279 An index can be in contravariant or in covariant form; in the latter
1280 case it is represented prepending a ``-`` to the index name. Adding
1281 ``-`` to a covariant (is_up=False) index makes it contravariant.
1283 Dummy indices have a name with head given by
1284 ``tensor_inde_type.dummy_name`` with underscore and a number.
1286 Similar to ``symbols`` multiple contravariant indices can be created
1287 at once using ``tensor_indices(s, typ)``, where ``s`` is a string
1288 of names.
1291 Examples
1292 ========
1294 >>> from sympy.tensor.tensor import TensorIndexType, TensorIndex, TensorHead, tensor_indices
1295 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1296 >>> mu = TensorIndex('mu', Lorentz, is_up=False)
1297 >>> nu, rho = tensor_indices('nu, rho', Lorentz)
1298 >>> A = TensorHead('A', [Lorentz, Lorentz])
1299 >>> A(mu, nu)
1300 A(-mu, nu)
1301 >>> A(-mu, -rho)
1302 A(mu, -rho)
1303 >>> A(mu, -mu)
1304 A(-L_0, L_0)
1305 """
1306 def __new__(cls, name, tensor_index_type, is_up=True):
1307 if isinstance(name, str):
1308 name_symbol = Symbol(name)
1309 elif isinstance(name, Symbol):
1310 name_symbol = name
1311 elif name is True:
1312 name = "_i{}".format(len(tensor_index_type._autogenerated))
1313 name_symbol = Symbol(name)
1314 tensor_index_type._autogenerated.append(name_symbol)
1315 else:
1316 raise ValueError("invalid name")
1318 is_up = sympify(is_up)
1319 return Basic.__new__(cls, name_symbol, tensor_index_type, is_up)
1321 @property
1322 def name(self):
1323 return self.args[0].name
1325 @property
1326 def tensor_index_type(self):
1327 return self.args[1]
1329 @property
1330 def is_up(self):
1331 return self.args[2]
1333 def _print(self):
1334 s = self.name
1335 if not self.is_up:
1336 s = '-%s' % s
1337 return s
1339 def __lt__(self, other):
1340 return ((self.tensor_index_type, self.name) <
1341 (other.tensor_index_type, other.name))
1343 def __neg__(self):
1344 t1 = TensorIndex(self.name, self.tensor_index_type,
1345 (not self.is_up))
1346 return t1
1349def tensor_indices(s, typ):
1350 """
1351 Returns list of tensor indices given their names and their types.
1353 Parameters
1354 ==========
1356 s : string of comma separated names of indices
1358 typ : ``TensorIndexType`` of the indices
1360 Examples
1361 ========
1363 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices
1364 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1365 >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz)
1366 """
1367 if isinstance(s, str):
1368 a = [x.name for x in symbols(s, seq=True)]
1369 else:
1370 raise ValueError('expecting a string')
1372 tilist = [TensorIndex(i, typ) for i in a]
1373 if len(tilist) == 1:
1374 return tilist[0]
1375 return tilist
1378class TensorSymmetry(Basic):
1379 """
1380 Monoterm symmetry of a tensor (i.e. any symmetric or anti-symmetric
1381 index permutation). For the relevant terminology see ``tensor_can.py``
1382 section of the combinatorics module.
1384 Parameters
1385 ==========
1387 bsgs : tuple ``(base, sgs)`` BSGS of the symmetry of the tensor
1389 Attributes
1390 ==========
1392 ``base`` : base of the BSGS
1393 ``generators`` : generators of the BSGS
1394 ``rank`` : rank of the tensor
1396 Notes
1397 =====
1399 A tensor can have an arbitrary monoterm symmetry provided by its BSGS.
1400 Multiterm symmetries, like the cyclic symmetry of the Riemann tensor
1401 (i.e., Bianchi identity), are not covered. See combinatorics module for
1402 information on how to generate BSGS for a general index permutation group.
1403 Simple symmetries can be generated using built-in methods.
1405 See Also
1406 ========
1408 sympy.combinatorics.tensor_can.get_symmetric_group_sgs
1410 Examples
1411 ========
1413 Define a symmetric tensor of rank 2
1415 >>> from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, get_symmetric_group_sgs, TensorHead
1416 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1417 >>> sym = TensorSymmetry(get_symmetric_group_sgs(2))
1418 >>> T = TensorHead('T', [Lorentz]*2, sym)
1420 Note, that the same can also be done using built-in TensorSymmetry methods
1422 >>> sym2 = TensorSymmetry.fully_symmetric(2)
1423 >>> sym == sym2
1424 True
1425 """
1426 def __new__(cls, *args, **kw_args):
1427 if len(args) == 1:
1428 base, generators = args[0]
1429 elif len(args) == 2:
1430 base, generators = args
1431 else:
1432 raise TypeError("bsgs required, either two separate parameters or one tuple")
1434 if not isinstance(base, Tuple):
1435 base = Tuple(*base)
1436 if not isinstance(generators, Tuple):
1437 generators = Tuple(*generators)
1439 return Basic.__new__(cls, base, generators, **kw_args)
1441 @property
1442 def base(self):
1443 return self.args[0]
1445 @property
1446 def generators(self):
1447 return self.args[1]
1449 @property
1450 def rank(self):
1451 return self.generators[0].size - 2
1453 @classmethod
1454 def fully_symmetric(cls, rank):
1455 """
1456 Returns a fully symmetric (antisymmetric if ``rank``<0)
1457 TensorSymmetry object for ``abs(rank)`` indices.
1458 """
1459 if rank > 0:
1460 bsgs = get_symmetric_group_sgs(rank, False)
1461 elif rank < 0:
1462 bsgs = get_symmetric_group_sgs(-rank, True)
1463 elif rank == 0:
1464 bsgs = ([], [Permutation(1)])
1465 return TensorSymmetry(bsgs)
1467 @classmethod
1468 def direct_product(cls, *args):
1469 """
1470 Returns a TensorSymmetry object that is being a direct product of
1471 fully (anti-)symmetric index permutation groups.
1473 Notes
1474 =====
1476 Some examples for different values of ``(*args)``:
1477 ``(1)`` vector, equivalent to ``TensorSymmetry.fully_symmetric(1)``
1478 ``(2)`` tensor with 2 symmetric indices, equivalent to ``.fully_symmetric(2)``
1479 ``(-2)`` tensor with 2 antisymmetric indices, equivalent to ``.fully_symmetric(-2)``
1480 ``(2, -2)`` tensor with the first 2 indices commuting and the last 2 anticommuting
1481 ``(1, 1, 1)`` tensor with 3 indices without any symmetry
1482 """
1483 base, sgs = [], [Permutation(1)]
1484 for arg in args:
1485 if arg > 0:
1486 bsgs2 = get_symmetric_group_sgs(arg, False)
1487 elif arg < 0:
1488 bsgs2 = get_symmetric_group_sgs(-arg, True)
1489 else:
1490 continue
1491 base, sgs = bsgs_direct_product(base, sgs, *bsgs2)
1493 return TensorSymmetry(base, sgs)
1495 @classmethod
1496 def riemann(cls):
1497 """
1498 Returns a monotorem symmetry of the Riemann tensor
1499 """
1500 return TensorSymmetry(riemann_bsgs)
1502 @classmethod
1503 def no_symmetry(cls, rank):
1504 """
1505 TensorSymmetry object for ``rank`` indices with no symmetry
1506 """
1507 return TensorSymmetry([], [Permutation(rank+1)])
1510@deprecated(
1511 """
1512 The tensorsymmetry() function is deprecated. Use the TensorSymmetry
1513 constructor instead.
1514 """,
1515 deprecated_since_version="1.5",
1516 active_deprecations_target="deprecated-tensorsymmetry",
1517)
1518def tensorsymmetry(*args):
1519 """
1520 Returns a ``TensorSymmetry`` object. This method is deprecated, use
1521 ``TensorSymmetry.direct_product()`` or ``.riemann()`` instead.
1523 Explanation
1524 ===========
1526 One can represent a tensor with any monoterm slot symmetry group
1527 using a BSGS.
1529 ``args`` can be a BSGS
1530 ``args[0]`` base
1531 ``args[1]`` sgs
1533 Usually tensors are in (direct products of) representations
1534 of the symmetric group;
1535 ``args`` can be a list of lists representing the shapes of Young tableaux
1537 Notes
1538 =====
1540 For instance:
1541 ``[[1]]`` vector
1542 ``[[1]*n]`` symmetric tensor of rank ``n``
1543 ``[[n]]`` antisymmetric tensor of rank ``n``
1544 ``[[2, 2]]`` monoterm slot symmetry of the Riemann tensor
1545 ``[[1],[1]]`` vector*vector
1546 ``[[2],[1],[1]`` (antisymmetric tensor)*vector*vector
1548 Notice that with the shape ``[2, 2]`` we associate only the monoterm
1549 symmetries of the Riemann tensor; this is an abuse of notation,
1550 since the shape ``[2, 2]`` corresponds usually to the irreducible
1551 representation characterized by the monoterm symmetries and by the
1552 cyclic symmetry.
1553 """
1554 from sympy.combinatorics import Permutation
1556 def tableau2bsgs(a):
1557 if len(a) == 1:
1558 # antisymmetric vector
1559 n = a[0]
1560 bsgs = get_symmetric_group_sgs(n, 1)
1561 else:
1562 if all(x == 1 for x in a):
1563 # symmetric vector
1564 n = len(a)
1565 bsgs = get_symmetric_group_sgs(n)
1566 elif a == [2, 2]:
1567 bsgs = riemann_bsgs
1568 else:
1569 raise NotImplementedError
1570 return bsgs
1572 if not args:
1573 return TensorSymmetry(Tuple(), Tuple(Permutation(1)))
1575 if len(args) == 2 and isinstance(args[1][0], Permutation):
1576 return TensorSymmetry(args)
1577 base, sgs = tableau2bsgs(args[0])
1578 for a in args[1:]:
1579 basex, sgsx = tableau2bsgs(a)
1580 base, sgs = bsgs_direct_product(base, sgs, basex, sgsx)
1581 return TensorSymmetry(Tuple(base, sgs))
1583@deprecated(
1584 "TensorType is deprecated. Use tensor_heads() instead.",
1585 deprecated_since_version="1.5",
1586 active_deprecations_target="deprecated-tensortype",
1587)
1588class TensorType(Basic):
1589 """
1590 Class of tensor types. Deprecated, use tensor_heads() instead.
1592 Parameters
1593 ==========
1595 index_types : list of ``TensorIndexType`` of the tensor indices
1596 symmetry : ``TensorSymmetry`` of the tensor
1598 Attributes
1599 ==========
1601 ``index_types``
1602 ``symmetry``
1603 ``types`` : list of ``TensorIndexType`` without repetitions
1604 """
1605 is_commutative = False
1607 def __new__(cls, index_types, symmetry, **kw_args):
1608 assert symmetry.rank == len(index_types)
1609 obj = Basic.__new__(cls, Tuple(*index_types), symmetry, **kw_args)
1610 return obj
1612 @property
1613 def index_types(self):
1614 return self.args[0]
1616 @property
1617 def symmetry(self):
1618 return self.args[1]
1620 @property
1621 def types(self):
1622 return sorted(set(self.index_types), key=lambda x: x.name)
1624 def __str__(self):
1625 return 'TensorType(%s)' % ([str(x) for x in self.index_types])
1627 def __call__(self, s, comm=0):
1628 """
1629 Return a TensorHead object or a list of TensorHead objects.
1631 Parameters
1632 ==========
1634 s : name or string of names.
1636 comm : Commutation group.
1638 see ``_TensorManager.set_comm``
1639 """
1640 if isinstance(s, str):
1641 names = [x.name for x in symbols(s, seq=True)]
1642 else:
1643 raise ValueError('expecting a string')
1644 if len(names) == 1:
1645 return TensorHead(names[0], self.index_types, self.symmetry, comm)
1646 else:
1647 return [TensorHead(name, self.index_types, self.symmetry, comm) for name in names]
1650@deprecated(
1651 """
1652 The tensorhead() function is deprecated. Use tensor_heads() instead.
1653 """,
1654 deprecated_since_version="1.5",
1655 active_deprecations_target="deprecated-tensorhead",
1656)
1657def tensorhead(name, typ, sym=None, comm=0):
1658 """
1659 Function generating tensorhead(s). This method is deprecated,
1660 use TensorHead constructor or tensor_heads() instead.
1662 Parameters
1663 ==========
1665 name : name or sequence of names (as in ``symbols``)
1667 typ : index types
1669 sym : same as ``*args`` in ``tensorsymmetry``
1671 comm : commutation group number
1672 see ``_TensorManager.set_comm``
1673 """
1674 if sym is None:
1675 sym = [[1] for i in range(len(typ))]
1676 with ignore_warnings(SymPyDeprecationWarning):
1677 sym = tensorsymmetry(*sym)
1678 return TensorHead(name, typ, sym, comm)
1681class TensorHead(Basic):
1682 """
1683 Tensor head of the tensor.
1685 Parameters
1686 ==========
1688 name : name of the tensor
1689 index_types : list of TensorIndexType
1690 symmetry : TensorSymmetry of the tensor
1691 comm : commutation group number
1693 Attributes
1694 ==========
1696 ``name``
1697 ``index_types``
1698 ``rank`` : total number of indices
1699 ``symmetry``
1700 ``comm`` : commutation group
1702 Notes
1703 =====
1705 Similar to ``symbols`` multiple TensorHeads can be created using
1706 ``tensorhead(s, typ, sym=None, comm=0)`` function, where ``s``
1707 is the string of names and ``sym`` is the monoterm tensor symmetry
1708 (see ``tensorsymmetry``).
1710 A ``TensorHead`` belongs to a commutation group, defined by a
1711 symbol on number ``comm`` (see ``_TensorManager.set_comm``);
1712 tensors in a commutation group have the same commutation properties;
1713 by default ``comm`` is ``0``, the group of the commuting tensors.
1715 Examples
1716 ========
1718 Define a fully antisymmetric tensor of rank 2:
1720 >>> from sympy.tensor.tensor import TensorIndexType, TensorHead, TensorSymmetry
1721 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1722 >>> asym2 = TensorSymmetry.fully_symmetric(-2)
1723 >>> A = TensorHead('A', [Lorentz, Lorentz], asym2)
1725 Examples with ndarray values, the components data assigned to the
1726 ``TensorHead`` object are assumed to be in a fully-contravariant
1727 representation. In case it is necessary to assign components data which
1728 represents the values of a non-fully covariant tensor, see the other
1729 examples.
1731 >>> from sympy.tensor.tensor import tensor_indices
1732 >>> from sympy import diag
1733 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1734 >>> i0, i1 = tensor_indices('i0:2', Lorentz)
1736 Specify a replacement dictionary to keep track of the arrays to use for
1737 replacements in the tensorial expression. The ``TensorIndexType`` is
1738 associated to the metric used for contractions (in fully covariant form):
1740 >>> repl = {Lorentz: diag(1, -1, -1, -1)}
1742 Let's see some examples of working with components with the electromagnetic
1743 tensor:
1745 >>> from sympy import symbols
1746 >>> Ex, Ey, Ez, Bx, By, Bz = symbols('E_x E_y E_z B_x B_y B_z')
1747 >>> c = symbols('c', positive=True)
1749 Let's define `F`, an antisymmetric tensor:
1751 >>> F = TensorHead('F', [Lorentz, Lorentz], asym2)
1753 Let's update the dictionary to contain the matrix to use in the
1754 replacements:
1756 >>> repl.update({F(-i0, -i1): [
1757 ... [0, Ex/c, Ey/c, Ez/c],
1758 ... [-Ex/c, 0, -Bz, By],
1759 ... [-Ey/c, Bz, 0, -Bx],
1760 ... [-Ez/c, -By, Bx, 0]]})
1762 Now it is possible to retrieve the contravariant form of the Electromagnetic
1763 tensor:
1765 >>> F(i0, i1).replace_with_arrays(repl, [i0, i1])
1766 [[0, -E_x/c, -E_y/c, -E_z/c], [E_x/c, 0, -B_z, B_y], [E_y/c, B_z, 0, -B_x], [E_z/c, -B_y, B_x, 0]]
1768 and the mixed contravariant-covariant form:
1770 >>> F(i0, -i1).replace_with_arrays(repl, [i0, -i1])
1771 [[0, E_x/c, E_y/c, E_z/c], [E_x/c, 0, B_z, -B_y], [E_y/c, -B_z, 0, B_x], [E_z/c, B_y, -B_x, 0]]
1773 Energy-momentum of a particle may be represented as:
1775 >>> from sympy import symbols
1776 >>> P = TensorHead('P', [Lorentz], TensorSymmetry.no_symmetry(1))
1777 >>> E, px, py, pz = symbols('E p_x p_y p_z', positive=True)
1778 >>> repl.update({P(i0): [E, px, py, pz]})
1780 The contravariant and covariant components are, respectively:
1782 >>> P(i0).replace_with_arrays(repl, [i0])
1783 [E, p_x, p_y, p_z]
1784 >>> P(-i0).replace_with_arrays(repl, [-i0])
1785 [E, -p_x, -p_y, -p_z]
1787 The contraction of a 1-index tensor by itself:
1789 >>> expr = P(i0)*P(-i0)
1790 >>> expr.replace_with_arrays(repl, [])
1791 E**2 - p_x**2 - p_y**2 - p_z**2
1792 """
1793 is_commutative = False
1795 def __new__(cls, name, index_types, symmetry=None, comm=0):
1796 if isinstance(name, str):
1797 name_symbol = Symbol(name)
1798 elif isinstance(name, Symbol):
1799 name_symbol = name
1800 else:
1801 raise ValueError("invalid name")
1803 if symmetry is None:
1804 symmetry = TensorSymmetry.no_symmetry(len(index_types))
1805 else:
1806 assert symmetry.rank == len(index_types)
1808 obj = Basic.__new__(cls, name_symbol, Tuple(*index_types), symmetry)
1809 obj.comm = TensorManager.comm_symbols2i(comm)
1810 return obj
1812 @property
1813 def name(self):
1814 return self.args[0].name
1816 @property
1817 def index_types(self):
1818 return list(self.args[1])
1820 @property
1821 def symmetry(self):
1822 return self.args[2]
1824 @property
1825 def rank(self):
1826 return len(self.index_types)
1828 def __lt__(self, other):
1829 return (self.name, self.index_types) < (other.name, other.index_types)
1831 def commutes_with(self, other):
1832 """
1833 Returns ``0`` if ``self`` and ``other`` commute, ``1`` if they anticommute.
1835 Returns ``None`` if ``self`` and ``other`` neither commute nor anticommute.
1836 """
1837 r = TensorManager.get_comm(self.comm, other.comm)
1838 return r
1840 def _print(self):
1841 return '%s(%s)' %(self.name, ','.join([str(x) for x in self.index_types]))
1843 def __call__(self, *indices, **kw_args):
1844 """
1845 Returns a tensor with indices.
1847 Explanation
1848 ===========
1850 There is a special behavior in case of indices denoted by ``True``,
1851 they are considered auto-matrix indices, their slots are automatically
1852 filled, and confer to the tensor the behavior of a matrix or vector
1853 upon multiplication with another tensor containing auto-matrix indices
1854 of the same ``TensorIndexType``. This means indices get summed over the
1855 same way as in matrix multiplication. For matrix behavior, define two
1856 auto-matrix indices, for vector behavior define just one.
1858 Indices can also be strings, in which case the attribute
1859 ``index_types`` is used to convert them to proper ``TensorIndex``.
1861 Examples
1862 ========
1864 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorSymmetry, TensorHead
1865 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
1866 >>> a, b = tensor_indices('a,b', Lorentz)
1867 >>> A = TensorHead('A', [Lorentz]*2, TensorSymmetry.no_symmetry(2))
1868 >>> t = A(a, -b)
1869 >>> t
1870 A(a, -b)
1872 """
1874 updated_indices = []
1875 for idx, typ in zip(indices, self.index_types):
1876 if isinstance(idx, str):
1877 idx = idx.strip().replace(" ", "")
1878 if idx.startswith('-'):
1879 updated_indices.append(TensorIndex(idx[1:], typ,
1880 is_up=False))
1881 else:
1882 updated_indices.append(TensorIndex(idx, typ))
1883 else:
1884 updated_indices.append(idx)
1886 updated_indices += indices[len(updated_indices):]
1888 tensor = Tensor(self, updated_indices, **kw_args)
1889 return tensor.doit()
1891 # Everything below this line is deprecated
1893 def __pow__(self, other):
1894 deprecate_data()
1895 with ignore_warnings(SymPyDeprecationWarning):
1896 if self.data is None:
1897 raise ValueError("No power on abstract tensors.")
1898 from .array import tensorproduct, tensorcontraction
1899 metrics = [_.data for _ in self.index_types]
1901 marray = self.data
1902 marraydim = marray.rank()
1903 for metric in metrics:
1904 marray = tensorproduct(marray, metric, marray)
1905 marray = tensorcontraction(marray, (0, marraydim), (marraydim+1, marraydim+2))
1907 return marray ** (other * S.Half)
1909 @property
1910 def data(self):
1911 deprecate_data()
1912 with ignore_warnings(SymPyDeprecationWarning):
1913 return _tensor_data_substitution_dict[self]
1915 @data.setter
1916 def data(self, data):
1917 deprecate_data()
1918 with ignore_warnings(SymPyDeprecationWarning):
1919 _tensor_data_substitution_dict[self] = data
1921 @data.deleter
1922 def data(self):
1923 deprecate_data()
1924 if self in _tensor_data_substitution_dict:
1925 del _tensor_data_substitution_dict[self]
1927 def __iter__(self):
1928 deprecate_data()
1929 with ignore_warnings(SymPyDeprecationWarning):
1930 return self.data.__iter__()
1932 def _components_data_full_destroy(self):
1933 """
1934 EXPERIMENTAL: do not rely on this API method.
1936 Destroy components data associated to the ``TensorHead`` object, this
1937 checks for attached components data, and destroys components data too.
1938 """
1939 # do not garbage collect Kronecker tensor (it should be done by
1940 # ``TensorIndexType`` garbage collection)
1941 deprecate_data()
1942 if self.name == "KD":
1943 return
1945 # the data attached to a tensor must be deleted only by the TensorHead
1946 # destructor. If the TensorHead is deleted, it means that there are no
1947 # more instances of that tensor anywhere.
1948 if self in _tensor_data_substitution_dict:
1949 del _tensor_data_substitution_dict[self]
1952def tensor_heads(s, index_types, symmetry=None, comm=0):
1953 """
1954 Returns a sequence of TensorHeads from a string `s`
1955 """
1956 if isinstance(s, str):
1957 names = [x.name for x in symbols(s, seq=True)]
1958 else:
1959 raise ValueError('expecting a string')
1961 thlist = [TensorHead(name, index_types, symmetry, comm) for name in names]
1962 if len(thlist) == 1:
1963 return thlist[0]
1964 return thlist
1967class TensExpr(Expr, ABC):
1968 """
1969 Abstract base class for tensor expressions
1971 Notes
1972 =====
1974 A tensor expression is an expression formed by tensors;
1975 currently the sums of tensors are distributed.
1977 A ``TensExpr`` can be a ``TensAdd`` or a ``TensMul``.
1979 ``TensMul`` objects are formed by products of component tensors,
1980 and include a coefficient, which is a SymPy expression.
1983 In the internal representation contracted indices are represented
1984 by ``(ipos1, ipos2, icomp1, icomp2)``, where ``icomp1`` is the position
1985 of the component tensor with contravariant index, ``ipos1`` is the
1986 slot which the index occupies in that component tensor.
1988 Contracted indices are therefore nameless in the internal representation.
1989 """
1991 _op_priority = 12.0
1992 is_commutative = False
1994 def __neg__(self):
1995 return self*S.NegativeOne
1997 def __abs__(self):
1998 raise NotImplementedError
2000 def __add__(self, other):
2001 return TensAdd(self, other).doit()
2003 def __radd__(self, other):
2004 return TensAdd(other, self).doit()
2006 def __sub__(self, other):
2007 return TensAdd(self, -other).doit()
2009 def __rsub__(self, other):
2010 return TensAdd(other, -self).doit()
2012 def __mul__(self, other):
2013 """
2014 Multiply two tensors using Einstein summation convention.
2016 Explanation
2017 ===========
2019 If the two tensors have an index in common, one contravariant
2020 and the other covariant, in their product the indices are summed
2022 Examples
2023 ========
2025 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensor_heads
2026 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
2027 >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz)
2028 >>> g = Lorentz.metric
2029 >>> p, q = tensor_heads('p,q', [Lorentz])
2030 >>> t1 = p(m0)
2031 >>> t2 = q(-m0)
2032 >>> t1*t2
2033 p(L_0)*q(-L_0)
2034 """
2035 return TensMul(self, other).doit()
2037 def __rmul__(self, other):
2038 return TensMul(other, self).doit()
2040 def __truediv__(self, other):
2041 other = _sympify(other)
2042 if isinstance(other, TensExpr):
2043 raise ValueError('cannot divide by a tensor')
2044 return TensMul(self, S.One/other).doit()
2046 def __rtruediv__(self, other):
2047 raise ValueError('cannot divide by a tensor')
2049 def __pow__(self, other):
2050 deprecate_data()
2051 with ignore_warnings(SymPyDeprecationWarning):
2052 if self.data is None:
2053 raise ValueError("No power without ndarray data.")
2054 from .array import tensorproduct, tensorcontraction
2055 free = self.free
2056 marray = self.data
2057 mdim = marray.rank()
2058 for metric in free:
2059 marray = tensorcontraction(
2060 tensorproduct(
2061 marray,
2062 metric[0].tensor_index_type.data,
2063 marray),
2064 (0, mdim), (mdim+1, mdim+2)
2065 )
2066 return marray ** (other * S.Half)
2068 def __rpow__(self, other):
2069 raise NotImplementedError
2071 @property
2072 @abstractmethod
2073 def nocoeff(self):
2074 raise NotImplementedError("abstract method")
2076 @property
2077 @abstractmethod
2078 def coeff(self):
2079 raise NotImplementedError("abstract method")
2081 @abstractmethod
2082 def get_indices(self):
2083 raise NotImplementedError("abstract method")
2085 @abstractmethod
2086 def get_free_indices(self) -> list[TensorIndex]:
2087 raise NotImplementedError("abstract method")
2089 @abstractmethod
2090 def _replace_indices(self, repl: dict[TensorIndex, TensorIndex]) -> TensExpr:
2091 raise NotImplementedError("abstract method")
2093 def fun_eval(self, *index_tuples):
2094 deprecate_fun_eval()
2095 return self.substitute_indices(*index_tuples)
2097 def get_matrix(self):
2098 """
2099 DEPRECATED: do not use.
2101 Returns ndarray components data as a matrix, if components data are
2102 available and ndarray dimension does not exceed 2.
2103 """
2104 from sympy.matrices.dense import Matrix
2105 deprecate_data()
2106 with ignore_warnings(SymPyDeprecationWarning):
2107 if 0 < self.rank <= 2:
2108 rows = self.data.shape[0]
2109 columns = self.data.shape[1] if self.rank == 2 else 1
2110 if self.rank == 2:
2111 mat_list = [] * rows
2112 for i in range(rows):
2113 mat_list.append([])
2114 for j in range(columns):
2115 mat_list[i].append(self[i, j])
2116 else:
2117 mat_list = [None] * rows
2118 for i in range(rows):
2119 mat_list[i] = self[i]
2120 return Matrix(mat_list)
2121 else:
2122 raise NotImplementedError(
2123 "missing multidimensional reduction to matrix.")
2125 @staticmethod
2126 def _get_indices_permutation(indices1, indices2):
2127 return [indices1.index(i) for i in indices2]
2129 def expand(self, **hints):
2130 return _expand(self, **hints).doit()
2132 def _expand(self, **kwargs):
2133 return self
2135 def _get_free_indices_set(self):
2136 indset = set()
2137 for arg in self.args:
2138 if isinstance(arg, TensExpr):
2139 indset.update(arg._get_free_indices_set())
2140 return indset
2142 def _get_dummy_indices_set(self):
2143 indset = set()
2144 for arg in self.args:
2145 if isinstance(arg, TensExpr):
2146 indset.update(arg._get_dummy_indices_set())
2147 return indset
2149 def _get_indices_set(self):
2150 indset = set()
2151 for arg in self.args:
2152 if isinstance(arg, TensExpr):
2153 indset.update(arg._get_indices_set())
2154 return indset
2156 @property
2157 def _iterate_dummy_indices(self):
2158 dummy_set = self._get_dummy_indices_set()
2160 def recursor(expr, pos):
2161 if isinstance(expr, TensorIndex):
2162 if expr in dummy_set:
2163 yield (expr, pos)
2164 elif isinstance(expr, (Tuple, TensExpr)):
2165 for p, arg in enumerate(expr.args):
2166 yield from recursor(arg, pos+(p,))
2168 return recursor(self, ())
2170 @property
2171 def _iterate_free_indices(self):
2172 free_set = self._get_free_indices_set()
2174 def recursor(expr, pos):
2175 if isinstance(expr, TensorIndex):
2176 if expr in free_set:
2177 yield (expr, pos)
2178 elif isinstance(expr, (Tuple, TensExpr)):
2179 for p, arg in enumerate(expr.args):
2180 yield from recursor(arg, pos+(p,))
2182 return recursor(self, ())
2184 @property
2185 def _iterate_indices(self):
2186 def recursor(expr, pos):
2187 if isinstance(expr, TensorIndex):
2188 yield (expr, pos)
2189 elif isinstance(expr, (Tuple, TensExpr)):
2190 for p, arg in enumerate(expr.args):
2191 yield from recursor(arg, pos+(p,))
2193 return recursor(self, ())
2195 @staticmethod
2196 def _contract_and_permute_with_metric(metric, array, pos, dim):
2197 # TODO: add possibility of metric after (spinors)
2198 from .array import tensorcontraction, tensorproduct, permutedims
2200 array = tensorcontraction(tensorproduct(metric, array), (1, 2+pos))
2201 permu = list(range(dim))
2202 permu[0], permu[pos] = permu[pos], permu[0]
2203 return permutedims(array, permu)
2205 @staticmethod
2206 def _match_indices_with_other_tensor(array, free_ind1, free_ind2, replacement_dict):
2207 from .array import permutedims
2209 index_types1 = [i.tensor_index_type for i in free_ind1]
2211 # Check if variance of indices needs to be fixed:
2212 pos2up = []
2213 pos2down = []
2214 free2remaining = free_ind2[:]
2215 for pos1, index1 in enumerate(free_ind1):
2216 if index1 in free2remaining:
2217 pos2 = free2remaining.index(index1)
2218 free2remaining[pos2] = None
2219 continue
2220 if -index1 in free2remaining:
2221 pos2 = free2remaining.index(-index1)
2222 free2remaining[pos2] = None
2223 free_ind2[pos2] = index1
2224 if index1.is_up:
2225 pos2up.append(pos2)
2226 else:
2227 pos2down.append(pos2)
2228 else:
2229 index2 = free2remaining[pos1]
2230 if index2 is None:
2231 raise ValueError("incompatible indices: %s and %s" % (free_ind1, free_ind2))
2232 free2remaining[pos1] = None
2233 free_ind2[pos1] = index1
2234 if index1.is_up ^ index2.is_up:
2235 if index1.is_up:
2236 pos2up.append(pos1)
2237 else:
2238 pos2down.append(pos1)
2240 if len(set(free_ind1) & set(free_ind2)) < len(free_ind1):
2241 raise ValueError("incompatible indices: %s and %s" % (free_ind1, free_ind2))
2243 # Raise indices:
2244 for pos in pos2up:
2245 index_type_pos = index_types1[pos]
2246 if index_type_pos not in replacement_dict:
2247 raise ValueError("No metric provided to lower index")
2248 metric = replacement_dict[index_type_pos]
2249 metric_inverse = _TensorDataLazyEvaluator.inverse_matrix(metric)
2250 array = TensExpr._contract_and_permute_with_metric(metric_inverse, array, pos, len(free_ind1))
2251 # Lower indices:
2252 for pos in pos2down:
2253 index_type_pos = index_types1[pos]
2254 if index_type_pos not in replacement_dict:
2255 raise ValueError("No metric provided to lower index")
2256 metric = replacement_dict[index_type_pos]
2257 array = TensExpr._contract_and_permute_with_metric(metric, array, pos, len(free_ind1))
2259 if free_ind1:
2260 permutation = TensExpr._get_indices_permutation(free_ind2, free_ind1)
2261 array = permutedims(array, permutation)
2263 if hasattr(array, "rank") and array.rank() == 0:
2264 array = array[()]
2266 return free_ind2, array
2268 def replace_with_arrays(self, replacement_dict, indices=None):
2269 """
2270 Replace the tensorial expressions with arrays. The final array will
2271 correspond to the N-dimensional array with indices arranged according
2272 to ``indices``.
2274 Parameters
2275 ==========
2277 replacement_dict
2278 dictionary containing the replacement rules for tensors.
2279 indices
2280 the index order with respect to which the array is read. The
2281 original index order will be used if no value is passed.
2283 Examples
2284 ========
2286 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices
2287 >>> from sympy.tensor.tensor import TensorHead
2288 >>> from sympy import symbols, diag
2290 >>> L = TensorIndexType("L")
2291 >>> i, j = tensor_indices("i j", L)
2292 >>> A = TensorHead("A", [L])
2293 >>> A(i).replace_with_arrays({A(i): [1, 2]}, [i])
2294 [1, 2]
2296 Since 'indices' is optional, we can also call replace_with_arrays by
2297 this way if no specific index order is needed:
2299 >>> A(i).replace_with_arrays({A(i): [1, 2]})
2300 [1, 2]
2302 >>> expr = A(i)*A(j)
2303 >>> expr.replace_with_arrays({A(i): [1, 2]})
2304 [[1, 2], [2, 4]]
2306 For contractions, specify the metric of the ``TensorIndexType``, which
2307 in this case is ``L``, in its covariant form:
2309 >>> expr = A(i)*A(-i)
2310 >>> expr.replace_with_arrays({A(i): [1, 2], L: diag(1, -1)})
2311 -3
2313 Symmetrization of an array:
2315 >>> H = TensorHead("H", [L, L])
2316 >>> a, b, c, d = symbols("a b c d")
2317 >>> expr = H(i, j)/2 + H(j, i)/2
2318 >>> expr.replace_with_arrays({H(i, j): [[a, b], [c, d]]})
2319 [[a, b/2 + c/2], [b/2 + c/2, d]]
2321 Anti-symmetrization of an array:
2323 >>> expr = H(i, j)/2 - H(j, i)/2
2324 >>> repl = {H(i, j): [[a, b], [c, d]]}
2325 >>> expr.replace_with_arrays(repl)
2326 [[0, b/2 - c/2], [-b/2 + c/2, 0]]
2328 The same expression can be read as the transpose by inverting ``i`` and
2329 ``j``:
2331 >>> expr.replace_with_arrays(repl, [j, i])
2332 [[0, -b/2 + c/2], [b/2 - c/2, 0]]
2333 """
2334 from .array import Array
2336 indices = indices or []
2337 remap = {k.args[0] if k.is_up else -k.args[0]: k for k in self.get_free_indices()}
2338 for i, index in enumerate(indices):
2339 if isinstance(index, (Symbol, Mul)):
2340 if index in remap:
2341 indices[i] = remap[index]
2342 else:
2343 indices[i] = -remap[-index]
2345 replacement_dict = {tensor: Array(array) for tensor, array in replacement_dict.items()}
2347 # Check dimensions of replaced arrays:
2348 for tensor, array in replacement_dict.items():
2349 if isinstance(tensor, TensorIndexType):
2350 expected_shape = [tensor.dim for i in range(2)]
2351 else:
2352 expected_shape = [index_type.dim for index_type in tensor.index_types]
2353 if len(expected_shape) != array.rank() or (not all(dim1 == dim2 if
2354 dim1.is_number else True for dim1, dim2 in zip(expected_shape,
2355 array.shape))):
2356 raise ValueError("shapes for tensor %s expected to be %s, "\
2357 "replacement array shape is %s" % (tensor, expected_shape,
2358 array.shape))
2360 ret_indices, array = self._extract_data(replacement_dict)
2362 last_indices, array = self._match_indices_with_other_tensor(array, indices, ret_indices, replacement_dict)
2363 return array
2365 def _check_add_Sum(self, expr, index_symbols):
2366 from sympy.concrete.summations import Sum
2367 indices = self.get_indices()
2368 dum = self.dum
2369 sum_indices = [ (index_symbols[i], 0,
2370 indices[i].tensor_index_type.dim-1) for i, j in dum]
2371 if sum_indices:
2372 expr = Sum(expr, *sum_indices)
2373 return expr
2375 def _expand_partial_derivative(self):
2376 # simply delegate the _expand_partial_derivative() to
2377 # its arguments to expand a possibly found PartialDerivative
2378 return self.func(*[
2379 a._expand_partial_derivative()
2380 if isinstance(a, TensExpr) else a
2381 for a in self.args])
2384class TensAdd(TensExpr, AssocOp):
2385 """
2386 Sum of tensors.
2388 Parameters
2389 ==========
2391 free_args : list of the free indices
2393 Attributes
2394 ==========
2396 ``args`` : tuple of addends
2397 ``rank`` : rank of the tensor
2398 ``free_args`` : list of the free indices in sorted order
2400 Examples
2401 ========
2403 >>> from sympy.tensor.tensor import TensorIndexType, tensor_heads, tensor_indices
2404 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
2405 >>> a, b = tensor_indices('a,b', Lorentz)
2406 >>> p, q = tensor_heads('p,q', [Lorentz])
2407 >>> t = p(a) + q(a); t
2408 p(a) + q(a)
2410 Examples with components data added to the tensor expression:
2412 >>> from sympy import symbols, diag
2413 >>> x, y, z, t = symbols("x y z t")
2414 >>> repl = {}
2415 >>> repl[Lorentz] = diag(1, -1, -1, -1)
2416 >>> repl[p(a)] = [1, 2, 3, 4]
2417 >>> repl[q(a)] = [x, y, z, t]
2419 The following are: 2**2 - 3**2 - 2**2 - 7**2 ==> -58
2421 >>> expr = p(a) + q(a)
2422 >>> expr.replace_with_arrays(repl, [a])
2423 [x + 1, y + 2, z + 3, t + 4]
2424 """
2426 def __new__(cls, *args, **kw_args):
2427 args = [_sympify(x) for x in args if x]
2428 args = TensAdd._tensAdd_flatten(args)
2429 args.sort(key=default_sort_key)
2430 if not args:
2431 return S.Zero
2432 if len(args) == 1:
2433 return args[0]
2435 return Basic.__new__(cls, *args, **kw_args)
2437 @property
2438 def coeff(self):
2439 return S.One
2441 @property
2442 def nocoeff(self):
2443 return self
2445 def get_free_indices(self) -> list[TensorIndex]:
2446 return self.free_indices
2448 def _replace_indices(self, repl: dict[TensorIndex, TensorIndex]) -> TensExpr:
2449 newargs = [arg._replace_indices(repl) if isinstance(arg, TensExpr) else arg for arg in self.args]
2450 return self.func(*newargs)
2452 @memoize_property
2453 def rank(self):
2454 if isinstance(self.args[0], TensExpr):
2455 return self.args[0].rank
2456 else:
2457 return 0
2459 @memoize_property
2460 def free_args(self):
2461 if isinstance(self.args[0], TensExpr):
2462 return self.args[0].free_args
2463 else:
2464 return []
2466 @memoize_property
2467 def free_indices(self):
2468 if isinstance(self.args[0], TensExpr):
2469 return self.args[0].get_free_indices()
2470 else:
2471 return set()
2473 def doit(self, **hints):
2474 deep = hints.get('deep', True)
2475 if deep:
2476 args = [arg.doit(**hints) for arg in self.args]
2477 else:
2478 args = self.args
2480 # if any of the args are zero (after doit), drop them. Otherwise, _tensAdd_check will complain about non-matching indices, even though the TensAdd is correctly formed.
2481 args = [arg for arg in args if arg != S.Zero]
2483 if len(args) == 0:
2484 return S.Zero
2485 elif len(args) == 1:
2486 return args[0]
2488 # now check that all addends have the same indices:
2489 TensAdd._tensAdd_check(args)
2491 # Collect terms appearing more than once, differing by their coefficients:
2492 args = TensAdd._tensAdd_collect_terms(args)
2494 # collect canonicalized terms
2495 def sort_key(t):
2496 if not isinstance(t, TensExpr):
2497 return [], [], []
2498 if hasattr(t, "_index_structure") and hasattr(t, "components"):
2499 x = get_index_structure(t)
2500 return t.components, x.free, x.dum
2501 return [], [], []
2502 args.sort(key=sort_key)
2504 if not args:
2505 return S.Zero
2506 # it there is only a component tensor return it
2507 if len(args) == 1:
2508 return args[0]
2510 obj = self.func(*args)
2511 return obj
2513 @staticmethod
2514 def _tensAdd_flatten(args):
2515 # flatten TensAdd, coerce terms which are not tensors to tensors
2516 a = []
2517 for x in args:
2518 if isinstance(x, (Add, TensAdd)):
2519 a.extend(list(x.args))
2520 else:
2521 a.append(x)
2522 args = [x for x in a if x.coeff]
2523 return args
2525 @staticmethod
2526 def _tensAdd_check(args):
2527 # check that all addends have the same free indices
2529 def get_indices_set(x: Expr) -> set[TensorIndex]:
2530 if isinstance(x, TensExpr):
2531 return set(x.get_free_indices())
2532 return set()
2534 indices0 = get_indices_set(args[0])
2535 list_indices = [get_indices_set(arg) for arg in args[1:]]
2536 if not all(x == indices0 for x in list_indices):
2537 raise ValueError('all tensors must have the same indices')
2539 @staticmethod
2540 def _tensAdd_collect_terms(args):
2541 # collect TensMul terms differing at most by their coefficient
2542 terms_dict = defaultdict(list)
2543 scalars = S.Zero
2544 if isinstance(args[0], TensExpr):
2545 free_indices = set(args[0].get_free_indices())
2546 else:
2547 free_indices = set()
2549 for arg in args:
2550 if not isinstance(arg, TensExpr):
2551 if free_indices != set():
2552 raise ValueError("wrong valence")
2553 scalars += arg
2554 continue
2555 if free_indices != set(arg.get_free_indices()):
2556 raise ValueError("wrong valence")
2557 # TODO: what is the part which is not a coeff?
2558 # needs an implementation similar to .as_coeff_Mul()
2559 terms_dict[arg.nocoeff].append(arg.coeff)
2561 new_args = [TensMul(Add(*coeff), t).doit() for t, coeff in terms_dict.items() if Add(*coeff) != 0]
2562 if isinstance(scalars, Add):
2563 new_args = list(scalars.args) + new_args
2564 elif scalars != 0:
2565 new_args = [scalars] + new_args
2566 return new_args
2568 def get_indices(self):
2569 indices = []
2570 for arg in self.args:
2571 indices.extend([i for i in get_indices(arg) if i not in indices])
2572 return indices
2574 def _expand(self, **hints):
2575 return TensAdd(*[_expand(i, **hints) for i in self.args])
2577 def __call__(self, *indices):
2578 deprecate_call()
2579 free_args = self.free_args
2580 indices = list(indices)
2581 if [x.tensor_index_type for x in indices] != [x.tensor_index_type for x in free_args]:
2582 raise ValueError('incompatible types')
2583 if indices == free_args:
2584 return self
2585 index_tuples = list(zip(free_args, indices))
2586 a = [x.func(*x.substitute_indices(*index_tuples).args) for x in self.args]
2587 res = TensAdd(*a).doit()
2588 return res
2590 def canon_bp(self):
2591 """
2592 Canonicalize using the Butler-Portugal algorithm for canonicalization
2593 under monoterm symmetries.
2594 """
2595 expr = self.expand()
2596 args = [canon_bp(x) for x in expr.args]
2597 res = TensAdd(*args).doit()
2598 return res
2600 def equals(self, other):
2601 other = _sympify(other)
2602 if isinstance(other, TensMul) and other.coeff == 0:
2603 return all(x.coeff == 0 for x in self.args)
2604 if isinstance(other, TensExpr):
2605 if self.rank != other.rank:
2606 return False
2607 if isinstance(other, TensAdd):
2608 if set(self.args) != set(other.args):
2609 return False
2610 else:
2611 return True
2612 t = self - other
2613 if not isinstance(t, TensExpr):
2614 return t == 0
2615 else:
2616 if isinstance(t, TensMul):
2617 return t.coeff == 0
2618 else:
2619 return all(x.coeff == 0 for x in t.args)
2621 def __getitem__(self, item):
2622 deprecate_data()
2623 with ignore_warnings(SymPyDeprecationWarning):
2624 return self.data[item]
2626 def contract_delta(self, delta):
2627 args = [x.contract_delta(delta) for x in self.args]
2628 t = TensAdd(*args).doit()
2629 return canon_bp(t)
2631 def contract_metric(self, g):
2632 """
2633 Raise or lower indices with the metric ``g``.
2635 Parameters
2636 ==========
2638 g : metric
2640 contract_all : if True, eliminate all ``g`` which are contracted
2642 Notes
2643 =====
2645 see the ``TensorIndexType`` docstring for the contraction conventions
2646 """
2648 args = [contract_metric(x, g) for x in self.args]
2649 t = TensAdd(*args).doit()
2650 return canon_bp(t)
2652 def substitute_indices(self, *index_tuples):
2653 new_args = []
2654 for arg in self.args:
2655 if isinstance(arg, TensExpr):
2656 arg = arg.substitute_indices(*index_tuples)
2657 new_args.append(arg)
2658 return TensAdd(*new_args).doit()
2660 def _print(self):
2661 a = []
2662 args = self.args
2663 for x in args:
2664 a.append(str(x))
2665 s = ' + '.join(a)
2666 s = s.replace('+ -', '- ')
2667 return s
2669 def _extract_data(self, replacement_dict):
2670 from sympy.tensor.array import Array, permutedims
2671 args_indices, arrays = zip(*[
2672 arg._extract_data(replacement_dict) if
2673 isinstance(arg, TensExpr) else ([], arg) for arg in self.args
2674 ])
2675 arrays = [Array(i) for i in arrays]
2676 ref_indices = args_indices[0]
2677 for i in range(1, len(args_indices)):
2678 indices = args_indices[i]
2679 array = arrays[i]
2680 permutation = TensMul._get_indices_permutation(indices, ref_indices)
2681 arrays[i] = permutedims(array, permutation)
2682 return ref_indices, sum(arrays, Array.zeros(*array.shape))
2684 @property
2685 def data(self):
2686 deprecate_data()
2687 with ignore_warnings(SymPyDeprecationWarning):
2688 return _tensor_data_substitution_dict[self.expand()]
2690 @data.setter
2691 def data(self, data):
2692 deprecate_data()
2693 with ignore_warnings(SymPyDeprecationWarning):
2694 _tensor_data_substitution_dict[self] = data
2696 @data.deleter
2697 def data(self):
2698 deprecate_data()
2699 with ignore_warnings(SymPyDeprecationWarning):
2700 if self in _tensor_data_substitution_dict:
2701 del _tensor_data_substitution_dict[self]
2703 def __iter__(self):
2704 deprecate_data()
2705 if not self.data:
2706 raise ValueError("No iteration on abstract tensors")
2707 return self.data.flatten().__iter__()
2709 def _eval_rewrite_as_Indexed(self, *args):
2710 return Add.fromiter(args)
2712 def _eval_partial_derivative(self, s):
2713 # Evaluation like Add
2714 list_addends = []
2715 for a in self.args:
2716 if isinstance(a, TensExpr):
2717 list_addends.append(a._eval_partial_derivative(s))
2718 # do not call diff if s is no symbol
2719 elif s._diff_wrt:
2720 list_addends.append(a._eval_derivative(s))
2722 return self.func(*list_addends)
2725class Tensor(TensExpr):
2726 """
2727 Base tensor class, i.e. this represents a tensor, the single unit to be
2728 put into an expression.
2730 Explanation
2731 ===========
2733 This object is usually created from a ``TensorHead``, by attaching indices
2734 to it. Indices preceded by a minus sign are considered contravariant,
2735 otherwise covariant.
2737 Examples
2738 ========
2740 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorHead
2741 >>> Lorentz = TensorIndexType("Lorentz", dummy_name="L")
2742 >>> mu, nu = tensor_indices('mu nu', Lorentz)
2743 >>> A = TensorHead("A", [Lorentz, Lorentz])
2744 >>> A(mu, -nu)
2745 A(mu, -nu)
2746 >>> A(mu, -mu)
2747 A(L_0, -L_0)
2749 It is also possible to use symbols instead of inidices (appropriate indices
2750 are then generated automatically).
2752 >>> from sympy import Symbol
2753 >>> x = Symbol('x')
2754 >>> A(x, mu)
2755 A(x, mu)
2756 >>> A(x, -x)
2757 A(L_0, -L_0)
2759 """
2761 is_commutative = False
2763 _index_structure = None # type: _IndexStructure
2764 args: tuple[TensorHead, Tuple]
2766 def __new__(cls, tensor_head, indices, *, is_canon_bp=False, **kw_args):
2767 indices = cls._parse_indices(tensor_head, indices)
2768 obj = Basic.__new__(cls, tensor_head, Tuple(*indices), **kw_args)
2769 obj._index_structure = _IndexStructure.from_indices(*indices)
2770 obj._free = obj._index_structure.free[:]
2771 obj._dum = obj._index_structure.dum[:]
2772 obj._ext_rank = obj._index_structure._ext_rank
2773 obj._coeff = S.One
2774 obj._nocoeff = obj
2775 obj._component = tensor_head
2776 obj._components = [tensor_head]
2777 if tensor_head.rank != len(indices):
2778 raise ValueError("wrong number of indices")
2779 obj.is_canon_bp = is_canon_bp
2780 obj._index_map = Tensor._build_index_map(indices, obj._index_structure)
2781 return obj
2783 @property
2784 def free(self):
2785 return self._free
2787 @property
2788 def dum(self):
2789 return self._dum
2791 @property
2792 def ext_rank(self):
2793 return self._ext_rank
2795 @property
2796 def coeff(self):
2797 return self._coeff
2799 @property
2800 def nocoeff(self):
2801 return self._nocoeff
2803 @property
2804 def component(self):
2805 return self._component
2807 @property
2808 def components(self):
2809 return self._components
2811 @property
2812 def head(self):
2813 return self.args[0]
2815 @property
2816 def indices(self):
2817 return self.args[1]
2819 @property
2820 def free_indices(self):
2821 return set(self._index_structure.get_free_indices())
2823 @property
2824 def index_types(self):
2825 return self.head.index_types
2827 @property
2828 def rank(self):
2829 return len(self.free_indices)
2831 @staticmethod
2832 def _build_index_map(indices, index_structure):
2833 index_map = {}
2834 for idx in indices:
2835 index_map[idx] = (indices.index(idx),)
2836 return index_map
2838 def doit(self, **hints):
2839 args, indices, free, dum = TensMul._tensMul_contract_indices([self])
2840 return args[0]
2842 @staticmethod
2843 def _parse_indices(tensor_head, indices):
2844 if not isinstance(indices, (tuple, list, Tuple)):
2845 raise TypeError("indices should be an array, got %s" % type(indices))
2846 indices = list(indices)
2847 for i, index in enumerate(indices):
2848 if isinstance(index, Symbol):
2849 indices[i] = TensorIndex(index, tensor_head.index_types[i], True)
2850 elif isinstance(index, Mul):
2851 c, e = index.as_coeff_Mul()
2852 if c == -1 and isinstance(e, Symbol):
2853 indices[i] = TensorIndex(e, tensor_head.index_types[i], False)
2854 else:
2855 raise ValueError("index not understood: %s" % index)
2856 elif not isinstance(index, TensorIndex):
2857 raise TypeError("wrong type for index: %s is %s" % (index, type(index)))
2858 return indices
2860 def _set_new_index_structure(self, im, is_canon_bp=False):
2861 indices = im.get_indices()
2862 return self._set_indices(*indices, is_canon_bp=is_canon_bp)
2864 def _set_indices(self, *indices, is_canon_bp=False, **kw_args):
2865 if len(indices) != self.ext_rank:
2866 raise ValueError("indices length mismatch")
2867 return self.func(self.args[0], indices, is_canon_bp=is_canon_bp).doit()
2869 def _get_free_indices_set(self):
2870 return {i[0] for i in self._index_structure.free}
2872 def _get_dummy_indices_set(self):
2873 dummy_pos = set(itertools.chain(*self._index_structure.dum))
2874 return {idx for i, idx in enumerate(self.args[1]) if i in dummy_pos}
2876 def _get_indices_set(self):
2877 return set(self.args[1].args)
2879 @property
2880 def free_in_args(self):
2881 return [(ind, pos, 0) for ind, pos in self.free]
2883 @property
2884 def dum_in_args(self):
2885 return [(p1, p2, 0, 0) for p1, p2 in self.dum]
2887 @property
2888 def free_args(self):
2889 return sorted([x[0] for x in self.free])
2891 def commutes_with(self, other):
2892 """
2893 :param other:
2894 :return:
2895 0 commute
2896 1 anticommute
2897 None neither commute nor anticommute
2898 """
2899 if not isinstance(other, TensExpr):
2900 return 0
2901 elif isinstance(other, Tensor):
2902 return self.component.commutes_with(other.component)
2903 return NotImplementedError
2905 def perm2tensor(self, g, is_canon_bp=False):
2906 """
2907 Returns the tensor corresponding to the permutation ``g``.
2909 For further details, see the method in ``TIDS`` with the same name.
2910 """
2911 return perm2tensor(self, g, is_canon_bp)
2913 def canon_bp(self):
2914 if self.is_canon_bp:
2915 return self
2916 expr = self.expand()
2917 g, dummies, msym = expr._index_structure.indices_canon_args()
2918 v = components_canon_args([expr.component])
2919 can = canonicalize(g, dummies, msym, *v)
2920 if can == 0:
2921 return S.Zero
2922 tensor = self.perm2tensor(can, True)
2923 return tensor
2925 def split(self):
2926 return [self]
2928 def _expand(self, **kwargs):
2929 return self
2931 def sorted_components(self):
2932 return self
2934 def get_indices(self) -> list[TensorIndex]:
2935 """
2936 Get a list of indices, corresponding to those of the tensor.
2937 """
2938 return list(self.args[1])
2940 def get_free_indices(self) -> list[TensorIndex]:
2941 """
2942 Get a list of free indices, corresponding to those of the tensor.
2943 """
2944 return self._index_structure.get_free_indices()
2946 def _replace_indices(self, repl: dict[TensorIndex, TensorIndex]) -> TensExpr:
2947 # TODO: this could be optimized by only swapping the indices
2948 # instead of visiting the whole expression tree:
2949 return self.xreplace(repl)
2951 def as_base_exp(self):
2952 return self, S.One
2954 def substitute_indices(self, *index_tuples):
2955 """
2956 Return a tensor with free indices substituted according to ``index_tuples``.
2958 ``index_types`` list of tuples ``(old_index, new_index)``.
2960 Examples
2961 ========
2963 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensor_heads, TensorSymmetry
2964 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
2965 >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz)
2966 >>> A, B = tensor_heads('A,B', [Lorentz]*2, TensorSymmetry.fully_symmetric(2))
2967 >>> t = A(i, k)*B(-k, -j); t
2968 A(i, L_0)*B(-L_0, -j)
2969 >>> t.substitute_indices((i, k),(-j, l))
2970 A(k, L_0)*B(-L_0, l)
2971 """
2972 indices = []
2973 for index in self.indices:
2974 for ind_old, ind_new in index_tuples:
2975 if (index.name == ind_old.name and index.tensor_index_type ==
2976 ind_old.tensor_index_type):
2977 if index.is_up == ind_old.is_up:
2978 indices.append(ind_new)
2979 else:
2980 indices.append(-ind_new)
2981 break
2982 else:
2983 indices.append(index)
2984 return self.head(*indices)
2986 def __call__(self, *indices):
2987 deprecate_call()
2988 free_args = self.free_args
2989 indices = list(indices)
2990 if [x.tensor_index_type for x in indices] != [x.tensor_index_type for x in free_args]:
2991 raise ValueError('incompatible types')
2992 if indices == free_args:
2993 return self
2994 t = self.substitute_indices(*list(zip(free_args, indices)))
2996 # object is rebuilt in order to make sure that all contracted indices
2997 # get recognized as dummies, but only if there are contracted indices.
2998 if len({i if i.is_up else -i for i in indices}) != len(indices):
2999 return t.func(*t.args)
3000 return t
3002 # TODO: put this into TensExpr?
3003 def __iter__(self):
3004 deprecate_data()
3005 with ignore_warnings(SymPyDeprecationWarning):
3006 return self.data.__iter__()
3008 # TODO: put this into TensExpr?
3009 def __getitem__(self, item):
3010 deprecate_data()
3011 with ignore_warnings(SymPyDeprecationWarning):
3012 return self.data[item]
3014 def _extract_data(self, replacement_dict):
3015 from .array import Array
3016 for k, v in replacement_dict.items():
3017 if isinstance(k, Tensor) and k.args[0] == self.args[0]:
3018 other = k
3019 array = v
3020 break
3021 else:
3022 raise ValueError("%s not found in %s" % (self, replacement_dict))
3024 # TODO: inefficient, this should be done at root level only:
3025 replacement_dict = {k: Array(v) for k, v in replacement_dict.items()}
3026 array = Array(array)
3028 dum1 = self.dum
3029 dum2 = other.dum
3031 if len(dum2) > 0:
3032 for pair in dum2:
3033 # allow `dum2` if the contained values are also in `dum1`.
3034 if pair not in dum1:
3035 raise NotImplementedError("%s with contractions is not implemented" % other)
3036 # Remove elements in `dum2` from `dum1`:
3037 dum1 = [pair for pair in dum1 if pair not in dum2]
3038 if len(dum1) > 0:
3039 indices1 = self.get_indices()
3040 indices2 = other.get_indices()
3041 repl = {}
3042 for p1, p2 in dum1:
3043 repl[indices2[p2]] = -indices2[p1]
3044 for pos in (p1, p2):
3045 if indices1[pos].is_up ^ indices2[pos].is_up:
3046 metric = replacement_dict[indices1[pos].tensor_index_type]
3047 if indices1[pos].is_up:
3048 metric = _TensorDataLazyEvaluator.inverse_matrix(metric)
3049 array = self._contract_and_permute_with_metric(metric, array, pos, len(indices2))
3050 other = other.xreplace(repl).doit()
3051 array = _TensorDataLazyEvaluator.data_contract_dum([array], dum1, len(indices2))
3053 free_ind1 = self.get_free_indices()
3054 free_ind2 = other.get_free_indices()
3056 return self._match_indices_with_other_tensor(array, free_ind1, free_ind2, replacement_dict)
3058 @property
3059 def data(self):
3060 deprecate_data()
3061 with ignore_warnings(SymPyDeprecationWarning):
3062 return _tensor_data_substitution_dict[self]
3064 @data.setter
3065 def data(self, data):
3066 deprecate_data()
3067 # TODO: check data compatibility with properties of tensor.
3068 with ignore_warnings(SymPyDeprecationWarning):
3069 _tensor_data_substitution_dict[self] = data
3071 @data.deleter
3072 def data(self):
3073 deprecate_data()
3074 with ignore_warnings(SymPyDeprecationWarning):
3075 if self in _tensor_data_substitution_dict:
3076 del _tensor_data_substitution_dict[self]
3077 if self.metric in _tensor_data_substitution_dict:
3078 del _tensor_data_substitution_dict[self.metric]
3080 def _print(self):
3081 indices = [str(ind) for ind in self.indices]
3082 component = self.component
3083 if component.rank > 0:
3084 return ('%s(%s)' % (component.name, ', '.join(indices)))
3085 else:
3086 return ('%s' % component.name)
3088 def equals(self, other):
3089 if other == 0:
3090 return self.coeff == 0
3091 other = _sympify(other)
3092 if not isinstance(other, TensExpr):
3093 assert not self.components
3094 return S.One == other
3096 def _get_compar_comp(self):
3097 t = self.canon_bp()
3098 r = (t.coeff, tuple(t.components), \
3099 tuple(sorted(t.free)), tuple(sorted(t.dum)))
3100 return r
3102 return _get_compar_comp(self) == _get_compar_comp(other)
3104 def contract_metric(self, g):
3105 # if metric is not the same, ignore this step:
3106 if self.component != g:
3107 return self
3108 # in case there are free components, do not perform anything:
3109 if len(self.free) != 0:
3110 return self
3112 #antisym = g.index_types[0].metric_antisym
3113 if g.symmetry == TensorSymmetry.fully_symmetric(-2):
3114 antisym = 1
3115 elif g.symmetry == TensorSymmetry.fully_symmetric(2):
3116 antisym = 0
3117 elif g.symmetry == TensorSymmetry.no_symmetry(2):
3118 antisym = None
3119 else:
3120 raise NotImplementedError
3121 sign = S.One
3122 typ = g.index_types[0]
3124 if not antisym:
3125 # g(i, -i)
3126 sign = sign*typ.dim
3127 else:
3128 # g(i, -i)
3129 sign = sign*typ.dim
3131 dp0, dp1 = self.dum[0]
3132 if dp0 < dp1:
3133 # g(i, -i) = -D with antisymmetric metric
3134 sign = -sign
3136 return sign
3138 def contract_delta(self, metric):
3139 return self.contract_metric(metric)
3141 def _eval_rewrite_as_Indexed(self, tens, indices):
3142 from sympy.tensor.indexed import Indexed
3143 # TODO: replace .args[0] with .name:
3144 index_symbols = [i.args[0] for i in self.get_indices()]
3145 expr = Indexed(tens.args[0], *index_symbols)
3146 return self._check_add_Sum(expr, index_symbols)
3148 def _eval_partial_derivative(self, s): # type: (Tensor) -> Expr
3150 if not isinstance(s, Tensor):
3151 return S.Zero
3152 else:
3154 # @a_i/@a_k = delta_i^k
3155 # @a_i/@a^k = g_ij delta^j_k
3156 # @a^i/@a^k = delta^i_k
3157 # @a^i/@a_k = g^ij delta_j^k
3158 # TODO: if there is no metric present, the derivative should be zero?
3160 if self.head != s.head:
3161 return S.Zero
3163 # if heads are the same, provide delta and/or metric products
3164 # for every free index pair in the appropriate tensor
3165 # assumed that the free indices are in proper order
3166 # A contravariante index in the derivative becomes covariant
3167 # after performing the derivative and vice versa
3169 kronecker_delta_list = [1]
3171 # not guarantee a correct index order
3173 for (count, (iself, iother)) in enumerate(zip(self.get_free_indices(), s.get_free_indices())):
3174 if iself.tensor_index_type != iother.tensor_index_type:
3175 raise ValueError("index types not compatible")
3176 else:
3177 tensor_index_type = iself.tensor_index_type
3178 tensor_metric = tensor_index_type.metric
3179 dummy = TensorIndex("d_" + str(count), tensor_index_type,
3180 is_up=iself.is_up)
3181 if iself.is_up == iother.is_up:
3182 kroneckerdelta = tensor_index_type.delta(iself, -iother)
3183 else:
3184 kroneckerdelta = (
3185 TensMul(tensor_metric(iself, dummy),
3186 tensor_index_type.delta(-dummy, -iother))
3187 )
3188 kronecker_delta_list.append(kroneckerdelta)
3189 return TensMul.fromiter(kronecker_delta_list).doit()
3190 # doit necessary to rename dummy indices accordingly
3193class TensMul(TensExpr, AssocOp):
3194 """
3195 Product of tensors.
3197 Parameters
3198 ==========
3200 coeff : SymPy coefficient of the tensor
3201 args
3203 Attributes
3204 ==========
3206 ``components`` : list of ``TensorHead`` of the component tensors
3207 ``types`` : list of nonrepeated ``TensorIndexType``
3208 ``free`` : list of ``(ind, ipos, icomp)``, see Notes
3209 ``dum`` : list of ``(ipos1, ipos2, icomp1, icomp2)``, see Notes
3210 ``ext_rank`` : rank of the tensor counting the dummy indices
3211 ``rank`` : rank of the tensor
3212 ``coeff`` : SymPy coefficient of the tensor
3213 ``free_args`` : list of the free indices in sorted order
3214 ``is_canon_bp`` : ``True`` if the tensor in in canonical form
3216 Notes
3217 =====
3219 ``args[0]`` list of ``TensorHead`` of the component tensors.
3221 ``args[1]`` list of ``(ind, ipos, icomp)``
3222 where ``ind`` is a free index, ``ipos`` is the slot position
3223 of ``ind`` in the ``icomp``-th component tensor.
3225 ``args[2]`` list of tuples representing dummy indices.
3226 ``(ipos1, ipos2, icomp1, icomp2)`` indicates that the contravariant
3227 dummy index is the ``ipos1``-th slot position in the ``icomp1``-th
3228 component tensor; the corresponding covariant index is
3229 in the ``ipos2`` slot position in the ``icomp2``-th component tensor.
3231 """
3232 identity = S.One
3234 _index_structure = None # type: _IndexStructure
3236 def __new__(cls, *args, **kw_args):
3237 is_canon_bp = kw_args.get('is_canon_bp', False)
3238 args = list(map(_sympify, args))
3240 """
3241 If the internal dummy indices in one arg conflict with the free indices
3242 of the remaining args, we need to rename those internal dummy indices.
3243 """
3244 free = [get_free_indices(arg) for arg in args]
3245 free = set(itertools.chain(*free)) #flatten free
3246 newargs = []
3247 for arg in args:
3248 dum_this = set(get_dummy_indices(arg))
3249 dum_other = [get_dummy_indices(a) for a in newargs]
3250 dum_other = set(itertools.chain(*dum_other)) #flatten dum_other
3251 free_this = set(get_free_indices(arg))
3252 if len(dum_this.intersection(free)) > 0:
3253 exclude = free_this.union(free, dum_other)
3254 newarg = TensMul._dedupe_indices(arg, exclude)
3255 else:
3256 newarg = arg
3257 newargs.append(newarg)
3259 args = newargs
3261 # Flatten:
3262 args = [i for arg in args for i in (arg.args if isinstance(arg, (TensMul, Mul)) else [arg])]
3264 args, indices, free, dum = TensMul._tensMul_contract_indices(args, replace_indices=False)
3266 # Data for indices:
3267 index_types = [i.tensor_index_type for i in indices]
3268 index_structure = _IndexStructure(free, dum, index_types, indices, canon_bp=is_canon_bp)
3270 obj = TensExpr.__new__(cls, *args)
3271 obj._indices = indices
3272 obj._index_types = index_types[:]
3273 obj._index_structure = index_structure
3274 obj._free = index_structure.free[:]
3275 obj._dum = index_structure.dum[:]
3276 obj._free_indices = {x[0] for x in obj.free}
3277 obj._rank = len(obj.free)
3278 obj._ext_rank = len(obj._index_structure.free) + 2*len(obj._index_structure.dum)
3279 obj._coeff = S.One
3280 obj._is_canon_bp = is_canon_bp
3281 return obj
3283 index_types = property(lambda self: self._index_types)
3284 free = property(lambda self: self._free)
3285 dum = property(lambda self: self._dum)
3286 free_indices = property(lambda self: self._free_indices)
3287 rank = property(lambda self: self._rank)
3288 ext_rank = property(lambda self: self._ext_rank)
3290 @staticmethod
3291 def _indices_to_free_dum(args_indices):
3292 free2pos1 = {}
3293 free2pos2 = {}
3294 dummy_data = []
3295 indices = []
3297 # Notation for positions (to better understand the code):
3298 # `pos1`: position in the `args`.
3299 # `pos2`: position in the indices.
3301 # Example:
3302 # A(i, j)*B(k, m, n)*C(p)
3303 # `pos1` of `n` is 1 because it's in `B` (second `args` of TensMul).
3304 # `pos2` of `n` is 4 because it's the fifth overall index.
3306 # Counter for the index position wrt the whole expression:
3307 pos2 = 0
3309 for pos1, arg_indices in enumerate(args_indices):
3311 for index_pos, index in enumerate(arg_indices):
3312 if not isinstance(index, TensorIndex):
3313 raise TypeError("expected TensorIndex")
3314 if -index in free2pos1:
3315 # Dummy index detected:
3316 other_pos1 = free2pos1.pop(-index)
3317 other_pos2 = free2pos2.pop(-index)
3318 if index.is_up:
3319 dummy_data.append((index, pos1, other_pos1, pos2, other_pos2))
3320 else:
3321 dummy_data.append((-index, other_pos1, pos1, other_pos2, pos2))
3322 indices.append(index)
3323 elif index in free2pos1:
3324 raise ValueError("Repeated index: %s" % index)
3325 else:
3326 free2pos1[index] = pos1
3327 free2pos2[index] = pos2
3328 indices.append(index)
3329 pos2 += 1
3331 free = [(i, p) for (i, p) in free2pos2.items()]
3332 free_names = [i.name for i in free2pos2.keys()]
3334 dummy_data.sort(key=lambda x: x[3])
3335 return indices, free, free_names, dummy_data
3337 @staticmethod
3338 def _dummy_data_to_dum(dummy_data):
3339 return [(p2a, p2b) for (i, p1a, p1b, p2a, p2b) in dummy_data]
3341 @staticmethod
3342 def _tensMul_contract_indices(args, replace_indices=True):
3343 replacements = [{} for _ in args]
3345 #_index_order = all(_has_index_order(arg) for arg in args)
3347 args_indices = [get_indices(arg) for arg in args]
3348 indices, free, free_names, dummy_data = TensMul._indices_to_free_dum(args_indices)
3350 cdt = defaultdict(int)
3352 def dummy_name_gen(tensor_index_type):
3353 nd = str(cdt[tensor_index_type])
3354 cdt[tensor_index_type] += 1
3355 return tensor_index_type.dummy_name + '_' + nd
3357 if replace_indices:
3358 for old_index, pos1cov, pos1contra, pos2cov, pos2contra in dummy_data:
3359 index_type = old_index.tensor_index_type
3360 while True:
3361 dummy_name = dummy_name_gen(index_type)
3362 if dummy_name not in free_names:
3363 break
3364 dummy = TensorIndex(dummy_name, index_type, True)
3365 replacements[pos1cov][old_index] = dummy
3366 replacements[pos1contra][-old_index] = -dummy
3367 indices[pos2cov] = dummy
3368 indices[pos2contra] = -dummy
3369 args = [
3370 arg._replace_indices(repl) if isinstance(arg, TensExpr) else arg
3371 for arg, repl in zip(args, replacements)]
3373 dum = TensMul._dummy_data_to_dum(dummy_data)
3374 return args, indices, free, dum
3376 @staticmethod
3377 def _get_components_from_args(args):
3378 """
3379 Get a list of ``Tensor`` objects having the same ``TIDS`` if multiplied
3380 by one another.
3381 """
3382 components = []
3383 for arg in args:
3384 if not isinstance(arg, TensExpr):
3385 continue
3386 if isinstance(arg, TensAdd):
3387 continue
3388 components.extend(arg.components)
3389 return components
3391 @staticmethod
3392 def _rebuild_tensors_list(args, index_structure):
3393 indices = index_structure.get_indices()
3394 #tensors = [None for i in components] # pre-allocate list
3395 ind_pos = 0
3396 for i, arg in enumerate(args):
3397 if not isinstance(arg, TensExpr):
3398 continue
3399 prev_pos = ind_pos
3400 ind_pos += arg.ext_rank
3401 args[i] = Tensor(arg.component, indices[prev_pos:ind_pos])
3403 def doit(self, **hints):
3404 is_canon_bp = self._is_canon_bp
3405 deep = hints.get('deep', True)
3406 if deep:
3407 args = [arg.doit(**hints) for arg in self.args]
3409 """
3410 There may now be conflicts between dummy indices of different args
3411 (each arg's doit method does not have any information about which
3412 dummy indices are already used in the other args), so we
3413 deduplicate them.
3414 """
3415 rule = dict(zip(self.args, args))
3416 rule = self._dedupe_indices_in_rule(rule)
3417 args = [rule[a] for a in self.args]
3419 else:
3420 args = self.args
3422 args = [arg for arg in args if arg != self.identity]
3424 # Extract non-tensor coefficients:
3425 coeff = reduce(lambda a, b: a*b, [arg for arg in args if not isinstance(arg, TensExpr)], S.One)
3426 args = [arg for arg in args if isinstance(arg, TensExpr)]
3428 if len(args) == 0:
3429 return coeff
3431 if coeff != self.identity:
3432 args = [coeff] + args
3433 if coeff == 0:
3434 return S.Zero
3436 if len(args) == 1:
3437 return args[0]
3439 args, indices, free, dum = TensMul._tensMul_contract_indices(args)
3441 # Data for indices:
3442 index_types = [i.tensor_index_type for i in indices]
3443 index_structure = _IndexStructure(free, dum, index_types, indices, canon_bp=is_canon_bp)
3445 obj = self.func(*args)
3446 obj._index_types = index_types
3447 obj._index_structure = index_structure
3448 obj._ext_rank = len(obj._index_structure.free) + 2*len(obj._index_structure.dum)
3449 obj._coeff = coeff
3450 obj._is_canon_bp = is_canon_bp
3451 return obj
3453 # TODO: this method should be private
3454 # TODO: should this method be renamed _from_components_free_dum ?
3455 @staticmethod
3456 def from_data(coeff, components, free, dum, **kw_args):
3457 return TensMul(coeff, *TensMul._get_tensors_from_components_free_dum(components, free, dum), **kw_args).doit()
3459 @staticmethod
3460 def _get_tensors_from_components_free_dum(components, free, dum):
3461 """
3462 Get a list of ``Tensor`` objects by distributing ``free`` and ``dum`` indices on the ``components``.
3463 """
3464 index_structure = _IndexStructure.from_components_free_dum(components, free, dum)
3465 indices = index_structure.get_indices()
3466 tensors = [None for i in components] # pre-allocate list
3468 # distribute indices on components to build a list of tensors:
3469 ind_pos = 0
3470 for i, component in enumerate(components):
3471 prev_pos = ind_pos
3472 ind_pos += component.rank
3473 tensors[i] = Tensor(component, indices[prev_pos:ind_pos])
3474 return tensors
3476 def _get_free_indices_set(self):
3477 return {i[0] for i in self.free}
3479 def _get_dummy_indices_set(self):
3480 dummy_pos = set(itertools.chain(*self.dum))
3481 return {idx for i, idx in enumerate(self._index_structure.get_indices()) if i in dummy_pos}
3483 def _get_position_offset_for_indices(self):
3484 arg_offset = [None for i in range(self.ext_rank)]
3485 counter = 0
3486 for i, arg in enumerate(self.args):
3487 if not isinstance(arg, TensExpr):
3488 continue
3489 for j in range(arg.ext_rank):
3490 arg_offset[j + counter] = counter
3491 counter += arg.ext_rank
3492 return arg_offset
3494 @property
3495 def free_args(self):
3496 return sorted([x[0] for x in self.free])
3498 @property
3499 def components(self):
3500 return self._get_components_from_args(self.args)
3502 @property
3503 def free_in_args(self):
3504 arg_offset = self._get_position_offset_for_indices()
3505 argpos = self._get_indices_to_args_pos()
3506 return [(ind, pos-arg_offset[pos], argpos[pos]) for (ind, pos) in self.free]
3508 @property
3509 def coeff(self):
3510 # return Mul.fromiter([c for c in self.args if not isinstance(c, TensExpr)])
3511 return self._coeff
3513 @property
3514 def nocoeff(self):
3515 return self.func(*[t for t in self.args if isinstance(t, TensExpr)]).doit()
3517 @property
3518 def dum_in_args(self):
3519 arg_offset = self._get_position_offset_for_indices()
3520 argpos = self._get_indices_to_args_pos()
3521 return [(p1-arg_offset[p1], p2-arg_offset[p2], argpos[p1], argpos[p2]) for p1, p2 in self.dum]
3523 def equals(self, other):
3524 if other == 0:
3525 return self.coeff == 0
3526 other = _sympify(other)
3527 if not isinstance(other, TensExpr):
3528 assert not self.components
3529 return self.coeff == other
3531 return self.canon_bp() == other.canon_bp()
3533 def get_indices(self):
3534 """
3535 Returns the list of indices of the tensor.
3537 Explanation
3538 ===========
3540 The indices are listed in the order in which they appear in the
3541 component tensors.
3542 The dummy indices are given a name which does not collide with
3543 the names of the free indices.
3545 Examples
3546 ========
3548 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensor_heads
3549 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
3550 >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz)
3551 >>> g = Lorentz.metric
3552 >>> p, q = tensor_heads('p,q', [Lorentz])
3553 >>> t = p(m1)*g(m0,m2)
3554 >>> t.get_indices()
3555 [m1, m0, m2]
3556 >>> t2 = p(m1)*g(-m1, m2)
3557 >>> t2.get_indices()
3558 [L_0, -L_0, m2]
3559 """
3560 return self._indices
3562 def get_free_indices(self) -> list[TensorIndex]:
3563 """
3564 Returns the list of free indices of the tensor.
3566 Explanation
3567 ===========
3569 The indices are listed in the order in which they appear in the
3570 component tensors.
3572 Examples
3573 ========
3575 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensor_heads
3576 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
3577 >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz)
3578 >>> g = Lorentz.metric
3579 >>> p, q = tensor_heads('p,q', [Lorentz])
3580 >>> t = p(m1)*g(m0,m2)
3581 >>> t.get_free_indices()
3582 [m1, m0, m2]
3583 >>> t2 = p(m1)*g(-m1, m2)
3584 >>> t2.get_free_indices()
3585 [m2]
3586 """
3587 return self._index_structure.get_free_indices()
3589 def _replace_indices(self, repl: dict[TensorIndex, TensorIndex]) -> TensExpr:
3590 return self.func(*[arg._replace_indices(repl) if isinstance(arg, TensExpr) else arg for arg in self.args])
3592 def split(self):
3593 """
3594 Returns a list of tensors, whose product is ``self``.
3596 Explanation
3597 ===========
3599 Dummy indices contracted among different tensor components
3600 become free indices with the same name as the one used to
3601 represent the dummy indices.
3603 Examples
3604 ========
3606 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensor_heads, TensorSymmetry
3607 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
3608 >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz)
3609 >>> A, B = tensor_heads('A,B', [Lorentz]*2, TensorSymmetry.fully_symmetric(2))
3610 >>> t = A(a,b)*B(-b,c)
3611 >>> t
3612 A(a, L_0)*B(-L_0, c)
3613 >>> t.split()
3614 [A(a, L_0), B(-L_0, c)]
3615 """
3616 if self.args == ():
3617 return [self]
3618 splitp = []
3619 res = 1
3620 for arg in self.args:
3621 if isinstance(arg, Tensor):
3622 splitp.append(res*arg)
3623 res = 1
3624 else:
3625 res *= arg
3626 return splitp
3628 def _expand(self, **hints):
3629 # TODO: temporary solution, in the future this should be linked to
3630 # `Expr.expand`.
3631 args = [_expand(arg, **hints) for arg in self.args]
3632 args1 = [arg.args if isinstance(arg, (Add, TensAdd)) else (arg,) for arg in args]
3633 return TensAdd(*[
3634 TensMul(*i) for i in itertools.product(*args1)]
3635 )
3637 def __neg__(self):
3638 return TensMul(S.NegativeOne, self, is_canon_bp=self._is_canon_bp).doit()
3640 def __getitem__(self, item):
3641 deprecate_data()
3642 with ignore_warnings(SymPyDeprecationWarning):
3643 return self.data[item]
3645 def _get_args_for_traditional_printer(self):
3646 args = list(self.args)
3647 if self.coeff.could_extract_minus_sign():
3648 # expressions like "-A(a)"
3649 sign = "-"
3650 if args[0] == S.NegativeOne:
3651 args = args[1:]
3652 else:
3653 args[0] = -args[0]
3654 else:
3655 sign = ""
3656 return sign, args
3658 def _sort_args_for_sorted_components(self):
3659 """
3660 Returns the ``args`` sorted according to the components commutation
3661 properties.
3663 Explanation
3664 ===========
3666 The sorting is done taking into account the commutation group
3667 of the component tensors.
3668 """
3669 cv = [arg for arg in self.args if isinstance(arg, TensExpr)]
3670 sign = 1
3671 n = len(cv) - 1
3672 for i in range(n):
3673 for j in range(n, i, -1):
3674 c = cv[j-1].commutes_with(cv[j])
3675 # if `c` is `None`, it does neither commute nor anticommute, skip:
3676 if c not in (0, 1):
3677 continue
3678 typ1 = sorted(set(cv[j-1].component.index_types), key=lambda x: x.name)
3679 typ2 = sorted(set(cv[j].component.index_types), key=lambda x: x.name)
3680 if (typ1, cv[j-1].component.name) > (typ2, cv[j].component.name):
3681 cv[j-1], cv[j] = cv[j], cv[j-1]
3682 # if `c` is 1, the anticommute, so change sign:
3683 if c:
3684 sign = -sign
3686 coeff = sign * self.coeff
3687 if coeff != 1:
3688 return [coeff] + cv
3689 return cv
3691 def sorted_components(self):
3692 """
3693 Returns a tensor product with sorted components.
3694 """
3695 return TensMul(*self._sort_args_for_sorted_components()).doit()
3697 def perm2tensor(self, g, is_canon_bp=False):
3698 """
3699 Returns the tensor corresponding to the permutation ``g``
3701 For further details, see the method in ``TIDS`` with the same name.
3702 """
3703 return perm2tensor(self, g, is_canon_bp=is_canon_bp)
3705 def canon_bp(self):
3706 """
3707 Canonicalize using the Butler-Portugal algorithm for canonicalization
3708 under monoterm symmetries.
3710 Examples
3711 ========
3713 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorHead, TensorSymmetry
3714 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
3715 >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz)
3716 >>> A = TensorHead('A', [Lorentz]*2, TensorSymmetry.fully_symmetric(-2))
3717 >>> t = A(m0,-m1)*A(m1,-m0)
3718 >>> t.canon_bp()
3719 -A(L_0, L_1)*A(-L_0, -L_1)
3720 >>> t = A(m0,-m1)*A(m1,-m2)*A(m2,-m0)
3721 >>> t.canon_bp()
3722 0
3723 """
3724 if self._is_canon_bp:
3725 return self
3726 expr = self.expand()
3727 if isinstance(expr, TensAdd):
3728 return expr.canon_bp()
3729 if not expr.components:
3730 return expr
3731 t = expr.sorted_components()
3732 g, dummies, msym = t._index_structure.indices_canon_args()
3733 v = components_canon_args(t.components)
3734 can = canonicalize(g, dummies, msym, *v)
3735 if can == 0:
3736 return S.Zero
3737 tmul = t.perm2tensor(can, True)
3738 return tmul
3740 def contract_delta(self, delta):
3741 t = self.contract_metric(delta)
3742 return t
3744 def _get_indices_to_args_pos(self):
3745 """
3746 Get a dict mapping the index position to TensMul's argument number.
3747 """
3748 pos_map = {}
3749 pos_counter = 0
3750 for arg_i, arg in enumerate(self.args):
3751 if not isinstance(arg, TensExpr):
3752 continue
3753 assert isinstance(arg, Tensor)
3754 for i in range(arg.ext_rank):
3755 pos_map[pos_counter] = arg_i
3756 pos_counter += 1
3757 return pos_map
3759 def contract_metric(self, g):
3760 """
3761 Raise or lower indices with the metric ``g``.
3763 Parameters
3764 ==========
3766 g : metric
3768 Notes
3769 =====
3771 See the ``TensorIndexType`` docstring for the contraction conventions.
3773 Examples
3774 ========
3776 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensor_heads
3777 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
3778 >>> m0, m1, m2 = tensor_indices('m0,m1,m2', Lorentz)
3779 >>> g = Lorentz.metric
3780 >>> p, q = tensor_heads('p,q', [Lorentz])
3781 >>> t = p(m0)*q(m1)*g(-m0, -m1)
3782 >>> t.canon_bp()
3783 metric(L_0, L_1)*p(-L_0)*q(-L_1)
3784 >>> t.contract_metric(g).canon_bp()
3785 p(L_0)*q(-L_0)
3786 """
3787 expr = self.expand()
3788 if self != expr:
3789 expr = canon_bp(expr)
3790 return contract_metric(expr, g)
3791 pos_map = self._get_indices_to_args_pos()
3792 args = list(self.args)
3794 #antisym = g.index_types[0].metric_antisym
3795 if g.symmetry == TensorSymmetry.fully_symmetric(-2):
3796 antisym = 1
3797 elif g.symmetry == TensorSymmetry.fully_symmetric(2):
3798 antisym = 0
3799 elif g.symmetry == TensorSymmetry.no_symmetry(2):
3800 antisym = None
3801 else:
3802 raise NotImplementedError
3804 # list of positions of the metric ``g`` inside ``args``
3805 gpos = [i for i, x in enumerate(self.args) if isinstance(x, Tensor) and x.component == g]
3806 if not gpos:
3807 return self
3809 # Sign is either 1 or -1, to correct the sign after metric contraction
3810 # (for spinor indices).
3811 sign = 1
3812 dum = self.dum[:]
3813 free = self.free[:]
3814 elim = set()
3815 for gposx in gpos:
3816 if gposx in elim:
3817 continue
3818 free1 = [x for x in free if pos_map[x[1]] == gposx]
3819 dum1 = [x for x in dum if pos_map[x[0]] == gposx or pos_map[x[1]] == gposx]
3820 if not dum1:
3821 continue
3822 elim.add(gposx)
3823 # subs with the multiplication neutral element, that is, remove it:
3824 args[gposx] = 1
3825 if len(dum1) == 2:
3826 if not antisym:
3827 dum10, dum11 = dum1
3828 if pos_map[dum10[1]] == gposx:
3829 # the index with pos p0 contravariant
3830 p0 = dum10[0]
3831 else:
3832 # the index with pos p0 is covariant
3833 p0 = dum10[1]
3834 if pos_map[dum11[1]] == gposx:
3835 # the index with pos p1 is contravariant
3836 p1 = dum11[0]
3837 else:
3838 # the index with pos p1 is covariant
3839 p1 = dum11[1]
3841 dum.append((p0, p1))
3842 else:
3843 dum10, dum11 = dum1
3844 # change the sign to bring the indices of the metric to contravariant
3845 # form; change the sign if dum10 has the metric index in position 0
3846 if pos_map[dum10[1]] == gposx:
3847 # the index with pos p0 is contravariant
3848 p0 = dum10[0]
3849 if dum10[1] == 1:
3850 sign = -sign
3851 else:
3852 # the index with pos p0 is covariant
3853 p0 = dum10[1]
3854 if dum10[0] == 0:
3855 sign = -sign
3856 if pos_map[dum11[1]] == gposx:
3857 # the index with pos p1 is contravariant
3858 p1 = dum11[0]
3859 sign = -sign
3860 else:
3861 # the index with pos p1 is covariant
3862 p1 = dum11[1]
3864 dum.append((p0, p1))
3866 elif len(dum1) == 1:
3867 if not antisym:
3868 dp0, dp1 = dum1[0]
3869 if pos_map[dp0] == pos_map[dp1]:
3870 # g(i, -i)
3871 typ = g.index_types[0]
3872 sign = sign*typ.dim
3874 else:
3875 # g(i0, i1)*p(-i1)
3876 if pos_map[dp0] == gposx:
3877 p1 = dp1
3878 else:
3879 p1 = dp0
3881 ind, p = free1[0]
3882 free.append((ind, p1))
3883 else:
3884 dp0, dp1 = dum1[0]
3885 if pos_map[dp0] == pos_map[dp1]:
3886 # g(i, -i)
3887 typ = g.index_types[0]
3888 sign = sign*typ.dim
3890 if dp0 < dp1:
3891 # g(i, -i) = -D with antisymmetric metric
3892 sign = -sign
3893 else:
3894 # g(i0, i1)*p(-i1)
3895 if pos_map[dp0] == gposx:
3896 p1 = dp1
3897 if dp0 == 0:
3898 sign = -sign
3899 else:
3900 p1 = dp0
3901 ind, p = free1[0]
3902 free.append((ind, p1))
3903 dum = [x for x in dum if x not in dum1]
3904 free = [x for x in free if x not in free1]
3906 # shift positions:
3907 shift = 0
3908 shifts = [0]*len(args)
3909 for i in range(len(args)):
3910 if i in elim:
3911 shift += 2
3912 continue
3913 shifts[i] = shift
3914 free = [(ind, p - shifts[pos_map[p]]) for (ind, p) in free if pos_map[p] not in elim]
3915 dum = [(p0 - shifts[pos_map[p0]], p1 - shifts[pos_map[p1]]) for i, (p0, p1) in enumerate(dum) if pos_map[p0] not in elim and pos_map[p1] not in elim]
3917 res = sign*TensMul(*args).doit()
3918 if not isinstance(res, TensExpr):
3919 return res
3920 im = _IndexStructure.from_components_free_dum(res.components, free, dum)
3921 return res._set_new_index_structure(im)
3923 def _set_new_index_structure(self, im, is_canon_bp=False):
3924 indices = im.get_indices()
3925 return self._set_indices(*indices, is_canon_bp=is_canon_bp)
3927 def _set_indices(self, *indices, is_canon_bp=False, **kw_args):
3928 if len(indices) != self.ext_rank:
3929 raise ValueError("indices length mismatch")
3930 args = list(self.args)[:]
3931 pos = 0
3932 for i, arg in enumerate(args):
3933 if not isinstance(arg, TensExpr):
3934 continue
3935 assert isinstance(arg, Tensor)
3936 ext_rank = arg.ext_rank
3937 args[i] = arg._set_indices(*indices[pos:pos+ext_rank])
3938 pos += ext_rank
3939 return TensMul(*args, is_canon_bp=is_canon_bp).doit()
3941 @staticmethod
3942 def _index_replacement_for_contract_metric(args, free, dum):
3943 for arg in args:
3944 if not isinstance(arg, TensExpr):
3945 continue
3946 assert isinstance(arg, Tensor)
3948 def substitute_indices(self, *index_tuples):
3949 new_args = []
3950 for arg in self.args:
3951 if isinstance(arg, TensExpr):
3952 arg = arg.substitute_indices(*index_tuples)
3953 new_args.append(arg)
3954 return TensMul(*new_args).doit()
3956 def __call__(self, *indices):
3957 deprecate_call()
3958 free_args = self.free_args
3959 indices = list(indices)
3960 if [x.tensor_index_type for x in indices] != [x.tensor_index_type for x in free_args]:
3961 raise ValueError('incompatible types')
3962 if indices == free_args:
3963 return self
3964 t = self.substitute_indices(*list(zip(free_args, indices)))
3966 # object is rebuilt in order to make sure that all contracted indices
3967 # get recognized as dummies, but only if there are contracted indices.
3968 if len({i if i.is_up else -i for i in indices}) != len(indices):
3969 return t.func(*t.args)
3970 return t
3972 def _extract_data(self, replacement_dict):
3973 args_indices, arrays = zip(*[arg._extract_data(replacement_dict) for arg in self.args if isinstance(arg, TensExpr)])
3974 coeff = reduce(operator.mul, [a for a in self.args if not isinstance(a, TensExpr)], S.One)
3975 indices, free, free_names, dummy_data = TensMul._indices_to_free_dum(args_indices)
3976 dum = TensMul._dummy_data_to_dum(dummy_data)
3977 ext_rank = self.ext_rank
3978 free.sort(key=lambda x: x[1])
3979 free_indices = [i[0] for i in free]
3980 return free_indices, coeff*_TensorDataLazyEvaluator.data_contract_dum(arrays, dum, ext_rank)
3982 @property
3983 def data(self):
3984 deprecate_data()
3985 with ignore_warnings(SymPyDeprecationWarning):
3986 dat = _tensor_data_substitution_dict[self.expand()]
3987 return dat
3989 @data.setter
3990 def data(self, data):
3991 deprecate_data()
3992 raise ValueError("Not possible to set component data to a tensor expression")
3994 @data.deleter
3995 def data(self):
3996 deprecate_data()
3997 raise ValueError("Not possible to delete component data to a tensor expression")
3999 def __iter__(self):
4000 deprecate_data()
4001 with ignore_warnings(SymPyDeprecationWarning):
4002 if self.data is None:
4003 raise ValueError("No iteration on abstract tensors")
4004 return self.data.__iter__()
4006 @staticmethod
4007 def _dedupe_indices(new, exclude):
4008 """
4009 exclude: set
4010 new: TensExpr
4012 If ``new`` has any dummy indices that are in ``exclude``, return a version
4013 of new with those indices replaced. If no replacements are needed,
4014 return None
4016 """
4017 exclude = set(exclude)
4018 dums_new = set(get_dummy_indices(new))
4019 free_new = set(get_free_indices(new))
4021 conflicts = dums_new.intersection(exclude)
4022 if len(conflicts) == 0:
4023 return None
4025 """
4026 ``exclude_for_gen`` is to be passed to ``_IndexStructure._get_generator_for_dummy_indices()``.
4027 Since the latter does not use the index position for anything, we just
4028 set it as ``None`` here.
4029 """
4030 exclude.update(dums_new)
4031 exclude.update(free_new)
4032 exclude_for_gen = [(i, None) for i in exclude]
4033 gen = _IndexStructure._get_generator_for_dummy_indices(exclude_for_gen)
4034 repl = {}
4035 for d in conflicts:
4036 if -d in repl.keys():
4037 continue
4038 newname = gen(d.tensor_index_type)
4039 new_d = d.func(newname, *d.args[1:])
4040 repl[d] = new_d
4041 repl[-d] = -new_d
4043 if len(repl) == 0:
4044 return None
4046 new_renamed = new._replace_indices(repl)
4047 return new_renamed
4049 def _dedupe_indices_in_rule(self, rule):
4050 """
4051 rule: dict
4053 This applies TensMul._dedupe_indices on all values of rule.
4055 """
4056 index_rules = {k:v for k,v in rule.items() if isinstance(k, TensorIndex)}
4057 other_rules = {k:v for k,v in rule.items() if k not in index_rules.keys()}
4058 exclude = set(self.get_indices())
4060 newrule = {}
4061 newrule.update(index_rules)
4062 exclude.update(index_rules.keys())
4063 exclude.update(index_rules.values())
4064 for old, new in other_rules.items():
4065 new_renamed = TensMul._dedupe_indices(new, exclude)
4066 if old == new or new_renamed is None:
4067 newrule[old] = new
4068 else:
4069 newrule[old] = new_renamed
4070 exclude.update(get_indices(new_renamed))
4071 return newrule
4073 def _eval_rewrite_as_Indexed(self, *args):
4074 from sympy.concrete.summations import Sum
4075 index_symbols = [i.args[0] for i in self.get_indices()]
4076 args = [arg.args[0] if isinstance(arg, Sum) else arg for arg in args]
4077 expr = Mul.fromiter(args)
4078 return self._check_add_Sum(expr, index_symbols)
4080 def _eval_partial_derivative(self, s):
4081 # Evaluation like Mul
4082 terms = []
4083 for i, arg in enumerate(self.args):
4084 # checking whether some tensor instance is differentiated
4085 # or some other thing is necessary, but ugly
4086 if isinstance(arg, TensExpr):
4087 d = arg._eval_partial_derivative(s)
4088 else:
4089 # do not call diff is s is no symbol
4090 if s._diff_wrt:
4091 d = arg._eval_derivative(s)
4092 else:
4093 d = S.Zero
4094 if d:
4095 terms.append(TensMul.fromiter(self.args[:i] + (d,) + self.args[i + 1:]))
4096 return TensAdd.fromiter(terms)
4099class TensorElement(TensExpr):
4100 """
4101 Tensor with evaluated components.
4103 Examples
4104 ========
4106 >>> from sympy.tensor.tensor import TensorIndexType, TensorHead, TensorSymmetry
4107 >>> from sympy import symbols
4108 >>> L = TensorIndexType("L")
4109 >>> i, j, k = symbols("i j k")
4110 >>> A = TensorHead("A", [L, L], TensorSymmetry.fully_symmetric(2))
4111 >>> A(i, j).get_free_indices()
4112 [i, j]
4114 If we want to set component ``i`` to a specific value, use the
4115 ``TensorElement`` class:
4117 >>> from sympy.tensor.tensor import TensorElement
4118 >>> te = TensorElement(A(i, j), {i: 2})
4120 As index ``i`` has been accessed (``{i: 2}`` is the evaluation of its 3rd
4121 element), the free indices will only contain ``j``:
4123 >>> te.get_free_indices()
4124 [j]
4125 """
4127 def __new__(cls, expr, index_map):
4128 if not isinstance(expr, Tensor):
4129 # remap
4130 if not isinstance(expr, TensExpr):
4131 raise TypeError("%s is not a tensor expression" % expr)
4132 return expr.func(*[TensorElement(arg, index_map) for arg in expr.args])
4133 expr_free_indices = expr.get_free_indices()
4134 name_translation = {i.args[0]: i for i in expr_free_indices}
4135 index_map = {name_translation.get(index, index): value for index, value in index_map.items()}
4136 index_map = {index: value for index, value in index_map.items() if index in expr_free_indices}
4137 if len(index_map) == 0:
4138 return expr
4139 free_indices = [i for i in expr_free_indices if i not in index_map.keys()]
4140 index_map = Dict(index_map)
4141 obj = TensExpr.__new__(cls, expr, index_map)
4142 obj._free_indices = free_indices
4143 return obj
4145 @property
4146 def free(self):
4147 return [(index, i) for i, index in enumerate(self.get_free_indices())]
4149 @property
4150 def dum(self):
4151 # TODO: inherit dummies from expr
4152 return []
4154 @property
4155 def expr(self):
4156 return self._args[0]
4158 @property
4159 def index_map(self):
4160 return self._args[1]
4162 @property
4163 def coeff(self):
4164 return S.One
4166 @property
4167 def nocoeff(self):
4168 return self
4170 def get_free_indices(self):
4171 return self._free_indices
4173 def _replace_indices(self, repl: dict[TensorIndex, TensorIndex]) -> TensExpr:
4174 # TODO: can be improved:
4175 return self.xreplace(repl)
4177 def get_indices(self):
4178 return self.get_free_indices()
4180 def _extract_data(self, replacement_dict):
4181 ret_indices, array = self.expr._extract_data(replacement_dict)
4182 index_map = self.index_map
4183 slice_tuple = tuple(index_map.get(i, slice(None)) for i in ret_indices)
4184 ret_indices = [i for i in ret_indices if i not in index_map]
4185 array = array.__getitem__(slice_tuple)
4186 return ret_indices, array
4189class WildTensorHead(TensorHead):
4190 """
4191 A wild object that is used to create ``WildTensor`` instances
4193 Explanation
4194 ===========
4196 Examples
4197 ========
4198 >>> from sympy.tensor.tensor import TensorHead, TensorIndex, WildTensorHead, TensorIndexType
4199 >>> R3 = TensorIndexType('R3', dim=3)
4200 >>> p = TensorIndex('p', R3)
4201 >>> q = TensorIndex('q', R3)
4203 A WildTensorHead can be created without specifying a ``TensorIndexType``
4205 >>> W = WildTensorHead("W")
4207 Calling it with a ``TensorIndex`` creates a ``WildTensor`` instance.
4209 >>> type(W(p))
4210 <class 'sympy.tensor.tensor.WildTensor'>
4212 The ``TensorIndexType`` is automatically detected from the index that is passed
4214 >>> W(p).component
4215 W(R3)
4217 Calling it with no indices returns an object that can match tensors with any number of indices.
4219 >>> K = TensorHead('K', [R3])
4220 >>> Q = TensorHead('Q', [R3, R3])
4221 >>> W().matches(K(p))
4222 {W: K(p)}
4223 >>> W().matches(Q(p,q))
4224 {W: Q(p, q)}
4226 If you want to ignore the order of indices while matching, pass ``unordered_indices=True``.
4228 >>> U = WildTensorHead("U", unordered_indices=True)
4229 >>> W(p,q).matches(Q(q,p))
4230 >>> U(p,q).matches(Q(q,p))
4231 {U(R3,R3): _WildTensExpr(Q(q, p))}
4233 Parameters
4234 ==========
4235 name : name of the tensor
4236 unordered_indices : whether the order of the indices matters for matching
4237 (default: False)
4239 See also
4240 ========
4241 ``WildTensor``
4242 ``TensorHead``
4244 """
4245 def __new__(cls, name, index_types=None, symmetry=None, comm=0, unordered_indices=False):
4246 if isinstance(name, str):
4247 name_symbol = Symbol(name)
4248 elif isinstance(name, Symbol):
4249 name_symbol = name
4250 else:
4251 raise ValueError("invalid name")
4253 if index_types is None:
4254 index_types = []
4256 if symmetry is None:
4257 symmetry = TensorSymmetry.no_symmetry(len(index_types))
4258 else:
4259 assert symmetry.rank == len(index_types)
4261 if symmetry != TensorSymmetry.no_symmetry(len(index_types)):
4262 raise NotImplementedError("Wild matching based on symmetry is not implemented.")
4264 obj = Basic.__new__(cls, name_symbol, Tuple(*index_types), sympify(symmetry), sympify(comm), sympify(unordered_indices))
4265 obj.comm = TensorManager.comm_symbols2i(comm)
4266 obj.unordered_indices = unordered_indices
4268 return obj
4270 def __call__(self, *indices, **kwargs):
4271 tensor = WildTensor(self, indices, **kwargs)
4272 return tensor.doit()
4275class WildTensor(Tensor):
4276 """
4277 A wild object which matches ``Tensor`` instances
4279 Explanation
4280 ===========
4281 This is instantiated by attaching indices to a ``WildTensorHead`` instance.
4283 Examples
4284 ========
4285 >>> from sympy.tensor.tensor import TensorHead, TensorIndex, WildTensorHead, TensorIndexType
4286 >>> W = WildTensorHead("W")
4287 >>> R3 = TensorIndexType('R3', dim=3)
4288 >>> p = TensorIndex('p', R3)
4289 >>> q = TensorIndex('q', R3)
4290 >>> K = TensorHead('K', [R3])
4291 >>> Q = TensorHead('Q', [R3, R3])
4293 Matching also takes the indices into account
4294 >>> W(p).matches(K(p))
4295 {W(R3): _WildTensExpr(K(p))}
4296 >>> W(p).matches(K(q))
4297 >>> W(p).matches(K(-p))
4299 If you want to match objects with any number of indices, just use a ``WildTensor`` with no indices.
4300 >>> W().matches(K(p))
4301 {W: K(p)}
4302 >>> W().matches(Q(p,q))
4303 {W: Q(p, q)}
4305 See Also
4306 ========
4307 ``WildTensorHead``
4308 ``Tensor``
4310 """
4311 def __new__(cls, tensor_head, indices, **kw_args):
4312 is_canon_bp = kw_args.pop("is_canon_bp", False)
4314 if tensor_head.func == TensorHead:
4315 """
4316 If someone tried to call WildTensor by supplying a TensorHead (not a WildTensorHead), return a normal tensor instead. This is helpful when using subs on an expression to replace occurrences of a WildTensorHead with a TensorHead.
4317 """
4318 return Tensor(tensor_head, indices, is_canon_bp=is_canon_bp, **kw_args)
4319 elif tensor_head.func == _WildTensExpr:
4320 return tensor_head(*indices)
4322 indices = cls._parse_indices(tensor_head, indices)
4323 index_types = [ind.tensor_index_type for ind in indices]
4324 tensor_head = tensor_head.func(
4325 tensor_head.name,
4326 index_types,
4327 symmetry=None,
4328 comm=tensor_head.comm,
4329 unordered_indices=tensor_head.unordered_indices,
4330 )
4332 obj = Basic.__new__(cls, tensor_head, Tuple(*indices))
4333 obj.name = tensor_head.name
4334 obj._index_structure = _IndexStructure.from_indices(*indices)
4335 obj._free = obj._index_structure.free[:]
4336 obj._dum = obj._index_structure.dum[:]
4337 obj._ext_rank = obj._index_structure._ext_rank
4338 obj._coeff = S.One
4339 obj._nocoeff = obj
4340 obj._component = tensor_head
4341 obj._components = [tensor_head]
4342 if tensor_head.rank != len(indices):
4343 raise ValueError("wrong number of indices")
4344 obj.is_canon_bp = is_canon_bp
4345 obj._index_map = obj._build_index_map(indices, obj._index_structure)
4347 return obj
4350 def matches(self, expr, repl_dict=None, old=False):
4351 if not isinstance(expr, TensExpr) and expr != S(1):
4352 return None
4354 if repl_dict is None:
4355 repl_dict = {}
4356 else:
4357 repl_dict = repl_dict.copy()
4359 if len(self.indices) > 0:
4360 if not hasattr(expr, "get_free_indices"):
4361 return None
4362 expr_indices = expr.get_free_indices()
4363 if len(expr_indices) != len(self.indices):
4364 return None
4365 if self._component.unordered_indices:
4366 m = self._match_indices_ignoring_order(expr)
4367 if m is None:
4368 return None
4369 else:
4370 repl_dict.update(m)
4371 else:
4372 for i in range(len(expr_indices)):
4373 m = self.indices[i].matches(expr_indices[i])
4374 if m is None:
4375 return None
4376 else:
4377 repl_dict.update(m)
4379 repl_dict[self.component] = _WildTensExpr(expr)
4380 else:
4381 #If no indices were passed to the WildTensor, it may match tensors with any number of indices.
4382 repl_dict[self] = expr
4384 return repl_dict
4386 def _match_indices_ignoring_order(self, expr, repl_dict=None, old=False):
4387 """
4388 Helper method for matches. Checks if the indices of self and expr
4389 match disregarding index ordering.
4390 """
4391 if repl_dict is None:
4392 repl_dict = {}
4393 else:
4394 repl_dict = repl_dict.copy()
4396 def siftkey(ind):
4397 if isinstance(ind, WildTensorIndex):
4398 if ind.ignore_updown:
4399 return "wild, updown"
4400 else:
4401 return "wild"
4402 else:
4403 return "nonwild"
4405 indices_sifted = sift(self.indices, siftkey)
4407 matched_indices = []
4408 expr_indices_remaining = expr.get_indices()
4409 for ind in indices_sifted["nonwild"]:
4410 matched_this_ind = False
4411 for e_ind in expr_indices_remaining:
4412 if e_ind in matched_indices:
4413 continue
4414 m = ind.matches(e_ind)
4415 if m is not None:
4416 matched_this_ind = True
4417 repl_dict.update(m)
4418 matched_indices.append(e_ind)
4419 break
4420 if not matched_this_ind:
4421 return None
4423 expr_indices_remaining = [i for i in expr_indices_remaining if i not in matched_indices]
4424 for ind in indices_sifted["wild"]:
4425 matched_this_ind = False
4426 for e_ind in expr_indices_remaining:
4427 m = ind.matches(e_ind)
4428 if m is not None:
4429 if -ind in repl_dict.keys() and -repl_dict[-ind] != m[ind]:
4430 return None
4431 matched_this_ind = True
4432 repl_dict.update(m)
4433 matched_indices.append(e_ind)
4434 break
4435 if not matched_this_ind:
4436 return None
4438 expr_indices_remaining = [i for i in expr_indices_remaining if i not in matched_indices]
4439 for ind in indices_sifted["wild, updown"]:
4440 matched_this_ind = False
4441 for e_ind in expr_indices_remaining:
4442 m = ind.matches(e_ind)
4443 if m is not None:
4444 if -ind in repl_dict.keys() and -repl_dict[-ind] != m[ind]:
4445 return None
4446 matched_this_ind = True
4447 repl_dict.update(m)
4448 matched_indices.append(e_ind)
4449 break
4450 if not matched_this_ind:
4451 return None
4453 if len(matched_indices) < len(self.indices):
4454 return None
4455 else:
4456 return repl_dict
4458class WildTensorIndex(TensorIndex):
4459 """
4460 A wild object that matches TensorIndex instances.
4462 Examples
4463 ========
4464 >>> from sympy.tensor.tensor import TensorIndex, TensorIndexType, WildTensorIndex
4465 >>> R3 = TensorIndexType('R3', dim=3)
4466 >>> p = TensorIndex("p", R3)
4468 By default, covariant indices only match with covariant indices (and
4469 similarly for contravariant)
4471 >>> q = WildTensorIndex("q", R3)
4472 >>> (q).matches(p)
4473 {q: p}
4474 >>> (q).matches(-p)
4476 If you want matching to ignore whether the index is co/contra-variant, set
4477 ignore_updown=True
4479 >>> r = WildTensorIndex("r", R3, ignore_updown=True)
4480 >>> (r).matches(-p)
4481 {r: -p}
4482 >>> (r).matches(p)
4483 {r: p}
4485 Parameters
4486 ==========
4487 name : name of the index (string), or ``True`` if you want it to be
4488 automatically assigned
4489 tensor_index_type : ``TensorIndexType`` of the index
4490 is_up : flag for contravariant index (is_up=True by default)
4491 ignore_updown : bool, Whether this should match both co- and contra-variant
4492 indices (default:False)
4493 """
4494 def __new__(cls, name, tensor_index_type, is_up=True, ignore_updown=False):
4495 if isinstance(name, str):
4496 name_symbol = Symbol(name)
4497 elif isinstance(name, Symbol):
4498 name_symbol = name
4499 elif name is True:
4500 name = "_i{}".format(len(tensor_index_type._autogenerated))
4501 name_symbol = Symbol(name)
4502 tensor_index_type._autogenerated.append(name_symbol)
4503 else:
4504 raise ValueError("invalid name")
4506 is_up = sympify(is_up)
4507 ignore_updown = sympify(ignore_updown)
4508 return Basic.__new__(cls, name_symbol, tensor_index_type, is_up, ignore_updown)
4510 @property
4511 def ignore_updown(self):
4512 return self.args[3]
4514 def __neg__(self):
4515 t1 = WildTensorIndex(self.name, self.tensor_index_type,
4516 (not self.is_up), self.ignore_updown)
4517 return t1
4519 def matches(self, expr, repl_dict=None, old=False):
4520 if not isinstance(expr, TensorIndex):
4521 return None
4522 if self.tensor_index_type != expr.tensor_index_type:
4523 return None
4524 if not self.ignore_updown:
4525 if self.is_up != expr.is_up:
4526 return None
4528 if repl_dict is None:
4529 repl_dict = {}
4530 else:
4531 repl_dict = repl_dict.copy()
4533 repl_dict[self] = expr
4534 return repl_dict
4537class _WildTensExpr(Basic):
4538 """
4539 INTERNAL USE ONLY
4541 This is an object that helps with replacement of WildTensors in expressions.
4542 When this object is set as the tensor_head of a WildTensor, it replaces the
4543 WildTensor by a TensExpr (passed when initializing this object).
4545 Examples
4546 ========
4547 >>> from sympy.tensor.tensor import WildTensorHead, TensorIndex, TensorHead, TensorIndexType
4548 >>> W = WildTensorHead("W")
4549 >>> R3 = TensorIndexType('R3', dim=3)
4550 >>> p = TensorIndex('p', R3)
4551 >>> q = TensorIndex('q', R3)
4552 >>> K = TensorHead('K', [R3])
4553 >>> print( ( K(p) ).replace( W(p), W(q)*W(-q)*W(p) ) )
4554 K(R_0)*K(-R_0)*K(p)
4556 """
4557 def __init__(self, expr):
4558 if not isinstance(expr, TensExpr):
4559 raise TypeError("_WildTensExpr expects a TensExpr as argument")
4560 self.expr = expr
4562 def __call__(self, *indices):
4563 return self.expr._replace_indices(dict(zip(self.expr.get_free_indices(), indices)))
4565 def __neg__(self):
4566 return self.func(self.expr*S.NegativeOne)
4568 def __abs__(self):
4569 raise NotImplementedError
4571 def __add__(self, other):
4572 if other.func != self.func:
4573 raise TypeError(f"Cannot add {self.func} to {other.func}")
4574 return self.func(self.expr+other.expr)
4576 def __radd__(self, other):
4577 if other.func != self.func:
4578 raise TypeError(f"Cannot add {self.func} to {other.func}")
4579 return self.func(other.expr+self.expr)
4581 def __sub__(self, other):
4582 return self + (-other)
4584 def __rsub__(self, other):
4585 return other + (-self)
4587 def __mul__(self, other):
4588 raise NotImplementedError
4590 def __rmul__(self, other):
4591 raise NotImplementedError
4593 def __truediv__(self, other):
4594 raise NotImplementedError
4596 def __rtruediv__(self, other):
4597 raise NotImplementedError
4599 def __pow__(self, other):
4600 raise NotImplementedError
4602 def __rpow__(self, other):
4603 raise NotImplementedError
4606def canon_bp(p):
4607 """
4608 Butler-Portugal canonicalization. See ``tensor_can.py`` from the
4609 combinatorics module for the details.
4610 """
4611 if isinstance(p, TensExpr):
4612 return p.canon_bp()
4613 return p
4616def tensor_mul(*a):
4617 """
4618 product of tensors
4619 """
4620 if not a:
4621 return TensMul.from_data(S.One, [], [], [])
4622 t = a[0]
4623 for tx in a[1:]:
4624 t = t*tx
4625 return t
4628def riemann_cyclic_replace(t_r):
4629 """
4630 replace Riemann tensor with an equivalent expression
4632 ``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)``
4634 """
4635 free = sorted(t_r.free, key=lambda x: x[1])
4636 m, n, p, q = [x[0] for x in free]
4637 t0 = t_r*Rational(2, 3)
4638 t1 = -t_r.substitute_indices((m,m),(n,q),(p,n),(q,p))*Rational(1, 3)
4639 t2 = t_r.substitute_indices((m,m),(n,p),(p,n),(q,q))*Rational(1, 3)
4640 t3 = t0 + t1 + t2
4641 return t3
4643def riemann_cyclic(t2):
4644 """
4645 Replace each Riemann tensor with an equivalent expression
4646 satisfying the cyclic identity.
4648 This trick is discussed in the reference guide to Cadabra.
4650 Examples
4651 ========
4653 >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorHead, riemann_cyclic, TensorSymmetry
4654 >>> Lorentz = TensorIndexType('Lorentz', dummy_name='L')
4655 >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz)
4656 >>> R = TensorHead('R', [Lorentz]*4, TensorSymmetry.riemann())
4657 >>> t = R(i,j,k,l)*(R(-i,-j,-k,-l) - 2*R(-i,-k,-j,-l))
4658 >>> riemann_cyclic(t)
4659 0
4660 """
4661 t2 = t2.expand()
4662 if isinstance(t2, (TensMul, Tensor)):
4663 args = [t2]
4664 else:
4665 args = t2.args
4666 a1 = [x.split() for x in args]
4667 a2 = [[riemann_cyclic_replace(tx) for tx in y] for y in a1]
4668 a3 = [tensor_mul(*v) for v in a2]
4669 t3 = TensAdd(*a3).doit()
4670 if not t3:
4671 return t3
4672 else:
4673 return canon_bp(t3)
4676def get_lines(ex, index_type):
4677 """
4678 Returns ``(lines, traces, rest)`` for an index type,
4679 where ``lines`` is the list of list of positions of a matrix line,
4680 ``traces`` is the list of list of traced matrix lines,
4681 ``rest`` is the rest of the elements of the tensor.
4682 """
4683 def _join_lines(a):
4684 i = 0
4685 while i < len(a):
4686 x = a[i]
4687 xend = x[-1]
4688 xstart = x[0]
4689 hit = True
4690 while hit:
4691 hit = False
4692 for j in range(i + 1, len(a)):
4693 if j >= len(a):
4694 break
4695 if a[j][0] == xend:
4696 hit = True
4697 x.extend(a[j][1:])
4698 xend = x[-1]
4699 a.pop(j)
4700 continue
4701 if a[j][0] == xstart:
4702 hit = True
4703 a[i] = reversed(a[j][1:]) + x
4704 x = a[i]
4705 xstart = a[i][0]
4706 a.pop(j)
4707 continue
4708 if a[j][-1] == xend:
4709 hit = True
4710 x.extend(reversed(a[j][:-1]))
4711 xend = x[-1]
4712 a.pop(j)
4713 continue
4714 if a[j][-1] == xstart:
4715 hit = True
4716 a[i] = a[j][:-1] + x
4717 x = a[i]
4718 xstart = x[0]
4719 a.pop(j)
4720 continue
4721 i += 1
4722 return a
4724 arguments = ex.args
4725 dt = {}
4726 for c in ex.args:
4727 if not isinstance(c, TensExpr):
4728 continue
4729 if c in dt:
4730 continue
4731 index_types = c.index_types
4732 a = []
4733 for i in range(len(index_types)):
4734 if index_types[i] is index_type:
4735 a.append(i)
4736 if len(a) > 2:
4737 raise ValueError('at most two indices of type %s allowed' % index_type)
4738 if len(a) == 2:
4739 dt[c] = a
4740 #dum = ex.dum
4741 lines = []
4742 traces = []
4743 traces1 = []
4744 #indices_to_args_pos = ex._get_indices_to_args_pos()
4745 # TODO: add a dum_to_components_map ?
4746 for p0, p1, c0, c1 in ex.dum_in_args:
4747 if arguments[c0] not in dt:
4748 continue
4749 if c0 == c1:
4750 traces.append([c0])
4751 continue
4752 ta0 = dt[arguments[c0]]
4753 ta1 = dt[arguments[c1]]
4754 if p0 not in ta0:
4755 continue
4756 if ta0.index(p0) == ta1.index(p1):
4757 # case gamma(i,s0,-s1) in c0, gamma(j,-s0,s2) in c1;
4758 # to deal with this case one could add to the position
4759 # a flag for transposition;
4760 # one could write [(c0, False), (c1, True)]
4761 raise NotImplementedError
4762 # if p0 == ta0[1] then G in pos c0 is mult on the right by G in c1
4763 # if p0 == ta0[0] then G in pos c1 is mult on the right by G in c0
4764 ta0 = dt[arguments[c0]]
4765 b0, b1 = (c0, c1) if p0 == ta0[1] else (c1, c0)
4766 lines1 = lines[:]
4767 for line in lines:
4768 if line[-1] == b0:
4769 if line[0] == b1:
4770 n = line.index(min(line))
4771 traces1.append(line)
4772 traces.append(line[n:] + line[:n])
4773 else:
4774 line.append(b1)
4775 break
4776 elif line[0] == b1:
4777 line.insert(0, b0)
4778 break
4779 else:
4780 lines1.append([b0, b1])
4782 lines = [x for x in lines1 if x not in traces1]
4783 lines = _join_lines(lines)
4784 rest = []
4785 for line in lines:
4786 for y in line:
4787 rest.append(y)
4788 for line in traces:
4789 for y in line:
4790 rest.append(y)
4791 rest = [x for x in range(len(arguments)) if x not in rest]
4793 return lines, traces, rest
4796def get_free_indices(t):
4797 if not isinstance(t, TensExpr):
4798 return ()
4799 return t.get_free_indices()
4802def get_indices(t):
4803 if not isinstance(t, TensExpr):
4804 return ()
4805 return t.get_indices()
4807def get_dummy_indices(t):
4808 if not isinstance(t, TensExpr):
4809 return ()
4810 inds = t.get_indices()
4811 free = t.get_free_indices()
4812 return [i for i in inds if i not in free]
4814def get_index_structure(t):
4815 if isinstance(t, TensExpr):
4816 return t._index_structure
4817 return _IndexStructure([], [], [], [])
4820def get_coeff(t):
4821 if isinstance(t, Tensor):
4822 return S.One
4823 if isinstance(t, TensMul):
4824 return t.coeff
4825 if isinstance(t, TensExpr):
4826 raise ValueError("no coefficient associated to this tensor expression")
4827 return t
4829def contract_metric(t, g):
4830 if isinstance(t, TensExpr):
4831 return t.contract_metric(g)
4832 return t
4835def perm2tensor(t, g, is_canon_bp=False):
4836 """
4837 Returns the tensor corresponding to the permutation ``g``
4839 For further details, see the method in ``TIDS`` with the same name.
4840 """
4841 if not isinstance(t, TensExpr):
4842 return t
4843 elif isinstance(t, (Tensor, TensMul)):
4844 nim = get_index_structure(t).perm2tensor(g, is_canon_bp=is_canon_bp)
4845 res = t._set_new_index_structure(nim, is_canon_bp=is_canon_bp)
4846 if g[-1] != len(g) - 1:
4847 return -res
4849 return res
4850 raise NotImplementedError()
4853def substitute_indices(t, *index_tuples):
4854 if not isinstance(t, TensExpr):
4855 return t
4856 return t.substitute_indices(*index_tuples)
4859def _expand(expr, **kwargs):
4860 if isinstance(expr, TensExpr):
4861 return expr._expand(**kwargs)
4862 else:
4863 return expr.expand(**kwargs)