adapt to "symbols are not callable" change in sympy-1.3

pySecDec/code_writer/make_package.py

@@ -33,6 +33,9 @@
 _sympy_zero = sp.sympify(0)
 _sympy_one = sp.sympify(1)
 
+# sympy symbols are no longer callable starting from version 1.3
+_to_function = lambda x: sp.Function(str(x))
+
 # define the internal names to be used in FORM
 internal_prefix = 'SecDecInternal'
 FORM_names = dict(
@@ -1915,7 +1918,7 @@ def secondary_decomposition_with_splitting(sector, indices, split_sectors=None):
         # remove `polynomial_names` from the `remainder_expression`
         this_primary_sector_remainder_expression = remainder_expression
         for poly_name in reversed_polynomial_names:
-            this_primary_sector_remainder_expression = this_primary_sector_remainder_expression.replace(-1, poly_name(*symbols_polynomials_to_decompose), remove=True)
+            this_primary_sector_remainder_expression = this_primary_sector_remainder_expression.replace(-1, _to_function(poly_name)(*symbols_polynomials_to_decompose), remove=True)
 
         # If there is a nontrivial primary decomposition, remove the integration variables from the `remainder_expression`
         for i,var in enumerate(all_integration_variables):

pySecDec/decomposition/test_common.py

@@ -109,7 +109,7 @@ def setUp(self):
 
     #@attr('active')
     def test_sector2array(self):
-        SecDecInternalCast, SecDecInternalOther = sp.symbols('SecDecInternalCast SecDecInternalOther')
+        SecDecInternalCast, SecDecInternalOther = sp.symbols('SecDecInternalCast SecDecInternalOther', cls=sp.Function)
 
         combined_expolists, combined_coeffs = _sector2array(self.sector)
 
@@ -139,7 +139,7 @@ def test_sector2array(self):
 
     #@attr('active')
     def test_sector2array_cancelling(self):
-        SecDecInternalCast = sp.symbols('SecDecInternalCast')
+        SecDecInternalCast = sp.symbols('SecDecInternalCast', cls=sp.Function)
         a = sp.symbols('a')
 
         mono = Polynomial([(0,1)],[1])
@@ -167,8 +167,8 @@ def test_sector2array_cancelling(self):
     #@attr('active')
     def test_sector2array_other_exponent(self):
 
-        SecDecInternalCast, SecDecInternalOther = sp.symbols('SecDecInternalCast SecDecInternalOther')
-        SecDecInternalExponent = sp.symbols('SecDecInternalExponent')
+        SecDecInternalCast, SecDecInternalOther = sp.symbols('SecDecInternalCast SecDecInternalOther', cls=sp.Function)
+        SecDecInternalExponent = sp.symbols('SecDecInternalExponent', cls=sp.Function)
 
         combined_expolists, combined_coeffs = _sector2array(self.sector_other_exponent)
 

pySecDec/loop_integral/from_propagators.py

@@ -7,6 +7,9 @@
 import sympy as sp
 import numpy as np
 
+# sympy symbols are no longer callable starting from version 1.3
+_to_function = lambda x: sp.Function(str(x))
+
 class LoopIntegralFromPropagators(LoopIntegral):
     __doc__ = r'''
     Construct the Feynman parametrization of a
@@ -278,8 +281,8 @@ def _numerator_tensors(self):
                     for i,momentum in enumerate(momenta):
 
                         # search for ``momentum(index)``
-                        while term.subs(momentum(index), 0) == 0: # expression has a factor `momentum(index)`
-                            term /= momentum(index)
+                        while term.subs(_to_function(momentum)(index), 0) == 0: # expression has a factor `momentum(index)`
+                            term /= _to_function(momentum)(index)
                             lst.append((i,index))
                             index_count += 1
 
@@ -326,7 +329,7 @@ def preliminary_numerator(self):
 
         aM = self.aM
         Q = self.Q
-        g = self.metric_tensor
+        g = _to_function(self.metric_tensor)
         D = self.dimensionality
         L = self.L
         Feynman_parameters_U_F = self.Feynman_parameters + sp.sympify(['U', 'F'])
@@ -498,7 +501,7 @@ def reduce_gamma(r_over_two):
                         this_tensor_P_factor = aM[external_momentum_index].dot(Q)
                         for j, coeff in enumerate(this_tensor_P_factor.coeffs):
                             if isinstance(coeff, sp.Expr):
-                                this_tensor_P_factor.coeffs[j] = coeff.subs((p, p(Lorentz_index_P)) for p in self.external_momenta)
+                                this_tensor_P_factor.coeffs[j] = coeff.subs((p, _to_function(p)(Lorentz_index_P)) for p in self.external_momenta)
 
                         # must append ``F`` and ``U`` to the parameters of ``this_tensor_P_factor``
                         this_tensor_P_factor.expolist = np.hstack([this_tensor_P_factor.expolist, np.zeros((len(this_tensor_P_factor.expolist), 2), dtype=int)])
@@ -518,14 +521,14 @@ def reduce_gamma(r_over_two):
                             if metric_unmatched:
                                 if Lorentz_index_1 == Lorentz_index_external:
                                     contractions_unmatched[i] = external_momentum_unmatched = False
-                                    this_numerator_summand *= self.external_momenta[external_momentum_index](Lorentz_index_2)
+                                    this_numerator_summand *= _to_function(self.external_momenta[external_momentum_index])(Lorentz_index_2)
                                     break
                                 if Lorentz_index_2 == Lorentz_index_external:
                                     contractions_unmatched[i] = external_momentum_unmatched = False
-                                    this_numerator_summand *= self.external_momenta[external_momentum_index](Lorentz_index_1)
+                                    this_numerator_summand *= _to_function(self.external_momenta[external_momentum_index])(Lorentz_index_1)
                                     break
                         if external_momentum_unmatched:
-                            this_numerator_summand *= self.external_momenta[external_momentum_index](Lorentz_index_external)
+                            this_numerator_summand *= _to_function(self.external_momenta[external_momentum_index])(Lorentz_index_external)
 
                     # ------------------------------------------------------------------------------
 
@@ -558,8 +561,8 @@ def replacement_rules_with_Lorentz_indices(self):
                 pattern_with_index = pattern
                 replacement_with_index = replacement
                 for p in self.external_momenta:
-                    pattern_with_index = pattern_with_index.subs(p, p(mu))
-                    replacement_with_index = replacement_with_index.subs(p, p(mu))
+                    pattern_with_index = pattern_with_index.subs(p, _to_function(p)(mu))
+                    replacement_with_index = replacement_with_index.subs(p, _to_function(p)(mu))
                 replacement_rules.append( (pattern_with_index, replacement_with_index) )
 
         return replacement_rules

pySecDec/loop_integral/test_loop_integral.py

@@ -614,7 +614,8 @@ def test_2L_box_with_numerator(self):
     #@attr('active')
     @attr('slow')
     def test_rank3_numerator_2L_double_triangle(self):
-        k1, k2, p1, p2, p3, p4, s, t, mu, nu, scalar_factor, eps, F = sp.symbols('k1 k2 p1 p2 p3 p4 s t mu nu scalar_factor eps F')
+        k1,  k2,  p1,  p2,  p3,  p4, s, t, mu, nu, scalar_factor, eps, F = sp.symbols('k1 k2 p1 p2 p3 p4 s t mu nu scalar_factor eps F')
+        k1f, k2f, p1f, p2f, p3f, p4f = sp.symbols('k1 k2 p1 p2 p3 p4', cls=sp.Function)
         sp11, sp12, sp13, sp22, sp23, sp33 = sp.symbols('sp11 sp12 sp13 sp22 sp23 sp33')
 
         z = sp.sympify(['z%i'%i for i in range(6+1)]) # Feynman parameters (only use z[1], z[2], ..., z[6] but not z[0])
@@ -630,7 +631,7 @@ def test_rank3_numerator_2L_double_triangle(self):
                              ('p1*p3', 'sp13')]
 
         propagators = [k1**2,(k1+p1)**2,(k1-k2-p3-p4)**2,(k2)**2,(k2+p2)**2,(k2+p3)**2]
-        numerator = k1(mu)*k2(nu)*k1(mu)*p2(nu)*2
+        numerator = k1f(mu)*k2f(nu)*k1f(mu)*p2f(nu)*2
 
         li = LoopIntegralFromPropagators(propagators, loop_momenta, external_momenta,
                                            numerator=numerator, Feynman_parameters=z[1:],
@@ -720,16 +721,17 @@ def test_rank3_numerator_2L_double_triangle(self):
     #@attr('active')
     @attr('slow')
     def test_rank2_numerator_2L_box(self):
-        k1, k2, p1, p2, p3, p4, s, t, mu, scalar_factor, D, eps = sp.symbols('k1 k2 p1 p2 p3 p4 s t mu scalar_factor D eps')
+        k1,  k2,  p1,  p2,  p3,  p4, s, t, mu, scalar_factor, D, eps = sp.symbols('k1 k2 p1 p2 p3 p4 s t mu scalar_factor D eps')
+        k1f, k2f, p1f, p2f, p3f, p4f = sp.symbols('k1 k2 p1 p2 p3 p4', cls=sp.Function)
 
         z = sp.sympify(['z%i'%i for i in range(7+1)]) # Feynman parameters (only use z[1], z[2], ..., z[7] but not z[0])
         loop_momenta = [k1, k2]
         external_momenta = [p1, p2, p3, p4]
         indices = [mu]
 
         propagators = [k1**2,(k1+p2)**2,(k1-p1)**2,(k1-k2)**2,(k2+p2)**2,(k2-p1)**2,(k2+p2+p3)**2]
-        numerator1 = 2*k1(mu)*k2(mu)
-        numerator2 = 2*k1(mu)*k1(mu)
+        numerator1 = 2*k1f(mu)*k2f(mu)
+        numerator2 = 2*k1f(mu)*k1f(mu)
         replacement_rules = [(p1**2, 0),
                              (p2**2, 0),
                              (p3**2, 0),
@@ -878,7 +880,7 @@ def test_rank2_numerator_2L_box(self):
             Feynman_parametrized_numerator = - Feynman_parametrized_numerator
 
             # only one scalar product per term --> can safely remove dummy index
-            Feynman_parametrized_numerator = Feynman_parametrized_numerator.subs(p1(mu), p1).subs(p2(mu), p2).subs(p3(mu), p3).subs(p4(mu), p4)
+            Feynman_parametrized_numerator = Feynman_parametrized_numerator.subs(p1f(mu), p1).subs(p2f(mu), p2).subs(p3f(mu), p3).subs(p4f(mu), p4)
 
             if i == 1:
                 Feynman_parametrized_numerator = Feynman_parametrized_numerator.subs(D, 4-2*eps)

pySecDec/test_algebra.py

@@ -1350,7 +1350,8 @@ def test_polynomial(self):
 
     #@attr('active')
     def test_function(self):
-        a, x, my_function = sp.symbols('a x my_function')
+        a, x = sp.symbols('a x')
+        my_function = sp.Function('my_function')
         sympy_function = my_function(a*x)
         function = Expression(sympy_function, ['x'])