From 9d55b94980eedad367b67e07a253337df16edc20 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 19 Nov 2023 22:24:35 -0500 Subject: [PATCH] wrap mcintegration --- Project.toml | 2 ++ docs/src/solvers/IntegralSolvers.md | 5 ++- src/Integrals.jl | 48 +++++++++++++++++++++++++++-- src/algorithms.jl | 18 +++++++++++ test/interface_tests.jl | 3 +- 5 files changed, 72 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 1b901a7a..6d819a67 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "4.1.0" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167" MonteCarloIntegration = "4886b29c-78c9-11e9-0a6e-41e1f4161f7b" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -41,6 +42,7 @@ FastGaussQuadrature = "0.5" ForwardDiff = "0.10" HCubature = "1.4" LinearAlgebra = "1.9" +MCIntegration = "0.4" MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3, 0.1" QuadGK = "2.5" Reexport = "0.2, 1.0" diff --git a/docs/src/solvers/IntegralSolvers.md b/docs/src/solvers/IntegralSolvers.md index 3997566c..241d4d29 100644 --- a/docs/src/solvers/IntegralSolvers.md +++ b/docs/src/solvers/IntegralSolvers.md @@ -4,7 +4,9 @@ The following algorithms are available: - `QuadGKJL`: Uses QuadGK.jl. Requires `nout=1` and `batch=0`, in-place is not allowed. - `HCubatureJL`: Uses HCubature.jl. Requires `batch=0`. - - `VEGAS`: Uses MonteCarloIntegration.jl. Requires `nout=1`. Works only for `>1`-dimensional integrations. + - `VEGAS`: Uses MonteCarloIntegration.jl. Requires `nout=1`. Works only for + `>1`-dimensional integrations. + - `VEGASMC`: Uses MCIntegration.jl. Requires `nout=1`. Works only for `>1`-dimensional integrations. - `CubatureJLh`: h-Cubature from Cubature.jl. Requires `using Cubature`. - `CubatureJLp`: p-Cubature from Cubature.jl. Requires `using Cubature`. - `CubaVegas`: Vegas from Cuba.jl. Requires `using Cuba`, `nout=1`. @@ -20,6 +22,7 @@ The following algorithms are available: QuadGKJL HCubatureJL VEGAS +VEGASMC CubaVegas CubaSUAVE CubaDivonne diff --git a/src/Integrals.jl b/src/Integrals.jl index e6e3de74..7b50cf8c 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -4,7 +4,7 @@ if !isdefined(Base, :get_extension) using Requires end -using Reexport, MonteCarloIntegration, QuadGK, HCubature +using Reexport, MonteCarloIntegration, QuadGK, HCubature, MCIntegration @reexport using SciMLBase using LinearAlgebra @@ -174,7 +174,51 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success) end -export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, QuadratureRule, TrapezoidalRule + +function __solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, domain, p; + reltol = nothing, abstol = nothing, maxiters = 1000) + lb, ub = domain + mid = collect((lb + ub) / 2) + vars = Continuous(vec([tuple(a,b) for (a,b) in zip(lb, ub)])) + + if prob.f isa BatchIntegralFunction + error("VEGASMC doesn't support batching. See https://github.com/numericalEFT/MCIntegration.jl/issues/29") + else + if isinplace(prob) + f0 = similar(prob.f.integrand_prototype) + f_ = (x, f, c) -> begin + n = 0 + for v in x + mid[n+=1] = first(v) + end + prob.f(f0, mid, p) + f .= vec(f0) + end + else + f0 = prob.f(mid, p) + f_ = (x, c) -> begin + n = 0 + for v in x + mid[n+=1] = first(v) + end + prob.f(mid, p) + end + end + dof = f0 isa Number ? 1 : ones(Int, length(f0)) + res = integrate(f_, inplace=isinplace(prob), var=vars, dof=dof, solver=:vegasmc, + neval=alg.neval, niter=min(maxiters,alg.niter), block=alg.block, adapt=alg.adapt, + gamma=alg.gamma, verbose=alg.verbose, debug=alg.debug) + out, err, chi = if f0 isa Number + map(only, (res.mean, res.stdev, res.chi2)) + else + map(a -> reshape(a, size(f0)), (res.mean, res.stdev, res.chi2)) + end + SciMLBase.build_solution(prob, VEGASMC(), out, err, chi=chi, retcode = ReturnCode.Success) + end +end + + +export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, TrapezoidalRule export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre export CubatureJLh, CubatureJLp export ArblibJL diff --git a/src/algorithms.jl b/src/algorithms.jl index d40dceaa..c294fff6 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -87,6 +87,24 @@ struct VEGAS <: SciMLBase.AbstractIntegralAlgorithm end VEGAS(; nbins = 100, ncalls = 1000, debug = false) = VEGAS(nbins, ncalls, debug) + +""" + VEGASMC() + +Markov-chain based Vegas algorithm from MCIntegration.jl +""" +struct VEGASMC <: SciMLBase.AbstractIntegralAlgorithm + neval::Int + niter::Int + block::Int + adapt::Bool + gamma::Float64 + verbose::Int + debug::Bool +end +VEGASMC(; neval=10^4, niter=20, block=16, adapt=true, gamma=1.0, verbose=-2, debug=false) = + VEGASMC(neval, niter, block, adapt, gamma, verbose, debug) + """ GaussLegendre{C, N, W} diff --git a/test/interface_tests.jl b/test/interface_tests.jl index 2b61abd9..f0d7eaff 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -7,7 +7,7 @@ max_nout_test = 2 reltol = 1e-3 abstol = 1e-3 -algs = [QuadGKJL, HCubatureJL, CubatureJLh, CubatureJLp, VEGAS, #CubaVegas, +algs = [QuadGKJL, HCubatureJL, CubatureJLh, CubatureJLp, VEGAS, VEGASMC, #CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre, ArblibJL] alg_req = Dict(QuadGKJL => (nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, @@ -16,6 +16,7 @@ alg_req = Dict(QuadGKJL => (nout = 1, allows_batch = false, min_dim = 1, max_dim max_dim = Inf, allows_iip = true), VEGAS => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, allows_iip = true), + VEGASMC => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = Inf, allows_iip = true), CubatureJLh => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, allows_iip = true), CubatureJLp => (nout = Inf, allows_batch = true, min_dim = 1,