Tensor refactory

[F. B.]

I came to the conclusion that the best solution to implement Quantum Field Theory in SymPy is to use the unify module combined with the new tensor module.

Unfortunately the new tensor module has a complicated arg-tree structure. The main problems concerning the implementation of QFT (let's ignore General Relativity for now, which would make everything even more complicated) are:

• Operators on tensors
• Unification/rewrite rules on tensors

Operators are the most urgent issue. The current arg-tree and internal structure of TensMul and TensAdd do not make it easy to add new objects containing tensor expressions. My current idea is to refactor both TensMul and TensAdd, and relocate all index and symmetry manipulation logic to either a TIDS object or to a new object (let's say TensorContainer).

The TensorContainer option would be the best one, so that TensMul and possible future operators would all inherit that object. All index symmetry operations and indices contractions would be handled by TensorContainer, so that TensMul is no longer be the only such container. TensAdd is then a collection of TensorContainer. A derivative operator would then contain the derivative variable and a TensMul object, but I still have to make up my mind on how to exactly represent operators.

What about not storing index symmetry information inside the arg-tree? I mean, it could work in a way similar to the new assumptions, with a global dictionary storing such info. I am currently proposing such an approach for components data in a PR.

[F. B.]

The second point concerns unification of tensor, there are some issues on contracted indices, with some boilerplate

from sympy.tensor.tensor import *
L = TensorIndexType('L')
i0, i1, i2 = tensor_indices('i0:3', L)
A = tensorhead('A', [L]*2, [[1]*2])

The arguments

>>> A(i0, i1).args
(1, (A(L,L),), (i0, i1))

>>> A(i0, -i0).args
(1, (A(L,L),), (dummy_index_0, -dummy_index_0))

That is, the third argument is the list of indices, and in the contracted form the index i0 gets replaced by dummy_index_0, which makes unification hard:

>>> from sympy.unify import *

>>> list(unify(A(i0, i1), A(i0, i2), {}, variables=[i2]))
[{i2: i1}]
>>> list(unify(A(i0, -i0), A(i1, -i1), {}, variables=[i1]))
[{}]

The second unify expression does not return anything because the arguments i0 and i1 have been substituted. A possible solution is to keep the arguments of a TensMul as they have been created, and store the dummy_index_0 inside the corresponding TIDS object.

A possible implementation of TensorDerivative (by which I mean the derivatives usually employed to get the equation of motions from the Lagrangian in QFT) can simply rely on rewriterule contained in the unify module.

The ability to define new objects such as TensorDerivative, PartialDerivative, CovariantDerivative is the primary reason why I want to make the tensor expression objects just containers, and possible transfer all index manipulation algorithms to either TIDS or a new object.

[Matthew Rocklin]

Maybe tensor indices should be identified as variables/wilds in unification ?  Should A(i, j, -j) unify to A(a, b, -b)?

--
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/0bcab04e-6d84-465e-8fa8-2592b28b8b3c%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Ondřej Čertík]

Do you want to do a G+ hangout about this? Looks like there are couple
options to go about this, it might be the fastest to brainstorm this
face to face.

Ondrej

> https://groups.google.com/d/msgid/sympy/CAJ8oX-E4QCqars8wyu1fnckpW1Xfaiyk3-FPmvpfX-Sy9KZiHw%40mail.gmail.com.

[Matthew Rocklin]

I'm off of work this week so I'm generally around.  As a warning I've been playing around with logpy recently and will probably push that a little bit.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVC2dOYvobCDcjMBySQx8TS-j7peBsxk2NmqKBrMWqRBQQ%40mail.gmail.com.

[Stefan Krastanov]

> I came to the conclusion that the best solution to implement Quantum Field
> Theory in SymPy is to use the unify module combined with the new tensor
> module.

I can not agree more with statement. I would be the happiest person on
Earth if you/the community succeeds at such a refactoring. I actually
tried to do it on one or two occasions (it actually happened exactly
when you joined us), but such a task is beyond me and regrettably I
will not have free time in the near future to help.

Sorry for this rather nonconstructive post, I just wanted to say that
I am happy that somebody is thinking about this.

[Matthew Rocklin]

Stefan, can you give a brief explanation of your experience?  In particular what difficulties did you run into?

--
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/CAD8szLyr1qqFoKwUcOOCrcfs%2B%3DeretJ8g1bNAnb5CTJ0HrcY-w%40mail.gmail.com.

[Stefan Krastanov]

> Stefan, can you give a brief explanation of your experience? In particular
> what difficulties did you run into?

The current tensor canonicalization and representation for tensors
does not adhere to most of sympy's class invariants.

- eq is occasionally about mathematical equality, not structural one,
and it is even raising errors in some cases, which break dictionary
lookups and so on

- you can not recreate objects using func and args

- substitution does not work (nothing related to tree traversal works actually)

- many others in the same style

In one sentence: the implemented algorithms are very smart, but they
completely disregard the invariants that sympy uses to work with the
symbolic trees.

I tried to go step by step and reinstate these invariants, but my
knowledge of the module was insufficient. (actually, going in the
other direction might be easier: the module when I was looking at it
was not using any of sympy's methods, so it was not really necessary
to subclass Basic or Expr).

In any case Franz is more aware of these issues than me and he was
working on fixing some of them (I was not following the work so I do
not know how much of this is still relevant).

> https://groups.google.com/d/msgid/sympy/CAJ8oX-EoyBb3XcdCSsEeqBV-9su%2B7sMB8etZZmfwmsdeGm3XPg%40mail.gmail.com.

[Aaron Meurer]

On Fri, Mar 28, 2014 at 9:13 PM, Stefan Krastanov
<stefan.k...@yale.edu> wrote:
>> Stefan, can you give a brief explanation of your experience? In particular
>> what difficulties did you run into?
>
> The current tensor canonicalization and representation for tensors
> does not adhere to most of sympy's class invariants.
>
> - eq is occasionally about mathematical equality, not structural one,
> and it is even raising errors in some cases, which break dictionary
> lookups and so on
>
> - you can not recreate objects using func and args
>
> - substitution does not work (nothing related to tree traversal works actually)
>
> - many others in the same style
>
> In one sentence: the implemented algorithms are very smart, but they
> completely disregard the invariants that sympy uses to work with the
> symbolic trees.
>
> I tried to go step by step and reinstate these invariants, but my
> knowledge of the module was insufficient. (actually, going in the
> other direction might be easier: the module when I was looking at it
> was not using any of sympy's methods, so it was not really necessary
> to subclass Basic or Expr).

If you're using .args you might as well subclass from Basic, because
you get a bunch of useful methods like xreplace, __eq__, and __hash__
for free. There are maybe still some methods on Basic that don't
belong there, but we should just work on moving them out.

And if you're not using .args then you don't adhere to the basic
invariants of SymPy.

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAD8szLw1M6E-OLrhsqqW0_RsocXywtF8Ost2%3D3kD7Vp4-3O_Kw%40mail.gmail.com.

[F. B.]

Actually I refactored the tensor module the last autumn and it now should have invariant argument trees.

I don't know whether __eq__ overloading may be an issue to unify, as well as the standard __hash__ for tensor expressions.

The issue now is more about how to represent operators on tensors in a frozen state, eg partial derivative by an index or by a tensor.

Note that now some of the logic is also handled by an object called TIDS, whose purpose is just to contain the calculation data pertaining to TensMul.

[F. B.]

 Currently, __eq__ and __hash__ are overwritten, in order to make the == operator correspond the tha mathematical equality. I see this is not standard in SymPy, since math equality is expressed by Eq(  ,  ) , hope this won't cause much trouble to unify/rewriterule, otherwise the == behavior will have to be changed as well.

 

And if you're not using .args then you don't adhere to the basic
invariants of SymPy.

This issue has already been solved.

[F. B.]

By the way, when I was referring to an external tensor index symmetry container, I meant something like what I'm doing here for the tensor components:

https://github.com/sympy/sympy/pull/2805

This PR solves the issues of associating tensor components data to tensor objects. Previously they were injected into the tensor classes, this created problems as that info was lost when the tensor was rebuilt through .func . Now there is an external dictionary mapping the tensor objects to the components data (which are numpy ndarrays). In any case, what about this PR? Any further suggestions/issues?

The tensor index symmetry could as well be stored into a global dictionary, that's an idea. But in reality, I think the unify/rewriterule would seamlessly work even if the tensor index symmetry is still stored inside the arg-tree.

Anyway, please review #2805, so that we can finally solve some urgent issues about components data of tensors, as well as some bugs previously reported on this mailing list.

[Stefan Krastanov]

Here is an example why `__eq__` should correspond to structural equality:

>>> class a(object):
... def __eq__(self, other):
... raise NotImplementedError()
...

>>> i = a()

>>> list = [1,2,3,i]

>>> list.index(i)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "<stdin>", line 3, in __eq__
NotImplementedError



__eq__ is part of basic python protocols, we should not make it do
other things.

[F. B.]

OK, so as soon as the currently pending PR is merged, I'll try to address this issue.

I think that the problem can be solved by renaming __eq__ to _eval_Eq and then solving all issues of confusion around.

[F. B.]

By the way, I missed another point in my previous posts, we also need a representation for tensor fields, i.e. tensors depending on another scalar/tensor, like $$ A^\mu (x^\nu) $$.

[F. B.]

OK, I addressed the issue pointed out by Krastanov, namely that __eq__ should be a structural equality, and not a mathematical equality:

https://github.com/sympy/sympy/pull/7357

[F. B.]

OK, after the last PR has been merged, I did some tests on the indices and there are still problems with .match( ) and .replace( ), while .subs( ), .replace(..., simultaneous=False), and .xreplace( ) seem to work. The current problem is that wilds are not allowed in the expression tree, because the object containing the tensor indices tries to access internal properties of the indices, which Wild does not possess, thus raising an error.

I did a grep of the unify module, and it looks like .match and .replace are not used, so this should not be an urgent issue.

I don't know whether it's worth to address this issue, but I envisioned an attempt to leave the current procedural programming style to approach a more functional style. There should also be a refactory in order to allow operators such as the partial derivative to act on tensor expressions, possible ideas:

• create a DummyTensorIndex type, similar to TensorIndex, dummy indices would be stored into the arg-tree as instances of DummyTensorIndex instead of TensorIndex.
• it makes no sense anymore to consider the index order until the operators are unevaluated, e.g.
  PartialDerivative(A(i)*B(j, k)*A(-k) + C(-k, j, i)*A(k), A(h))
  
  has a (i, j, -h) index signature, but it makes no sense to consider the index order. Yet one could wish to multiply the unevaluated operator by a tensor and perform a contraction, e.g.
  D(-h)*PartialDerivative(A(i)*B(j, k)*A(-k) + C(-k, j, i)*A(k), A(h))
  
  This should be a TensMul object, but it cannot possess positional info. Solving this issue would also pave the way to unexpanded tensor expressions.
• Should DummyTensorIndex and TensorIndex have a wild=True/False parameter to solve the issue of matching contracted indices?

In any case, introducing operators requires some refactory of the current index algorithms. Fortunately there is now a separate TIDS class used to stored all such algorithms.

[Matthew Rocklin]

Yeah, Wild subclasses from Expr, alas.  The logpy solution is to have a dispatched isvar function.

It seems to me that the all tensor indices are in some sense wild.  Shouldn't A[i, -i] unify to A[j, -j] ?  

--

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/eb0fdec5-18e6-47c3-b631-e3a791ab2f6e%40googlegroups.com.

[F. B.]

On Wednesday, April 2, 2014 7:35:45 PM UTC+2, Matthew wrote:

Yeah, Wild subclasses from Expr, alas.  The logpy solution is to have a dispatched isvar function.

It seems to me that the all tensor indices are in some sense wild.  Shouldn't A[i, -i] unify to A[j, -j] ?  

Actually not, that was a particular case because A(i, -i) is A(j, -j), and this is also equivalent to a scalar (j and -j are summed over). A(i, j) instead is equivalent to a matrix. If A has associated components data, A(i, j).data is the N x N ndarray, and A(i, -i) is the trace.

Consider A(i, -i, j) and A(j, -i, i), these are not equal and should not unify, not even if both i and j are wild. Of course a tensor constructor would raise an exception if an index is repeated more than twice.

[Matthew Rocklin]

Right, but ideally A(i, -i, j) and A(j, -i, i) wouldn't unify.  Actually in this sort of case I suppose that they would because it could be that i == j.  Generally speaking though unification variables do need match consistently within a term.  (a, a) does not match to (1, 2).  Perhaps we could consider all tensor indices on one side to be wild?

--

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/d7b6882e-d224-48b9-8261-bd5ff81cb84c%40googlegroups.com.

[F. B.]

On Wednesday, April 2, 2014 9:56:29 PM UTC+2, Matthew wrote:

Right, but ideally A(i, -i, j) and A(j, -i, i) wouldn't unify.  Actually in this sort of case I suppose that they would because it could be that i == j.

Mmm, you're right. I didn't consider this problem. I believe that some precautions should be taken when using the unification, in order to avoid nonsense expressions, i.e. unification results should first be verified. I didn't test this on rewriterule, but there could be an uncaught exception.
 

Generally speaking though unification variables do need match consistently within a term.  (a, a) does not match to (1, 2).  Perhaps we could consider all tensor indices on one side to be wild?

Well, maybe it's better to continue this discussion as soon as there will be some code attempt to apply rewriterule on tensors. There are still some issues to solve, after that I would like to try to implement a partial derivative operator on tensors with rewriterule.
 

[Aaron Meurer]

Yes, exactly. In Python, __eq__ really means equality, and things like sets and lists will assume objects that return True with __eq__ can be interchanged with one another. SymPy follows this convention. If a == b, any SymPy algorithm might interchange a with b and consider it to be valid.

Also take a read of https://github.com/sympy/sympy/wiki/Automatic-Simplification. It's best to allow different things ways of writing mathematically the same thing to be written down differently, and deal with them being mathematically equal in algorithms. 

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/CAD8szLxfMiWvhFxkELqSwFdJTauce8jr0%3DHz-AmiCA3UnGs%2BSw%40mail.gmail.com.

[F. B.]

On Wednesday, April 9, 2014 1:36:56 AM UTC+2, Aaron Meurer wrote:

Yes, exactly. In Python, __eq__ really means equality, and things like sets and lists will assume objects that return True with __eq__ can be interchanged with one another. SymPy follows this convention. If a == b, any SymPy algorithm might interchange a with b and consider it to be valid.

Also take a read of https://github.com/sympy/sympy/wiki/Automatic-Simplification. It's best to allow different things ways of writing mathematically the same thing to be written down differently, and deal with them being mathematically equal in algorithms. 

I also share the feeling that canonicalization upon construction is wrong, especially as it is computationally expensive for tensors. TensAdd is currently canonicalizing everything.

I tried to make TensMul and TensAdd inherit Mul and Add, it looks great to use all polynomial manipulation stuff with tensors, except that methods inherited by Mul and Add completely ignore indices (this is bad). For example, if TensMul inherits Mul, it is possible to use tens_mul_instance.expand(), unfortunately the resulting expression is not a TensAdd, it's an Add object. Is there any way to overcome this issue? Tensor expressions should really be handled to some degree as polynomials.

[Aaron Meurer]

There are probably little ways around these things, but nothing clean without dispatching in the core.

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/e9a07348-9bb1-48af-a93c-7e458bb5f0f2%40googlegroups.com.

[F. B.]

On Friday, April 11, 2014 1:29:11 AM UTC+2, Aaron Meurer wrote:

There are probably little ways around these things, but nothing clean without dispatching in the core.

Another point in favor of multiple dispatching.

By the way, tensor expressions should just become ordinary expression with the addition of an index management mechanism, as well as other features such as components data association.

Of course some precautions should be taken, for example all indices have to be contracted if you take the exponential of a tensor expression (I am wondering, did anyone ever define a unique and consistent way to generalize the matrix exponential to tensors of any rank?).

[Aaron Meurer]

This has always been the original motivation for multiple dispatch, at
least in my mind. You can't make objects that do their own thing in
Add or Mul or whatever. There are a dozen example of this throughout
SymPy, and a dozen more in user code. There are many hacky ways around
it, but none are satisfactory.

The problem is, how do you dispatch Add(*args). Any argument of the
Add might want to do anything with any other argument. You don't want
to require that arguments be next to each other, because then
something as simple as Add(yourobject, 0, yourobject) wouldn't do the
right thing. You can do the n**2 passes, but does it remain efficient
at that point?

I think we should just start to play with this, especially now that we
have a decent implementation of multiple dispatch. I'd personally
rather play with this with a module that I can understand (so, e.g.,
matrix expressions rather than tensores), but anything is better than
nothing.

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/b7592a77-00ef-4c82-a94b-f0e904581c63%40googlegroups.com.

[Matthew Rocklin]

Sets is a good play module for this purpose.  Union and Intersection simplify pairwise much like Add and Mul would.  I like the pairwise methods we have but the strategy to orchestrate them could use some cleaning up.  Might be a good place to build intuition.  

See https://github.com/sympy/sympy/pull/2979

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6L%2B0wx_h1YFT83swB90WwYq_fL3W6q-61frXxr-rUo0YQ%40mail.gmail.com.

[Aaron Meurer]

Oh, that's even better. Then we don't have to touch the core, and we
can merge even experimental things faster.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-Fh-a3sJ9NFRG%3DNt1bxcbkqcTLzcOPTMjGU4Hy%2BxEmdAg%40mail.gmail.com.

[Matthew Rocklin]

Well, sets is in the core.  Interval, Union, and FiniteSet are all used in various parts.  Fortunately though I think we refactor without changing functionality.  So far I've just replaced a custom double dispatch system with multipledispatch.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6%2BOOW5iHNsOCbFW26A6zh7j%3DnCm57pNTu8hTCv8m81i5Q%40mail.gmail.com.

[Matthew Rocklin]

The benefit of sets is that everyone can understand it but it still has enough complex elements (ProductSet, EmptySet) to capture the complexity of a lot of the more challenging systems (though probably not quite as challenging as tensors.)

[Aaron Meurer]

I don't mean the core as in sympy.core. I mean the core as in the base
classes Basic, Expr, Add, Mul (and to some extend the numbers). Stuff
that *everything* uses. Not everything uses the sets. Unless I'm
missing some very core use of the sets.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-FgYuCqRHmwEPL%2Bt-5YPTmB8XPwzcZevOpUNgAAcUFp7g%40mail.gmail.com.

[Matthew Rocklin]

All I know is that Interval and Union were there before I got here.  Maybe integrals?

I've changed things in FiniteSet and had some core things break.

But yes, Set is not Expr or Basic.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6KynwwTuTnPtEzKL45vgUg5Q00zdxXc-jSQb9K_L8LtNA%40mail.gmail.com.

[Aaron Meurer]

I think Interval was just added for its own sake.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-Hrb-Uq0T5fqauTDJyva0eQzMg4vhHyM29xf%2B%3D_FiK%2BEw%40mail.gmail.com.