wrap mcintegration

Author: lxvm

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"

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

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

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}
 

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,