Status of Tensor module?

[Jamie Titus]

Hello all,

I was wondering if anybody knew what the status of the Tensor module is - i.e., is work being done on it, what is the known functionality, and what is the intended functionality. For the past week I have been attempting to use Sympy for basic General Relativity (homework problems really), and it has not gone well. My naive review of the documentation suggests that I should be able to use Einstein Notation, since I am completely unfamiliar with Penrose's abstract index notation (and it doesn't seem to be as useful computationally, frankly, though the wikipedia article on it is very short). However, I keep encountering errors when I try to do basic operations such as add or multiply tensors, and I am never quite sure if they are bugs or features. A minimal example is that when you attempt to symmetrize a tensor you just get back double the tensor components, i.e., given a TensorHead A, and indices i0, i1, A(i0, i1) + A(i1, i0) should be manifestly symmetric, but it just results in 2*A(i0, i1). If necessary I can provide minimal working examples in code. I wasn't sure if this is the right place for that.

If there is somebody who is acting as a contact for the Tensor module, I'd very much like to know whether I'm encountering bugs or 'features.' If the former, I'm willing to try to submit patches to fix some of the problems that keep cropping up, or at least document them. I honestly found the Tensor module very confusing from the standpoint of somebody who is interested in solving simple GR problems, so I'd like to try to understand the rationale behind it and make things clearer for other users.

[Ondřej Čertík]

Hi Jamie,

On Fri, Oct 31, 2014 at 10:53 AM, Jamie Titus <sen...@gmail.com> wrote:
> Hello all,
>
> I was wondering if anybody knew what the status of the Tensor module is -
> i.e., is work being done on it, what is the known functionality, and what is
> the intended functionality. For the past week I have been attempting to use
> Sympy for basic General Relativity (homework problems really), and it has
> not gone well. My naive review of the documentation suggests that I should
> be able to use Einstein Notation, since I am completely unfamiliar with
> Penrose's abstract index notation (and it doesn't seem to be as useful
> computationally, frankly, though the wikipedia article on it is very short).
> However, I keep encountering errors when I try to do basic operations such
> as add or multiply tensors, and I am never quite sure if they are bugs or
> features. A minimal example is that when you attempt to symmetrize a tensor
> you just get back double the tensor components, i.e., given a TensorHead A,
> and indices i0, i1, A(i0, i1) + A(i1, i0) should be manifestly symmetric,
> but it just results in 2*A(i0, i1). If necessary I can provide minimal
> working examples in code. I wasn't sure if this is the right place for that.

Thanks for your feedback, yes, this is the right place to ask. Would
you mind sending
few snippets that you found confusing?

>
> If there is somebody who is acting as a contact for the Tensor module, I'd
> very much like to know whether I'm encountering bugs or 'features.' If the
> former, I'm willing to try to submit patches to fix some of the problems
> that keep cropping up, or at least document them. I honestly found the
> Tensor module very confusing from the standpoint of somebody who is
> interested in solving simple GR problems, so I'd like to try to understand
> the rationale behind it and make things clearer for other users.

GR was one of the main applications that I wanted to use SymPy for
when I started it. However, it was a lot of work to get the basic
functionality working, so I never got to actually write the tensor
module. Some people who worked on the tensor module are:

Francesco Bonazzi
Björn Dahlgren
Mario Pernici
...

Essentially if you go into the file and do "git blame", it tells you
who wrote it.

I would be very much interested in fixing it, so that it is easy to
use for users. If you could send us all the things that surprised you,
that would be a start. Then we can figure out how to best fix that.

Many thanks,
Ondrej

[Jamie Titus]

I think the most instructive thing that I found confusing is was my attempt to symmetrize a tensor. This minimal snippet, taken almost entirely from the examples in documentation, should demonstrate my point.

from sympy import *
from sympy.tensor.tensor import TensorIndexType, tensorsymmetry, TensorType

Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
sym2 = tensorsymmetry([1]*2)
S2 = TensorType([Lorentz]*2, sym2)
A = S2('A')
from sympy.tensor.tensor import tensor_indices, tensorhead
Lorentz.data = [1, -1, -1, -1]
i0, i1 = tensor_indices('i0:2', Lorentz)
A.data = [[j+2*i for j in range(2)] for i in range(2)]
A.data
A_symm = A(i0, i1) + A(i1, i0)
A_symm.data

A.data is [[0, 1], [2, 3]], One would want A_symm.data to be [[0, 3],[3, 6]], a symmetric tensor. However, Sympy is computing it as 2*A. Indeed, that is exactly how sympy represents A_symm.

I start off with this example because I am genuinely not sure if this is the intended behaviour, since it is unclear to me whether or not the Tensor module is supposed to support typical Einstein notation. If this is a "bug," then there are many things which must be tested and fixed. If this is a "feature," then I probably need to look for a different module to do GR with.

[Francesco Bonazzi]

On Friday, October 31, 2014 5:53:23 PM UTC+1, Jamie Titus wrote:

Hello all,

I was wondering if anybody knew what the status of the Tensor module is - i.e., is work being done on it, what is the known functionality, and what is the intended functionality. For the past week I have been attempting to use Sympy for basic General Relativity (homework problems really), and it has not gone well.

Yes, unfortunately there's still a lot to do before it can be used, just consider that it is completely unable to express derivatives.
 

My naive review of the documentation suggests that I should be able to use Einstein Notation, since I am completely unfamiliar with Penrose's abstract index notation (and it doesn't seem to be as useful computationally, frankly, though the wikipedia article on it is very short).

The idea behind the abstract index notation is simple: the structure of contracted indices is basis-independent, A(i0, -i0) has the same value in all coordinate systems, so someone may want to reason about the index structure (and contractions) without working on a coordinate basis.

However, I keep encountering errors when I try to do basic operations such as add or multiply tensors, and I am never quite sure if they are bugs or features. A minimal example is that when you attempt to symmetrize a tensor you just get back double the tensor components, i.e., given a TensorHead A, and indices i0, i1, A(i0, i1) + A(i1, i0) should be manifestly symmetric, but it just results in 2*A(i0, i1). If necessary I can provide minimal working examples in code. I wasn't sure if this is the right place for that.

Whenever you add tensor expressions, the canonicalizer is automatically called. Given the set of all permissible permutations of indices (and anti commutations by inverting the sign of the tensor), the canonicalizer picks the one which would come first alphabetically, that is A(i1, i0) becomes A(i0, i1). Thus A(i0, i1) + A(i0, i1) gets simplified to 2*A(i0, i1).

It is not clear whether this is the desired behaviour, also considering that such an excessive use of the canonicalizer seriously slows down performance.
 

I honestly found the Tensor module very confusing from the standpoint of somebody who is interested in solving simple GR problems, so I'd like to try to understand the rationale behind it and make things clearer for other users.

Currently, it simply handles basic polynomials of tensors in expanded form. I added the .data property to do some reasoning with component data, but I'm not very glad with it. I think that the tensor module need an extension to support indices representing coordinate systems and basis elements of those coordinate systems. This would require to add support for repeated indices to the canonicalizer.

Anyways, if you would like to do some calculations, I suggest you have a look at sympy.diffgeom. Here an example of how to derive the Schwarzschild metric with diffgeom:

http://krastanov.files.wordpress.com/2012/07/schwarzschild.pdf

[Francesco Bonazzi]

On Friday, October 31, 2014 8:16:35 PM UTC+1, Jamie Titus wrote:

I think the most instructive thing that I found confusing is was my attempt to symmetrize a tensor. This minimal snippet, taken almost entirely from the examples in documentation, should demonstrate my point.

from sympy import *
from sympy.tensor.tensor import TensorIndexType, tensorsymmetry, TensorType

Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
sym2 = tensorsymmetry([1]*2)
S2 = TensorType([Lorentz]*2, sym2)
A = S2('A')
from sympy.tensor.tensor import tensor_indices, tensorhead
Lorentz.data = [1, -1, -1, -1]
i0, i1 = tensor_indices('i0:2', Lorentz)
A.data = [[j+2*i for j in range(2)] for i in range(2)]
A.data
A_symm = A(i0, i1) + A(i1, i0)
A_symm.data

A.data is [[0, 1], [2, 3]], One would want A_symm.data to be [[0, 3],[3, 6]], a symmetric tensor. However, Sympy is computing it as 2*A. Indeed, that is exactly how sympy represents A_symm.

You declare A to be symmetric, but you add non-symmetric data to it. The code should raise an error on the assignment, I am aware of this problem, the solution of which requires the canonicalizer to support repeated indices.

If you call tensorsymmetry([1]*2) which is equivalent to tensorsymmetry([1,1]), you declare the indices to be symmetric, so A(i0, i1) = A(i1, i0). Use tensorsymmetry([1], [1]) if you don't want that.

In [28]: tensorsymmetry([1,1]).generators
Out[28]: (Permutation(3)(0, 1),)

In [29]: tensorsymmetry([1],[1]).generators
Out[29]: (Permutation(3),)

You see, in Out[28] there is a permutation (0, 1), meaning that the first and second indices can be permuted. In Out[29] there is no permutation, meaning that A(i0, i1) would be unequal to A(i1, i0).

I start off with this example because I am genuinely not sure if this is the intended behaviour, since it is unclear to me whether or not the Tensor module is supposed to support typical Einstein notation. If this is a "bug," then there are many things which must be tested and fixed. If this is a "feature," then I probably need to look for a different module to do GR with.

The approach to the tensor module is very ambitious, support both symbolic tensor expressions and component data, be coordinate-independent yet allow to operate on coordinate bases. A more challenging point would be to add support to covariant and partial derivatives. Unfortunately this requires a lot of time. The algorithm to handle the symbolic index structure is quite complicated.

In sympy.diffgeom you simply work on the components of a tensor using the TensorProduct object. You don't have a compact expression as A(i0, i1), but you can do partial, covariant and Lie derivatives, which are currently unsupported in sympy.tensor.tensor.

[Alan Bromborsky]

The following might be of interest. Attached is the documentation for the new geometric algebra/calculus module and a paper on how to apply geometric algebra to classical general relativity (not the gauge theory of gravity).  The module (which uses sympy) is at github.com/brombo/galgebra.  Also attached (manifold.py and manifold.pdf) are examples of tensors on manifolds.  The tensor instantiator is "Mlt" which stands for multilinear transfomation which is away of doing tensors without index notation.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b1b7eb06-13ad-4789-a721-11ec17abd0b2%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Aaron Meurer]

When in doubt, don't do automatic simplification. You can always have
simplify() or doit() or whatever canonicalize the expression. If an
expression doesn't auto-simplify it is always possible to force it to
simplify later, but if it always auto-simplifies, it is impossible to
make it not do so.

I wrote some thoughts on this at
https://github.com/sympy/sympy/wiki/Automatic-Simplification a while
back. This case clearly violates the rule that it should be fast, and
the rule that it should be something that no one would ever not want
to do.

Aaron Meurer

>
>>
>> I honestly found the Tensor module very confusing from the standpoint of
>> somebody who is interested in solving simple GR problems, so I'd like to try
>> to understand the rationale behind it and make things clearer for other
>> users.
>
>
> Currently, it simply handles basic polynomials of tensors in expanded form.
> I added the .data property to do some reasoning with component data, but I'm
> not very glad with it. I think that the tensor module need an extension to
> support indices representing coordinate systems and basis elements of those
> coordinate systems. This would require to add support for repeated indices
> to the canonicalizer.
>
> Anyways, if you would like to do some calculations, I suggest you have a
> look at sympy.diffgeom. Here an example of how to derive the Schwarzschild
> metric with diffgeom:
>
> http://krastanov.files.wordpress.com/2012/07/schwarzschild.pdf
>

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/d7d0eaa7-e616-407e-a27b-a1c4ff2ff6c1%40googlegroups.com.

[Jamie Titus]

Thank you for your reply.

I think I have a few suggestions:

You are correct that I have defined a symmetric tensor and added asymmetric data! However, I took this example almost verbatim from the code examples in the documentation on the website for TensExpr. I believe it would be good to update the examples under TensExpr so that the declared symmetry matches the data. I would also have very much appreciated a simple explanation in the documentation to the BSGS function from the perspective of a non-mathematician.

Second - as a student approaching GR from a physics perspective, the need to declare a symmetry for the indices is very non-intuitive, especially because I have never studied permutations formally, and am therefore completely unfamiliar with the BSGS formalism. Would it be possible to create tensors which have only the identity permutation by default? That way, those of us who are relatively less concerned with canonicalization can use tensors without worrying about conflicts between the data we put in and the symbolic representation.

Third - I tried running my code with the symmetry definition you provided. The symbolic representation was correct, because I got
>>> A_symm

>>> A(i0, i1) + A(i1, i0)

However, the data had the same problem as before, i.e.,
>>> A_symm.data
>>> array([[0, 2], [4, 6]], dtype=object)

The _TensorDataLazyEvaluator does not attempt to permute the indices in the _get function when working with a TensAdd expression. Instead, it just adds the indices in alphabetical order.  So if you attempt

>>> T(i0, i1, i2) = B(i0, i1, i2) + C(i1, i0, i2) + D(i2, i0, i1)
>>> T(i0, i1, i2).data

The result will be incorrect, but more specifically, assuming i0 = 'i0', etc, the _get function will return you the data for
B(i0, i1, i2) + C(i0, i1, i2) + D(i0, i1, i2)

I actually know how to fix this (indeed, when hacking around on my own computer I did fix it) by checking the index positions and judicious use of the np.swapaxes function.

Covariant and partial differentiation I have actually implemented, in a hacky sort of way, but I'm running into bugs in the multiplication that I haven't tracked down yet.

I guess what I want to know is - Are you interested in taking the current sympy.tensor.tensor module and fixing these bugs? I am willing to spend some time implementing/fixing the problems with the _TensorDataLazyEvaluator class, I just don't know what is planned.

With regards to the diffgeom module, I found that to be unintuitive and even less documented than tensor. It was not really clear to me how I was supposed to define arbitrary tensors (i guess build them up from basevector fields? I didn't look into it much), and it is very very nice to have the Einstein index notation. Keep in mind that I am not any kind of GR theorist or Differential Geometer, I'm an experimentalist who just wants an easy way to convert equations into code. Therefore, if you can tell me what needs to be done, I'm willing to work on supporting Einstein Index Notation in sympy.tensor.tensor.

[Francesco Bonazzi]

You are correct that I have defined a symmetric tensor and added asymmetric data! However, I took this example almost verbatim from the code examples in the documentation on the website for TensExpr. I believe it would be good to update the examples under TensExpr so that the declared symmetry matches the data.

Yes, that is definitely an error.

Second - as a student approaching GR from a physics perspective, the need to declare a symmetry for the indices is very non-intuitive, especially because I have never studied permutations formally, and am therefore completely unfamiliar with the BSGS formalism. Would it be possible to create tensors which have only the identity permutation by default? That way, those of us who are relatively less concerned with canonicalization can use tensors without worrying about conflicts between the data we put in and the symbolic representation.

The symmetry requirement could be made optional by filling a default value. Anyway, I would consider this to be low-priority, as there are more serious issues.

BSGS is not a formalism of tensors, it's a generator of the permutation group which is used by an algorithm employed to simplify expressions according to their symmetries.

Consider this example:

In [43]: A = tensorhead('A', [Lorentz]*2, [[2]])

In [44]: B = tensorhead('B', [Lorentz]*2, [[1]*2])

In [45]: expr = A(i0, i1)*B(-i0, -i1)

In [46]: expr.canon_bp()
Out[46]: 0

It clearly recognizes that expr is zero (double contraction between symmetric and antisymmetric tensors).

Third - I tried running my code with the symmetry definition you provided. The symbolic representation was correct, because I got
>>> A_symm

>>> A(i0, i1) + A(i1, i0)

However, the data had the same problem as before, i.e.,
>>> A_symm.data
>>> array([[0, 2], [4, 6]], dtype=object)

The _TensorDataLazyEvaluator does not attempt to permute the indices in the _get function when working with a TensAdd expression. Instead, it just adds the indices in alphabetical order.  So if you attempt

>>> T(i0, i1, i2) = B(i0, i1, i2) + C(i1, i0, i2) + D(i2, i0, i1)
>>> T(i0, i1, i2).data

The result will be incorrect, but more specifically, assuming i0 = 'i0', etc, the _get function will return you the data for
B(i0, i1, i2) + C(i0, i1, i2) + D(i0, i1, i2)

I actually know how to fix this (indeed, when hacking around on my own computer I did fix it) by checking the index positions and judicious use of the np.swapaxes function.

That is clearly a bug. The module sympy.tensor.tensor has undergone very little testing, so it's likely to have various kinds of problems.

By the way, I am not very happy with the current way component data are added to the tensor. The abstract index notation is meant to be basis-free, and defining components on it does not make any sense. I believe there is a need of a new kind of indices, say BasisTensorIndex, which would take a sympy.diffgeom.CoordSys object as a parameter. Obviously TensorIndexType should correspondingly be assigned a sympy.diffgeom.Patch instance. After this, the .data assignment should be permissible only if the indices are concrete symbols, that is, they have a coordinate system basis. One could then read the component data in a new coordinate system by defining a new CoordSys instance, connecting it to the other one as defined in sympy.diffgeom, then creating other "BasisTensorIndex" referred to this new coord system.

Such kind of indices would allow repeated indices, e.g. A(-i, i, i), in the abstract index notation this is not allowed. Obviously this requires support to the canonicalization of repeated indices, or disabling canonicalization for them.

Covariant and partial differentiation I have actually implemented, in a hacky sort of way, but I'm running into bugs in the multiplication that I haven't tracked down yet.

My idea is to create CovarDerivative and PartialDerivative objects. One would use CovarDerivative(i0, A(i1, -i0)), noting that the covariant derivative is OK if expressed with the abstract index notation, while the partial derivative needs the coord-system based index notation.

The problem of this idea is that it is incompatible with the current algorithm that detects the repeated indices.

I guess what I want to know is - Are you interested in taking the current sympy.tensor.tensor module and fixing these bugs? I am willing to spend some time implementing/fixing the problems with the _TensorDataLazyEvaluator class, I just don't know what is planned.

I would like to have a major rewrite of that module, but I don't have much time now. One point is to make TensAdd and TensMul inherit Add and Mul, and adapt the index contraction tools to handle unexpanded tensor expressions, e.g. A(i) * ( B(-i, j) + C(j, -i) ) , and unexpanded derivatives, e.g. PartialDerivative(i, A(-i, j) + B(-i, j)).

By the way, any contribution is welcome. Discussion about the source code usually takes place on github.

With regards to the diffgeom module, I found that to be unintuitive and even less documented than tensor. It was not really clear to me how I was supposed to define arbitrary tensors (i guess build them up from basevector fields? I didn't look into it much), and it is very very nice to have the Einstein index notation.

diffgeom only supports fully covariant tensors by creating a TensorProduct of one-forms, and adding those tensor products. For example, the matrix [[1, 2], [3, 4]] would be represented as

TP(dx, dx) + 2*TP(dx, dy) + 3*TP(dy, dx) + 4*TP(dy, dy)

where dx, dy are the basis one-forms, obviously TP(dx, dy, dz) would be the component of a rank-3 tensor. It's another way to deal with them, by explicitly showing all components.
 

Keep in mind that I am not any kind of GR theorist or Differential Geometer, I'm an experimentalist who just wants an easy way to convert equations into code. Therefore, if you can tell me what needs to be done, I'm willing to work on supporting Einstein Index Notation in sympy.tensor.tensor.

The main problem is to reason about the existing code. I tried to clean it by separating some of the index algorithms from the tensor expression classes, and about one year ago I put it in TIDS, a class you can find on that file.

You could fix the bugs you reported if you like, that's easy. The best is to discuss about code on github, which is better suited for that, feel free to create PR or gists.

[Jamie Titus]

It sounds like the module needs a bit of planning and organization before I just go hacking away and changing things (aside from one or two bugs).

Where on GitHub is the best place to discuss that? I'm completely unfamiliar with GitHub, though I have some experience using Git. It seems to me that it
might be immature to start submitting Pull Requests because we first need to think very carefully about how to implement coordinate dependence
in a rational and efficient way. But I am unfamiliar with the GitHub workflow, so if a Pull Request is the place to discuss that, just let me know.

[Francesco Bonazzi]

Apart from the bug fixes, it's quite hard to go on rigorously.

You wrote you have some hackish stuff about derivatives, on github you can create a gist, which is usually a piece of code unrelated to a project, which people can see and comment.

If you realized, the tensor module only allows expanded polynomial expressions, that is difficultto change because the algorithm was written for that, and supporting something different would require a llot of rewrite.

I have a pending PR on the tensor module, which is a major rewrite of the code, but I'm not very satisfied of it.

An alternative is to expose special objects such as partial derivatives to the index structure of a tensor.

If you like to do something, try to see how to create a covariant derivative object that exposes the index structure.

Say, CovarDerivative(i, A(j, k)) should be treated as a tensor of indices i j k. You will realize that the way the tensor module was written makes this task quite hard.

[Jamie Titus]

On Saturday, November 1, 2014 5:01:04 PM UTC-7, Francesco Bonazzi wrote:

Apart from the bug fixes, it's quite hard to go on rigorously.

You wrote you have some hackish stuff about derivatives, on github you can create a gist, which is usually a piece of code unrelated to a project, which people can see and comment.

I have made a gist with some code I wrote to do derivatives. https://gist.github.com/seniosh/ef0e1550110eec0c6b41#file-grtensor-py where can that then be discussed? I already realize several problems (partial differentiation gets treated as a Tensor, which is a problem if you want to switch coordinates), but it's a
starting point for discussion at least.
 

If you realized, the tensor module only allows expanded polynomial expressions, that is difficultto change because the algorithm was written for that, and supporting something different would require a llot of rewrite.

I have a pending PR on the tensor module, which is a major rewrite of the code, but I'm not very satisfied of it.

Should I move discussion to that PR? Or should I make a new one?
 

An alternative is to expose special objects such as partial derivatives to the index structure of a tensor.

If you like to do something, try to see how to create a covariant derivative object that exposes the index structure.

Say, CovarDerivative(i, A(j, k)) should be treated as a tensor of indices i j k. You will realize that the way the tensor module was written makes this task quite hard.

When I implemented it, my CovarDerivative function operated on data, and then returned a new TensMul object, so CovarDerivative(i, A(j, k)) returns dA(i, j, k). 

[Joachim Durchholz]

Am 02.11.2014 um 01:16 schrieb Jamie Titus:
> I have made a gist with some code I wrote to do derivatives.
> https://gist.github.com/seniosh/ef0e1550110eec0c6b41#file-grtensor-py where
> can that then be discussed?

For pull requests, you can click on a line number and get an input field
to write your comment for that line.
I'd have thought that that's possible in gists, too, but I haven't found
a way to do so, so maybe you still need to set up a Github repo (create
a Github account, fork sympy, git clone the fork to your local disk,
commit your changes to the clone and push them to the fork on Github -
the whole process is a bit lengthy, but Github comes with a very
extensive blow-by-blow tutorial for it so that should work for you).

[Aaron Meurer]

Gists are a bad way to discuss code, because they live in isolation,
and you can't make line comments.

If you have code, you can use a pull request. It is OK to submit a
pull request that shouldn't be merged, but is just for discussion.
Just make that clear in the title of the pull request. Ideally the
code will evolve into something that can be merged, although it
doesn't always.

If there is no code, just discuss either here on the mailing list, or
on a GitHub issue.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/5c27afa1-d279-40d8-bea2-a72aaf4f3d18%40googlegroups.com.

[Francesco Bonazzi]

On Sunday, November 2, 2014 1:16:35 AM UTC+1, Jamie Titus wrote:

On Saturday, November 1, 2014 5:01:04 PM UTC-7, Francesco Bonazzi wrote:

Apart from the bug fixes, it's quite hard to go on rigorously.

You wrote you have some hackish stuff about derivatives, on github you can create a gist, which is usually a piece of code unrelated to a project, which people can see and comment.

I have made a gist with some code I wrote to do derivatives. https://gist.github.com/seniosh/ef0e1550110eec0c6b41#file-grtensor-py where can that then be discussed? I already realize several problems (partial differentiation gets treated as a Tensor, which is a problem if you want to switch coordinates), but it's a
starting point for discussion at least.

Component operations only are already partly implemented in sympy.diffgeom, tensors there are fully-covariant only, unfortunately. It would be clean to share the same algorithm between those two modules, as the operations are practically the same when acting on components, sympy.diffgeom simply has a different way to show the same objects to the end-user.

Furthermore, in sympy.tensor.tensor one may simply wish to manipulate tensor formulae without any knowledge of the components of the tensors being handled. Maybe you would like to express the covariant derivative of a tensor whose components are unknown to you.

If you want to operate on components of tensors I suggest you to look closer to sympy.diffgeom, you may be unfamiliar with its syntax, but it's not that complicated to learn it. It's just a different syntax to express the same things concerning components of tensors. The function twoform_to_matrix clearly illustrates this kind of correspondence, where you transform a summation of tensor products to a matrix (matrix = 2-rank tensor without valence markings). Try to play with that.

 

If you realized, the tensor module only allows expanded polynomial expressions, that is difficultto change because the algorithm was written for that, and supporting something different would require a llot of rewrite.

I have a pending PR on the tensor module, which is a major rewrite of the code, but I'm not very satisfied of it.

Should I move discussion to that PR? Or should I make a new one?

My pending PR is https://github.com/sympy/sympy/pull/7762

I am not happy with it, and I think I am going to close it, because the performance loss is significant.

The ideas I introduced with it are still valid. The current tensor module handles polynomial expressions of tensors, that is TensMul and TensAdd. TensMul contains a list of TensorHead and TensorIndex objects, while TensAdd is a list of TensMul objects.

All operations, that is, multiplication, addition, canonicalization, assume this strict structure, TensMul is supposed to have a defined index order while TensAdd does have it. All operations assume this strictly, so introducing objects such as CovarDerivative requires major rewrites, and that's the reason why my PR has been pending for a long time.

My idea is to introduce the following API which all classes inheriting TensExpr should have:

• .has_index_order()  specifying whether the expression has an index order, A(i0, i1) has it, A(i0, i1) + B(i1, i0) does not have it.
• .free_indices_list() returns a list of the free indices if the order exists, otherwise raises an error.
• .free_indices_set() returns a set (=collection without order) of indices, should always work.
• equivalent expressions for dummy indices.
• indices symmetry information and methods should raise an error if called on expressions without index order. Currently they are all implemented on TensMul, but future objects such as CovarDerivative may also need all those methods.
• no abstract class related to index order: TensMul should be expressible in forms without index order, such as A(i)*(B(j, k) + C(k, j)), it is not a class-feature.
  

All current algorithms assume that TensMul always has index order and that it's the only one allowed to possess it. Obviously, if we introduce CovarDerivative(i0, A(i1, i1)), this should be exposed as a tensor with index order with indices [i0, i1, i2], which is not a TensMul.

 

An alternative is to expose special objects such as partial derivatives to the index structure of a tensor.

If you like to do something, try to see how to create a covariant derivative object that exposes the index structure.

Say, CovarDerivative(i, A(j, k)) should be treated as a tensor of indices i j k. You will realize that the way the tensor module was written makes this task quite hard.

When I implemented it, my CovarDerivative function operated on data, and then returned a new TensMul object, so CovarDerivative(i, A(j, k)) returns dA(i, j, k). 

Operation on data would fit more on sympy.diffgeom, which already supports partial, covariant and Lie derivatives. If you wish to produce a fast implementation of tensor objects as multidimensional arrays with valence markings, I would suggest to create a new class, say TensorArray, inside sympy.diffgeom, providing a compact way to represent tensor components and information on the coordinate system you're working on, and valence (i.e. covariantness/contravariantness) markings.

Supplying sympy.tensor.tensor with all of its missing features is much harder and requires a lot of time, but if you wish to play with that, I can support you.

[Jamie Titus]

I think I understand the intended structure of the tensor.tensor module a bit more now. A lot of the things I have fooled around with are quick hacks that don't really fit the use of each class (i.e., just returning a TensMul from a covariant derivative function). A quick question - why is it true that TensAdd objects are not supposed to support index order? Is this a function of the other algorithms like canonicalization that I'm not as interested in?

I'd like to ask more generally what use a tensor is if the indices are unordered, but I suppose that just reveals my ignorance as a novice physicist - I think of Tensors in a pretty concrete matrix-based way, and when I'm not thinking of them as matrices, I'm thinking of them as functions, which requires an argument order.

I will take your advice and look more at sympy.diffgeom for the time being - do you know if it's possible to implement mixed tensors there? Perhaps it would be easier to build a TensorArray class that provides an interface more like sympy.tensor.tensor while using the implementation of sympy.diffgeom.

Thank you for taking the time to explain everything about the structure of the module.

[Francesco Bonazzi]

On Monday, November 3, 2014 5:10:11 PM UTC+1, Jamie Titus wrote:

I think I understand the intended structure of the tensor.tensor module a bit more now. A lot of the things I have fooled around with are quick hacks that don't really fit the use of each class (i.e., just returning a TensMul from a covariant derivative function). A quick question - why is it true that TensAdd objects are not supposed to support index order? Is this a function of the other algorithms like canonicalization that I'm not as interested in?

If you have an expression like A(i, j) + B(j, i) - what is the index order? [i, j] or [j, i]? It exposes a set of free (i.e. uncontracted) indices, {i, j}. The single terms of the addition have an order, [i, j] and [j, i] respectively.
 

I'd like to ask more generally what use a tensor is if the indices are unordered, but I suppose that just reveals my ignorance as a novice physicist - I think of Tensors in a pretty concrete matrix-based way, and when I'm not thinking of them as matrices, I'm thinking of them as functions, which requires an argument order.

A tensor expression may not have a clear order. When you think of tensors as multidimensional arrays, you have to pick an index order, so in A(i, j) + B(j, i) you may pick the order [i, j], transpose the matrix of B's components and add it to A. It currently fails to transpose, this is the bug you reported. If you picked [j, i] instead of [i, j], you would operate on the transpose.

I will take your advice and look more at sympy.diffgeom for the time being - do you know if it's possible to implement mixed tensors there? Perhaps it would be easier to build a TensorArray class that provides an interface more like sympy.tensor.tensor while using the implementation of sympy.diffgeom.

Yes, I tried to draft a TensorArray some months ago, but then ran out of time. In sympy.diffgeom there is a wedge product, which makes sense only on differential forms, which transform covariantly.

Anyways, creating a TensorArray is the fastest approach to GR. I would devise it something like

class TensorArray(Expr):
    def __new__(cls, ndarray, shape, valence, coord_sys=None):
        return Expr.__new__(ndarray, shape, coord_sys)

def numpy_to_TensorArray(ndarr, valence, coord_sys=None):
    ...

Thank you for taking the time to explain everything about the structure of the module.

You're welcome.

[Francesco Bonazzi]

The documentation bug has been corrected on the master branch.

Concerning the components data of A(i, j) + A(j, i), I thought that the faster solution it to remap the axes of the numpy.ndarray to the lexicographically sorted indices. That is, the components data of both A(i, j) and A(j, i) would both be sorted as (i, j). This would break compatibility with A(j, i).data, which would get transposed on future versions of SymPy.

An alternative is to make indexing mandatory when there is no evident index order. That is, A(i, j).data is OK, while (A(i, j) + A(j, i)).data would raise an error. Maybe a new method could be introduced, say .get_data_by_index_order(i, j)