Fix a bunch of doc builds

docs/Project.toml

@@ -24,6 +24,7 @@ Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 LuxCUDA = "d0bbae9a-e099-4d5b-a835-1c6931763bda"
 LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
+MLDataDevices = "7e8f7934-dd98-4c1a-8fe8-92b47a384d40"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"

docs/make.jl

@@ -26,7 +26,13 @@ makedocs(sitename = "Overview of Julia's SciML",
         "https://bkamins.github.io/julialang/2020/12/24/minilanguage.html",
         "https://arxiv.org/abs/2109.06786",
         "https://arxiv.org/abs/2001.04385",
-        "https://code.visualstudio.com/"],
+        "https://code.visualstudio.com/",
+        "https://www.biorxiv.org/content/10.1101/2020.11.28.402297v2",
+        "https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/docs/src/dev/how_it_works.md",
+        "https://github.com/SciML/Optimization.jl/blob/master/lib/OptimizationPolyalgorithms/src/OptimizationPolyalgorithms.jl",
+        "http://www.scholarpedia.org/article/Differential-algebraic_equations",
+        "https://computing.llnl.gov/projects/sundials"
+        ],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/stable/",
         mathengine = mathengine),
@@ -79,7 +85,6 @@ makedocs(sitename = "Overview of Julia's SciML",
                 "highlevels/developer_documentation.md"],
             "Extra Learning Resources" => ["highlevels/learning_resources.md"]
         ]],
-    warnonly = true,
 )
 
 deploydocs(repo = "github.com/SciML/SciMLDocs")

docs/src/getting_started/first_simulation.md

@@ -269,8 +269,7 @@ import DifferentialEquations as DE
 import ModelingToolkit as MTK
 import ModelingToolkit: t_nounits as t, D_nounits as D, @variables, @parameters, @named
 @parameters α=1.5 β=1.0 γ=3.0 δ=1.0
-@variables t 🐰(t)=1.0 🐺(t)=1.0
-D = MTK.Differential(t)
+@variables 🐰(t)=1.0 🐺(t)=1.0
 eqs = [D(🐰) ~ α * 🐰 - β * 🐰 * 🐺,
     D(🐺) ~ -γ * 🐺 + δ * 🐰 * 🐺]
 

docs/src/highlevels/modeling_languages.md

@@ -80,7 +80,7 @@ making it a solid foundation for any agent-based model.
 ## Unitful.jl: A Julia package for physical units
 
 Supports not only SI units, but also any other unit system.
-[Unitful.jl](https://painterqubits.github.io/Unitful.jl/stable/) has minimal run-time penalty of units.
+[Unitful.jl](https://juliaphysics.github.io/Unitful.jl/stable/) has minimal run-time penalty of units.
 Includes facilities for dimensional analysis, and integrates easily with the usual mathematical operations and collections that are defined in Julia.
 
 ## ReactionMechanismSimulator.jl: Simulation and Analysis of Large Chemical Reaction Systems

docs/src/showcase/brusselator.md

@@ -65,7 +65,7 @@ We wish to obtain the solution to this PDE on a timespan of ``t \in [0,11.5]``.
 
 ## Defining the symbolic PDEsystem with ModelingToolkit.jl
 
-With `ModelingToolkit.jl`, we first symbolically define the system, see also the docs for [`PDESystem`](https://docs.sciml.ai/ModelingToolkit/stable/systems/PDESystem/):
+With `ModelingToolkit.jl`, we first symbolically define the system, see also the docs for [`PDESystem`](https://docs.sciml.ai/ModelingToolkit/stable/API/PDESystem/):
 
 ```@example bruss
 import ModelingToolkit as MTK

docs/src/showcase/missing_physics.md

@@ -410,11 +410,11 @@ regressions:
 ```@example ude
 options = DataDrivenDiffEq.DataDrivenCommonOptions(maxiters = 10_000,
     normalize = DataDrivenDiffEq.DataNormalization(DataDrivenDiffEq.ZScoreTransform),
-    selector = bic, digits = 1,
+    selector = DataDrivenDiffEq.bic, digits = 1,
     data_processing = DataDrivenDiffEq.DataProcessing(split = 0.9,
         batchsize = 30,
         shuffle = true,
-        rng = StableRNG(1111)))
+        rng = StableRNGs.StableRNG(1111)))
 
 full_res = DataDrivenDiffEq.solve(full_problem, basis, opt, options = options)
 full_eqs = DataDrivenDiffEq.get_basis(full_res)
@@ -424,11 +424,11 @@ println(full_res)
 ```@example ude
 options = DataDrivenDiffEq.DataDrivenCommonOptions(maxiters = 10_000,
     normalize = DataDrivenDiffEq.DataNormalization(DataDrivenDiffEq.ZScoreTransform),
-    selector = bic, digits = 1,
+    selector = DataDrivenDiffEq.bic, digits = 1,
     data_processing = DataDrivenDiffEq.DataProcessing(split = 0.9,
         batchsize = 30,
         shuffle = true,
-        rng = StableRNG(1111)))
+        rng = StableRNGs.StableRNG(1111)))
 
 ideal_res = DataDrivenDiffEq.solve(ideal_problem, basis, opt, options = options)
 ideal_eqs = DataDrivenDiffEq.get_basis(ideal_res)
@@ -438,11 +438,11 @@ println(ideal_res)
 ```@example ude
 options = DataDrivenDiffEq.DataDrivenCommonOptions(maxiters = 10_000,
     normalize = DataDrivenDiffEq.DataNormalization(DataDrivenDiffEq.ZScoreTransform),
-    selector = bic, digits = 1,
+    selector = DataDrivenDiffEq.bic, digits = 1,
     data_processing = DataDrivenDiffEq.DataProcessing(split = 0.9,
         batchsize = 30,
         shuffle = true,
-        rng = StableRNG(1111)))
+        rng = StableRNGs.StableRNG(1111)))
 
 nn_res = DataDrivenDiffEq.solve(nn_problem, basis, opt, options = options)
 nn_eqs = DataDrivenDiffEq.get_basis(nn_res)
@@ -525,14 +525,14 @@ p1 = Plots.plot(t, abs.(Array(solution) .- estimate)' .+ eps(Float32),
     legend = :topright)
 
 # Plot L₂
-p2 = plot3d(X̂[1, :], X̂[2, :], Ŷ[2, :], lw = 3,
+p2 = Plots.plot3d(X̂[1, :], X̂[2, :], Ŷ[2, :], lw = 3,
     title = "Neural Network Fit of U2(t)", color = c1,
     label = "Neural Network", xaxis = "x", yaxis = "y",
     titlefont = "Helvetica", legendfont = "Helvetica",
     legend = :bottomright)
 Plots.plot!(X̂[1, :], X̂[2, :], Ȳ[2, :], lw = 3, label = "True Missing Term", color = c2)
 
-p3 = scatter(solution, color = [c1 c2], label = ["x data" "y data"],
+p3 = Plots.scatter(solution, color = [c1 c2], label = ["x data" "y data"],
     title = "Extrapolated Fit From Short Training Data",
     titlefont = "Helvetica", legendfont = "Helvetica",
     markersize = 5)
@@ -542,8 +542,8 @@ Plots.plot!(p3, true_solution_long, color = [c1 c2], linestyle = :dot, lw = 5,
 Plots.plot!(p3, estimate_long, color = [c3 c4], lw = 1,
     label = ["Estimated x(t)" "Estimated y(t)"])
 Plots.plot!(p3, [2.99, 3.01], [0.0, 10.0], lw = 1, color = :black, label = nothing)
-annotate!([(1.5, 13, text("Training \nData", 10, :center, :top, :black, "Helvetica"))])
-l = @layout [grid(1, 2)
-             grid(1, 1)]
+Plots.annotate!([(1.5, 13, Plots.text("Training \nData", 10, :center, :top, :black, "Helvetica"))])
+l = Plots.@layout [Plots.grid(1, 2)
+                   Plots.grid(1, 1)]
 Plots.plot(p1, p2, p3, layout = l)
 ```

docs/src/showcase/ode_types.md

@@ -11,7 +11,7 @@ few useful/interesting types that can be used:
 | ArbFloat                 | [ArbNumerics.jl](https://github.com/JeffreySarnoff/ArbNumerics.jl)                         | More efficient higher precision solutions                       |
 | Measurement              | [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl)                    | Uncertainty propagation                                         |
 | Particles                | [MonteCarloMeasurements.jl](https://github.com/baggepinnen/MonteCarloMeasurements.jl) | Uncertainty propagation                                         |
-| Unitful                  | [Unitful.jl](https://painterqubits.github.io/Unitful.jl/stable/)                      | Unit-checked arithmetic                                         |
+| Unitful                  | [Unitful.jl](https://juliaphysics.github.io/Unitful.jl/stable/)                      | Unit-checked arithmetic                                         |
 | Quaternion               | [Quaternions.jl](https://juliageometry.github.io/Quaternions.jl/stable/)              | Quaternions, duh.                                               |
 | Fun                      | [ApproxFun.jl](https://juliaapproximation.github.io/ApproxFun.jl/latest/)             | Representing PDEs as ODEs in function spaces                    |
 | AbstractOrthoPoly        | [PolyChaos.jl](https://docs.sciml.ai/PolyChaos/stable/)                               | Polynomial Chaos Expansion (PCE) for uncertainty quantification |

docs/src/showcase/optimal_data_gathering_for_missing_physics.md

@@ -11,7 +11,8 @@ To this end, we will rely on the following packages:
 using Random; Random.seed!(984519674645)
 using StableRNGs; rng = StableRNG(845652695)
 import ModelingToolkit as MTK
-import ModelingToolkit: t_nounits as t, D_nounits as D, @mtkcompile, mtkcompile
+import ModelingToolkit: t_nounits as t, D_nounits as D, @mtkmodel, @mtkcompile, mtkcompile
+using ModelingToolkit
 import ModelingToolkitNeuralNets
 import OrdinaryDiffEqRosenbrock as ODE
 import SymbolicIndexingInterface
@@ -58,11 +59,11 @@ This can be implemented in MTK as:
         y_x_s = 0.777
         m = 0.0
     end
-    MTK.@parameters begin
+    @parameters begin
         controls[1:length(optimization_state)-1] = optimization_state[2:end], [tunable = false] # optimization_state is defined further below
         Q_in = optimization_initial, [tunable = false] # similar for optimization state
     end
-    MTK.@variables begin
+    @variables begin
         C_s(t) = 1.0
         C_x(t) = 1.0
         V(t) = 7.0
@@ -106,7 +107,7 @@ We thus extend the bioreactor MTK model with this equation:
 ```@example DoE
 @mtkmodel TrueBioreactor begin
     @extend Bioreactor()
-    MTK.@parameters begin
+    @parameters begin
         μ_max = 0.421
         K_s = 0.439*10
     end
@@ -127,11 +128,11 @@ Similarly, we can extend the bioreactor with a neural network to represent this
                           Lux.Dense(5, 1, x->1*sigmoid(x)))
     end
     @components begin
-        nn = NeuralNetworkBlock(; n_input=1, n_output=1, chain, rng)
+        nn = ModelingToolkitNeuralNets.NeuralNetworkBlock(; n_input=1, n_output=1, chain, rng)
     end
     @equations begin
-        nn.output.u[1] ~ μ
-        nn.input.u[1] ~ C_s
+        nn.outputs[1] ~ μ
+        nn.inputs[1] ~ C_s
     end
 end
 nothing # hide
@@ -195,7 +196,7 @@ function loss(x, (probs, get_varss, datas))
     loss
 end
 of = OPT.OptimizationFunction{true}(loss, SMS.AutoZygote())
-x0 = reduce(vcat, getindex.((default_values(ude_bioreactor),), tunable_parameters(ude_bioreactor)))
+x0 = reduce(vcat, getindex.((MTK.default_values(ude_bioreactor),), MTK.tunable_parameters(ude_bioreactor)))
 get_vars = getu(ude_bioreactor, [ude_bioreactor.C_s])
 ps = ([ude_prob], [get_vars], [data]);
 op = OPT.OptimizationProblem(of, x0, ps)
@@ -411,7 +412,7 @@ plot(ude_sol2[3,:])
 ude_prob_remake = remake(ude_prob, p=ude_prob2.p)
 sol_remake = ODE.solve(ude_prob_remake, ODE.Rodas5P())
 plot(sol_remake[3,:])
-x0 = reduce(vcat, getindex.((default_values(ude_bioreactor),), tunable_parameters(ude_bioreactor)))
+x0 = reduce(vcat, getindex.((MTK.default_values(ude_bioreactor),), MTK.tunable_parameters(ude_bioreactor)))
 
 get_vars2 = getu(ude_bioreactor2, [ude_bioreactor2.C_s])
 
@@ -522,7 +523,7 @@ sol3 = ODE.solve(prob3, ODE.Rodas5P())
 @mtkcompile ude_bioreactor3 = UDEBioreactor()
 ude_prob3 = ODE.ODEProblem(ude_bioreactor3, [], (0.0, 15.0), tstops=0:15, save_everystep=false)
 
-x0 = reduce(vcat, getindex.((default_values(ude_bioreactor3),), tunable_parameters(ude_bioreactor3)))
+x0 = reduce(vcat, getindex.((MTK.default_values(ude_bioreactor3),), MTK.tunable_parameters(ude_bioreactor3)))
 
 get_vars3 = getu(ude_bioreactor3, [ude_bioreactor3.C_s])
 
@@ -556,4 +557,4 @@ like with a double division.
 This is because symbolic regression considers multiplication and division to have the same complexity.
 
 In this tutorial, we have shown that experimental design can be used to explore the state space of a dynamic system in a thoughtful way,
-such that missing physics can be recovered in an efficient manner.
+such that missing physics can be recovered in an efficient manner.

docs/src/showcase/optimization_under_uncertainty.md

@@ -9,6 +9,7 @@ First let's consider a 2D bouncing ball, where the states are the horizontal pos
 ```@example control
 import DifferentialEquations as DE
 import Plots
+import Statistics
 
 function ball!(du, u, p, t)
     du[1] = u[2]
@@ -19,7 +20,7 @@ end
 
 ground_condition(u, t, integrator) = u[3]
 ground_affect!(integrator) = integrator.u[4] = -integrator.p[2] * integrator.u[4]
-ground_cb = DE.DE.ContinuousCallback(ground_condition, ground_affect!)
+ground_cb = DE.ContinuousCallback(ground_condition, ground_affect!)
 
 u0 = [0.0, 2.0, 50.0, 0.0]
 tspan = (0.0, 50.0)
@@ -34,8 +35,8 @@ For this particular problem, we wish to measure the impact distance from a point
 
 ```@example control
 stop_condition(u, t, integrator) = u[1] - 25.0
-stop_cb = DE.DE.ContinuousCallback(stop_condition, DE.terminate!)
-cbs = DE.DE.CallbackSet(ground_cb, stop_cb)
+stop_cb = DE.ContinuousCallback(stop_condition, DE.terminate!)
+cbs = DE.CallbackSet(ground_cb, stop_cb)
 
 tspan = (0.0, 1500.0)
 prob = DE.ODEProblem(ball!, u0, tspan, p)
@@ -46,7 +47,7 @@ Plots.plot(sol, vars = (1, 3), label = nothing, xlabel = "x", ylabel = "y")
 To help visualize this problem, we plot as follows, where the star indicates a desired impact location
 
 ```@example control
-rectangle(xc, yc, w, h) = Shape(xc .+ [-w, w, w, -w] ./ 2.0, yc .+ [-h, -h, h, h] ./ 2.0)
+rectangle(xc, yc, w, h) = Plots.Shape(xc .+ [-w, w, w, -w] ./ 2.0, yc .+ [-h, -h, h, h] ./ 2.0)
 
 begin
     Plots.plot(sol, vars = (1, 3), label = nothing, lw = 3, c = :black)
@@ -68,9 +69,9 @@ import Distributions
 cor_dist = Distributions.truncated(Distributions.Normal(0.9, 0.02), 0.9 - 3 * 0.02, 1.0)
 trajectories = 100
 
-prob_func(prob, i, repeat) = DE.DE.remake(prob, p = [p[1], rand(cor_dist)])
-ensemble_prob = DE.DE.EnsembleProblem(prob, prob_func = prob_func)
-ensemblesol = DE.solve(ensemble_prob, DE.Tsit5(), DE.DE.EnsembleThreads(), trajectories = trajectories,
+prob_func(prob, i, repeat) = DE.remake(prob, p = [p[1], rand(cor_dist)])
+ensemble_prob = DE.EnsembleProblem(prob, prob_func = prob_func)
+ensemblesol = DE.solve(ensemble_prob, DE.Tsit5(), DE.EnsembleThreads(), trajectories = trajectories,
     callback = cbs)
 
 begin # plot
@@ -97,10 +98,10 @@ obs(sol, p) = abs2(sol[3, end] - 25)
 With the observable defined, we can compute the expected squared miss distance from our Monte Carlo simulation results as
 
 ```@example control
-mean_ensemble = mean([obs(sol, p) for sol in ensemblesol])
+mean_ensemble = Statistics.mean([obs(sol, p) for sol in ensemblesol])
 ```
 
-Alternatively, we can use the `Koopman()` algorithm in SciMLExpectations.jl to compute this expectation much more efficiently as
+Alternatively, we can use the `SciMLExpectations.Koopman()` algorithm in SciMLExpectations.jl to compute this expectation much more efficiently as
 
 ```@example control
 import SciMLExpectations
@@ -230,9 +231,9 @@ end
 Using the previously computed optimal initial conditions, let's compute the probability of hitting this wall
 
 ```@example control
-sm = SystemMap(DE.remake(prob, u0 = make_u0(minx)), DE.Tsit5(), callback = cbs)
-exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
-sol = DE.solve(exprob, Koopman(), ireltol = 1e-5)
+sm = SciMLExpectations.SystemMap(DE.remake(prob, u0 = make_u0(minx)), DE.Tsit5(), callback = cbs)
+exprob = SciMLExpectations.ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
+sol = DE.solve(exprob, SciMLExpectations.Koopman(), ireltol = 1e-5)
 sol.u
 ```
 
@@ -241,17 +242,17 @@ We then set up the constraint function for NLopt just as before.
 ```@example control
 function 𝔼_constraint(res, θ, pars)
     prob = DE.ODEProblem(ball!, make_u0(θ), tspan, p)
-    sm = SystemMap(prob, DE.Tsit5(), callback = cbs)
-    exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
-    sol = DE.solve(exprob, Koopman(), ireltol = 1e-5)
+    sm = SciMLExpectations.SystemMap(prob, DE.Tsit5(), callback = cbs)
+    exprob = SciMLExpectations.ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
+    sol = DE.solve(exprob, SciMLExpectations.Koopman(), ireltol = 1e-5)
     res .= sol.u
 end
 opt_lcons = [-Inf]
 opt_ucons = [0.01]
-optimizer = OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
+optimizer = OptimizationMOI.MOI.OptimizerWithAttributes(OptimizationNLopt.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,
+opt_f = OPT.OptimizationFunction(𝔼_loss, OPT.AutoForwardDiff(), cons = 𝔼_constraint)
+opt_prob = OPT.OptimizationProblem(opt_f, opt_ini; lb = opt_lb, ub = opt_ub, lcons = opt_lcons,
     ucons = opt_ucons)
 opt_sol = DE.solve(opt_prob, optimizer)
 minx2 = opt_sol.u

docs/src/showcase/pinngpu.md

@@ -22,7 +22,8 @@ our packages look like:
 # High Level Interface
 import NeuralPDE
 import ModelingToolkit as MTK
-using MTK: @parameters, @variables, Differential, Interval, PDESystem
+using ModelingToolkit: @parameters, @variables, Differential, PDESystem
+using DomainSets: Interval
 
 # Optimization Libraries
 import Optimization as OPT
@@ -32,7 +33,8 @@ import OptimizationOptimisers
 import Lux
 import LuxCUDA
 import ComponentArrays
-const gpud = LuxCUDA.gpu_device() # allocate a GPU device
+import MLDataDevices
+const gpud = MLDataDevices.gpu_device() # allocate a GPU device
 
 # Standard Libraries
 import Printf
@@ -100,11 +102,11 @@ bcs = [u(t_min, x, y) ~ analytic_sol_func(t_min, x, y),
     u(t, x, y_max) ~ analytic_sol_func(t, x, y_max)]
 
 # Space and time domains
-domains = [t ∈ MTK.Interval(t_min, t_max),
-    x ∈ MTK.Interval(x_min, x_max),
-    y ∈ MTK.Interval(y_min, y_max)]
+domains = [t ∈ Interval(t_min, t_max),
+    x ∈ Interval(x_min, x_max),
+    y ∈ Interval(y_min, y_max)]
 
-@named pde_system = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
+MTK.@named pde_system = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
 ```
 
 !!! note