Migrate to SciMLBase v2

.github/workflows/CI.yml

@@ -19,7 +19,6 @@ jobs:
  - Core
  version:
  - '1'
- - '1.6'
  steps:
  - uses: actions/checkout@v4
  - uses: julia-actions/setup-julia@v1
@@ -39,7 +38,7 @@ jobs:
  - uses: julia-actions/julia-runtest@v1
  - uses: julia-actions/julia-processcoverage@v1
  with:
- directories: src,ext,lib/IntegralsCuba/src,lib/IntegralsCubature/src
+ directories: src,ext
  - uses: codecov/codecov-action@v3
  with:
  files: lcov.info

.github/workflows/CompatHelper.yml

@@ -23,4 +23,4 @@ jobs:
  - name: CompatHelper.main()
  env:
  GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
- run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs", "lib/IntegralsCubature", "lib/IntegralsCuba"])'
+ run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs"])'

.github/workflows/Documentation.yml

@@ -25,7 +25,7 @@ jobs:
  run: julia --project=docs/ --code-coverage=user docs/make.jl
  - uses: julia-actions/julia-processcoverage@v1
  with:
- directories: src,lib/IntegralsCuba/src,lib/IntegralsCubature/src
+ directories: src
  - uses: codecov/codecov-action@v3
  with:
  files: lcov.info

Project.toml

@@ -15,32 +15,41 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
+Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
+IntegralsCubaExt = "Cuba"
+IntegralsCubatureExt = "Cubature"
 IntegralsFastGaussQuadratureExt = "FastGaussQuadrature"
 IntegralsForwardDiffExt = "ForwardDiff"
 IntegralsZygoteExt = ["Zygote", "ChainRulesCore"]
 
 [compat]
 ChainRulesCore = "0.10.7, 1"
 CommonSolve = "0.2"
+Cuba = "2"
+Cubature = "1"
 Distributions = "0.23, 0.24, 0.25"
 FastGaussQuadrature = "0.5"
 ForwardDiff = "0.10"
 HCubature = "1.4"
-MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3"
+LinearAlgebra = "1.9"
+MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3, 0.1"
 QuadGK = "2.5"
 Reexport = "0.2, 1.0"
 Requires = "1"
-SciMLBase = "1.98"
+SciMLBase = "2.6"
 Zygote = "0.4.22, 0.5, 0.6"
-julia = "1.6"
+julia = "1.9"
 
 [extras]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
+Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
@@ -53,4 +62,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature"]
+test = ["SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"]

docs/Project.toml

@@ -1,19 +1,19 @@
 [deps]
+Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
+Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
-IntegralsCuba = "e00cd5f1-6337-4131-8b37-28b2fe4cd6cb"
-IntegralsCubature = "c31f79ba-6e32-46d4-a52f-182a8ac42a54"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+Cuba = "2"
+Cubature = "1"
 Distributions = "0.25"
 Documenter = "1"
 FiniteDiff = "2"
 ForwardDiff = "0.10"
 Integrals = "3"
-IntegralsCuba = "0.3"
-IntegralsCubature = "0.2"
 Zygote = "0.6"

docs/src/solvers/IntegralSolvers.md

@@ -5,12 +5,12 @@ 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.
- - `CubatureJLh`: h-Cubature from Cubature.jl. Requires `using IntegralsCubature`.
- - `CubatureJLp`: p-Cubature from Cubature.jl. Requires `using IntegralsCubature`.
- - `CubaVegas`: Vegas from Cuba.jl. Requires `using IntegralsCuba`, `nout=1`.
- - `CubaSUAVE`: SUAVE from Cuba.jl. Requires `using IntegralsCuba`.
- - `CubaDivonne`: Divonne from Cuba.jl. Requires `using IntegralsCuba`. Works only for `>1`-dimensional integrations.
- - `CubaCuhre`: Cuhre from Cuba.jl. Requires `using IntegralsCuba`. 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`.
+ - `CubaSUAVE`: SUAVE from Cuba.jl. Requires `using Cuba`.
+ - `CubaDivonne`: Divonne from Cuba.jl. Requires `using Cuba`. Works only for `>1`-dimensional integrations.
+ - `CubaCuhre`: Cuhre from Cuba.jl. Requires `using Cuba`. Works only for `>1`-dimensional integrations.
  - `GaussLegendre`: Uses Gauss-Legendre quadrature with nodes and weights from FastGaussQuadrature.jl.
  - `QuadratureRule`: Accepts a user-defined function that returns nodes and weights.
 

docs/src/tutorials/differentiating_integrals.md

@@ -6,7 +6,7 @@ calls. It integrates with ForwardDiff.jl for forward-mode automatic differentiat
 and Zygote.jl for reverse-mode automatic differentiation. For example:
 
 ```@example AD
-using Integrals, ForwardDiff, FiniteDiff, Zygote, IntegralsCuba
+using Integrals, ForwardDiff, FiniteDiff, Zygote, Cuba
 f(x, p) = sum(sin.(x .* p))
 lb = ones(2)
 ub = 3ones(2)

docs/src/tutorials/numerical_integrals.md

@@ -53,7 +53,7 @@ This means that a new output vector is created every time the function `f` is ca
 If we do not want these allocations, we can also define `f` in-position.
 
 ```@example integrate3
-using Integrals, IntegralsCubature
+using Integrals, Cubature
 function f(y, u, p)
  y[1] = sum(sin.(u))
  y[2] = sum(cos.(u))
@@ -74,7 +74,7 @@ The batch interface allows us to compute multiple points at once.
 For example, here we do allocation-free multithreading with Cubature.jl:
 
 ```@example integrate4
-using Integrals, IntegralsCubature, Base.Threads
+using Integrals, Cubature, Base.Threads
 function f(y, u, p)
  Threads.@threads for i in 1:size(u, 2)
  y[1, i] = sum(sin.(@view(u[:, i])))
@@ -97,7 +97,7 @@ the change is a one-argument change:
 
 ```@example integrate5
 using Integrals
-using IntegralsCuba
+using Cuba
 f(u, p) = sum(sin.(u))
 prob = IntegralProblem(f, ones(3), 3ones(3))
 sol = solve(prob, CubaCuhre(); reltol = 1e-3, abstol = 1e-3)

ext/IntegralsCubaExt.jl

@@ -0,0 +1,119 @@
+module IntegralsCubaExt
+
+using Integrals, Cuba
+import Integrals: transformation_if_inf, scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
+ CubaSUAVE, CubaDivonne, CubaCuhre
+
+function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm,
+ sensealg,
+ domain, p;
+ reltol = 1e-8, abstol = 1e-8,
+ maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int))
+ @assert maxiters>=1000 "maxiters for $alg should be larger than 1000"
+ lb, ub = domain
+ mid = (lb + ub) / 2
+ ndim = length(mid)
+ (vol = prod(map(-, ub, lb))) isa Real ||
+ throw(ArgumentError("Cuba.jl only supports real-valued integrands"))
+ # we could support other types by multiplying by the jacobian determinant at the end
+
+ if prob.f isa BatchIntegralFunction
+ # nvec == 1 in Cuba will change vectors to matrices, so we won't support it when
+ # batching
+ (nvec = prob.f.max_batch) > 1 ||
+ throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))
+
+ if mid isa Real
+ _x = zeros(typeof(mid), prob.f.max_batch)
+ scale = x -> scale_x!(resize!(_x, length(x)), ub, lb, vec(x))
+ else
+ _x = zeros(eltype(mid), length(mid), prob.f.max_batch)
+ scale = x -> scale_x!(view(_x, :, 1:size(x, 2)), ub, lb, x)
+ end
+
+ if isinplace(prob)
+ fsize = size(prob.f.integrand_prototype)[begin:(end - 1)]
+ y = similar(prob.f.integrand_prototype, fsize..., nvec)
+ ax = map(_ -> (:), fsize)
+ f = function (x, dx)
+ dy = @view(y[ax..., begin:(begin + size(dx, 2) - 1)])
+ prob.f(dy, scale(x), p)
+ dx .= reshape(dy, :, size(dx, 2)) .* vol
+ end
+ else
+ y = mid isa Number ? prob.f(typeof(mid)[], p) :
+ prob.f(Matrix{typeof(mid)}(undef, length(mid), 0), p)
+ fsize = size(y)[begin:(end - 1)]
+ f = (x, dx) -> dx .= reshape(prob.f(scale(x), p), :, size(dx, 2)) .* vol
+ end
+ ncomp = prod(fsize)
+ else
+ nvec = 1
+
+ if mid isa Real
+ scale = x -> scale_x(ub, lb, only(x))
+ else
+ _x = zeros(eltype(mid), length(mid))
+ scale = x -> scale_x!(_x, ub, lb, x)
+ end
+
+ if isinplace(prob)
+ y = similar(prob.f.integrand_prototype)
+ f = (x, dx) -> dx .= vec(prob.f(y, scale(x), p)) .* vol
+ else
+ y = prob.f(mid, p)
+ f = (x, dx) -> dx .= Iterators.flatten(prob.f(scale(x), p)) .* vol
+ end
+ ncomp = length(y)
+ end
+
+ if alg isa CubaVegas
+ out = Cuba.vegas(f, ndim, ncomp; rtol = reltol,
+ atol = abstol, nvec = nvec,
+ maxevals = maxiters,
+ flags = alg.flags, seed = alg.seed, minevals = alg.minevals,
+ nstart = alg.nstart, nincrease = alg.nincrease,
+ gridno = alg.gridno)
+ elseif alg isa CubaSUAVE
+ out = Cuba.suave(f, ndim, ncomp; rtol = reltol,
+ atol = abstol, nvec = nvec,
+ maxevals = maxiters,
+ flags = alg.flags, seed = alg.seed, minevals = alg.minevals,
+ nnew = alg.nnew, nmin = alg.nmin, flatness = alg.flatness)
+ elseif alg isa CubaDivonne
+ out = Cuba.divonne(f, ndim, ncomp; rtol = reltol,
+ atol = abstol, nvec = nvec,
+ maxevals = maxiters,
+ flags = alg.flags, seed = alg.seed, minevals = alg.minevals,
+ key1 = alg.key1, key2 = alg.key2, key3 = alg.key3,
+ maxpass = alg.maxpass, border = alg.border,
+ maxchisq = alg.maxchisq, mindeviation = alg.mindeviation)
+ elseif alg isa CubaCuhre
+ out = Cuba.cuhre(f, ndim, ncomp; rtol = reltol,
+ atol = abstol, nvec = nvec,
+ maxevals = maxiters,
+ flags = alg.flags, minevals = alg.minevals, key = alg.key)
+ end
+
+ # out.integral is a Vector{Float64}, but we want to return it to the shape of the integrand
+ if prob.f isa BatchIntegralFunction
+ if y isa AbstractVector
+ val = out.integral[1]
+ else
+ val = reshape(out.integral, fsize)
+ end
+ else
+ if y isa Real
+ val = out.integral[1]
+ elseif y isa AbstractVector
+ val = out.integral
+ else
+ val = reshape(out.integral, size(y))
+ end
+ end
+
+ SciMLBase.build_solution(prob, alg, val, out.error,
+ chi = out.probability, retcode = ReturnCode.Success)
+end
+
+end