[sympy] Tensor expression with canonicalization: part1 (#1700)
pernici
It is implemented a tensor class with monoterm canonicalization
and algebraic operations on tensors.
As an application,
check of identities involving the cyclic symmetry of the Riemann tensor
```
>>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorS
>>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
>>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz)
>>> symr = TensorSymmetry(riemann_bsgs)
>>> R4 = TensorType([Lorentz]*4, symr)
>>> R = R4('R')
>>> t = R(i,j,k,l)*(R(-i,-j,-k,-l) - 2*R(-i,-k,-j,-l))
>>> riemann_cyclic(t)
0
```
You can merge this Pull Request by running:
git pull https://github.com/pernici/sympy tensor1
Or view, comment on, or merge it at:
https://github.com/sympy/sympy/pull/1700
Commit Summary
- pep8 tensor_can
- reformat for sphinx
- typo edit
- Fixed a bug in canonicalize; corrected comments in test_tensor_can.py
- Tensor expressions with canonicalization
File Changes
- A doc/src/modules/tensor/tensor.rst (41)
- M sympy/combinatorics/tensor_can.py (209)
- M sympy/combinatorics/tests/test_tensor_can.py (14)
- M sympy/core/basic.py (1)
- M sympy/printing/str.py (6)
- A sympy/tensor/tensor.py (1260)
- A sympy/tensor/tests/test_tensor.py (665)
Patch Links
—
Reply to this email directly or view it on GitHub.
Julien Rioux
SymPy Bot Summary: 
Failed after merging pernici/tensor1 (38e15aa) into master (d503614).
@pernici: Please fix the test failures.
PyPy 2.0.0-beta-1; 2.7.3-final-42: fail
Python 2.7.2-final-0: fail
Python 3.2.1-final-0: fail
Sphinx 1.1.3: fail
Docs build command: make clean && make html-errors && make latex && cd _build/latex && xelatex sympy-*.tex
Matthew Rocklin
In doc/src/modules/tensor/tensor.rst:
> @@ -0,0 +1,41 @@ > +.. _tensor-tensor:
I love that you included a docs page with the first PR. We should encourage this more in the future.
Matthew Rocklin
In sympy/tensor/tensor.py:
> @@ -0,0 +1,1277 @@ > +from collections import defaultdict > +from sympy.core import Basic, sympify, Add, Mul, S > +from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, bsgs_direct_product, canonicalize, riemann_bsgs > + > +class TensorIndexType(object):
You have a few objects that don't inherit from Basic. If these contain mathematical information about an expression then they should be Basic so that other sympy utilities can manipulate them (through subs, replace, free_symbols, rules, etc...)
Matthew Rocklin
In sympy/tensor/tensor.py:
> +class TensorIndexType(object): > + """ > + A TensorIndexType is characterized by its name, by metric_sym, > + giving the symmetry of its metric. > + > + If a dimension ``dim`` is defined, it can be a symbol or an integer. > + """ > + def __init__(self, name, metric_sym=0, dim=None, eps_dim = None, > + dummy_fmt=None): > + """ > + name name of the tensor type > + > + ``metric_sym``: > + 0 symmetric > + 1 antisymmetric > + None no symmetry
Is there a mathematical backing behind this mapping? If not how about using strings or constants. For example why not
TensorIndexType("some name", "symmetric" ... )
rather than
TensorIndexType("some name", 1 ... )
Matthew Rocklin
In sympy/tensor/tensor.py:
> + ========
> +
> + >>> from sympy.tensor.tensor import TensorIndexType, TensorIndex, TensorSymmetry, TensorType, get_symmetric_group_sgs
> + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
> + >>> i = TensorIndex('i', Lorentz); i
> + i
> + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1))
> + >>> S1 = TensorType([Lorentz], sym1)
> + >>> A, B = S1('A,B')
> + >>> A(i)*B(-i)
> + A(L_0)*B(-L_0)
> + """
> + def __init__(self, name, tensortype, is_contravariant=True):
> + self.name = name
> + self.tensortype = tensortype
> + self.is_contravariant = is_contravariant
Why is contravariant vs is covariant?
Matthew Rocklin
In sympy/tensor/tensor.py:
> + ``s`` string of comma separated names of indices
> +
> + ``typ`` list of TensorIndexType of the indices
> +
> + Examples
> + ========
> +
> + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices
> + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
> + >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz)
> + """
> + a = s.split(',')
> + return [TensorIndex(i, typ) for i in a]
> +
> +
> +class TensorSymmetry(object):
This is a very small class. Does it need to be a class?
Matthew Rocklin
In sympy/tensor/tensor.py:
> + >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz)
> + >>> sym2 = TensorSymmetry(get_symmetric_group_sgs(2))
> + >>> S2 = TensorType([Lorentz]*2, sym2)
> + >>> V = S2('V')
> + """
> + is_commutative = False
> +
> + def __new__(cls, index_types, symmetry, **kw_args):
> + obj = Basic.__new__(cls, **kw_args)
> + obj.index_types = []
> + obj.index_types = index_types
> + obj.types = list(set(obj.index_types))
> + obj.types.sort(key=lambda x: x.name)
> + obj.symmetry = symmetry
> + assert symmetry.rank == len(index_types) + 2
> + return obj
This doesn't follow standard use of SymPy args.
Each sympy Basic object is a node in a tree. That node has a type (TensorType here) and a bunch of children which we call args. Mathematical properties are then determined through computations on the args. In this style your TensorType class would be written something like this (skipping kw_args out of laziness).
class TensorType(Basic): def __new__(cls, index_types, symmetry): assert symmetry.rank = len(index_types) + 2 return Basic.__new__(cls, index_types, symmetry) index_types = property(lambda self: self.args[0]) symmetry = property(lambda self: self.args[1]) types = property(lambda self: sorted(set(self.index_types), key=lambda x: x.name))
I'm sure I've made some errors above but I hope that it conveys the general message.
By doing this general SymPy utilities can more effectively inter-operate with your code in the future.
Matthew Rocklin
In sympy/tensor/tensor.py:
> + obj.index_types = index_types > + obj.types = list(set(obj.index_types)) > + obj.types.sort(key=lambda x: x.name
) > + obj.symmetry = symmetry > + assert symmetry.rank == len(index_types) + 2 > + return obj > + > + def __str__(self): > + return 'TensorType(%s)' %([str(x) for x in self.index_types]) > + > + def __call__(self, s, commuting=0): > + """ > + commuting: > + None no commutation rule > + 0 commutes > + 1 anticommutes
I would choose a representation which more closely matched mathematics. I would use "commutes" rather than 0.
Matthew Rocklin
In sympy/tensor/tensor.py:
> + >>> sym2 = TensorSymmetry(get_symmetric_group_sgs(2))
> + >>> S2 = TensorType([Lorentz]*2, sym2)
> + >>> A = S2('A')
> + """
> + assert isinstance(name, basestring)
> +
> + obj = Basic.__new__(cls, **kw_args)
> + obj.name = name
> + obj.index_types = typ.index_types
> + obj.rank = len(obj.index_types)
> + obj.types = typ.types
> + obj.symmetry = typ.symmetry
> + obj.commuting = commuting
> + return obj
> +
> + def __eq__(self, other):
You shouldn't have to implement __eq__. Basic should take care of it for you if you follow the standard args style I mentioned above. We actually prefer that subclasses of Basic don't implement __eq__
Matthew Rocklin
In sympy/tensor/tensor.py:
> + obj.index_types = index_types > + obj.types = list(set(obj.index_types)) > + obj.types.sort(key=lambda x: x.name
) > + obj.symmetry = symmetry > + assert symmetry.rank == len(index_types) + 2 > + return obj > + > + def __str__(self): > + return 'TensorType(%s)' %([str(x) for x in self.index_types]) > + > + def __call__(self, s, commuting=0): > + """ > + commuting: > + None no commutation rule > + 0 commutes > + 1 anticommutes
Or how about commutes = True/False/None this follows SymPy's standard three value logic.
Matthew Rocklin
In sympy/tensor/tensor.py:
> + if self.is_TensAdd: > + args = self.args + (other,) > + return TensAdd(*args) > + return TensAdd(self, other) > + > + def __radd__(self, other): > + return TensAdd(other, self) > + > + def __sub__(self, other): > + other = sympify(other) > + return TensAdd(self, -other) > + > + def __rsub__(self, other): > + return TensAdd(other, -self) > + > + def __mul__(self, other):
Is there a way to have a tensor expression that isn't entirely distributed?
Matthew Rocklin
In sympy/tensor/tensor.py:
> + >>> a, b = tensor_indices('a,b', Lorentz)
> + >>> sym = TensorSymmetry(get_symmetric_group_sgs(1))
> + >>> S1 = TensorType([Lorentz], sym)
> + >>> p, q = S1('p,q')
> + >>> t1 = p(a)
> + >>> t2 = q(a)
> + >>> t1 + t2
> + p(a) + q(a)
> + """
> + is_Tensor = True
> + is_TensAdd = True
> +
> + def __new__(cls, *args, **kw_args):
> + """
> + args tuple of addends
> + """
This function is huge and scares me. Can it be factored down or simplified at all?
Matthew Rocklin
In sympy/tensor/tensor.py:
> + t = TensAdd(*args)
> + return t.canon_bp()
> +
> + def _pretty(self):
> + a = []
> + args = self.args
> + for x in args:
> + a.append(str(x))
> + a.sort()
> + s = ' + '.join(a)
> + s = s.replace('+ -', '- ')
> + return s
> +
> +class Tensor(TensExpr):
> + """
> + Product of tensors
Should the name be TensorProduct or something similar? Why does this class get the simplest name?
Matthew Rocklin
I'm unable to competently review the mathematics but I have left some comments on general organization. My thoughts are as follows.
Its clear that this PR contains a substantial amount of mathematical logic that is highly appropriate for SymPy (we want this). However I think that the current implementation is isolated. It will work well with itself but will not extend well to future work within the SymPy ecosystem. I've made a few direct comments on specific issues (e.g. use args and Basic so that we have a well formed expression tree.)
In general I would also like to suggest that as much of the mathematical knowledge you've encoded here be as far removed from the expression tree as possible. If you can pull logic from a method out to a function I suggest doing so (although this is only my preference, not necessarily SymPy wisdom.) For example I really like that the Butler-Portugal algorithm is separate from this implementation of tensor expressions. This ensures that it can be easily used in the future should another tensor implementation emerge (a possibility we need to plan for given this topic's history in SymPy). You have some large methods in your classes, I suspect that these encode a fair amount of math. Is it possible to pull out large chunks of them to general functions as is done with BP?
I suspect that @Krastanov will be more able to judge the mathematical content than I.
Aaron Meurer
In sympy/tensor/tensor.py:
> + obj.index_types = index_types > + obj.types = list(set(obj.index_types)) > + obj.types.sort(key=lambda x: x.name) > + obj.symmetry = symmetry > + assert symmetry.rank == len(index_types) + 2 > + return obj > + > + def __str__(self): > + return 'TensorType(%s)' %([str(x) for x in self.index_types]) > + > + def __call__(self, s, commuting=0): > + """ > + commuting: > + None no commutation rule > + 0 commutes > + 1 anticommutes
The problem is that "anticommutes" doesn't really mean "not commutes", so commutes=False for it is a little confusing for it (the former implies the latter, but not the reverse). A way would be to just have both "commutes" and "anticommutes", and have it be an error if they are both true (or it implies that it is zero or something). Probably the assumptions system could be leveraged somehow to do this automatically. Ditto for symmetric/antisymmetric.
Aaron Meurer
I also like having the algorithms separate from the user-level classes. This is how the polys work. There, there are actually several layers of abstraction. Aside from maintaining good abstraction, which is always a goal, it allows to keep the user-level code high-level (readable and usable via things like syntactic sugar and metaclassery), but the core algorithm fast and readable. I use readable in both places because they apply to different things. For user-level code (or more generally, the higher up in the abstraction layers you are), the actual use of the code should be readable. For low-level algorithmic code, the implementation itself should be readable, because it is the system critical part.
pernici
In sympy/tensor/tensor.py:
> + ``s`` string of comma separated names of indices
> +
> + ``typ`` list of TensorIndexType of the indices
> +
> + Examples
> + ========
> +
> + >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices
> + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
> + >>> a, b, c, d = tensor_indices('a,b,c,d', Lorentz)
> + """
> + a = s.split(',')
> + return [TensorIndex(i, typ) for i in a]
> +
> +
> +class TensorSymmetry(object):
The tensor symmetry might be specified using Young tableaux instead of generators; this class should deal with that
in the future.
pernici
In sympy/tensor/tensor.py:
> + >>> sym2 = TensorSymmetry(get_symmetric_group_sgs(2))
> + >>> S2 = TensorType([Lorentz]*2, sym2)
> + >>> A = S2('A')
> + """
> + assert isinstance(name, basestring)
> +
> + obj = Basic.__new__(cls, **kw_args)
> + obj.name= name > + obj.index_types = typ.index_types > + obj.rank = len(obj.index_types) > + obj.types = typ.types > + obj.symmetry = typ.symmetry > + obj.commuting = commuting > + return obj > + > + def __eq__(self, other):
Nice; I verified that doing this __eq__ is not necessary anymore in TensorIndex and TensorHead
pernici
In sympy/tensor/tensor.py:
> + obj.index_types = index_types > + obj.types = list(set(obj.index_types)) > + obj.types.sort(key=lambda x: x.name
) > + obj.symmetry = symmetry > + assert symmetry.rank == len(index_types) + 2 > + return obj > + > + def __str__(self): > + return 'TensorType(%s)' %([str(x) for x in self.index_types]) > + > + def __call__(self, s, commuting=0): > + """ > + commuting: > + None no commutation rule > + 0 commutes > + 1 anticommutes
For the moment I will not change this; the rules for commutation should be more detailed; for instance two tensors
can be non commuting but commuting one with the other, so commuting=None is not satisfactory in this case.
For instance Cadabra has this feature.
pernici
In sympy/tensor/tensor.py:
> + if self.is_TensAdd: > + args = self.args + (other,) > + return TensAdd(*args) > + return TensAdd(self, other) > + > + def __radd__(self, other): > + return TensAdd(other, self) > + > + def __sub__(self, other): > + other = sympify(other) > + return TensAdd(self, -other) > + > + def __rsub__(self, other): > + return TensAdd(other, -self) > + > + def __mul__(self, other):
It would be nice, but I do not know enough the internals of SymPy to attempt to do it.
pernici
In sympy/tensor/tensor.py:
> + >>> a, b = tensor_indices('a,b', Lorentz)
> + >>> sym = TensorSymmetry(get_symmetric_group_sgs(1))
> + >>> S1 = TensorType([Lorentz], sym)
> + >>> p, q = S1('p,q')
> + >>> t1 = p(a)
> + >>> t2 = q(a)
> + >>> t1 + t2
> + p(a) + q(a)
> + """
> + is_Tensor = True
> + is_TensAdd = True
> +
> + def __new__(cls, *args, **kw_args):
> + """
> + args tuple of addends
> + """
If factored two functions out of it; it will be in next commit.
pernici
However, I think that this class should be renamed to something like "Manifold/TangentSpace" for the reasons that I already stated.
The indices can be formal objects, like in the case of the Lorentz index in dimensional regularization, so
I think that TensorIndexType is fine.
pernici
Is TensorSymmetry used for anything else besides defining a TensorType?
No, but as I wrote in another comment, TensorSymmetry should in the future allow to specify the symmetry of the
tensor using Young tableaux instead of writing directly the BSGS.
pernici
In sympy/tensor/tensor.py:
> +class TensorIndexType(object): > + """ > + A TensorIndexType is characterized by its name, by metric_sym, > + giving the symmetry of its metric. > + > + If a dimension ``dim`` is defined, it can be a symbol or an integer. > + """ > + def __init__(self, name, metric_sym=0, dim=None, eps_dim = None, > + dummy_fmt=None): > + """ > + name name of the tensor type > + > + ``metric_sym``: > + 0 symmetric > + 1 antisymmetric > + None no symmetry
I am changing metric_sym to metric_antisym; if it is True the metric is antisymmetric, if False it is symmetric,
if None there is no metric.
Similarly commuting to anticommuting in TensorType.__call__, sym to antisym in get_symmetric_group_sgs
pernici
In sympy/tensor/tensor.py:
> + ========
> +
> + >>> from sympy.tensor.tensor import TensorIndexType, TensorIndex, TensorSymmetry, TensorType, get_symmetric_group_sgs
> + >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
> + >>> i = TensorIndex('i', Lorentz); i
> + i
> + >>> sym1 = TensorSymmetry(get_symmetric_group_sgs(1))
> + >>> S1 = TensorType([Lorentz], sym1)
> + >>> A, B = S1('A,B')
> + >>> A(i)*B(-i)
> + A(L_0)*B(-L_0)
> + """
> + def __init__(self, name, tensortype, is_contravariant=True):
> + self.name = name
> + self.tensortype = tensortype
> + self.is_contravariant = is_contravariant
With the representation i for contravariant index and -i for covariant index, the contravariant index "comes first" so is_contravariant seems to me more natural.
pernici
In sympy/tensor/tensor.py:
> + t = TensAdd(*args)
> + return t.canon_bp()
> +
> + def _pretty(self):
> + a = []
> + args = self.args
> + for x in args:
> + a.append(str(x))
> + a.sort()
> + s = ' + '.join(a)
> + s = s.replace('+ -', '- ')
> + return s
> +
> +class Tensor(TensExpr):
> + """
> + Product of tensors
I have replaced Tensor with TensMul
pernici
I think that it makes sense to merge TensorHead and TensorType. TensorType seems to be used as a shortcut for creating similar tensors (like symbols, vars, numbered_symbols within the rest of sympy).
A TensorType creates objects with given tensor symmetry, e.g. a symmetric tensor of rank 2; from it one
can generate several TensorHead with that symmetry; therefore it seems to me convenient to keep both.
pernici
An important detail is that we have no way to account for the difference between a TangentSpace and its dual, but we do not care because this is a metric manifold and we have TheCanonicalWayToSpecifyADual.
In the cases without metric, like in the case of the defining representation of SU(N) for N > 2 there is
no metric converting it to the conjugate representation; the defining representation is represented by a
covariant vector, the conjugate representation by a contravariant vector.
In this PR they have the same TensorType S1 = TensorType([Color], sym1);q, qbar = S1('q, qbar')
One cannot raise the index of.Group Theory` by Cvitanovic, which has a tensorial approach
I have started to look at
to group theory, to see how well this PR works in that case.
Stefan Krastanov
However, I think that this class should be renamed to something like "Manifold/TangentSpace" for the reasons that I already stated.
The indices can be formal objects, like in the case of the Lorentz index in dimensional regularization, so
I think that TensorIndexType is fine.
This is a neat detail. Maybe it can be documented, but anyway, the current name is ok.
Stefan Krastanov
An important detail is that we have no way to account for the difference between a TangentSpace and its dual, but we do not care because this is a metric manifold and we have TheCanonicalWayToSpecifyADual.
In the cases without metric, like in the case of the defining representation of SU(N) for N > 2 there is
no metric converting it to the conjugate representation; the defining representation is represented by a
covariant vector, the conjugate representation by a contravariant vector.
I do not see a problem...
Maybe I am just not understanding the issue, but I think that there is a metric in this case. The usual representation of SU(N) is a Hilbert space which by definition is a vector space with a scalar product. This scalar product plays the role of a metric.
And more generally, I was left with the impression that for any vector space (without a scalar product structure) there is always a canonical way to define a scalar product and hence get a metric...
Stefan Krastanov
I think that it makes sense to merge TensorHead and TensorType. TensorType seems to be used as a shortcut for creating similar tensors (like symbols, vars, numbered_symbols within the rest of sympy).
A TensorType creates objects with given tensor symmetry, e.g. a symmetric tensor of rank 2; from it one
can generate several TensorHead with that symmetry; therefore it seems to me convenient to keep both.
It would be convenient to have a utility function that produces TensorHeads with the same symmetry, however I do not think it makes sense for it to be:
- a whole class (As far as I understand all the important stuff is done within the
__call__method. Why not just use a function) - a class that subclasses Basic (You do not create expression trees containing
TensorType, nor compareTensorType. Basically,TensorTypeis not something that is symbolically manipulated (it is used to create an object that can be symbolically manipulted))
Stefan Krastanov
Summary: These two questions seem important to me:
- pernici@38e15aa#sympy-tensor-tensor-py-P60 (having a special class for the metric)
-
#1700 (
TensorType)
Stefan Krastanov
Summary (Addendum): These two questions seem important to me:
- https://github.com/pernici/sympy/commit/38e15aa3a510ca9b1fb389361e86e91121b2e096#sympy-tensor-tensor-py-P60 (having a special class for the metric)
- refactoring `TensorType`
- https://github.com/pernici/sympy/commit/38e15aa3a510ca9b1fb389361e86e91121b2e096#commitcomment-2351505 (The use of `is_Tensor`)
Aaron Meurer
SymPy Bot Summary: Failed after merging pernici/tensor1 (9a52a76) into master (2abef4f).
@pernici: Please fix the test failures.
Python 2.5.0-final-0: pass
Python 2.6.6-final-0: fail
Python 2.7.2-final-0: pass
Python 2.6.8-final-0: pass
Python 2.7.3-final-0: fail
PyPy 2.0.0-beta-1; 2.7.3-final-42: pass
Python 3.2.2-final-0: pass
Python 3.3.0-final-0: pass
Python 3.2.3-final-0: pass
Python 3.3.0-final-0: pass
Python 3.3.0-final-0: pass
**Sphinx 1.1.3:** fail
Aaron Meurer
And more generally, I was left with the impression that for any vector space (without a scalar product structure) there is always a canonical way to define a scalar product and hence get a metric...
I've never heard of this. I know that you can do it for a normed space iff the norm satisfies the parallelogram identity.
pernici
Maybe I am just not understanding the issue, but I think that there is a metric in this case. The usual representation of SU(N) is a Hilbert space which by definition is a vector space with a scalar product. This scalar product plays the role of a metric.
In SU(N) for N > 2 qbar(i) cannot be equal to g(i, j)*q(-j) for some metric g(i, j) because
the conjugate defining representation is inequivalent to the defining representation, so there cannot be a
linear relation among these two representations.
The scalar product in a Hilbert space is sesquilinear, that is it involves complex conjugation, not only linear
relations.
pernici
Summary (Addendum): These three questions seem important to me:
having a special class for the metric
refactoring TensorType
The use of is_Tensor
The third point is done; I made a change in TensorIndexType for the first point; I do not know if it is what you wanted.
I have still to look at the second point.
Do you have an idea how to implement tensor fields, derivatives, covariant derivatives?