pySecDec¶
pySecDec [PSD17], [PSD18] is a toolbox for the calculation of dimensionally regulated parameter integrals using the sector decomposition approach [BH00]; see also [Hei08], [BHJ+15].
 Please cite the following references if you use pySecDec for a scientific publication:
Installation¶
Download the Program and Install¶
pySecDec should run fine with both, python 2.7 and python 3 on unixlike systems.
Before you install pySecDec, make sure that you have recent versions of numpy (http://www.numpy.org/) and sympy (http://www.sympy.org/) installed. The version of sympy should be 0.7.6 or higher, the version of numpy should be 1.6 or higher. Type
$ python c "import numpy"
$ python c "import sympy"
to check for their availability.
In case either numpy or sympy are missing on your machine, it is easiest to install them from your package repository. Alternatively, and in particular if you do not have administrator rights, pip (https://pip.pypa.io/en/stable/) may be used to perform the installation.
To install pySecDec download and upack the tarball from http://secdec.hepforge.org/. The tarball contains a distribution of pySecDec and the additional dependencies listed below. Typing
$ make
should build all redistributed packages and display two commands
to be added to your .bashrc
or .profile
.
The Geomethod and Normaliz¶
Note
If you are not urgently interested in using the
geometric decomposition
, you
can ignore this section for the beginning. The instructions below are
not essential for a pySecDec installation. You can still install
normaliz after installing pySecDec. All but the
geometric decomposition
routines work without normaliz.
If you want to use the geometric decomposition
module, you need the normaliz [BIR] command line executable.
The geometric decomposition
module is
designed for normaliz version 3  currently versions
3.0.0
, 3.1.0
, 3.1.1
, 3.3.0
, 3.4.0
, 3.5.4
, 3.6.0
, 3.6.2
, and 3.7.3
are known to work. We recommend to set your $PATH
such that the
normaliz executable is found. Alternatively, you can pass the path to the normaliz
executable directly to the functions that need it.
Drawing Feynman Diagrams with neato¶
In order to use plot_diagram()
, the command line tool
neato must be available. The function loop_package()
tries
to call plot_diagram()
if given a
LoopIntegralFromGraph
and issues a warning on failure. That
warning can be safely ignored if you are not interested in the drawing.
neato is part of the graphviz package. It is available in many package repositories and at http://www.graphviz.org.
Additional Dependencies for Generated c++ Packages¶
The intended main usage of pySecDec is to make it write c++ packages using the functions
pySecDec.code_writer.make_package()
and pySecDec.loop_integral.loop_package()
.
In order to build these c++ packages, the following additional nonpythonbased libraries
and programs are required:
 CUBA (http://www.feynarts.de/cuba/)
 QMC (https://github.com/mppmu/qmc)
 FORM (http://www.nikhef.nl/~form/)
 SecDecUtil (part of pySecDec, see SedDecUtil), depends on:
The functions pySecDec.code_writer.make_package()
and pySecDec.loop_integral.loop_package()
can use the external program nauty [MP+14] to find all sector symmetries and therefore reduce the number of
sectors:
 NAUTY (http://pallini.di.uniroma1.it/)
These packages are redistributed with the pySecDec tarball; i.e. you don’t have to install any of them yourself.
Getting Started¶
After installation, you should have a folder examples in your main pySecDec directory. Here we describe a few of the examples available in the examples directory. A full list of examples is given in List of Examples.
A Simple Example¶
We first show how to compute a simple dimensionally regulated integral:
To run the example change to the easy directory and run the commands:
$ python generate_easy.py
$ make C easy
$ python integrate_easy.py
Additional build options are discussed in the next section. This will evaluate and print the result of the integral:
Numerical Result: + (1.00015897181235158e+00 +/ 4.03392522752491021e03)*eps^1 + (3.06903035514056399e01 +/ 2.82319349818329918e03) + O(eps)
Analytic Result: + (1.000000)*eps^1 + (0.306853) + O(eps)
The file generate_easy.py
defines the integral and calls pySecDec to perform the sector decomposition.
When run it produces the directory easy which contains the code required to numerically evaluate the integral.
The make command builds this code and produces a library.
The file integrate_easy.py
loads the integral library and evaluates the integral.
The user is encouraged to copy and adapt these files to evaluate their own integrals.
Note
If the user is interested in evaluating a loop integral there are many convenience functions that make this much easier. Please see Evaluating a Loop Integral for more details.
In generate_easy.py
we first import make_package
, a function which can decompose, subtract and expand regulated integrals and write a C++ package to evaluate them.
To define our integral we give it a name which will be used as the name of the output directory and C++ namespace.
The integration_variables are declared along with a list of the name of the regulators.
We must specify a list of the requested_orders to which pySecDec should expand our integral in each regulator.
Here we specify requested_orders = [0]
which instructs make_package
to expand the integral up to and including .
Next, we declare the polynomials_to_decompose, here sympy syntax should be used.
from pySecDec import make_package
make_package(
name = 'easy',
integration_variables = ['x','y'],
regulators = ['eps'],
requested_orders = [0],
polynomials_to_decompose = ['(x+y)^(2+eps)'],
)
Once the C++ library has been written and built we run integrate_easy.py
.
Here the library is loaded using IntegralLibrary
.
Calling the instance of IntegralLibrary
with easy_integral()
numerically evaluates the integral and returns the result.
from pySecDec.integral_interface import IntegralLibrary
from math import log
# load c++ library
easy = IntegralLibrary('easy/easy_pylink.so')
# integrate
_, _, result = easy()
# print result
print('Numerical Result:' + result)
print('Analytic Result:' + ' + (%f)*eps^1 + (%f) + O(eps)' % (1.0,1.0log(2.0)))
Evaluating a Loop Integral¶
A simple example of the evaluation of a loop integral with pySecDec is box1L.
This example computes a oneloop box with one offshell leg (with offshellness s1
) and one internal massive line (with mass squared msq
), it is shown in Fig. 2.1.
To run the example change to the box1L directory and run the commands:
$ python generate_box1L.py
$ make C box1L
$ python integrate_box1L.py
This will print the result of the integral evaluated with Mandelstam invariants s=4.0
, t=0.75
and s1=1.25
, msq=1.0
:
eps^2: 0.142868356275422825  1.63596224151119965e6*I +/ ( 0.00118022544307414272 + 0.000210769456586696187*I )
eps^1: 0.639405625715768089 + 1.34277036689902802e6*I +/ ( 0.00650722394065588166 + 0.000971496627153705891*I )
eps^0 : 0.425514350373418893 + 1.86892487760861536*I +/ ( 0.00706834403694714484 + 0.0186497890361357298*I )
The file generate_box1L.py
defines the loop integral and calls pySecDec to perform the sector decomposition. When run it produces the directory box1L which contains the code required to numerically evaluate the integral. The make command builds this code and produces a library. The file integrate_box1L.py
loads the integral library and evaluates the integral for a specified numerical point.
The content of the python files is described in detail in the following sections. The user is encouraged to copy and adapt these files to evaluate their own loop integrals.
Defining a Loop Integral¶
To explain the input format, let us look at generate_box1L.py
from the oneloop box example. The first two lines read
import pySecDec as psd
from pySecDec.loop_integral import loop_package
They say that the module pySecDec should be imported with the alias psd, and that the
function loop_package
from the module loop_integral
is needed.
The following part contains the definition of the loop integral li
:
li = psd.loop_integral.LoopIntegralFromGraph(
# give adjacency list and indicate whether the propagator connecting the numbered vertices is massive or massless in the first entry of each list item.
internal_lines = [['m',[1,2]],[0,[2,3]],[0,[3,4]],[0,[4,1]]],
# contains the names of the external momenta and the label of the vertex they are attached to
external_lines = [['p1',1],['p2',2],['p3',3],['p4',4]],
# define the kinematics and the names for the kinematic invariants
replacement_rules = [
('p1*p1', 's1'),
('p2*p2', 0),
('p3*p3', 0),
('p4*p4', 0),
('p3*p2', 't/2'),
('p1*p2', 's/2s1/2'),
('p1*p4', 't/2s1/2'),
('p2*p4', 's1/2t/2s/2'),
('p3*p4', 's/2'),
('m**2', 'msq')
]
)
Here the class LoopIntegralFromGraph
is used to Feynman parametrize the loop integral given the adjacency list. Alternatively, the class LoopIntegralFromPropagators
can be used to construct the Feynman integral given the momentum representation.
The symbols for the kinematic invariants and the masses also need to be given as an ordered list. The ordering is important as the numerical values assigned to these list elements at the numerical evaluation stage should have the same order.
Mandelstam_symbols = ['s','t','s1']
mass_symbols = ['msq']
Next, the function loop_package
is called. It will create a folder called box1L.
It performs the algebraic sector decomposition steps and writes a package containing the C++ code for the numerical evaluation.
The argument requested_order specifies the order in the regulator to which the integral should be expanded.
For a complete list of possible options see loop_package
.
loop_package(
name = 'box1L',
loop_integral = li,
real_parameters = Mandelstam_symbols + mass_symbols,
# the highest order of the final epsilon expansion > change this value to whatever you think is appropriate
requested_order = 0,
# the optimization level to use in FORM (can be 0, 1, 2, 3, 4)
form_optimization_level = 2,
# the WorkSpace parameter for FORM
form_work_space = '100M',
# the method to be used for the sector decomposition
# valid values are ``iterative`` or ``geometric`` or ``geometric_ku``
decomposition_method = 'iterative',
# if you choose ``geometric[_ku]`` and 'normaliz' is not in your
# $PATH, you can set the path to the 'normaliz' commandline
# executable here
#normaliz_executable='/path/to/normaliz',
)
Building the C++ Library¶
After running the python script generate_box1L.py the folder box1L is created and should contain the following files and subdirectories
Makefile Makefile.conf README box1L.hpp codegen integrate_box1L.cpp cuda_integrate_box1L.cpp pylink src
in the folder box1L, typing
$ make
will create the static library libbox1L.a
and box1L_pylink.so
which can be linked to external programs.
The make
command can also be run in parallel by using the j
option. The number of threads each instance of tform
uses can be
set via the environment variable FORMTHREADS.
New in version 1.4: The environment variable FORMOPT sets FORM’s code optimization level. If not set, the value that was passed to make_package
or loop_package
is used.
To build the dynamic library libbox1L.so
set dynamic
as build target:
$ make dynamic
The code generation with FORM without subsequent compilation can be run by setting source
as build target.
To build the library with nvcc for GPU support, type
$ CXX=nvcc SECDEC_WITH_CUDA=sm_XX make
where sm_XX
must be replaced by the target GPU architechtures, see the arch option of NVCC.
To evaluate the integral numerically a program can call one of these libraries. How to do this interactively or via a python script is explained in the section Python Interface. Alternatively, a C++ program can be produced as explained in the section C++ Interface.
Python Interface (basic)¶
To evaluate the integral for a given numerical point we can use integrate_box1L.py
.
First it imports the necessary python packages and loads the C++ library.
from __future__ import print_function
from pySecDec.integral_interface import IntegralLibrary
import sympy as sp
# load c++ library
box1L = IntegralLibrary('box1L/box1L_pylink.so')
Next, an integrator is configured for the numerical integration. The full list of available integrators and their options is given in integral_interface
.
# choose integrator
box.use_Vegas(flags=2) # ``flags=2``: verbose > see Cuba manual
If you want to use GPUs, change to the CudaQmc
integrator. For example, to run on all available GPUs and CPU cores
using the Korobov transform with weight 3, change the above lines to
# choose integrator
box.use_Qmc(transform='Korobov3')
Calling the box
library numerically evaluates the integral.
Note that the order of the real parameters must match that specified in generate_box1L.py
.
A list of possible settings for the library, in particular details of how to set the contour deformation parameters, is given in IntegralLibrary
.
To change the accuracy settings of the integration, the most important parameters are epsrel
, epsabs
and maxeval
, which
can be added to the integrator argument list:
# choose integrator
box.use_Vegas(flags=2,epsrel=0.01, epsabs=1e07, maxeval=1000000)
# integrate
str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor = box1L(real_parameters=[4.0, 0.75, 1.25, 1.0])
In case of a sign check error (sign_check_error), the arguments number_of_presamples
, deformation_parameters_maximum
, and deformation_parameters_minimum
as described in
IntegralLibrary
can be used to modify the contour.
At this point the string str_integral_with_prefactor
contains the full result of the integral and can be manipulated as required.
In the integrate_box1L.py
an example is shown how to parse the expression with sympy and access individual orders of the regulator.
Note
Instead of parsing the result, it can simply be printed with the line print(str_integral_with_prefactor)
.
# convert complex numbers from c++ to sympy notation
str_integral_with_prefactor = str_integral_with_prefactor.replace(',','+I*')
str_prefactor = str_prefactor.replace(',','+I*')
str_integral_without_prefactor = str_integral_without_prefactor.replace(',','+I*')
# convert result to sympy expressions
integral_with_prefactor = sp.sympify(str_integral_with_prefactor.replace('+/','*value+error*'))
integral_with_prefactor_err = sp.sympify(str_integral_with_prefactor.replace('+/','*value+error*'))
prefactor = sp.sympify(str_prefactor)
integral_without_prefactor = sp.sympify(str_integral_without_prefactor.replace('+/','*value+error*'))
integral_without_prefactor_err = sp.sympify(str_integral_without_prefactor.replace('+/','*value+error*'))
# examples how to access individual orders
print('Numerical Result')
print('eps^2:', integral_with_prefactor.coeff('eps',2).coeff('value'), '+/ (', integral_with_prefactor_err.coeff('eps',2).coeff('error'), ')')
print('eps^1:', integral_with_prefactor.coeff('eps',1).coeff('value'), '+/ (', integral_with_prefactor_err.coeff('eps',1).coeff('error'), ')')
print('eps^0 :', integral_with_prefactor.coeff('eps',0).coeff('value'), '+/ (', integral_with_prefactor_err.coeff('eps',0).coeff('error'), ')')
An example of how to loop over several kinematic points is shown in the example multiple_kinematic_points.py.
C++ Interface (advanced)¶
Usually it is easier to obtain a numerical result using the Python Interface.
However, the library can also be used directly from C++.
Inside the generated box1L folder the file integrate_box1L.cpp
demonstrates this.
The function print_integral_info
shows how to access the important variables of the integral library.
In the main
function a kinematic point must be specified by setting the real_parameters
variable, for example:
int main()
{
// User Specified Phasespace point
const std::vector<box1L::real_t> real_parameters = {4.0, 0.75, 1.25, 1.0}; // EDIT: kinematic point specified here
const std::vector<box1L::complex_t> complex_parameters = { };
The name::make_integrands()
function returns an secdecutil::IntegrandContainer
for each sector and regulator order:
// Generate the integrands (optimization of the contour if applicable)
const std::vector<box1L::nested_series_t<box1L::integrand_t>> sector_integrands = box1L::make_integrands(real_parameters, complex_parameters);
The contour deformation has to be adjusted in case of a sign check error (sign_check_error). This can be done via additional arguments to name::make_integrands()
.
The sectors can be added before integration:
// Add integrands of sectors (together flag)
const box1L::nested_series_t<box1L::integrand_t> all_sectors = std::accumulate(++sector_integrands.begin(), sector_integrands.end(), *sector_integrands.begin() );
An secdecutil::Integrator
is constructed and its parameters are set:
// Integrate
secdecutil::cuba::Vegas<box1L::integrand_return_t> integrator;
integrator.flags = 2; // verbose output > see cuba manual
To numerically integrate the functions the secdecutil::Integrator::integrate()
function is applied to each secdecutil::IntegrandContainer
using secdecutil::deep_apply()
:
const box1L::nested_series_t<secdecutil::UncorrelatedDeviation<box1L::integrand_return_t>> result_all = secdecutil::deep_apply( all_sectors, integrator.integrate );
The remaining lines print the result:
std::cout << "" << std::endl << std::endl;
std::cout << " integral info  " << std::endl;
print_integral_info();
std::cout << std::endl;
std::cout << " integral without prefactor  " << std::endl;
std::cout << result_all << std::endl << std::endl;
std::cout << " prefactor  " << std::endl;
const box1L::nested_series_t<box1L::integrand_return_t> prefactor = box1L::prefactor(real_parameters, complex_parameters);
std::cout << prefactor << std::endl << std::endl;
std::cout << " full result (prefactor*integral)  " << std::endl;
std::cout << prefactor*result_all << std::endl;
return 0;
}
After editing the real_parameters
as described above the C++ program can be built and executed with the commands
$ make integrate_box1L
$ ./integrate_box1L
New in version 1.4.
The similar template file cuda_integrate_box1L.cpp
provides an example to run on GPUs. The main differences are in the lines that generate, add, and integrate the integrands.
Rather than name::make_integrands()
, name::make_cuda_integrands()
is called:
// Generate the integrands (optimization of the contour if applicable)
const std::vector<box1L::nested_series_t<box1L::cuda_integrand_t>> sector_integrands = box1L::make_cuda_integrands(real_parameters, complex_parameters);
If the integrands are added together before integration, the sum command is as follows:
// Add integrands of sectors (together flag)
const box1L::nested_series_t<box1L::cuda_together_integrand_t> all_sectors =
std::accumulate(++sector_integrands.begin(), sector_integrands.end(), box1L::cuda_together_integrand_t()+*sector_integrands.begin());
Note the conversion from name::cuda_integrand_t
to name::cuda_together_integrand_t
. The CUDAcapable version of the Qmc
integrator takes additional the template arguments box1L::maximal_number_of_integration_variables
, integrators::transforms::Korobov<3>::type
,
and name::cuda_integrand_t
:
// Integrate
secdecutil::integrators::Qmc<
box1L::integrand_return_t,
box1L::maximal_number_of_integration_variables,
integrators::transforms::Korobov<3>::type, // EDIT: integral transform specified to "Korobov<3>"
box1L::cuda_together_integrand_t
> integrator;
integrator.verbosity = 1;
const box1L::nested_series_t<secdecutil::UncorrelatedDeviation<box1L::integrand_return_t>> result_all = secdecutil::deep_apply( all_sectors, integrator.integrate );
If the integrands are integrated separately, name::cuda_together_integrand_t
should be changed to name::cuda_integrand_t
. If your integral
is higher than seven dimensional, changing the integral transform to integrators::transforms::Baker::type
may improve the accuracy of the result. For further
options of the integrator we refer to Section 4.5.2.
List of Examples¶
Here we list the available examples. For more details regarding each example see [PSD17] and [PSD18].
easy:  a simple parametric integral, described in Section 2.1 
easy_cuda:  the same integral as in easy but computed on GPUs with CUDA 
box1L:  a simple 1loop, 4point, 4propagator integral, described in Section 2.2 
triangle2L:  a 2loop, 3point, 6propagator diagram, also known as P126 
box2L_numerator:  a massless planar onshell 2loop, 4point, 7propagator box with a numerator, either defined as an inverse propagator
box2L_invprop.py or in terms of contracted Lorentz vectors box2L_contracted_tensor.py 
pentabox_fin:  a 2loop, 5point, 8propagator diagram, evaluated in dimensions where it is finite 
triangle3L:  a 2loop, 3point, 7propagator integral, demonstrates that the symmetry finder can significantly reduce the number of sectors 
formfactor4L:  a singlescale 4loop 3point integral in dimensions 
bubble6L:  a singlescale 6loop 2point integral, evaluated at a Euclidean phasespace point 
elliptic2L_euclidean:  an integral known to contain elliptic functions, evaluated at a Euclidean phasespace point 
elliptic2L_physical:  an integral known to contain elliptic functions, evaluated at a physical phasespace point 
banana_3mass:  a 3loop 2point integral with three different internal masses known to contain hyperelliptic functions, evaluated at a physical phasespace point 
hyperelliptic:  a 2loop 4point nonplanar integral known to contain hyperelliptic functions, evaluated at a physical phasespace point 
triangle2L_split:  a 2loop, 3point, 6propagator integral without a Euclidean region due to special kinematics 
Nbox2L_split:  three 2loop, 4point, 5propagator integrals that need split=True due to special kinematics 
hypergeo5F4:  a general dimensionally regulated parameter integral 
4photon1L_amplitude:  calculation of the 4photon amplitude, showing how to use pySecDec as an integral library in a larger context 
two_regulators:  an integral involving poles in two different regulators. 
userdefined_cpp:  a collection of examples demonstrating how to combine polynomials to be decomposed with other userdefined functions 
Overview¶
pySecDec consists of several modules that provide functions and classes for specific purposes. In this overview, we present only the most important aspects of selected modules. These are exactly the modules necessary to set up the algebraic computation of a Feynman loop integral requisite for the numerical evaluation. For detailed instruction of a specific function or class, please be referred to the reference guide.
The Algebra Module¶
The algebra module implements a very basic computer algebra system. pySecDec uses both sympy and numpy. Although sympy in principle provides everything we need, it is way too slow for typical applications. That is because sympy is completely written in python without making use of any precompiled functions. pySecDec’s algebra module uses the in general faster numpy function wherever possible.
Polynomials¶
Since sector decomposition is an algorithm that acts on polynomials, we start with
the key class Polynomial
.
As the name suggests, the Polynomial
class
is a container for multivariate polynomials, i.e. functions of the form:
A multivariate polynomial is completely determined by its coefficients and
the exponents . The Polynomial
class stores these in two arrays:
>>> from pySecDec.algebra import Polynomial
>>> poly = Polynomial([[1,0], [0,2]], ['A', 'B'])
>>> poly
+ (A)*x0 + (B)*x1**2
>>> poly.expolist
array([[1, 0],
[0, 2]])
>>> poly.coeffs
array([A, B], dtype=object)
It is also possible to instantiate the Polynomial
by its algebraic representation:
>>> poly2 = Polynomial.from_expression('A*x0 + B*x1**2', ['x0','x1'])
>>> poly2
+ (A)*x0 + (B)*x1**2
>>> poly2.expolist
array([[1, 0],
[0, 2]])
>>> poly2.coeffs
array([A, B], dtype=object)
Note that the second argument of
Polynomial.from_expression()
defines the variables .
Within the Polynomial
class, basic operations are implemented:
>>> poly + 1
+ (1) + (B)*x1**2 + (A)*x0
>>> 2 * poly
+ (2*A)*x0 + (2*B)*x1**2
>>> poly + poly
+ (2*B)*x1**2 + (2*A)*x0
>>> poly * poly
+ (B**2)*x1**4 + (2*A*B)*x0*x1**2 + (A**2)*x0**2
>>> poly ** 2
+ (B**2)*x1**4 + (2*A*B)*x0*x1**2 + (A**2)*x0**2
General Expressions¶
In order to perform the pySecDec.subtraction
and pySecDec.expansion
,
we have to introduce more complex algebraic constructs.
General expressions can be entered in a straightforward way:
>>> from pySecDec.algebra import Expression
>>> log_of_x = Expression('log(x)', ['x'])
>>> log_of_x
log( + (1)*x)
All expressions in the context of this algebra module are based
on extending or combining the Polynomials
introduced above.
In the example above, log_of_x
is a
LogOfPolynomial
, which
is a derived class from Polynomial
:
>>> type(log_of_x)
<class 'pySecDec.algebra.LogOfPolynomial'>
>>> isinstance(log_of_x, Polynomial)
True
>>> log_of_x.expolist
array([[1]])
>>> log_of_x.coeffs
array([1], dtype=object)
We have seen an extension to the
Polynomial
class, now let us consider
a combination:
>>> more_complex_expression = log_of_x * log_of_x
>>> more_complex_expression
(log( + (1)*x)) * (log( + (1)*x))
We just introduced the Product
of two LogOfPolynomials
:
>>> type(more_complex_expression)
<class 'pySecDec.algebra.Product'>
As suggested before, the Product
combines two Polynomials
. They
are accessible through the factors
:
>>> more_complex_expression.factors[0]
log( + (1)*x)
>>> more_complex_expression.factors[1]
log( + (1)*x)
>>> type(more_complex_expression.factors[0])
<class 'pySecDec.algebra.LogOfPolynomial'>
>>> type(more_complex_expression.factors[1])
<class 'pySecDec.algebra.LogOfPolynomial'>
Important
When working with this algebra module, it is important to understand that
everything is based on the class
Polynomial
.
To emphasize the importance of the above statement, consider the following code:
>>> expression1 = Expression('x*y', ['x', 'y'])
>>> expression2 = Expression('x*y', ['x'])
>>> type(expression1)
<class 'pySecDec.algebra.Polynomial'>
>>> type(expression2)
<class 'pySecDec.algebra.Polynomial'>
>>> expression1
+ (1)*x*y
>>> expression2
+ (y)*x
Although expression1
and expression2
are mathematically identical,
they are treated differently by the algebra module. In expression1
, both,
x
and y
, are considered as variables of the
Polynomial
. In contrast, y
is treated
as coefficient in expression2
:
>>> expression1.expolist
array([[1, 1]])
>>> expression1.coeffs
array([1], dtype=object)
>>> expression2.expolist
array([[1]])
>>> expression2.coeffs
array([y], dtype=object)
The second argument of the function Expression
controls how the variables are distributed among the coefficients and the variables
in the underlying Polynomials
.
Keep that in mind in order to avoid confusion. One can always check which symbols are
considered as variables by asking for the symbols
:
>>> expression1.symbols
[x, y]
>>> expression2.symbols
[x]
Feynman Parametrization of Loop Integrals¶
The primary purpose of pySecDec is the numerical calculation of loop integrals as they arise in fixed order calculations in quantum field theories. In the first step of our approach, the loop integral is converted from the momentum representation to the Feynman parameter representation, see for example [Hei08] (Chapter 3).
The module pySecDec.loop_integral
implements exactly that conversion.
The most basic use is to calculate the first and the second
Symanzik polynomial U
and F
, respectively, from the propagators of a loop integral.
One Loop Bubble¶
To calculate U
and F
of the one loop bubble, type the following
commands:
>>> from pySecDec.loop_integral import LoopIntegralFromPropagators
>>> propagators = ['k**2', '(k  p)**2']
>>> loop_momenta = ['k']
>>> one_loop_bubble = LoopIntegralFromPropagators(propagators, loop_momenta)
>>> one_loop_bubble.U
+ (1)*x0 + (1)*x1
>>> one_loop_bubble.F
+ (p**2)*x0*x1
The example above among other useful features is also stated in the full documenation of
LoopIntegralFromPropagators()
in the reference guide.
Two Loop Planar Box with Numerator¶
Consider the propagators of the two loop planar box:
>>> propagators = ['k1**2','(k1+p2)**2',
... '(k1p1)**2','(k1k2)**2',
... '(k2+p2)**2','(k2p1)**2',
... '(k2+p2+p3)**2']
>>> loop_momenta = ['k1','k2']
We could now instantiate the LoopIntegral
just like before. However, let us consider an additional numerator:
>>> numerator = 'k1(mu)*k1(mu) + 2*k1(mu)*p3(mu) + p3(mu)*p3(mu)' # (k1 + p3) ** 2
In order to unambiguously define the loop integral, we must state which
symbols denote the Lorentz_indices
(just mu
in this case here) and the external_momenta
:
>>> external_momenta = ['p1','p2','p3','p4']
>>> Lorentz_indices=['mu']
With that, we can Feynman parametrize the two loop box with a numerator:
>>> box = LoopIntegralFromPropagators(propagators, loop_momenta, external_momenta,
... numerator=numerator, Lorentz_indices=Lorentz_indices)
>>> box.U
+ (1)*x3*x6 + (1)*x3*x5 + (1)*x3*x4 + (1)*x2*x6 + (1)*x2*x5 + (1)*x2*x4 + (1)*x2*x3 + (1)*x1*x6 + (1)*x1*x5 + (1)*x1*x4 + (1)*x1*x3 + (1)*x0*x6 + (1)*x0*x5 + (1)*x0*x4 + (1)*x0*x3
>>> box.F
+ (p1**2  2*p1*p2  2*p1*p3  p2**2  2*p2*p3  p3**2)*x3*x5*x6 + (p3**2)*x3*x4*x6 + (p1**2  2*p1*p2  p2**2)*x3*x4*x5 + (p1**2  2*p1*p2  2*p1*p3  p2**2  2*p2*p3  p3**2)*x2*x5*x6 + (p3**2)*x2*x4*x6 + (p1**2  2*p1*p2  p2**2)*x2*x4*x5 + (p1**2  2*p1*p2  2*p1*p3  p2**2  2*p2*p3  p3**2)*x2*x3*x6 + (p1**2  2*p1*p2  p2**2)*x2*x3*x4 + (p1**2  2*p1*p2  2*p1*p3  p2**2  2*p2*p3  p3**2)*x1*x5*x6 + (p3**2)*x1*x4*x6 + (p1**2  2*p1*p2  p2**2)*x1*x4*x5 + (p3**2)*x1*x3*x6 + (p1**2  2*p1*p2  p2**2)*x1*x3*x5 + (p1**2  2*p1*p2  p2**2)*x1*x2*x6 + (p1**2  2*p1*p2  p2**2)*x1*x2*x5 + (p1**2  2*p1*p2  p2**2)*x1*x2*x4 + (p1**2  2*p1*p2  p2**2)*x1*x2*x3 + (p1**2  2*p1*p2  2*p1*p3  p2**2  2*p2*p3  p3**2)*x0*x5*x6 + (p3**2)*x0*x4*x6 + (p1**2  2*p1*p2  p2**2)*x0*x4*x5 + (p2**2  2*p2*p3  p3**2)*x0*x3*x6 + (p1**2)*x0*x3*x5 + (p2**2)*x0*x3*x4 + (p1**2)*x0*x2*x6 + (p1**2)*x0*x2*x5 + (p1**2)*x0*x2*x4 + (p1**2)*x0*x2*x3 + (p2**2)*x0*x1*x6 + (p2**2)*x0*x1*x5 + (p2**2)*x0*x1*x4 + (p2**2)*x0*x1*x3
>>> box.numerator
+ (2*eps*p3(mu)**2 + 2*p3(mu)**2)*U**2 + (eps  2)*x6*F + (eps  2)*x5*F + (eps  2)*x4*F + (eps  2)*x3*F + (4*eps*p2(mu)*p3(mu)  4*eps*p3(mu)**2  4*p2(mu)*p3(mu)  4*p3(mu)**2)*x3*x6*U + (4*eps*p1(mu)*p3(mu) + 4*p1(mu)*p3(mu))*x3*x5*U + (4*eps*p2(mu)*p3(mu)  4*p2(mu)*p3(mu))*x3*x4*U + (2*eps*p2(mu)**2 + 4*eps*p2(mu)*p3(mu) + 2*eps*p3(mu)**2 + 2*p2(mu)**2 + 4*p2(mu)*p3(mu) + 2*p3(mu)**2)*x3**2*x6**2 + (4*eps*p1(mu)*p2(mu)  4*eps*p1(mu)*p3(mu)  4*p1(mu)*p2(mu)  4*p1(mu)*p3(mu))*x3**2*x5*x6 + (2*eps*p1(mu)**2 + 2*p1(mu)**2)*x3**2*x5**2 + (4*eps*p2(mu)**2 + 4*eps*p2(mu)*p3(mu) + 4*p2(mu)**2 + 4*p2(mu)*p3(mu))*x3**2*x4*x6 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x3**2*x4*x5 + (2*eps*p2(mu)**2 + 2*p2(mu)**2)*x3**2*x4**2 + (4*eps*p1(mu)*p3(mu) + 4*p1(mu)*p3(mu))*x2*x6*U + (4*eps*p1(mu)*p3(mu) + 4*p1(mu)*p3(mu))*x2*x5*U + (4*eps*p1(mu)*p3(mu) + 4*p1(mu)*p3(mu))*x2*x4*U + (4*eps*p1(mu)*p3(mu) + 4*p1(mu)*p3(mu))*x2*x3*U + (4*eps*p1(mu)*p2(mu)  4*eps*p1(mu)*p3(mu)  4*p1(mu)*p2(mu)  4*p1(mu)*p3(mu))*x2*x3*x6**2 + (4*eps*p1(mu)**2  4*eps*p1(mu)*p2(mu)  4*eps*p1(mu)*p3(mu) + 4*p1(mu)**2  4*p1(mu)*p2(mu)  4*p1(mu)*p3(mu))*x2*x3*x5*x6 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2*x3*x5**2 + (8*eps*p1(mu)*p2(mu)  4*eps*p1(mu)*p3(mu)  8*p1(mu)*p2(mu)  4*p1(mu)*p3(mu))*x2*x3*x4*x6 + (4*eps*p1(mu)**2  4*eps*p1(mu)*p2(mu) + 4*p1(mu)**2  4*p1(mu)*p2(mu))*x2*x3*x4*x5 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x2*x3*x4**2 + (4*eps*p1(mu)*p2(mu)  4*eps*p1(mu)*p3(mu)  4*p1(mu)*p2(mu)  4*p1(mu)*p3(mu))*x2*x3**2*x6 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2*x3**2*x5 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x2*x3**2*x4 + (2*eps*p1(mu)**2 + 2*p1(mu)**2)*x2**2*x6**2 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2**2*x5*x6 + (2*eps*p1(mu)**2 + 2*p1(mu)**2)*x2**2*x5**2 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2**2*x4*x6 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2**2*x4*x5 + (2*eps*p1(mu)**2 + 2*p1(mu)**2)*x2**2*x4**2 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2**2*x3*x6 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2**2*x3*x5 + (4*eps*p1(mu)**2 + 4*p1(mu)**2)*x2**2*x3*x4 + (2*eps*p1(mu)**2 + 2*p1(mu)**2)*x2**2*x3**2 + (4*eps*p2(mu)*p3(mu)  4*p2(mu)*p3(mu))*x1*x6*U + (4*eps*p2(mu)*p3(mu)  4*p2(mu)*p3(mu))*x1*x5*U + (4*eps*p2(mu)*p3(mu)  4*p2(mu)*p3(mu))*x1*x4*U + (4*eps*p2(mu)*p3(mu)  4*p2(mu)*p3(mu))*x1*x3*U + (4*eps*p2(mu)**2 + 4*eps*p2(mu)*p3(mu) + 4*p2(mu)**2 + 4*p2(mu)*p3(mu))*x1*x3*x6**2 + (4*eps*p1(mu)*p2(mu) + 4*eps*p2(mu)**2 + 4*eps*p2(mu)*p3(mu)  4*p1(mu)*p2(mu) + 4*p2(mu)**2 + 4*p2(mu)*p3(mu))*x1*x3*x5*x6 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x1*x3*x5**2 + (8*eps*p2(mu)**2 + 4*eps*p2(mu)*p3(mu) + 8*p2(mu)**2 + 4*p2(mu)*p3(mu))*x1*x3*x4*x6 + (4*eps*p1(mu)*p2(mu) + 4*eps*p2(mu)**2  4*p1(mu)*p2(mu) + 4*p2(mu)**2)*x1*x3*x4*x5 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1*x3*x4**2 + (4*eps*p2(mu)**2 + 4*eps*p2(mu)*p3(mu) + 4*p2(mu)**2 + 4*p2(mu)*p3(mu))*x1*x3**2*x6 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x1*x3**2*x5 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1*x3**2*x4 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x1*x2*x6**2 + (8*eps*p1(mu)*p2(mu)  8*p1(mu)*p2(mu))*x1*x2*x5*x6 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x1*x2*x5**2 + (8*eps*p1(mu)*p2(mu)  8*p1(mu)*p2(mu))*x1*x2*x4*x6 + (8*eps*p1(mu)*p2(mu)  8*p1(mu)*p2(mu))*x1*x2*x4*x5 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x1*x2*x4**2 + (8*eps*p1(mu)*p2(mu)  8*p1(mu)*p2(mu))*x1*x2*x3*x6 + (8*eps*p1(mu)*p2(mu)  8*p1(mu)*p2(mu))*x1*x2*x3*x5 + (8*eps*p1(mu)*p2(mu)  8*p1(mu)*p2(mu))*x1*x2*x3*x4 + (4*eps*p1(mu)*p2(mu)  4*p1(mu)*p2(mu))*x1*x2*x3**2 + (2*eps*p2(mu)**2 + 2*p2(mu)**2)*x1**2*x6**2 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1**2*x5*x6 + (2*eps*p2(mu)**2 + 2*p2(mu)**2)*x1**2*x5**2 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1**2*x4*x6 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1**2*x4*x5 + (2*eps*p2(mu)**2 + 2*p2(mu)**2)*x1**2*x4**2 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1**2*x3*x6 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1**2*x3*x5 + (4*eps*p2(mu)**2 + 4*p2(mu)**2)*x1**2*x3*x4 + (2*eps*p2(mu)**2 + 2*p2(mu)**2)*x1**2*x3**2
We can also generate the output in terms of Mandelstam invariants:
>>> replacement_rules = [
... ('p1*p1', 0),
... ('p2*p2', 0),
... ('p3*p3', 0),
... ('p4*p4', 0),
... ('p1*p2', 's/2'),
... ('p2*p3', 't/2'),
... ('p1*p3', 's/2t/2')
... ]
>>> box = LoopIntegralFromPropagators(propagators, loop_momenta, external_momenta,
... numerator=numerator, Lorentz_indices=Lorentz_indices,
... replacement_rules=replacement_rules)
>>> box.U
+ (1)*x3*x6 + (1)*x3*x5 + (1)*x3*x4 + (1)*x2*x6 + (1)*x2*x5 + (1)*x2*x4 + (1)*x2*x3 + (1)*x1*x6 + (1)*x1*x5 + (1)*x1*x4 + (1)*x1*x3 + (1)*x0*x6 + (1)*x0*x5 + (1)*x0*x4 + (1)*x0*x3
>>> box.F
+ (s)*x3*x4*x5 + (s)*x2*x4*x5 + (s)*x2*x3*x4 + (s)*x1*x4*x5 + (s)*x1*x3*x5 + (s)*x1*x2*x6 + (s)*x1*x2*x5 + (s)*x1*x2*x4 + (s)*x1*x2*x3 + (s)*x0*x4*x5 + (t)*x0*x3*x6
>>> box.numerator
+ (eps  2)*x6*F + (eps  2)*x5*F + (eps  2)*x4*F + (eps  2)*x3*F + (2*eps*t  2*t)*x3*x6*U + (4*eps*(s/2  t/2)  2*s  2*t)*x3*x5*U + (2*eps*t  2*t)*x3*x4*U + (2*eps*t + 2*t)*x3**2*x6**2 + (2*eps*s  4*eps*(s/2  t/2) + 2*t)*x3**2*x5*x6 + (2*eps*t + 2*t)*x3**2*x4*x6 + (2*eps*s  2*s)*x3**2*x4*x5 + (4*eps*(s/2  t/2)  2*s  2*t)*x2*x6*U + (4*eps*(s/2  t/2)  2*s  2*t)*x2*x5*U + (4*eps*(s/2  t/2)  2*s  2*t)*x2*x4*U + (4*eps*(s/2  t/2)  2*s  2*t)*x2*x3*U + (2*eps*s  4*eps*(s/2  t/2) + 2*t)*x2*x3*x6**2 + (2*eps*s  4*eps*(s/2  t/2) + 2*t)*x2*x3*x5*x6 + (4*eps*s  4*eps*(s/2  t/2)  2*s + 2*t)*x2*x3*x4*x6 + (2*eps*s  2*s)*x2*x3*x4*x5 + (2*eps*s  2*s)*x2*x3*x4**2 + (2*eps*s  4*eps*(s/2  t/2) + 2*t)*x2*x3**2*x6 + (2*eps*s  2*s)*x2*x3**2*x4 + (2*eps*t  2*t)*x1*x6*U + (2*eps*t  2*t)*x1*x5*U + (2*eps*t  2*t)*x1*x4*U + (2*eps*t  2*t)*x1*x3*U + (2*eps*t + 2*t)*x1*x3*x6**2 + (2*eps*s + 2*eps*t  2*s + 2*t)*x1*x3*x5*x6 + (2*eps*s  2*s)*x1*x3*x5**2 + (2*eps*t + 2*t)*x1*x3*x4*x6 + (2*eps*s  2*s)*x1*x3*x4*x5 + (2*eps*t + 2*t)*x1*x3**2*x6 + (2*eps*s  2*s)*x1*x3**2*x5 + (2*eps*s  2*s)*x1*x2*x6**2 + (4*eps*s  4*s)*x1*x2*x5*x6 + (2*eps*s  2*s)*x1*x2*x5**2 + (4*eps*s  4*s)*x1*x2*x4*x6 + (4*eps*s  4*s)*x1*x2*x4*x5 + (2*eps*s  2*s)*x1*x2*x4**2 + (4*eps*s  4*s)*x1*x2*x3*x6 + (4*eps*s  4*s)*x1*x2*x3*x5 + (4*eps*s  4*s)*x1*x2*x3*x4 + (2*eps*s  2*s)*x1*x2*x3**2
Sector Decomposition¶
The sector decomposition algorithm aims to factorize the polynomials as products of a monomial and a polynomial with nonzero constant term:
Factorizing polynomials in that way by expoliting integral transformations is the first step in an algorithm for solving dimensionally regulated integrals of the form
The iterative sector decomposition splits the integral and remaps the integration domain until all polynomials in all arising integrals (called sectors) have the desired form . An introduction to the sector decomposition approach can be found in [Hei08].
To demonstrate the pySecDec.decomposition
module, we decompose the polynomials
>>> p1 = Polynomial.from_expression('x + A*y', ['x','y','z'])
>>> p2 = Polynomial.from_expression('x + B*y*z', ['x','y','z'])
Let us first focus on the iterative decomposition of p1
. In the pySecDec
framework, we first have to pack p1
into a Sector
:
>>> from pySecDec.decomposition import Sector
>>> initial_sector = Sector([p1])
>>> print(initial_sector)
Sector:
Jacobian= + (1)
cast=[( + (1)) * ( + (1)*x + (A)*y)]
other=[]
We can now run the iterative decomposition and take a look at the decomposed sectors:
>>> from pySecDec.decomposition.iterative import iterative_decomposition
>>> decomposed_sectors = iterative_decomposition(initial_sector)
>>> for sector in decomposed_sectors:
... print(sector)
... print('\n')
...
Sector:
Jacobian= + (1)*x
cast=[( + (1)*x) * ( + (1) + (A)*y)]
other=[]
Sector:
Jacobian= + (1)*y
cast=[( + (1)*y) * ( + (1)*x + (A))]
other=[]
The decomposition of p2
needs two iterations and yields three sectors:
>>> initial_sector = Sector([p2])
>>> decomposed_sectors = iterative_decomposition(initial_sector)
>>> for sector in decomposed_sectors:
... print(sector)
... print('\n')
...
Sector:
Jacobian= + (1)*x
cast=[( + (1)*x) * ( + (1) + (B)*y*z)]
other=[]
Sector:
Jacobian= + (1)*x*y
cast=[( + (1)*x*y) * ( + (1) + (B)*z)]
other=[]
Sector:
Jacobian= + (1)*y*z
cast=[( + (1)*y*z) * ( + (1)*x + (B))]
other=[]
Note that we declared z
as a variable for sector p1
evne though it does not depend on it.
This declaration is necessary if we want to simultaneously decompose p1
and p2
:
>>> initial_sector = Sector([p1, p2])
>>> decomposed_sectors = iterative_decomposition(initial_sector)
>>> for sector in decomposed_sectors:
... print(sector)
... print('\n')
...
Sector:
Jacobian= + (1)*x
cast=[( + (1)*x) * ( + (1) + (A)*y), ( + (1)*x) * ( + (1) + (B)*y*z)]
other=[]
Sector:
Jacobian= + (1)*x*y
cast=[( + (1)*y) * ( + (1)*x + (A)), ( + (1)*x*y) * ( + (1) + (B)*z)]
other=[]
Sector:
Jacobian= + (1)*y*z
cast=[( + (1)*y) * ( + (1)*x*z + (A)), ( + (1)*y*z) * ( + (1)*x + (B))]
other=[]
We just fully decomposed p1
and p2
. In some cases, one may want to bring
one polynomial, say p1
, into standard form, but not neccessarily the other.
For that purpose, the Sector
can take
a second argument. In the following code example, we bring p1
into standard
form, apply all transformations to p2
as well, but stop before p2
is fully
decomposed:
>>> initial_sector = Sector([p1], [p2])
>>> decomposed_sectors = iterative_decomposition(initial_sector)
>>> for sector in decomposed_sectors:
... print(sector)
... print('\n')
...
Sector:
Jacobian= + (1)*x
cast=[( + (1)*x) * ( + (1) + (A)*y)]
other=[ + (1)*x + (B)*x*y*z]
Sector:
Jacobian= + (1)*y
cast=[( + (1)*y) * ( + (1)*x + (A))]
other=[ + (1)*x*y + (B)*y*z]
Subtraction¶
In the subtraction, we want to perform those integrations that lead to divergencies. The master formula for one integration variables is
where is denotes the pth derivative of with respect to . The equation above effectively defines the remainder term . All terms on the right hand side of the equation above are constructed to be free of divergencies. For more details and the generalization to multiple variables, we refer the reader to [Hei08]. In the following, we show how to use the implementation in pySecDec.
To initialize the subtraction, we first define a factorized expression of the form :
>>> from pySecDec.algebra import Expression
>>> symbols = ['x','y','eps']
>>> x_monomial = Expression('x**(1  b_x*eps)', symbols)
>>> y_monomial = Expression('y**(2  b_y*eps)', symbols)
>>> cal_I = Expression('cal_I(x, y, eps)', symbols)
We must pack the monomials into a pySecDec.algebra.Product
:
>>> from pySecDec.algebra import Product
>>> monomials = Product(x_monomial, y_monomial)
Although this seems to be to complete input according to the equation
above, we are still missing a structure to store poles in. The function
pySecDec.subtraction.integrate_pole_part()
is designed to return
an iterable of the same type as the input. That is particularly important
since the output of the subtraction of one variable is the input for the
subtraction of the next variable. We will see this iteration later. Initially,
we do not have poles yet, therefore we define a one of the required type:
>>> from pySecDec.algebra import Pow
>>> import numpy as np
>>> polynomial_one = Polynomial(np.zeros([1,len(symbols)], dtype=int), np.array([1]), symbols, copy=False)
>>> pole_part_initializer = Pow(polynomial_one, polynomial_one)
pole_part_initializer
is of type pySecDec.algebra.Pow
and has polynomial_one
in the exponent. We initialize the base with polynomial_one
; i.e. a one packed into
a polynomial. The function pySecDec.subtraction.integrate_pole_part()
populates the
base with factors of when poles arise.
We are now ready to build the subtraction_initializer
 the pySecDec.algebra.Product
to be passed into pySecDec.subtraction.integrate_pole_part()
.
>>> from pySecDec.subtraction import integrate_pole_part
>>> subtraction_initializer = Product(monomials, pole_part_initializer, cal_I)
>>> x_subtracted = integrate_pole_part(subtraction_initializer, 0)
The second argument of pySecDec.subtraction.integrate_pole_part()
specifies
to which variable we want to apply the master formula, here we choose .
First, remember that the x monomial is a dimensionally regulated .
Therefore, the sum collapses to only one term and we have two terms in total.
Each term corresponds to one entry in the list x_subtracted
:
>>> len(x_subtracted)
2
x_subtracted
has the same structure as our input. The first factor of each term
stores the remaining monomials:
>>> x_subtracted[0].factors[0]
(( + (1))**( + (b_x)*eps + (1))) * (( + (1)*y)**( + (b_y)*eps + (2)))
>>> x_subtracted[1].factors[0]
(( + (1)*x)**( + (b_x)*eps + (1))) * (( + (1)*y)**( + (b_y)*eps + (2)))
The second factor stores the poles. There is an epsilon pole in the first term, but still none in the second:
>>> x_subtracted[0].factors[1]
( + (b_x)*eps) ** ( + (1))
>>> x_subtracted[1].factors[1]
( + (1)) ** ( + (1))
The last factor catches everything that is not covered by the first two fields:
>>> x_subtracted[0].factors[2]
(cal_I( + (0), + (1)*y, + (1)*eps))
>>> x_subtracted[1].factors[2]
(cal_I( + (1)*x, + (1)*y, + (1)*eps)) + (( + (1)) * (cal_I( + (0), + (1)*y, + (1)*eps)))
We have now performed the subtraction for . Because in and output have a similar structure, we can easily perform the subtraction for as well:
>>> x_and_y_subtracted = []
>>> for s in x_subtracted:
... x_and_y_subtracted.extend( integrate_pole_part(s,1) )
Alternatively, we can directly instruct pySecDec.subtraction.integrate_pole_part()
to perform both subtractions:
>>> alternative_x_and_y_subtracted = integrate_pole_part(subtraction_initializer,0,1)
In both cases, the result is a list of the terms appearing on the right hand side of the master equation.
Expansion¶
The purpose of the expansion
module is,
as the name suggests, to provide routines to perform a series expansion.
The module basically implements two routines  the Taylor expansion
(pySecDec.expansion.expand_Taylor()
) and an expansion of polyrational
functions supporting singularities in the expansion variable
(pySecDec.expansion.expand_singular()
).
Taylor Expansion¶
The function pySecDec.expansion.expand_Taylor()
implements the ordinary
Taylor expansion. It takes an algebraic expression (in the sense of the
algebra module, the index of the expansion variable
and the order to which the expression shall be expanded:
>>> from pySecDec.algebra import Expression
>>> from pySecDec.expansion import expand_Taylor
>>> expression = Expression('x**eps', ['eps'])
>>> expand_Taylor(expression, 0, 2).simplify()
+ (1) + (log( + (x)))*eps + ((log( + (x))) * (log( + (x))) * ( + (1/2)))*eps**2
It is also possible to expand an expression in multiple variables simultaneously:
>>> expression = Expression('x**(eps + alpha)', ['eps', 'alpha'])
>>> expand_Taylor(expression, [0,1], [2,0]).simplify()
+ (1) + (log( + (x)))*eps + ((log( + (x))) * (log( + (x))) * ( + (1/2)))*eps**2
The command above instructs pySecDec.expansion.expand_Taylor()
to expand
the expression
to the second order in the variable indexed 0
(eps
)
and to the zeroth order in the variable indexed 1
(alpha
).
Laurent Expansion¶
pySecDec.expansion.expand_singular()
Laurent expands polyrational functions.
Its input is more restrictive than for the Taylor expansion.
It expects a Product
where the factors are either
Polynomials
or
ExponentiatedPolynomials
with exponent = 1
:
>>> from pySecDec.expansion import expand_singular
>>> expression = Expression('1/(eps + alpha)', ['eps', 'alpha']).simplify()
>>> expand_singular(expression, 0, 1)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/pcl340a/sjahn/Projects/pySecDec/pySecDec/expansion.py", line 241, in expand_singular
return _expand_and_flatten(product, indices, orders, _expand_singular_step)
File "/home/pcl340a/sjahn/Projects/pySecDec/pySecDec/expansion.py", line 209, in _expand_and_flatten
expansion = recursive_expansion(expression, indices, orders)
File "/home/pcl340a/sjahn/Projects/pySecDec/pySecDec/expansion.py", line 198, in recursive_expansion
expansion = expansion_one_variable(expression, index, order)
File "/home/pcl340a/sjahn/Projects/pySecDec/pySecDec/expansion.py", line 82, in _expand_singular_step
raise TypeError('`product` must be a `Product`')
TypeError: `product` must be a `Product`
>>> expression # ``expression`` is indeed a polyrational function.
( + (1)*alpha + (1)*eps)**(1)
>>> type(expression) # It is just not packed in a ``Product`` as ``expand_singular`` expects.
<class 'pySecDec.algebra.ExponentiatedPolynomial'>
>>> from pySecDec.algebra import Product
>>> expression = Product(expression)
>>> expand_singular(expression, 0, 1)
+ (( + (1)) * (( + (1)*alpha)**(1))) + (( + (1)) * (( + (1)*alpha**2)**(1)))*eps
Like in the Taylor expansion, we can expand simultaneously in
multiple parameters. Note, however, that the result of the Laurent expansion depends
on the ordering of the expansion variables. The second argument of pySecDec.expansion.expand_singular()
determines the order of the expansion:
>>> expression = Expression('1/(2*eps) * 1/(eps + alpha)', ['eps', 'alpha']).simplify()
>>> eps_first = expand_singular(expression, [0,1], [1,1])
>>> eps_first
+ (( + (1/2)) * (( + (1))**(1)))*eps**1*alpha**1 + (( + (1/2)) * (( + (1))**(1)))*alpha**2 + (( + (1)) * (( + (2))**(1)))*eps*alpha**3
>>> alpha_first = expand_singular(expression, [1,0], [1,1])
>>> alpha_first
+ (( + (1/2)) * (( + (1))**(1)))*eps**2 + (( + (1/2)) * (( + (1))**(1)))*eps**3*alpha
The expression printed out by our algebra module are quite messy. In order to obtain nicer output, we can convert these expressions to the slower but more high level sympy:
>>> import sympy as sp
>>> eps_first = expand_singular(expression, [0,1], [1,1])
>>> alpha_first = expand_singular(expression, [1,0], [1,1])
>>> sp.sympify(eps_first)
1/(2*alpha*eps)  1/(2*alpha**2) + eps/(2*alpha**3)
>>> sp.sympify(alpha_first)
alpha/(2*eps**3) + 1/(2*eps**2)
SecDecUtil¶
SecDecUtil is a standalone autotoolsc++ package, that collects common helper classes and
functions needed by the c++ code generated using loop_package
or make_package
. Everything defined by the SecDecUtil
is put into the c++ namepace secdecutil.
Series¶
A class template for containing (optionally truncated) Laurent series. Multivariate series can be represented as series of series.
This class overloads the arithmetic operators (+
, 
, *
, /
) and the comparator operators (==
, !=
).
A string representation can be obtained using the <<
operator.
The at(i)
and [i]
operators return the coefficient of the i
^{th} power of the expansion parameter. Otherwise elements can be accessed identically to std::vector
.
 template<typename
T
>
classSeries
¶
 std::string
expansion_parameter
¶A string representing the expansion parameter of the series (default
x
)
 int
get_order_min
() const¶Returns the lowest order in the series.
 int
get_order_max
() const¶Returns the highest order in the series.
 bool
get_truncated_above
() const¶Checks whether the series is truncated from above.
 bool
has_term
(int order) const¶Checks whether the series has a term at order
order
.
Example:
#include <iostream>
#include <secdecutil/series.hpp>
int main()
{
secdecutil::Series<int> exact(2,1,{1,2,3,4},false,"eps");
secdecutil::Series<int> truncated(2,1,{1,2,3,4},true,"eps");
secdecutil::Series<secdecutil::Series<int>> multivariate(1,2,
{
{2,1,{1,2},false,"alpha"},
{2,1,{3,4},false,"alpha"},
},false,"eps"
);
std::cout << "exact: " << exact << std::endl;
std::cout << "truncated: " << truncated << std::endl;
std::cout << "multivariate: " << multivariate << std::endl << std::endl;
std::cout << "exact + 1: " << exact + 1 << std::endl;
std::cout << "exact * exact: " << exact * exact << std::endl;
std::cout << "exact * truncated: " << exact * truncated << std::endl;
std::cout << "exact.at(2): " << exact.at(2) << std::endl;
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include std=c++11 example.cpp o example lm && ./example
Output:
exact: + (1)*eps^2 + (2)*eps^1 + (3) + (4)*eps
truncated: + (1)*eps^2 + (2)*eps^1 + (3) + (4)*eps + O(eps^2)
multivariate: + ( + (1)*alpha^2 + (2)*alpha^1)*eps + ( + (3)*alpha^2 + (4)*alpha^1)*eps^2
exact + 1: + (1)*eps^2 + (2)*eps^1 + (4) + (4)*eps
exact * exact: + (1)*eps^4 + (4)*eps^3 + (10)*eps^2 + (20)*eps^1 + (25) + (24)*eps + (16)*eps^2
exact * truncated: + (1)*eps^4 + (4)*eps^3 + (10)*eps^2 + (20)*eps^1 + O(eps^0)
exact.at(2): 1
Deep Apply¶
A general concept to apply a std::function
to a nested data structure. If the applied std::function
is not void then deep_apply()
returns a nested data structure of the return values. Currently secdecutil implements this for std::vector
and Series
.
This concept allows, for example, the elements of a nested series to be edited without knowing the depth of the nested structure.
Example (complex conjugate a Series
):
#include <iostream>
#include <complex>
#include <secdecutil/series.hpp>
#include <secdecutil/deep_apply.hpp>
int main()
{
std::function<std::complex<double>(std::complex<double>)> conjugate =
[] (std::complex<double> element)
{
return std::conj(element);
};
secdecutil::Series<std::complex<double>> u(1,0,{{1,2},{3,4}},false,"eps");
secdecutil::Series<secdecutil::Series<std::complex<double>>> m(1,1,{{1,1,{{1,2}},false,"alpha"},},false,"eps");
std::cout << "u: " << u << std::endl;
std::cout << "m: " << m << std::endl << std::endl;
std::cout << "conjugated u: " << secdecutil::deep_apply(u, conjugate) << std::endl;
std::cout << "conjugated m: " << secdecutil::deep_apply(m, conjugate) << std::endl;
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include std=c++11 example.cpp o example lm && ./example
Output:
u: + ((1,2))*eps^1 + ((3,4))
m: + ( + ((1,2))*alpha)*eps
conjugated u: + ((1,2))*eps^1 + ((3,4))
conjugated m: + ( + ((1,2))*alpha)*eps
Uncertainties¶
A class template which implements uncertainty propagation for uncorrelated random variables by overloads of the +
, 
, *
and partially /
.
Division by UncorrelatedDeviation
is not implemented as it is not always defined. It has special overloads for std::complex<T>
.
Note
Division by UncorrelatedDeviation
is not implemented as this operation is not always well defined.
Specifically, it is ill defined in the case that the errors are Gaussian distributed as the expectation value,
where
is undefined in the Riemann or Lebesgue sense. The rule can not be derived from the first principles of probability theory.
The rules implemented for real valued error propagation are:
For complex numbers the above rules are implemented for the real and imaginary parts individually.
The expectation value.
The standard deviation.
Example:
#include <iostream>
#include <complex>
#include <secdecutil/uncertainties.hpp>
int main()
{
secdecutil::UncorrelatedDeviation<double> r(1.,0.5);
secdecutil::UncorrelatedDeviation<std::complex<double>> c({2.,3.},{0.6,0.7});
std::cout << "r: " << r << std::endl;
std::cout << "c: " << c << std::endl << std::endl;
std::cout << "r.value: " << r.value << std::endl;
std::cout << "r.uncertainty: " << r.uncertainty << std::endl;
std::cout << "r + c: " << r + c << std::endl;
std::cout << "r * c: " << r * c << std::endl;
std::cout << "r / 3.0: " << r / 3. << std::endl;
// std::cout << "1. / r: " << 1. / r << std::endl; // ERROR
// std::cout << "c / r: " << c / r << std::endl; // ERROR
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include std=c++11 example.cpp o example lm && ./example
Output:
r: 1 +/ 0.5
c: (2,3) +/ (0.6,0.7)
r.value: 1
r.uncertainty: 0.5
r + c: (3,3) +/ (0.781025,0.7)
r * c: (2,3) +/ (1.20416,1.69189)
r / 3.0: 0.333333 +/ 0.166667
Integrand Container¶
A class template for containing integrands. It stores the number of integration variables and the integrand as a std::function
.
This class overloads the arithmetic operators (+
, 
, *
, /
) and the call operator (()
).
Example (add two IntegrandContainer
and evaluate one point):
#include <iostream>
#include <secdecutil/integrand_container.hpp>
int main()
{
using input_t = const double * const;
using return_t = double;
std::function<return_t(input_t)> f1 = [] (input_t x) { return 2*x[0]; };
secdecutil::IntegrandContainer<return_t,input_t> c1(1,f1);
std::function<return_t(input_t)> f2 = [] (input_t x) { return x[0]*x[1]; };
secdecutil::IntegrandContainer<return_t,input_t> c2(2,f2);
secdecutil::IntegrandContainer<return_t,input_t> c3 = c1 + c2;
const double point[]{1.0,2.0};
std::cout << "c1.number_of_integration_variables: " << c1.number_of_integration_variables << std::endl;
std::cout << "c2.number_of_integration_variables: " << c2.number_of_integration_variables << std::endl << std::endl;
std::cout << "c3.number_of_integration_variables: " << c3.number_of_integration_variables << std::endl;
std::cout << "c3.integrand(point): " << c3.integrand(point) << std::endl;
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include std=c++11 example.cpp o example lm && ./example
Output:
c1.number_of_integration_variables: 1
c2.number_of_integration_variables: 2
c3.number_of_integration_variables: 2
c3.integrand(point): 4
Integrator¶
A base class template from which integrator implementations inherit. It defines the minimal API available for all integrators.
 template<typename
return_t
, typenameinput_t
, typenamecontainer_t
= secdecutil::IntegrandContainer<return_t, input_t const *const >>
classIntegrator
¶
 bool
together
¶(Only available if
return_t
is astd::complex
type) Iftrue
after each call of the function both the real and imaginary parts are passed to the underlying integrator. Iffalse
after each call of the function only the real or imaginary part is passed to the underlying integrator. For some adaptive integrators considering the real and imaginary part of a complex function separately can improve the sampling. Default:false
.
 UncorrelatedDeviation<return_t>
integrate
(const IntegrandContainer<return_t, input_t const *const >&)¶Integrates the
IntegrandContainer
and returns the value and uncertainty as anUncorrelatedDeviation
.
An integrator that chooses another integrator based on the dimension of the integrand.
 template<typename
return_t
, typenameinput_t
>
classMultiIntegrator
¶
 Integrator<return_t, input_t> &
low_dimensional_integrator
¶Reference to the integrator to be used if the integrand has a lower dimension than
critical_dim
.
 Integrator<return_t, input_t> &
high_dimensional_integrator
¶Reference to the integrator to be used if the integrand has dimension
critical_dim
or higher.
 int
critical_dim
¶The dimension below which the
low_dimensional_integrator
is used.
CQuad¶
For one dimensional integrals, we wrap the cquad integrator form the GNU scientifc library (gsl).
 CQuad takes the following options:
epsrel
 The desired relative accuracy for the numerical evaluation. Default:0.01
.epsabs
 The desired absolute accuracy for the numerical evaluation. Default:1e7
.n
 The size of the workspace. This value can only be set in the constructor. Changing this attribute of an instance is not possible. Default:100
.verbose
 Whether or not to print status information. Default:false
.zero_border
 The minimal value an integration variable can take. Default:0.0
. (new in version 1.3)
Qmc¶
The quasimonte carlo integrator as described in [PSD18]. Using a quasimonte integrator to compute sector decomposed integrals was pioneered in [LWY+15].

template<typename
return_t
, integrators::Umaxdim
, template<typename, typename, integrators::U> classtransform_t
, typenamecontainer_t
= secdecutil::IntegrandContainer<return_t, typename remove_complex<return_t>::type const *const >, template<typename, typename, integrators::U> classfitfunction_t
= void_template>
classQmc
: Integrator<return_t, return_t, container_t>, public integrators::Qmc<return_t, return_t, maxdim, transform_t, fitfunction_t>¶ Derived from
secdecutil::Integrator
and::integrators::Qmc
 the underlying standalone implementation of the Qmc.
 The most important fields and template argments of
Qmc
are: minn
 The minimal number of points in the Qmc lattice. Will be augmented to the next larger availablen
.minm
 The minimal number of random shifts.maxeval
 The maximal number of integrand evaluations.epsrel
 The desired relative accuracy for the numerical evaluation.epsabs
 The desired absolute accuracy for the numerical evaluation.maxdim
 The highest dimension theQmc
instance can be used for.transform_t
 The periodizing transform to apply prior to integration.fitfunction_t
 The fit function transform to apply for adaptive integration.verbosity
 Controls the amount of status messages during integration. Can be0
,1
,2
, or3
.devices
 Astd::set
of devices to run on.1
denotes the CPU, positive integers refer to GPUs.
Refer to the documentation of the standalone Qmc for the default values and additional information.
An integral transform has to be chosen by setting the template argument transform_t
. Available transforms
are e.g. Korobov<r0,r1>
and Sidi<r0>
, please refer to the underlying Qmc implementation for a complete list.
The fit function for adaptive integration can be set by the fitfunction_t
, e.g. PolySingular
. If not set,
the default of the underlying Qmc implementation is used.
Examples how to use the Qmc on the CPU and on both, CPU and GPU are shown below.
Cuba¶
 Currently we wrap the following Cuba integrators:
Vegas
Suave
Divonne
Cuhre
 The Cuba integrators all implement:
epsrel
 The desired relative accuracy for the numerical evaluation. Default:0.01
.epsabs
 The desired absolute accuracy for the numerical evaluation. Default:1e7
.flags
 Sets the Cuba verbosity flags. Theflags=2
means that the Cuba input parameters and the result after each iteration are written to the log file of the numerical integration. Default:0
.seed
 The seed used to generate random numbers for the numerical integration with Cuba. Default:0
.mineval
 The number of evaluations which should at least be done before the numerical integrator returns a result. Default:0
.maxeval
 The maximal number of evaluations to be performed by the numerical integrator. Default:1000000
.zero_border
 The minimal value an integration variable can take. Default:0.0
. (new in version 1.3)
The available integrator specific parameters and their default values are:
Vegas  Suave  Divonne  Cuhre 

nstart (1000 ) 
nnew (1000 ) 
key1 (2000 ) 
key (0 ) 
nincrease (500 ) 
nmin (10 ) 
key2 (1 ) 

nbatch (500 ) 
flatness (25.0 ) 
key3 (1 ) 

maxpass (4 ) 

border (0.0 ) 

maxchisq (1.0 ) 

mindeviation (0.15 ) 
For the description of these more specific parameters we refer to the Cuba manual.
Examples¶
Integrate Real Function with Cuba Vegas¶
Example:
#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/cuba.hpp>
int main()
{
using input_t = const double * const;
using return_t = double;
secdecutil::cuba::Vegas<return_t> integrator;
integrator.epsrel = 1e4;
integrator.maxeval = 1e7;
secdecutil::IntegrandContainer<return_t,input_t> c(2, [] (input_t x) { return x[0]*x[1]; });
secdecutil::UncorrelatedDeviation<return_t> result = integrator.integrate(c);
std::cout << "result: " << result << std::endl;
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include L${SECDEC_CONTRIB}/lib std=c++11 example.cpp o example lcuba lm && ./example
Output:
result: 0.250002 +/ 2.4515e05
Integrate Complex Function with Cuba Vegas¶
Example:
#include <iostream>
#include <complex>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/cuba.hpp>
int main()
{
using input_t = const double * const;
using return_t = std::complex<double>;
secdecutil::cuba::Vegas<return_t> integrator;
std::function<return_t(input_t)> f = [] (input_t x) { return return_t{x[0],x[1]}; };
secdecutil::IntegrandContainer<return_t,input_t> c(2,f);
integrator.together = false; // integrate real and imaginary part separately (default)
secdecutil::UncorrelatedDeviation<return_t> result_separate = integrator.integrate(c);
integrator.together = true; // integrate real and imaginary part simultaneously
secdecutil::UncorrelatedDeviation<return_t> result_together = integrator.integrate(c);
std::cout << "result_separate: " << result_separate << std::endl;
std::cout << "result_together: " << result_together << std::endl;
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include L${SECDEC_CONTRIB}/lib std=c++11 example.cpp o example lcuba lm && ./example
Output:
result_separate: (0.499889,0.500284) +/ (0.00307225,0.00305688)
result_together: (0.499924,0.500071) +/ (0.00357737,0.00357368)
Integrate Real Function with Cuba Vegas or CQuad¶
Example:
#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/integrator.hpp>
#include <secdecutil/integrators/cuba.hpp>
#include <secdecutil/integrators/cquad.hpp>
int main()
{
using input_base_t = double;
using input_t = const input_base_t * const;
using return_t = double;
secdecutil::cuba::Vegas<return_t> vegas;
vegas.epsrel = 1e5;
vegas.maxeval = 1e7;
secdecutil::gsl::CQuad<return_t> cquad;
cquad.epsrel = 1e10;
cquad.epsabs = 1e13;
secdecutil::MultiIntegrator<return_t,input_base_t> integrator(cquad,vegas,2);
secdecutil::IntegrandContainer<return_t,input_t> one_dimensional(1, [] (input_t x) { return x[0]; });
secdecutil::IntegrandContainer<return_t,input_t> two_dimensional(2, [] (input_t x) { return x[0]*x[1]; });
secdecutil::UncorrelatedDeviation<return_t> result_1d = integrator.integrate(one_dimensional); // uses cquad
secdecutil::UncorrelatedDeviation<return_t> result_2d = integrator.integrate(two_dimensional); // uses vegas
std::cout << "result_1d: " << result_1d << std::endl;
std::cout << "result_2d: " << result_2d << std::endl;
}
Compile/Run:
$ c++ I${SECDEC_CONTRIB}/include L${SECDEC_CONTRIB}/lib std=c++11 example.cpp o example lcuba lgsl lgslcblas lm && ./example
Output:
result_1d: 0.5 +/ 9.58209e17
result_2d: 0.25 +/ 5.28953e06
Set the integral transform of the Qmc¶
Example:
#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/qmc.hpp>
using input_base_t = double;
using input_t = input_base_t const * const;
using return_t = double;
using container_t = secdecutil::IntegrandContainer<return_t,input_t>;
using result_t = secdecutil::UncorrelatedDeviation<return_t>;
const int seed = 12345, maxdim = 4;
int main()
{
/*
* minimal instantiation
*/
secdecutil::integrators::Qmc
<
return_t, // the return type of the integrand
maxdim, // the highest dimension this integrator will be used for
::integrators::transforms::Baker::type // the integral transform
> integrator_baker;
integrator_baker.randomgenerator.seed(seed);
/*
* disable adaptation
*/
secdecutil::integrators::Qmc
<
return_t, // the return type of the integrand
maxdim, // the highest dimension this integrator will be used for
::integrators::transforms::Korobov<4,1>::type, // the integral transform
container_t, // the functor type to be passed to this integrator
::integrators::fitfunctions::None::type // the fit funtion
> integrator_korobov4x1;
integrator_korobov4x1.randomgenerator.seed(seed);
/*
* enable adaptation
*/
secdecutil::integrators::Qmc
<
return_t, // the return type of the integrand
maxdim, // the highest dimension this integrator will be used for
::integrators::transforms::Sidi<3>::type, // the integral transform
container_t, // the functor type to be passed to this integrator
::integrators::fitfunctions::PolySingular::type // the fit funtion
> integrator_sidi3_adaptive;
integrator_sidi3_adaptive.randomgenerator.seed(seed);
// define the integrand as a functor
container_t integrand(
4, // dimension
[] (input_t x) { return x[0]*x[1]*x[2]*x[3]; } // integrand function
);
// compute the integral with different settings
result_t result_baker = integrator_baker.integrate(integrand);
result_t result_korobov4x1 = integrator_korobov4x1.integrate(integrand);
result_t result_sidi3_adaptive = integrator_sidi3_adaptive.integrate(integrand);
// print the results
std::cout << "baker: " << result_baker << std::endl;
std::cout << "Korobov (weights 4, 1): " << result_korobov4x1 << std::endl;
std::cout << "Sidi (weight 3, adaptive): " << result_sidi3_adaptive << std::endl;
}
Compile/Run:
c++ I${SECDEC_CONTRIB}/include pthread L${SECDEC_CONTRIB}/lib std=c++11 example.cpp o example lm && ./example
Output:
baker: 0.0625 +/ 7.93855e08
Korobov (weights 4, 1): 0.0625108 +/ 2.97931e05
Sidi (weight 3, adaptive): 0.0625 +/ 4.33953e09
Run the Qmc on GPUs¶
Example:
#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/qmc.hpp>
using input_base_t = double;
using input_t = input_base_t const * const;
using return_t = double;
using container_t = secdecutil::IntegrandContainer<return_t,input_t>;
using result_t = secdecutil::UncorrelatedDeviation<return_t>;
/*
* `container_t` cannot be used on the GPU > define a different container type
*/
struct cuda_integrand_t
{
const int number_of_integration_variables = 4;
// integrand function
#ifdef __CUDACC__
__host__ __device__
#endif
return_t operator()(input_t x)
{
return x[0]*x[1]*x[2]*x[3];
};
} cuda_integrand;
const int seed = 12345, maxdim = 4;
int main()
{
/*
* Qmc capable of sampling on the GPU
*/
secdecutil::integrators::Qmc
<
return_t, // the return type of the integrand
maxdim, // the highest dimension this integrator will be used for
::integrators::transforms::Sidi<3>::type, // the integral transform
cuda_integrand_t, // the functor type to be passed to this integrator
::integrators::fitfunctions::PolySingular::type // the fit funtion (optional)
> integrator_sidi3_adaptive_gpu;
integrator_sidi3_adaptive_gpu.randomgenerator.seed(seed);
// compute the integral with different settings
result_t result_sidi3_adaptive_gpu = integrator_sidi3_adaptive_gpu.integrate(cuda_integrand);
// print the results
std::cout << "Sidi (weight 3, adaptive): " << result_sidi3_adaptive_gpu << std::endl;
}
Compile/Run:
nvcc x cu I${SECDEC_CONTRIB}/include L${SECDEC_CONTRIB}/lib std=c++11 example.cpp o example lgsl lgslcblas lm && ./example # with GPU
c++ I${SECDEC_CONTRIB}/include pthread L${SECDEC_CONTRIB}/lib std=c++11 example.cpp o example lgsl lgslcblas lm && ./example # without GPU
Output:
Sidi (weight 3, adaptive): 0.0625 +/ 4.33953e09
Reference Guide¶
This section describes all public functions and classes in pySecDec.
Algebra¶
Implementation of a simple computer algebra system.

class
pySecDec.algebra.
ExponentiatedPolynomial
(expolist, coeffs, exponent=1, polysymbols='x', copy=True)¶ Like
Polynomial
, but with a global exponent.Parameters:  expolist – iterable of iterables; The variable’s powers for each term.
 coeffs – iterable; The coefficients of the polynomial.
 exponent – object, optional; The global exponent.
 polysymbols –
iterable or string, optional; The symbols to be used for the polynomial variables when converted to string. If a string is passed, the variables will be consecutively numbered.
 For example: expolist=[[2,0],[1,1]] coeffs=[“A”,”B”]
 polysymbols=’x’ (default) <> “A*x0**2 + B*x0*x1”
 polysymbols=[‘x’,’y’] <> “A*x**2 + B*x*y”
 copy –
bool; Whether or not to copy the expolist, the coeffs, and the exponent.
Note
If copy is
False
, it is assumed that the expolist, the coeffs and the exponent have the correct type.

copy
()¶ Return a copy of a
Polynomial
or a subclass.

derive
(index)¶ Generate the derivative by the parameter indexed index.
Parameters: index – integer; The index of the paramater to derive by.

simplify
()¶ Apply the identity <something>**0 = 1 or <something>**1 = <something> or 1**<something> = 1 if possible, otherwise call the simplify method of the base class. Convert
exponent
to symbol if possible.

pySecDec.algebra.
Expression
(expression, polysymbols, follow_functions=False)¶ Convert a sympy expression to an expression in terms of this module.
Parameters:  expression – string or sympy expression; The expression to be converted
 polysymbols – iterable of strings or sympy symbols;
The symbols to be stored as
expolists
(seePolynomial
) where possible.  follow_functions – bool, optional (default =
False
); If true, return the converted expression and a list ofFunction
that occur in the expression.

class
pySecDec.algebra.
Function
(symbol, *arguments, **kwargs)¶ Symbolic function that can take care of parameter transformations. It keeps track of all taken derivatives: When
derive()
is called, save the multiindex of the taken derivative.The derivative multiindices are the keys in the dictionary
self.derivative_tracks
. The values are lists with two elements: Its first element is the index to derive the derivative indicated by the multiindex in the second element by, in order to abtain the derivative indicated by the key:>>> from pySecDec.algebra import Polynomial, Function >>> x = Polynomial.from_expression('x', ['x','y']) >>> y = Polynomial.from_expression('y', ['x','y']) >>> poly = x**2*y + y**2 >>> func = Function('f', x, y) >>> ddfuncd0d1 = func.derive(0).derive(1) >>> func Function(f( + (1)*x, + (1)*y), derivative_tracks = {(1, 0): [0, (0, 0)], (1, 1): [1, (1, 0)]}) >>> func.derivative_tracks {(1, 0): [0, (0, 0)], (1, 1): [1, (1, 0)]} >>> func.compute_derivatives(poly) {(1, 0): + (2)*x*y, (1, 1): + (2)*x}
Parameters:  symbol – string; The symbol to be used to represent the Function.
 arguments – arbitrarily many
_Expression
; The arguments of the Function.  copy – bool; Whether or not to copy the arguments.

compute_derivatives
(expression=None)¶ Compute all derivatives of
expression
that are mentioned inself.derivative_tracks
. The purpose of this function is to avoid computing the same derivatives multiple times.Parameters: expression – _Expression
, optional; The expression to compute the derivatives of. If not provided, the derivatives are shown as in terms of the function’s derivativesdfd<index>
.

derive
(index)¶ Generate the derivative by the parameter indexed index. The derivative of a function with symbol
f
by some index is denoted asdfd<index>
.Parameters: index – integer; The index of the paramater to derive by.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –

simplify
()¶ Simplify the arguments.

class
pySecDec.algebra.
Log
(arg, copy=True)¶ The (natural) logarithm to base e (2.718281828459..). Store the expressions
log(arg)
.Parameters:  arg –
_Expression
; The argument of the logarithm.  copy – bool; Whether or not to copy the arg.

derive
(index)¶ Generate the derivative by the parameter indexed index.
Parameters: index – integer; The index of the paramater to derive by.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –

simplify
()¶ Apply
log(1) = 0
.
 arg –

class
pySecDec.algebra.
LogOfPolynomial
(expolist, coeffs, polysymbols='x', copy=True)¶ The natural logarithm of a
Polynomial
.Parameters:  expolist – iterable of iterables; The variable’s powers for each term.
 coeffs – iterable; The coefficients of the polynomial.
 exponent – object, optional; The global exponent.
 polysymbols –
iterable or string, optional; The symbols to be used for the polynomial variables when converted to string. If a string is passed, the variables will be consecutively numbered.
 For example: expolist=[[2,0],[1,1]] coeffs=[“A”,”B”]
 polysymbols=’x’ (default) <> “A*x0**2 + B*x0*x1”
 polysymbols=[‘x’,’y’] <> “A*x**2 + B*x*y”

derive
(index)¶ Generate the derivative by the parameter indexed index.
Parameters: index – integer; The index of the paramater to derive by.

static
from_expression
(expression, polysymbols)¶ Alternative constructor. Construct the
LogOfPolynomial
from an algebraic expression.Parameters:  expression – string or sympy expression; The algebraic representation of the polynomial, e.g. “5*x1**2 + x1*x2”
 polysymbols – iterable of strings or sympy symbols; The symbols to be interpreted as the polynomial variables, e.g. “[‘x1’,’x2’]”.

simplify
()¶ Apply the identity
log(1) = 0
, otherwise call the simplify method of the base class.

class
pySecDec.algebra.
Polynomial
(expolist, coeffs, polysymbols='x', copy=True)¶ Container class for polynomials. Store a polynomial as list of lists counting the powers of the variables. For example the polynomial “x1**2 + x1*x2” is stored as [[2,0],[1,1]].
Coefficients are stored in a separate list of strings, e.g. “A*x0**2 + B*x0*x1” <> [[2,0],[1,1]] and [“A”,”B”].
Parameters:  expolist –
iterable of iterables; The variable’s powers for each term.
Hint
Negative powers are allowed.
 coeffs – 1d arraylike with numerical or sympysymbolic (see http://www.sympy.org/) content, e.g. [x,1,2] where x is a sympy symbol; The coefficients of the polynomial.
 polysymbols –
iterable or string, optional; The symbols to be used for the polynomial variables when converted to string. If a string is passed, the variables will be consecutively numbered.
 For example: expolist=[[2,0],[1,1]] coeffs=[“A”,”B”]
 polysymbols=’x’ (default) <> “A*x0**2 + B*x0*x1”
 polysymbols=[‘x’,’y’] <> “A*x**2 + B*x*y”
 copy –
bool; Whether or not to copy the expolist and the coeffs.
Note
If copy is
False
, it is assumed that the expolist and the coeffs have the correct type.

becomes_zero_for
(zero_params)¶ Return True if the polynomial becomes zero if the parameters passed in zero_params are set to zero. Otherwise, return False.
Parameters: zero_params – iterable of integers; The indices of the parameters to be checked.

copy
()¶ Return a copy of a
Polynomial
or a subclass.

derive
(index)¶ Generate the derivative by the parameter indexed index.
Parameters: index – integer; The index of the paramater to derive by.

static
from_expression
(expression, polysymbols)¶ Alternative constructor. Construct the polynomial from an algebraic expression.
Parameters:  expression – string or sympy expression; The algebraic representation of the polynomial, e.g. “5*x1**2 + x1*x2”
 polysymbols – iterable of strings or sympy symbols; The symbols to be interpreted as the polynomial variables, e.g. “[‘x1’,’x2’]”.

has_constant_term
(indices=None)¶ Return True if the polynomial can be written as:
Otherwise, return False.
Parameters: indices – list of integers or None; The indices of the polysymbols to consider. If None
(default) all indices are taken into account.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –

simplify
(deep=True)¶ Combine terms that have the same exponents of the variables.
Parameters: deep – bool; If True
(default) call the simplify method of the coefficients if they are of type_Expression
.
 expolist –

class
pySecDec.algebra.
Pow
(base, exponent, copy=True)¶ Exponential. Store two expressions
A
andB
to be interpreted as the exponentialA**B
.Parameters:  base –
_Expression
; The baseA
of the exponential.  exponent –
_Expression
; The exponentB
.  copy – bool; Whether or not to copy base and exponent.

derive
(index)¶ Generate the derivative by the parameter indexed index.
Parameters: index – integer; The index of the paramater to derive by.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –

simplify
()¶ Apply the identity <something>**0 = 1 or <something>**1 = <something> or 1**<something> = 1 if possible. Convert to
ExponentiatedPolynomial
orPolynomial
if possible.
 base –

class
pySecDec.algebra.
Product
(*factors, **kwargs)¶ Product of polynomials. Store one or polynomials to be interpreted as product .
Parameters:  factors – arbitrarily many instances of
Polynomial
; The factors .  copy – bool; Whether or not to copy the factors.
can be accessed with
self.factors[i]
.Example:
p = Product(p0, p1) p0 = p.factors[0] p1 = p.factors[1]

derive
(index)¶ Generate the derivative by the parameter indexed index. Return an instance of the optimized
ProductRule
.Parameters: index – integer; The index of the paramater to derive by.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –
 factors – arbitrarily many instances of

class
pySecDec.algebra.
ProductRule
(*expressions, **kwargs)¶ Store an expression of the form
The main reason for introducing this class is a speedup when calculating derivatives. In particular, this class implements simplifications such that the number of terms grows less than exponentially (scaling of the naive implementation of the product rule) with the number of derivatives.
Parameters: expressions – arbitrarily many expressions; The expressions . 
copy
()¶ Return a copy of a
ProductRule
.

derive
(index)¶ Generate the derivative by the parameter indexed index. Note that this class is particularly designed to hold derivatives of a product.
Parameters: index – integer; The index of the paramater to derive by.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –

simplify
()¶ Combine terms that have the same derivatives of the expressions.

to_sum
()¶ Convert the
ProductRule
toSum


class
pySecDec.algebra.
Sum
(*summands, **kwargs)¶ Sum of polynomials. Store one or polynomials to be interpreted as product .
Parameters:  summands – arbitrarily many instances of
Polynomial
; The summands .  copy – bool; Whether or not to copy the summands.
can be accessed with
self.summands[i]
.Example:
p = Sum(p0, p1) p0 = p.summands[0] p1 = p.summands[1]

derive
(index)¶ Generate the derivative by the parameter indexed index.
Parameters: index – integer; The index of the paramater to derive by.

replace
(index, value, remove=False)¶ Replace a variable in an expression by a number or a symbol. The entries in all
expolist
of the underlyingPolynomial
are set to zero. The coefficients are modified according to value and the powers indicated in theexpolist
.Parameters:  expression –
_Expression
; The expression to replace the variable.  index – integer; The index of the variable to be replaced.
 value – number or sympy expression; The value to insert for the chosen variable.
 remove – bool;
Whether or not to remove the replaced
parameter from the
parameters
in the expression.
 expression –
 summands – arbitrarily many instances of
Loop Integral¶
This module defines routines to Feynman parametrize a loop integral and build a c++ package that numerically integrates over the sector decomposed integrand.
Feynman Parametrization¶
Routines to Feynman parametrize a loop integral.

class
pySecDec.loop_integral.
LoopIntegral
(*args, **kwargs)¶ Container class for loop integrals. The main purpose of this class is to convert a loop integral from the momentum representation to the Feynman parameter representation.
It is possible to provide either the graph of the loop integrals as adjacency list, or the propagators.
The Feynman parametrized integral is a product of the following expressions that are accessible as member properties:
self.regulator ** self.regulator_power
self.Gamma_factor
self.exponentiated_U
self.exponentiated_F
self.numerator
self.measure
,
where
self
is an instance of eitherLoopIntegralFromGraph
orLoopIntegralFromPropagators
.When inverse propagators or nonnumerical propagator powers are present (see powerlist), some Feynman_parameters drop out of the integral. The variables to integrate over can be accessed as
self.integration_variables
.While
self.numerator
describes the numerator polynomial generated by tensor numerators or inverse propagators,self.measure
contains the monomial associated with the integration measure in the case of propagator powers . The Gamma functions in the denominator belonging to the measure, however, are multiplied to the overall Gamma factor given byself.Gamma_factor
.Changed in version 1.2.2: The overall sign is included in
self.Gamma_factor
.See also
 input as graph:
LoopIntegralFromGraph
 input as list of propagators:
LoopIntegralFromPropagators

class
pySecDec.loop_integral.
LoopIntegralFromGraph
(internal_lines, external_lines, replacement_rules=[], Feynman_parameters='x', regulator='eps', regulator_power=0, dimensionality='42*eps', powerlist=[])¶ Construct the Feynman parametrization of a loop integral from the graph using the cut construction method.
Example:
>>> from pySecDec.loop_integral import * >>> internal_lines = [['0',[1,2]], ['m',[2,3]], ['m',[3,1]]] >>> external_lines = [['p1',1],['p2',2],['p12',3]] >>> li = LoopIntegralFromGraph(internal_lines, external_lines) >>> li.exponentiated_U ( + (1)*x0 + (1)*x1 + (1)*x2)**(2*eps  1) >>> li.exponentiated_F ( + (m**2)*x2**2 + (2*m**2  p12**2)*x1*x2 + (m**2)*x1**2 + (m**2  p1**2)*x0*x2 + (m**2  p2**2)*x0*x1)**(eps  1)
Parameters:  internal_lines – iterable of internal line specification, consisting
of string or sympy expression for mass and a pair
of strings or numbers for the vertices, e.g.
[['m', [1,2]], ['0', [2,1]]]
.  external_lines – iterable of external line specification, consisting
of string or sympy expression for external momentum
and a strings or number for the vertex, e.g.
[['p1', 1], ['p2', 2]]
.  replacement_rules – iterable of iterables with two strings or sympy
expressions, optional;
Symbolic replacements to be made for the external
momenta, e.g. definition of Mandelstam variables.
Example:
[('p1*p2', 's'), ('p1**2', 0)]
wherep1
andp2
are external momenta. It is also possible to specify vector replacements, for example[('p4', '(p1+p2+p3)')]
.  Feynman_parameters – iterable or string, optional; The symbols to be used for the Feynman parameters. If a string is passed, the Feynman parameter variables will be consecutively numbered starting from zero.
 regulator –
string or sympy symbol, optional; The symbol to be used for the dimensional regulator (typically or )
Note
If you change this symbol, you have to adapt the dimensionality accordingly.
 regulator_power –
integer; An additional factor to the numerator.
See also
 dimensionality – string or sympy expression, optional; The dimensionality; typically , which is the default value.
 powerlist – iterable, optional; The powers of the propergators, possibly dependent on the regulator. In case of negative powers, the numerator is constructed by taking derivatives with respect to the corresponding Feynman parameters as explained in Section 3.2.4 of Ref. [BHJ+15]. If negative powers are combined with a tensor numerator, the derivatives act on the Feynmanparametrized tensor numerator as well, which leads to a consistent result.
 internal_lines – iterable of internal line specification, consisting
of string or sympy expression for mass and a pair
of strings or numbers for the vertices, e.g.

class
pySecDec.loop_integral.
LoopIntegralFromPropagators
(propagators, loop_momenta, external_momenta=[], Lorentz_indices=[], numerator=1, metric_tensor='g', replacement_rules=[], Feynman_parameters='x', regulator='eps', regulator_power=0, dimensionality='42*eps', powerlist=[])¶ Construct the Feynman parametrization of a loop integral from the algebraic momentum representation.
Example:
>>> from pySecDec.loop_integral import * >>> propagators = ['k**2', '(k  p)**2'] >>> loop_momenta = ['k'] >>> li = LoopIntegralFromPropagators(propagators, loop_momenta) >>> li.exponentiated_U ( + (1)*x0 + (1)*x1)**(2*eps  2) >>> li.exponentiated_F ( + (p**2)*x0*x1)**(eps)
The 1st (U) and 2nd (F) Symanzik polynomials and their exponents can also be accessed independently:
>>> li.U + (1)*x0 + (1)*x1 >>> li.F + (p**2)*x0*x1 >>> >>> li.exponent_U 2*eps  2 >>> li.exponent_F eps
Parameters:  propagators – iterable of strings or sympy expressions;
The propagators, e.g.
['k1**2', '(k1k2)**2  m1**2']
.  loop_momenta – iterable of strings or sympy expressions;
The loop momenta, e.g.
['k1','k2']
.  external_momenta –
iterable of strings or sympy expressions, optional; The external momenta, e.g.
['p1','p2']
. Specifying the external_momenta is only required when a numerator is to be constructed.See also
parameter numerator
 Lorentz_indices –
iterable of strings or sympy expressions, optional; Symbols to be used as Lorentz indices in the numerator.
See also
parameter numerator
 numerator –
string or sympy expression, optional; The numerator of the loop integral. Scalar products must be passed in index notation e.g.
k1(mu)*k2(mu)
. The numerator must be a sum of products of exclusively: numbers
 scalar products (e.g.
p1(mu)*k1(mu)*p1(nu)*k2(nu)
)  symbols (e.g.
s
,eps
)
 Examples:
p1(mu)*k1(mu)*p1(nu)*k2(nu) + 4*s*eps*k1(mu)*k1(mu)
p1(mu)*(k1(mu) + k2(mu))*p1(nu)*k2(nu)
p1(mu)*k1(mu)
Note
In order to use the resulting LoopIntegral as an argument to the function
pySecDec.loop_integral.loop_package()
, the resulting Feynman parametrizedself.numerator
must be expressible as apySecDec.algebra.Polynomial
such that all coefficients are purely numeric. In addition, all scalar products of the external momenta must be expressed in terms of Mandelstam variables using the replacement_rules.Warning
All Lorentz indices (including the contracted ones and also including the numbers that have been used) must be explicitly defined using the parameter Lorentz_indices.
Hint
It is possible to use numbers as indices, for example
p1(mu)*p2(mu)*k1(nu)*k2(nu) = p1(1)*p2(1)*k1(2)*k2(2)
.Hint
The numerator may have uncontracted indices, e.g.
k1(mu)*k2(nu)
. If indices are left open, however, the LoopIntegral cannot be used with the package generatorpySecDec.loop_integral.loop_package()
.  metric_tensor – string or sympy symbol, optional; The symbol to be used for the (Minkowski) metric tensor .
 replacement_rules – iterable of iterables with two strings or sympy
expressions, optional;
Symbolic replacements to be made for the external
momenta, e.g. definition of Mandelstam variables.
Example:
[('p1*p2', 's'), ('p1**2', 0)]
wherep1
andp2
are external momenta. It is also possible to specify vector replacements, for example[('p4', '(p1+p2+p3)')]
.  Feynman_parameters – iterable or string, optional; The symbols to be used for the Feynman parameters. If a string is passed, the Feynman parameter variables will be consecutively numbered starting from zero.
 regulator –
string or sympy symbol, optional; The symbol to be used for the dimensional regulator (typically or )
Note
If you change this symbol, you have to adapt the dimensionality accordingly.
 regulator_power –
integer; An additional factor to the numerator.
See also
 dimensionality – string or sympy expression, optional; The dimensionality; typically , which is the default value.
 powerlist – iterable, optional; The powers of the propergators, possibly dependent on the regulator. In case of negative powers, the numerator is constructed by taking derivatives with respect to the corresponding Feynman parameters as explained in Section 3.2.4 of Ref. [BHJ+15]. If negative powers are combined with a tensor numerator, the derivatives act on the Feynmanparametrized tensor numerator as well, which leads to a consistent result.
 propagators – iterable of strings or sympy expressions;
The propagators, e.g.
Loop Package¶
This module contains the function that generates a c++ package.

pySecDec.loop_integral.
loop_package
(name, loop_integral, requested_order, real_parameters=[], complex_parameters=[], contour_deformation=True, additional_prefactor=1, form_optimization_level=2, form_work_space='500M', decomposition_method='iterative', normaliz_executable='normaliz', enforce_complex=False, split=False, ibp_power_goal=1, use_iterative_sort=True, use_light_Pak=True, use_dreadnaut=False, use_Pak=True, processes=None)¶ Decompose, subtract and expand a Feynman parametrized loop integral. Return it as c++ package.
See also
This function is a wrapper around
pySecDec.code_writer.make_package()
.See also
The generated library is described in Generated C++ Libraries.
Parameters:  name – string; The name of the c++ namespace and the output directory.
 loop_integral –
pySecDec.loop_integral.LoopIntegral
; The loop integral to be computed.  requested_orders – integer; Compute the expansion in the regulator to this order.
 real_parameters – iterable of strings or sympy symbols, optional; Parameters to be interpreted as real numbers, e.g. Mandelstam invariants and masses.
 complex_parameters – iterable of strings or sympy symbols, optional; Parameters to be interpreted as complex numbers. To use the complex mass scheme, define the masses as complex parameters.
 contour_deformation – bool, optional;
Whether or not to produce code for contour
deformation.
Default:
True
.  additional_prefactor – string or sympy expression, optional; An additional factor to be multiplied to the loop integral. It may depend on the regulator, the real_parameters, and the complex_parameters.
 form_optimization_level – integer out of the interval [0,4], optional;
The optimization level to be used in FORM.
Default:
2
.  form_work_space – string, optional;
The FORM WorkSpace. Default:
'500M'
.  decomposition_method –
string, optional; The strategy for decomposing the polynomials. The following strategies are available:
 ’iterative’ (default)
 ’geometric’
 ’geometric_ku’
Note
For ‘geometric’ and ‘geometric_ku’, the thirdparty program “normaliz” is needed. See The Geomethod and Normaliz.
 normaliz_executable – string, optional; The command to run normaliz. normaliz is only required if decomposition_method is set to ‘geometric’ or ‘geometric_ku’. Default: ‘normaliz’
 enforce_complex – bool, optional;
Whether or not the generated integrand functions
should have a complex return type even though
they might be purely real.
The return type of the integrands is automatically
complex if contour_deformation is
True
or if there are complex_parameters. In other cases, the calculation can typically be kept purely real. Most commonly, this flag is needed iflog(<negative real>)
occurs in one of the integrand functions. However, pySecDec will suggest setting this flag toTrue
in that case. Default:False
 split – bool, optional;
Whether or not to split the integration domain in
order to map singularities from to
. Set this option to
True
if you have singularties when one or more integration variables are one. Default:False
 ibp_power_goal –
number or iterable of number, optional; The power_goal that is forwarded to
integrate_by_parts()
.This option controls how the subtraction terms are generated. Setting it to
numpy.inf
disablesintegrate_by_parts()
, while0
disablesintegrate_pole_part()
.See also
To generate the subtraction terms, this function first calls
integrate_by_parts()
for each integration variable with the give ibp_power_goal. Thenintegrate_pole_part()
is called.Default:
1
 use_iterative_sort – bool;
Whether or not to use
squash_symmetry_redundant_sectors_sort()
withiterative_sort()
to find sector symmetries. Default:True
 use_light_Pak – bool;
Whether or not to use
squash_symmetry_redundant_sectors_sort()
withlight_Pak_sort()
to find sector symmetries. Default:True
 use_dreadnaut – bool or string, optional;
Whether or not to use
squash_symmetry_redundant_sectors_dreadnaut()
to find sector symmetries. If given a string, interpret that string as the command line executable dreadnaut. IfTrue
, try$SECDEC_CONTRIB/bin/dreadnaut
and, if the environment variable$SECDEC_CONTRIB
is not set,dreadnaut
. Default:False
 use_Pak – bool;
Whether or not to use
squash_symmetry_redundant_sectors_sort()
withPak_sort()
to find sector symmetries. Default:True
 processes – integer or None, optional;
The maximal number of processes to be used. If
None
, the number of CPUsmultiprocessing.cpu_count()
is used. New in version 1.3. Default:None
Drawing Feynman Diagrams¶
Use the following function to draw Feynman diagrams.

pySecDec.loop_integral.draw.
plot_diagram
(internal_lines, external_lines, filename, powerlist=None, neato='neato', extension='pdf', Gstart=0)¶ Draw a Feynman diagram using Graphviz (neato).
Thanks to Viktor Papara <papara@mpp.mpg.de> for his major contributions to this function.
Note
This function requires the command line tool neato. See also Drawing Feynman Diagrams with neato.
Warning
The target is overwritten without prompt if it exists already.
Parameters:  internal_lines – list;
Adjacency list of internal lines,
e.g.
[['m',['a',4]],['m',[4,5]], ['m',['a',5]],[0,[1,2]],[0,[4,1]],[0,[2,5]]]
 external_lines – list; Adjacency list of external lines, e.g. [[‘p1’,1],[‘p2’,2],[‘p3’,’a’]]
 filename – string; The name of the output file. The generated file gets this name plus the file extension.
 powerlist – list, optional; The powers of the propagators defined by the internal_lines.
 neato – string, default: “neato”; The shell command to call “neato”.
 extension – string, default: “pdf”; The file extension. This also defines the output format.
 Gstart – nonnegative int; The is value is passed to “neato” with the “Gstart” option. Try changing this value if the visualization looks bad.
 internal_lines – list;
Adjacency list of internal lines,
e.g.
Polytope¶
The polytope class as required by
pySecDec.decomposition.geometric
.

class
pySecDec.polytope.
Polytope
(vertices=None, facets=None)¶ Representation of a polytope defined by either its vertices or its facets. Call
complete_representation()
to translate from one to the other representation.Parameters:  vertices – two dimensional array; The polytope in vertex representation. Each row is interpreted as one vertex.
 facets –
two dimensional array; The polytope in facet representation. Each row represents one facet . A row in facets is interpreted as one normal vector with additionally the constant in the last column. The points of the polytope obey

complete_representation
(normaliz='normaliz', workdir='normaliz_tmp', keep_workdir=False)¶ Transform the vertex representation of a polytope to the facet representation or the other way round. Remove surplus entries in
self.facets
orself.vertices
.Note
This function calls the command line executable of normaliz [BIR]. See The Geomethod and Normaliz for installation and a list of tested versions.
Parameters:  normaliz – string; The shell command to run normaliz.
 workdir –
string; The directory for the communication with normaliz. A directory with the specified name will be created in the current working directory. If the specified directory name already exists, an
OSError
is raised.Note
The communication with normaliz is done via files.
 keep_workdir – bool; Whether or not to delete the workdir after execution.

vertex_incidence_lists
()¶ Return for each vertex the list of facets it lies in (as dictonary). The keys of the output dictonary are the vertices while the values are the indices of the facets in
self.facets
.

pySecDec.polytope.
convex_hull
(*polynomials)¶ Calculate the convex hull of the Minkowski sum of all polynomials in the input. The algorithm sets all coefficients to one first and then only keeps terms of the polynomial product that have coefficient 1. Return the list of these entries in the expolist of the product of all input polynomials.
Parameters: polynomials – abritrarily many instances of Polynomial
where all of these have an equal number of variables; The polynomials to calculate the convex hull for.

pySecDec.polytope.
triangulate
(cone, normaliz='normaliz', workdir='normaliz_tmp', keep_workdir=False, switch_representation=False)¶ Split a cone into simplicial cones; i.e. cones defined by exactly rays where is the dimensionality.
Note
This function calls the command line executable of normaliz [BIR]. See The Geomethod and Normaliz for installation and a list of tested versions.
Parameters:  cone – two dimensional array; The defining rays of the cone.
 normaliz – string; The shell command to run normaliz.
 workdir –
string; The directory for the communication with normaliz. A directory with the specified name will be created in the current working directory. If the specified directory name already exists, an
OSError
is raised.Note
The communication with normaliz is done via files.
 keep_workdir – bool; Whether or not to delete the workdir after execution.
 switch_representation – bool; Whether or not to switch between facet and vertex/ray representation.
Decomposition¶
The core of sector decomposition. This module implements the actual decomposition routines.
Common¶
This module collects routines that are used by multiple decompition modules.

class
pySecDec.decomposition.
Sector
(cast, other=[], Jacobian=None)¶ Container class for sectors that arise during the sector decomposition.
Parameters:  cast – iterable of
algebra.Product
or ofalgebra.Polynomial
; The polynomials to be cast to the form <monomial> * (const + …)  other – iterable of
algebra.Polynomial
, optional; All variable transformations are applied to these polynomials but it is not attempted to achieve the form <monomial> * (const + …)  Jacobian –
algebra.Polynomial
with one term, optional; The Jacobian determinant of this sector. If not provided, the according unit monomial (1*x0^0*x1^0…) is assumed.
 cast – iterable of

pySecDec.decomposition.
squash_symmetry_redundant_sectors_sort
(sectors, sort_function, indices=None)¶ Reduce a list of sectors by squashing duplicates with equal integral.
If two sectors only differ by a permutation of the polysymbols (to be interpreted as integration variables over some inteval), then the two sectors integrate to the same value. Thus we can drop one of them and count the other twice. The multiple counting of a sector is accounted for by increasing the coefficient of the Jacobian by one.
Equivalence up to permutation is established by applying the sort_function to each sector, this brings them into a canonical form. Sectors with identical canonical forms differ only by a permutation.
Note: whether all symmetries are found depends on the choice of sort_function. The sort function
pySecDec.matrix_sort.Pak_sort()
should find all symmetries whilst the sort functionspySecDec.matrix_sort.iterative_sort()
andpySecDec.matrix_sort.light_Pak_sort()
are faster but do not identify all symmetries.See also:
squash_symmetry_redundant_sectors_dreadnaut()
Example:
>>> from pySecDec.algebra import Polynomial >>> from pySecDec.decomposition import Sector >>> from pySecDec.decomposition import squash_symmetry_redundant_sectors_sort >>> from pySecDec.matrix_sort import Pak_sort >>> >>> poly = Polynomial([(0,1),(1,0)], ['a','b']) >>> swap = Polynomial([(1,0),(0,1)], ['a','b']) >>> Jacobian_poly = Polynomial([(1,0)], [3]) # three >>> Jacobian_swap = Polynomial([(0,1)], [5]) # five >>> sectors = ( ... Sector([poly],Jacobian=Jacobian_poly), ... Sector([swap],Jacobian=Jacobian_swap) ... ) >>> >>> reduced_sectors = squash_symmetry_redundant_sectors_sort(sectors, ... Pak_sort) >>> len(reduced_sectors) # symmetry x0 <> x1 1 >>> # The Jacobians are added together to account >>> # for the double counting of the sector. >>> reduced_sectors[0].Jacobian + (8)*x0
Parameters:  sectors – iterable of
Sector
; the sectors to be reduced.  sort_function –
pySecDec.matrix_sort.iterative_sort()
,pySecDec.matrix_sort.light_Pak_sort()
, orpySecDec.matrix_sort.Pak_sort()
; The function to be used for finding a canonical form of the sectors.  indices – iterable of integers, optional; The indices of the variables to consider. If not provided, all indices are taken into account.
 sectors – iterable of

pySecDec.decomposition.
squash_symmetry_redundant_sectors_dreadnaut
(sectors, indices=None, dreadnaut='dreadnaut', workdir='dreadnaut_tmp', keep_workdir=False)¶ Reduce a list of sectors by squashing duplicates with equal integral.
Each
Sector
is converted to aPolynomial
which is represented as a graph following the example of [MP+14] (v2.6 Figure 7, Isotopy of matrices).We first multiply each polynomial in the sector by a unique tag then sum the polynomials of the sector, this converts a sector to a polynomial. Next, we convert the expolist of the resulting polynomial to a graph where each unique exponent in the expolist is considered to be a different symbol. Each unique coefficient in the polynomial`s coeffs is assigned a vertex and connected to the row vertex of any term it multiplies. The external program dreadnaut is then used to bring the graph into a canonical form and provide a hash. Sectors with equivalent hashes may be identical, their canonical graphs are compared and if they are identical the sectors are combined.
Note
This function calls the command line executable of dreadnaut [MP+14]. It has been tested with dreadnaut version nauty26r7.
See also:
squash_symmetry_redundant_sectors_sort()
Parameters:  sectors – iterable of
Sector
; the sectors to be reduced.  indices – iterable of integers, optional; The indices of the variables to consider. If not provided, all indices are taken into account.
 dreadnaut – string; The shell command to run dreadnaut.
 workdir –
string; The directory for the communication with dreadnaut. A directory with the specified name will be created in the current working directory. If the specified directory name already exists, an
OSError
is raised.Note
The communication with dreadnaut is done via files.
 keep_workdir – bool; Whether or not to delete the workdir after execution.
 sectors – iterable of
Iterative¶
The iterative sector decomposition routines.

exception
pySecDec.decomposition.iterative.
EndOfDecomposition
¶ This exception is raised if the function
iteration_step()
is called although the sector is already in standard form.

pySecDec.decomposition.iterative.
find_singular_set
(sector, indices=None)¶ Function within the iterative sector decomposition procedure which heuristically chooses an optimal decomposition set. The strategy was introduced in arXiv:hepph/0004013 [BH00] and is described in 4.2.2 of arXiv:1410.7939 [Bor14]. Return a list of indices.
Parameters:  sector –
Sector
; The sector to be decomposed.  indices – iterable of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.
 sector –

pySecDec.decomposition.iterative.
iteration_step
(sector, indices=None)¶ Run a single step of the iterative sector decomposition as described in chapter 3.2 (part II) of arXiv:0803.4177v2 [Hei08]. Return an iterator of
Sector
 the arising subsectors.Parameters:  sector –
Sector
; The sector to be decomposed.  indices – iterable of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.
 sector –

pySecDec.decomposition.iterative.
iterative_decomposition
(sector, indices=None)¶ Run the iterative sector decomposition as described in chapter 3.2 (part II) of arXiv:0803.4177v2 [Hei08]. Return an iterator of
Sector
 the arising subsectors.Parameters:  sector –
Sector
; The sector to be decomposed.  indices – iterable of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.
 sector –

pySecDec.decomposition.iterative.
primary_decomposition
(sector, indices=None)¶ Perform the primary decomposition as described in chapter 3.2 (part I) of arXiv:0803.4177v2 [Hei08]. Return a list of
Sector
 the primary sectors. For N Feynman parameters, there are N primary sectors where the ith Feynman parameter is set to 1 in sector i.See also
Parameters:  sector –
Sector
; The container holding the polynomials (typically and ) to eliminate the Dirac delta from.  indices – iterable of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.
 sector –

pySecDec.decomposition.iterative.
primary_decomposition_polynomial
(polynomial, indices=None)¶ Perform the primary decomposition on a single polynomial.
See also
Parameters:  polynomial –
algebra.Polynomial
; The polynomial to eliminate the Dirac delta from.  indices – iterable of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.
 polynomial –

pySecDec.decomposition.iterative.
remap_parameters
(singular_parameters, Jacobian, *polynomials)¶ Remap the Feynman parameters according to eq. (16) of arXiv:0803.4177v2 [Hei08]. The parameter whose index comes first in singular_parameters is kept fix.
The remapping is done in place; i.e. the polynomials are NOT copied.
Parameters:  singular_parameters – list of integers; The indices such that at least one of polynomials becomes zero if all .
 Jacobian –
Polynomial
; The Jacobian determinant is multiplied to this polynomial.  polynomials – abritrarily many instances of
algebra.Polynomial
where all of these have an equal number of variables; The polynomials of Feynman parameters to be remapped. These are typically and .
Example:
remap_parameters([1,2], Jacobian, F, U)
Geometric¶
The geometric sector decomposition routines.

pySecDec.decomposition.geometric.
Cheng_Wu
(sector, index=1)¶ Replace one Feynman parameter by one. This means integrating out the Dirac delta according to the ChengWu theorem.
Parameters:  sector –
Sector
; The container holding the polynomials (typically and ) to eliminate the Dirac delta from.  index – integer, optional; The index of the Feynman parameter to eliminate. Default: 1 (the last Feynman parameter)
 sector –

pySecDec.decomposition.geometric.
generate_fan
(*polynomials)¶ Calculate the fan of the polynomials in the input. The rays of a cone are given by the exponent vectors after factoring out a monomial together with the standard basis vectors. Each choice of factored out monomials gives a different cone. Only full () dimensional cones in need to be considered.
Parameters: polynomials – abritrarily many instances of Polynomial
where all of these have an equal number of variables; The polynomials to calculate the fan for.

pySecDec.decomposition.geometric.
geometric_decomposition
(sector, indices=None, normaliz='normaliz', workdir='normaliz_tmp')¶ Run the sector decomposition using the geomethod as described in [BHJ+15].
Note
This function calls the command line executable of normaliz [BIR]. See The Geomethod and Normaliz for installation and a list of tested versions.
Parameters:  sector –
Sector
; The sector to be decomposed.  indices – list of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.  normaliz – string; The shell command to run normaliz.
 workdir –
string; The directory for the communication with normaliz. A directory with the specified name will be created in the current working directory. If the specified directory name already exists, an
OSError
is raised.Note
The communication with normaliz is done via files.
 sector –

pySecDec.decomposition.geometric.
geometric_decomposition_ku
(sector, indices=None, normaliz='normaliz', workdir='normaliz_tmp')¶ Run the sector decomposition using the original geometric decomposition strategy by Kaneko and Ueda as described in [KU10].
Note
This function calls the command line executable of normaliz [BIR]. See The Geomethod and Normaliz for installation and a list of tested versions.
Parameters:  sector –
Sector
; The sector to be decomposed.  indices – list of integers or None;
The indices of the parameters to be considered as
integration variables. By default (
indices=None
), all parameters are considered as integration variables.  normaliz – string; The shell command to run normaliz.
 workdir –
string; The directory for the communication with normaliz. A directory with the specified name will be created in the current working directory. If the specified directory name already exists, an
OSError
is raised.Note
The communication with normaliz is done via files.
 sector –

pySecDec.decomposition.geometric.
transform_variables
(polynomial, transformation, polysymbols='y')¶ Transform the parameters of a
pySecDec.algebra.Polynomial
,, where is the transformation matrix.
Parameters:  polynomial –
pySecDec.algebra.Polynomial
; The polynomial to transform the variables in.  transformation – two dimensional array; The transformation matrix .
 polysymbols – string or iterable of strings;
The symbols for the new variables. This argument
is passed to the default constructor of
pySecDec.algebra.Polynomial
. Refer to the documentation ofpySecDec.algebra.Polynomial
for further details.
 polynomial –
Splitting¶
Routines to split the integration between and . This maps singularities from to .

pySecDec.decomposition.splitting.
find_singular_sets_at_one
(polynomial)¶ Find all possible sets of parameters such that the polynomial’s constant term vanishes if these parameters are set to one.
Example:
>>> from pySecDec.algebra import Polynomial >>> from pySecDec.decomposition.splitting import find_singular_sets_at_one >>> polysymbols = ['x0', 'x1'] >>> poly = Polynomial.from_expression('1  10*x0  x1', polysymbols) >>> find_singular_sets_at_one(poly) [(1,)]
Parameters: polynomial – Polynomial
; The polynomial to search in.

pySecDec.decomposition.splitting.
remap_one_to_zero
(polynomial, *indices)¶ Apply the transformation to polynomial for the parameters of the given indices.
Parameters:  polynomial –
Polynomial
; The polynomial to apply the transformation to.  indices – arbitrarily many int;
The indices of the
polynomial.polysymbols
to apply the transformation to.
Example:
>>> from pySecDec.algebra import Polynomial >>> from pySecDec.decomposition.splitting import remap_one_to_zero >>> polysymbols = ['x0'] >>> polynomial = Polynomial.from_expression('x0', polysymbols) >>> remap_one_to_zero(polynomial, 0) + (1) + (1)*x0
 polynomial –

pySecDec.decomposition.splitting.
split
(sector, seed, *indices)¶ Split the integration interval for the parameters given by indices. The splitting point is fixed using numpy’s random number generator.
Return an iterator of
Sector
 the arising subsectors.Parameters: sector – Sector
; The sector to be split. :param seed;
 integer; The seed for the random number generator that is used to fix the splitting point.
Parameters: indices – arbitrarily many integers; The indices of the variables to be split.

pySecDec.decomposition.splitting.
split_singular
(sector, seed, indices=[])¶ Split the integration interval for the parameters that can lead to singularities at one for the polynomials in
sector.cast
.Return an iterator of
Sector
 the arising subsectors.Parameters:  sector –
Sector
; The sector to be split.  seed – integer; The seed for the random number generator that is used to fix the splitting point.
 indices – iterables of integers; The indices of the variables to be split if required. An empty iterator means that all variables may potentially be split.
 sector –
Matrix Sort¶
Algorithms to sort a matrix when column and row permutations are allowed.

pySecDec.matrix_sort.
Pak_sort
(matrix, *indices)¶ Inplace modify the matrix to some canonical ordering, when permutations of rows and columns are allowed.
The indices parameter can contain a list of lists of column indices. Only the columns present in the same list are swapped with each other.
The implementation of this function is described in chapter 2 of [Pak11].
Note
If not all indices are considered the resulting matrix may not be canonical.
See also
Parameters:  matrix – 2D arraylike; The matrix to be canonicalized.
 indices – arbitrarily many iterables of nonnegative integers;
The groups of columns to permute.
Default:
range(1,matrix.shape[1])

pySecDec.matrix_sort.
iterative_sort
(matrix)¶ Inplace modify the matrix to some ordering, when permutations of rows and columns (excluding the first) are allowed.
Note
This function may result in different orderings depending on the initial ordering.
See also
Parameters: matrix – 2D arraylike; The matrix to be canonicalized.

pySecDec.matrix_sort.
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 minimized second column.
Note
This function may result in different orderings depending on the initial ordering.
See also
Parameters: matrix – 2D arraylike; The matrix to be canonicalized.
Subtraction¶
Routines to isolate the divergencies in an expansion.

pySecDec.subtraction.
integrate_by_parts
(polyprod, power_goals, indices)¶ Repeatedly apply integration by parts,
, where denotes the derivative of with respect to . The iteration stops, when power_goal_j.
See also
This function provides an alternative to
integrate_pole_part()
.Parameters:  polyprod –
algebra.Product
of the form<product of <monomial>**(a_j + ...)> * <regulator poles of cal_I> * <cal_I>
; The input product as decribed above. The <product of <monomial>**(a_j + …)> should be apySecDec.algebra.Product
of <monomial>**(a_j + …). as described below. The <monomial>**(a_j + …) should be anpySecDec.algebra.ExponentiatedPolynomial
withexponent
being aPolynomial
of the regulators . Although no dependence on the Feynman parameters is expected in theexponent
, the polynomial variables should be the Feynman parameters and the regulators. The constant term of the exponent should be numerical. The polynomial variables ofmonomial
and the other factors (interpreted as ) are interpreted as the Feynman parameters and the epsilon regulators. Make sure that the last factor (<cal_I>
) is defined and finite for . All poles for should be made explicit by putting them into<regulator poles of cal_I>
aspySecDec.algebra.Pow
withexponent = 1
and thebase
of typepySecDec.algebra.Polynomial
.  power_goals – number or iterable of numbers, e.g. float, integer, …; The stopping criterion for the iteration.
 indices – iterable of integers; The index/indices of the parameter(s) to partially integrate. in the formulae above.
Return the pole part and the numerically integrable remainder as a list. Each returned list element has the same structure as the input polyprod.
 polyprod –

pySecDec.subtraction.
integrate_pole_part
(polyprod, *indices)¶ Transform an integral of the form
into the form
, where denotes the pth derivative of with respect to . The equations above are to be understood schematically.
See also
This function implements the transformation from equation (19) to (21) as described in arXiv:0803.4177v2 [Hei08].
Parameters:  polyprod –
algebra.Product
of the form<product of <monomial>**(a_j + ...)> * <regulator poles of cal_I> * <cal_I>
; The input product as decribed above. The <product of <monomial>**(a_j + …)> should be apySecDec.algebra.Product
of <monomial>**(a_j + …). as described below. The <monomial>**(a_j + …) should be anpySecDec.algebra.ExponentiatedPolynomial
withexponent
being aPolynomial
of the regulators . Although no dependence on the Feynman parameters is expected in theexponent
, the polynomial variables should be the Feynman parameters and the regulators. The constant term of the exponent should be numerical. The polynomial variables ofmonomial
and the other factors (interpreted as ) are interpreted as the Feynman parameters and the epsilon regulators. Make sure that the last factor (<cal_I>
) is defined and finite for . All poles for should be made explicit by putting them into<regulator poles of cal_I>
aspySecDec.algebra.Pow
withexponent = 1
and thebase
of typepySecDec.algebra.Polynomial
.  indices – arbitrarily many integers; The index/indices of the parameter(s) to partially integrate. in the formulae above.
Return the pole part and the numerically integrable remainder as a list. That is the sum and the integrand of equation (21) in arXiv:0803.4177v2 [Hei08]. Each returned list element has the same structure as the input polyprod.
 polyprod –

pySecDec.subtraction.
pole_structure
(monomial_product, *indices)¶ Return a list of the unregulated exponents of the parameters specified by indices in monomial_product.
Parameters:  monomial_product –
pySecDec.algebra.ExponentiatedPolynomial
withexponent
being aPolynomial
; The monomials of the subtraction to extract the pole structure from.  indices – arbitrarily many integers; The index/indices of the parameter(s) to partially investigate.
 monomial_product –
Expansion¶
Routines to series expand singular and nonsingular expressions.

exception
pySecDec.expansion.
OrderError
¶ This exception is raised if an expansion to a lower than the lowest order of an expression is requested.

pySecDec.expansion.
expand_Taylor
(expression, indices, orders)¶ Series/Taylor expand a nonsingular expression around zero.
Return a
algebra.Polynomial
 the series expansion.Parameters:  expression – an expression composed of the types defined in
the module
algebra
; The expression to be series expanded.  indices – integer or iterable of integers; The indices of the parameters to expand. The ordering of the indices defines the ordering of the expansion.
 order – integer or iterable of integers; The order to which the expansion is to be calculated.
 expression – an expression composed of the types defined in
the module

pySecDec.expansion.
expand_singular
(product, indices, orders)¶ Series expand a potentially singular expression of the form
Return a
algebra.Polynomial
 the series expansion.See also
To expand more general expressions use
expand_sympy()
.Parameters:  product –
algebra.Product
with factors of the form<polynomial>
and<polynomial> ** 1
; The expression to be series expanded.  indices – integer or iterable of integers; The indices of the parameters to expand. The ordering of the indices defines the ordering of the expansion.
 order – integer or iterable of integers; The order to which the expansion is to be calculated.
 product –

pySecDec.expansion.
expand_sympy
(expression, variables, orders)¶ Expand a sympy expression in the variables to given orders. Return the expansion as nested
pySecDec.algebra.Polynomial
.See also
This function is a generalization of
expand_singular()
.Parameters:  expression – string or sympy expression; The expression to be expanded
 variables – iterable of strings or sympy symbols; The variables to expand the expression in.
 orders – iterable of integers; The orders to expand to.
Code Writer¶
This module collects routines to create a c++ library.
Make Package¶
This is the main function of pySecDec.

pySecDec.code_writer.
make_package
(name, integration_variables, regulators, requested_orders, polynomials_to_decompose, polynomial_names=[], other_polynomials=[], prefactor=1, remainder_expression=1, functions=[], real_parameters=[], 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_iterative_sort=True, use_light_Pak=True, use_dreadnaut=False, use_Pak=True, processes=None)¶ Decompose, subtract and expand an expression. Return it as c++ package.
See also
In order to decompose a loop integral, use the function
pySecDec.loop_integral.loop_package()
.See also
The generated library is described in Generated C++ Libraries.
Parameters:  name – string; The name of the c++ namepace and the output directory.
 integration_variables – iterable of strings or sympy symbols; The variables that are to be integrated from 0 to 1.
 regulators – iterable of strings or sympy symbols; The UV/IR regulators of the integral.
 requested_orders – iterable of integers; Compute the expansion in the regulators to these orders.
 polynomials_to_decompose – iterable of strings or sympy expressions or
pySecDec.algebra.ExponentiatedPolynomial
orpySecDec.algebra.Polynomial
; The polynomials to be decomposed.  polynomial_names – iterable of strings; Assign symbols for the polynomials_to_decompose. These can be referenced in the other_polynomials; see other_polynomials for details.
 other_polynomials –
iterable of strings or sympy expressions or
pySecDec.algebra.ExponentiatedPolynomial
orpySecDec.algebra.Polynomial
; Additional polynomials where no decomposition is attempted. The symbols defined in polynomial_names can be used to reference the polynomials_to_decompose. This is particularly useful when computing loop integrals where the “numerator” can depend on the first and second Symanzik polynomials.Example (1loop bubble with numerator):
>>> polynomials_to_decompose = ["(x0 + x1)**(2*eps  4)", ... "(p**2*x0*x1)**(eps))"] >>> polynomial_names = ["U", "F"] >>> other_polynomials = [""" (eps  1)*s*U**2 ... + (eps  2)*F ...  (eps  1)*2*s*x0*U ... + (eps  1)*s*x0**2"""]
See also
Note that the polynomial_names refer to the polynomials_to_decompose without their exponents.
 prefactor – string or sympy expression, optional; A factor that does not depend on the integration variables.
 remainder_expression –
string or sympy expression or
pySecDec.algebra._Expression
, optional; An additional factor.Dummy function must be provided with all arguments, e.g.
remainder_expression='exp(eps)*f(x0,x1)'
. In addition, all dummy function must be listed in functions.  functions –
iterable of strings or sympy symbols, optional; Function symbols occuring in remainder_expression, e.g.``[‘f’]``.
Note
Only userdefined functions that are provided as c++callable code should be mentioned here. Listing basic mathematical functions (e.g.
log
,pow
,exp
,sqrt
, …) is not required and considered an error to avoid name conflicts.Note
The power function pow and the logarithm log use the nonstandard continuation with an infinitesimal negative imaginary part on the negative real axis (e.g.
log(1) = i*pi
).  real_parameters – iterable of strings or sympy symbols, optional; Symbols to be interpreted as real variables.
 complex_parameters – iterable of strings or sympy symbols, optional; Symbols to be interpreted as complex variables.
 form_optimization_level – integer out of the interval [0,4], optional;
The optimization level to be used in FORM.
Default:
2
.  form_work_space – string, optional;
The FORM WorkSpace. Default:
'500M'
.  form_insertion_depth – nonnegative integer, optional;
How deep FORM should try to resolve nested function
calls. Default:
5
.  contour_deformation_polynomial – string or sympy symbol, optional;
The name of the polynomial in polynomial_names
that is to be continued to the complex plane
according to a prescription.
For loop integrals, this is the second Symanzik
polynomial
F
. If not provided, no code for contour deformation is created.  positive_polynomials – iterable of strings or sympy symbols, optional;
The names of the polynomials in polynomial_names
that should always have a positive real part.
For loop integrals, this applies to the first Symanzik
polynomial
U
. If not provided, no polynomial is checked for positiveness. If contour_deformation_polynomial isNone
, this parameter is ignored.  decomposition_method –
string, optional; The strategy to decompose the polynomials. The following strategies are available:
 ’iterative_no_primary’ (default)
 ’geometric_no_primary’
 ’iterative’
 ’geometric’
 ’geometric_ku’
’iterative’, ‘geometric’, and ‘geometric_ku’ are only valid for loop integrals. An end user should always use ‘iterative_no_primary’ or ‘geometric_no_primary’ here. In order to compute loop integrals, please use the function
pySecDec.loop_integral.loop_package()
.  normaliz_executable – string, optional; The command to run normaliz. normaliz is only required if decomposition_method starts with ‘geometric’. Default: ‘normaliz’
 enforce_complex – bool, optional;
Whether or not the generated integrand functions
should have a complex return type even though
they might be purely real.
The return type of the integrands is automatically
complex if contour_deformation is
True
or if there are complex_parameters. In other cases, the calculation can typically be kept purely real. Most commonly, this flag is needed iflog(<negative real>)
occurs in one of the integrand functions. However, pySecDec will suggest setting this flag toTrue
in that case. Default:False
 split – bool or integer, optional;
Whether or not to split the integration domain
in order to map singularities from to
. Set this option to
True
if you have singularties when one or more integration variables are one. If an integer is passed, that integer is used as seed to generate the splitting point. Default:False
 ibp_power_goal –
number or iterable of number, optional; The power_goal that is forwarded to
integrate_by_parts()
.This option controls how the subtraction terms are generated. Setting it to
numpy.inf
disablesintegrate_by_parts()
, while0
disablesintegrate_pole_part()
.See also
To generate the subtraction terms, this function first calls
integrate_by_parts()
for each integration variable with the give ibp_power_goal. Thenintegrate_pole_part()
is called.Default:
1
 use_iterative_sort – bool;
Whether or not to use
squash_symmetry_redundant_sectors_sort()
withiterative_sort()
to find sector symmetries. Default:True
 use_light_Pak – bool;
Whether or not to use
squash_symmetry_redundant_sectors_sort()
withlight_Pak_sort()
to find sector symmetries. Default:True
 use_dreadnaut – bool or string, optional;
Whether or not to use
squash_symmetry_redundant_sectors_dreadnaut()
to find sector symmetries. If given a string, interpret that string as the command line executable dreadnaut. IfTrue
, try$SECDEC_CONTRIB/bin/dreadnaut
and, if the environment variable$SECDEC_CONTRIB
is not set,dreadnaut
. Default:False
 use_Pak – bool;
Whether or not to use
squash_symmetry_redundant_sectors_sort()
withPak_sort()
to find sector symmetries. Default:True
 processes – integer or None, optional;
The maximal number of processes to be used. If
None
, the number of CPUsmultiprocessing.cpu_count()
is used. New in version 1.3. Default:None
Template Parser¶
Functions to generate c++ sources from template files.

pySecDec.code_writer.template_parser.
parse_template_file
(src, dest, replacements={})¶ Copy a file from src to dest replacing
%(...)
instructions in the standard python way.Warning
If the file specified in dest exists, it is overwritten without prompt.
See also
Parameters:  src – str; The path to the template file.
 dest – str; The path to the destination file.
 replacements –
dict; The replacements to be performed. The standard python replacement rules apply:
>>> '%(var)s = %(value)i' % dict( ... var = 'my_variable', ... value = 5) 'my_variable = 5'

pySecDec.code_writer.template_parser.
parse_template_tree
(src, dest, replacements_in_files={}, filesystem_replacements={})¶ Copy a directory tree from src to dest using
parse_template_file()
for each file and replacing the filenames according to filesystem_replacements.See also
Parameters:  src – str; The path to the template directory.
 dest – str; The path to the destination directory.
 replacements_in_files –
dict; The replacements to be performed in the files. The standard python replacement rules apply:
>>> '%(var)s = %(value)i' % dict( ... var = 'my_variable', ... value = 5) 'my_variable = 5'
 filesystem_replacements – dict;
Renaming rules for the destination files. and
directories. If a file or directory name in
the source tree src matches a key in this
dictionary, it is renamed to the corresponding
value. If the value is
None
, the corresponding file is ignored.
Generated C++ Libraries¶
A C++ Library to numerically compute a given integral (loop integral) can be generated by the make_package()
(loop_package()
) functions.
The name passed to the make_package()
or loop_package()
function will be used as the C++ namespace of the generated library.
A program demonstrating the use of the C++ library is generated for each integral and written to name/integrate_name.cpp
. Here we document the C++ library API.
See also

typedef double
real_t
¶ The real type used by the library.

type
integrand_return_t
¶ The return type of the integrand function. If the integral has complex parameters or uses contour deformation or if enforce_complex is set to
True
in the call tomake_package()
orloop_package()
then integrand_return_t is complex_t. Otherwise integrand_return_t is real_t.

template<typename
T
>
usingnested_series_t
= secdecutil::Series<secdecutil::Series<...<T>>>¶ A potentially nested
secdecutil::Series
representing the series expansion in each of the regulators. If the integral depends on only one regulator (for example, a loop integral generated withloop_package()
) this type will be asecdecutil::Series
. For integrals that depend on multiple regulators then this will be a series of series representing the multivariate series. This type can be used to write code that can handle integrals depending on arbitrarily many regulators.See also

typedef secdecutil::IntegrandContainer<integrand_return_t, real_t const *const >
integrand_t
¶ The type of the integrand. Within the generated C++ library integrands are stored in a container along with the number of integration variables upon which they depend. These containers can be passed to an integrator for numerical integration.
See also

type
cuda_integrand_t
¶ New in version 1.4.
The type of a single integrand (sector) usable on a CUDA device (GPU). This container can be passed to an integrator for numerical integration.

type
cuda_together_integrand_t
¶ New in version 1.4.
The type of a sum of integrands (sectors) usable on a CUDA device (GPU). This container can be passed to an integrator for numerical integration.

const unsigned long long
number_of_sectors
¶ The number of sectors generated by the sector decomposition.
Changed in version 1.3.1: Type was
unsigned int
in earlier versions of pySecDec.

const unsigned int
maximal_number_of_integration_variables
¶ The number of integration variables after primary decomposition. This provides an upper bound in the number of integration variables for all integrand functions. The actual number of integration variables may be lower for a given integrand.

const unsigned int
number_of_regulators
¶ The number of regulators on which the integral depends.

const unsigned int
number_of_real_parameters
¶ The number of real parameters on which the integral depends.

const std::vector<std::string>
names_of_real_parameters
¶ An ordered vector of string representations of the names of the real parameters.

const unsigned int
number_of_complex_parameters
¶ The number of complex parameters on which the integral depends.

const std::vector<std::string>
names_of_complex_parameters
¶ An ordered vector of string representations of the names of the complex parameters.

const std::vector<int>
lowest_orders
¶ A vector of the lowest order of each regulator which appears in the integral, not including the prefactor.

const std::vector<int>
highest_orders
¶ A vector of the highest order of each regulator which appears in the integral, not including the prefactor. This depends on the requested_orders and prefactor/additional_prefactor parameter passed to
make_package()
orloop_package()
. In the case ofloop_package()
it also depends on the function prefactor of the integral which appears upon Feynman parametrization.

const std::vector<int>
lowest_prefactor_orders
¶ A vector of the lowest order of each regulator which appears in the prefactor of the integral.

const std::vector<int>
highest_prefactor_orders
¶ A vector of the highest order of each regulator which appears in the prefactor of the integral.

const std::vector<int>
requested_orders
¶ A vector of the requested orders of each regulator used to generate the C++ library, i.e. the requested_orders parameter passed to
make_package()
orloop_package()
.

const std::vector<nested_series_t<sector_container_t>> &
get_sectors
()¶ Changed in version 1.3.1: The variable
sectors
has been replaced by this function.A low level interface for obtaining the underlying integrand C++ functions.
Warning
The precise definition and usage of
get_sectors()
is likely to change in future versions of pySecDec.

nested_series_t<integrand_return_t>
prefactor
(const std::vector<real_t> &real_parameters, const std::vector<complex_t> &complex_parameters)¶ The series expansion of the integral prefactor evaluated with the given parameters. If the library was generated using
make_package()
it will be equal to the prefactor passed tomake_package()
. If the library was generated withloop_package()
it will be the product of the additional_prefactor passed toloop_package()
and the function prefactor of the integral which appears upon Feynman parametrization.

const std::vector<std::vector<real_t>>
pole_structures
¶ A vector of the powers of the monomials that can be factored out of each sector of the polynomial during the decomposition.
Example: an integral depending on variables and may have two sectors, the first may have a monomial factored out and the second may have a monomial factored out during the decomposition. The resulting pole_structures would read
{ {1,2}, {1,0} }
. Poles of type are known as logarithmic poles, poles of type are known as linear poles.

std::vector<nested_series_t<integrand_t>>
make_integrands
(const std::vector<real_t> &real_parameters, const std::vector<complex_t> &complex_parameters)¶ (without contour deformation)

std::vector<nested_series_t<cuda_integrand_t>>
make_cuda_integrands
(const std::vector<real_t> &real_parameters, const std::vector<complex_t> &complex_parameters)¶ New in version 1.4.
(without contour deformation) (CUDA only)

std::vector<nested_series_t<integrand_t>>
make_integrands
(const std::vector<real_t> &real_parameters, const std::vector<complex_t> &complex_parameters, unsigned number_of_presamples = 100000, real_t deformation_parameters_maximum = 1., real_t deformation_parameters_minimum = 1.e5, real_t deformation_parameters_decrease_factor = 0.9)¶ (with contour deformation)

std::vector<nested_series_t<cuda_integrand_t>>
make_cuda_integrands
(const std::vector<real_t> &real_parameters, const std::vector<complex_t> &complex_parameters, unsigned number_of_presamples = 100000, real_t deformation_parameters_maximum = 1., real_t deformation_parameters_minimum = 1.e5, real_t deformation_parameters_decrease_factor = 0.9)¶ New in version 1.4.
(with contour deformation) (CUDA only)
Gives a vector containing the series expansions of individual sectors of the integrand after sector decomposition with the specified real_paraemters and complex_parameters bound. Each element of the vector contains the series expansion of an individual sector. The series consists of instances of
integrand_t
(cuda_integrand_t
) which contain the integrand functions and the number of integration variables upon which they depend. The real and complex parameters are bound to the values passed in real_parameters and complex_parameters. If enabled, contour deformation is controlled by the parameters number_of_presamples, deformation_parameters_maximum, deformation_parameters_minimum, deformation_parameters_decrease_factor which are documented inpySecDec.integral_interface.IntegralLibrary
. In case of a sign check error (sign_check_error), manually set number_of_presamples, deformation_parameters_maximum, and deformation_parameters_minimum.Passing the integrand_t to the
secdecutil::Integrator::integrate()
function of an instance of a particularsecdecutil::Integrator
will return the numerically evaluated integral. To integrate all orders of all sectorssecdecutil::deep_apply()
can be used.Note
This is the recommended way to access the integrand functions.
Integral Interface¶
An interface to libraries generated by
pySecDec.code_writer.make_package()
or
pySecDec.loop_integral.loop_package()
.

class
pySecDec.integral_interface.
CPPIntegrator
¶ Abstract base class for integrators to be used with an
IntegralLibrary
. This class holds a pointer to the c++ integrator and defines the destructor.

class
pySecDec.integral_interface.
CQuad
(integral_library, epsrel=0.01, epsabs=1e07, n=100, verbose=False, zero_border=0.0)¶ Wrapper for the cquad integrator defined in the gsl library.
Parameters: integral_library – IntegralLibrary
; The integral to be computed with this integrator.The other options are defined in Section 4.5.1 and in the gsl manual.

class
pySecDec.integral_interface.
CudaQmc
(integral_library, transform, fitfunction='default', generatingvectors='default', epsrel=0.0, epsabs=0.0, maxeval=0, errormode='default', evaluateminn=0, minn=0, minm=0, maxnperpackage=0, maxmperpackage=0, cputhreads=0, cudablocks=0, cudathreadsperblock=0, verbosity=0, seed=0, devices=[])¶ Wrapper for the Qmc integrator defined in the integrators library for GPU use.
Parameters:  integral_library –
IntegralLibrary
; The integral to be computed with this integrator.  errormode – string;
The errormode parameter of the Qmc, can be
"default"
,"all"
, and"largest"
."default"
takes the default from the Qmc library. See the Qmc docs for details on the other settings.  transform – string;
An integral transform related to the parameter
P of the Qmc. The possible choices
correspond to the integral transforms of the
underlying Qmc implementation. Possible values
are,
"none"
,"baker"
,sidi#
,"korobov#"
, andkorobov#x#
where any#
(the rank of the Korobov/Sidi transform) must be an integer between 1 and 6.  fitfunction – string;
An integral transform related to the parameter
F of the Qmc. The possible choices
correspond to the integral transforms of the
underlying Qmc implementation. Possible values
are
"default"
,"none"
,"polysingular"
.  generatingvectors – string;
The name of a set of generating vectors.
The possible choices correspond to the available generating
vectors of the underlying Qmc implementation. Possible values
are
"default"
,"cbcpt_dn1_100"
,"cbcpt_dn2_6"
and"cbcpt_cfftw1_6"
.
The other options are defined in the Qmc docs. If an argument is set to 0 then the default of the underlying Qmc implementation is used.
 integral_library –

class
pySecDec.integral_interface.
Cuhre
(integral_library, epsrel=0.01, epsabs=1e07, flags=0, mineval=0, maxeval=1000000, zero_border=0.0, key=0, real_complex_together=False)¶ Wrapper for the Cuhre integrator defined in the cuba library.
Parameters: integral_library – IntegralLibrary
; The integral to be computed with this integrator.The other options are defined in Section 4.5.3 and in the cuba manual.

class
pySecDec.integral_interface.
Divonne
(integral_library, epsrel=0.01, epsabs=1e07, flags=0, seed=0, mineval=0, maxeval=1000000, zero_border=0.0, key1=2000, key2=1, key3=1, maxpass=4, border=0.0, maxchisq=1.0, mindeviation=0.15, real_complex_together=False)¶ Wrapper for the Divonne integrator defined in the cuba library.
Parameters: integral_library – IntegralLibrary
; The integral to be computed with this integrator.The other options are defined in Section 4.5.3 and in the cuba manual.

class
pySecDec.integral_interface.
IntegralLibrary
(shared_object_path)¶ Interface to a c++ library produced by
make_package()
orloop_package()
.Parameters: shared_object_path – str; The path to the file “<name>_pylink.so” that can be built by the command
$ make pylink
in the root directory of the c++ library.
Instances of this class can be called with the following arguments:
Parameters:  real_parameters – iterable of float; The real_parameters of the library.
 complex_parameters – iterable of complex; The complex parameters of the library.
 together – bool, optional;
Whether to integrate the sum of all sectors
or to integrate the sectors separately.
Default:
True
.  number_of_presamples – unsigned int, optional;
The number of samples used for the
contour optimization.
A larger value here may resolve a sign
check error (sign_check_error).
This option is ignored if the integral
library was created without deformation.
Default:
100000
.  deformation_parameters_maximum – float, optional;
The maximal value the deformation parameters
can obtain.
Lower this value if you get a sign check
error (sign_check_error).
If
number_of_presamples=0
, all are set to this value. This option is ignored if the integral library was created without deformation. Default:1.0
.  deformation_parameters_minimum – float, optional;
The minimal value the deformation parameters
can obtain.
Lower this value if you get a sign check
error (sign_check_error).
If
number_of_presamples=0
, all are set to this value. This option is ignored if the integral library was created without deformation. Default:1e5
.  deformation_parameters_decrease_factor – float, optional;
If the sign check with the optimized
fails during the presampling
stage, all are multiplied
by this value until the sign check passes.
We recommend to rather change
number_of_presamples
,deformation_parameters_maximum
, anddeformation_parameters_minimum
in case of a sign check error. This option is ignored if the integral library was created without deformation. Default:0.9
.
The call operator returns three strings: * The integral without its prefactor * The prefactor * The integral multiplied by the prefactor
The integrator can be configured by calling the member methods
use_Vegas()
,use_Suave()
,use_Divonne()
,use_Cuhre()
,use_CQuad()
, anduse_Qmc()
. The available options are listed in the documentation ofVegas
,Suave
,Divonne
,Cuhre
,CQuad
,Qmc
(CudaQmc
for GPU version), respectively.CQuad
can only be used for one dimensional integrals. A call touse_CQuad()
configures the integrator to useCQuad
if possible (1D) and the previously defined integrator otherwise. By default,CQuad
(1D only) andVegas
are used with their default arguments. For details about the options, refer to the cuba and the gsl manual.Further information about the library is stored in the member variable info of type
dict
.

class
pySecDec.integral_interface.
MultiIntegrator
(integral_library, low_dim_integrator, high_dim_integrator, critical_dim)¶ New in version 1.3.1.
Wrapper for the
secdecutil::MultiIntegrator
.Parameters:  integral_library –
IntegralLibrary
; The integral to be computed with this integrator.  low_dim_integrator –
CPPIntegrator
; The integrator to be used if the integrand is lower dimensional than critical_dim.  high_dim_integrator –
CPPIntegrator
; The integrator to be used if the integrand has dimension critical_dim or higher.  critical_dim – integer; The dimension below which the low_dimensional_integrator is used.
Use this class to switch between integrators based on the dimension of the integrand when integrating the integral_ibrary. For example, “
CQuad
for 1D andVegas
otherwise” is implemented as:integral_library.integrator = MultiIntegrator(integral_library,CQuad(integral_library),Vegas(integral_library),2)
MultiIntegrator
can be nested to implement multiple critical dimensions. To use e.g.CQuad
for 1D,Cuhre
for 2D and 3D, andVegas
otherwise, do:integral_library.integrator = MultiIntegrator(integral_library,CQuad(integral_library),MultiIntegrator(integral_library,Cuhre(integral_library),Vegas(integral_library),4),2)
Warning
The integral_library passed to the integrators must be the same for all of them. Furthermore, an integrator can only be used to integrate the integral_library it has beeen constructed with.
Warning
The
MultiIntegrator
cannot be used withCudaQmc
. integral_library –

class
pySecDec.integral_interface.
Qmc
(integral_library, transform, fitfunction='default', generatingvectors='default', epsrel=0.0, epsabs=0.0, maxeval=0, errormode='default', evaluateminn=0, minn=0, minm=0, maxnperpackage=0, maxmperpackage=0, cputhreads=0, cudablocks=0, cudathreadsperblock=0, verbosity=0, seed=0, devices=[])¶ Wrapper for the Qmc integrator defined in the integrators library.
Parameters:  integral_library –
IntegralLibrary
; The integral to be computed with this integrator.  errormode – string;
The errormode parameter of the Qmc, can be
"default"
,"all"
, and"largest"
."default"
takes the default from the Qmc library. See the Qmc docs for details on the other settings.  transform – string;
An integral transform related to the parameter
P of the Qmc. The possible choices
correspond to the integral transforms of the
underlying Qmc implementation. Possible values
are,
"none"
,"baker"
,sidi#
,"korobov#"
, andkorobov#x#
where any#
(the rank of the Korobov/Sidi transform) must be an integer between 1 and 6.  fitfunction – string;
An integral transform related to the parameter
F of the Qmc. The possible choices
correspond to the integral transforms of the
underlying Qmc implementation. Possible values
are
"default"
,"none"
,"polysingular"
.  generatingvectors – string;
The name of a set of generating vectors.
The possible choices correspond to the available generating
vectors of the underlying Qmc implementation. Possible values
are
"default"
,"cbcpt_dn1_100"
,"cbcpt_dn2_6"
and"cbcpt_cfftw1_6"
.
See also
The most important options are described in Section 4.5.2.
The other options are defined in the Qmc docs. If an argument is set to 0 then the default of the underlying Qmc implementation is used.
 integral_library –

class
pySecDec.integral_interface.
Suave
(integral_library, epsrel=0.01, epsabs=1e07, flags=0, seed=0, mineval=0, maxeval=1000000, zero_border=0.0, nnew=1000, nmin=10, flatness=25.0, real_complex_together=False)¶ Wrapper for the Suave integrator defined in the cuba library.
Parameters: integral_library – IntegralLibrary
; The integral to be computed with this integrator.The other options are defined in Section 4.5.3 and in the cuba manual.

class
pySecDec.integral_interface.
Vegas
(integral_library, epsrel=0.01, epsabs=1e07, flags=0, seed=0, mineval=0, maxeval=1000000, zero_border=0.0, nstart=1000, nincrease=500, nbatch=1000, real_complex_together=False)¶ Wrapper for the Vegas integrator defined in the cuba library.
Parameters: integral_library – IntegralLibrary
; The integral to be computed with this integrator.The other options are defined in Section 4.5.3 and in the cuba manual.
Miscellaneous¶
Collection of generalpurpose helper functions.

pySecDec.misc.
adjugate
(M)¶ Calculate the adjugate of a matrix.
Parameters: M – a squarematrixlike array;

pySecDec.misc.
all_pairs
(iterable)¶ Return all possible pairs of a given set.
all_pairs([1,2,3,4]) > [(1,2),(3,4)] [(1,3),(2,4)] [(1,4),(2,3)]
Parameters: iterable – iterable; The set to be split into all possible pairs.

pySecDec.misc.
argsort_2D_array
(array)¶ Sort a 2D array according to its row entries. The idea is to bring identical rows together.
See also
If your array is not two dimesional use
argsort_ND_array()
. Example:
input sorted 1 2 3 1 2 3 2 3 4 1 2 3 1 2 3 2 3 4
Return the indices like numpy’s
argsort()
would.Parameters: array – 2D array; The array to be argsorted.

pySecDec.misc.
argsort_ND_array
(array)¶ Like
argsort_2D_array()
, this function groups identical entries in an array with any dimensionality greater than (or equal to) two together.Return the indices like numpy’s
argsort()
would.See also
Parameters: array – ND array, ; The array to be argsorted.

pySecDec.misc.
assert_degree_at_most_max_degree
(expression, variables, max_degree, error_message)¶ Assert that expression is a polynomial of degree less or equal max_degree in the variables.

pySecDec.misc.
cached_property
(method)¶ Like the builtin property to be used as decorator but the method is only called once per instance.
Example:
class C(object): 'Sum up the numbers from one to `N`.' def __init__(self, N): self.N = N @cached_property def sum(self): result = 0 for i in range(1, self.N + 1): result += i return result

pySecDec.misc.
det
(M)¶ Calculate the determinant of a matrix.
Parameters: M – a squarematrixlike array;

pySecDec.misc.
doc
(docstring)¶ Decorator that replaces a function’s docstring with docstring.
Example:
@doc('documentation of `some_funcion`') def some_function(*args, **kwargs): pass

pySecDec.misc.
flatten
(polynomial, depth=inf)¶ Convert nested polynomials; i.e. polynomials that have polynomials in their coefficients to one single polynomial.
Parameters:  polynomial –
pySecDec.algebra.Polynomial
; The polynomial to “flatten”.  depth – integer;
The maximum number of recursion steps.
If not provided, stop if the coefficient
is not a
pySecDec.algebra.Polynomial
.
 polynomial –

pySecDec.misc.
lowest_order
(expression, variable)¶ Find the lowest order of expression’s series expansion in variable.
Example:
>>> from pySecDec.misc import lowest_order >>> lowest_order('exp(eps)', 'eps') 0 >>> lowest_order('gamma(eps)', 'eps') 1
Parameters:  expression – string or sympy expression; The expression to compute the lowest expansion order of.
 variable – string or sympy expression; The variable in which to expand.

pySecDec.misc.
missing
(full, part)¶ Return the elements in full that are not contained in part. Raise ValueError if an element is in part but not in full.
missing([1,2,3], [1]) > [2,3]
missing([1,2,3,1], [1,2]) > [3,1]
missing([1,2,3], [1,'a']) > ValueError
Parameters:  full – iterable; The set of elements to complete part with.
 part – iterable; The set to be completed to a superset of full.

pySecDec.misc.
parallel_det
(M, pool)¶ Calculate the determinant of a matrix in parallel.
Parameters:  M – a squarematrixlike array;
 pool –
multiprocessing.Pool
; The pool to be used.
Example:
>>> from pySecDec.misc import parallel_det >>> from multiprocessing import Pool >>> from sympy import sympify >>> M = [['m11','m12','m13','m14'], ... ['m21','m22','m23','m24'], ... ['m31','m32','m33','m34'], ... ['m41','m42','m43','m44']] >>> M = sympify(M) >>> parallel_det(M, Pool(2)) # 2 processes m11*(m22*(m33*m44  m34*m43)  m23*(m32*m44  m34*m42) + m24*(m32*m43  m33*m42))  m12*(m21*(m33*m44  m34*m43)  m23*(m31*m44  m34*m41) + m24*(m31*m43  m33*m41)) + m13*(m21*(m32*m44  m34*m42)  m22*(m31*m44  m34*m41) + m24*(m31*m42  m32*m41))  m14*(m21*(m32*m43  m33*m42)  m22*(m31*m43  m33*m41) + m23*(m31*m42  m32*m41))

pySecDec.misc.
powerset
(iterable, min_length=0, stride=1)¶ Return an iterator over the powerset of a given set.
powerset([1,2,3]) > () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)
Parameters:  iterable – iterable; The set to generate the powerset for.
 min_length – integer, optional;
Only generate sets with minimal given length.
Default:
0
.  stride – integer;
Only generate sets that have a multiple of
stride elements.
powerset([1,2,3], stride=2) > () (1,2) (1,3) (2,3)

pySecDec.misc.
rangecomb
(low, high)¶ Return an iterator over the occuring orders in a multivariate series expansion between low and high.
Parameters:  low – vectorlike array; The lowest orders.
 high – vectorlike array; The highest orders.
Example:
>>> from pySecDec.misc import rangecomb >>> all_orders = rangecomb([1,2], [0,0]) >>> list(all_orders) [(1, 2), (1, 1), (1, 0), (0, 2), (0, 1), (0, 0)]

pySecDec.misc.
sympify_symbols
(iterable, error_message, allow_number=False)¶ sympify each item in iterable and assert that it is a symbol.
Frequently Asked Questions¶
How can I adjust the integrator parameters?¶
If the python interface is used for the numerical integration, i.e. a python script like examples/integrate_box1L.py
, the integration parameters can be specified in the argument list of the integrator call.
For example, using Vegas as integrator:
box1L.use_Vegas(flags=2, epsrel=1e3, epsabs=1e12, nstart=5000, nincrease=10000, maxeval=10000000, real_complex_together=True)
Or, using Divonne as integrator:
box1L.use_Divonne(flags=2, epsrel=1e3, epsabs=1e12, maxeval=10000000, border=1e8, real complex together=True)
The parameter real complex together tells the integrator to integrate real and imaginary parts simultaneously. A complete list of possible options for the integrators can be found in integral_interface
.
If the C++ interface is used, the options can be specified as fields of the integrator.
For example, after running examples/generate_box1L.py
, in the file examples/box1L/integrate_box1L.cpp
, you can modify the corresponding block to e.g.:
// Integrate
secdecutil::cuba::Vegas<box1L::integrand_return_t> integrator;
integrator.flags = 2; // verbose output > see cuba manual
integrator.epsrel = 1e2;
integrator.epsabs = 1e12;
integrator.nstart = 5000;
integrator.nincrease = 10000;
integrator.maxeval = 10000000;
integrator.together = true;
In order to set the Divonne integrator with the same parameters as above, do:
// Integrate
secdecutil::cuba::Divonne<box1L::integrand_return_t> integrator;
integrator.flags = 2; // verbose output > see cuba manual
integrator.epsrel = 1e2;
integrator.epsabs = 1e12;
integrator.maxeval = 10000000;
integrator.border = 1e8;
integrator.together = true;
More information about the C++ integrator class can be found in Section 4.5.
How can I request a higher numerical accuracy?¶
The integrator stops if any of the following conditions is fulfilled: (1) epsrel
is reached, (2) epsabs
is reached, (3) maxeval
is reached.
Therefore, setting these parameters accordingly will cause the integrator to make more iterations and reach a more accurate result.
How can I tune the contour deformation parameters?¶
You can specify the parameters in the argument of the integral call in the python script for the integration, see e.g. line 12 of examples/integrate box1L.py
:
str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor=box1L(real_parameters=[4.,0.75,1.25,1.],number_of_presamples=10**6,deformation_parameters_maximum=0.5)
This sets the number of presampling points to 10**6
(default: 10**5
) and the maximum value for the contour deformation parameter deformation_parameters_maximum
to 0.5
(default: 1
). The user should make sure that deformation parameters maximum is always larger than deformation_parameters_minimum (default: 1e5
). These parameters are described in IntegralLibrary
.
What can I do if the program stops with an error message containing sign_check_error?¶
This error occurs if the contour deformation leads to a wrong sign of the Feynman prescription, usually due to the fact that the deformation parameter is too large. Choose a larger value for number_of_presamples
and a smaller value (e.g. 0.5
) for deformation_parameters_maximum
(see item above). If that does not help, you can try 0.1
instead of 0.5
for deformation_parameters_maximum
. The relevant parameters are described in IntegralLibrary
.
What does additional_prefactor mean exactly?¶
We should first point out that the conventions for additional prefactors defined by the user have been changed between SecDec 3 and pySecDec. The prefactor specified by the user will now be included in the numerical result.
To make clear what is meant by “additional”, we repeat our conventions for Feynman integrals here.
A scalar Feynman graph in dimensions at loops with propagators, where the propagators can have arbitrary, not necessarily integer powers , has the following representation in momentum space:
where the are linear combinations of external momenta and loop momenta .
Introducing Feynman parameters leads to:
The prefactor coming from the Feynman parametrisation will always be included in the numerical result, corresponding to additional_prefactor=1 (default), i.e. the program will return the numerical value for . If the user defines additional_prefactor=’gamma(32*eps)’, this prefactor will be expanded in and included in the numerical result returned by pySecDec, in addition to the one coming from the Feynman parametrisation.
For general polynomials not related to loop integrals, i.e. in make_package
, the prefactor provided by the user is the only prefactor, as there is no prefactor coming from a Feynman parametrisation in this case. This is the reason why in make_package
the keyword for the prefactor defined by the user is prefactor
, while in loop_package
it is additional_prefactor
.
What can I do if I get nan?¶
This means that the integral does not converge which can have several reasons. When Divonne is used as an integrator, it is important to use a nonzero value for border, e.g. border=1e8
. Vegas is in general the most robust integrator. When using Vegas, try to increase the values for nstart
and nincrease
, for example nstart=10000
(default: 1000
) and nincrease=5000
(default: 500
).
If the integral is nonEuclidean, make sure that contour_deformation=True is set.
Another reason for getting nan can be that the integral has singularities at and therefore needs usage of the split
option, see item below.
What can I use as numerator of a loop integral?¶
The numerator must be a sum of products of numbers, scalar products (e.g. p1(mu)*k1(mu)*p1(nu)*k2(nu)
and/or symbols (e.g. m
). The numerator can also be an inverse propagator.
In addition, the numerator must be finite in the limit . The default numerator is 1
.
Examples:
p1(mu)*k1(mu)*p1(nu)*k2(nu) + 4*s*eps*k1(mu)*k1(mu)
p1(mu)*(k1(mu) + k2(mu))*p1(nu)*k2(nu)
p1(mu)*k1(mu)
More details can be found in LoopIntegralFromPropagators
.
How can I integrate just one coefficient of a particular order in the regulator?¶
You can pick a certain order in the C++ interface (see C++ Interface (advanced)). To integrate only one order, for example the finite part, change the line:
const box1L::nested_series_t<secdecutil::UncorrelatedDeviation<box1L::integrand_return_t>> result_all = secdecutil::deep_apply( all_sectors, integrator.integrate );
to:
int order = 0; // compute finite part only
const secdecutil::UncorrelatedDeviation<box1L::integrand_return_t> result_order = secdecutil::deep_apply(all_sectors.at(order), integrator.integrate );
where box1L
is to be replaced by the name of your integral. In addition, you should change the lines:
std::cout << " integral without prefactor  " << std::endl;
std::cout << result_all << std::endl << std::endl;
to:
std::cout << " integral without prefactor  " << std::endl;
std::cout << result_order << std::endl << std::endl;
and remove the lines:
std::cout << " prefactor  " << std::endl;
const box1L::nested_series_t<box1L::integrand_return_t> prefactor = box1L::prefactor(real_parameters, complex_parameters);
std::cout << prefactor << std::endl << std::endl;
std::cout << " full result (prefactor*integral)  " << std::endl;
std::cout << prefactor*result_all << std::endl;
because the expansion of the prefactor will in general mix with the pole coefficients and thus affect the finite part. We should point out however that deleting these lines also means that the result will not contain any prefactor, not even the one coming from the Feynman parametrisation.
How can I use complex masses?¶
In the python script generating the expressions for the integral, define mass symbols in the same way as for real masses, e.g:
Mandelstam_symbols=['s']
mass_symbols=['msq']
Then, in loop_package
define:
real parameters = Mandelstam_symbols,
complex parameters = mass_symbols,
In the integration script (using the python interface), the numerical values for the complex parameters are given after the ones for the real parameters:
str_integral_without_prefactor, str_prefactor, str_integral_with_prefactor = integral(real_parameters=[4.],complex_parameters=[1.0.0038j])
Note that in python the letter j
is used rather than i
for the imaginary part.
In the C++ interface, you can set (for the example triangle2L):
const std::vector<triangle2L::real_t> real_parameters = { 4. };
const std::vector<triangle2L::complex_t> complex_parameters = { {1.,0.0038} };
When should I use the “split” option?¶
The modules loop_package
and make_package
have the option to split the integration domain (split=True
). This option can be useful for integrals which do not have a Euclidean region. If certain kinematic conditions are fulfilled, for example if the integral contains massive onshell lines, it can happen that singularities at remain in the polynomial after the decomposition. The split option remaps these singularities to the origin of parameter space. If your integral is of this type, and with the standard approach the numerical integration does not seem to converge, try the split
option. It produces a lot more sectors, so it should not be used without need. We also would like to mention that very often a change of basis to increase the (negative) power of the polynomial can be beneficial if integrals of this type occur in the calculation.
References¶
[BH00]  T. Binoth and G. Heinrich,
An automatized algorithm to compute infrared divergent
multiloop integrals, Nucl. Phys. B 585 (2000) 741,

[BHJ+15]  S. Borowka, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk, T. Zirke,
SecDec3.0: numerical evaluation of multiscale integrals beyond one loop, 2015, Comput.Phys.Comm.196,

[BIR]  W. Bruns and B. Ichim and T. Römer and R. Sieg and C. Söger,
Normaliz. Algorithms for rational cones and affine monoids,
available at https://www.normaliz.uniosnabrueck.de

[BIS16]  W. Bruns, B. Ichim, C. Söger,
The power of pyramid decomposition in Normaliz, 2016, J.Symb.Comp.74, 513–536,

[Bor14]  S. Borowka,
Evaluation of multiloop multiscale integrals and phenomenological twoloop applications, 2014, PhD Thesis  Technische Universität München

[GKR+11]  J. Gluza, K. Kajda, T. Riemann, V. Yundin,
Numerical Evaluation of Tensor Feynman Integrals in Euclidean Kinematics, 2011, Eur.Phys.J.C71,

[GSL]  M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, P. Alken, M. Booth, F. Rossi,
GNU Scientific Library Reference Manual  Third Edition, 2009, Network Theory Ltd.,
ISBN: 0954612078 (ISBN13: 9780954612078),
available at http://www.gnu.org/software/gsl/

[Hah05]  T. Hahn,
CUBA: A Library for multidimensional numerical integration, 2005, Comput.Phys.Comm.168, 7895,

[Hah16]  T. Hahn,
Concurrent Cuba, 2016, Comput.Phys.Comm.207, 341349,

[Hei08]  G. Heinrich,
Sector Decomposition, 2008, Int.J.Mod.Phys.A23,

[KU10]  T. Kaneko and T. Ueda,
A Geometric method of sector decomposition, 2010, Comput.Phys.Comm.181,

[KUV13]  J. Kuipers, T. Ueda, J. A. M. Vermaseren,
Code Optimization in FORM, 2015, Comput.Phys.Comm.189, 119,

[LWY+15]  Z. Li, J. Wang, Q.S. Yan, X. Zhao,
Efficient Numerical Evaluation of Feynman Integrals, 2016, Chin.Phys.C40 No. 3, 033103,

[MP+14]  B. D. McKay and A. Piperno,
Practical graph isomorphism, II, 2014, Journal of Symbolic Computation, 60, 94112,

[Pak11]  A. Pak,
The toolbox of modern multiloop calculations: novel
analytic and semianalytic techniques, 2012, J. Phys.: Conf. Ser. 368 012049,

[PSD17]  S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk, T. Zirke,
pySecDec: A toolbox for the numerical evaluation of multiscale integrals, 2018, Comput.Phys.Comm.222,

[PSD18]  S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk,
A GPU compatible quasiMonte Carlo integrator interfaced to pySecDec,

[RUV17]  B. Ruijl, T. Ueda, J. Vermaseren,
FORM version 4.2,

[Ver00]  J. A. M. Vermaseren,
New features of FORM,
