From c58f22ab5a1da4e06eab9856238a59240e01ab6c Mon Sep 17 00:00:00 2001 From: Kun Chen Date: Thu, 22 Sep 2022 12:06:44 -0400 Subject: [PATCH] add support for userdata and the corresponding tests (#39) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 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 Co-authored-by: Mosè Giordano --- .gitignore | 3 ++ docs/src/index.md | 67 +++++++++++++++++++------- src/Cuba.jl | 38 ++++++++++++++- src/cuhre.jl | 17 ++++--- src/divonne.jl | 17 ++++--- src/suave.jl | 17 ++++--- src/vegas.jl | 17 ++++--- test/runtests.jl | 119 +++++++++++++++++++++++++++++++++------------- 8 files changed, 220 insertions(+), 75 deletions(-) diff --git a/.gitignore b/.gitignore index 58f509a..7c2cd9a 100644 --- 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 index 55ec4fb..ef66c0e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -25,7 +25,7 @@ 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) | " | " | @@ -33,7 +33,7 @@ All algorithms provided by Cuba library are supported in `Cuba.jl`: * 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 | " | " | @@ -41,7 +41,7 @@ All algorithms provided by Cuba library are supported in `Cuba.jl`: * 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 index d186ddf..0158435 100644 --- 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 index 054e56c..57b8746 100644 --- 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 index 4fb7b98..0146b6a 100644 --- 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 index 7353a18..38590ed 100644 --- 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)) diff --git a/src/vegas.jl b/src/vegas.jl index 4aebbef..ea3c946 100644 --- a/src/vegas.jl +++ b/src/vegas.jl @@ -1,7 +1,8 @@ ### vegas.jl --- Integrate with Vegas method -struct Vegas{T} <: Integrand{T} +struct Vegas{T, D} <: Integrand{T} func::T + userdata::D ndim::Int ncomp::Int nvec::Int64 @@ -19,8 +20,11 @@ struct Vegas{T} <: Integrand{T} spin::Ptr{Cvoid} end -@inline function dointegrate!(x::Vegas{T}, integrand, integral, - error, prob, neval, fail, nregions) where {T} +@inline function dointegrate!(x::Vegas{T, D}, integrand, integral, + error, prob, neval, fail, nregions) where {T, D} + + userdata = ismissing(x.userdata) ? x.func : (x.func, x.userdata) + ccall((:llVegas, libcuba), Cdouble, (Cint, # ndim Cint, # ncomp @@ -45,7 +49,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.nstart, x.nincrease, x.nbatch, x.gridno, x.statefile, x.spin, # Output @@ -61,6 +65,7 @@ components. `ndim` and `ncomp` default to 1. Accepted keywords: +* `userdata` * `nvec` * `rtol` * `atol` @@ -82,9 +87,9 @@ function vegas(integrand::T, ndim::Integer=1, ncomp::Integer=1; nstart::Integer=NSTART, nincrease::Integer=NINCREASE, nbatch::Integer=NBATCH, gridno::Integer=GRIDNO, 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(Vegas(integrand, ndim, ncomp, Int64(nvec), Cdouble(rtol_), + return dointegrate(Vegas(integrand, userdata, ndim, ncomp, Int64(nvec), Cdouble(rtol_), Cdouble(atol_), flags, seed, trunc(Int64, minevals), trunc(Int64, maxevals), Int64(nstart), Int64(nincrease), Int64(nbatch), gridno, diff --git a/test/runtests.jl b/test/runtests.jl index 9628d43..b09b145 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,9 +4,9 @@ using Cuba using Test using Distributed, LinearAlgebra -f1(x,y,z) = sin(x)*cos(y)*exp(z) -f2(x,y,z) = exp(-(x*x + y*y + z*z)) -f3(x,y,z) = 1/(1 - x*y*z) +f1(x, y, z) = sin(x) * cos(y) * exp(z) +f2(x, y, z) = exp(-(x * x + y * y + z * z)) +f3(x, y, z) = 1 / (1 - x * y * z) function integrand1(x, f) f[1] = f1(x[1], x[2], x[3]) @@ -17,73 +17,125 @@ end # Make sure using "addprocs" doesn't make the program segfault. addprocs(1) Cuba.cores(0, 10000) -Cuba.accel(0,1000) +Cuba.accel(0, 1000) # Test results and make sure the estimation of error is exact. ncomp = 3 @testset "$alg" for (alg, atol) in ((vegas, 1e-4), (suave, 1e-3), - (divonne, 1e-4), (cuhre, 1e-8)) + (divonne, 1e-4), (cuhre, 1e-8)) # Make sure that using maxevals > typemax(Int32) doesn't result into InexactError. if alg == divonne result = @inferred alg(integrand1, 3, ncomp, atol=atol, rtol=1e-8, - flags=0, border = 1e-5, maxevals = 3000000000) + flags=0, border=1e-5, maxevals=3000000000) else result = @inferred alg(integrand1, 3, ncomp, atol=atol, rtol=1e-8, - flags=0, maxevals = 3e9) + flags=0, maxevals=3e9) end # Analytic expressions: ((e-1)*(1-cos(1))*sin(1), (sqrt(pi)*erf(1)/2)^3, zeta(3)) for (i, answer) in enumerate((0.6646696797813771, 0.41653838588663805, 1.2020569031595951)) - @test result[1][i] ≈ answer atol=result[2][i] + @test result[1][i] ≈ answer atol = result[2][i] + end +end + +struct UserData + para::Float64 +end +function integrand2(x, f, userdata) + f[1] = f1(x[1], x[2], x[3]) + userdata.para + f[2] = f2(x[1], x[2], x[3]) + userdata.para + f[3] = f3(x[1], x[2], x[3]) + userdata.para +end +@testset "$alg with userdata" for (alg, atol) in ((vegas, 1e-4), (suave, 1e-3), + (divonne, 1e-4), (cuhre, 1e-8)) + # Make sure that using maxevals > typemax(Int32) doesn't result into InexactError. + if alg == divonne + result = @inferred alg(integrand2, 3, ncomp, atol=atol, rtol=1e-8, + flags=0, border=1e-5, maxevals=3000000000, userdata=UserData(0.1)) + else + result = @inferred alg(integrand2, 3, ncomp, atol=atol, rtol=1e-8, + flags=0, maxevals=3e9, userdata=UserData(0.1)) + end + # Analytic expressions: ((e-1)*(1-cos(1))*sin(1)+1.0, (sqrt(pi)*erf(1)/2)^3+1.0, zeta(3)+1.0) + for (i, answer) in enumerate((0.7646696797813771, 0.51653838588663805, 1.3020569031595951)) + @test result[1][i] ≈ answer atol = result[2][i] end end @testset "ndim" begin func(x, f) = (f[] = norm(x)) - answer_1d = 1/2 # integral of abs(x) in 1D - answer_2d = (8 * asinh(1) + 2^(7/2))/24 # integral of sqrt(x^2 + y^2) in 2D - @test @inferred(vegas(func))[1][1] ≈ answer_1d rtol = 1e-4 - @test @inferred(suave(func))[1][1] ≈ answer_1d rtol = 1e-2 + answer_1d = 1 / 2 # integral of abs(x) in 1D + answer_2d = (8 * asinh(1) + 2^(7 / 2)) / 24 # integral of sqrt(x^2 + y^2) in 2D + @test @inferred(vegas(func))[1][1] ≈ answer_1d rtol = 1e-4 + @test @inferred(suave(func))[1][1] ≈ answer_1d rtol = 1e-2 @test @inferred(divonne(func))[1][1] ≈ answer_2d rtol = 1e-4 - @test @inferred(cuhre(func))[1][1] ≈ answer_2d rtol = 1e-4 + @test @inferred(cuhre(func))[1][1] ≈ answer_2d rtol = 1e-4 @test_throws ArgumentError cuhre(func, 1) @test_throws ArgumentError divonne(func, 1) end @testset "Integral over infinite domain" begin - func(x) = log(1 + x^2)/(1 + x^2) - result, rest = @inferred cuhre((x, f) -> f[1] = func(x[1]/(1 - x[1]))/(1 - x[1])^2, - atol = 1e-12, rtol = 1e-10) + func(x) = log(1 + x^2) / (1 + x^2) + result, rest = @inferred cuhre((x, f) -> f[1] = func(x[1] / (1 - x[1])) / (1 - x[1])^2, + atol=1e-12, rtol=1e-10) @test result[1] ≈ pi * log(2) atol = 3e-12 end @testset "Vectorization" begin for alg in (vegas, suave, divonne, cuhre) - result1, err1, _ = @inferred alg((x,f) -> f[1] = x[1] + cos(x[2]) - exp(x[3]), 3) - result2, err2, _ = @inferred alg((x,f) -> f[1,:] .= x[1,:] .+ cos.(x[2,:]) .- exp.(x[3,:]), - 3, nvec = 10) + result1, err1, _ = @inferred alg((x, f) -> f[1] = x[1] + cos(x[2]) - exp(x[3]), 3) + result2, err2, _ = @inferred alg((x, f) -> f[1, :] .= x[1, :] .+ cos.(x[2, :]) .- exp.(x[3, :]), + 3, nvec=10) + @test result1 == result2 + @test err1 == err2 + result1, err1, _ = @inferred alg((x, f) -> begin + f[1] = sin(x[1]) + f[2] = sqrt(x[2]) + end, 2, 2) + result2, err2, _ = @inferred alg((x, f) -> begin + f[1, :] .= sin.(x[1, :]) + f[2, :] .= sqrt.(x[2, :]) + end, + 2, 2, nvec=10) + @test result1 == result2 + @test err1 == err2 + end +end + +@testset "Vectorization with userdata" begin + for alg in (vegas, suave, divonne, cuhre) + result1, err1, _ = @inferred alg((x, f, u) -> f[1] = u.para + x[1] + cos(x[2]) - exp(x[3]), 3; userdata=UserData(0.1)) + result2, err2, _ = @inferred alg((x, f, u) -> f[1, :] .= u.para .+ x[1, :] .+ cos.(x[2, :]) .- exp.(x[3, :]), + 3, nvec=10, userdata=UserData(0.1)) @test result1 == result2 - @test err1 == err2 - result1, err1, _ = @inferred alg((x,f) -> begin f[1] = sin(x[1]); f[2] = sqrt(x[2]) end, 2, 2) - result2, err2, _ = @inferred alg((x,f) -> begin f[1,:] .= sin.(x[1,:]); f[2,:] .= sqrt.(x[2,:]) end, - 2, 2, nvec = 10) + @test err1 == err2 + result1, err1, _ = @inferred alg((x, f, u) -> begin + f[1] = sin(x[1]) + u.para + f[2] = sqrt(x[2]) + u.para + end, 2, 2; userdata=UserData(0.1)) + result2, err2, _ = @inferred alg((x, f, u) -> begin + f[1, :] .= sin.(x[1, :]) .+ u.para + f[2, :] .= sqrt.(x[2, :]) .+ u.para + end, + 2, 2, nvec=10, userdata=UserData(0.1)) @test result1 == result2 - @test err1 == err2 + @test err1 == err2 end end + @testset "Show" begin @test occursin("Note: The desired accuracy was reached", - repr(vegas((x, f) -> f[1] = x[1]))) + repr(vegas((x, f) -> f[1] = x[1]))) @test occursin("Note: The accuracy was not met", - repr(suave((x,f) -> f[1] = x[1], atol = 1e-12, rtol = 1e-12))) + repr(suave((x, f) -> f[1] = x[1], atol=1e-12, rtol=1e-12))) @test occursin("Try increasing `maxevals` to", - repr(divonne((x, f) -> f[1] = exp(x[1])*cos(x[1]), - atol = 1e-9, rtol = 1e-9))) + repr(divonne((x, f) -> f[1] = exp(x[1]) * cos(x[1]), + atol=1e-9, rtol=1e-9))) @test occursin("Note: Dimension out of range", - repr(Cuba.dointegrate(Cuba.Cuhre((x, f) -> f[1] = x[1], 1, 1, Int64(Cuba.NVEC), - Cuba.RTOL, Cuba.ATOL, Cuba.FLAGS, - Int64(Cuba.MINEVALS), Int64(Cuba.MAXEVALS), - Cuba.KEY, Cuba.STATEFILE, Cuba.SPIN)))) + repr(Cuba.dointegrate(Cuba.Cuhre((x, f) -> f[1] = x[1], missing, 1, 1, Int64(Cuba.NVEC), + Cuba.RTOL, Cuba.ATOL, Cuba.FLAGS, + Int64(Cuba.MINEVALS), Int64(Cuba.MAXEVALS), + Cuba.KEY, Cuba.STATEFILE, Cuba.SPIN)))) end @@ -93,3 +145,6 @@ Cuba.exit(C_NULL, C_NULL) # Dummy call just to increase code coverage Cuba.integrand_ptr(Cuba.generic_integrand!) +Cuba.integrand_ptr_userdata(Cuba.generic_integrand!, missing) +Cuba.integrand_ptr_nvec(Cuba.generic_integrand!) +Cuba.integrand_ptr_nvec_userdata(Cuba.generic_integrand!, missing)