Documenter 1.0 upgrade

.github/workflows/CI.yml

@@ -3,9 +3,13 @@ on:
   pull_request:
     branches:
       - master
+    paths-ignore:
+      - 'docs/**'
   push:
     branches:
       - master
+    paths-ignore:
+      - 'docs/**'
 jobs:
   test:
     runs-on: ubuntu-latest

.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()'
+        run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs"])'

docs/Project.toml

@@ -1,13 +1,16 @@
 [deps]
 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"
-NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SciMLExpectations = "afe9f18d-7609-4d0e-b7b7-af0cb72b8ea8"
@@ -17,4 +20,4 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
-Documenter = "0.27"
+Documenter = "1"

docs/make.jl

@@ -8,18 +8,9 @@ include("pages.jl")
 makedocs(sitename = "SciMLExpectations.jl",
          authors = "Chris Rackauckas",
          modules = [SciMLExpectations],
-         linkcheck = true,
-         strict = [
-             :doctest,
-             :linkcheck,
-             :parse_error,
-             :example_block,
-             # Other available options are
-             # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
-         ],
-         clean = true, doctest = false,
-         format = Documenter.HTML(analytics = "UA-90474609-3",
-                                  assets = ["assets/favicon.ico"],
+         clean = true, doctest = false, linkcheck = true,
+         warnonly = [:missing_docs],
+         format = Documenter.HTML(assets = ["assets/favicon.ico"],
                                   canonical = "https://docs.sciml.ai/SciMLExpectations/stable/"),
          pages = pages)
 

docs/src/index.md

@@ -75,32 +75,19 @@ Pkg.status(; mode = PKGMODE_MANIFEST) # hide
 </details>
 ```
 
-```@raw html
-You can also download the 
-<a href="
-```
-
 ```@eval
 using TOML
+using Markdown
 version = TOML.parse(read("../../Project.toml", String))["version"]
 name = TOML.parse(read("../../Project.toml", String))["name"]
-link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
-       "/assets/Manifest.toml"
-```
-
-```@raw html
-">manifest</a> file and the
-<a href="
-```
-
-```@eval
-using TOML
-version = TOML.parse(read("../../Project.toml", String))["version"]
-name = TOML.parse(read("../../Project.toml", String))["name"]
-link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
-       "/assets/Project.toml"
-```
-
-```@raw html
-">project</a> file.
+link_manifest = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
+                "/assets/Manifest.toml"
+link_project = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
+               "/assets/Project.toml"
+Markdown.parse("""You can also download the
+[manifest]($link_manifest)
+file and the
+[project]($link_project)
+file.
+""")
 ```

docs/src/manual/solve.md

@@ -1,5 +1,6 @@
 # Solving Expectation Problems
 
 ```@docs
-solve(prob::ExpectationProblem, expalg, args...)
+solve(prob::ExpectationProblem, expalg::MonteCarlo)
+solve(prob::ExpectationProblem, expalg::Koopman)
 ```

docs/src/tutorials/optimization_under_uncertainty.md

@@ -1,13 +1,13 @@
-# Optimization Under Uncertainty with SciMLExpectations.jl
+# [Optimization Under Uncertainty](@id optimization_under_uncertainty)
 
-This tutorial gives an overview of how to leverage the efficient Koopman expectation method from SciMLExpectations to perform optimization under uncertainty. We demonstrate this by using a bouncing ball model with an uncertain model parameter. We also demonstrate its application to problems with probabilistic constraints, in particular a special class of constraints called chance constraints.
+This tutorial showcases how to leverage the efficient Koopman expectation method from SciMLExpectations to perform optimization under uncertainty. We demonstrate this by using a bouncing ball model with an uncertain model parameter. We also demonstrate its application to problems with probabilistic constraints, in particular a special class of constraints called chance constraints.
 
 ## System Model
 
 First let's consider a 2D bouncing ball, where the states are the horizontal position $x$, horizontal velocity $\dot{x}$, vertical position $y$, and vertical velocity $\dot{y}$. This model has two system parameters, acceleration due to gravity and coefficient of restitution (models energy loss when the ball impacts the ground). We can simulate such a system using `ContinuousCallback` as
 
 ```@example control
-using OrdinaryDiffEq, Plots
+using DifferentialEquations, Plots
 
 function ball!(du, u, p, t)
     du[1] = u[2]
@@ -27,7 +27,6 @@ p = [9.807, 0.9]
 prob = ODEProblem(ball!, u0, tspan, p)
 sol = solve(prob, Tsit5(), callback = ground_cb)
 plot(sol, vars = (1, 3), label = nothing, xlabel = "x", ylabel = "y")
-plot!(background_color = :transparent)
 ```
 
 For this particular problem, we wish to measure the impact distance from a point $y=25$ on a wall at $x=25$. So, we introduce an additional callback that terminates the simulation on wall impact.
@@ -71,7 +70,7 @@ trajectories = 100
 prob_func(prob, i, repeat) = remake(prob, p = [p[1], rand(cor_dist)])
 ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
 ensemblesol = solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories = trajectories,
-                    callback = cbs)
+    callback = cbs)
 
 begin # plot
     plot(ensemblesol, vars = (1, 3), lw = 1)
@@ -116,53 +115,38 @@ sol.u
 
 We now wish to optimize the initial position ($x_0,y_0$) and horizontal velocity ($\dot{x}_0$) of the system to minimize the expected squared miss distance from the star, where $x_0\in\left[-100,0\right]$, $y_0\in\left[1,3\right]$, and $\dot{x}_0\in\left[10,50\right]$. We will demonstrate this using a gradient-based optimization approach from NLopt.jl using `ForwardDiff.jl` AD through the expectation calculation.
 
-First, we load the required packages and define our loss function
-
 ```@example control
-using NLopt, SciMLSensitivity, ForwardDiff
-
+using Optimization, OptimizationNLopt, OptimizationMOI
 make_u0(θ) = [θ[1], θ[2], θ[3], 0.0]
-
-function 𝔼_loss(θ)
+function 𝔼_loss(θ, pars)
     prob = ODEProblem(ball!, make_u0(θ), tspan, p)
     sm = SystemMap(prob, Tsit5(), callback = cbs)
     exprob = ExpectationProblem(sm, obs, h, gd; nout = 1)
     sol = solve(exprob, Koopman(), ireltol = 1e-5)
     sol.u
 end
-```
-
-NLopt requires that this loss function return the loss as above, but also do an inplace update of the gradient. So, we wrap this function to put it in the form required by NLopt.
-
-```@example control
-function 𝔼_loss_nlopt(θ, ∇)
-    length(∇) > 0 ? ForwardDiff.gradient!(∇, 𝔼_loss, θ) : nothing
-    𝔼_loss(θ)
-end
-```
-
-We then optimize using the [Method of Moving Asymptotes](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#mma-method-of-moving-asymptotes-and-ccsa) algorithm (`:LD_MMA`)
-
-```@example control
-opt = Opt(:LD_MMA, 3)
-opt.lower_bounds = [-100.0, 1.0, 10.0]
-opt.upper_bounds = [0.0, 3.0, 50.0]
-opt.xtol_rel = 1e-3
-opt.min_objective = 𝔼_loss_nlopt
-(minf, minx, ret) = NLopt.optimize(opt, [-1.0, 2.0, 50.0])
+opt_f = OptimizationFunction(𝔼_loss, Optimization.AutoForwardDiff())
+opt_ini = [-1.0, 2.0, 50.0]
+opt_lb = [-100.0, 1.0, 10.0]
+opt_ub = [0.0, 3.0, 50.0]
+opt_prob = OptimizationProblem(opt_f, opt_ini; lb = opt_lb, ub = opt_ub)
+optimizer = OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
+    "algorithm" => :LD_MMA)
+opt_sol = solve(opt_prob, optimizer)
+minx = opt_sol.u
 ```
 
 Let's now visualize 100 Monte Carlo simulations
 
 ```@example control
 ensembleprob = EnsembleProblem(remake(prob, u0 = make_u0(minx)), prob_func = prob_func)
 ensemblesol = solve(ensembleprob, Tsit5(), EnsembleThreads(), trajectories = 100,
-                    callback = cbs)
+    callback = cbs)
 
 begin
     plot(ensemblesol, vars = (1, 3), lw = 1, alpha = 0.1)
     plot!(solve(remake(prob, u0 = make_u0(minx)), Tsit5(), callback = cbs),
-          vars = (1, 3), label = nothing, c = :black, lw = 3, ls = :dash)
+        vars = (1, 3), label = nothing, c = :black, lw = 3, ls = :dash)
     xlabel!("x [m]")
     ylabel!("y [m]")
     plot!(rectangle(27.5, 25, 5, 50), c = :red, label = nothing)
@@ -175,7 +159,7 @@ end
 Looks pretty good! But, how long did it take? Let's benchmark.
 
 ```@example control
-@time NLopt.optimize(opt, [-1.0, 2.0, 50.0])
+@time solve(opt_prob, optimizer)
 ```
 
 Not bad for bound constrained optimization under uncertainty of a hybrid system!
@@ -191,7 +175,7 @@ begin
     xlabel!("x [m]")
     ylabel!("y [m]")
     plot!([constraint[1], constraint[1]], [0.0, constraint[2]], lw = 5, c = :black,
-          label = nothing)
+        label = nothing)
     scatter!([25], [25], marker = :star, ms = 10, label = nothing, c = :green)
     ylims!(0.0, 50.0)
     xlims!(minx[1], 27.5)
@@ -208,16 +192,16 @@ function constraint_affect!(integrator)
     integrator.u[3] < constraint[2] ? terminate!(integrator) : nothing
 end
 constraint_cb = ContinuousCallback(constraint_condition, constraint_affect!,
-                                   save_positions = (true, false));
+    save_positions = (true, false));
 constraint_cbs = CallbackSet(ground_cb, stop_cb, constraint_cb)
 
 ensemblesol = solve(ensembleprob, Tsit5(), EnsembleThreads(), trajectories = 500,
-                    callback = constraint_cbs)
+    callback = constraint_cbs)
 
 begin
     plot(ensemblesol, vars = (1, 3), lw = 1, alpha = 0.1)
     plot!(solve(remake(prob, u0 = make_u0(minx)), Tsit5(), callback = constraint_cbs),
-          vars = (1, 3), label = nothing, c = :black, lw = 3, ls = :dash)
+        vars = (1, 3), label = nothing, c = :black, lw = 3, ls = :dash)
 
     xlabel!("x [m]")
     ylabel!("y [m]")
@@ -234,13 +218,16 @@ That doesn't look good!
 We now need a second observable for the system. To compute a probability of impact, we use an indicator function for if a trajectory impacts the wall. In other words, this functions returns 1 if the trajectory hits the wall and 0 otherwise.
 
 ```@example control
-constraint_obs(sol, p) = sol[1, end] ≈ constraint[1] ? one(sol[1, end]) : zero(sol[1, end])
+function constraint_obs(sol, p)
+    sol((constraint[1] - sol[1, 1]) / sol[2, 1])[3] <= constraint[2] ? one(sol[1, end]) :
+    zero(sol[1, end])
+end
 ```
 
 Using the previously computed optimal initial conditions, let's compute the probability of hitting this wall
 
 ```@example control
-sm = SystemMap(remake(prob, u0 = make_u0(minx)), Tsit5(), callback = constraint_cbs)
+sm = SystemMap(remake(prob, u0 = make_u0(minx)), Tsit5(), callback = cbs)
 exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
 sol = solve(exprob, Koopman(), ireltol = 1e-5)
 sol.u
@@ -249,56 +236,49 @@ sol.u
 We then set up the constraint function for NLopt just as before.
 
 ```@example control
-function 𝔼_constraint(θ)
+function 𝔼_constraint(res, θ, pars)
     prob = ODEProblem(ball!, make_u0(θ), tspan, p)
-    sm = SystemMap(prob, Tsit5(), callback = constraint_cbs)
+    sm = SystemMap(prob, Tsit5(), callback = cbs)
     exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
     sol = solve(exprob, Koopman(), ireltol = 1e-5)
-    sol.u
-end
-function 𝔼_constraint_nlopt(x, ∇)
-    length(∇) > 0 ? ForwardDiff.gradient!(∇, 𝔼_constraint, x) : nothing
-    𝔼_constraint(x) - 0.01
+    res .= sol.u
 end
-```
-
-Note that NLopt requires the constraint function to be of the form $g(x) \leq 0$. Hence, why we return `𝔼_constraint(x) - 0.01` for the 1% chance constraint.
-
-The rest of the NLopt setup looks the same as before, except for adding the inequality constraint
-
-```@example control
-opt = Opt(:LD_MMA, 3)
-opt.lower_bounds = [-100.0, 1.0, 10.0]
-opt.upper_bounds = [0.0, 3.0, 50.0]
-opt.xtol_rel = 1e-3
-opt.min_objective = 𝔼_loss_nlopt
-inequality_constraint!(opt, 𝔼_constraint_nlopt, 1e-5)
-(minf2, minx2, ret2) = NLopt.optimize(opt, [-1.0, 2.0, 50.0])
+opt_lcons = [-Inf]
+opt_ucons = [0.01]
+optimizer = OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
+    "algorithm" => :LD_MMA)
+opt_f = OptimizationFunction(𝔼_loss, Optimization.AutoForwardDiff(), cons = 𝔼_constraint)
+opt_prob = OptimizationProblem(opt_f, opt_ini; lb = opt_lb, ub = opt_ub, lcons = opt_lcons,
+    ucons = opt_ucons)
+opt_sol = solve(opt_prob, optimizer)
+minx2 = opt_sol.u
 ```
 
 The probability of impacting the wall is now
 
 ```@example control
-λ = 𝔼_constraint(minx2)
+container = zeros(1)
+𝔼_constraint(container, minx2, nothing)
+λ = container[1]
 ```
 
 We can check if this is within tolerance by
 
 ```@example control
-λ - 0.01 <= 1e-5
+λ <= 0.01 + 1e-5
 ```
 
 Again, we plot some Monte Carlo simulations from this result as follows
 
 ```@example control
 ensembleprob = EnsembleProblem(remake(prob, u0 = make_u0(minx2)), prob_func = prob_func)
 ensemblesol = solve(ensembleprob, Tsit5(), EnsembleThreads(),
-                    trajectories = 500, callback = constraint_cbs)
+    trajectories = 500, callback = constraint_cbs)
 
 begin
     plot(ensemblesol, vars = (1, 3), lw = 1, alpha = 0.1)
     plot!(solve(remake(prob, u0 = make_u0(minx2)), Tsit5(), callback = constraint_cbs),
-          vars = (1, 3), label = nothing, c = :black, lw = 3, ls = :dash)
+        vars = (1, 3), label = nothing, c = :black, lw = 3, ls = :dash)
     plot!([constraint[1], constraint[1]], [0.0, constraint[2]], lw = 5, c = :black)
 
     xlabel!("x [m]")