Fixes

Project.toml

@@ -30,7 +30,7 @@ DiffEqNoiseProcess = "5"
 Distributions = "0.23, 0.24, 0.25"
 FiniteDiff = "2"
 ForwardDiff = "0.10"
-Integrals = "4"
+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"

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/src/tutorials/gpu_bayesian.md

@@ -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]
@@ -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

@@ -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