Merge branch 'selective_ibp'

high_level_tests/selective_ibp/Makefile

@@ -0,0 +1,29 @@
+include ./Makefile.conf
+include ../test_suite/Makefile.conf
+
+$(INTEGRAL_NAME) : generate_$(INTEGRAL_NAME).py
+	$(PYTHON) $<
+
+.PHONY : static-library
+static-library : $(INTEGRAL_NAME)
+	$(MAKE) -C $(INTEGRAL_NAME) static
+
+.PHONY : pylink-library
+pylink-library : $(INTEGRAL_NAME) static-library
+	$(MAKE) -C $(INTEGRAL_NAME) pylink
+
+.PHONY : test
+test : static-library pylink-library
+	echo $(EXAMPLES_DIR)
+	$(MAKE) -C test INTEGRAL_NAME=$(INTEGRAL_NAME) $@
+
+.PHONY : clean
+clean :
+	$(MAKE) -C $(INTEGRAL_NAME) $@
+	$(MAKE) -C test INTEGRAL_NAME=$(INTEGRAL_NAME) $@
+
+.PHONY : very-clean
+very-clean :
+	rm -rf $(INTEGRAL_NAME)
+	$(MAKE) -C test INTEGRAL_NAME=$(INTEGRAL_NAME) $@
+

high_level_tests/selective_ibp/Makefile.conf

@@ -0,0 +1 @@
+INTEGRAL_NAME = userdefined_cpp

high_level_tests/selective_ibp/functions_implementation.hpp

@@ -0,0 +1,53 @@
+#ifndef userdefined_cpp_src_functions_hpp_included
+#define userdefined_cpp_src_functions_hpp_included
+
+#include "userdefined_cpp.hpp"
+#include "SecDecInternalFunctions.hpp"
+
+#include <cmath>
+#include <complex>
+
+namespace userdefined_cpp
+{
+
+    /*
+     * Declarations of the `functions` and their required
+     * derivatives are declared here. The derivative of a function
+     * 'f' by its i-th argument is denoted as 'dfdi'. To implement
+     * the functions listed below, you can either add "inline"
+     * keywords and define these functions here, or you define the
+     * functions in a separate '.cpp' file. If you decide for a
+     * separate file, the file name can be arbitrary up to the
+     * '.cpp' suffix. Furthermore, the '.cpp' file must be located
+     * in this directory ('src/'). If you want to link against
+     * an external library (e.g. the gsl), you should add the
+     * corresponding compiler and linker flags to the "Makefile.conf"
+     * in the top level directory.
+     *
+     * Note: Not all functions listed here may actually be needed.
+     *       This file lists all derivatives that occurred in the
+     *       calculation. It is possible that some dropped out due
+     *       to algebraic simplifications after this list was
+     *       generated.
+     */
+
+    template<typename T>
+    integrand_return_t func(T arg)
+    {
+        return 1.0;
+    }
+
+    template<typename T>
+    integrand_return_t dfuncd0(T arg)
+    {
+        return 0.0;
+    }
+
+    inline
+    integrand_return_t HeavisideTheta(real_t x)
+    {
+        return x > 0.0 ? 1.0 : 0.0;
+    }
+
+};
+#endif

high_level_tests/selective_ibp/functions_template.hpp

@@ -0,0 +1,44 @@
+#ifndef userdefined_cpp_src_functions_hpp_included
+#define userdefined_cpp_src_functions_hpp_included
+
+#include "userdefined_cpp.hpp"
+#include "SecDecInternalFunctions.hpp"
+
+#include <cmath>
+#include <complex>
+
+namespace userdefined_cpp
+{
+
+    /*
+     * Declarations of the `functions` and their required
+     * derivatives are declared here. The derivative of a function
+     * 'f' by its i-th argument is denoted as 'dfdi'. To implement
+     * the functions listed below, you can either add "inline"
+     * keywords and define these functions here, or you define the
+     * functions in a separate '.cpp' file. If you decide for a
+     * separate file, the file name can be arbitrary up to the
+     * '.cpp' suffix. Furthermore, the '.cpp' file must be located
+     * in this directory ('src/'). If you want to link against
+     * an external library (e.g. the gsl), you should add the
+     * corresponding compiler and linker flags to the "Makefile.conf"
+     * in the top level directory.
+     *
+     * Note: Not all functions listed here may actually be needed.
+     *       This file lists all derivatives that occurred in the
+     *       calculation. It is possible that some dropped out due
+     *       to algebraic simplifications after this list was
+     *       generated.
+     */
+
+/*
+%s
+
+%s
+
+%s
+
+*/
+
+};
+#endif

high_level_tests/selective_ibp/generate_userdefined_cpp.py

@@ -0,0 +1,61 @@
+#! /usr/bin/env python
+import shutil
+from itertools import permutations
+from pySecDec import make_package
+
+name = 'userdefined_cpp'
+
+make_package(
+
+name=name,
+integration_variables = ['x','y','z'],
+
+# use integration by parts rather than subtraction for ``x`` and ``y``
+# generate subtraction term for ``z`` to avoid derivative of HeavisideTheta
+ibp_power_goal = [0,0,-1],
+
+# the order here defines the order of the expansion
+regulators = ['eps'],
+functions = ['HeavisideTheta','func'],
+
+polynomials_to_decompose = ['(x+y)**(-2+eps)','z**(-1+eps)'],
+remainder_expression = 'HeavisideTheta(1/4-z)*func(y)',
+# full analytic result for func(y)=1:
+#   (4^-eps (-2 + 2^eps))/((-1 + eps) eps^2)
+# = 1/eps^2 + (1 - Log[2] - Log[4])/eps +  1/2 (2 - 2 Log[2] - Log[2]^2 - 2 Log[4] + 2 Log[2] Log[4] + Log[4]^2) + O[eps]
+# = 1.0000000000000000000/eps^2 - 1.0794415416798359283/eps + 0.60214400703386905808 + O[eps]
+
+# the highest orders of the final regulator expansion
+# the order here matches the order of ``regulators``
+requested_orders = [0],
+
+)
+
+# check the generated "functions.hpp"
+cpp_function_declarations = (
+'''\
+    template<typename T0>
+    integrand_return_t HeavisideTheta(T0 arg0);\
+''',
+'''\
+    template<typename T0>
+    integrand_return_t dfuncd0(T0 arg0);\
+''',
+'''\
+    template<typename T0>
+    integrand_return_t func(T0 arg0);\
+'''
+)
+with open(name + '/src/functions.hpp') as generated_header_file:
+    generated_header = generated_header_file.read()
+with open('functions_template.hpp') as target_header_template_file:
+    target_header_template = target_header_template_file.read()
+matched_header = False
+for ordering in permutations(cpp_function_declarations):
+    if generated_header == target_header_template % ordering:
+        matched_header = True
+        break
+assert matched_header, 'mismatch between generated and expected "functions.hpp"'
+
+# copy 'functions.hpp' (predefined for this example) to required directory
+shutil.copy('functions_implementation.hpp',name+'/src/functions.hpp')

high_level_tests/selective_ibp/test/Makefile

@@ -0,0 +1,15 @@
+include ../Makefile.conf
+include ../../test_suite/Makefile.conf
+
+.PHONY : test
+test : test.log
+test.log : test.py
+	(echo "PYTHON test:" && $(PYTHON) test.py -v 2>&1) 2>&1 >test.log \
+	&& printf "\n@@@ SUCCESS @@@" >> test.log || printf "\n@@@ FAILURE @@@" >> test.log
+
+.PHONY : clean very-clean
+clean :
+	rm -rf *.o *.exe *.log
+
+very-clean  : clean
+

high_level_tests/selective_ibp/test/test.py

@@ -0,0 +1,81 @@
+from __future__ import print_function
+from pySecDec.integral_interface import IntegralLibrary
+import sympy as sp
+import unittest
+
+class CheckLib(unittest.TestCase):
+    def setUp(self):
+        # load c++ library
+        self.lib = IntegralLibrary('../userdefined_cpp/userdefined_cpp_pylink.so')
+
+        # full analytic result for func(y)=1:
+        #   (4^-eps (-2 + 2^eps))/((-1 + eps) eps^2)
+        # = 1/eps^2 + (1 - Log[2] - Log[4])/eps +  1/2 (2 - 2 Log[2] - Log[2]^2 - 2 Log[4] + 2 Log[2] Log[4] + Log[4]^2) + O[eps]
+        # = 1.0000000000000000000/eps^2 - 1.0794415416798359283/eps + 0.60214400703386905808 + O[eps]
+        self.target_result = {
+                                 -2:   1.0,
+                                 -1:  -1.0794415416798359283,
+                                  0:   0.60214400703386905808
+                             }
+        self.epsrel = 1e-4
+        self.maxeval = 10**8
+
+    def check_result(self, computed_series, epsrel):
+        # convert result to sympy expressions
+        integral_with_prefactor = sp.sympify(  computed_series.replace(',','+I*').replace('+/-','*value+error*')  )
+
+        for order in range(-2,0):
+            value = integral_with_prefactor.coeff('eps',order).coeff('value')
+            error = integral_with_prefactor.coeff('eps',order).coeff('error')
+
+            # check that the uncertainties are reasonable
+            self.assertLessEqual(error, abs(4*epsrel * self.target_result[order]))
+
+            # check that the desired uncertainties are reached
+            self.assertLessEqual(error, abs(epsrel * value) )
+
+            # check integral value
+            self.assertAlmostEqual(  value, self.target_result[order], delta=epsrel*abs(self.target_result[order])  )
+
+    def test_Vegas(self):
+        # choose integrator
+        self.lib.use_Vegas(flags=2, epsrel=self.epsrel, maxeval=self.maxeval) # ``flags=2``: verbose --> see Cuba manual
+
+        # integrate
+        str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor = self.lib()
+
+        # check
+        self.check_result(str_integral_with_prefactor, self.epsrel)
+
+    def test_Suave(self):
+        # choose integrator
+        self.lib.use_Suave(flags=2, epsrel=self.epsrel, maxeval=self.maxeval, nnew=10**5) # ``flags=2``: verbose --> see Cuba manual
+
+        # integrate
+        str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor = self.lib()
+
+        # check
+        self.check_result(str_integral_with_prefactor, self.epsrel)
+
+    def test_Divonne(self):
+        # choose integrator
+        self.lib.use_Divonne(flags=2, epsrel=self.epsrel, border=1e-12, maxeval=self.maxeval) # ``flags=2``: verbose --> see Cuba manual
+
+        # integrate
+        str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor = self.lib()
+
+        # check
+        self.check_result(str_integral_with_prefactor, self.epsrel)
+
+    def test_Cuhre(self):
+        # choose integrator
+        self.lib.use_Cuhre(flags=2, epsrel=self.epsrel, maxeval=self.maxeval) # ``flags=2``: verbose --> see Cuba manual
+
+        # integrate
+        str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor = self.lib()
+
+        # check
+        self.check_result(str_integral_with_prefactor, self.epsrel)
+
+if __name__ == '__main__':
+    unittest.main()