Repository: sympy/sympy
Pull request: hep: Fixes bug in gamma_trace of sum of products containing GammaMatrix mixed with other symbols
#### References to other Issues or PRs

Fixes #13636 

#### Brief description of what is fixed or changed

<!-- BEGIN RELEASE NOTES -->
* physics.hep
  * Fixed a bug in finding traces of sums of products of GammaMatrix mixed with other factors. (#13636)
<!-- END RELEASE NOTES -->


Problem statement:
Strangely behaved gamma_trace (gamma matrix trace)
I am trying to use the high energy physics module of SymPy to compute traces of gamma matrices using the dedicated function `sympy.physics.hep.gamma_matrices.gamma_trace`.  I have come accross the following strange behaviour: in some cases, the trace of a sum is not computed to be the sum of the traces. Here is an example:

```python
from sympy import Symbol
from sympy.tensor.tensor import tensor_indices, tensorhead
from sympy.physics.hep.gamma_matrices import GammaMatrix as G, gamma_trace, LorentzIndex

pi, ki, pf = tensorhead("pi, ki, pf", [LorentzIndex], [[1]])
i0, i1, i2, i3, i4 = tensor_indices("i0:5", LorentzIndex)
x = Symbol("x")
pis = pi(i2) * G(-i2)
kis = ki(i3) * G(-i3)
pfs = pf(i4) * G(-i4)

# Show the trace of A + B
print gamma_trace(pfs * G(i0) * kis * G(i1) * pis * G(-i1) * kis * G(-i0) + pfs * G(i0) * kis * G(i1) * pis * x * G(-i0) * pi(-i1))
# Show the trace of A
print gamma_trace(pfs * G(i0) * kis * G(i1) * pis * G(-i1) * kis * G(-i0))
# Show the trace of B
print gamma_trace(pfs * G(i0) * kis * G(i1) * pis * x * G(-i0) * pi(-i1))
```

It may be a normal behaviour, but I couldn't find any hint in the docs dedicated to the module.

Hints from issue discussion:
`gamma_trace(B)` differs from that part of `gamma_trace(A + B)` which comes from `B` because it is processed by `_simplify_single_line`. It calls `extract_type_tens` in order to collect all gamma matrices, but only the explicit gamma matrix factors are found, not those embedded in `pis, kis, pfs`. Hence it seems that the result of `gamma_trace(B)` is wrong.
The result of `gamma_trace(B)` is actually right and similarly for `gamma_trace(A)`. Only `gamma_trace(A + B)` is problematic. Also, I am not sure I understand the argument about the "explicit gamma factors". Do you mean there is a problem with using the intermediary variables `pis, kis, pfs` ? Yet, it seems the documentation gives similar examples.

All in all, should this behaviour be considered a bug? And in any case, what is the correct way to compute the trace of a tensor expression involving gamma matrices without worrying too much about this kind of pitfalls?
Confirmed this is still a bug in sympy 1.10. This should be considered a bug.
The `sympy.physics.hep` module is unmaintained right now. If you want to contribute improvements or take over maintenance then go ahead.
OK, I'm pretty sure I have now found the right fix for this. Passes all tests including newly added test for this bug. I will try to create a pull request tomorrow after I get some sleep, run flake8 and black, make sure I have set up a branch in my fork correctly, etc.
You don't need to run black. Since black is not currently used if you run it then it will probably rewrite everything so it is best not to do that.

Commit message: Add myself to .mailmap

Patch to review:
diff --git a/.mailmap b/.mailmap
index e832c7be1ef2..03eca177f6ff 100644
--- a/.mailmap
+++ b/.mailmap
@@ -588,6 +588,8 @@ Gilbert Gede <gilbertgede@gmail.com> <ggede@ucdavis.edu>
 Gilles Schintgen <gschintgen@hambier.lu>
 Gina <Dr-G@users.noreply.github.com>
 Gleb Siroki <g.shiroki@gmail.com>
+Glenn Horton-Smith <glenn.hortonsmith@gmail.com>
+Glenn Horton-Smith <glenn.hortonsmith@gmail.com> Glenn Horton-Smith <gahs@phys.ksu.edu>
 GolimarOurHero <metalera94@hotmail.com>
 Goutham Lakshminarayan <dl.goutham@gmail.com> Goutham <devnull@localhost>
 Govind Sahai <gsiitbhu@gmail.com>
diff --git a/sympy/physics/hep/gamma_matrices.py b/sympy/physics/hep/gamma_matrices.py
index 23284e4a9297..3009550888f7 100644
--- a/sympy/physics/hep/gamma_matrices.py
+++ b/sympy/physics/hep/gamma_matrices.py
@@ -190,7 +190,7 @@ def gamma_trace(t):
 
     """
     if isinstance(t, TensAdd):
-        res = TensAdd(*[_trace_single_line(x) for x in t.args])
+        res = TensAdd(*[gamma_trace(x) for x in t.args])
         return res
     t = _simplify_single_line(t)
     res = _trace_single_line(t)
diff --git a/sympy/physics/hep/tests/test_gamma_matrices.py b/sympy/physics/hep/tests/test_gamma_matrices.py
index 27509803f24f..47e4abbac664 100644
--- a/sympy/physics/hep/tests/test_gamma_matrices.py
+++ b/sympy/physics/hep/tests/test_gamma_matrices.py
@@ -3,6 +3,7 @@
     TensExpr, canon_bp
 from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex, \
     kahane_simplify, gamma_trace, _simplify_single_line, simplify_gamma_expression
+from sympy import Symbol
 
 
 def _is_tensor_eq(arg1, arg2):
@@ -399,3 +400,25 @@ def test_gamma_matrix_trace():
     t = ps*ps*ps*ps*ps*ps*ps*ps
     r = gamma_trace(t)
     assert r.equals(4*p2*p2*p2*p2)
+
+def test_bug_13636():
+    """Test issue 13636 regarding handling traces of sums of products
+    of GammaMatrix mixed with other factors."""
+    pi, ki, pf = tensor_heads("pi, ki, pf", [LorentzIndex])
+    i0, i1, i2, i3, i4 = tensor_indices("i0:5", LorentzIndex)
+    x = Symbol("x")
+    pis = pi(i2) * G(-i2)
+    kis = ki(i3) * G(-i3)
+    pfs = pf(i4) * G(-i4)
+
+    a = pfs * G(i0) * kis * G(i1) * pis * G(-i1) * kis * G(-i0)
+    b = pfs * G(i0) * kis * G(i1) * pis * x * G(-i0) * pi(-i1)
+    ta = gamma_trace(a)
+    tb = gamma_trace(b)
+    t_a_plus_b = gamma_trace(a + b)
+    assert ta.equals(
+        -16 * ki(i0) * ki(-i0) * pf(i1) * pi(-i1)
+        + 32 * ki(i0) * ki(i1) * pf(-i0) * pi(-i1)
+    )
+    assert tb.equals(-8 * x * ki(i0) * pf(-i0) * pi(i1) * pi(-i1))
+    assert t_a_plus_b.equals(ta + tb)


Code review comments:
diff_hunk: @@ -399,3 +400,25 @@ def test_gamma_matrix_trace():
     t = ps*ps*ps*ps*ps*ps*ps*ps
     r = gamma_trace(t)
     assert r.equals(4*p2*p2*p2*p2)
+
+def test_bug_13636():
+    """Test issue 13636 regarding handling traces of sums of products
+    of GammaMatrix mixed with other factors."""
+    pi, ki, pf = tensor_heads("pi, ki, pf", [LorentzIndex])
+    i0, i1, i2, i3, i4 = tensor_indices("i0:5", LorentzIndex)
+    x = Symbol("x")
+    pis = pi(i2) * G(-i2)
+    kis = ki(i3) * G(-i3)
+    pfs = pf(i4) * G(-i4)
+
+    a = pfs * G(i0) * kis * G(i1) * pis * G(-i1) * kis * G(-i0)
+    b = pfs * G(i0) * kis * G(i1) * pis * x * G(-i0) * pi(-i1)
+    ta = gamma_trace(a)
+    tb = gamma_trace(b)
+    t_a_plus_b = gamma_trace(a + b)
+    assert ta.equals(
+        -16 * ki(i0) * ki(-i0) * pf(i1) * pi(-i1)
+        + 32 * ki(i0) * ki(i1) * pf(-i0) * pi(-i1)
+    )
+    assert tb.equals(-8 * x * ki(i0) * pf(-i0) * pi(i1) * pi(-i1))
+    assert t_a_plus_b.equals(ta + tb)
path: sympy/physics/hep/tests/test_gamma_matrices.py
text: @user1:
Does `==` not work here or does it need to use `.equals`?

@author:
Good question. I just checked, and indeed == does not work here. The .equals() function for tensor objects applies sympify and .cannon_bp to the tensor expression before comparing, while ._eq_ is just Basic._eq_, and is false for this comparison because the expression trees are not exactly identical. The other test functions such as test_gamma_matrix_trace() also use .equals().


@user1:
Usually it is better to construct the exact expression for comparison in tests so that they can be compared with `==`. Sometimes that is not be best option because the exact expression is confusing and does not display the intent of the test well.

I'm not sure what the preferred answer would be here. What exactly is the difference between the compared expressions? You can use `srepr` as one way to view the expression structure explicitly.

@author:
ta == `4*((-4)*ki(L_0)*ki(-L_0)*pf(L_1)*pi(-L_1) + 8*ki(L_0)*ki(L_1)*pf(-L_0)*pi(-L_1))` which is the same as the expression being tested against except for distribution of the factor of 4 over the coefficients. I can use the exact expression for the test, I don't mind.

tb and ta_plus_tb are the same as the expressions they are being compared against, so we could use == there with no problem.

@author:
This is fixed now.