add support for userdata and the corresponding tests (#39)

* add support for userdata and the corresponding tests * make userdata type concrete * working on index.html * add docs for userdata keyword in index.md * add back doctests * change version X.X.X to 2.3.0 * add tests for vectorization with userdata Co-authored-by: Kun Chen <[email protected]> Co-authored-by: Mosè Giordano <[email protected]>
giordano · Sep 22, 2022 · c58f22a · c58f22a
1 parent 01a5f4c
commit c58f22a
Show file tree

Hide file tree

Showing 8 changed files with 220 additions and 75 deletions.
diff --git a/.gitignore b/.gitignore
@@ -6,3 +6,6 @@
 /benchmarks/benchmark-c
 Manifest.toml
 /docs/build/
+
+*.history
+*.vscode
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -25,23 +25,23 @@ All algorithms provided by Cuba library are supported in `Cuba.jl`:
 * [Vegas](https://en.wikipedia.org/wiki/VEGAS_algorithm):
 
   | Basic integration method                                                                | Type                                                                 | [Variance reduction](https://en.wikipedia.org/wiki/Variance_reduction)   |
-  |-----------------------------------------------------------------------------------------|----------------------------------------------------------------------|--------------------------------------------------------------------------|
+  | --------------------------------------------------------------------------------------- | -------------------------------------------------------------------- | ------------------------------------------------------------------------ |
   | [Sobol quasi-random sample](https://en.wikipedia.org/wiki/Sobol_sequence)               | [Monte Carlo](https://en.wikipedia.org/wiki/Monte_Carlo_integration) | [importance sampling](https://en.wikipedia.org/wiki/Importance_sampling) |
   | [Mersenne Twister pseudo-random sample](https://en.wikipedia.org/wiki/Mersenne_Twister) | "                                                                    | "                                                                        |
   | [Ranlux pseudo-random sample](http://arxiv.org/abs/hep-lat/9309020)                     | "                                                                    | "                                                                        |
 
 * Suave
 
   | Basic integration method              | Type        | Variance reduction                                                                                         |
-  |---------------------------------------|-------------|------------------------------------------------------------------------------------------------------------|
+  | ------------------------------------- | ----------- | ---------------------------------------------------------------------------------------------------------- |
   | Sobol quasi-random sample             | Monte Carlo | globally [adaptive subdivision](https://en.wikipedia.org/wiki/Adaptive_quadrature) and importance sampling |
   | Mersenne Twister pseudo-random sample | "           | "                                                                                                          |
   | Ranlux pseudo-random sample           | "           | "                                                                                                          |
 
 * Divonne
 
   | Basic integration method              | Type          | Variance reduction                                                                                                    |
-  |---------------------------------------|---------------|-----------------------------------------------------------------------------------------------------------------------|
+  | ------------------------------------- | ------------- | --------------------------------------------------------------------------------------------------------------------- |
   | Korobov quasi-random sample           | Monte Carlo   | [stratified sampling](https://en.wikipedia.org/wiki/Stratified_sampling) aided by methods from numerical optimization |
   | Sobol quasi-random sample             | "             | "                                                                                                                     |
   | Mersenne Twister pseudo-random sample | "             | "                                                                                                                     |
@@ -51,7 +51,7 @@ All algorithms provided by Cuba library are supported in `Cuba.jl`:
 * Cuhre
 
   | Basic integration method | Type          | Variance reduction            |
-  |--------------------------|---------------|-------------------------------|
+  | ------------------------ | ------------- | ----------------------------- |
   | cubature rules           | deterministic | globally adaptive subdivision |
 
 For more details on the algorithms see the manual included in Cuba
@@ -341,7 +341,7 @@ These are optional keywords common to all functions:
   follows:
 
   | `seed`   | `level` (bits 8--31 of `flags`) | Generator                        |
-  |----------|---------------------------------|----------------------------------|
+  | -------- | ------------------------------- | -------------------------------- |
   | zero     | N/A                             | Sobol (quasi-random)             |
   | non-zero | zero                            | Mersenne Twister (pseudo-random) |
   | non-zero | non-zero                        | Ranlux (pseudo-random)           |
@@ -394,6 +394,13 @@ These are optional keywords common to all functions:
     support parallelization, so this keyword should not have a value
     different from `C_NULL`.
 
+-   `userdata` (arbitrary julia type, default: `missing`): user data 
+    passed to the integrand. See [Passing data to the integrand function](@ref) for a usage example.
+
+!!! note "Compatibility"
+
+	The keyword `userdata` is only supported in `Cuba.jl` after the version 2.3.0.
+
 #### Vegas-Specific Keywords
 
 These optional keywords can be passed only to [`vegas`](@ref):
@@ -903,9 +910,7 @@ Exact result:   1.0 + 1.0im
 
 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.
+example, when the integrand function depends on changing parameters.
 
 For example, the [cumulative distribution
 function](https://en.wikipedia.org/wiki/Cumulative_distribution_function)
@@ -915,23 +920,51 @@ defined by
 
 ```math
 F(x; k) = \int_{0}^{x} \frac{t^{k/2 - 1}\exp(-t/2)}{2^{k/2}\Gamma(k/2)}
+\,\mathrm{d}t = \frac{x}{2^{k/2}\Gamma(k/2)} \int_{0}^{1} (xt)^{k/2 - 1}\exp(-xt/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](https://docs.julialang.org/en/v1/manual/variables-and-scoping/).
-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
+The integrand depends on user-defined parameters ``x`` and ``k``. One option is
+passing a tuple `(x, k)` to the integrand using the `userdata` keyword argument
+in [`vegas`](@ref), [`suave`](@ref), [`divonne`](@ref) or [`cuhre`](@ref).
+The following julia script uses this trick to compute ``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](https://github.com/jiahao/GSL.jl) package.
 
 ```julia
 julia> using Cuba, GSL, Printf, SpecialFunctions
 
+julia> function integrand(t, f, userdata)
+           # Chi-squared probability density function, without constant denominator.
+           # The result of integration will be divided by that factor.
+           # userdata is a tuple (x, k/2), see below
+           x, k = userdata
+           f[1] = (t[1]*x)^(k/2 - 1.0)*exp(-(t[1]*x)/2)
+       end
+
+julia> chi2cdf(x::Real, k::Real) = x*cuhre(integrand; userdata = (x, k))[1][1]/(2^(k/2)*gamma(k/2))
+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
+```
+
+An alternative solution to pass the user-defined data is implementing the integrand 
+as a nested inner function. The inner function can access any variable visible
+in its [scope](https://docs.julialang.org/en/v1/manual/variables-and-scoping/).
+
+```julia
+julia> using Cuba, GSL, Printf, SpecialFunctions
+
 julia> function chi2cdf(x::Real, k::Real)
            k2 = k/2
            # Chi-squared probability density function, without constant denominator.

diff --git a/src/Cuba.jl b/src/Cuba.jl
@@ -72,6 +72,15 @@ function generic_integrand!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
     func!(x, f)
     return Cint(0)
 end
+function generic_integrand_userdata!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
+                            f_::Ptr{Cdouble}, func!_and_userdata)
+    func!, userdata = func!_and_userdata
+    # Get arrays from "x_" and "f_" pointers.
+    x = unsafe_wrap(Array, x_, (ndim,))
+    f = unsafe_wrap(Array, f_, (ncomp,))
+    func!(x, f, userdata)
+    return Cint(0)
+end
 function generic_integrand!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
                             f_::Ptr{Cdouble}, func!, nvec::Cint)
     # Get arrays from "x_" and "f_" pointers.
@@ -80,6 +89,15 @@ function generic_integrand!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
     func!(x, f)
     return Cint(0)
 end
+function generic_integrand_userdata!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
+                            f_::Ptr{Cdouble}, func!_and_userdata, nvec::Cint)
+    func!, userdata = func!_and_userdata
+    # Get arrays from "x_" and "f_" pointers.
+    x = unsafe_wrap(Array, x_, (ndim, nvec))
+    f = unsafe_wrap(Array, f_, (ncomp, nvec))
+    func!(x, f, userdata)
+    return Cint(0)
+end
 
 # Return pointer for "integrand", to be passed as "integrand" argument to Cuba functions.
 integrand_ptr(integrand::T)  where {T} = @cfunction(generic_integrand!, Cint,
@@ -88,6 +106,14 @@ integrand_ptr(integrand::T)  where {T} = @cfunction(generic_integrand!, Cint,
                                                      Ref{Cint}, # ncomp
                                                      Ptr{Cdouble}, # f
                                                      Ref{T})) # userdata
+
+integrand_ptr_userdata(integrand::T, userdata::D)  where {T, D} = @cfunction(generic_integrand_userdata!, Cint,
+                                                    (Ref{Cint}, # ndim
+                                                     Ptr{Cdouble}, # x
+                                                     Ref{Cint}, # ncomp
+                                                     Ptr{Cdouble}, # f
+                                                     Ref{Tuple{T, D}})) # userdata
+
 integrand_ptr_nvec(integrand::T) where {T} = @cfunction(generic_integrand!, Cint,
                                                         (Ref{Cint}, # ndim
                                                          Ptr{Cdouble}, # x
@@ -96,6 +122,14 @@ integrand_ptr_nvec(integrand::T) where {T} = @cfunction(generic_integrand!, Cint
                                                          Ref{T}, # userdata
                                                          Ref{Cint})) # nvec
 
+integrand_ptr_nvec_userdata(integrand::T, userdata::D) where {T, D} = @cfunction(generic_integrand_userdata!, Cint,
+                                                        (Ref{Cint}, # ndim
+                                                         Ptr{Cdouble}, # x
+                                                         Ref{Cint}, # ncomp
+                                                         Ptr{Cdouble}, # f
+                                                         Ref{Tuple{T, D}}, # userdata
+                                                         Ref{Cint})) # nvec
+
 abstract type Integrand{T} end
 
 function __init__()
@@ -162,9 +196,9 @@ end
     # Choose the integrand function wrapper based on the value of `nvec`.  This function is
     # called only once, so the overhead of the following if should be negligible.
     if x.nvec == 1
-        integrand = integrand_ptr(x.func)
+        integrand = ismissing(x.userdata) ? integrand_ptr(x.func) : integrand_ptr_userdata(x.func, x.userdata)
     else
-        integrand = integrand_ptr_nvec(x.func)
+        integrand = ismissing(x.userdata) ? integrand_ptr_nvec(x.func) : integrand_ptr_nvec_userdata(x.func, x.userdata)
     end
     integral  = Vector{Cdouble}(undef, x.ncomp)
     error     = Vector{Cdouble}(undef, x.ncomp)

diff --git a/src/cuhre.jl b/src/cuhre.jl
@@ -1,7 +1,8 @@
 ### cuhre.jl --- Integrate with Cuhre method
 
-struct Cuhre{T} <: Integrand{T}
+struct Cuhre{T, D} <: Integrand{T}
     func::T
+    userdata::D
     ndim::Int
     ncomp::Int
     nvec::Int64
@@ -15,8 +16,11 @@ struct Cuhre{T} <: Integrand{T}
     spin::Ptr{Cvoid}
 end
 
-@inline function dointegrate!(x::Cuhre{T}, integrand, integral,
-                              error, prob, neval, fail, nregions) where {T}
+@inline function dointegrate!(x::Cuhre{T, D}, integrand, integral,
+                              error, prob, neval, fail, nregions) where {T, D}
+
+    userdata = ismissing(x.userdata) ? x.func : (x.func, x.userdata)
+
     ccall((:llCuhre, libcuba), Cdouble,
           (Cint, # ndim
            Cint, # ncomp
@@ -38,7 +42,7 @@ end
            Ptr{Cdouble}, # error
            Ptr{Cdouble}),# prob
           # Input
-          x.ndim, x.ncomp, integrand, x.func, x.nvec, x.rtol, x.atol,
+          x.ndim, x.ncomp, integrand, userdata, x.nvec, x.rtol, x.atol,
           x.flags, x.minevals, x.maxevals, x.key, x.statefile, x.spin,
           # Output
           nregions, neval, fail, integral, error, prob)
@@ -53,6 +57,7 @@ components.  `ncomp` defaults to 1, `ndim` defaults to 2 and must be ≥ 2.
 
 Accepted keywords:
 
+* `userdata`
 * `nvec`
 * `rtol`
 * `atol`
@@ -68,13 +73,13 @@ function cuhre(integrand::T, ndim::Integer=2, ncomp::Integer=1;
                flags::Integer=FLAGS, minevals::Real=MINEVALS,
                maxevals::Real=MAXEVALS, key::Integer=KEY,
                statefile::AbstractString=STATEFILE, spin::Ptr{Cvoid}=SPIN,
-               abstol=missing, reltol=missing) where {T}
+               abstol=missing, reltol=missing, userdata=missing) where {T}
     atol_,rtol_ = tols(atol,rtol,abstol,reltol)
     # Cuhre requires "ndim" to be at least 2, even for an integral over a one
     # dimensional domain.  Instead, we don't prevent users from setting wrong
     # "ndim" values like 0 or negative ones.
     ndim >=2 || throw(ArgumentError("In Cuhre ndim must be ≥ 2"))
-    return dointegrate(Cuhre(integrand, ndim, ncomp, Int64(nvec), Cdouble(rtol_),
+    return dointegrate(Cuhre(integrand, userdata, ndim, ncomp, Int64(nvec), Cdouble(rtol_),
                              Cdouble(atol_), flags, trunc(Int64, minevals),
                              trunc(Int64, maxevals), key, String(statefile), spin))
 end
diff --git a/src/divonne.jl b/src/divonne.jl
@@ -1,7 +1,8 @@
 ### divonne.jl --- Integrate with Divonne method
 
-struct Divonne{T} <: Integrand{T}
+struct Divonne{T, D} <: Integrand{T}
     func::T
+    userdata::D
     ndim::Int
     ncomp::Int
     nvec::Int64
@@ -27,8 +28,11 @@ struct Divonne{T} <: Integrand{T}
     spin::Ptr{Cvoid}
 end
 
-@inline function dointegrate!(x::Divonne{T}, integrand, integral,
-                              error, prob, neval, fail, nregions) where {T}
+@inline function dointegrate!(x::Divonne{T, D}, integrand, integral,
+                              error, prob, neval, fail, nregions) where {T, D}
+
+    userdata = ismissing(x.userdata) ? x.func : (x.func, x.userdata)
+
     ccall((:llDivonne, libcuba), Cdouble,
           (Cint, # ndim
            Cint, # ncomp
@@ -62,7 +66,7 @@ end
            Ptr{Cdouble}, # error
            Ptr{Cdouble}),# prob
           # Input
-          x.ndim, x.ncomp, integrand, x.func, x.nvec, x.rtol,
+          x.ndim, x.ncomp, integrand, userdata, x.nvec, x.rtol,
           x.atol, x.flags, x.seed, x.minevals, x.maxevals, x.key1, x.key2,
           x.key3, x.maxpass, x.border, x.maxchisq, x.mindeviation, x.ngiven,
           x.ldxgiven, x.xgiven, x.nextra, x.peakfinder, x.statefile, x.spin,
@@ -79,6 +83,7 @@ components. `ncomp` defaults to 1, `ndim` defaults to 2 and must be ≥ 2.
 
 Accepted keywords:
 
+* `userdata`
 * `nvec`
 * `rtol`
 * `atol`
@@ -116,13 +121,13 @@ function divonne(integrand::T, ndim::Integer=2, ncomp::Integer=1;
                  nextra::Integer=NEXTRA,
                  peakfinder::Ptr{Cvoid}=PEAKFINDER,
                  statefile::AbstractString=STATEFILE,
-                 spin::Ptr{Cvoid}=SPIN, reltol=missing, abstol=missing) where {T}
+                 spin::Ptr{Cvoid}=SPIN, reltol=missing, abstol=missing, userdata=missing) where {T}
     atol_,rtol_ = tols(atol,rtol,abstol,reltol)
     # Divonne requires "ndim" to be at least 2, even for an integral over a one
     # dimensional domain.  Instead, we don't prevent users from setting wrong
     # "ndim" values like 0 or negative ones.
     ndim >=2 || throw(ArgumentError("In Divonne ndim must be ≥ 2"))
-    return dointegrate(Divonne(integrand, ndim, ncomp, Int64(nvec), Cdouble(rtol_),
+    return dointegrate(Divonne(integrand, userdata, ndim, ncomp, Int64(nvec), Cdouble(rtol_),
                                Cdouble(atol_), flags, seed, trunc(Int64, minevals),
                                trunc(Int64, maxevals), key1, key2, key3, maxpass,
                                Cdouble(border), Cdouble(maxchisq), Cdouble(mindeviation),

diff --git a/src/suave.jl b/src/suave.jl
@@ -1,7 +1,8 @@
 ### suave.jl --- Integrate with Suave method
 
-struct Suave{T} <: Integrand{T}
+struct Suave{T, D} <: Integrand{T}
     func::T
+    userdata::D
     ndim::Int
     ncomp::Int
     nvec::Int64
@@ -18,8 +19,11 @@ struct Suave{T} <: Integrand{T}
     spin::Ptr{Cvoid}
 end
 
-@inline function dointegrate!(x::Suave{T}, integrand, integral,
-                              error, prob, neval, fail, nregions) where {T}
+@inline function dointegrate!(x::Suave{T, D}, integrand, integral,
+                              error, prob, neval, fail, nregions) where {T, D}
+
+    userdata = ismissing(x.userdata) ? x.func : (x.func, x.userdata)
+
     ccall((:llSuave, libcuba), Cdouble,
           (Cint, # ndim
            Cint, # ncomp
@@ -44,7 +48,7 @@ end
            Ptr{Cdouble}, # error
            Ptr{Cdouble}),# prob
           # Input
-          x.ndim, x.ncomp, integrand, x.func, x.nvec,
+          x.ndim, x.ncomp, integrand, userdata, x.nvec,
           x.rtol, x.atol, x.flags, x.seed, x.minevals, x.maxevals,
           x.nnew, x.nmin, x.flatness, x.statefile, x.spin,
           # Output
@@ -60,6 +64,7 @@ components. `ndim` and `ncomp` default to 1.
 
 Accepted keywords:
 
+* `userdata`
 * `nvec`
 * `rtol`
 * `atol`
@@ -79,9 +84,9 @@ function suave(integrand::T, ndim::Integer=1, ncomp::Integer=1;
                minevals::Real=MINEVALS, maxevals::Real=MAXEVALS,
                nnew::Integer=NNEW, nmin::Integer=NMIN, flatness::Real=FLATNESS,
                statefile::AbstractString=STATEFILE, spin::Ptr{Cvoid}=SPIN,
-               reltol=missing, abstol=missing) where {T}
+               reltol=missing, abstol=missing, userdata=missing) where {T}
     atol_,rtol_ = tols(atol,rtol,abstol,reltol)
-    return dointegrate(Suave(integrand, ndim, ncomp, Int64(nvec), Cdouble(rtol_),
+    return dointegrate(Suave(integrand, userdata, ndim, ncomp, Int64(nvec), Cdouble(rtol_),
                              Cdouble(atol_), flags, seed, trunc(Int64, minevals),
                              trunc(Int64, maxevals), Int64(nnew), Int64(nmin),
                              Cdouble(flatness), String(statefile), spin))