Skip to content

Commit d817af7

Browse files
TorkelETorkelE
committed
init
1 parent 91bf817 commit d817af7

File tree

2 files changed

+186
-20
lines changed

2 files changed

+186
-20
lines changed

‎docs/src/inverse_problems/optimization_ode_param_fitting.md

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
# [Parameter Fitting for ODEs using SciML/Optimization.jl and DiffEqParamEstim.jl](@id optimization_parameter_fitting)
2+
Fitting parameters to data involves solving an optimisation problem (that is, finding the parameter set that optimally fits you model to your data, typically by minimising a cost function). The SciML ecosystem's primary package for solving optimisation problems is [Optimization.jl](https://github.com/SciML/Optimization.jl). It provides access to a variety of solvers from a single common interface, wrapping a large number of optimisation methods that have been implemented in Julia into a single common interface.
3+
4+
This tutorial both demonstrate how to create parameter fitting cost functions using the [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl) package, and how to use Optimization.jl to minimise these. Optimization.jl can also be used in other contexts, such as finding parameter sets that maximises the magnitude of some system behaviour. More details on how to use these packages can be found in their [respective](https://docs.sciml.ai/Optimization/stable/) [documentations](https://docs.sciml.ai/DiffEqParamEstim/stable/).
5+
6+
## Basic example
7+
8+
Let us consider a simple catalysis network, where an enzyme (*E*) turns a substrate (*S*) into a product (*P*):
9+
```@example diffeq_param_estim_1
10+
using Catalyst
11+
rn = @reaction_network begin
12+
kB, S + E --> SE
13+
kD, SE --> S + E
14+
kP, SE --> P + E
15+
end
16+
```
17+
From some known initial condition, and a true parameter set (which we later want to recover from the data) we generate synthetic data (on which we will demonstrate the fitting process).
18+
```@example diffeq_param_estim_1
19+
# Define initial conditions and parameters.
20+
u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
21+
p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
22+
23+
# Generate synthetic data.
24+
using OrdinaryDiffEq
25+
oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true)
26+
true_sol = solve(oprob_true, Tsit5())
27+
data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
28+
data_ts = data_sol.t[2:end]
29+
data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
30+
31+
# Plots the true solutions and the (synthetic data) measurements.
32+
using Plots
33+
plot(true_sol; idxs=:P, label="True solution", lw=8)
34+
plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue)
35+
```
36+
37+
Next, we will use DiffEqParamEstim.jl to build a loss function that describing how well our model's simulations fit the data.
38+
```@example diffeq_param_estim_1
39+
using DiffEqParamEstim, Optimization
40+
p_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0]
41+
oprob = ODEProblem(rn, u0, (0.0, 10.0), p_dummy)
42+
loss_function = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); maxiters=10000, verbose=false, save_idxs=4)
43+
nothing # hide
44+
```
45+
To `build_loss_objective` we provide the following arguments:
46+
- `oprob`: The `ODEProblem` with which we simulate our model (using some dummy parameter values, since we do not know these).
47+
- `Tsit5()`: The numeric integrator we wish to simulate our model with.
48+
- `L2Loss(data_ts, data_vals)`: Defines the loss function. While [other alternatives](https://docs.sciml.ai/DiffEqParamEstim/stable/getting_started/#Alternative-Cost-Functions-for-Increased-Robustness) are available, `L2Loss` is the simplest one (measuring the sum of squared distances between model simulations and data measurements). Its first argument is the time points at which the data is collected, and the second is the data's values.
49+
- `Optimization.AutoForwardDiff()`: Our choice of [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) framework.
50+
51+
Furthermore, we can pass any number of additional optional arguments, these are then passed to the internal `solve()` function (that are used to simulate our ODE). Here we provide the following additional arguments:
52+
- `maxiters=10000`: If the ODE integrator takes a very large number of steps, that is typically the sign of a very poor fit. Reducing the `maxiters` threshold reduces the time we waste on evaluating such models.
53+
- `verbose=false`: The simulation of models with highly unsuitable parameter sets typically generate various warnings (such as simulation terminations due to reaching the `maxiters` value). To avoid an overflow of such (here unnecessary) warnings as we evaluate a large number of parameter sets, we turn warnings off.
54+
- `save_idxs=4`: The measured species (*P*) is the 4th species in our species vector (`species(rn)`). To ensure that the concentration of the right species is evaluated against the data, we set the numeric integrator to only save the value of this species.
55+
56+
Now we can create an `OptimizationProblem` using our `loss_function` and some initial guess of parameter values from which the optimiser will start:
57+
```@example diffeq_param_estim_1
58+
optprob = OptimizationProblem(loss_function, [1.0, 1.0, 1.0])
59+
nothing # hide
60+
```
61+
62+
!!! note
63+
`OptimizationProblem` cannot currently accept parameter values in the form of a map (e.g. `[:kB => 1.0, :kD => 1.0, :kP => 1.0]`). These must be provided as individual values (using the same order as the parameters occur in in the `parameters(rs)` vector).
64+
65+
Finally, we can optimise `optprob` to find the parameter set that best fits our data. Optimization.jl does not provide any optimisation methods by default. Instead, for each supported optimisation package, it provides a corresponding wrapper-package to import that optimisation package for using with Optimization. E.g., if we wish to [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) method, we must install and import the OptimizationOptimJL package. A summary of all, by Optimization.jl, supported optimisation packages can be found [here](https://docs.sciml.ai/Optimization/stable/#Overview-of-the-Optimizers). Here, we import the Optim.jl package and uses it to minimise our cost function (thus finding a parameter set that fits the data):
66+
```@example diffeq_param_estim_1
67+
using OptimizationOptimJL
68+
optsol = solve(optprob, Optim.NelderMead())
69+
nothing # hide
70+
```
71+
72+
We can now simulate our model for the corresponding parameter set, checking that it fits our data.
73+
```@example diffeq_param_estim_1
74+
oprob_fitted = remake(oprob; p=optsol.u)
75+
fitted_sol = solve(oprob_fitted, Tsit5())
76+
plot!(fitted_sol; idxs=:P, label="Fitted solution", linestyle=:dash, lw=6, color=:lightblue)
77+
```
78+
79+
!!! note
80+
Here, a good exercise is to check the resulting parameter set and note that, while it creates a good fit to the data, it does not actually correspond to the original parameter set. [Identifiability](https://www.sciencedirect.com/science/article/pii/S1364815218307278) is a concept that studies how to deal with this problem.
81+
82+
Say that we instead would like to use the [Broyden–Fletcher–Goldfarb–Shannon](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) algorithm, as implemented by the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl) package. In this case we would run:
83+
```@example diffeq_param_estim_1
84+
using OptimizationNLopt
85+
sol = solve(optprob, NLopt.LD_LBFGS())
86+
nothing # hide
87+
```
88+
89+
## Optimisation problems with data for multiple species
90+
Imagine that, in our previous example, we had measurements the concentration of both *S* and *P*:
91+
```@example diffeq_param_estim_1
92+
data_vals_S = (0.8 .+ 0.4*rand(10)) .* data_sol[:S][2:end]
93+
data_vals_P = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
94+
95+
plot(true_sol; idxs=[:S, :P], label=["True S" "True P"], lw=8)
96+
plot!(data_ts, data_vals_S; label="Measured S", seriestype=:scatter, ms=6, color=:blue)
97+
plot!(data_ts, data_vals_P; label="Measured P", seriestype=:scatter, ms=6, color=:red)
98+
```
99+
100+
In this case we would have to use the `L2Loss(data_ts, hcat(data_vals_S, data_vals_P))` and `save_idxs=[1,4]` arguments in `loss_function`:
101+
```@example diffeq_param_estim_1
102+
loss_function_S_P = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, Array(hcat(data_vals_S, data_vals_P)')), Optimization.AutoForwardDiff(); maxiters=10000, verbose=false, save_idxs=[1,4])
103+
nothing # hide
104+
```
105+
Here, `Array(hcat(data_vals_S, data_vals_P)')` is required to pu the data in the right form (in this case, a 2x10 matrix).
106+
107+
We can now fit our model to data and plot the results:
108+
```@example diffeq_param_estim_1
109+
optprob_S_P = OptimizationProblem(loss_function_S_P, [1.0,1.0, 1.0])
110+
optsol_S_P = solve(optprob_S_P, Optim.NelderMead())
111+
oprob_fitted_S_P = remake(oprob; p=optsol_S_P.u)
112+
fitted_sol_S_P = solve(oprob_fitted_S_P, Tsit5())
113+
plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle=:dash, lw=6, color=[:lightblue :pink])
114+
```
115+
116+
## Setting parameter constraints and boundaries
117+
Sometimes, it is desired to set boundaries on parameter values. Indeed, this can speed up the optimisation process (by preventing searching through unfeasible parts of parameter space) or be a prerequisite of some optimisation methods. This can be done by passing the `lb` (lower bounds) and `up` (upper bounds) arguments to `OptimizationProblem`. These are vectors (of the same length as the number of parameters), with each argument corresponding to the boundary value of the parameter with the same index (as used in the `parameters(rn)` vector). If we wish to constrain each parameter to the interval *(0.1, 10.0)* this can be done through:
118+
```@example diffeq_param_estim_1
119+
optprob = OptimizationProblem(loss_function, [1.0, 1.0, 1.0]; lb = [0.1, 0.1, 0.1], ub = [10.0, 10.0, 10.0])
120+
nothing # hide
121+
```
122+
123+
In addition to boundaries, Optimization.jl also supports setting [linear and non-linear constraints](https://docs.sciml.ai/Optimization/stable/tutorials/constraints/#constraints) on its output solution.
124+
125+
## Parameter fitting with known parameters
126+
If we from previous knowledge know that *kD = 0.1*, and only would like to fit the values of *kD* and *kP*, this can be achieved through `build_loss_objective`'s `prob_generator` argument. First, we create a function (`fixed_p_prob_generator`) that modifies our `ODEProblem` to incorporate this knowledge:
127+
```@example diffeq_param_estim_1
128+
fixed_p_prob_generator(prob, p) = remake(prob; p = vcat(p[1], 0.1, p[2]))
129+
nothing # hide
130+
```
131+
Here, it takes the `ODEProblem` (`prob`) we simulates, and the parameter set used (`p`), during the optimisation process, and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_remake)). Now we create our modified loss function:
132+
```@example diffeq_param_estim_1
133+
loss_function_fixed_kD = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); prob_generator = fixed_p_prob_generator, maxiters=10000, verbose=false, save_idxs=4)
134+
nothing # hide
135+
```
136+
137+
We can create an optimisation problem from this one like previously, but keeping in mind that it (and its output results) only contains two parameter values (*kB* and *kP):
138+
```@example diffeq_param_estim_1
139+
optprob_fixed_kD = OptimizationProblem(loss_function_fixed_kD, [1.0, 1.0])
140+
optsol_fixed_kD = solve(optprob_fixed_kD, Optim.NelderMead())
141+
nothing # hide
142+
```
143+
144+
## Fitting parameters on the logarithmic scale
145+
Often it can be advantageous to fit parameters on a [logarithmic, rather than linear, scale](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646). The best way to proceed is to simply replace each parameters in the model definition by its logarithmic version:
146+
```@example diffeq_param_estim_2
147+
using Catalyst
148+
rn = @reaction_network begin
149+
10^kB, S + E --> SE
150+
10^kD, SE --> S + E
151+
10^kP, SE --> P + E
152+
end
153+
```
154+
And then proceeding, by keeping in mind that parameter values are logarithmic. Here, setting
155+
```@example diffeq_param_estim_2
156+
p_true = [:kB => 0.0, :kD => -1.0, :kP => 10^(0.5)]
157+
nothing # hide
158+
```
159+
corresponds to the same true parameter values as used previously (`[:kB => 1.0, :kD => 0.1, :kP => 0.5]`).
160+
161+
## Parameter fitting to multiple experiments
162+
Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref advanced_simulations_ensemble_problems). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
163+
164+
## Optimisation solver options
165+
Optimization.jl supports various [optimisation solver options](https://docs.sciml.ai/Optimization/stable/API/solve/) that can be supplied to the `solve` command. For example, to set a maximum number of seconds (after which the optimisation process is terminated), you can use the `maxtime` argument:
166+
```@example diffeq_param_estim_1
167+
optsol_fixed_kD = solve(optprob, Optim.NelderMead(); maxtime=100)
168+
nothing # hide
169+
```
170+
171+
---
172+
## References
173+
[^1]: [Alejandro F. Villaverde, Dilan Pathirana, Fabian Fröhlich, Jan Hasenauer, Julio R. Banga, *A protocol for dynamic model calibration*, Briefings in Bioinformatics (2023).](https://academic.oup.com/bib/article/23/1/bbab387/6383562?login=false)

‎docs/src/inverse_problems/parameter_estimation.md

Lines changed: 13 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,24 +1,16 @@
1-
# [Parameter Estimation](@id parameter_estimation)
2-
The parameters of a model, generated by Catalyst, can be estimated using various
3-
packages available in the Julia ecosystem. Refer
4-
[here](https://docs.sciml.ai/Overview/stable/highlevels/inverse_problems/) for
5-
more extensive information. Below follows a quick tutorial of how
6-
[Optimization.jl](https://docs.sciml.ai/Optimization/stable/) can be used to fit
7-
a parameter set to data.
1+
# [Fitting Parameters for an Oscillatory System](@id parameter_estimation)
2+
In this examples we will use [Optimization.jl](https://github.com/SciML/Optimization.jl) to fit the parameters of an oscillatory system (the Brusselator). 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. As we complete this procedure (which can be advantageous for oscillatory systems) we reach a global optimum.
83

94
First, we fetch the required packages.
105
```@example pe1
116
using Catalyst
12-
using DifferentialEquations
13-
using SciMLSensitivity
7+
using OrdinaryDiffEq
148
using Optimization
15-
16-
# for the ADAM optimizer
17-
using OptimizationOptimisers
9+
using OptimizationOptimisers # Required for the ADAM optimizer.
10+
using SciMLSensitivity # Required for `Optimization.AutoZygote()` automatic differentiation option.
1811
```
1912

20-
Next, we declare our model. For our example, we will use the Brusselator, a
21-
simple oscillator.
13+
Next, we declare our model, the Brusselator oscillator.
2214
```@example pe1
2315
brusselator = @reaction_network begin
2416
A, ∅ --> X
@@ -27,10 +19,11 @@ brusselator = @reaction_network begin
2719
1, X --> ∅
2820
end
2921
p_real = [:A => 1., :B => 2.]
22+
nothing # hide
3023
```
3124

3225
We simulate our model, and from the simulation generate sampled data points
33-
(with added noise), to which we will attempt to fit a parameter et.
26+
(to which we add noise). We will use this data to fit the parameters of our model.
3427
```@example pe1
3528
u0 = [:X => 1.0, :Y => 1.0]
3629
tspan = (0.0, 30.0)
@@ -54,8 +47,8 @@ scatter!(sample_times, sample_vals'; color = [:blue :red], legend = nothing)
5447

5548
Next, we create a function to fit the parameters using the `ADAM` optimizer. For
5649
a given initial estimate of the parameter values, `pinit`, this function will
57-
fit parameter values, `p`, to our data samples. `tend` is used to indicate the
58-
time interval over which to fit to the ODE solution.
50+
fit parameter values, `p`, to our data samples. We use `tend` to indicate the
51+
time interval over which we fit the model.
5952
```@example pe1
6053
function optimise_p(pinit, tend)
6154
function loss(p, _)
@@ -124,10 +117,10 @@ The final parameter estimate is then
124117
```@example pe1
125118
p_estimate
126119
```
127-
which is close to the actual parameters of `[1.0, 2.0]`.
120+
which is close to the actual parameter set of `[1.0, 2.0]`.
128121

129-
## Why we fit the parameters in iterations.
130-
The reason we chose to fit the model on a smaller interval to begin with, and
122+
## Why we fit the parameters in iterations
123+
As previously mentioned, the reason we chose to fit the model on a smaller interval to begin with, and
131124
then extend the interval, is to avoid getting stuck in a local minimum. Here
132125
specifically, we chose our initial interval to be smaller than a full cycle of
133126
the oscillation. If we had chosen to fit a parameter set on the full interval

0 commit comments

Comments
 (0)