diff --git a/README.md b/README.md index e713985a75a..947e6d4115a 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,3 @@ [![Documentation and Deployment](https://badge.buildkite.com/b4720f5f3a8d2d1ac8bbaebf5e0e38aaae28b6c8fb436ba05e.svg)](https://buildkite.com/julialang/scimldocs?branch=main) Attempts at a single documentation for the full SciML interface. - diff --git a/docs/make.jl b/docs/make.jl index b17f51ce158..04b614b9370 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,76 +8,74 @@ ENV["GKSwstype"] = "100" using Plots mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/require", "[tex]/mathtools"]), - :tex => Dict("inlineMath" => [["\$", "\$"], ["\\(", "\\)"]], - "packages" => [ - "base", - "ams", - "autoload", - "mathtools", - "require", - ]))) + :tex => Dict("inlineMath" => [["\$", "\$"], ["\\(", "\\)"]], + "packages" => [ + "base", + "ams", + "autoload", + "mathtools", + "require" + ]))) makedocs(sitename = "Overview of Julia's SciML", - authors = "Chris Rackauckas", - modules = Module[], - clean = true, doctest = false, linkcheck = true, - linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1477274812460449793", - "https://epubs.siam.org/doi/10.1137/0903023", - "https://bkamins.github.io/julialang/2020/12/24/minilanguage.html", - "https://arxiv.org/abs/2109.06786", - "https://arxiv.org/abs/2001.04385", - - ], - format = Documenter.HTML(assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/stable/", - mathengine = mathengine), - pages = [ - "SciML: Open Source Software for Scientific Machine Learning with Julia" => "index.md", - "Getting Started" => [ - "getting_started/getting_started.md", - "New User Tutorials" => [ - "getting_started/installation.md", - "getting_started/first_simulation.md", - "getting_started/first_optimization.md", - "getting_started/fit_simulation.md", - "getting_started/find_root.md", - ], - "Comparison With Other Tools" => [ - "comparisons/python.md", - "comparisons/matlab.md", - "comparisons/r.md", - "comparisons/cppfortran.md", - ], - ], - "Showcase of Cool Examples" => Any["showcase/showcase.md", - "Automated Model Discovery" => Any["showcase/missing_physics.md", - "showcase/bayesian_neural_ode.md", - "showcase/blackhole.md"], - "Solving Difficult Equations Efficiently" => Any["showcase/brusselator.md", - "showcase/pinngpu.md", - "showcase/massively_parallel_gpu.md", - "showcase/gpu_spde.md"], - "Useful Cool Wonky Things" => Any["showcase/ode_types.md", - "showcase/symbolic_analysis.md", - "showcase/optimization_under_uncertainty.md"]], - "What is SciML?" => ["overview.md", - "Solvers" => ["highlevels/equation_solvers.md", - "highlevels/inverse_problems.md", - "highlevels/partial_differential_equation_solvers.md"], - "Modeling Tools" => ["highlevels/modeling_languages.md", - "highlevels/model_libraries_and_importers.md", - "highlevels/symbolic_tools.md", - "highlevels/array_libraries.md"], - "Simulation Analysis" => ["highlevels/parameter_analysis.md", - "highlevels/uncertainty_quantification.md", - "highlevels/plots_visualization.md"], - "Machine Learning" => ["highlevels/function_approximation.md", - "highlevels/implicit_layers.md", - "highlevels/symbolic_learning.md"], - "Developer Tools" => ["highlevels/numerical_utilities.md", - "highlevels/interfaces.md", - "highlevels/developer_documentation.md"], - "Extra Learning Resources" => ["highlevels/learning_resources.md"], - ]]) + authors = "Chris Rackauckas", + modules = Module[], + clean = true, doctest = false, linkcheck = true, + linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1477274812460449793", + "https://epubs.siam.org/doi/10.1137/0903023", + "https://bkamins.github.io/julialang/2020/12/24/minilanguage.html", + "https://arxiv.org/abs/2109.06786", + "https://arxiv.org/abs/2001.04385"], + format = Documenter.HTML(assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/stable/", + mathengine = mathengine), + pages = [ + "SciML: Open Source Software for Scientific Machine Learning with Julia" => "index.md", + "Getting Started" => [ + "getting_started/getting_started.md", + "New User Tutorials" => [ + "getting_started/installation.md", + "getting_started/first_simulation.md", + "getting_started/first_optimization.md", + "getting_started/fit_simulation.md", + "getting_started/find_root.md" + ], + "Comparison With Other Tools" => [ + "comparisons/python.md", + "comparisons/matlab.md", + "comparisons/r.md", + "comparisons/cppfortran.md" + ] + ], + "Showcase of Cool Examples" => Any["showcase/showcase.md", + "Automated Model Discovery" => Any["showcase/missing_physics.md", + "showcase/bayesian_neural_ode.md", + "showcase/blackhole.md"], + "Solving Difficult Equations Efficiently" => Any["showcase/brusselator.md", + "showcase/pinngpu.md", + "showcase/massively_parallel_gpu.md", + "showcase/gpu_spde.md"], + "Useful Cool Wonky Things" => Any["showcase/ode_types.md", + "showcase/symbolic_analysis.md", + "showcase/optimization_under_uncertainty.md"]], + "What is SciML?" => ["overview.md", + "Solvers" => ["highlevels/equation_solvers.md", + "highlevels/inverse_problems.md", + "highlevels/partial_differential_equation_solvers.md"], + "Modeling Tools" => ["highlevels/modeling_languages.md", + "highlevels/model_libraries_and_importers.md", + "highlevels/symbolic_tools.md", + "highlevels/array_libraries.md"], + "Simulation Analysis" => ["highlevels/parameter_analysis.md", + "highlevels/uncertainty_quantification.md", + "highlevels/plots_visualization.md"], + "Machine Learning" => ["highlevels/function_approximation.md", + "highlevels/implicit_layers.md", + "highlevels/symbolic_learning.md"], + "Developer Tools" => ["highlevels/numerical_utilities.md", + "highlevels/interfaces.md", + "highlevels/developer_documentation.md"], + "Extra Learning Resources" => ["highlevels/learning_resources.md"] + ]]) deploydocs(repo = "github.com/SciML/SciMLDocs") diff --git a/docs/make_aggregate.jl b/docs/make_aggregate.jl index 7017a45bb34..5872636cdf3 100644 --- a/docs/make_aggregate.jl +++ b/docs/make_aggregate.jl @@ -28,8 +28,9 @@ docsmodules = [ =# ], "Solvers" => [ - "Equation Solvers" => ["LinearSolve", "NonlinearSolve", "DiffEqDocs", "Integrals", "DifferenceEquations", - "Optimization", "JumpProcesses"], + "Equation Solvers" => [ + "LinearSolve", "NonlinearSolve", "DiffEqDocs", "Integrals", + "DifferenceEquations", "Optimization", "JumpProcesses", "LineSearch"], #= "Third-Party Equation Solvers" => [ "LowRankIntegrators", diff --git a/docs/src/comparisons/matlab.md b/docs/src/comparisons/matlab.md index f57cc778d84..6bd73707f59 100644 --- a/docs/src/comparisons/matlab.md +++ b/docs/src/comparisons/matlab.md @@ -74,10 +74,10 @@ The following chart will help you get quickly acquainted with Julia's SciML Tool |:------------------- |:------------------------------------------------------------------------------------- | | `plot` | [Plots](https://docs.juliaplots.org/stable/), [Makie](https://docs.makie.org/stable/) | | `sparse` | [SparseArrays](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#Sparse-Arrays) | -| `interp1` | [DataInterpolations](https://docs.sciml.ai/DataInterpolations/) | -| `\`, `gmres`, `cg` | [LinearSolve](http://linearsolve.sciml.ai/dev/) | +| `interp1` | [DataInterpolations](https://docs.sciml.ai/DataInterpolations/) | +| `\`, `gmres`, `cg` | [LinearSolve](https://linearsolve.sciml.ai/dev/) | | `fsolve` | [NonlinearSolve](https://nonlinearsolve.sciml.ai/) | -| `quad` | [Integrals](https://docs.sciml.ai/Integrals/stable/) | +| `quad` | [Integrals](https://docs.sciml.ai/Integrals/stable/) | | `fmincon` | [Optimization](https://optimization.sciml.ai/) | | `odeXX` | [DifferentialEquations](https://diffeq.sciml.ai/latest/) | | `ode45` | `Tsit5` | diff --git a/docs/src/comparisons/python.md b/docs/src/comparisons/python.md index 706f855e167..9c527057bc9 100644 --- a/docs/src/comparisons/python.md +++ b/docs/src/comparisons/python.md @@ -61,7 +61,7 @@ The following chart will help you get quickly acquainted with Julia's SciML Tool |:---------------------------- |:----------------------------------------------------------------------------------------------------------------------------------------------------------------------- | | Matplotlib | [Plots](https://docs.juliaplots.org/stable/), [Makie](https://docs.makie.org/stable/) | | `scipy.special` | [SpecialFunctions](https://github.com/JuliaMath/SpecialFunctions.jl) | -| `scipy.linalg.solve` | [LinearSolve](http://linearsolve.sciml.ai/dev/) | +| `scipy.linalg.solve` | [LinearSolve](https://linearsolve.sciml.ai/dev/) | | `scipy.integrate` | [Integrals](https://docs.sciml.ai/Integrals/stable/) | | `scipy.optimize` | [Optimization](https://optimization.sciml.ai/) | | `scipy.optimize.fsolve` | [NonlinearSolve](https://nonlinearsolve.sciml.ai/) | @@ -71,7 +71,7 @@ The following chart will help you get quickly acquainted with Julia's SciML Tool | `scipy.sparse` | [SparseArrays](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#Sparse-Arrays), [ARPACK](https://github.com/JuliaLinearAlgebra/Arpack.jl) | | `odeint`/`solve_ivp` | [DifferentialEquations](https://diffeq.sciml.ai/latest/) | | `scipy.integrate.solve_bvp` | [Boundary-value problem](https://diffeq.sciml.ai/latest/tutorials/bvp_example/#Boundary-Value-Problems) | -| PyTorch | [Flux](https://fluxml.ai/), [Lux](http://lux.csail.mit.edu/stable/) | +| PyTorch | [Flux](https://fluxml.ai/), [Lux](https://lux.csail.mit.edu/stable/) | | gillespy2 | [Catalyst](https://catalyst.sciml.ai/dev/), [JumpProcesses](https://github.com/SciML/JumpProcesses.jl) | | scipy.optimize.approx_fprime | [FiniteDiff](https://github.com/JuliaDiff/FiniteDiff.jl) | | autograd | [ForwardDiff\*](https://github.com/JuliaDiff/ForwardDiff.jl), [Enzyme\*](https://github.com/EnzymeAD/Enzyme.jl), [DiffEqSensitivity](https://sensitivity.sciml.ai/dev/) | diff --git a/docs/src/getting_started/find_root.md b/docs/src/getting_started/find_root.md index 377ebe7a479..2f501229eaf 100644 --- a/docs/src/getting_started/find_root.md +++ b/docs/src/getting_started/find_root.md @@ -56,7 +56,7 @@ prob = NonlinearProblem(ns, []) sol = solve(prob, NewtonRaphson()) # Analyze the solution -@show sol[[x,y,z]], sol.resid +@show sol[[x, y, z]], sol.resid ``` ## Step-by-Step Solution @@ -178,5 +178,5 @@ We can check it as follows: ```@example first_rootfind # Analyze the solution -@show sol[[x,y,z]], sol.resid +@show sol[[x, y, z]], sol.resid ``` diff --git a/docs/src/getting_started/first_optimization.md b/docs/src/getting_started/first_optimization.md index 903dba81007..5012623a119 100644 --- a/docs/src/getting_started/first_optimization.md +++ b/docs/src/getting_started/first_optimization.md @@ -15,11 +15,11 @@ Let's solve our first optimization problem! The following parts of the SciML Ecosystem will be used in this tutorial: -| Module | Description | -|:---------------------------------------------------------------------------------------------- |:---------------------------------- | -| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The numerical optimization package | -| [OptimizationNLopt.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/nlopt/) | The NLopt optimizers we will use | -| [ForwardDiff.jl](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Optimization.AutoForwardDiff) | The automatic differentiation library for gradients | +| Module | Description | +|:------------------------------------------------------------------------------------------------------------------- |:--------------------------------------------------- | +| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The numerical optimization package | +| [OptimizationNLopt.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/nlopt/) | The NLopt optimizers we will use | +| [ForwardDiff.jl](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Optimization.AutoForwardDiff) | The automatic differentiation library for gradients | ## Problem Setup @@ -68,9 +68,9 @@ To do this tutorial, we will need a few components: - [Optimization.jl](https://docs.sciml.ai/Optimization/stable/), the optimization interface. - [OptimizationNLopt.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/nlopt/), the optimizers we will use. - - [ForwardDiff.jl](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Optimization.AutoForwardDiff), + - [ForwardDiff.jl](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Optimization.AutoForwardDiff), the automatic differentiation library for gradients - + Note that Optimization.jl is an interface for optimizers, and thus we always have to choose which optimizer we want to use. Here we choose to demonstrate `OptimizationNLopt` because of its efficiency and versatility. But there are many other possible choices. Check out @@ -102,6 +102,7 @@ parameters, and write out the loss function on a vector-defined state as follows # Define the problem to optimize L(u, p) = (p[1] - u[1])^2 + p[2] * (u[2] - u[1]^2)^2 ``` + Next we need to create an `OptimizationFunction` where we tell Optimization.jl to use the ForwardDiff.jl package for creating the gradient and other derivatives required by the optimizer. diff --git a/docs/src/getting_started/first_simulation.md b/docs/src/getting_started/first_simulation.md index 957d4672443..2a291e2aa01 100644 --- a/docs/src/getting_started/first_simulation.md +++ b/docs/src/getting_started/first_simulation.md @@ -173,9 +173,9 @@ to represent an `ODESystem` with the following: @mtkbuild sys = ODESystem(eqs, t) ``` -Notice that in our equations we have an algebraic equation `z ~ x + y`. This is not a -differential equation but an algebraic equation, and thus we call this set of equations a -Differential-Algebraic Equation (DAE). The symbolic system of ModelingToolkit can eliminate +Notice that in our equations we have an algebraic equation `z ~ x + y`. This is not a +differential equation but an algebraic equation, and thus we call this set of equations a +Differential-Algebraic Equation (DAE). The symbolic system of ModelingToolkit can eliminate such equations to return simpler forms to numerically approximate. Notice that what is returned is an `ODESystem`, but now with the simplified set of diff --git a/docs/src/getting_started/fit_simulation.md b/docs/src/getting_started/fit_simulation.md index 7929c7115c0..d5631c28f65 100644 --- a/docs/src/getting_started/fit_simulation.md +++ b/docs/src/getting_started/fit_simulation.md @@ -9,12 +9,12 @@ stable and efficient. Let's see this in action. The following parts of the SciML Ecosystem will be used in this tutorial: -| Module | Description | -|:---------------------------------------------------------------------------------------------------------------- |:--------------------------------------------------------- | -| [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) | The differential equation solvers | -| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The numerical optimization package | +| Module | Description | +|:------------------------------------------------------------------------------------------------------------------------------------------------------ |:--------------------------------------------------------- | +| [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) | The differential equation solvers | +| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The numerical optimization package | | [OptimizationPolyalgorithms.jl](https://github.com/SciML/Optimization.jl/blob/master/lib/OptimizationPolyalgorithms/src/OptimizationPolyalgorithms.jl) | The optimizers we will use | -| [SciMLSensitivity.jl](https://docs.sciml.ai/SciMLSensitivity/dev/) | The connection of the SciML ecosystems to differentiation | +| [SciMLSensitivity.jl](https://docs.sciml.ai/SciMLSensitivity/dev/) | The connection of the SciML ecosystems to differentiation | Along with the following general ecosystem packages: @@ -53,6 +53,7 @@ over a time span of ``t_0 = 0`` to ``t_f = 10`` at every ``\Delta t = 1``. Can we figure out what the parameter values should be directly from the data? ## Solution as Copy-Pastable Code + ```@example using DifferentialEquations, Optimization, OptimizationPolyalgorithms, SciMLSensitivity using ForwardDiff, Plots @@ -64,7 +65,7 @@ y_data = [1.000 0.259 2.015 1.908 0.323 0.629 3.458 0.508 0.314 4.547 0.906] xy_data = vcat(x_data, y_data) # Plot the provided data -scatter(t_data, xy_data', label=["x Data" "y Data"]) +scatter(t_data, xy_data', label = ["x Data" "y Data"]) # Setup the ODE function function lotka_volterra!(du, u, p, t) @@ -118,9 +119,9 @@ optprob = OptimizationProblem(optf, pguess) # Optimize the ODE parameters for best fit to our data pfinal = solve(optprob, PolyOpt(), - callback = callback, - maxiters = 200) -α, β, γ, δ = round.(pfinal, digits=1) + callback = callback, + maxiters = 200) +α, β, γ, δ = round.(pfinal, digits = 1) ``` ## Step-by-Step Solution @@ -132,13 +133,13 @@ To do this tutorial, we will need a few components. This is done using the Julia ```julia using Pkg Pkg.add([ - "DifferentialEquations", - "Optimization", - "OptimizationPolyalgorithms", - "SciMLSensitivity", - "ForwardDiff", - "Plots", - ]) + "DifferentialEquations", + "Optimization", + "OptimizationPolyalgorithms", + "SciMLSensitivity", + "ForwardDiff", + "Plots" +]) ``` Now we're ready. Let's load in these packages: @@ -162,11 +163,11 @@ y_data = [1.000 0.259 2.015 1.908 0.323 0.629 3.458 0.508 0.314 4.547 0.906] xy_data = vcat(x_data, y_data) # Plot the provided data -scatter(t_data, xy_data', label=["x Data" "y Data"]) +scatter(t_data, xy_data', label = ["x Data" "y Data"]) ``` !!! note - + The `Array` `xy_data` above has been oriented with time instances as columns so that it can be directly compared with an `ODESolution` object. (See [Solution Handling](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#solution) @@ -239,7 +240,7 @@ scatter!(plt, t_data, xy_data', label = ["x Data" "y Data"]) to fit our parameters to our training data. !!! note - + For more details on using DifferentialEquations.jl, check out the [getting started with DifferentialEquations.jl tutorial](https://docs.sciml.ai/DiffEqDocs/stable/getting_started/). @@ -263,7 +264,7 @@ For example, if we wanted to change the initial condition `u0` of our ODE, we co For our case, we want to change around just the parameters, so we can do `remake(prob, p = newp)`. It is faster to `remake` an existing `SciMLProblem` than to create a new problem every iteration. !!! note - + `remake` can change multiple items at once by passing more keyword arguments, i.e., `remake(prob, u0 = newu0, p = newp)`. This can be used to extend the example to simultaneously learn the initial conditions and parameters! @@ -347,54 +348,51 @@ optprob = OptimizationProblem(optf, pguess) # Optimize the ODE parameters for best fit to our data pfinal = solve(optprob, - PolyOpt(), - callback = callback, - maxiters = 200) -α, β, γ, δ = round.(pfinal, digits=1) + PolyOpt(), + callback = callback, + maxiters = 200) +α, β, γ, δ = round.(pfinal, digits = 1) ``` !!! note - + When referencing the documentation for DifferentialEquations.jl and Optimization.jl simultaneously, note that the variables `f`, `u`, and `p` will refer to different quantities. - + DifferentialEquations.jl: - + ```math \frac{du}{dt} = f(u,p,t) ``` - - - `f` in `ODEProblem` is the function defining the derivative `du` in the ODE. - + + - `f` in `ODEProblem` is the function defining the derivative `du` in the ODE. + Here: `lotka_volterra!` - - - `u` in `ODEProblem` contains the state variables of `f`. - + + - `u` in `ODEProblem` contains the state variables of `f`. + Here: `x` and `y` - - - `p` in `ODEProblem` contains the parameter variables of `f`. - + - `p` in `ODEProblem` contains the parameter variables of `f`. + Here: `\alpha`, `\beta`, `\gamma`, and `\delta` - - - `t` is the independent (time) variable. - + - `t` is the independent (time) variable. + Here: indirectly defined with `tspan` in `ODEProblem` and `saveat` in `solve` - + Optimization.jl: - + ```math \min_{u} f(u,p) ``` - - - `f` in `OptimizationProblem` is the function to minimize (optimize). - + + - `f` in `OptimizationProblem` is the function to minimize (optimize). + Here: the [anonymous function](https://docs.julialang.org/en/v1/manual/functions/#man-anonymous-functions) `(x, _) -> loss(x)` - - - `u` in `OptimizationProblem` contains the state variables of `f` to be optimized. - + + - `u` in `OptimizationProblem` contains the state variables of `f` to be optimized. + Here: the ODE parameters `\alpha`, `\beta`, `\gamma`, and `\delta` stored in `p` - - - `p` in `OptimizationProblem` contains any fixed - [hyperparameters](https://en.wikipedia.org/wiki/Hyperparameter_(machine_learning)) of `f`. - + - `p` in `OptimizationProblem` contains any fixed + [hyperparameters](https://en.wikipedia.org/wiki/Hyperparameter_(machine_learning)) of `f`. + Here: our `loss` function does not require any hyperparameters, so we pass `_` for this `p`. diff --git a/docs/src/getting_started/installation.md b/docs/src/getting_started/installation.md index b6b4c3dd16f..07e3c630d2b 100644 --- a/docs/src/getting_started/installation.md +++ b/docs/src/getting_started/installation.md @@ -26,7 +26,7 @@ You can run SciML with Julia in any development environment you please, but our environment is VS Code. For more information on using Julia with VS Code, check out the [Julia VS Code Extension website](https://www.julia-vscode.org/). Let's install it! -First [download VS Code from the official website](https://code.visualstudio.com/download). +First [download VS Code from the official website](https://code.visualstudio.com/). Next, open Visual Studio Code and click Extensions. diff --git a/docs/src/highlevels/learning_resources.md b/docs/src/highlevels/learning_resources.md index 71ed6f05a66..4f95119e2d7 100644 --- a/docs/src/highlevels/learning_resources.md +++ b/docs/src/highlevels/learning_resources.md @@ -33,6 +33,7 @@ classic SIR epidemic model. ## Other Books Featuring SciML - [Nonlinear Dynamics: A Concise Introduction Interlaced with Code](https://link.springer.com/book/10.1007/978-3-030-91032-7) + - [Numerical Methods for Scientific Computing: The Definitive Manual for Math Geeks](https://www.equalsharepress.com/) - [Fundamentals of Numerical Computation](https://tobydriscoll.net/fnc-julia/frontmatter.html) - [Statistics with Julia](https://statisticswithjulia.org/) diff --git a/docs/src/highlevels/model_libraries_and_importers.md b/docs/src/highlevels/model_libraries_and_importers.md index 1b61e2a8af2..4818fd5d6fb 100644 --- a/docs/src/highlevels/model_libraries_and_importers.md +++ b/docs/src/highlevels/model_libraries_and_importers.md @@ -29,7 +29,7 @@ quickly building up complex differential equation models. It includes: ## SBMLToolkit.jl: SBML Import [SBMLToolkit.jl](https://github.com/SciML/SBMLToolkit.jl) is a library for reading -[SBML files](https://synonym.caltech.edu/#:%7E:text=What%20is%20SBML%3F,field%20of%20the%20life%20sciences.) +[SBML files](https://sbml.org/) into the standard formats for Catalyst.jl and ModelingToolkit.jl. There are well over one thousand biological models available in the [BioModels Repository](https://www.ebi.ac.uk/biomodels/). diff --git a/docs/src/highlevels/modeling_languages.md b/docs/src/highlevels/modeling_languages.md index 2ac6d2c2c56..576759575b6 100644 --- a/docs/src/highlevels/modeling_languages.md +++ b/docs/src/highlevels/modeling_languages.md @@ -64,7 +64,7 @@ This image that went viral is actually runnable code from [ParameterizedFunction ## MomentClosure.jl: Automated Generation of Moment Closure Equations -[MomentClosure.jl](https://docs.sciml.ai/MomentClosure/dev/) is a library for generating the moment +[MomentClosure.jl](https://github.com/augustinas1/MomentClosure.jl) is a library for generating the moment closure equations for a given chemical master equation or stochastic differential equation. Thus instead of solving a stochastic model thousands of times to find the mean and variance, this library can generate the deterministic equations for how the mean and variance evolve in order to be solved in a single run. MomentClosure.jl @@ -74,18 +74,18 @@ processes. ## Agents.jl: Agent-Based Modeling Framework in Julia If one wants to do agent-based modeling in Julia, -[Agents.jl](https://docs.sciml.ai/Agents/stable/) is the go-to library. It's fast and flexible, +[Agents.jl](https://github.com/JuliaDynamics/Agents.jl) is the go-to library. It's fast and flexible, 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://docs.sciml.ai/Unitful/stable/) has minimal run-time penalty of units. +[Unitful.jl](https://painterqubits.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 -[ReactionMechanismSimulator.jl](https://docs.sciml.ai/ReactionMechanismSimulator/stable/) +[ReactionMechanismSimulator.jl](https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl) is a tool for simulating and analyzing large chemical reaction mechanisms. It interfaces with the ReactionMechanismGenerator suite for automatically constructing reaction pathways from chemical components to quickly build realistic models of chemical systems. @@ -99,7 +99,7 @@ allow for directly solving the weak form of the stochastic model. ## AlgebraicPetri.jl: Applied Category Theory of Modeling -[AlgebraicPetri.jl](https://docs.sciml.ai/AlgebraicPetri/stable/) is a library for automating the intuitive +[AlgebraicPetri.jl](https://github.com/AlgebraicJulia/AlgebraicPetri.jl) is a library for automating the intuitive generation of dynamical models using a Category theory-based approach. ## QuantumOptics.jl: Simulating quantum systems. diff --git a/docs/src/highlevels/numerical_utilities.md b/docs/src/highlevels/numerical_utilities.md index 7fec4992842..6771cae71b3 100644 --- a/docs/src/highlevels/numerical_utilities.md +++ b/docs/src/highlevels/numerical_utilities.md @@ -37,6 +37,7 @@ schemes which are composable with automatic differentiation and the SciML ecosys methods and regression techniques for handling noisy data. Its methods include: - `ConstantInterpolation(u,t)` - A piecewise constant interpolation. + - `LinearInterpolation(u,t)` - A linear interpolation. - `QuadraticInterpolation(u,t)` - A quadratic interpolation. - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`. @@ -44,7 +45,7 @@ methods and regression techniques for handling noisy data. Its methods include: - `CubicSpline(u,t)` - A cubic spline interpolation. - `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are: - + + `d` - degree of B-spline + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method. diff --git a/docs/src/index.md b/docs/src/index.md index ed17b8485ae..c99ba0f66de 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -78,4 +78,4 @@ file and the [project]($link_project) file. """) -``` \ No newline at end of file +``` diff --git a/docs/src/overview.md b/docs/src/overview.md index 3c3891fb435..45235dd942a 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -26,36 +26,36 @@ and much more. For an overview of the broad goals of the SciML organization, wat Below is a simplification of the user-facing packages for use in scientific computing and SciML workflows. -| Workflow Element | SciML-Supported Julia packages | -|:------------------------------------------------------------------------------------------------------- |:----------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| Plotting and Visualization | [Plots\*](https://docs.juliaplots.org/stable/), [Makie\*](https://docs.makie.org/stable/) | -| Sparse matrix | [SparseArrays\*](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#Sparse-Arrays) | -| Interpolation/approximation | [DataInterpolations\*](https://docs.sciml.ai/DataInterpolations/stable), [ApproxFun\*](https://juliaapproximation.github.io/ApproxFun.jl/stable/) | -| Linear system / least squares | [LinearSolve](https://docs.sciml.ai/LinearSolve/stable) | -| Nonlinear system / rootfinding | [NonlinearSolve](https://docs.sciml.ai/NonlinearSolve/stable/) | -| Polynomial roots | [Polynomials\*](https://juliamath.github.io/Polynomials.jl/stable/#Root-finding-1) | -| Integration | [Integrals](https://docs.sciml.ai/Integrals/stable/) | -| Nonlinear Optimization | [Optimization](https://docs.sciml.ai/Optimization/stable/) | -| Other Optimization (linear, quadratic, convex, etc.) | [JuMP\*](https://github.com/jump-dev/JuMP.jl) | -| [Initial-value problem](https://docs.sciml.ai/DiffEqDocs/latest/getting_started/) | [DifferentialEquations](https://docs.sciml.ai/DiffEqDocs/stable/) | -| [Boundary-value problem](https://docs.sciml.ai/DiffEqDocs/latest/tutorials/bvp_example/) | [DifferentialEquations](https://docs.sciml.ai/DifferenceEquations/stable/) | -| Continuous-Time Markov Chains (Poisson Jumps), Jump Diffusions | [JumpProcesses](https://docs.sciml.ai/JumpProcesses/stable/) | -| Finite differences | [FiniteDifferences\*](https://juliadiff.org/FiniteDifferences.jl/latest/), [FiniteDiff\*](https://github.com/JuliaDiff/FiniteDiff.jl) | -| Automatic Differentiation | [ForwardDiff\*](https://github.com/JuliaDiff/ForwardDiff.jl), [Enzyme\*](https://github.com/EnzymeAD/Enzyme.jl), [DiffEqSensitivity](https://docs.sciml.ai/SciMLSensitivity/stable/) | -| Bayesian Modeling | [Turing\*](https://turinglang.org/docs/) | -| Deep Learning | [Flux\*](https://fluxml.ai/) | -| Acausal Modeling / DAEs | [ModelingToolkit](https://docs.sciml.ai/ModelingToolkit/stable/) | -| Chemical Reaction Networks | [Catalyst](https://docs.sciml.ai/Catalyst/stable/) | -| Symbolic Computing | [Symbolics](https://docs.sciml.ai/Symbolics/stable/) | -| Fast Fourier Transform | [FFTW\*](https://github.com/JuliaMath/FFTW.jl) | -| Partial Differential Equation Discretizations | Associated Julia packages | -| --- | --- | -| Finite Differences | [MethodOfLines](https://docs.sciml.ai/MethodOfLines/stable/) | -| Discontinuous Galerkin | [Trixi\*](https://github.com/trixi-framework/Trixi.jl) | -| Finite Element | [Gridap\*](https://github.com/gridap/Gridap.jl) | -| Physics-Informed Neural Networks | [NeuralPDE](https://docs.sciml.ai/NeuralPDE/stable/) | -| Neural Operators | [NeuralOperators](https://docs.sciml.ai/NeuralOperators/stable/) | -| High Dimensional Deep Learning | [HighDimPDE](https://docs.sciml.ai/HighDimPDE/stable/) | +| Workflow Element | SciML-Supported Julia packages | +|:---------------------------------------------------------------------------------------- |:------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | +| Plotting and Visualization | [Plots\*](https://docs.juliaplots.org/stable/), [Makie\*](https://docs.makie.org/stable/) | +| Sparse matrix | [SparseArrays\*](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#Sparse-Arrays) | +| Interpolation/approximation | [DataInterpolations\*](https://docs.sciml.ai/DataInterpolations/stable), [ApproxFun\*](https://juliaapproximation.github.io/ApproxFun.jl/stable/) | +| Linear system / least squares | [LinearSolve](https://docs.sciml.ai/LinearSolve/stable) | +| Nonlinear system / rootfinding | [NonlinearSolve](https://docs.sciml.ai/NonlinearSolve/stable/) | +| Polynomial roots | [Polynomials\*](https://juliamath.github.io/Polynomials.jl/stable/#Root-finding-1) | +| Integration | [Integrals](https://docs.sciml.ai/Integrals/stable/) | +| Nonlinear Optimization | [Optimization](https://docs.sciml.ai/Optimization/stable/) | +| Other Optimization (linear, quadratic, convex, etc.) | [JuMP\*](https://github.com/jump-dev/JuMP.jl) | +| [Initial-value problem](https://docs.sciml.ai/DiffEqDocs/latest/getting_started/) | [DifferentialEquations](https://docs.sciml.ai/DiffEqDocs/stable/) | +| [Boundary-value problem](https://docs.sciml.ai/DiffEqDocs/latest/tutorials/bvp_example/) | [DifferentialEquations](https://docs.sciml.ai/DifferenceEquations/stable/) | +| Continuous-Time Markov Chains (Poisson Jumps), Jump Diffusions | [JumpProcesses](https://docs.sciml.ai/JumpProcesses/stable/) | +| Finite differences | [FiniteDifferences\*](https://juliadiff.org/FiniteDifferences.jl/latest/), [FiniteDiff\*](https://github.com/JuliaDiff/FiniteDiff.jl) | +| Automatic Differentiation | [ForwardDiff\*](https://github.com/JuliaDiff/ForwardDiff.jl), [Enzyme\*](https://github.com/EnzymeAD/Enzyme.jl), [DiffEqSensitivity](https://docs.sciml.ai/SciMLSensitivity/stable/) | +| Bayesian Modeling | [Turing\*](https://turinglang.org/docs/) | +| Deep Learning | [Flux\*](https://fluxml.ai/) | +| Acausal Modeling / DAEs | [ModelingToolkit](https://docs.sciml.ai/ModelingToolkit/stable/) | +| Chemical Reaction Networks | [Catalyst](https://docs.sciml.ai/Catalyst/stable/) | +| Symbolic Computing | [Symbolics](https://docs.sciml.ai/Symbolics/stable/) | +| Fast Fourier Transform | [FFTW\*](https://github.com/JuliaMath/FFTW.jl) | +| Partial Differential Equation Discretizations | Associated Julia packages | +| --- | --- | +| Finite Differences | [MethodOfLines](https://docs.sciml.ai/MethodOfLines/stable/) | +| Discontinuous Galerkin | [Trixi\*](https://github.com/trixi-framework/Trixi.jl) | +| Finite Element | [Gridap\*](https://github.com/gridap/Gridap.jl) | +| Physics-Informed Neural Networks | [NeuralPDE](https://docs.sciml.ai/NeuralPDE/stable/) | +| Neural Operators | [NeuralOperators](https://docs.sciml.ai/NeuralOperators/stable/) | +| High Dimensional Deep Learning | [HighDimPDE](https://docs.sciml.ai/HighDimPDE/stable/) | \* Denotes a non-SciML package that is heavily tested against as part of SciML workflows and has frequent collaboration with the SciML developers. diff --git a/docs/src/showcase/bayesian_neural_ode.md b/docs/src/showcase/bayesian_neural_ode.md index 87cb997fbb2..54a9ab0c1ac 100644 --- a/docs/src/showcase/bayesian_neural_ode.md +++ b/docs/src/showcase/bayesian_neural_ode.md @@ -1,174 +1,174 @@ -# [Uncertainty Quantified Deep Bayesian Model Discovery](@id bnode) - -In this tutorial, we show how SciML can combine the differential equation solvers seamlessly -with Bayesian estimation libraries like AdvancedHMC.jl and Turing.jl. This enables -converting Neural ODEs to Bayesian Neural ODEs, which enables us to estimate the error in -the Neural ODE estimation and forecasting. In this tutorial, a working example of the -Bayesian Neural ODE: NUTS sampler is shown. - -!!! note - - For more details, have a look at this paper: https://arxiv.org/abs/2012.07244 - -## Step 1: Import Libraries - -For this example, we will need the following libraries: - -```@example bnode -# SciML Libraries -using SciMLSensitivity, DifferentialEquations - -# ML Tools -using Lux, Zygote - -# External Tools -using Random, Plots, AdvancedHMC, MCMCChains, StatsPlots, ComponentArrays -``` - -## Setup: Get the data from the Spiral ODE example - -We will also need data to fit against. As a demonstration, we will generate our data -using a simple cubic ODE `u' = A*u^3` as follows: - -```@example bnode -u0 = [2.0; 0.0] -datasize = 40 -tspan = (0.0, 1) -tsteps = range(tspan[1], tspan[2], length = datasize) -function trueODEfunc(du, u, p, t) - true_A = [-0.1 2.0; -2.0 -0.1] - du .= ((u .^ 3)'true_A)' -end -prob_trueode = ODEProblem(trueODEfunc, u0, tspan) -ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps)) -``` - -We will want to train a neural network to capture the dynamics that fit `ode_data`. - -## Step 2: Define the Neural ODE architecture. - -Note that this step potentially offers a lot of flexibility in the number of layers/ number -of units in each layer. It may not necessarily be true that a 100 units architecture is -better at prediction/forecasting than a 50 unit architecture. On the other hand, a -complicated architecture can take a huge computational time without increasing performance. - -```@example bnode -dudt2 = Lux.Chain(x -> x .^ 3, - Lux.Dense(2, 50, tanh), - Lux.Dense(50, 2)) - -rng = Random.default_rng() -p, st = Lux.setup(rng, dudt2) -const _st = st -function neuralodefunc(u, p, t) - dudt2(u, p, _st)[1] -end -function prob_neuralode(u0, p) - prob = ODEProblem(neuralodefunc, u0, tspan, p) - sol = solve(prob, Tsit5(), saveat = tsteps) -end -p = ComponentArray{Float64}(p) -const _p = p -``` - -Note that the `f64` is required to put the Lux neural network into Float64 precision. - -## Step 3: Define the loss function for the Neural ODE. - -```@example bnode -function predict_neuralode(p) - p = p isa ComponentArray ? p : convert(typeof(_p),p) - Array(prob_neuralode(u0, p)) -end -function loss_neuralode(p) - pred = predict_neuralode(p) - loss = sum(abs2, ode_data .- pred) - return loss, pred -end -``` - -## Step 4: Now we start integrating the Bayesian estimation workflow as prescribed by the AdvancedHMC interface with the NeuralODE defined above - -The AdvancedHMC interface requires us to specify: (a) the Hamiltonian log density and its gradient , (b) the sampler and (c) the step size adaptor function. - -For the Hamiltonian log density, we use the loss function. The θ*θ term denotes the use of Gaussian priors. - -The user can make several modifications to Step 4. The user can try different acceptance ratios, warmup samples and posterior samples. One can also use the Variational Inference (ADVI) framework, which doesn't work quite as well as NUTS. The SGLD (Stochastic Gradient Langevin Descent) sampler is seen to have a better performance than NUTS. Have a look at https://sebastiancallh.github.io/post/langevin/ for a brief introduction to SGLD. - -```@example bnode -l(θ) = -sum(abs2, ode_data .- predict_neuralode(θ)) - sum(θ .* θ) -function dldθ(θ) - x, lambda = Zygote.pullback(l, θ) - grad = first(lambda(1)) - return x, grad -end - -metric = DiagEuclideanMetric(length(p)) -h = Hamiltonian(metric, l, dldθ) -``` - -We use the NUTS sampler with an acceptance ratio of δ= 0.45 in this example. In addition, we use Nesterov Dual Averaging for the Step Size adaptation. - -We sample using 500 warmup samples and 500 posterior samples. - -```@example bnode -integrator = Leapfrog(find_good_stepsize(h, p)) -kernel = HMCKernel(Trajectory{MultinomialTS}(integrator, GeneralisedNoUTurn())) -adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.45, integrator)) -samples, stats = sample(h, kernel, p, 500, adaptor, 500; progress = true) -``` - -## Step 5: Plot diagnostics - -Now let's make sure the fit is good. This can be done by looking at the chain mixing plot -and the autocorrelation plot. First, let's create the chain mixing plot using the plot -recipes from ???? - -```@example bnode -samples = hcat(samples...) -samples_reduced = samples[1:5, :] -samples_reshape = reshape(samples_reduced, (500, 5, 1)) -Chain_Spiral = Chains(samples_reshape) -plot(Chain_Spiral) -``` - -Now we check the autocorrelation plot: - -```@example bnode -autocorplot(Chain_Spiral) -``` - -As another diagnostic, let's check the result on retrodicted data. To do this, we generate -solutions of the Neural ODE on samples of the neural network parameters, and check the -results of the predictions against the data. Let's start by looking at the time series: - -```@example bnode -pl = scatter(tsteps, ode_data[1, :], color = :red, label = "Data: Var1", xlabel = "t", - title = "Spiral Neural ODE") -scatter!(tsteps, ode_data[2, :], color = :blue, label = "Data: Var2") -for k in 1:300 - resol = predict_neuralode(samples[:, 100:end][:, rand(1:400)]) - plot!(tsteps, resol[1, :], alpha = 0.04, color = :red, label = "") - plot!(tsteps, resol[2, :], alpha = 0.04, color = :blue, label = "") -end - -losses = map(x -> loss_neuralode(x)[1], eachcol(samples)) -idx = findmin(losses)[2] -prediction = predict_neuralode(samples[:, idx]) -plot!(tsteps, prediction[1, :], color = :black, w = 2, label = "") -plot!(tsteps, prediction[2, :], color = :black, w = 2, label = "Best fit prediction", - ylims = (-2.5, 3.5)) -``` - -That showed the time series form. We can similarly do a phase-space plot: - -```@example bnode -pl = scatter(ode_data[1, :], ode_data[2, :], color = :red, label = "Data", xlabel = "Var1", - ylabel = "Var2", title = "Spiral Neural ODE") -for k in 1:300 - resol = predict_neuralode(samples[:, 100:end][:, rand(1:400)]) - plot!(resol[1, :], resol[2, :], alpha = 0.04, color = :red, label = "") -end -plot!(prediction[1, :], prediction[2, :], color = :black, w = 2, - label = "Best fit prediction", ylims = (-2.5, 3)) -``` +# [Uncertainty Quantified Deep Bayesian Model Discovery](@id bnode) + +In this tutorial, we show how SciML can combine the differential equation solvers seamlessly +with Bayesian estimation libraries like AdvancedHMC.jl and Turing.jl. This enables +converting Neural ODEs to Bayesian Neural ODEs, which enables us to estimate the error in +the Neural ODE estimation and forecasting. In this tutorial, a working example of the +Bayesian Neural ODE: NUTS sampler is shown. + +!!! note + + For more details, have a look at this paper: https://arxiv.org/abs/2012.07244 + +## Step 1: Import Libraries + +For this example, we will need the following libraries: + +```@example bnode +# SciML Libraries +using SciMLSensitivity, DifferentialEquations + +# ML Tools +using Lux, Zygote + +# External Tools +using Random, Plots, AdvancedHMC, MCMCChains, StatsPlots, ComponentArrays +``` + +## Setup: Get the data from the Spiral ODE example + +We will also need data to fit against. As a demonstration, we will generate our data +using a simple cubic ODE `u' = A*u^3` as follows: + +```@example bnode +u0 = [2.0; 0.0] +datasize = 40 +tspan = (0.0, 1) +tsteps = range(tspan[1], tspan[2], length = datasize) +function trueODEfunc(du, u, p, t) + true_A = [-0.1 2.0; -2.0 -0.1] + du .= ((u .^ 3)'true_A)' +end +prob_trueode = ODEProblem(trueODEfunc, u0, tspan) +ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps)) +``` + +We will want to train a neural network to capture the dynamics that fit `ode_data`. + +## Step 2: Define the Neural ODE architecture. + +Note that this step potentially offers a lot of flexibility in the number of layers/ number +of units in each layer. It may not necessarily be true that a 100 units architecture is +better at prediction/forecasting than a 50 unit architecture. On the other hand, a +complicated architecture can take a huge computational time without increasing performance. + +```@example bnode +dudt2 = Lux.Chain(x -> x .^ 3, + Lux.Dense(2, 50, tanh), + Lux.Dense(50, 2)) + +rng = Random.default_rng() +p, st = Lux.setup(rng, dudt2) +const _st = st +function neuralodefunc(u, p, t) + dudt2(u, p, _st)[1] +end +function prob_neuralode(u0, p) + prob = ODEProblem(neuralodefunc, u0, tspan, p) + sol = solve(prob, Tsit5(), saveat = tsteps) +end +p = ComponentArray{Float64}(p) +const _p = p +``` + +Note that the `f64` is required to put the Lux neural network into Float64 precision. + +## Step 3: Define the loss function for the Neural ODE. + +```@example bnode +function predict_neuralode(p) + p = p isa ComponentArray ? p : convert(typeof(_p), p) + Array(prob_neuralode(u0, p)) +end +function loss_neuralode(p) + pred = predict_neuralode(p) + loss = sum(abs2, ode_data .- pred) + return loss, pred +end +``` + +## Step 4: Now we start integrating the Bayesian estimation workflow as prescribed by the AdvancedHMC interface with the NeuralODE defined above + +The AdvancedHMC interface requires us to specify: (a) the Hamiltonian log density and its gradient , (b) the sampler and (c) the step size adaptor function. + +For the Hamiltonian log density, we use the loss function. The θ*θ term denotes the use of Gaussian priors. + +The user can make several modifications to Step 4. The user can try different acceptance ratios, warmup samples and posterior samples. One can also use the Variational Inference (ADVI) framework, which doesn't work quite as well as NUTS. The SGLD (Stochastic Gradient Langevin Descent) sampler is seen to have a better performance than NUTS. Have a look at https://sebastiancallh.github.io/post/langevin/ for a brief introduction to SGLD. + +```@example bnode +l(θ) = -sum(abs2, ode_data .- predict_neuralode(θ)) - sum(θ .* θ) +function dldθ(θ) + x, lambda = Zygote.pullback(l, θ) + grad = first(lambda(1)) + return x, grad +end + +metric = DiagEuclideanMetric(length(p)) +h = Hamiltonian(metric, l, dldθ) +``` + +We use the NUTS sampler with an acceptance ratio of δ= 0.45 in this example. In addition, we use Nesterov Dual Averaging for the Step Size adaptation. + +We sample using 500 warmup samples and 500 posterior samples. + +```@example bnode +integrator = Leapfrog(find_good_stepsize(h, p)) +kernel = HMCKernel(Trajectory{MultinomialTS}(integrator, GeneralisedNoUTurn())) +adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.45, integrator)) +samples, stats = sample(h, kernel, p, 500, adaptor, 500; progress = true) +``` + +## Step 5: Plot diagnostics + +Now let's make sure the fit is good. This can be done by looking at the chain mixing plot +and the autocorrelation plot. First, let's create the chain mixing plot using the plot +recipes from ???? + +```@example bnode +samples = hcat(samples...) +samples_reduced = samples[1:5, :] +samples_reshape = reshape(samples_reduced, (500, 5, 1)) +Chain_Spiral = Chains(samples_reshape) +plot(Chain_Spiral) +``` + +Now we check the autocorrelation plot: + +```@example bnode +autocorplot(Chain_Spiral) +``` + +As another diagnostic, let's check the result on retrodicted data. To do this, we generate +solutions of the Neural ODE on samples of the neural network parameters, and check the +results of the predictions against the data. Let's start by looking at the time series: + +```@example bnode +pl = scatter(tsteps, ode_data[1, :], color = :red, label = "Data: Var1", xlabel = "t", + title = "Spiral Neural ODE") +scatter!(tsteps, ode_data[2, :], color = :blue, label = "Data: Var2") +for k in 1:300 + resol = predict_neuralode(samples[:, 100:end][:, rand(1:400)]) + plot!(tsteps, resol[1, :], alpha = 0.04, color = :red, label = "") + plot!(tsteps, resol[2, :], alpha = 0.04, color = :blue, label = "") +end + +losses = map(x -> loss_neuralode(x)[1], eachcol(samples)) +idx = findmin(losses)[2] +prediction = predict_neuralode(samples[:, idx]) +plot!(tsteps, prediction[1, :], color = :black, w = 2, label = "") +plot!(tsteps, prediction[2, :], color = :black, w = 2, label = "Best fit prediction", + ylims = (-2.5, 3.5)) +``` + +That showed the time series form. We can similarly do a phase-space plot: + +```@example bnode +pl = scatter(ode_data[1, :], ode_data[2, :], color = :red, label = "Data", xlabel = "Var1", + ylabel = "Var2", title = "Spiral Neural ODE") +for k in 1:300 + resol = predict_neuralode(samples[:, 100:end][:, rand(1:400)]) + plot!(resol[1, :], resol[2, :], alpha = 0.04, color = :red, label = "") +end +plot!(prediction[1, :], prediction[2, :], color = :black, w = 2, + label = "Best fit prediction", ylims = (-2.5, 3)) +``` diff --git a/docs/src/showcase/blackhole.md b/docs/src/showcase/blackhole.md index b2d1225af0a..d3b90dc892d 100644 --- a/docs/src/showcase/blackhole.md +++ b/docs/src/showcase/blackhole.md @@ -39,7 +39,7 @@ And external libraries: | Module | Description | |:---------------------------------------------------------------------------- |:--------------------------------------------------- | -| [Lux.jl](http://lux.csail.mit.edu/stable/) | The deep learning (neural network) framework | +| [Lux.jl](https://lux.csail.mit.edu/stable/) | The deep learning (neural network) framework | | [ComponentArrays.jl](https://jonniedie.github.io/ComponentArrays.jl/stable/) | For the `ComponentArray` type to match Lux to SciML | | [LineSearches.jl](https://github.com/JuliaNLSolvers/LineSearches.jl) | Allows for setting a line search for optimization | | [DataFrames.jl](https://dataframes.juliadata.org/stable/) | A nice and easy data handling format | @@ -100,7 +100,6 @@ function NewtonianOrbitModel(u, model_params, t) ϕ̇ = numer / denom return [χ̇, ϕ̇] - end function RelativisticOrbitModel(u, model_params, t) @@ -123,10 +122,9 @@ function RelativisticOrbitModel(u, model_params, t) ϕ̇ = numer / (M * (p^(3 / 2)) * denom) return [χ̇, ϕ̇] - end -function AbstractNNOrbitModel(u, model_params, t; NN=nothing, NN_params=nothing) +function AbstractNNOrbitModel(u, model_params, t; NN = nothing, NN_params = nothing) #= Defines system of odes which describes motion of point like particle with Newtonian physics, uses @@ -152,12 +150,11 @@ function AbstractNNOrbitModel(u, model_params, t; NN=nothing, NN_params=nothing) ϕ̇ = (numer / denom) * nn[2] return [χ̇, ϕ̇] - end function AbstractNROrbitModel(u, model_params, t; - NN_chiphi=nothing, NN_chiphi_params=nothing, - NN_pe=nothing, NN_pe_params=nothing) + NN_chiphi = nothing, NN_chiphi_params = nothing, + NN_pe = nothing, NN_pe_params = nothing) #= Defines system of odes which describes motion of point like particle with Newtonian physics, uses @@ -205,7 +202,7 @@ end =# using DelimitedFiles -function soln2orbit(soln, model_params=nothing) +function soln2orbit(soln, model_params = nothing) #= Performs change of variables: (χ(t),ϕ(t)) ↦ (x(t),y(t)) @@ -235,7 +232,7 @@ function soln2orbit(soln, model_params=nothing) return orbit end -function orbit2tensor(orbit, component, mass=1.0) +function orbit2tensor(orbit, component, mass = 1.0) #= Construct trace-free moment tensor Ι(t) for orbit from BH orbit (x(t),y(t)) @@ -264,20 +261,20 @@ end function d_dt(v::AbstractVector, dt) # uses second-order one-sided difference stencils at the endpoints; see https://doi.org/10.1090/S0025-5718-1988-0935077-0 a = -3 / 2 * v[1] + 2 * v[2] - 1 / 2 * v[3] - b = (v[3:end] .- v[1:end-2]) / 2 - c = 3 / 2 * v[end] - 2 * v[end-1] + 1 / 2 * v[end-2] + b = (v[3:end] .- v[1:(end - 2)]) / 2 + c = 3 / 2 * v[end] - 2 * v[end - 1] + 1 / 2 * v[end - 2] return [a; b; c] / dt end function d2_dt2(v::AbstractVector, dt) # uses second-order one-sided difference stencils at the endpoints; see https://doi.org/10.1090/S0025-5718-1988-0935077-0 a = 2 * v[1] - 5 * v[2] + 4 * v[3] - v[4] - b = v[1:end-2] .- 2 * v[2:end-1] .+ v[3:end] - c = 2 * v[end] - 5 * v[end-1] + 4 * v[end-2] - v[end-3] + b = v[1:(end - 2)] .- 2 * v[2:(end - 1)] .+ v[3:end] + c = 2 * v[end] - 5 * v[end - 1] + 4 * v[end - 2] - v[end - 3] return [a; b; c] / (dt^2) end -function h_22_quadrupole_components(dt, orbit, component, mass=1.0) +function h_22_quadrupole_components(dt, orbit, component, mass = 1.0) #= x(t) and y(t) inputs are the trajectory of the orbiting BH. @@ -292,7 +289,7 @@ function h_22_quadrupole_components(dt, orbit, component, mass=1.0) return 2 * mtensor_ddot end -function h_22_quadrupole(dt, orbit, mass=1.0) +function h_22_quadrupole(dt, orbit, mass = 1.0) h11 = h_22_quadrupole_components(dt, orbit, (1, 1), mass) h22 = h_22_quadrupole_components(dt, orbit, (2, 2), mass) h12 = h_22_quadrupole_components(dt, orbit, (1, 2), mass) @@ -300,7 +297,6 @@ function h_22_quadrupole(dt, orbit, mass=1.0) end function h_22_strain_one_body(dt, orbit) - h11, h12, h22 = h_22_quadrupole(dt, orbit) h₊ = h11 - h22 @@ -322,7 +318,7 @@ end function h_22_strain_two_body(dt, orbit1, mass1, orbit2, mass2) # compute (2,2) mode strain from orbits of BH 1 of mass1 and BH2 of mass 2 - @assert abs(mass1 + mass2 - 1.0) < 1e-12 "Masses do not sum to unity" + @assert abs(mass1 + mass2 - 1.0)<1e-12 "Masses do not sum to unity" h11, h12, h22 = h_22_quadrupole_two_body(dt, orbit1, mass1, orbit2, mass2) @@ -349,10 +345,9 @@ function one2two(path, m1, m2) return r1, r2 end -function compute_waveform(dt, soln, mass_ratio, model_params=nothing) - - @assert mass_ratio <= 1.0 "mass_ratio must be <= 1" - @assert mass_ratio >= 0.0 "mass_ratio must be non-negative" +function compute_waveform(dt, soln, mass_ratio, model_params = nothing) + @assert mass_ratio<=1.0 "mass_ratio must be <= 1" + @assert mass_ratio>=0.0 "mass_ratio must be non-negative" orbit = soln2orbit(soln, model_params) if mass_ratio > 0 @@ -368,15 +363,14 @@ function compute_waveform(dt, soln, mass_ratio, model_params=nothing) end function interpolate_time_series(tsteps, tdata, fdata) - - @assert length(tdata) == length(fdata) "lengths of tdata and fdata must match" + @assert length(tdata)==length(fdata) "lengths of tdata and fdata must match" interp_fdata = zeros(length(tsteps)) - for j = 1:length(tsteps) - for i = 1:length(tdata)-1 - if tdata[i] <= tsteps[j] < tdata[i+1] - weight = (tsteps[j] - tdata[i]) / (tdata[i+1] - tdata[i]) - interp_fdata[j] = (1 - weight) * fdata[i] + weight * fdata[i+1] + for j in 1:length(tsteps) + for i in 1:(length(tdata) - 1) + if tdata[i] <= tsteps[j] < tdata[i + 1] + weight = (tsteps[j] - tdata[i]) / (tdata[i + 1] - tdata[i]) + interp_fdata[j] = (1 - weight) * fdata[i] + weight * fdata[i + 1] break end end @@ -385,7 +379,7 @@ function interpolate_time_series(tsteps, tdata, fdata) return interp_fdata end -function file2waveform(tsteps, filename="waveform.txt") +function file2waveform(tsteps, filename = "waveform.txt") # read in file f = open(filename, "r") @@ -399,7 +393,7 @@ function file2waveform(tsteps, filename="waveform.txt") return waveform end -function file2trajectory(tsteps, filename="trajectoryA.txt") +function file2trajectory(tsteps, filename = "trajectoryA.txt") # read in file f = open(filename, "r") @@ -430,7 +424,7 @@ u0 = Float64[pi, 0.0] # initial conditions datasize = 250 tspan = (0.0f0, 6.0f4) # timespace for GW waveform tsteps = range(tspan[1], tspan[2], length = datasize) # time at each timestep -dt_data = tsteps[2] - tsteps[1] +dt_data = tsteps[2] - tsteps[1] dt = 100.0 model_params = [100.0, 1.0, 0.5]; # p, M, e ``` @@ -439,15 +433,15 @@ and demonstrate the gravitational waveform: ```@example ude prob = ODEProblem(RelativisticOrbitModel, u0, tspan, model_params) -soln = Array(solve(prob, RK4(), saveat = tsteps, dt = dt, adaptive=false)) +soln = Array(solve(prob, RK4(), saveat = tsteps, dt = dt, adaptive = false)) waveform = compute_waveform(dt_data, soln, mass_ratio, model_params)[1] plt = plot(tsteps, waveform, - markershape=:circle, markeralpha = 0.25, - linewidth = 2, alpha = 0.5, - label="waveform data", xlabel="Time", ylabel="Waveform") + markershape = :circle, markeralpha = 0.25, + linewidth = 2, alpha = 0.5, + label = "waveform data", xlabel = "Time", ylabel = "Waveform") ``` -Looks great! +Looks great! ## Automating the Discovery of Relativistic Equations from Newtonian Physics @@ -463,25 +457,26 @@ p, st = Lux.setup(rng, NN) NN_params = ComponentArray{Float64}(p) function ODE_model(u, NN_params, t) - du = AbstractNNOrbitModel(u, model_params, t, NN=NN, NN_params=NN_params) + du = AbstractNNOrbitModel(u, model_params, t, NN = NN, NN_params = NN_params) return du end ``` -Next, we can compute the orbital trajectory and gravitational waveform using the neural network with its initial weights. +Next, we can compute the orbital trajectory and gravitational waveform using the neural network with its initial weights. ```@example ude prob_nn = ODEProblem(ODE_model, u0, tspan, NN_params) -soln_nn = Array(solve(prob_nn, RK4(), u0 = u0, p = NN_params, saveat = tsteps, dt = dt, adaptive=false)) +soln_nn = Array(solve( + prob_nn, RK4(), u0 = u0, p = NN_params, saveat = tsteps, dt = dt, adaptive = false)) waveform_nn = compute_waveform(dt_data, soln_nn, mass_ratio, model_params)[1] plot!(plt, tsteps, waveform_nn, - markershape=:circle, markeralpha = 0.25, - linewidth = 2, alpha = 0.5, - label="waveform NN") + markershape = :circle, markeralpha = 0.25, + linewidth = 2, alpha = 0.5, + label = "waveform NN") display(plt) ``` -This is the model before training. +This is the model before training. Next, we define the objective (loss) function to be minimized when training the neural differential equations. @@ -489,12 +484,15 @@ Next, we define the objective (loss) function to be minimized when training the function loss(NN_params) first_obs_to_use_for_training = 1 last_obs_to_use_for_training = length(waveform) - obs_to_use_for_training = first_obs_to_use_for_training:last_obs_to_use_for_training - - pred = Array(solve(prob_nn, RK4(), u0 = u0, p = NN_params, saveat = tsteps, dt = dt, adaptive=false)) + obs_to_use_for_training = first_obs_to_use_for_training:last_obs_to_use_for_training + + pred = Array(solve( + prob_nn, RK4(), u0 = u0, p = NN_params, saveat = tsteps, dt = dt, adaptive = false)) pred_waveform = compute_waveform(dt_data, pred, mass_ratio, model_params)[1] - loss = ( sum(abs2, view(waveform,obs_to_use_for_training) .- view(pred_waveform,obs_to_use_for_training) ) ) + loss = (sum(abs2, + view(waveform, obs_to_use_for_training) .- + view(pred_waveform, obs_to_use_for_training))) return loss, pred_waveform end ``` @@ -510,7 +508,7 @@ We'll use the following callback to save the history of the loss values. ```@example ude losses = [] -callback(θ,l,pred_waveform; doplot = true) = begin +callback(θ, l, pred_waveform; doplot = true) = begin push!(losses, l) #= Disable plotting as it trains since in docs display(l) @@ -539,30 +537,36 @@ The next cell initializes the weights of the neural network and then trains the Training uses the BFGS optimizers. This seems to give good results because the Newtonian model seems to give a very good initial guess. ```@example ude -NN_params = NN_params .* 0 + Float64(1e-4) * randn(StableRNG(2031), eltype(NN_params), size(NN_params)) +NN_params = NN_params .* 0 + + Float64(1e-4) * randn(StableRNG(2031), eltype(NN_params), size(NN_params)) adtype = Optimization.AutoZygote() optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype) optprob = Optimization.OptimizationProblem(optf, ComponentVector{Float64}(NN_params)) -res1 = Optimization.solve(optprob, OptimizationOptimisers.Adam(0.001f0), callback=callback, maxiters=100) +res1 = Optimization.solve( + optprob, OptimizationOptimisers.Adam(0.001f0), callback = callback, maxiters = 100) optprob = Optimization.OptimizationProblem(optf, res1.u) -res2 = Optimization.solve(optprob, BFGS(initial_stepnorm=0.01, linesearch=LineSearches.BackTracking()), callback=callback, maxiters=20) +res2 = Optimization.solve( + optprob, BFGS(initial_stepnorm = 0.01, linesearch = LineSearches.BackTracking()), + callback = callback, maxiters = 20) ``` ## Result Analysis -Now, we'll plot the learned solutions of the neural ODE and compare them to our full physical model and the Newtonian model. +Now, we'll plot the learned solutions of the neural ODE and compare them to our full physical model and the Newtonian model. ```@example ude -reference_solution = solve(remake(prob, p = model_params, saveat = tsteps, tspan=tspan), - RK4(), dt = dt, adaptive=false) +reference_solution = solve(remake(prob, p = model_params, saveat = tsteps, tspan = tspan), + RK4(), dt = dt, adaptive = false) -optimized_solution = solve(remake(prob_nn, p = res2.minimizer, saveat = tsteps, tspan=tspan), - RK4(), dt = dt, adaptive=false) +optimized_solution = solve( + remake(prob_nn, p = res2.minimizer, saveat = tsteps, tspan = tspan), + RK4(), dt = dt, adaptive = false) Newtonian_prob = ODEProblem(NewtonianOrbitModel, u0, tspan, model_params) -Newtonian_solution = solve(remake(Newtonian_prob, p = model_params, saveat = tsteps, tspan=tspan), - RK4(), dt = dt, adaptive=false) +Newtonian_solution = solve( + remake(Newtonian_prob, p = model_params, saveat = tsteps, tspan = tspan), + RK4(), dt = dt, adaptive = false) true_orbit = soln2orbit(reference_solution, model_params) pred_orbit = soln2orbit(optimized_solution, model_params) @@ -575,44 +579,53 @@ Newt_waveform = compute_waveform(dt_data, Newtonian_solution, mass_ratio, model_ true_orbit = soln2orbit(reference_solution, model_params) pred_orbit = soln2orbit(optimized_solution, model_params) Newt_orbit = soln2orbit(Newtonian_solution, model_params) -plt = plot(true_orbit[1,:], true_orbit[2,:], linewidth = 2, label = "truth") -plot!(plt, pred_orbit[1,:], pred_orbit[2,:], linestyle = :dash, linewidth = 2, label = "prediction") -plot!(plt, Newt_orbit[1,:], Newt_orbit[2,:], linewidth = 2, label = "Newtonian") +plt = plot(true_orbit[1, :], true_orbit[2, :], linewidth = 2, label = "truth") +plot!(plt, pred_orbit[1, :], pred_orbit[2, :], + linestyle = :dash, linewidth = 2, label = "prediction") +plot!(plt, Newt_orbit[1, :], Newt_orbit[2, :], linewidth = 2, label = "Newtonian") ``` ```@example ude -plt = plot(tsteps,true_waveform, linewidth = 2, label = "truth", xlabel="Time", ylabel="Waveform") -plot!(plt,tsteps,pred_waveform, linestyle = :dash, linewidth = 2, label = "prediction") -plot!(plt,tsteps,Newt_waveform, linewidth = 2, label = "Newtonian") +plt = plot(tsteps, true_waveform, linewidth = 2, label = "truth", + xlabel = "Time", ylabel = "Waveform") +plot!(plt, tsteps, pred_waveform, linestyle = :dash, linewidth = 2, label = "prediction") +plot!(plt, tsteps, Newt_waveform, linewidth = 2, label = "Newtonian") ``` Now we'll do the same, but extrapolating the model out in time. ```@example ude -factor=5 - -extended_tspan = (tspan[1], factor*tspan[2]) -extended_tsteps = range(tspan[1], factor*tspan[2], length = factor*datasize) -reference_solution = solve(remake(prob, p = model_params, saveat = extended_tsteps, tspan=extended_tspan), - RK4(), dt = dt, adaptive=false) -optimized_solution = solve(remake(prob_nn, p = res2.minimizer, saveat = extended_tsteps, tspan=extended_tspan), - RK4(), dt = dt, adaptive=false) +factor = 5 + +extended_tspan = (tspan[1], factor * tspan[2]) +extended_tsteps = range(tspan[1], factor * tspan[2], length = factor * datasize) +reference_solution = solve( + remake(prob, p = model_params, saveat = extended_tsteps, tspan = extended_tspan), + RK4(), dt = dt, adaptive = false) +optimized_solution = solve( + remake(prob_nn, p = res2.minimizer, saveat = extended_tsteps, tspan = extended_tspan), + RK4(), dt = dt, adaptive = false) Newtonian_prob = ODEProblem(NewtonianOrbitModel, u0, tspan, model_params) -Newtonian_solution = solve(remake(Newtonian_prob, p = model_params, saveat = extended_tsteps, tspan=extended_tspan), - RK4(), dt = dt, adaptive=false) +Newtonian_solution = solve( + remake( + Newtonian_prob, p = model_params, saveat = extended_tsteps, tspan = extended_tspan), + RK4(), dt = dt, adaptive = false) true_orbit = soln2orbit(reference_solution, model_params) pred_orbit = soln2orbit(optimized_solution, model_params) Newt_orbit = soln2orbit(Newtonian_solution, model_params) -plt = plot(true_orbit[1,:], true_orbit[2,:], linewidth = 2, label = "truth") -plot!(plt, pred_orbit[1,:], pred_orbit[2,:], linestyle = :dash, linewidth = 2, label = "prediction") -plot!(plt, Newt_orbit[1,:], Newt_orbit[2,:], linewidth = 2, label = "Newtonian") +plt = plot(true_orbit[1, :], true_orbit[2, :], linewidth = 2, label = "truth") +plot!(plt, pred_orbit[1, :], pred_orbit[2, :], + linestyle = :dash, linewidth = 2, label = "prediction") +plot!(plt, Newt_orbit[1, :], Newt_orbit[2, :], linewidth = 2, label = "Newtonian") ``` ```@example ude true_waveform = compute_waveform(dt_data, reference_solution, mass_ratio, model_params)[1] pred_waveform = compute_waveform(dt_data, optimized_solution, mass_ratio, model_params)[1] Newt_waveform = compute_waveform(dt_data, Newtonian_solution, mass_ratio, model_params)[1] -plt = plot(extended_tsteps,true_waveform, linewidth = 2, label = "truth", xlabel="Time", ylabel="Waveform") -plot!(plt,extended_tsteps,pred_waveform, linestyle = :dash, linewidth = 2, label = "prediction") -plot!(plt,extended_tsteps,Newt_waveform, linewidth = 2, label = "Newtonian") +plt = plot(extended_tsteps, true_waveform, linewidth = 2, + label = "truth", xlabel = "Time", ylabel = "Waveform") +plot!(plt, extended_tsteps, pred_waveform, linestyle = :dash, + linewidth = 2, label = "prediction") +plot!(plt, extended_tsteps, Newt_waveform, linewidth = 2, label = "Newtonian") ``` diff --git a/docs/src/showcase/brusselator.md b/docs/src/showcase/brusselator.md index 53bcdd2068d..d86719a4496 100644 --- a/docs/src/showcase/brusselator.md +++ b/docs/src/showcase/brusselator.md @@ -123,7 +123,7 @@ partial differential equation transforms it into a new numerical problem. For ex | Finite Difference, Finite Volume, Finite Element, discretizing all variables | `NonlinearProblem` | | Finite Difference, Finite Volume, Finite Element, discretizing all variables except time | `ODEProblem`/`DAEProblem` | | Physics-Informed Neural Network | `OptimizationProblem` | -| Feynman-Kac Formula | `SDEProblem` | +| Feynman-Kac Formula | `SDEProblem` | | Universal Stochastic Differential Equation ([High dimensional PDEs](https://docs.sciml.ai/HighDimPDE/stable/)) | `OptimizationProblem` inverse problem over `SDEProblem` | Thus the process of solving a PDE is fundamentally about transforming its symbolic form @@ -147,7 +147,7 @@ dy = (y_max - y_min) / N order = 2 discretization = MOLFiniteDifference([x => dx, y => dy], t, approx_order = order, - grid_align = center_align) + grid_align = center_align) ``` Next, we `discretize` the system, converting the `PDESystem` in to an `ODEProblem`: @@ -278,8 +278,8 @@ function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata) end @btime solve(prob_sparse, - TRBDF2(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true), - save_everystep = false); + TRBDF2(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true), + save_everystep = false); ``` And now we're zooming! For more information on these performance improvements, check out diff --git a/docs/src/showcase/gpu_spde.md b/docs/src/showcase/gpu_spde.md index 580159cd605..f6d3db65d05 100644 --- a/docs/src/showcase/gpu_spde.md +++ b/docs/src/showcase/gpu_spde.md @@ -108,7 +108,7 @@ const Y = reshape([j for i in 1:100 for j in 1:100], N, N) const α₁ = 1.0 .* (X .>= 80) const Mx = Tridiagonal([1.0 for i in 1:(N - 1)], [-2.0 for i in 1:N], - [1.0 for i in 1:(N - 1)]) + [1.0 for i in 1:(N - 1)]) const My = copy(Mx) # Do the reflections, different for x and y operators Mx[2, 1] = 2.0 @@ -302,14 +302,14 @@ These last two ways enclose the pointer to our cache arrays locally but still pr function f(du,u,p,t) to the ODE solver. Now, since PDEs are large, many times we don't care about getting the whole timeseries. Using -the [output controls from DifferentialEquations.jl](http://diffeq.sciml.ai/latest/basics/common_solver_opts.html#Output-Control-1), we can make it only output the final timepoint. +the [output controls from DifferentialEquations.jl](https://diffeq.sciml.ai/latest/basics/common_solver_opts.html#Output-Control-1), we can make it only output the final timepoint. ```julia prob = ODEProblem(f, u0, (0.0, 100.0)) @time sol = solve(prob, ROCK2(), progress = true, save_everystep = false, - save_start = false); + save_start = false); @time sol = solve(prob, ROCK2(), progress = true, save_everystep = false, - save_start = false); + save_start = false); ``` Around 0.4 seconds. Much better. Also, if you're using VS Code, this'll give you a nice @@ -352,7 +352,7 @@ const Y = reshape([j for i in 1:100 for j in 1:100], N, N) const α₁ = 1.0 .* (X .>= 80) const Mx = Array(Tridiagonal([1.0 for i in 1:(N - 1)], [-2.0 for i in 1:N], - [1.0 for i in 1:(N - 1)])) + [1.0 for i in 1:(N - 1)])) const My = copy(Mx) Mx[2, 1] = 2.0 Mx[end - 1, end] = 2.0 @@ -440,12 +440,12 @@ end prob2 = ODEProblem(gf, gu0, (0.0, 100.0)) CUDA.allowscalar(false) # makes sure none of the slow fallbacks are used @time sol = solve(prob2, ROCK2(), progress = true, dt = 0.003, save_everystep = false, - save_start = false); + save_start = false); ``` ```@example spde @time sol = solve(prob2, ROCK2(), progress = true, dt = 0.003, save_everystep = false, - save_start = false); + save_start = false); ``` Go have fun. diff --git a/docs/src/showcase/massively_parallel_gpu.md b/docs/src/showcase/massively_parallel_gpu.md index 4d4ade85bb4..93c7ab6fb3d 100644 --- a/docs/src/showcase/massively_parallel_gpu.md +++ b/docs/src/showcase/massively_parallel_gpu.md @@ -90,14 +90,16 @@ sol = solve(monteprob, Tsit5(), EnsembleThreads(), trajectories = 10_000, saveat Now uhh, we just change `EnsembleThreads()` to `EnsembleGPUArray()` ```@example diffeqgpu -sol = solve(monteprob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()), trajectories = 10_000, saveat = 1.0f0) +sol = solve(monteprob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()), + trajectories = 10_000, saveat = 1.0f0) ``` Or for a more efficient version, `EnsembleGPUKernel()`. But that requires special solvers, so we also change to `GPUTsit5()`. ```@example diffeqgpu -sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()), trajectories = 10_000) +sol = solve( + monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()), trajectories = 10_000) ``` Okay, so that was anticlimactic, but that's the point: if it were harder than that, it diff --git a/docs/src/showcase/missing_physics.md b/docs/src/showcase/missing_physics.md index 57ac59ea8da..f10ff8b94cc 100644 --- a/docs/src/showcase/missing_physics.md +++ b/docs/src/showcase/missing_physics.md @@ -18,14 +18,14 @@ correct model and auto-complete it to find the missing physics. There are many packages which are used as part of this showcase. Let's detail what they are and how they are used. For the neural network training: -| Module | Description | -|:-------------------------------------------------------------------------------------------------------- |:----------------------------------------------------- | -| [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/) (DifferentialEquations.jl) | The numerical differential equation solvers | -| [SciMLSensitivity.jl](https://docs.sciml.ai/SciMLSensitivity/stable/) | The adjoint methods, defines gradients of ODE solvers | -| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The optimization library | -| [OptimizationOptimisers.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/optimisers/) | The optimization solver package with `Adam` | -| [OptimizationOptimJL.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/optim/) | The optimization solver package with `LBFGS` | -| [LineSearches.jl](https://julianlsolvers.github.io/LineSearches.jl/latest/index.html) | Line search algorithms package to be used with `LBFGS`| +| Module | Description | +|:-------------------------------------------------------------------------------------------------------- |:------------------------------------------------------ | +| [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/) (DifferentialEquations.jl) | The numerical differential equation solvers | +| [SciMLSensitivity.jl](https://docs.sciml.ai/SciMLSensitivity/stable/) | The adjoint methods, defines gradients of ODE solvers | +| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The optimization library | +| [OptimizationOptimisers.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/optimisers/) | The optimization solver package with `Adam` | +| [OptimizationOptimJL.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/optim/) | The optimization solver package with `LBFGS` | +| [LineSearches.jl](https://julianlsolvers.github.io/LineSearches.jl/latest/index.html) | Line search algorithms package to be used with `LBFGS` | For the symbolic model discovery: @@ -47,7 +47,7 @@ And external libraries: | Module | Description | |:---------------------------------------------------------------------------- |:--------------------------------------------------- | -| [Lux.jl](http://lux.csail.mit.edu/stable/) | The deep learning (neural network) framework | +| [Lux.jl](https://lux.csail.mit.edu/stable/) | The deep learning (neural network) framework | | [ComponentArrays.jl](https://jonniedie.github.io/ComponentArrays.jl/stable/) | For the `ComponentArray` type to match Lux to SciML | | [Plots.jl](https://docs.juliaplots.org/stable/) | The plotting and visualization library | | [StableRNGs.jl](https://docs.juliaplots.org/stable/) | Stable random seeding | @@ -138,7 +138,7 @@ rbf(x) = exp.(-(x .^ 2)) # Multilayer FeedForward const U = Lux.Chain(Lux.Dense(2, 5, rbf), Lux.Dense(5, 5, rbf), Lux.Dense(5, 5, rbf), - Lux.Dense(5, 2)) + Lux.Dense(5, 2)) # Get the initial parameters and state variables of the model p, st = Lux.setup(rng, U) const _st = st @@ -196,12 +196,12 @@ Knowing this, our `predict` function looks like: function predict(θ, X = Xₙ[:, 1], T = t) _prob = remake(prob_nn, u0 = X, tspan = (T[1], T[end]), p = θ) Array(solve(_prob, Vern7(), saveat = T, - abstol = 1e-6, reltol = 1e-6, - sensealg=QuadratureAdjoint(autojacvec=ReverseDiffVJP(true)))) + abstol = 1e-6, reltol = 1e-6, + sensealg = QuadratureAdjoint(autojacvec = ReverseDiffVJP(true)))) end ``` -There are many choices for the combination of sensitivity algorithm and automatic differentiation library (see [Choosing a Sensitivity Algorithm](https://docs.sciml.ai/SciMLSensitivity/stable/manual/differential_equation_sensitivities/#Choosing-a-Sensitivity-Algorithm). For example, you could have used `sensealg=ForwardDiffSensitivity()`. +There are many choices for the combination of sensitivity algorithm and automatic differentiation library (see [Choosing a Sensitivity Algorithm](https://docs.sciml.ai/SciMLSensitivity/stable/manual/differential_equation_sensitivities/#Choosing-a-Sensitivity-Algorithm). For example, you could have used `sensealg=ForwardDiffSensitivity()`. Now, for our loss function, we solve the ODE at our new parameters and check its L2 loss against the dataset. Using our `predict` function, this looks like: @@ -236,7 +236,7 @@ end Now we're ready to train! To run the training process, we will need to build an [`OptimizationProblem`](https://docs.sciml.ai/Optimization/stable/API/optimization_problem/). Because we have a lot of parameters, we will use -[Zygote.jl](https://docs.sciml.ai/Zygote/stable/). Optimization.jl makes the choice of +[Zygote.jl](https://fluxml.ai/Zygote.jl/dev/). Optimization.jl makes the choice of automatic differentiation easy, just by specifying an `adtype` in the [`OptimizationFunction` construction](https://docs.sciml.ai/Optimization/stable/API/optimization_function/) @@ -258,7 +258,8 @@ Thus we first solve the optimization problem with ADAM. Choosing a learning rate (tuned to be as high as possible that doesn't tend to make the loss shoot up), we see: ```@example ude -res1 = Optimization.solve(optprob, OptimizationOptimisers.Adam(), callback = callback, maxiters = 5000) +res1 = Optimization.solve( + optprob, OptimizationOptimisers.Adam(), callback = callback, maxiters = 5000) println("Training loss after $(length(losses)) iterations: $(losses[end])") ``` @@ -267,7 +268,8 @@ second optimization, and run it with BFGS. This looks like: ```@example ude optprob2 = Optimization.OptimizationProblem(optf, res1.u) -res2 = Optimization.solve(optprob2, LBFGS(linesearch = BackTracking()), callback = callback, maxiters = 1000) +res2 = Optimization.solve( + optprob2, LBFGS(linesearch = BackTracking()), callback = callback, maxiters = 1000) println("Final training loss after $(length(losses)) iterations: $(losses[end])") # Rename the best candidate @@ -283,9 +285,9 @@ How well did our neural network do? Let's take a look: ```@example ude # Plot the losses pl_losses = plot(1:5000, losses[1:5000], yaxis = :log10, xaxis = :log10, - xlabel = "Iterations", ylabel = "Loss", label = "ADAM", color = :blue) + xlabel = "Iterations", ylabel = "Loss", label = "ADAM", color = :blue) plot!(5001:length(losses), losses[5001:end], yaxis = :log10, xaxis = :log10, - xlabel = "Iterations", ylabel = "Loss", label = "LBFGS", color = :red) + xlabel = "Iterations", ylabel = "Loss", label = "LBFGS", color = :red) ``` Next, we compare the original data to the output of the UDE predictor. Note that we can even create more samples from the underlying model by simply adjusting the time steps! @@ -297,7 +299,7 @@ ts = first(solution.t):(mean(diff(solution.t)) / 2):last(solution.t) X̂ = predict(p_trained, Xₙ[:, 1], ts) # Trained on noisy data vs real solution pl_trajectory = plot(ts, transpose(X̂), xlabel = "t", ylabel = "x(t), y(t)", color = :red, - label = ["UDE Approximation" nothing]) + label = ["UDE Approximation" nothing]) scatter!(solution.t, transpose(Xₙ), color = :black, label = ["Measurements" nothing]) ``` @@ -310,7 +312,7 @@ Ȳ = [-p_[2] * (X̂[1, :] .* X̂[2, :])'; p_[3] * (X̂[1, :] .* X̂[2, :])'] Ŷ = U(X̂, p_trained, st)[1] pl_reconstruction = plot(ts, transpose(Ŷ), xlabel = "t", ylabel = "U(x,y)", color = :red, - label = ["UDE Approximation" nothing]) + label = ["UDE Approximation" nothing]) plot!(ts, transpose(Ȳ), color = :black, label = ["True Interaction" nothing]) ``` @@ -319,7 +321,7 @@ And have a nice look at all the information: ```@example ude # Plot the error pl_reconstruction_error = plot(ts, norm.(eachcol(Ȳ - Ŷ)), yaxis = :log, xlabel = "t", - ylabel = "L2-Error", label = nothing, color = :red) + ylabel = "L2-Error", label = nothing, color = :red) pl_missing = plot(pl_reconstruction, pl_reconstruction_error, layout = (2, 1)) pl_overall = plot(pl_trajectory, pl_missing) @@ -395,12 +397,12 @@ regressions: ```@example ude options = DataDrivenCommonOptions(maxiters = 10_000, - normalize = DataNormalization(ZScoreTransform), - selector = bic, digits = 1, - data_processing = DataProcessing(split = 0.9, - batchsize = 30, - shuffle = true, - rng = StableRNG(1111))) + normalize = DataNormalization(ZScoreTransform), + selector = bic, digits = 1, + data_processing = DataProcessing(split = 0.9, + batchsize = 30, + shuffle = true, + rng = StableRNG(1111))) full_res = solve(full_problem, basis, opt, options = options) full_eqs = get_basis(full_res) @@ -409,12 +411,12 @@ println(full_res) ```@example ude options = DataDrivenCommonOptions(maxiters = 10_000, - normalize = DataNormalization(ZScoreTransform), - selector = bic, digits = 1, - data_processing = DataProcessing(split = 0.9, - batchsize = 30, - shuffle = true, - rng = StableRNG(1111))) + normalize = DataNormalization(ZScoreTransform), + selector = bic, digits = 1, + data_processing = DataProcessing(split = 0.9, + batchsize = 30, + shuffle = true, + rng = StableRNG(1111))) ideal_res = solve(ideal_problem, basis, opt, options = options) ideal_eqs = get_basis(ideal_res) @@ -423,12 +425,12 @@ println(ideal_res) ```@example ude options = DataDrivenCommonOptions(maxiters = 10_000, - normalize = DataNormalization(ZScoreTransform), - selector = bic, digits = 1, - data_processing = DataProcessing(split = 0.9, - batchsize = 30, - shuffle = true, - rng = StableRNG(1111))) + normalize = DataNormalization(ZScoreTransform), + selector = bic, digits = 1, + data_processing = DataProcessing(split = 0.9, + batchsize = 30, + shuffle = true, + rng = StableRNG(1111))) nn_res = solve(nn_problem, basis, opt, options = options) nn_eqs = get_basis(nn_res) @@ -504,29 +506,29 @@ c3 = :blue # RGBA(255/255,90/255,0,1) # Orange c4 = :purple # RGBA(153/255,50/255,204/255,1) # Purple p1 = plot(t, abs.(Array(solution) .- estimate)' .+ eps(Float32), - lw = 3, yaxis = :log, title = "Timeseries of UODE Error", - color = [3 :orange], xlabel = "t", - label = ["x(t)" "y(t)"], - titlefont = "Helvetica", legendfont = "Helvetica", - legend = :topright) + lw = 3, yaxis = :log, title = "Timeseries of UODE Error", + color = [3 :orange], xlabel = "t", + label = ["x(t)" "y(t)"], + titlefont = "Helvetica", legendfont = "Helvetica", + legend = :topright) # Plot L₂ p2 = 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) + title = "Neural Network Fit of U2(t)", color = c1, + label = "Neural Network", xaxis = "x", yaxis = "y", + titlefont = "Helvetica", legendfont = "Helvetica", + legend = :bottomright) 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"], - title = "Extrapolated Fit From Short Training Data", - titlefont = "Helvetica", legendfont = "Helvetica", - markersize = 5) + title = "Extrapolated Fit From Short Training Data", + titlefont = "Helvetica", legendfont = "Helvetica", + markersize = 5) plot!(p3, true_solution_long, color = [c1 c2], linestyle = :dot, lw = 5, - label = ["True x(t)" "True y(t)"]) + label = ["True x(t)" "True y(t)"]) plot!(p3, estimate_long, color = [c3 c4], lw = 1, - label = ["Estimated x(t)" "Estimated y(t)"]) + label = ["Estimated x(t)" "Estimated y(t)"]) 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) diff --git a/docs/src/showcase/ode_types.md b/docs/src/showcase/ode_types.md index c586b6411d2..396db9764d0 100644 --- a/docs/src/showcase/ode_types.md +++ b/docs/src/showcase/ode_types.md @@ -10,7 +10,7 @@ few useful/interesting types that can be used: | BigFloat | Base Julia | Higher precision solutions | | ArbFloat | [ArbFloats.jl](https://github.com/JuliaArbTypes/ArbFloats.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 | +| Particles | [MonteCarloMeasurements.jl](https://github.com/baggepinnen/MonteCarloMeasurements.jl) | Uncertainty propagation | | Unitful | [Unitful.jl](https://painterqubits.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 | @@ -18,7 +18,7 @@ few useful/interesting types that can be used: | Num | [Symbolics.jl](https://symbolics.juliasymbolics.org/stable/) | Build symbolic expressions of ODE solution approximations | | Taylor | [TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl) | Build a Taylor series around a solution point | | Dual | [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/) | Perform forward-mode automatic differentiation | -| TrackedArray\TrackedReal | [ReverseDiff.jl](https://juliadiff.org/ReverseDiff.jl/) | Perform reverse-mode automatic differentiation | +| TrackedArray\TrackedReal | [ReverseDiff.jl](https://juliadiff.org/ReverseDiff.jl/) | Perform reverse-mode automatic differentiation | and on and on. That's only a subset of types people have effectively used on the SciML tools. @@ -33,13 +33,13 @@ time values in the same type as tspan, and the solution values in the same type initial condition. !!! note - + Support for this feature is restricted to the native algorithms of OrdinaryDiffEq.jl. The other solvers such as Sundials.jl, and ODEInterface.jl are incompatible with some number systems. !!! warn - + Adaptive timestepping requires that the time type is compatible with `sqrt` and `^` functions. Thus for example, `tspan` cannot be `Int` if adaptive timestepping is chosen. diff --git a/docs/src/showcase/pinngpu.md b/docs/src/showcase/pinngpu.md index 563149467cc..0875715cbcb 100644 --- a/docs/src/showcase/pinngpu.md +++ b/docs/src/showcase/pinngpu.md @@ -103,7 +103,7 @@ domains = [t ∈ Interval(t_min, t_max), ``` !!! note - + We used the wildcard form of the variable definition `@variables u(..)` which then requires that we always specify what the dependent variables of `u` are. This is because in the boundary conditions we change from using `u(t,x,y)` to more specific points and lines, like `u(t,x_max,y)`. @@ -117,10 +117,10 @@ We will use a simple multi-layer perceptron, like: using Lux inner = 25 chain = Chain(Dense(3, inner, Lux.σ), - Dense(inner, inner, Lux.σ), - Dense(inner, inner, Lux.σ), - Dense(inner, inner, Lux.σ), - Dense(inner, 1)) + Dense(inner, inner, Lux.σ), + Dense(inner, inner, Lux.σ), + Dense(inner, inner, Lux.σ), + Dense(inner, 1)) ps = Lux.setup(Random.default_rng(), chain)[1] ``` @@ -140,8 +140,8 @@ ps = ps |> ComponentArray |> gpud .|> Float64 ```@example pinn strategy = GridTraining(0.05) discretization = PhysicsInformedNN(chain, - strategy, - init_params = ps) + strategy, + init_params = ps) prob = discretize(pde_system, discretization) ``` @@ -179,9 +179,9 @@ function plot_(res) anim = @animate for (i, t) in enumerate(0:0.05:t_max) @info "Animating frame $i..." u_real = reshape([analytic_sol_func(t, x, y) for x in xs for y in ys], - (length(xs), length(ys))) + (length(xs), length(ys))) u_predict = reshape([Array(phi([t, x, y], res.u))[1] for x in xs for y in ys], - length(xs), length(ys)) + length(xs), length(ys)) u_error = abs.(u_predict .- u_real) title = @sprintf("predict, t = %.3f", t) p1 = plot(xs, ys, u_predict, st = :surface, label = "", title = title) diff --git a/docs/src/showcase/symbolic_analysis.md b/docs/src/showcase/symbolic_analysis.md index f21e2cdcbd4..b95fff2599e 100644 --- a/docs/src/showcase/symbolic_analysis.md +++ b/docs/src/showcase/symbolic_analysis.md @@ -118,7 +118,7 @@ Did you implement the DAE incorrectly? No. Is the solver broken? No. It turns out that this is a property of the DAE that we are attempting to solve. This kind of DAE is known as an index-3 DAE. For a complete discussion of DAE -index, see [this article](http://www.scholarpedia.org/article/Differential-algebraic_equations). +index, see [this article](https://www.scholarpedia.org/article/Differential-algebraic_equations). Essentially, the issue here is that we have 4 differential variables (``x``, ``v_x``, ``y``, ``v_y``) and one algebraic variable ``T`` (which we can know because there is no `D(T)` term in the equations). An index-1 DAE always satisfies that the Jacobian of @@ -171,7 +171,7 @@ we can avoid this by using methods like for automatically performing this reduction to index 1. While this requires the ModelingToolkit symbolic form, we use `modelingtoolkitize` to transform the numerical code into symbolic code, run `structural_simplify` to -simplify the system and lower the index, +simplify the system and lower the index, then transform back to numerical code with `ODEProblem`, and solve with a numerical solver. Let's try that out: @@ -192,7 +192,7 @@ plot and the plot is shown in the same order as the original numerical code. Note that we can even go a bit further. If we use the `ODEProblem` -constructor, we represent the mass matrix DAE of the index-reduced system, +constructor, we represent the mass matrix DAE of the index-reduced system, which can be solved via: ```@example indexred @@ -257,7 +257,7 @@ eqs = [ D(x4) ~ -k5 * x4 / (k6 + x4), D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6), D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10, - D(x7) ~ k9 * x6 * (k10 - x6) / k10, + D(x7) ~ k9 * x6 * (k10 - x6) / k10 ] # define the output functions (quantities that can be measured) @@ -273,7 +273,7 @@ After that, we are ready to check the system for local identifiability: # query local identifiability # we pass the ode-system local_id_all = assess_local_identifiability(de, measured_quantities = measured_quantities, - p = 0.99) + p = 0.99) # [ Info: Preproccessing `ModelingToolkit.ODESystem` object # 6-element Vector{Bool}: # 1 @@ -291,7 +291,7 @@ Let's try to check specific parameters and their combinations ```julia to_check = [k5, k7, k10 / k9, k5 + k6] local_id_some = assess_local_identifiability(de, measured_quantities = measured_quantities, - funcs_to_check = to_check, p = 0.99) + funcs_to_check = to_check, p = 0.99) # 4-element Vector{Bool}: # 1 # 1 @@ -333,7 +333,7 @@ eqs = [ D(x1) ~ -b * x1 + 1 / (c + x4), D(x2) ~ a * x1 - beta * x2, D(x3) ~ g * x2 - delta * x3, - D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3, + D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3 ] measured_quantities = [y ~ x1 + x2, y2 ~ x2] @@ -367,7 +367,7 @@ eqs = [ D(x1) ~ -b * x1 + 1 / (c + x4), D(x2) ~ a * x1 - beta * x2 - u1, D(x3) ~ g * x2 - delta * x3 + u2, - D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3, + D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3 ] measured_quantities = [y ~ x1 + x2, y2 ~ x2] @@ -377,7 +377,7 @@ to_check = [b, c] ode = ODESystem(eqs, t, name = :GoodwinOsc) global_id = assess_identifiability(ode, measured_quantities = measured_quantities, - funcs_to_check = to_check, p = 0.9) + funcs_to_check = to_check, p = 0.9) # Dict{Num, Symbol} with 2 entries: # b => :globally # c => :globally