2. 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.
2.1. 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:
$ python3 generate_easy.py
$ make C easy
$ python3 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.00000000000000022e+00,0.00000000000000000e+00) +/ (5.65352153979095401e17,0.00000000000000000e+00))*eps^1 + ((3.06852819440053548e01,0.00000000000000000e+00) +/ (1.18502493127591741e15,0.00000000000000000e+00)) + 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.
#!/usr/bin/env python3
from pySecDec import make_package
if __name__ == "__main__":
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.
#!/usr/bin/env python3
from pySecDec.integral_interface import IntegralLibrary
from math import log
if __name__ == "__main__":
# 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)))
2.2. 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:
$ python3 generate_box1L.py
$ make C box1L
$ python3 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.
2.2.1. 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
from pySecDec.loop_integral import loop_package
import pySecDec as psd
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_orders 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_orders = [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',
)
2.2.2. 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 README box1L.pdf box1L_integral integral_names.txt pylink
Makefile.conf box1L.hpp box1L_coefficients integrate_box1L.cpp src
in the folder box1L, typing
$ make
will create the static library box1L_integral/libbox1L_integral.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
To build the library with nvcc for GPU support, type
$ CXX=nvcc SECDEC_WITH_CUDA_FLAGS="arch=sm_XX" make
where sm_XX
must be replaced by the target GPU architechtures, see the arch option of NVCC.
The SECDEC_WITH_CUDA_FLAGS
environment variable, which enables GPU code compilation, contains flags which are passed to NVCC during code compilation and linking.
Multiple GPU architectures may be specified as described in the NVCC manual, for example
SECDEC_WITH_CUDA_FLAGS="gencode arch=compute_XX,code=sm_XX gencode arch=compute_YY,code=sm_YY"
where XX
and YY
are the target GPU architectures. The script
examples/easy/printcudaarch.sh
can be used to obtain the compute architecture of your current machine.
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.
2.2.3. 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
box1L.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
box1L.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
box1L.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 integrate_box1L_multiple_points.py.
2.2.4. 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.
After the lines parsing the input parameters, an secdecutil::Integrator
is constructed and its parameters are set:
// Set up Integrator
secdecutil::integrators::Qmc<
box1L::integrand_return_t,
box1L::maximal_number_of_integration_variables,
integrators::transforms::Korobov<3>::type,
box1L::user_integrand_t
> integrator;
integrator.verbosity = 1;
The amplitude is constructed via a call to name::make_amplitudes()
and packed into a name::handler_t
.
// Construct the amplitudes
std::vector<box1L::nested_series_t<box1L::sum_t>> unwrapped_amplitudes =
box1L::make_amplitudes(real_parameters, complex_parameters, "box1L_coefficients", integrator);
// Pack amplitudes into handler
box1L::handler_t<box1L::amplitudes_t> amplitudes
(
unwrapped_amplitudes,
integrator.epsrel, integrator.epsabs
// further optional arguments: maxeval, mineval, maxincreasefac, min_epsrel, min_epsabs, max_epsrel, max_epsabs
);
amplitudes.verbose = true;
If desired, the contour deformation can be adjusted via additional arguments to name::handler_t
.
See also
Section 4.1 and Section 5.9.1 for more detailed information about name::make_amplitudes()
and name::handler_t
.
To numerically integrate the sum of sectors, the name::handler_t::evaluate()
function is called:
// compute the amplitudes
const std::vector<box1L::nested_series_t<secdecutil::UncorrelatedDeviation<box1L::integrand_return_t>>> result = amplitudes.evaluate();
The remaining lines print the result:
// print the result
for (unsigned int amp_idx = 0; amp_idx < box1L::number_of_amplitudes; ++amp_idx)
std::cout << "amplitude" << amp_idx << " = " << result.at(amp_idx) << std::endl;
The C++ program can be built with the command:
$ make integrate_box1L
A kinematic point must be specified when calling the integrate_box1L
executable, the input format is:
$ ./integrate_box1L 4.0 0.75 1.25 1.0
where the arguments are the real_parameters
values for (s
, t
, s1
, msq
).
For integrals depending on complex_parameters
, their value is specified by a space separated pair of numbers representing the real and imaginary part.
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 QMC integrator we refer to Section 4.6.2.
2.3. Evaluating a Weighted Sum of Integrals¶
New in version 1.5.
Let us examine example easy_sum
, which demonstrates how two weighted sums of dimensionally regulated integrals can be evaluated.
The example computes the following two weighted sums:
where
First, we import the necessary python packages and open the if __name__ == "__main__"
guard, as required by multiprocessing.
#!/usr/bin/env python3
from pySecDec import Coefficient
from pySecDec import MakePackage
from pySecDec import sum_package
if __name__ == "__main__":
The common arguments for the integrals are collected in the common_args
dictionary.
common_args = {}
common_args['real_parameters'] = ['s']
common_args['regulators'] = ['eps']
common_args['requested_orders'] = [0]
Next, the coefficients of the integrals for each weighted sum are specified.
Each Coefficient
is specified as a list of numerator factors, list of denominator factors and a list of real or complex parameters on which the coefficient depends.
Coefficients can depend also depend the regulators, the sum_package
function will automatically determine the correct orders to which the coefficients and integrals should be expanded in order to obtain the requested_orders
.
coefficients = [
[ # sum1
Coefficient(['2*s'],['1'],['s']), # easy1
Coefficient(['3*s'],['1'],['s']) # easy2
],
[ # sum2
Coefficient(['s'],['2*eps'],['s']), # easy1
Coefficient(['s*eps'],['3'],['s']) # easy2
]
]
The integrals are specified using the MakePackage wrapper function (which has the same arguments as make_package
), for loop integrals the LoopPackage wrapper may be used (it has the same arguments as loop_package
).
integrals = [
MakePackage('easy1',
integration_variables = ['x','y'],
polynomials_to_decompose = ['(x+y)^(2+eps)'],
**common_args),
MakePackage('easy2',
integration_variables = ['x','y'],
polynomials_to_decompose = ['(2*x+3*y)^(1+eps)'],
**common_args)
]
Finally, the list of integrals and coefficients are passed to sum_package
. This will generate a C++ library which efficiently evaluates both weighted sums of integrals, sharing the results of the integrals between the different sums.
# generate code sum of (int * coeff)
sum_package('easy_sum', integrals,
coefficients = coefficients, **common_args)
The generated C++ library can be compiled and called via the python and/or C++ interface as described above.
2.4. Using Expansion By Regions (Generic Integral)¶
New in version 1.5.
The example make_regions_ebr
provides a simple introduction to the expansion by regions functionality within pySecDec.
For a more detailed discussion of expansion by regions see our paper [PSD21].
The necessary packages are loaded and the if __name__ == "__main__"
guard is opened.
#!/usr/bin/env python3
from pySecDec import sum_package, make_regions
if __name__ == "__main__":
Expansion by regions is applied to a generic integral using the make_regions
function.
regions_generators = make_regions(
name = 'make_regions_ebr',
integration_variables = ['x'],
regulators = ['delta'],
requested_orders = [0],
smallness_parameter = 't',
polynomials_to_decompose = ['(x)**(delta)','(t + x + x**2)**(1)'],
expansion_by_regions_order = 0,
real_parameters = ['t'],
complex_parameters = [],
decomposition_method = 'geometric_infinity_no_primary',
polytope_from_sum_of=[1]
)
The output of make_regions
can be passed to sum_package
in order to generate a C++ library suitable for evaluating the expanded integral.
sum_package(
'make_regions_ebr',
regions_generators,
regulators = ['delta'],
requested_orders = [0],
real_parameters = ['t']
)
The generated C++ library can be compiled and called via the python and/or C++ interface as described above.
2.5. Using Expansion By Regions (Loop Integral)¶
New in version 1.5.
The example generate_box1L_ebr
demonstrates how expansion by regions can be applied to loop integrals within pySecDec by applying it to the 1loop box integral as described in Section 4.2 of [Mis18].
For a more detailed discussion of expansion by regions see our paper [PSD21].
First, the necessary packages are loaded and the if __name__ == "__main__"
guard is opened.
#!/usr/bin/env python3
from pySecDec import sum_package, loop_regions
import pySecDec as psd
# This example is the oneloop box example in Go Mishima's paper arXiv:1812.04373
if __name__ == "__main__":
The loop integral can be constructed via the convenience functions in loop_integral
, here we use LoopintegralFromGraph
.
Note that powerlist=["1+n1","1+n1/2","1+n1/3","1+n1/5"]
, here n1
is an extra regulator required to regulate the singularities which appear when expanding this loop integral.
We use the “trick” of introducing only a single regulator divided by different prime numbers for each power, rather than unique regulators for each propagator (though this is also supported by pySecDec).
Poles in the extra regulator n1
may appear in individual regions but are expected to cancel when all regions are summed.
# here we define the Feynman diagram
li = psd.loop_integral.LoopIntegralFromGraph(
internal_lines = [['mt',[3,1]],['mt',[1,2]],['mt',[2,4]],['mt',[4,3]]],
external_lines = [['p1',1],['p2',2],['p3',3],['p4',4]],
powerlist=["1+n1","1+n1/2","1+n1/3","1+n1/5"],
regulators=["eps","n1"],
Feynman_parameters=["x%i" % i for i in range(1,5)], # renames the parameters to get the same polynomials as in 1812.04373
replacement_rules = [
# note that in those relations all momenta are incoming
# general relations:
('p1*p1', 'm1sq'),
('p2*p2', 'm2sq'),
('p3*p3', 'm3sq'),
('p4*p4', 'm4sq'),
('p1*p2', 's/2(m1sq+m2sq)/2'),
('p1*p3', 't/2(m1sq+m3sq)/2'),
('p1*p4', 'u/2(m1sq+m4sq)/2'),
('p2*p3', 'u/2(m2sq+m3sq)/2'),
('p2*p4', 't/2(m2sq+m4sq)/2'),
('p3*p4', 's/2(m3sq+m4sq)/2'),
('u', '(m1sq+m2sq+m3sq+m4sq)st'),
# relations for our specific case:
('mt**2', 'mtsq'),
('m1sq',0),
('m2sq',0),
('m3sq','mHsq'),
('m4sq','mHsq'),
('mHsq', 0),
])
Expansion by regions is applied to a loop integral using the loop_regions
function.
We expand around a small mass mtsq.
# find the regions
generators_args = loop_regions(
name = "box1L_ebr",
loop_integral=li,
smallness_parameter = "mtsq",
expansion_by_regions_order=0)
The output of loop_regions
can be passed to sum_package
in order to generate a C++ library suitable for evaluating the expanded integral.
# write the code to sum up the regions
sum_package("box1L_ebr",
generators_args,
li.regulators,
requested_orders = [0,0],
real_parameters = ['s','t','u','mtsq'],
complex_parameters = [])
The generated C++ library can be compiled and called via the python and/or C++ interface as described above.
2.6. List of Examples¶
Here we list the available examples. For more details regarding each example see [PSD17], [PSD18] and [PSD21].
easy: 
a simple parametric integral, described in Section 2.1 
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

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 
hypergeo5F4: 
a general dimensionally regulated parameter integral 
hz2L_nonplanar: 
a 2loop, 4point, 7propagator integral with internal and external masses 
box1L_ebr: 
uses expansion by regions to expand a 1loop box with a small internal mass, this integral is also considered in Section 4.2 of [Mis18] 
bubble1L_ebr: 
uses expansion by regions to expand a 1loop, 2point integral in various limits, demonstrates the use of an additional regulator as described in [PSD21] 
bubble1L_dotted_ebr: 
uses expansion by regions to expand a 1loop, 2point integral, demonstrates the and methods described in [PSD21] 
bubble2L_largem_ebr: 
uses expansion by regions to expand a 1loop, 2point integral with a large mass 
bubble2L_smallm_ebr: 
uses expansion by regions to expand a 1loop, 2point integral with a small mass 
formfactor1L_ebr: 
uses expansion by regions to compute various 1loop, 3point form factor integrals from the literature, demonstrates the use
of 
triangle2L_ebr: 
uses expansion by regions to compute a 2loop, 3point integral with a large mass 
make_regions_ebr: 
uses expansion by regions to compute a simple generic integral with a small parameter 
easy_sum: 
calculates the sum of two integrals with different coefficients, demonstrates the use of 
yyyy1L: 
calculates a 1loop 4photon helicity amplitude, demonstrates the use of 
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 
regions: 
prints a list of the regions obtained by applying expansion by regions to formfactor1L_massless 