Opened 9 months ago
Closed 8 months ago
#32355 closed defect (fixed)
Bug in multiple tensor contractions with scalar result
Reported by: | egourgoulhon | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-9.5 |
Component: | manifolds | Keywords: | tensor contraction |
Cc: | tscrim, gh-mjungmath, mkoeppe | Merged in: | |
Authors: | Eric Gourgoulhon | Reviewers: | Michael Jung |
Report Upstream: | N/A | Work issues: | |
Branch: | 48f5523 (Commits, GitHub, GitLab) | Commit: | 48f5523074fbbb835fa41be1d567d244e9244ae4 |
Dependencies: | Stopgaps: |
Description (last modified by )
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 aij bij 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 aij bji 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.
Change History (14)
comment:1 Changed 9 months ago by
- Description modified (diff)
- Summary changed from Bug in multiple tensor contractions to Bug in multiple tensor contractions with scalar result
comment:2 Changed 9 months ago by
comment:3 Changed 9 months ago by
- Branch set to public/manifolds/contraction_bug-32355
- Cc tscrim gh-mjungmath mkoeppe added
- Commit set to 26ba8f76d35623c499ce7c70a0b6a2afbb724929
- Status changed from new to needs_review
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)
|
comment:4 Changed 9 months ago by
comment:5 follow-up: ↓ 7 Changed 9 months ago by
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.
comment:6 Changed 9 months ago by
- Commit changed from 26ba8f76d35623c499ce7c70a0b6a2afbb724929 to 3e5ac7d48bc15509694e0b904b8595b38caa4e8a
Branch pushed to git repo; I updated commit sha1. New commits:
3e5ac7d | Improvement in the fix for #32355
|
comment:7 in reply to: ↑ 5 ; follow-up: ↓ 8 Changed 9 months ago by
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...
comment:8 in reply to: ↑ 7 ; follow-up: ↓ 10 Changed 9 months ago by
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.
comment:9 Changed 9 months ago by
- Commit changed from 3e5ac7d48bc15509694e0b904b8595b38caa4e8a to 48f5523074fbbb835fa41be1d567d244e9244ae4
Branch pushed to git repo; I updated commit sha1. New commits:
48f5523 | Yet another improvement in the fix for #32355
|
comment:10 in reply to: ↑ 8 Changed 9 months ago by
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.
comment:11 Changed 9 months ago by
- Reviewers set to Michael Jung
- Status changed from needs_review to positive_review
comment:12 Changed 9 months ago by
Thank you for the review!
comment:13 Changed 9 months ago by
- Milestone changed from sage-9.4 to sage-9.5
comment:14 Changed 8 months ago by
- Branch changed from public/manifolds/contraction_bug-32355 to 48f5523074fbbb835fa41be1d567d244e9244ae4
- Resolution set to fixed
- Status changed from positive_review to closed
The culprit is line 2187 of
src/sage/tensor/modules/comp.py
:ind
should be reordered before being passed toother[[...]]
.