Merge pull request #740 from SciML/new_SciMLOptimization_paramfitting_tutorial

New SciML/Optimization for fitting parameters of ODEs tutorial tutorial

Author: TorkelE

docs/Project.toml

@@ -2,6 +2,7 @@
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DiffEqParamEstim = "1130ab10-4a5a-5621-a13d-e4788d82bd4c"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -12,6 +13,8 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
@@ -30,6 +33,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 BifurcationKit = "0.3"
 Catalyst = "13"
 DataFrames = "1"
+DiffEqParamEstim = "2.1"
 DifferentialEquations = "7.7"
 Distributions = "0.25"
 Documenter = "0.27"
@@ -39,6 +43,8 @@ ModelingToolkit = "8.47"
 NonlinearSolve = "3.4.0"
 Optim = "1"
 Optimization = "3.19"
+OptimizationNLopt = "0.1.8"
+OptimizationOptimJL = "0.1.14"
 OptimizationOptimisers = "0.1.1"
 OrdinaryDiffEq = "6"
 PEtab = "2"

docs/pages.jl

@@ -16,8 +16,9 @@ pages = Any["Home" => "index.md",
                                            "catalyst_applications/homotopy_continuation.md",
                                            "catalyst_applications/nonlinear_solve.md",
                                            "catalyst_applications/bifurcation_diagrams.md"],
-            "Inverse Problems" => Any["inverse_problems/parameter_estimation.md",
-                                      "inverse_problems/petab_ode_param_fitting.md",
-                                      "inverse_problems/structural_identifiability.md"],
+            "Inverse Problems" => Any["inverse_problems/optimization_ode_param_fitting.md",
+                                        "inverse_problems/petab_ode_param_fitting.md",
+                                        "inverse_problems/structural_identifiability.md",
+                                        "Inverse problem examples" => Any["inverse_problems/examples/ode_fitting_oscillation.md"]],
             "FAQs" => "faqs.md",
             "API" => "api/catalyst_api.md"]

docs/src/catalyst_applications/advanced_simulations.md

@@ -9,7 +9,7 @@ DifferentialEquations package](https://docs.sciml.ai/DiffEqDocs/stable/), which
 Catalyst uses for all simulations.
 
 
-## Monte Carlo simulations using `EnsembleProblem`s
+## [Monte Carlo simulations using `EnsembleProblem`s](@id advanced_simulations_ensemble_problems)
 In many contexts one needs to run multiple simulations of a model, for example
 to collect statistics of SDE or jump process solutions, or to systematically
 vary parameter values within a model. While it is always possible to manually

docs/src/catalyst_applications/simulation_structure_interfacing.md

@@ -33,7 +33,7 @@ oprob[:X1] = 10.0
 ```
 with parameters using the same notation.
 
-#### Remaking problems using the `remake` function
+#### [Remaking problems using the `remake` function](@id simulation_structure_interfacing_remake)
 Typically, when modifying problems, it is recommended to use the `remake` function. Unlike when we do `oprob[:X1] = 10.0` (which modifies the problem in question), `remake` creates a new problem object. The `remake` function takes a problem as input, and any fields you wish to modify (and their new values) as optional inputs. Thus, we can do:
 ```@example ex1
 using DifferentialEquations

docs/src/inverse_problems/parameter_estimation.md renamed to docs/src/inverse_problems/examples/ode_fitting_oscillation.md

@@ -1,24 +1,16 @@
-# [Parameter Estimation](@id parameter_estimation)
-The parameters of a model, generated by Catalyst, can be estimated using various
-packages available in the Julia ecosystem. Refer
-[here](https://docs.sciml.ai/Overview/stable/highlevels/inverse_problems/) for
-more extensive information. Below follows a quick tutorial of how
-[Optimization.jl](https://docs.sciml.ai/Optimization/stable/) can be used to fit
-a parameter set to data.
+# [Fitting Parameters for an Oscillatory System](@id parameter_estimation)
+In this example we will use [Optimization.jl](https://github.com/SciML/Optimization.jl) to fit the parameters of an oscillatory system (the Brusselator) to data. Here, special consideration is taken to avoid reaching a local minimum. Instead of fitting the entire time series directly, we will start with fitting parameter values for the first period, and then use those as an initial guess for fitting the next (and then these to find the next one, and so on). Using this procedure is advantageous for oscillatory systems, and enables us to reach the global optimum.
 
 First, we fetch the required packages.
 ```@example pe1
 using Catalyst
-using DifferentialEquations
-using SciMLSensitivity
+using OrdinaryDiffEq
 using Optimization
-
-# for the ADAM optimizer
-using OptimizationOptimisers
+using OptimizationOptimisers # Required for the ADAM optimizer.
+using SciMLSensitivity # Required for `Optimization.AutoZygote()` automatic differentiation option.
 ```
 
-Next, we declare our model. For our example, we will use the Brusselator, a
-simple oscillator.
+Next, we declare our model, the Brusselator oscillator.
 ```@example pe1
 brusselator = @reaction_network begin
     A, ∅ --> X
@@ -27,10 +19,11 @@ brusselator = @reaction_network begin
     1, X --> ∅
 end
 p_real = [:A => 1., :B => 2.]
+nothing # hide
 ```
 
 We simulate our model, and from the simulation generate sampled data points
-(with added noise), to which we will attempt to fit a parameter et.
+(to which we add noise). We will use this data to fit the parameters of our model.
 ```@example pe1
 u0 = [:X => 1.0, :Y => 1.0]
 tspan = (0.0, 30.0)
@@ -54,8 +47,8 @@ scatter!(sample_times, sample_vals'; color = [:blue :red], legend = nothing)
 
 Next, we create a function to fit the parameters using the `ADAM` optimizer. For
 a given initial estimate of the parameter values, `pinit`, this function will
-fit parameter values, `p`, to our data samples. `tend` is used to indicate the
-time interval over which to fit to the ODE solution.
+fit parameter values, `p`, to our data samples. We use `tend` to indicate the
+time interval over which we fit the model.
 ```@example pe1
 function optimise_p(pinit, tend)
     function loss(p, _)
@@ -124,14 +117,14 @@ The final parameter estimate is then
 ```@example pe1
 p_estimate
 ```
-which is close to the actual parameters of `[1.0, 2.0]`.
+which is close to the actual parameter set of `[1.0, 2.0]`.
 
-## Why we fit the parameters in iterations.
-The reason we chose to fit the model on a smaller interval to begin with, and
+## Why we fit the parameters in iterations
+As previously mentioned, the reason we chose to fit the model on a smaller interval to begin with, and
 then extend the interval, is to avoid getting stuck in a local minimum. Here
 specifically, we chose our initial interval to be smaller than a full cycle of
 the oscillation. If we had chosen to fit a parameter set on the full interval
-immediately we would have received an inferior solution.
+immediately we would have obtained poor fit and an inaccurate estimate for the parameters.
 ```@example pe1
 p_estimate = optimise_p([5.0,5.0], 30.0)