Cuba

Cuba

Introduction

Cuba.jl is a Julia library for multidimensional numerical integration of real-valued functions of real arguments, using different algorithms.

This is just a Julia wrapper around the C Cuba library, version 4.2, by Thomas Hahn. All the credits goes to him for the underlying functions, blame me for any problem with the Julia interface.

All algorithms provided by Cuba library are supported in Cuba.jl:

For more details on the algorithms see the manual included in Cuba library and available in deps/usr/share/cuba.pdf after successful installation of Cuba.jl.

Integration is always performed on the $n$-dimensional unit hypercube $[0, 1]^{n}$.

Tip

If you want to compute an integral over a different set, you have to scale the integrand function in order to have an equivalent integral on $[0, 1]^{n}$ using substitution rules. For example, recall that in one dimension

\[\int_{a}^{b} f(x)\,\mathrm{d}x = \int_{0}^{1} f(a + (b - a) y) (b - a)\,\mathrm{d}y\]

where the final $(b - a)$ is the one-dimensional version of the Jacobian.

Integration over a semi-infinite or an inifite domain is a bit trickier, but you can follow this advice from Steven G. Johnson: to compute an integral over a semi-infinite interval, you can perform the change of variables $x=a+y/(1-y)$:

\[\int_{a}^{\infty} f(x)\,\mathrm{d}x = \int_{0}^{1} f\left(a + \frac{y}{1 - y}\right)\frac{1}{(1 - y)^2}\,\mathrm{d}y\]

For an infinite interval, you can perform the change of variables $x=(2y - 1)/((1 - y)y)$:

\[\int_{-\infty}^{\infty} f(x)\,\mathrm{d}x = \int_{0}^{1} f\left(\frac{2y - 1}{(1 - y)y}\right)\frac{2y^2 - 2y + 1}{(1 - y)^2y^2}\,\mathrm{d}y\]

In addition, recall that for an even function $\int_{-\infty}^{\infty} f(x)\,\mathrm{d}x = 2\int_{0}^{\infty}f(x)\,\mathrm{d}x$, while the integral of an odd function over the infinite interval $(-\infty, \infty)$ is zero.

All this generalizes straightforwardly to more than one dimension. In Examples section you can find the computation of a 3-dimensional integral with non-constant boundaries using Cuba.jl and two integrals over infinite domains.

Cuba.jl is available for GNU/Linux, FreeBSD, Mac OS, and Windows (i686 and x86_64 architectures).

Installation

The latest version of Cuba.jl is available for Julia 1.0 and later versions, and can be installed with Julia built-in package manager. In a Julia session run the commands

pkg> update
pkg> add Cuba

Older versions are also available for Julia 0.4-0.7.

Usage

After installing the package, run

using Cuba

or put this command into your Julia script.

Cuba.jl provides the following functions to integrate:

Cuba.vegas — Function.
vegas(integrand, ndim=1, ncomp=1[, keywords]) -> integral, error, probability, neval, fail, nregions

Calculate integral of integrand over the unit hypercube in ndim dimensions using Vegas algorithm. integrand is a vectorial function with ncomp components. ndim and ncomp default to 1.

Accepted keywords:

  • nvec
  • rtol
  • atol
  • flags
  • seed
  • minevals
  • maxevals
  • nstart
  • nincrease
  • nbatch
  • gridno
  • statefile
  • spin
source
Cuba.suave — Function.
suave(integrand, ndim=1, ncomp=1[, keywords]) -> integral, error, probability, neval, fail, nregions

Calculate integral of integrand over the unit hypercube in ndim dimensions using Suave algorithm. integrand is a vectorial function with ncomp components. ndim and ncomp default to 1.

Accepted keywords:

  • nvec
  • rtol
  • atol
  • flags
  • seed
  • minevals
  • maxevals
  • nnew
  • nmin
  • flatness
  • statefile
  • spin
source
Cuba.divonne — Function.
divonne(integrand, ndim=2, ncomp=1[, keywords]) -> integral, error, probability, neval, fail, nregions

Calculate integral of integrand over the unit hypercube in ndim dimensions using Divonne algorithm. integrand is a vectorial function with ncomp components. ncomp defaults to 1, ndim defaults to 2 and must be ≥ 2.

Accepted keywords:

  • nvec
  • rtol
  • atol
  • flags
  • seed
  • minevals
  • maxevals
  • key1
  • key2
  • key3
  • maxpass
  • border
  • maxchisq
  • mindeviation
  • ngiven
  • ldxgiven
  • xgiven
  • nextra
  • peakfinder
  • statefile
  • spin
source
Cuba.cuhre — Function.
cuhre(integrand, ndim=2, ncomp=1[, keywords]) -> integral, error, probability, neval, fail, nregions

Calculate integral of integrand over the unit hypercube in ndim dimensions using Cuhre algorithm. integrand is a vectorial function with ncomp components. ncomp defaults to 1, ndim defaults to 2 and must be ≥ 2.

Accepted keywords:

  • nvec
  • rtol
  • atol
  • flags
  • minevals
  • maxevals
  • key
  • statefile
  • spin
source

Large parts of the following sections are borrowed from Cuba manual. Refer to it for more information on the details.

Cuba.jl wraps the 64-bit integers functions of Cuba library, in order to push the range of certain counters to its full extent. In detail, the following arguments:

are passed to the Cuba library as 64-bit integers, so they are limited to be at most

julia> typemax(Int64)
9223372036854775807

There is no way to overcome this limit. See the following sections for the meaning of each argument.

Arguments

The only mandatory argument of integrator functions is:

Optional positional arguments are:

integrand should be a function integrand(x, f) taking two arguments:

x and f are matrices with dimensions (ndim, nvec) and (ncomp, nvec), respectively, when nvec > 1. See the Vectorization section below for more information.

Also anonymous functions are allowed as integrand. For those familiar with Cubature.jl package, this is the same syntax used for integrating vector-valued functions.

For example, the integral

\[\int_{0}^{1} \cos (x) \,\mathrm{d}x = \sin(1) = 0.8414709848078965\]

can be computed with one of the following commands

julia> vegas((x, f) -> f[1] = cos(x[1]))
Component:
 1: 0.8414910005259609 ± 5.2708169787733e-5 (prob.: 0.028607201257039333)
Integrand evaluations: 13500
Number of subregions:  0
Note: The desired accuracy was reached

julia> suave((x, f) -> f[1] = cos(x[1]))
Component:
 1: 0.8411523690658836 ± 8.357995611133613e-5 (prob.: 1.0)
Integrand evaluations: 22000
Number of subregions:  22
Note: The desired accuracy was reached

julia> divonne((x, f) -> f[1] = cos(x[1]))
Component:
 1: 0.841468071955942 ± 5.3955070531551656e-5 (prob.: 0.0)
Integrand evaluations: 1686
Number of subregions:  14
Note: The desired accuracy was reached

julia> cuhre((x, f) -> f[1] = cos(x[1]))
Component:
 1: 0.8414709848078966 ± 2.2204460420128823e-16 (prob.: 3.443539937576958e-5)
Integrand evaluations: 195
Number of subregions:  2
Note: The desired accuracy was reached

In section Examples you can find more complete examples. Note that x and f are both arrays with type Float64, so Cuba.jl can be used to integrate real-valued functions of real arguments. See how to work with a complex integrand.

Note: if you used Cuba.jl until version 0.0.4, be aware that the user interface has been reworked in version 0.0.5 in a backward incompatible way.

Optional Keywords

All other arguments required by Cuba integrator routines can be passed as optional keywords. Cuba.jl uses some reasonable default values in order to enable users to invoke integrator functions with a minimal set of arguments. Anyway, if you want to make sure future changes to some default values of keywords will not affect your current script, explicitely specify the value of the keywords.

Common Keywords

These are optional keywords common to all functions:

Vegas-Specific Keywords

These optional keywords can be passed only to vegas:

Suave-Specific Keywords

These optional keywords can be passed only to suave:

Divonne-Specific Keywords

These optional keywords can be passed only to divonne:

Cuhre-Specific Keyword

This optional keyword can be passed only to cuhre:

Output

The integrating functions vegas, suave, divonne, and cuhre return an Integral object whose fields are

integral :: Vector{Float64}
error    :: Vector{Float64}
probl    :: Vector{Float64}
neval    :: Int64
fail     :: Int32
nregions :: Int32

The first three fields are arrays with length ncomp, the last three ones are scalars. The Integral object can also be iterated over like a tuple. In particular, if you assign the output of integrator functions to the variable named result, you can access the value of the i-th component of the integral with result[1][i] or result.integral[i] and the associated error with result[2][i] or result.error[i].

Vectorization

Vectorization means evaluating the integrand function for several points at once. This is also known as Single Instruction Multiple Data (SIMD) paradigm and is different from ordinary parallelization where independent threads are executed concurrently. It is usually possible to employ vectorization on top of parallelization.

Cuba.jl cannot automatically vectorize the integrand function, of course, but it does pass (up to) nvec points per integrand call (Common Keywords). This value need not correspond to the hardware vector length –computing several points in one call can also make sense e.g. if the computations have significant intermediate results in common.

When nvec > 1, the input x is a matrix of dimensions (ndim, nvec), while the output f is a matrix with dimensions (ncomp, nvec). Vectorization can be used to evaluate more quickly the integrand function, for example by exploiting parallelism, thus speeding up computation of the integral. See the section Vectorized Function below for an example of a vectorized funcion.

Disambiguation

The nbatch argument of vegas is related in purpose but not identical to nvec. It internally partitions the sampling done by Vegas but has no bearing on the number of points given to the integrand. On the other hand, it it pointless to choose nvec > nbatch for Vegas.

Examples

One dimensional integral

The integrand of

\[\int_{0}^{1} \frac{\log(x)}{\sqrt{x}} \,\mathrm{d}x\]

has an algebraic-logarithmic divergence for $x = 0$, but the integral is convergent and its value is $-4$. Cuba.jl integrator routines can handle this class of functions and you can easily compute the numerical approximation of this integral using one of the following commands:

julia> vegas( (x,f) -> f[1] = log(x[1])/sqrt(x[1]))
Component:
 1: -3.9981623937128465 ± 0.00044066437168409464 (prob.: 0.2843052968819913)
Integrand evaluations: 1007500
Number of subregions:  0
Note: The accuracy was not met within the maximum number of evaluations

julia> suave( (x,f) -> f[1] = log(x[1])/sqrt(x[1]))
Component:
 1: -3.999976286717141 ± 0.00039504866661339003 (prob.: 1.0)
Integrand evaluations: 51000
Number of subregions:  51
Note: The desired accuracy was reached

julia> divonne( (x,f) -> f[1] = log(x[1])/sqrt(x[1]), atol = 1e-8, rtol = 1e-8)
Component:
 1: -3.999999899620808 ± 2.1865962888459237e-7 (prob.: 0.0)
Integrand evaluations: 1002059
Number of subregions:  1582
Note: The accuracy was not met within the maximum number of evaluations
Hint: Try increasing `maxevals` to 4884287

julia> cuhre( (x,f) -> f[1] = log(x[1])/sqrt(x[1]))
Component:
 1: -4.000000355067187 ± 0.0003395484028626406 (prob.: 0.0)
Integrand evaluations: 5915
Number of subregions:  46
Note: The desired accuracy was reached

Vector-valued integrand

Consider the integral

\[\int\limits_{\Omega} \boldsymbol{f}(x,y,z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\]

where $\Omega = [0, 1]^{3}$ and

\[\boldsymbol{f}(x,y,z) = \left(\sin(x)\cos(y)\exp(z), \,\exp(-(x^2 + y^2 + z^2)), \,\frac{1}{1 - xyz}\right)\]

In this case it is more convenient to write a simple Julia script to compute the above integral

julia> using Cuba, SpecialFunctions

julia> function integrand(x, f)
           f[1] = sin(x[1])*cos(x[2])*exp(x[3])
           f[2] = exp(-(x[1]^2 + x[2]^2 + x[3]^2))
           f[3] = 1/(1 - prod(x))
       end
integrand (generic function with 1 method)

julia> result, err = cuhre(integrand, 3, 3, atol=1e-12, rtol=1e-10);

julia> answer = ((ℯ-1)*(1-cos(1))*sin(1), (sqrt(pi)*erf(1)/2)^3, zeta(3));

julia> for i = 1:3
           println("Component ", i)
           println(" Result of Cuba: ", result[i], " ± ", err[i])
           println(" Exact result:   ", answer[i])
           println(" Actual error:   ", abs(result[i] - answer[i]))
       end
Component 1
 Result of Cuba: 0.6646696797813745 ± 1.0056262721114345e-13
 Exact result:   0.6646696797813771
 Actual error:   2.6645352591003757e-15
Component 2
 Result of Cuba: 0.41653838588064585 ± 2.932867102879894e-11
 Exact result:   0.41653838588663805
 Actual error:   5.992206730809357e-12
Component 3
 Result of Cuba: 1.2020569031649704 ± 1.1958521782293645e-10
 Exact result:   1.2020569031595951
 Actual error:   5.375255796025158e-12

Integral with non-constant boundaries

The integral

\[\int_{-y}^{y}\int_{0}^{z}\int_{0}^{\pi} \cos(x)\sin(y)\exp(z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\]

has non-constant boundaries. By applying the substitution rule repeatedly, you can scale the integrand function and get this equivalent integral over the fixed domain $\Omega = [0, 1]^{3}$

\[\int\limits_{\Omega} 2\pi^{3}yz^2 \cos(\pi yz(2x - 1)) \sin(\pi yz) \exp(\pi z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\]

that can be computed with Cuba.jl using the following Julia script

julia> using Cuba

julia> function integrand(x, f)
           f[1] = 2pi^3*x[2]*x[3]^2*cos(pi*x[2]*x[3]*(2*x[1] - 1.0))*
                  sin(pi*x[2]*x[3])*exp(pi*x[3])
       end
integrand (generic function with 1 method)

julia> result, err = cuhre(integrand, 3, 1, atol=1e-12, rtol=1e-10);

julia> answer = pi*ℯ^pi - (4ℯ^pi - 4)/5;

julia> begin
               println("Result of Cuba: ", result[1], " ± ", err[1])
               println("Exact result:   ", answer)
               println("Actual error:   ", abs(result[1] - answer))
       end
Result of Cuba: 54.98607586826155 ± 5.4606062189698135e-9
Exact result:   54.98607586789537
Actual error:   3.6617819887396763e-10

Integrals over Infinite Domains

Cuba.jl assumes always as integration domain the hypercube $[0, 1]^n$, but we have seen that using integration by substitution we can calculate integrals over different domains as well. In the Introduction we also proposed two useful substitutions that can be employed to change an infinite or semi-infinite domain into a finite one.

As a first example, consider the following integral with a semi-infinite domain:

\[\int_{0}^{\infty}\frac{\log(1 + x^2)}{1 + x^2}\,\mathrm{d}x\]

whose exact result is $\pi\log 2$. This can be computed as follows:

julia> using Cuba

julia> # The function we want to integrate over [0, ∞).

julia> func(x) = log(1 + x^2)/(1 + x^2)
func (generic function with 1 method)

julia> # Scale the function in order to integrate over [0, 1].

julia> function integrand(x, f)
           f[1] = func(x[1]/(1 - x[1]))/(1 - x[1])^2
       end
integrand (generic function with 1 method)

julia> result, err = cuhre(integrand, atol = 1e-12, rtol = 1e-10);

julia> answer = pi*log(2);

julia> begin
               println("Result of Cuba: ", result[1], " ± ", err[1])
               println("Exact result:   ", answer)
               println("Actual error:   ", abs(result[1] - answer))
       end
Result of Cuba: 2.1775860903056885 ± 2.150398850102772e-10
Exact result:   2.177586090303602
Actual error:   2.086331107875594e-12

Now we want to calculate this integral, over an infinite domain

\[\int_{-\infty}^{\infty} \frac{1 - \cos x}{x^2}\,\mathrm{d}x\]

which gives $\pi$. You can calculate the result with the code below. Note that integrand function has value $1/2$ for $x=0$, but you have to inform Julia about this.

julia> using Cuba

julia> # The function we want to integrate over (-∞, ∞).

julia> func(x) = x==0 ? 0.5*one(x) : (1 - cos(x))/x^2
func (generic function with 1 method)

julia> # Scale the function in order to integrate over [0, 1].

julia> function integrand(x, f)
           f[1] = func((2*x[1] - 1)/x[1]/(1 - x[1])) *
                   (2*x[1]^2 - 2*x[1] + 1)/x[1]^2/(1 - x[1])^2
       end
integrand (generic function with 1 method)

julia> result, err = cuhre(integrand, atol = 1e-7, rtol = 1e-7);

julia> answer = float(pi);

julia> begin
               println("Result of Cuba: ", result[1], " ± ", err[1])
               println("Exact result:   ", answer)
               println("Actual error:   ", abs(result[1] - answer))
       end
Result of Cuba: 3.1415928900555046 ± 2.050669142074607e-6
Exact result:   3.141592653589793
Actual error:   2.3646571145619077e-7

Complex integrand

As already explained, Cuba.jl operates on real quantities, so if you want to integrate a complex-valued function of complex arguments you have to treat complex quantities as 2-component arrays of real numbers. For example, if you do not remember Euler's formula, you can compute this simple integral

\[\int_{0}^{\pi/2} \exp(\mathrm{i} x)\,\mathrm{d}x\]

with the following code

julia> using Cuba

julia> function integrand(x, f)
           # Complex integrand, scaled to integrate in [0, 1].
           tmp = cis(x[1]*pi/2)*pi/2
           # Assign to two components of "f" the real
           # and imaginary part of the integrand.
           f[1], f[2] = reim(tmp)
       end
integrand (generic function with 1 method)

julia> result = cuhre(integrand, 2, 2);

julia> begin
           println("Result of Cuba: ", complex(result[1]...))
           println("Exact result:   ", complex(1.0, 1.0))
       end
Result of Cuba: 1.0 + 1.0im
Exact result:   1.0 + 1.0im

Passing data to the integrand function

Cuba Library allows program written in C and Fortran to pass extra data to the integrand function with userdata argument. This is useful, for example, when the integrand function depends on changing parameters. In Cuba.jl the userdata argument is not available, but you do not normally need it.

For example, the cumulative distribution function $F(x;k)$ of chi-squared distribution is defined by

\[F(x; k) = \int_{0}^{x} \frac{t^{k/2 - 1}\exp(-t/2)}{2^{k/2}\Gamma(k/2)} \,\mathrm{d}t\]

The cumulative distribution function depends on parameter $k$, but the function passed as integrand to Cuba.jl integrator routines accepts as arguments only the input and output vectors. However you can easily define a function to calculate a numerical approximation of $F(x; k)$ based on the above integral expression because the integrand can access any variable visible in its scope. The following Julia script computes $F(x = \pi; k)$ for different $k$ and compares the result with more precise values, based on the analytic expression of the cumulative distribution function, provided by GSL.jl package.

julia> using Cuba, GSL, Printf, SpecialFunctions

julia> function chi2cdf(x::Real, k::Real)
           k2 = k/2
           # Chi-squared probability density function, without constant denominator.
           # The result of integration will be divided by that factor.
           function chi2pdf(t::Float64)
               # "k2" is taken from the outside.
               return t^(k2 - 1.0)*exp(-t/2)
           end
           # Neither "x" is passed directly to the integrand function,
           # but is visible to it.  "x" is used to scale the function
           # in order to actually integrate in [0, 1].
           x*cuhre((t,f) -> f[1] = chi2pdf(t[1]*x))[1][1]/(2^k2*gamma(k2))
       end
chi2cdf (generic function with 1 method)

julia> x = float(pi);

julia> begin
            @printf("Result of Cuba: %.6f %.6f %.6f %.6f %.6f\n",
                    map((k) -> chi2cdf(x, k), collect(1:5))...)
            @printf("Exact result:   %.6f %.6f %.6f %.6f %.6f\n",
                    map((k) -> cdf_chisq_P(x, k), collect(1:5))...)
        end
Result of Cuba: 0.923681 0.792120 0.629694 0.465584 0.321833
Exact result:   0.923681 0.792120 0.629695 0.465584 0.321833

Vectorized Function

Consider the integral

\[\int\limits_{\Omega} \prod_{i=1}^{10} \cos(x_{i}) \,\mathrm{d}\boldsymbol{x} = \sin(1)^{10} = 0.1779883\dots\]

where $\Omega = [0, 1]^{10}$ and $\boldsymbol{x} = (x_{1}, \dots, x_{10})$ is a 10-dimensional vector. A simple way to compute this integral is the following:

julia> using Cuba, BenchmarkTools

julia> cuhre((x, f) -> f[] = prod(cos.(x)), 10)
Component:
 1: 0.1779870665870775 ± 1.0707995959536173e-6 (prob.: 0.2438374075714901)
Integrand evaluations: 7815
Fail:                  0
Number of subregions:  2

julia> @benchmark cuhre((x, f) -> f[] = prod(cos.(x)), 10)
BenchmarkTools.Trial:
  memory estimate:  2.62 MiB
  allocs estimate:  39082
  --------------
  minimum time:     1.633 ms (0.00% GC)
  median time:      1.692 ms (0.00% GC)
  mean time:        1.867 ms (8.62% GC)
  maximum time:     3.660 ms (45.54% GC)
  --------------
  samples:          2674
  evals/sample:     1

We can use vectorization in order to speed up evaluation of the integrand function.

julia> function fun_vec(x,f)
           f[1,:] .= 1.0
           for j in 1:size(x,2)
               for i in 1:size(x, 1)
                   f[1, j] *= cos(x[i, j])
               end
           end
       end
fun_vec (generic function with 1 method)

julia> cuhre(fun_vec, 10, nvec = 1000)
Component:
 1: 0.1779870665870775 ± 1.0707995959536173e-6 (prob.: 0.2438374075714901)
Integrand evaluations: 7815
Fail:                  0
Number of subregions:  2

julia> @benchmark cuhre(fun_vec, 10, nvec = 1000)
BenchmarkTools.Trial:
  memory estimate:  2.88 KiB
  allocs estimate:  54
  --------------
  minimum time:     949.976 μs (0.00% GC)
  median time:      954.039 μs (0.00% GC)
  mean time:        966.930 μs (0.00% GC)
  maximum time:     1.204 ms (0.00% GC)
  --------------
  samples:          5160
  evals/sample:     1

A further speed up can be gained by running the for loop in parallel with Threads.@threads. For example, running Julia with 4 threads:

julia> function fun_par(x,f)
           f[1,:] .= 1.0
           Threads.@threads for j in 1:size(x,2)
               for i in 1:size(x, 1)
                   f[1, j] *= cos(x[i, j])
               end
           end
       end
fun_par (generic function with 1 method)

julia> cuhre(fun_par, 10, nvec = 1000)
Component:
 1: 0.1779870665870775 ± 1.0707995959536173e-6 (prob.: 0.2438374075714901)
Integrand evaluations: 7815
Fail:                  0
Number of subregions:  2

julia> @benchmark cuhre(fun_par, 10, nvec = 1000)
BenchmarkTools.Trial:
  memory estimate:  3.30 KiB
  allocs estimate:  63
  --------------
  minimum time:     507.914 μs (0.00% GC)
  median time:      515.182 μs (0.00% GC)
  mean time:        520.667 μs (0.06% GC)
  maximum time:     3.801 ms (85.06% GC)
  --------------
  samples:          9565
  evals/sample:     1

Performance

Cuba.jl cannot (yet?) take advantage of parallelization capabilities of Cuba Library. Nonetheless, it has performances comparable with equivalent native C or Fortran codes based on Cuba library when CUBACORES environment variable is set to 0 (i.e., multithreading is disabled). The following is the result of running the benchmark present in test directory on a 64-bit GNU/Linux system running Julia 0.7.0-beta2.3 (commit 83ce9c7524) equipped with an Intel(R) Core(TM) i7-4700MQ CPU. The C and FORTRAN 77 benchmark codes have been compiled with GCC 7.3.0.

$ CUBACORES=0 julia -e 'using Pkg; cd(Pkg.dir("Cuba")); include("test/benchmark.jl")'
[ Info: Performance of Cuba.jl:
  0.257360 seconds (Vegas)
  0.682703 seconds (Suave)
  0.329552 seconds (Divonne)
  0.233190 seconds (Cuhre)
[ Info: Performance of Cuba Library in C:
  0.268249 seconds (Vegas)
  0.682682 seconds (Suave)
  0.319553 seconds (Divonne)
  0.234099 seconds (Cuhre)
[ Info: Performance of Cuba Library in Fortran:
  0.233532 seconds (Vegas)
  0.669809 seconds (Suave)
  0.284515 seconds (Divonne)
  0.195740 seconds (Cuhre)

Of course, native C and Fortran codes making use of Cuba Library outperform Cuba.jl when higher values of CUBACORES are used, for example:

$ CUBACORES=1 julia -e 'using Pkg; cd(Pkg.dir("Cuba")); include("test/benchmark.jl")'
[ Info: Performance of Cuba.jl:
  0.260080 seconds (Vegas)
  0.677036 seconds (Suave)
  0.342396 seconds (Divonne)
  0.233280 seconds (Cuhre)
[ Info: Performance of Cuba Library in C:
  0.096388 seconds (Vegas)
  0.574647 seconds (Suave)
  0.150003 seconds (Divonne)
  0.102817 seconds (Cuhre)
[ Info: Performance of Cuba Library in Fortran:
  0.094413 seconds (Vegas)
  0.556084 seconds (Suave)
  0.139606 seconds (Divonne)
  0.107335 seconds (Cuhre)

Cuba.jl internally fixes CUBACORES to 0 in order to prevent from forking julia processes that would only slow down calculations eating up the memory, without actually taking advantage of concurrency. Furthemore, without this measure, adding more Julia processes with addprocs() would only make the program segfault.

Related projects

There are other Julia packages for multidimenensional numerical integration:

Development

Cuba.jl is developed on GitHub: https://github.com/giordano/Cuba.jl. Feel free to report bugs and make suggestions at https://github.com/giordano/Cuba.jl/issues.

History

The ChangeLog of the package is available in NEWS.md file in top directory. There have been some breaking changes from time to time, beware of them when upgrading the package.

License

The Cuba.jl package is licensed under the GNU Lesser General Public License, the same as Cuba library. The original author is Mosè Giordano.

Credits

If you use this library for your work, please credit Thomas Hahn. Citable papers about Cuba Library: