Merge branch 'fix-Pak'

pySecDec/code_writer/make_package.py

@@ -11,7 +11,7 @@
                       ExponentiatedPolynomial, Pow, Product, \
                       ProductRule, Function, Sum
 from .. import decomposition
-from ..matrix_sort import iterative_sort, Pak_sort
+from ..matrix_sort import iterative_sort, Pak_sort, light_Pak_sort
 from ..subtraction import integrate_pole_part, integrate_by_parts, pole_structure as compute_pole_structure
 from ..expansion import expand_singular, expand_Taylor, expand_sympy, OrderError
 from ..misc import lowest_order, det
@@ -739,7 +739,8 @@ def make_package(name, integration_variables, regulators, requested_orders,
                  complex_parameters=[], form_optimization_level=2, form_work_space='500M',
                  form_insertion_depth=5, contour_deformation_polynomial=None, positive_polynomials=[],
                  decomposition_method='iterative_no_primary', normaliz_executable='normaliz',
-                 enforce_complex=False, split=False, ibp_power_goal=-1, use_dreadnaut=True):
+                 enforce_complex=False, split=False, ibp_power_goal=-1, use_dreadnaut=False,
+                 use_Pak=True):
     r'''
     Decompose, subtract and expand an expression.
     Return it as c++ package.
@@ -968,6 +969,13 @@ def make_package(name, integration_variables, regulators, requested_orders,
         ``$SECDEC_CONTRIB/bin/dreadnaut`` and, if the
         environment variable ``$SECDEC_CONTRIB`` is not set,
         ``dreadnaut``.
+        Default: ``False``
+
+    :param use_Pak:
+        bool;
+        Whether or not to use
+        :func:`.squash_symmetry_redundant_sectors_sort`
+        with :func:`.Pak_sort` to find sector symmetries.
         Default: ``True``
 
     '''
@@ -1054,11 +1062,14 @@ def reduce_sectors_by_symmetries(sectors, message):
         # find symmetries
         sectors = decomposition.squash_symmetry_redundant_sectors_sort(sectors, iterative_sort)
         print(message + ' after symmetry finding (iterative):', len(sectors))
-        sectors = decomposition.squash_symmetry_redundant_sectors_sort(sectors, Pak_sort)
-        print(message + ' after symmetry finding (iterative+Pak):', len(sectors))
+        sectors = decomposition.squash_symmetry_redundant_sectors_sort(sectors, light_Pak_sort)
+        print(message + ' after symmetry finding (light Pak):', len(sectors))
+        if use_Pak:
+            sectors = decomposition.squash_symmetry_redundant_sectors_sort(sectors, Pak_sort)
+            print(message + ' after symmetry finding (full Pak):', len(sectors))
         if use_dreadnaut:
             sectors = decomposition.squash_symmetry_redundant_sectors_dreadnaut(sectors, dreadnaut_executable, os.path.join(name,'dreadnaut_workdir'))
-            print(message + ' after symmetry finding (iterative+Pak+dreadnaut):', len(sectors))
+            print(message + ' after symmetry finding (dreadnaut):', len(sectors))
         return sectors
 
     # if splitting desired, implement it as additional primary decomposition

pySecDec/decomposition/common.py

@@ -5,7 +5,6 @@
 '''
 
 from ..algebra import Polynomial, ExponentiatedPolynomial, Product
-from ..matrix_sort import iterative_sort, Pak_sort
 from ..misc import argsort_ND_array
 import numpy as np
 import sympy as sp
@@ -606,7 +605,8 @@ def squash_symmetry_redundant_sectors_sort(sectors, sort_function):
         reduced.
 
     :param sort_function:
-        :func:`pySecDec.matrix_sort.iterative_sort` or
+        :func:`pySecDec.matrix_sort.iterative_sort`,
+        :func:`pySecDec.matrix_sort.light_Pak_sort`, or
         :func:`pySecDec.matrix_sort.Pak_sort`;
         The function to be used for finding a canonical
         form of the sectors.

pySecDec/decomposition/test_common.py

@@ -3,7 +3,7 @@
 from .common import *
 from .common import _sector2array, _collision_safe_hash
 from ..algebra import Polynomial, ExponentiatedPolynomial, Product
-from ..matrix_sort import iterative_sort, Pak_sort
+from ..matrix_sort import iterative_sort, Pak_sort, light_Pak_sort
 import unittest
 import sympy as sp
 import numpy as np
@@ -99,19 +99,6 @@ def setUp(self):
         self.sector_p0 = Sector([self.p0], Jacobian=self.Jacobian)
         self.sector_swapped_p0 = Sector([self.swapped_p0], Jacobian=self.swapped_Jacobian)
 
-        np.random.seed(0)
-        # Aside: Pak_sort passes with seed(10),low=0,high=3, size=(101,7)
-        self.p1_hard_expolist = np.random.randint(low=0,high=4, size=(501,7))
-        self.p1_hard_expolist_permuted=np.random.permutation(
-            np.transpose(np.random.permutation(np.transpose(self.p1_hard_expolist))))
-
-        self.Jacobian_hard = Polynomial([(1,1,1,1,1,1,1)], ['a'])
-        self.p1_hard = Polynomial(self.p1_hard_expolist, [1]*501)
-        self.sector_p1_hard = Sector([self.p1_hard],[], self.Jacobian_hard)
-
-        self.swapped_p1_hard = Polynomial(self.p1_hard_expolist_permuted, [1]*501)
-        self.sector_swapped_p1_hard = Sector([self.swapped_p1_hard],[], self.Jacobian_hard)
-
         self.a, self.b, self.c, self.d, self.e, self.f, self.g = sp.symbols('a b c d e f g')
 
     #@attr('active')
@@ -198,20 +185,69 @@ def test_squash_symmetry_redundant_sectors_2D(self):
         self.assertEqual((sp.sympify(reduced_sectors[0].cast[0]) - sp.sympify(self.p0.copy())).simplify(), 0)
 
     #@attr('active')
+    @attr('slow')
     def test_squash_symmetry_hard(self):
-        sectors = [self.sector_p1_hard.copy(), self.sector_swapped_p1_hard.copy()]
 
-        # test symmetry finding by sorting, fails
-        for sort_function in (iterative_sort, Pak_sort):
-            reduced_sectors = squash_symmetry_redundant_sectors_sort(sectors, sort_function)
-            self.assertNotEqual(len(reduced_sectors), 1)
-            self.assertNotEqual(reduced_sectors[0].Jacobian.coeffs[0], sp.sympify('2*a'))
+        # In python3, the hashes of symbols are random.
+        #     --> Select a random ``s`` until the fast
+        #         and incomplete algorithms miss symmetries.
+        too_easy = True
+        while too_easy:
+            # hard example: from "examples/triangle"
+            s = np.random.randint(2**63-1)
+            x0,x1,x2,x3,x4,x5 = symbols_hard = sp.sympify(['x%i'%i for i in range(6)])
+            Jacobian_hard = Polynomial([[1]*len(symbols_hard)], ['a'])
+            hard_p1 = (Polynomial.from_expression( + (1)*x1*x2 + (1)*x3 + (1)*x1*x3 + (1)*x1*x4 + (1)*x1*x5 + (1)*x1 + (1)*x4 + (1)*x2*x3 + (1)*x2*x4
+                                                   + (1)*x2*x5 + (1)*x5, symbols_hard) ** sp.sympify('3*eps')
+                      ).simplify()
+            hard_p2 = (Polynomial.from_expression( + (2*s)*x1*x2*x3 + (s)*x1*x2*x4 + (2*s)*x1*x2*x5 + (s)*x1*x2
+                                                   + (s)*x5 + (-s)*x2*x4*x5 + (s)*x2*x4 + (s)*x2*x5 + (s)*x3 + (-s)*x4*x5
+                                                   + (2*s)*x1*x3 + (s)*x1 + (-s)*x1*x4*x5 + (2*s)*x1*x4 + (s)*x1*x5
+                                                   + (s)*x2**2*x3 + (s)*x2**2*x4 + (s)*x2**2*x5 + (s)*x2*x3
+                                                   + (s)*x1**2*x2 + (s)*x1**2*x3 + (s)*x1**2*x4 + (s)*x1**2*x5 + (s)*x1**2
+                                                   + (s)*x1*x2**2 + (s)*x4, symbols_hard) ** sp.sympify('-2*eps - 2')
+                      ).simplify()
+            hard_p1_permuted = (Polynomial.from_expression( + (1)*x0*x1 + (1)*x4 + (1)*x0*x3 + (1)*x0*x4 + (1)*x5
+                                                            + (1)*x0*x5 + (1)*x1 + (1)*x1*x3 + (1)*x1*x4 + (1)*x1*x5
+                                                            + (1)*x3, symbols_hard) ** sp.sympify('3*eps')
+                               ).simplify()
+            hard_p2_permuted = (Polynomial.from_expression( + (s)*x0**2*x1 + (s)*x0**2*x3 + (s)*x0**2*x4 + (s)*x0**2*x5
+                                                            + (2*s)*x0*x1*x3 + (2*s)*x0*x1*x4 + (s)*x0*x1*x5 + (s)*x0*x1 + (s)*x0*x3
+                                                            + (-s)*x0*x4*x5 + (s)*x0*x4 + (s)*x0*x5 + (s)*x1**2*x3 + (s)*x1**2*x4
+                                                            + (s)*x1**2*x5 + (s)*x1**2 + (2*s)*x1*x3 + (-s)*x1*x4*x5 + (s)*x1*x4
+                                                            + (s)*x5 + (2*s)*x1*x5 + (s)*x1 + (s)*x3 + (-s)*x4*x5 + (s)*x4
+                                                            + (s)*x0*x1**2, symbols_hard) ** sp.sympify('-2*eps - 2')
+                               ).simplify()
+            sector_hard = Sector([hard_p1,hard_p2], [], Jacobian_hard)
+            sector_swapped_hard = Sector([hard_p1_permuted,hard_p2_permuted], [], Jacobian_hard)
+            sectors = [sector_hard.copy(), sector_swapped_hard.copy()]
+
+            try:
+                # test symmetry finding by iterative sorting, fails
+                reduced_sectors = squash_symmetry_redundant_sectors_sort(sectors, iterative_sort)
+                self.assertNotEqual(len(reduced_sectors), 1)
+                self.assertNotEqual(reduced_sectors[0].Jacobian.coeffs[0], sp.sympify('2*a'))
+
+                # test symmetry finding using light sorting, fails
+                reduced_sectors = squash_symmetry_redundant_sectors_sort(sectors, light_Pak_sort)
+                self.assertNotEqual(len(reduced_sectors), 1)
+                self.assertNotEqual(reduced_sectors[0].Jacobian.coeffs[0], sp.sympify('2*a'))
+
+                too_easy = False
+
+            except AssertionError:
+                pass
 
         # test symmetry finding by graph (using dreadnaut)
         reduced_sectors = squash_symmetry_redundant_sectors_dreadnaut(sectors, dreadnaut_executable, workdir='tmpdir_test_squash_symmetry_hard_python' + python_major_version)
         self.assertEqual(len(reduced_sectors), 1)
         self.assertEqual(reduced_sectors[0].Jacobian.coeffs[0], sp.sympify('2*a'))
 
+        # test symmetry finding using Pak's full algorithm
+        reduced_sectors = squash_symmetry_redundant_sectors_sort(sectors, Pak_sort)
+        self.assertEqual(len(reduced_sectors), 1)
+        self.assertEqual(reduced_sectors[0].Jacobian.coeffs[0], sp.sympify('2*a'))
+
     #@attr('active')
     def test_symmetry_4D(self):
         # sectors 0 and 2 are related by permutation, sector 1 is unrelated

pySecDec/loop_integral/loop_package.py

@@ -25,7 +25,8 @@ def loop_package(name, loop_integral, requested_order,
                  normaliz_executable='normaliz',
                  enforce_complex=False,
                  split=False, ibp_power_goal=-1,
-                 use_dreadnaut=True):
+                 use_dreadnaut=False,
+                 use_Pak=True):
     '''
     Decompose, subtract and expand a Feynman
     parametrized loop integral. Return it as
@@ -158,6 +159,13 @@ def loop_package(name, loop_integral, requested_order,
         ``$SECDEC_CONTRIB/bin/dreadnaut`` and, if the
         environment variable ``$SECDEC_CONTRIB`` is not set,
         ``dreadnaut``.
+        Default: ``False``
+
+    :param use_Pak:
+        bool;
+        Whether or not to use
+        :func:`.squash_symmetry_redundant_sectors_sort`
+        with :func:`.Pak_sort` to find sector symmetries.
         Default: ``True``
 
     '''

pySecDec/matrix_sort.py

@@ -21,7 +21,7 @@ def iterative_sort(matrix):
         orderings depending on the initial ordering.
 
     .. seealso::
-        :func:`.Pak_sort`
+        :func:`.Pak_sort`, :func:`.light_Pak_sort`
 
     :param matrix:
         2D array-like;
@@ -49,15 +49,72 @@ def Pak_sort(matrix):
     Inplace modify the `matrix` to some ordering,
     when permutations of rows and columns (excluding
     the first) are allowed. The implementation of
-    this function is describedin chapter 2 of
+    this function is described in chapter 2 of
     [Pak11]_.
 
     .. note::
         This function may result in different
         orderings depending on the initial ordering.
 
     .. seealso::
-        :func:`.iterative_sort`
+        :func:`.iterative_sort`, :func:`.light_Pak_sort`
+
+    :param matrix:
+        2D array-like;
+        The matrix to be canonicalized.
+
+    '''
+    options = [matrix]
+    for i in range(1, matrix.shape[1]):
+        permutations = []
+        for m in options:
+            for j in range(i, m.shape[1]):
+                permuted_matrix = m.copy()
+
+                # permute integration variables `column_to_swap` and `j`
+                permuted_matrix[:, i] = m[:, j]
+                permuted_matrix[:, j] = m[:, i]
+
+                # sort by rows
+                permuted_matrix[:] = permuted_matrix[argsort_2D_array(permuted_matrix)]
+
+                # transpose since we need column-wise ordering in the next step
+                permutations.append(permuted_matrix.T)
+
+        # sort the matrices from smallest to largest
+        sorted_matrix = argsort_ND_array(permutations)
+
+        # add all matrices that have the largest possible value for `i' to list of options to check
+        options = []
+        for k in range(len(sorted_matrix)-1,-1,-1):
+            if np.array_equal(permutations[sorted_matrix[k]][i], permutations[sorted_matrix[-1]][i]):
+                for mat in options:
+                    if np.array_equal(mat, permutations[sorted_matrix[k]].T):
+                        break
+                else:
+                    options.append(permutations[sorted_matrix[k]].T)
+            else:
+                break
+
+    # resulting options are all equivalent, take first
+    matrix[:] = options[0]
+
+def light_Pak_sort(matrix):
+    '''
+    Inplace modify the `matrix` to some ordering,
+    when permutations of rows and columns (excluding
+    the first) are allowed. The implementation of
+    this function is described in chapter 2 of
+    [Pak11]_. This function implements a lightweight
+    version: In step (v), we only consider one, not
+    all table copies with the maximized second column.
+
+    .. note::
+        This function may result in different
+        orderings depending on the initial ordering.
+
+    .. seealso::
+        :func:`.iterative_sort`, :func:`.Pak_sort`
 
     :param matrix:
         2D array-like;
@@ -66,7 +123,7 @@ def Pak_sort(matrix):
     '''
     for i in range(1,matrix.shape[1]):
         # sort all permutations of columns `i` and `j` (where `i`<=`j`) --> pick largest
-        permutaions = []
+        permutations = []
         for j in range(i,matrix.shape[1]):
             permuted_matrix = matrix.copy()
 
@@ -75,11 +132,11 @@ def Pak_sort(matrix):
             permuted_matrix[:,j] = matrix[:,i]
 
             # sort by rows
-            permuted_matrix[:] = permuted_matrix[argsort_2D_array(matrix)]
+            permuted_matrix[:] = permuted_matrix[argsort_2D_array(permuted_matrix)]
 
             # transpose since we need column-wise ordering in the next step
-            permutaions.append(permuted_matrix.T)
+            permutations.append(permuted_matrix.T)
 
-        # find the largest `i`th column in `permutaions` and keep that
-        index_of_largest = argsort_ND_array(permutaions)[-1]
-        matrix[:] = permutaions[index_of_largest].T # transpose back
+        # find the largest `i`th column in `permutations` and keep that
+        index_of_largest = argsort_ND_array(permutations)[-1]
+        matrix[:] = permutations[index_of_largest].T # transpose back

pySecDec/test_matrix_sort.py

@@ -51,7 +51,7 @@ def test_iterative_algorithm(self):
         self.check_algorithm(iterative_sort, self.matrix_3_by_4, expected_solutions)
 
     #@attr('active')
-    def test_Pak_algorithm(self):
+    def test_light_Pak_algorithm(self):
         expected_solutions = [
                                  [[2, 0, 2, 0],
                                   [2, 1, 2, 1],
@@ -62,3 +62,12 @@ def test_Pak_algorithm(self):
                                   [2, 2, 1, 1]]
                              ]
         self.check_algorithm(Pak_sort, self.matrix_3_by_4, expected_solutions)
+
+    #@attr('active')
+    def test_Pak_algorithm(self):
+        expected_solutions = [
+                                 [[2, 0, 2, 1],
+                                  [2, 2, 0, 0],
+                                  [2, 2, 1, 1]]
+                             ]
+        self.check_algorithm(Pak_sort, self.matrix_3_by_4, expected_solutions)