Bug in multiple tensor contractions with scalar result

[Eric Gourgoulhon]

Tensor contractions on multiple indices give wrong results when the outcome is a scalar and the index positions in each tensor do not have the same order. For instance, let us consider

sage: M = Manifold(2, 'M')                                                                          
sage: X.<x,y> = M.chart()                                                                           
sage: a = M.multivector_field(2)                                                                    
sage: a[0,1] = 1                                                                                    
sage: b = M.diff_form(2)                                                                            
sage: b[0,1] = 1                                                                                    

The double contraction a^ij b_ij is obtained via

sage: a.contract(0, 1, b, 0, 1).expr()                                                              
2
sage: (a['^ij']*b['_ij']).expr()  # equivalent to above                                                                 
2

which is correct. However, asking for the double contraction a^ij b_ji leads to

sage: a.contract(0, 1, b, 1, 0).expr()  
2
sage: (a['^ij']*b['_ji']).expr()  # equivalent to above                                                                     
2

which is obviously wrong: the result should be -2, since b is antisymmetric.

If the outcome is not a scalar field, there is no issue:

sage: c = M.tensor_field(3, 0, antisym=(0,1))                                                       
sage: c[0,1,0] = 1                                                                                  
sage: c.contract(0, 1, b, 0, 1).display()  # correct result                                                        
2 ∂/∂x
sage: c.contract(0, 1, b, 1, 0).display()  # correct result                                                           
-2 ∂/∂x

This bug has been reported in ​https://ask.sagemath.org/question/58379.

[Eric Gourgoulhon]

The culprit is line 2187 of src/sage/tensor/modules/comp.py:

    for ind in comp_for_contr.index_generator():
        res += self[[ind]] * other[[ind]]

ind should be reordered before being passed to other[[...]].

[Eric Gourgoulhon]

The bug is fixed in the above branch. Some slight clean up of the code for parallelized contractions has also been performed.

---

New commits:

​26ba8f7

Fix bug in Components.contract() (Trac #32355)

[Michael Jung]

Here are some comments.

• I think you don't need to initialize a new component in

          comp_for_contr = Components(self._ring, self._frame, ncontr,
                                      start_index=self._sindex)

and you can use self.index_generator() right away, no?

• What about this?

-            for ind_s in comp_for_contr.index_generator():
-                ind_o = [None for i in range(ncontr)]
-                for pos_s, pos_o in contractions:
-                    ind_o[pos_o] = ind_s[pos_s]
-                ind_pairs.append((ind_s, ind_o))
+            for ind_s in self.index_generator():
+                ind_o = [ind_s[pos_s] for pos_s, _ in sorted(contractions, key=itemgetter(1))]
+                ind_pairs.append((ind_s, ind_o))

You could also turn the outer for-loop into a list comprehension which would be faster but I reckon less well readable.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3e5ac7d

Improvement in the fix for #32355

[Eric Gourgoulhon]

Replying to gh-mjungmath:

Here are some comments.

Thanks for taking a look into this.

• I think you don't need to initialize a new component in

          comp_for_contr = Components(self._ring, self._frame, ncontr,
                                      start_index=self._sindex)

and you can use self.index_generator() right away, no?

Yes indeed! (this is needed for the generic case, but not for the specific case of a scalar output). The latest commit implement this change.

• What about this?

-            for ind_s in comp_for_contr.index_generator():
-                ind_o = [None for i in range(ncontr)]
-                for pos_s, pos_o in contractions:
-                    ind_o[pos_o] = ind_s[pos_s]
-                ind_pairs.append((ind_s, ind_o))
+            for ind_s in self.index_generator():
+                ind_o = [ind_s[pos_s] for pos_s, _ in sorted(contractions, key=itemgetter(1))]
+                ind_pairs.append((ind_s, ind_o))

I would prefer to keep the original version, since it is more efficient from a computational point of view: there is no need to sort a list. Granted: in most use cases, the list contains only a few elements and sorting shall be fast. So it is more a matter of taste...

[Michael Jung]

Replying to egourgoulhon:

I would prefer to keep the original version, since it is more efficient from a computational point of view: there is no need to sort a list. Granted: in most use cases, the list contains only a few elements and sorting shall be fast. So it is more a matter of taste...

You're right, this is faster! Then I still have one last comment:

-                ind_o = [None for i in range(ncontr)]
+                ind_o = [None] * ncontr

This is indeed faster.

If that is done, and patchbot is happy -> positive review.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​48f5523

Yet another improvement in the fix for #32355

[Eric Gourgoulhon]

Replying to gh-mjungmath: Then I still have one last comment:

-                ind_o = [None for i in range(ncontr)]
+                ind_o = [None] * ncontr

Done in the above commit.

[Eric Gourgoulhon]

Thank you for the review!