giordano / cuba.jl

Library for multidimensional numerical integration with four independent algorithms: Vegas, Suave, Divonne, and Cuhre.

Cuba.jl

Documentation	Build Status	Code Coverage

Introduction

Cuba.jl is a library for multidimensional numerical integration with different algorithms in Julia.

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. Feel free to report bugs and make suggestions at https://github.com/giordano/Cuba.jl/issues.

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

• vegas (type: Monte Carlo; variance reduction with importance sampling)
• suave (type: Monte Carlo; variance reduction with globally adaptive subdivision + importance sampling)
• divonne (type: Monte Carlo or deterministic; variance reduction with stratified sampling, aided by methods from numerical optimization)
• cuhre (type: deterministic; variance reduction with globally adaptive subdivision)

Integration is performed on the n-dimensional unit hypercube [0, 1]^n. 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.

Cuba.jl is available on all platforms supported by Julia.

Installation

The latest version of Cuba.jl is available for Julia 1.3 and later versions, and can be installed with Julia built-in package manager. In a Julia session, after entering the package manager mode with ], run the command

pkg> update
pkg> add Cuba

Older versions are also available for Julia 0.4-1.2.

Usage

After installing the package, run

using Cuba

or put this command into your Julia script.

Cuba.jl provides the following functions to integrate:

vegas(integrand, ndim, ncomp[; keywords...])
suave(integrand, ndim, ncomp[; keywords...])
divonne(integrand, ndim, ncomp[; keywords...])
cuhre(integrand, ndim, ncomp[; keywords...])

These functions wrap the 64-bit integers functions provided by the Cuba library.

The only mandatory argument is:

• function: the name of the function to be integrated

Optional positional arguments are:

• ndim: the number of dimensions of the integration domain. Defaults to 1 in vegas and suave, to 2 in divonne and cuhre. Note: ndim must be at least 2 with the latest two methods.
• ncomp: the number of components of the integrand. Defaults to 1

ndim and ncomp arguments must appear in this order, so you cannot omit ndim but not ncomp. integrand should be a function integrand(x, f) taking two arguments:

• the input vector x of length ndim
• the output vector f of length ncomp, used to set the value of each component of the integrand at point x

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

∫_0^1 cos(x) dx = 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

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 vectors 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 integration 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]. The details of other quantities can be found in Cuba manual.

All other arguments listed in Cuba documentation can be passed as optional keywords.

Documentation

A more detailed manual of Cuba.jl, with many complete examples, is available at https://giordano.github.io/Cuba.jl/stable/.

Related projects

There are other Julia packages for multidimenensional numerical integration:

• Cubature.jl
• HCubature.jl
• NIntegration.jl

License

The Cuba.jl package is released under the terms of the MIT "Expat" License. Note that the binary library Cuba is distributed with the GNU Lesser General Public License. The original author of Cuba.jl is Mosè Giordano. If you use this library for your work, please credit Thomas Hahn (citable papers about Cuba library: https://ui.adsabs.harvard.edu/abs/2005CoPhC.168...78H/abstract and https://ui.adsabs.harvard.edu/abs/2015JPhCS.608a2066H/abstract).

cuba.jl's Issues

Controlling Size of X Matrix using nvec

Hello,

I'm trying to use the nvec option to choose a large matrix for sampling. However, when I set nvec = 10000 the integration function seems to generate a matrix that's only Mx 500 when it should be M x nvec. Is there anyway to ensure that the nvec option is implemented exactly? Why does it deviate from the selected nvec?

Code below:

import Pkg
Pkg.activate("joint_timing")
Pkg.instantiate()
using Cuba, Distributions
using BenchmarkTools, Test, CUDA
using FLoops, FoldsCUDA
using SpecialFunctions
using Suppressor

# User Inputs
M= 5 # number of independent uniform random variables
atol=1e-10
rtol=1e-10
nvec=100000
maxevals=1000000

## Setting up functions 
function beta_pdf_cpu(x,a,b,dim)
    prod(x.^(a-1.0f0) .* (1.0f0 .- x).^(b-1.0f0)./(gamma(a)*gamma(b)/gamma(a+b)),dims=dim)
end

function int_cpu(x, f)
display(size(x)) 
  f[1] = beta_pdf_cpu(x,1.0,2.0,1)[1]
end


function int_thread_col_cpu(x, f)
    display(size(x))
    @floop for i in 1:size(x,2)
         f[i] = beta_pdf_cpu(x[:,i],1.0,2.0,1)[1]
    end
end

#Integrating   

cuhre(int_cpu, M, 1, atol=atol, rtol=rtol,nvec = nvec)
suave(int_thread_col_cpu, M, 1, atol=atol, rtol=rtol, nvec = nvec, maxevals = maxevals)

## You can see here that the size of X is only 5x500 when it should be 5xNvec = 5 x 100000

Regions -> True inside Cuba.jl?

In the Mathematica interface to Cuba, I can obtain a list of the regions used by option Regions -> True. Is there any way to get this information in Julia?

After skimming the web page, it's unclear to me whether the C++ library exports this data, so perhaps this cannot be done.

suave not handling large maxevals

Hello, the docs state:

nvec, minevals, maxevals, nnew, nmin, neval are passed to the Cuba library as 64-bit integers, so they are limited to be at most typemax(Int64)

However, suave() seems to have issues returning NaN for large maxevals < typemax(Int64). I have not observed this issue with the other integrators in Cuba.jl

MWE:

suave((x,f)->f.=3.14, 1, 1, 
        rtol=1e-3, atol=1e-3, 
        nvec=1, maxevals = typemax(Int64)*.499 #works
suave((x,f)->f.=3.14, 1, 1, 
        rtol=1e-3, atol=1e-3, 
        nvec=1, maxevals = typemax(Int64)*.5) # returns NaN or random small number

Run benchmark codes

Hi, Mosè. I have a problem in running the benchmark test on julia-1.5.2.

I follow the instructions in the documentation (https://giordano.github.io/Cuba.jl/latest/#Performance-1).

After running "CUBACORES=0 julia -e 'using Pkg; import Cuba; include(joinpath(dirname(dirname(pathof(Cuba))), "benchmarks", "benchmark.jl"))'
" in Terminal, it returns

"
[ Info: Performance of Cuba.jl:
0.383992 seconds (Vegas)
0.572846 seconds (Suave)
0.466027 seconds (Divonne)
0.336353 seconds (Cuhre)
[ Info: Performance of Cuba Library in C:
ERROR: LoadError: UndefVarError: LIBPATH_env not defined
"

I have tried to run the benchmark code on MacOS and Debian. Both of them fail with the same error message. I guess it might be related to the fact that the benchmark code was written for Julia 0.7.0-beta2.3 as mentioned in the documentation. If possible, can you check this problem?

divonne and cuhre pass vector of length 2 for one-dimensional integration

When integrating in one dimension, vegas and suave call the integrand with a vectir of length 1 (as expected), but divonne and cuhre pass a vector of length 2:

julia> vegas((x, f) -> ((@assert length(x) == 1); f[1] = norm(x)), 1)
...

julia> suave((x, f) -> ((@assert length(x) == 1); f[1] = norm(x)), 1)
...

julia> divonne((x, f) -> ((@assert length(x) == 1); f[1] = norm(x)), 1)
ERROR: AssertionError: length(x) == 1
julia> divonne((x, f) -> ((@assert length(x) == 2); f[1] = norm(x)), 1)
...

julia> cuhre((x, f) -> ((@assert length(x) == 1); f[1] = norm(x)), 1)
ERROR: AssertionError: length(x) == 1
julia> cuhre((x, f) -> ((@assert length(x) == 2); f[1] = norm(x)), 1)
...

Request: better user interface for integration over a non-unit cube

I think the interface would be much more intuitive if the users can directly specify the integration range, e.g. cuhre(f::Function, xmin, xmax), similar as the interface in Cubature.jl

Improve user interface

Currently, users have to write their own integrand functions using this awkward boilerplate:

function integrand(ndim::Cint, xx::Ptr{Cdouble}, ncomp::Cint, ff::Ptr{Cdouble},
                   userdata::Ptr{Void})
    # Take arrays from "xx" and "ff" pointers.
    x = pointer_to_array(xx, (ndim,))
    f = pointer_to_array(ff, (ncomp,))
    # Do calculations on "f" here
    #   ...
    # Store back the results to "ff"
    ff = pointer_from_objref(f)
return Cint(0)::Cint
end

It would be better to let users directly call integrator functions with something like this:

Vegas( (x)-> cos(x), ndim[; keywords])

Commit e61c7f1 in branch ui enables this feature, but that implementation (using eval) is really _slow_ ever for simple integrand functions

Combining Cuba.jl with GPU

Hello,

I'm trying to utilize gpu computation using Cuda.jl to speed up calculating integrals. Is it possible to do so? If so, how? Here is my example code:

using Cuba, Distributions
using BenchmarkTools, Test, CUDA
using FLoops, FoldsCUDA
using SpecialFunctions

# user input
M= 10  #number of independent uniform random variables
atol=1e-6
rtol=1e-3
nvec=1000000
maxevals=100000000

# setting up integration function
function int_thread_col_gpu_precompile_test(x, f)
    w = CuArray(x)
   @floop CUDAEx() for i in 1:size(x,2)
  #@floop for i in 1:size(x,2)
      val = reduce(*,@view(w[:,i]))
        f[i] = val
    end
    return f
end

cuhre(int_thread_col_gpu_precompile_test, M, 1, atol=atol, rtol=rtol)
suave(int_thread_col_gpu_precompile_test, M, 1, atol=atol, rtol=rtol, nvec = nvec, maxevals = maxevals)

If I use the CUDAEx() call, the code errors. If I don't the code works fine, but isn't exploiting the GPU effectively.

If I include the CUDAEx() call, the error message is

GPU compilation of kernel transduce_kernel!(Nothing, Transducers.Reduction{Transducers.Map{typeof(first)}, Transducers.BottomRF{Transducers.AdHocRF{var"#__##oninit_function#434#22", typeof(identity), InitialValues.AdjoinIdentity{var"#__##reducing_function#435#23"{Matrix{Float64}, CuDeviceMatrix{Float64, 1}}}, typeof(identity), var"#__##combine_function#436#24"}}}, Transducers.InitOf{Transducers.DefaultInitOf}, Int64, Base.OneTo{Int64}, UnitRange{Int64}) failed
KernelError: passing and using non-bitstype argument

Argument 3 to your kernel function is of type Transducers.Reduction{Transducers.Map{typeof(first)}, Transducers.BottomRF{Transducers.AdHocRF{var"#__##oninit_function#434#22", typeof(identity), InitialValues.AdjoinIdentity{var"#__##reducing_function#435#23"{Matrix{Float64}, CuDeviceMatrix{Float64, 1}}}, typeof(identity), var"#__##combine_function#436#24"}}}, which is not isbits:
  .inner is of type Transducers.BottomRF{Transducers.AdHocRF{var"#__##oninit_function#434#22", typeof(identity), InitialValues.AdjoinIdentity{var"#__##reducing_function#435#23"{Matrix{Float64}, CuDeviceMatrix{Float64, 1}}}, typeof(identity), var"#__##combine_function#436#24"}} which is not isbits.
    .inner is of type Transducers.AdHocRF{var"#__##oninit_function#434#22", typeof(identity), InitialValues.AdjoinIdentity{var"#__##reducing_function#435#23"{Matrix{Float64}, CuDeviceMatrix{Float64, 1}}}, typeof(identity), var"#__##combine_function#436#24"} which is not isbits.
      .next is of type InitialValues.AdjoinIdentity{var"#__##reducing_function#435#23"{Matrix{Float64}, CuDeviceMatrix{Float64, 1}}} which is not isbits.
        .op is of type var"#__##reducing_function#435#23"{Matrix{Float64}, CuDeviceMatrix{Float64, 1}} which is not isbits.
          .f is of type Matrix{Float64} which is not isbits.

Histogramming with vegas

i would like to build a histogram of some distribution while integrating it with vegas. Is there a way to access the "vegas-weight" of a point at each step in the integration. (with "vegas-weight" I mean the weight of the point determined by the importance function of vegas). If I am not mistaking the Fortran version of the code actually passes this weight to the integrand function.

`nvec` parameter not implemented?

I was trying to use the nvec parameter. (I think nvec could be used to implement parallelism, this is the way the Mathematica interface to cubalib implements parallelism.) However, I couldn't figure out how to do it. Looking at the code, it looks like the nvec parameter isn't passed to the integrand function:
https://github.com/giordano/Cuba.jl/blob/master/src/Cuba.jl#L86
and
https://github.com/giordano/Cuba.jl/blob/master/src/Cuba.jl#L96

Support Windows

Cuba.jl hasn't been tested on Windows and probably current build script won't work on that system. Also an AppVeyor recipe to test the package would be useful.

No Interface to Divonne Integrand Phase Argument

On page 14 of the Cuba paper (https://arxiv.org/abs/hep-ph/0404043) there is discussion about having a phase argument on the integrand so that an approximation with similar peak structure can be passed for the beginning of the algorithm. There doesn't seem to be support on the Julia interface for accessing this functionality

change license to MIT?

I was looking for LGPL packages, suspecting that there should not be many since for a Julia package there's no such thing as "linking" making LGPL equivalent to GPL. I found this package, which appears to be LGPL because it's a wrapper around an LGPL library (which is dynamically linked, so that makes sense). Wrapper code need not have the same license as the wrapped library, so this package could and arguably should be MIT licensed. As it stands, whereas Cuba can be used with/by non-open source software so long as modifications to Cuba are released in accordance with the LGPL, Cuba.jl cannot be used with/by non-open source software, which I don't think was the intention. Since there are still relatively few contributors, getting permission for a license change should be tractable.

Crash at maxevals=100

This gives a segmentation fault error:

using Cuba 
vegas((x, f) -> f[1] = cos(x[1]), maxevals=100)

It works fine at 1000, and also this works fine:

vegas((x, f) -> f[1] = cos(x[1]), nstart=99, nincrease=99, maxevals=100)

I get a similar error from other functions (suave, divonne, cuhre) at maxevals=100.

Cuba 0.4.0, Julia 0.6.0, on a mac.

Cannot differentiate through `cuhre`. Explicit type inference

I love the library and I don't any other better way in Julia to do >1 dim integrals.

Just discovered an issed trying calculate gradient of the integrated function

using Cuba
using ForwardDiff
#
b(cosθ,ϕ,c) = c[1]*cosθ*sin(ϕ) + c[2]*cosθ^2*cos(ϕ)^2
f(c) = cuhre((x,f)->f[1]=b(2x[1]-1,π*(2x[2]-1),c),2,1)
ForwardDiff.gradient(c->f(c).integral[1], [1,1.1])

giving the error

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(χ2),Float64},Float64,3})
Closest candidates are:
  Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:718
  Float64(!Matched::Int8) at float.jl:60
  ...

most likely at dointegral function

Cuba.jl/src/Cuba.jl

Line 181 in 3681ed8

@inline function dointegrate(x::T) where {T<:Integrand}

Interestingly a gradient of 1dim integral taken with QuadGK works (well, native Julia).

Support parallelization

I'm not sure it's actually feasible, but it would be great if Cuba.jl could take advantage of parallelization capability of Cuba Library. Concurrency is achieved using fork and wait, but trying to increase the number of Cuba cores in Cuba.jl, after having increased the number of Julia processes with addprocs (without this, Cuba will spawn useless Julia processes), results in an undefined reference to fork function.

Segmentation fault in `vegas` and `suave` with `maxevals` smaller than 256

I got segmentation error when trying to pass maxevals to vegas.
Minimal example: vegas((x, f) -> f[1] = cos(x[1]); atol=1e-4, maxevals=20)

All other kwargs seem to work fine.
Machine info: Pop!_OS 22.04 LTS and Apple M2 (Sonoma).
Cuba version: v2.3.0

Lots of allocations

Cuba produces a lot of allocations, 4 per iteration in the example below. Is there a way to avoid this?

julia> using Cuba

julia> f!(x,ret ) = ret
f! (generic function with 1 method)

julia> vegas(f!); # warmup

julia> @time vegas(f!)
  0.001090 seconds (4.01 k allocations: 187.984 KiB)
Component:
 1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000

Integration on general domains

I am wondering if it possible to integrate an n-dimensional function on a set that is not an n-dimensional box in a way similar to the Monte Carlo integration routines in Numerical Recipes. If the function returns 0 when x is not in the set maybe this will just work in the current implementation of Cuba.jl?

Fail to precompile Cuba

I have used the following code to install the Cuba package, and it works fine.

julia> Pkg.update()
julia> Pkg.add("Cuba")

However, i got an error 'ERROR: Failed to precompile Cuba', when I use the following code:

using Cuba

What is the solution for this please?

I have a windows 10 system, and just installed Julia 1.0.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

maxevals produces InexactError()

Simple increase in value for maxevals produces error. You can replicate problem with this simple code
(Debian Linux, julia 0.5.1, 64 bit kernel)

In[]:= vegas((x,f)->f[1]=cos(0.1*x[1]), maxevals = 1e9)
Out[]:= ([0.998335],[4.70961e-5],[-999.0],1000,0,0)

In[]:= vegas((x,f)->f[1]=cos(0.1*x[1]), maxevals = 3e9)
Out[]:= InexactError()
in dointegrate(::Symbol, ::Ptr{Void}, ::Int64, ::Int64, ::Ptr{Void}, ::Int64, ::Float64, ::Float64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Float64, ::Int64, ::Int64, ::Int64, ::Int64, ::Float64, ::Float64, ::Float64, ::Int64, ::Int64, ::Int64, ::Int64, ::Ptr{Void}, ::Int64, ::String, ::Ptr{Void}) at /home/sageusr/.julia/v0.5/Cuba/src/Cuba.jl:169
in #vegas#1(::Int64, ::Float64, ::Float64, ::Int64, ::Int64, ::Int64, ::Float64, ::Int64, ::Int64, ::Int64, ::Int64, ::String, ::Ptr{Void}, ::Cuba.#vegas, ::##33#34, ::Int64, ::Int64) at /home/sageusr/.julia/v0.5/Cuba/src/Cuba.jl:295
in (::Cuba.#kw##vegas)(::Array{Any,1}, ::Cuba.#vegas, ::Function, ::Int64, ::Int64) at ./:0 (repeats 2 times)

cuhre fails at high precisions (?)

The question mark in the title is because it seems to fail, but really I have no idea what the correct result is.

The attached code compares Cuba's cuhre with Cubature's hcubature. The results differ by about 1e-3, despite an abstol of 1e-6 and a reltol of 0, and no limit in maxevals. Both integrators report success. I'm filing this here (rather than at Cubature) because Cubature's results are more in line with what I expect, but I can't exclude an error in hcubature. Both integrators stick to their guns when increasing the precision. The integrand is an indicator function of a polygon

using Interpolations
using Cubature
using Cuba
using PyPlot

const xmax = 1.0
const reltol = 0.0
const abstol = 1e-6

const P1 = BSpline(Linear())

ε(x,y) = 3cos(2pi*x)*cos(2pi*y)+sin(4pi*x)*cos(4pi*y)
function test(p,N)
    x = 0:1/N:xmax-xmax/N
    y = x

    z = ε.(x,y')

    itp = interpolate(z, p, OnGrid())
    εp = scale(itp,x,y)

    function f(x,y,εF)
        εqxy = εp[x,y]
        εqxy < εF ? 1.0 : 0.0
    end

    function error(εF)
        res = cuhre(((x,out)->(out[1] = f(x[1],x[2],εF))), 2,1, reltol=reltol, abstol=abstol,maxevals=1e10)
        println()
        println(res[1][1])
        # println(res)
        res = hcubature(x->f(x[1],x[2],εF), (0,0), (xmax,xmax), reltol=reltol, abstol=abstol, maxevals=0)
        println(res[1])
    end

    error(1.95)
end

test(P1,10)

Cuhre doesn't work for ncomp > 1024

Hi Giordano,

I found that when the number of components of f is greater than 1024, cuhre fails to give the correct answer. Is this the limitation of the algorithm? Or is there some keyword to tune the maximum number of components? Thanks.

Jin