Merge pull request #148 from ArnoStrouwen/Integrals4

Integrals 4 upgrade

.github/workflows/CI.yml

@@ -12,14 +12,18 @@ on:
       - 'docs/**'
 jobs:
   test:
-    runs-on: ubuntu-latest
+    runs-on: ${{ matrix.os }}
     strategy:
+      fail-fast: false
       matrix:
         group:
           - Core
         version:
           - '1'
-          - '1.6'
+        os:
+          - ubuntu-latest
+          - macos-latest
+          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -23,14 +23,14 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 Aqua = "0.8"
 BenchmarkTools = "1"
 ComponentArrays = "0.15"
+Cuba = "2"
+Cubature = "1.5"
 DiffEqBase = "6.100"
 DiffEqNoiseProcess = "5"
 Distributions = "0.23, 0.24, 0.25"
 FiniteDiff = "2"
 ForwardDiff = "0.10"
-Integrals = "3.1"
-IntegralsCuba = "0.2"
-IntegralsCubature = "0.2"
+Integrals = "4.4"
 KernelDensity = "0.6"
 LinearAlgebra = "1"
 OrdinaryDiffEq = "6"
@@ -42,7 +42,7 @@ RecursiveArrayTools = "2, 3"
 Reexport = "0.2, 1.0"
 ReferenceTests = "0.10"
 SafeTestsets = "0.1"
-SciMLBase = "1, 2"
+SciMLBase = "2"
 StaticArrays = "1"
 Statistics = "1"
 StochasticDiffEq = "6"
@@ -56,11 +56,11 @@ julia = "1.6"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
+Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 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"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
@@ -74,4 +74,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 
 [targets]
-test = ["Test", "Aqua", "SafeTestsets", "TestExtras", "Pkg", "PkgBenchmark", "BenchmarkTools", "ReferenceTests", "SymPy", "Random", "OrdinaryDiffEq", "ComponentArrays", "ForwardDiff", "Integrals", "IntegralsCubature", "IntegralsCuba", "FiniteDiff", "RecursiveArrayTools", "StochasticDiffEq"]
+test = ["Test", "Aqua", "SafeTestsets", "TestExtras", "Pkg", "PkgBenchmark", "BenchmarkTools", "ReferenceTests", "SymPy", "Random", "OrdinaryDiffEq", "ComponentArrays", "ForwardDiff", "Integrals", "Cubature", "Cuba", "FiniteDiff", "RecursiveArrayTools", "StochasticDiffEq"]

README.md

@@ -26,7 +26,7 @@ the documentation, which contains the unreleased features.
 ### Example
 
 ```julia
-using SciMLExpectations, OrdinaryDiffEq, Distributions, IntegralsCubature
+using SciMLExpectations, OrdinaryDiffEq, Distributions, Cubature
 
 function eom!(du, u, p, t, A)
     du .= A * u
@@ -56,7 +56,7 @@ julia> analytical
 
 ```julia
 g(sol, p) = sol[:, end]
-exprob = ExpectationProblem(sm, g, cov, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, cov, gd)
 sol = solve(exprob, Koopman(); quadalg = CubatureJLh(),
             ireltol = 1e-3, iabstol = 1e-3)
 sol.u # Expectation of the states 1 and 2 at the final time point

docs/Project.toml

@@ -1,11 +1,11 @@
 [deps]
+Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
 DiffEqGPU = "071ae1c0-96b5-11e9-1965-c90190d839ea"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-IntegralsCuba = "e00cd5f1-6337-4131-8b37-28b2fe4cd6cb"
 KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
@@ -26,7 +26,6 @@ DifferentialEquations = "7"
 Distributions = "0.25"
 Documenter = "1"
 ForwardDiff = "0.10"
-IntegralsCuba = "0.3"
 KernelDensity = "0.6"
 MCMCChains = "6"
 Optimization = "3"

docs/src/tutorials/gpu_bayesian.md

@@ -22,7 +22,7 @@ Let's start by importing all the necessary libraries:
 using OrdinaryDiffEq
 using Turing, MCMCChains, Distributions
 using KernelDensity
-using SciMLExpectations, IntegralsCuba
+using SciMLExpectations, Cuba
 using DiffEqGPU
 using Plots, StatsPlots
 using Random;
@@ -166,7 +166,7 @@ ub = [p_kde_i.x[end] for p_kde_i in p_kde] #upper bound for integration
 gd = GenericDistribution(pdf_func, rand_func, lb, ub)
 h(x, u, p) = u, x
 sm = SystemMap(prob1, Tsit5())
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman(); reltol = 1e-2, abstol = 1e-2)
 sol.u
 ```

docs/src/tutorials/introduction.md

@@ -62,7 +62,7 @@ using SciMLExpectations
 gd = GenericDistribution(u0_dist...)
 h(x, u, p) = x, p
 sm = SystemMap(prob, Tsit5())
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, MonteCarlo(100000))
 sol.u
 ```
@@ -115,7 +115,7 @@ Changing the distribution, we arrive at
 ```@example introduction
 u0_dist = [Uniform(0.0, 10.0)]
 gd = GenericDistribution(u0_dist...)
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, MonteCarlo(100000))
 sol.u
 ```
@@ -125,7 +125,7 @@ and
 ```@example introduction
 u0_dist = [Uniform(0.0, 10.0)]
 gd = GenericDistribution(u0_dist...)
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -141,7 +141,7 @@ Note that the `Koopman()` algorithm doesn't currently support infinite or semi-i
 ```@example introduction
 u0_dist = [Normal(3.0, 2.0)]
 gd = GenericDistribution(u0_dist...)
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -157,7 +157,7 @@ Using a truncated distribution will alleviate this problem. However, there is an
 ```@example introduction
 u0_dist = [truncated(Normal(3.0, 2.0), -1000, 1000)]
 gd = GenericDistribution(u0_dist...)
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -167,7 +167,7 @@ whereas truncating at $\pm 4\sigma$ produces the correct result
 ```@example introduction
 u0_dist = [truncated(Normal(3.0, 2.0), -5, 11)]
 gd = GenericDistribution(u0_dist...)
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -177,20 +177,20 @@ If a large truncation is required, it is best practice to center the distributio
 ```@example introduction
 u0_dist = [truncated(Normal(3.0, 2.0), 3 - 1000, 3 + 1000)]
 gd = GenericDistribution(u0_dist...)
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(u0))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
 
 ## Vector-Valued Functions
 
-`ExpectationProblem` can also handle vector-valued functions. Simply pass the vector-valued function and set the `nout` kwarg to the length of the vector the function returns.
+`ExpectationProblem` can also handle vector-valued functions.
 
 Here, we demonstrate this by computing the expectation of `u` at `t=4.0s` and `t=6.0s`
 
 ```@example introduction
 g(sol, p) = [sol(4.0)[1], sol(6.0)[1]]
-exprob = ExpectationProblem(sm, g, h, gd; nout = 2)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -208,7 +208,7 @@ saveat = tspan[1]:0.5:tspan[2]
 g(sol, p) = sol[1, :]
 prob = ODEProblem(f, u0, tspan, p, saveat = saveat)
 sm = SystemMap(prob, Tsit5())
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(saveat))
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -234,21 +234,21 @@ end
 counter = [0, 0, 0]
 g(sol, p, power) = g(sol, p, power, counter)
 g(sol, p) = g(sol, p, 1)
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
 
 ```@example introduction
 g(sol, p) = g(sol, p, 2)
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
 
 ```@example introduction
 g(sol, p) = g(sol, p, 3)
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -269,7 +269,7 @@ function g(sol, p, counter)
 end
 counter = [0]
 g(sol, p) = g(sol, p, counter)
-exprob = ExpectationProblem(sm, g, h, gd; nout = 3)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 ```
@@ -305,7 +305,7 @@ function g(sol, p)
     x = sol(4.0)[1]
     [x, x^2]
 end
-exprob = ExpectationProblem(sm, g, h, gd; nout = 2)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 mean_g = sol.u[1]
@@ -327,7 +327,7 @@ saveat = tspan[1]:0.5:tspan[2]
 g(sol, p) = [sol[1, :]; sol[1, :] .^ 2]
 prob = ODEProblem(f, u0, tspan, p, saveat = saveat)
 sm = SystemMap(prob, Tsit5())
-exprob = ExpectationProblem(sm, g, h, gd; nout = length(saveat) * 2)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 
@@ -350,7 +350,7 @@ function g(sol, p)
     x = sol(4.0)[1]
     [x, x^2, x^3]
 end
-exprob = ExpectationProblem(sm, g, h, gd; nout = 3)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman())
 sol.u
 mean_g = sol.u[1]
@@ -367,7 +367,7 @@ It is also possible to solve the various simulations in parallel by using the `b
 The default quadrature algorithm to solve `ExpectationProblem` does not support batch-mode evaluation. So, we first load dependencies for additional quadrature algorithms
 
 ```@example introduction
-using IntegralsCuba
+using Cuba
 ```
 
 We then solve our expectation as before, using a `batch=10` multi-thread parallelization via `EnsembleThreads()` of Cuba's SUAVE algorithm. We also introduce additional uncertainty in the model parameter.
@@ -379,7 +379,7 @@ gd = GenericDistribution(u0_dist, p_dist)
 g(sol, p) = sol(6.0)[1]
 h(x, u, p) = [x[1]], [x[2]]
 sm = SystemMap(prob, Tsit5(), EnsembleThreads())
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 # batchmode = EnsembleThreads() #where to pass this?
 sol = solve(exprob, Koopman(), batch = 10, quadalg = CubaSUAVE())
 sol.u
@@ -392,8 +392,8 @@ Now, let's compare the performance of the batch and non-batch modes
 ```
 
 ```@example introduction
-solve(exprob, Koopman(), batch = 0, quadalg = CubaSUAVE())
-@time solve(exprob, Koopman(), batch = 0, quadalg = CubaSUAVE())
+solve(exprob, Koopman(), quadalg = CubaSUAVE())
+@time solve(exprob, Koopman(), quadalg = CubaSUAVE())
 ```
 
 It is also possible to parallelize across the GPU. However, one must be careful of the limitations of ensemble solutions with the GPU. Please refer to [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl) for details.
@@ -428,7 +428,7 @@ g(sol, p) = sol(6.0f0)[1]
 h(x, u, p) = [x[1]], [x[2]]
 prob = ODEProblem(f, u0, tspan, p)
 sm = SystemMap(prob, Tsit5(), EnsembleCPUArray())
-exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, g, h, gd)
 sol = solve(exprob, Koopman(), batch = 10, quadalg = CubaSUAVE())
 sol.u
 ```

docs/src/tutorials/optimization_under_uncertainty.md

@@ -106,7 +106,7 @@ using SciMLExpectations
 gd = GenericDistribution(cor_dist)
 h(x, u, p) = u, [p[1]; x[1]]
 sm = SystemMap(prob, Tsit5(), callback = cbs)
-exprob = ExpectationProblem(sm, obs, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, obs, h, gd)
 sol = solve(exprob, Koopman(), ireltol = 1e-5)
 sol.u
 ```
@@ -121,7 +121,7 @@ make_u0(θ) = [θ[1], θ[2], θ[3], 0.0]
 function 𝔼_loss(θ, pars)
     prob = ODEProblem(ball!, make_u0(θ), tspan, p)
     sm = SystemMap(prob, Tsit5(), callback = cbs)
-    exprob = ExpectationProblem(sm, obs, h, gd; nout = 1)
+    exprob = ExpectationProblem(sm, obs, h, gd)
     sol = solve(exprob, Koopman(), ireltol = 1e-5)
     sol.u
 end
@@ -228,7 +228,7 @@ Using the previously computed optimal initial conditions, let's compute the prob
 
 ```@example control
 sm = SystemMap(remake(prob, u0 = make_u0(minx)), Tsit5(), callback = cbs)
-exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
+exprob = ExpectationProblem(sm, constraint_obs, h, gd)
 sol = solve(exprob, Koopman(), ireltol = 1e-5)
 sol.u
 ```
@@ -239,7 +239,7 @@ We then set up the constraint function for NLopt just as before.
 function 𝔼_constraint(res, θ, pars)
     prob = ODEProblem(ball!, make_u0(θ), tspan, p)
     sm = SystemMap(prob, Tsit5(), callback = cbs)
-    exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
+    exprob = ExpectationProblem(sm, constraint_obs, h, gd)
     sol = solve(exprob, Koopman(), ireltol = 1e-5)
     res .= sol.u
 end

docs/src/tutorials/process_noise.md

@@ -5,7 +5,7 @@ This done by representing the Wiener process using the [Kosambi–Karhunen–Lo
 
 ```@example process_noise
 using SciMLExpectations
-using IntegralsCuba
+using Cuba
 using StochasticDiffEq
 using DiffEqNoiseProcess
 using Distributions
@@ -19,7 +19,7 @@ prob = SDEProblem(f, g, u0, (0.0, 1.0), noise = W)
 sm = ProcessNoiseSystemMap(prob, 8, LambaEM(), abstol = 1e-3, reltol = 1e-3)
 cov(x, u, p) = x, p
 observed(sol, p) = sol[:, end]
-exprob = ExpectationProblem(sm, observed, cov; nout = length(u0))
+exprob = ExpectationProblem(sm, observed, cov)
 sol1 = solve(exprob, Koopman(), ireltol = 1e-3, iabstol = 1e-3, batch = 64,
              quadalg = CubaDivonne())
 sol1.u