4. SecDecUtil

SecDecUtil is a standalone autotools-c++ 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.

4.1. 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 ith power of the expansion parameter. Otherwise elements can be accessed identically to std::vector.

template<typename T>
class Series
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.

Series(int order_min, int order_max, std::vector<T> content, bool truncated_above = true, const std::string expansion_parameter = "x")

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

4.2. 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.

template<typename Tout, typename Tin, template<typename...> class Tnest>
Tnest<Tout> deep_apply(Tnest<Tin> &nest, std::function<Tout(Tin)> &func)

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

4.3. 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,

\mathrm{E}\left[\frac{1}{X}\right] = \int_{-\infty}^{\infty} \frac{1}{X} p(X)\ \mathrm{d}X,

where

p(X) = \frac{1}{\sqrt{2 \pi \sigma^2 }} \exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right),

is undefined in the Riemann or Lebesgue sense. The rule \delta(a/b) = |a/b| \sqrt{ (\delta a/a)^2 + (\delta b/b)^2 } can not be derived from the first principles of probability theory.

The rules implemented for real valued error propagation are:

\delta(a+b) = \sqrt{(\delta a)^2 + (\delta b)^2},

\delta(a-b) = \sqrt{(\delta a)^2 + (\delta b)^2},

\delta(ab) = \sqrt{ (\delta a)^2 b^2 + (\delta b)^2 a^2 + (\delta a)^2 (\delta b)^2 }.

For complex numbers the above rules are implemented for the real and imaginary parts individually.

template<typename T>
class UncorrelatedDeviation
T value

The expectation value.

T uncertainty

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

4.4. 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 (+, -, *, /).

template<typename T, typename ...Args>
class IntegrandContainer
int number_of_integration_variables

The number of integration variables that the integrand depends on.

std::function<T(Args...)> integrand

The integrand function.

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

4.5. Integrator

A base class template from which integrator implementations inherit. It defines the minimal API available for all integrators.

template<typename return_t, typename input_t>
class Integrator
bool together

(Only available if return_t is a std::complex type) If true after each call of the function both the real and imaginary parts are passed to the underlying integrator. If false 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 an UncorrelatedDeviation.

4.5.1. 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: 1e-7.
  • flags - Sets the Cuba verbosity flags. The flags=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.

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.

4.5.2. Examples

4.5.2.1. 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 = 1e-4;
    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.4515e-05

4.5.2.2. 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)