Skip to content

Commit

Permalink
Merge pull request #1149 from SciML/change_optimization_parameter_fit…
Browse files Browse the repository at this point in the history
…ting_doc

Update Optimization/parameter fitting tutorial
TorkelE authored Dec 29, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
2 parents 1a70e28 + d4147d8 commit 4efe8ab
Showing 5 changed files with 182 additions and 159 deletions.
4 changes: 2 additions & 2 deletions docs/Project.toml
Original file line number Diff line number Diff line change
@@ -5,7 +5,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
DiffEqParamEstim = "1130ab10-4a5a-5621-a13d-e4788d82bd4c"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
@@ -23,6 +22,7 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
OptimizationEvolutionary = "cb963754-43f6-435e-8d4b-99009ff27753"
OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
@@ -50,7 +50,6 @@ CairoMakie = "0.12"
Catalyst = "14.4"
DataFrames = "1.6"
DiffEqBase = "6.159.0"
DiffEqParamEstim = "2.2"
Distributions = "0.25"
Documenter = "1.4.1"
DynamicalSystems = "3.3"
@@ -67,6 +66,7 @@ NonlinearSolve = "3.12, 4"
Optim = "1.9"
Optimization = "4"
OptimizationBBO = "0.4"
OptimizationEvolutionary = "0.4"
OptimizationNLopt = "0.3"
OptimizationOptimJL = "0.4"
OptimizationOptimisers = "0.3"
2 changes: 1 addition & 1 deletion docs/pages.jl
Original file line number Diff line number Diff line change
@@ -45,8 +45,8 @@ pages = Any[
"steady_state_functionality/dynamical_systems.md"
],
"Inverse problems" => Any[
"inverse_problems/optimization_ode_param_fitting.md",
"inverse_problems/petab_ode_param_fitting.md",
"inverse_problems/optimization_ode_param_fitting.md",
"inverse_problems/behaviour_optimisation.md",
# "inverse_problems/structural_identifiability.md",
"inverse_problems/global_sensitivity_analysis.md",
47 changes: 16 additions & 31 deletions docs/src/inverse_problems/behaviour_optimisation.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# [Optimization for non-data fitting purposes](@id behaviour_optimisation)
In previous tutorials we have described how to use PEtab.jl and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1].
In previous tutorials we have described how to use [PEtab.jl](@ref petab_parameter_fitting) and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom objective function, which minimizer can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. Many options used here are described in more detail in [the tutorial on using Optimization.jl for parameter fitting](@ref optimization_parameter_fitting). A more throughout description of how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1].

## [Maximising the pulse amplitude of an incoherent feed forward loop](@id behaviour_optimisation_IFFL_example)
Incoherent feedforward loops (network motifs where a single component both activates and deactivates a downstream component) are able to generate pulses in response to step inputs[^2]. In this tutorial we will consider such an incoherent feedforward loop, attempting to generate a system with as prominent a response pulse as possible.
@@ -17,40 +17,42 @@ end
```
To demonstrate this pulsing behaviour we will simulate the system for an example parameter set. We select an initial condition (`u0`) so the system begins in a steady state.
```@example behaviour_optimization
using OrdinaryDiffEqTsit5, Plots
using OrdinaryDiffEqDefault, Plots
example_p = [:pX => 0.1, :pY => 1.0, :pZ => 1.0]
tspan = (0.0, 50.0)
example_u0 = [:X => 0.1, :Y => 0.1, :Z => 1.0]
oprob = ODEProblem(incoherent_feed_forward, example_u0, tspan, example_p)
sol = solve(oprob, Tsit5())
sol = solve(oprob)
plot(sol)
```
Here we note that, while $X$ and $Y$ reach new steady state levels in response to the increase in $pX$, $Z$ resumes to its initial concentration after the pulse.

We will now attempt to find the parameter set $(pX,pY,pZ)$ which maximises the response pulse amplitude (defined by the maximum activity of $Z$ subtracted by its steady state activity). To do this, we create a custom cost function:
We will now attempt to find the parameter set $(pX,pY,pZ)$ which maximises the response pulse amplitude (defined by the maximum activity of $Z$ subtracted by its steady state activity). To do this, we create a custom objective function:
```@example behaviour_optimization
function pulse_amplitude(p, _)
ps = Dict([:pX => p[1], :pY => p[2], :pZ => p[2]])
u0_new = [:X => ps[:pX], :Y => ps[:pX]*ps[:pY], :Z => ps[:pZ]/ps[:pY]^2]
oprob_local = remake(oprob; u0= u0_new, p = ps)
sol = solve(oprob_local, Tsit5(); verbose = false, maxiters = 10000)
p = Dict([:pX => p[1], :pY => p[2], :pZ => p[2]])
u0 = [:X => p[:pX], :Y => p[:pX]*p[:pY], :Z => p[:pZ]/p[:pY]^2]
oprob_local = remake(oprob; u0, p)
sol = solve(oprob_local; verbose = false, maxiters = 10000)
SciMLBase.successful_retcode(sol) || return Inf
return -(maximum(sol[:Z]) - sol[:Z][1])
end
nothing # hide
```
This cost function takes two arguments (a parameter value `p`, and an additional one which we will ignore here but discuss later). It first calculates the new initial steady state concentration for the given parameter set. Next, it creates an updated `ODEProblem` using the steady state as initial conditions and the, to the cost function provided, input parameter set. While we could create a new `ODEProblem` within the cost function, cost functions are often called a large number of times during the optimisation process (making performance important). Here, using [`remake` on a previously created `ODEProblem`](@ref simulation_structure_interfacing_problems_remake) is more performant than creating a new one. Just like [when using Optimization.jl to fit parameters to data](@ref optimization_parameter_fitting), we use the `verbose = false` option to prevent unnecessary simulation printouts, and a reduced `maxiters` value to reduce time spent simulating (for the model) unsuitable parameter sets. We also use `SciMLBase.successful_retcode(sol)` to check whether the simulation return code indicates a successful simulation (and if it did not, returns a large cost function value). Finally, Optimization.jl finds the function's *minimum value*, so to find the *maximum* relative pulse amplitude, we make our cost function return the negative pulse amplitude.
This objective function takes two arguments (a parameter value `p`, and an additional one which we will ignore but is discussed in a note [here](@ref optimization_parameter_fitting_basics)). It first calculates the new initial steady state concentration for the given parameter set. Next, it creates an updated `ODEProblem` using the steady state as initial conditions and the, to the objective function provided, input parameter set. Finally, Optimization.jl finds the function's *minimum value*, so to find the *maximum* relative pulse amplitude, we make our objective function return the negative pulse amplitude.

Just like for [parameter fitting](@ref optimization_parameter_fitting), we create a `OptimizationProblem` using our cost function, and some initial guess of the parameter value. We also set upper and lower bounds for each parameter using the `lb` and `ub` optional arguments (in this case limiting each parameter's value to the interval $(0.1,10.0)$).
As described [in our tutorial on parameter fitting using Optimization.jl](@ref optimization_parameter_fitting_basics) we use `remake`, `verbose = false`, `maxiters = 10000`, and a check on the simulations return code, all providing various advantages to the optimisation procedure (as explained in that tutorial).

Just like for [parameter fitting](@ref optimization_parameter_fitting_basics), we create an `OptimizationProblem` using our objective function, and some initial guess of the parameter values. We also [set upper and lower bounds](@ref optimization_parameter_fitting_constraints) for each parameter using the `lb` and `ub` optional arguments (in this case limiting each parameter's value to the interval $(0.1,10.0)$).
```@example behaviour_optimization
using Optimization
initial_guess = [1.0, 1.0, 1.0]
opt_prob = OptimizationProblem(pulse_amplitude, initial_guess; lb = [1e-1, 1e-1, 1e-1], ub = [1e1, 1e1, 1e1])
nothing # hide
```
!!! note
As described in a [previous section on Optimization.jl](@ref optimization_parameter_fitting), `OptimizationProblem`s do not support setting parameter values using maps. We must instead set `initial_guess` values using a vector. Next, in the first line of our cost function, we reshape the parameter values to the common form used across Catalyst (e.g. `[:pX => p[1], :pY => p[2], :pZ => p[2]]`, however, here we use a dictionary to easier compute the steady state initial condition). We also note that the order used in this array corresponds to the order we give each parameter's bounds in `lb` and `ub`, and the order in which their values occur in the output solution.
As described in a [previous section on Optimization.jl](@ref optimization_parameter_fitting), `OptimizationProblem`s do not support setting parameter values using maps. We must instead set `initial_guess` values using a vector. Next, in the first line of our objective function, we reshape the parameter values to the common form used across Catalyst (e.g. `[:pX => p[1], :pY => p[2], :pZ => p[2]]`, however, here we use a dictionary to easier compute the steady state initial condition). We also note that the order used in this array corresponds to the order we give each parameter's bounds in `lb` and `ub`, and the order in which their values occur in the output solution.

As [previously described](@ref optimization_parameter_fitting), Optimization.jl supports a wide range of optimisation algorithms. Here we use one from [BlackBoxOptim.jl](https://github.com/robertfeldt/BlackBoxOptim.jl):
```@example behaviour_optimization
@@ -63,30 +65,13 @@ Finally, we plot a simulation using the found parameter set (stored in `opt_sol.
ps_res = Dict([:pX => opt_sol.u[1], :pY => opt_sol.u[2], :pZ => opt_sol.u[2]])
u0_res = [:X => ps_res[:pX], :Y => ps_res[:pX]*ps_res[:pY], :Z => ps_res[:pZ]/ps_res[:pY]^2]
oprob_res = remake(oprob; u0 = u0_res, p = ps_res)
sol_res = solve(oprob_res, Tsit5())
sol_res = solve(oprob_res)
plot(sol_res; idxs = :Z)
```
For this model, it turns out that $Z$'s maximum pulse amplitude is equal to twice its steady state concentration. Hence, the maximisation of its pulse amplitude is equivalent to maximising its steady state concentration.

!!! note
Especially if you check Optimization.jl's documentation, you will note that cost functions have the `f(u,p)` form. This is because `OptimizationProblem`s (like e.g. `ODEProblem`s) can take both variables (which can be varied in the optimisation problem), but also parameters that are fixed. In our case, the *optimisation variables* correspond to our *model parameters*. Hence, our model parameter values are the `u` input. This is also why we find the optimisation solution (our optimised parameter set) in `opt_sol`'s `u` field. Our optimisation problem does not actually have any parameters, hence, the second argument of `pulse_amplitude` is unused (that is why we call it `_`, a name commonly indicating unused function arguments).

There are several modifications to our problem where it would actually have parameters. E.g. our model might have had additional parameters (e.g. a degradation rate) which we would like to keep fixed throughout the optimisation process. If we then would like to run the optimisation process for several different values of these fixed parameters, we could have made them parameters to our `OptimizationProblem` (and their values provided as a third argument, after `initial_guess`).

## [Utilising automatic differentiation](@id behaviour_optimisation_AD)
Optimisation methods can be divided into differentiation-free and differentiation-based optimisation methods. E.g. consider finding the minimum of the function $f(x) = x^2$, given some initial guess of $x$. Here, we can simply compute the differential and descend along it until we find $x=0$ (admittedly, for this simple problem the minimum can be computed directly). This principle forms the basis of optimisation methods such as [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent), which utilises information of a function's differential to minimise it. When attempting to find a global minimum, to avoid getting stuck in local minimums, these methods are often augmented by additional routines. While the differentiation of most algebraic functions is trivial, it turns out that even complicated functions (such as the one we used above) can be differentiated computationally through the use of [*automatic differentiation* (AD)](https://en.wikipedia.org/wiki/Automatic_differentiation).

Through packages such as [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl), [ReverseDiff.jl](https://github.com/JuliaDiff/ReverseDiff.jl), and [Zygote.jl](https://github.com/FluxML/Zygote.jl), Julia supports AD for most code. Specifically for code including simulation of differential equations, differentiation is supported by [SciMLSensitivity.jl](https://github.com/SciML/SciMLSensitivity.jl). Generally, AD can be used without specific knowledge from the user, however, it requires an additional step in the construction of our `OptimizationProblem`. Here, we create a [specialised `OptimizationFunction` from our cost function](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#optfunction). To it, we will also provide our choice of AD method. There are [several alternatives](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Automatic-Differentiation-Construction-Choice-Recommendations), and in our case we will use `AutoForwardDiff()` (a good choice for small optimisation problems). We can then create a new `OptimizationProblem` using our updated cost function:
```@example behaviour_optimization
opt_func = OptimizationFunction(pulse_amplitude, AutoForwardDiff())
opt_prob = OptimizationProblem(opt_func, initial_guess; lb = [1e-1, 1e-1, 1e-1], ub = [1e1, 1e1, 1e1])
nothing # hide
```
Finally, we can find the optimum using some differentiation-based optimisation methods. Here we will use [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s `BFGS` method:
```julia
using OptimizationOptimJL
opt_sol = solve(opt_prob, OptimizationOptimJL.BFGS())
```
## [Other optimisation options](@id behaviour_optimisation_options)
How to use Optimization.jl is discussed in more detail in [this tutorial](@ref optimization_parameter_fitting). This includes options such as using [automatic differentiation](@ref optimization_parameter_fitting_AD), [setting constraints](@ref optimization_parameter_fitting_constraints), and setting [optimisation solver options](@ref optimization_parameter_fitting_solver_options). Finally, it discusses the advantages of [carrying out the fitting in logarithmic space](@ref optimization_parameter_fitting_log_scale), something which can be advantageous for the problem described above as well.

---
## [Citation](@id structural_identifiability_citation)
200 changes: 118 additions & 82 deletions docs/src/inverse_problems/optimization_ode_param_fitting.md

Large diffs are not rendered by default.

88 changes: 45 additions & 43 deletions docs/src/inverse_problems/petab_ode_param_fitting.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
# [Parameter Fitting for ODEs using PEtab.jl](@id petab_parameter_fitting)
The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data [^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data and is a good default choice when fitting reaction rate equation ODE models. This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation).
The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data [^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data and is a good default choice when fitting reaction rate equation ODE models. This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation).

While PEtab's interface generally is very flexible, there might be specific use-cases where it cannot create an appropriate cost-function. Here, it is recommended to instead look at using [Optimization.jl](@ref optimization_parameter_fitting).

## Introductory example
Let us consider a simple catalysis network, where an enzyme ($E$) turns a substrate ($S$) into a product ($P$):
@@ -19,18 +21,18 @@ u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
# Generate synthetic data.
using OrdinaryDiffEqTsit5
using OrdinaryDiffEqDefault
oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true)
true_sol = solve(oprob_true, Tsit5())
data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
true_sol = solve(oprob_true)
data_sol = solve(oprob_true; saveat = 1.0)
data_ts = data_sol.t[2:end]
data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
# Plots the true solutions and the (synthetic data) measurements.
using Plots
default(bottom_margin=4Plots.Measures.mm,left_margin=4Plots.Measures.mm) # hide
plot(true_sol; idxs=:P, label="True solution", lw=4)
plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue)
default(bottom_margin = 4Plots.Measures.mm, left_margin = 4Plots.Measures.mm) # hide
plot(true_sol; idxs = :P, label = "True solution", lw = 4)
plot!(data_ts, data_vals; label = "Measurements", seriestype = :scatter, ms = 6, color = :blue)
```

Generally, PEtab takes five different inputs to define an optimisation problem (the minimiser of which generates a fitted parameter set):
@@ -86,15 +88,15 @@ For cases where several simulation conditions are given, we also need to provide
Each individual measurement is provided as a row of a `DataFrame`. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three:
```@example petab1
using DataFrames
measurements = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals)
measurements = DataFrame(obs_id = "obs_P", time = data_ts, measurement = data_vals)
```
Since, in our example, all measurements are of the same observable, we can set `obs_id="obs_P"`. However, it is also possible to [include measurements from several different observables](@ref petab_simulation_measurements_several_observables).


### Creating a PEtabModel
Finally, we combine all inputs in a single `PEtabModel`. To it, we also pass the initial conditions of our simulations (using the `speciemap` argument). It is also possible to have [initial conditions with uncertainty](@ref petab_simulation_initial_conditions_uncertainty), [that vary between different simulations](@ref petab_simulation_conditions), or [that we attempt to fit to the data](@ref petab_simulation_initial_conditions_fitted).
```@example petab1
petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0)
petab_model = PEtabModel(rn, observables, measurements, params; speciemap = u0)
nothing # hide
```

@@ -116,9 +118,9 @@ res = calibrate(petab_problem, p0, IPNewton())

We can now simulate our model for the fitted parameter set, and compare the result to the measurements and true solution.
```@example petab1
oprob_fitted = remake(oprob_true; p=get_ps(res, petab_problem))
fitted_sol = solve(oprob_fitted, Tsit5())
plot!(fitted_sol; idxs=:P, label="Fitted solution", linestyle=:dash, lw=4, color=:lightblue)
oprob_fitted = remake(oprob_true; p = get_ps(res, petab_problem))
fitted_sol = solve(oprob_fitted)
plot!(fitted_sol; idxs = :P, label = "Fitted solution", linestyle = :dash, lw = 4, color = :lightblue)
```

Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/stable/API/), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note, that even if PEtab.jl has found the global optimum (which fits the data well), this does not actually correspond to the true parameter set. This phenomenon is related to the concept of *identifiability*, which is very important for parameter fitting.
@@ -155,7 +157,7 @@ It would also be possible to make $σ$ a *parameter that is fitted as a part of

PEtab.jl assumes that noise is normally distributed (with a standard deviation equal to the second argument provided to `PEtabObservable`). The only other (currently implemented) noise distribution is log-normally distributed noise, which is designated through the `transformation=:log` argument:
```@example petab1
obs_P = PEtabObservable(P, σ; transformation=:log)
obs_P = PEtabObservable(P, σ; transformation = :log)
```

## [Additional features: Parameters](@id petab_parameters)
@@ -170,14 +172,14 @@ nothing # hide
```
We then provide `parameter_map=[:kB => 1.0]` when we assembly our model:
```@example petab1
petab_model_known_param = PEtabModel(rn, observables, measurements, params; speciemap=u0, parametermap=[:kB => 1.0])
petab_model_known_param = PEtabModel(rn, observables, measurements, params; speciemap = u0, parametermap = [:kB => 1.0])
nothing # hide
```

### [Parameter bounds](@id petab_parameters_bounds)
By default, when fitted, potential parameter values are assumed to be in the interval $(1e-3, 1e3)$. When declaring a `PEtabParameter` it is possible to change these values through the `lb` and `ub` arguments. E.g. we could use
```@example petab1
par_kB = PEtabParameter(:kB; lb=1e-2, ub=1e2)
par_kB = PEtabParameter(:kB; lb = 1e-2, ub = 1e2)
```
to achieve the more conservative bound $(1e-2, 1e2)$ for the parameter $kB$.

@@ -186,19 +188,19 @@ to achieve the more conservative bound $(1e-2, 1e2)$ for the parameter $kB$.
Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous[^2]). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of $kB$ to linear:

```@example petab1
par_kB = PEtabParameter(:kB; scale=:lin)
par_kB = PEtabParameter(:kB; scale = :lin)
```

### [Parameter priors](@id petab_parameters_priors)
If we have prior knowledge about the distribution of a parameter, it is possible to incorporate this into the model. The prior can be any continuous, univariate, distribution from the [Distributions.jl package](https://github.com/JuliaStats/Distributions.jl). E.g we can use:

```@example petab1
using Distributions
par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2))
par_kB = PEtabParameter(:kB; prior = Normal(1.0,0.2))
```
to set a normally distributed prior (with mean `1.0` and standard deviation `0.2`) on the value of $kB$. By default, the prior is assumed to be on the linear scale of the parameter (before any potential log transform). To specify that the prior is on the logarithmic scale, the `prior_on_linear_scale=false` argument can be provided:
```@example petab1
par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2), prior_on_linear_scale=false)
par_kB = PEtabParameter(:kB; prior = Normal(1.0,0.2), prior_on_linear_scale = false)
```
In this example, setting `prior_on_linear_scale=false` makes sense as a (linear) normal distribution is non-zero for negative values (an alternative is to use a log-normal distribution, e.g. `prior=LogNormal(3.0, 3.0)`).

@@ -232,15 +234,15 @@ u0 = [:E => 1.0, :SE => 0.0, :P => 0.0]
p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
# Simulate data.
using OrdinaryDiffEqTsit5
using OrdinaryDiffEqDefault
t1, d1 = let
oprob_true = ODEProblem(rn, [:S=>1.0; u0], (0.0, 10.0), p_true)
data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
oprob_true = ODEProblem(rn, [:S => 1.0; u0], (0.0, 10.0), p_true)
data_sol = solve(oprob_true; saveat = 1.0)
data_sol.t[2:end], (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
end
t2, d2 = let
oprob_true = ODEProblem(rn, [:S=>0.5; u0], (0.0, 10.0), p_true)
data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
data_sol = solve(oprob_true; saveat = 1.0)
data_sol.t[2:end], (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
end
@@ -258,11 +260,11 @@ c2 = Dict(:S => 0.5)
simulation_conditions = Dict("c1" => c1, "c2" => c2)
using DataFrames
m1 = DataFrame(simulation_id="c1", obs_id="obs_P", time=t1, measurement=d1)
m2 = DataFrame(simulation_id="c2", obs_id="obs_P", time=t2, measurement=d2)
m1 = DataFrame(simulation_id = "c1", obs_id = "obs_P", time = t1, measurement = d1)
m2 = DataFrame(simulation_id = "c2", obs_id = "obs_P", time = t2, measurement = d2)
measurements = vcat(m1,m2)
petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0,
petab_model = PEtabModel(rn, observables, measurements, params; speciemap = u0,
simulation_conditions = simulation_conditions)
nothing # hide
```
@@ -288,8 +290,8 @@ nothing # hide
```
we are able to include all these measurements in the same `measurements` `DataFrame`:
```@example petab1
m1 = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals_S)
m2 = DataFrame(obs_id="obs_S", time=data_ts, measurement=data_vals_P)
m1 = DataFrame(obs_id = "obs_P", time = data_ts, measurement = data_vals_S)
m2 = DataFrame(obs_id = "obs_S", time = data_ts, measurement = data_vals_P)
measurements = vcat(m1,m2)
```
which then can be used as input to `PEtabModel`.
@@ -346,10 +348,10 @@ end
u0 = [:SE => 0.0, :P => 0.0]
p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5, :S0=>1.0, :E0 => 1.0]
using OrdinaryDiffEqTsit5
using OrdinaryDiffEqDefault
oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true)
true_sol = solve(oprob_true, Tsit5())
data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
true_sol = solve(oprob_true)
data_sol = solve(oprob_true; saveat = 1.0)
data_ts = data_sol.t[2:end]
data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
@@ -365,9 +367,9 @@ par_E0 = PEtabParameter(:E0; prior = Normal(1.0, 0.1))
params = [par_kB, par_kD, par_kP, par_S0, par_E0]
using DataFrames
measurements = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals)
measurements = DataFrame(obs_id = "obs_P", time = data_ts, measurement = data_vals)
petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0)
petab_model = PEtabModel(rn, observables, measurements, params; speciemap = u0)
nothing # hide
```
Here, when we fit our data we will also gain values for `S0` and `E0`, however, unless we are explicitly interested in these, they can be ignored.
@@ -380,9 +382,9 @@ Here is an example, adapted from the [more detailed PEtab.jl documentation](http
```@example petab1
using OrdinaryDiffEqRosenbrock
PEtabODEProblem(petab_model,
odesolver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8),
gradient_method=:ForwardDiff,
hessian_method=:ForwardDiff)
odesolver = ODESolver(Rodas5P(), abstol = 1e-8, reltol = 1e-8),
gradient_method = :ForwardDiff,
hessian_method = :ForwardDiff)
nothing # hide
```
where we simulate our ODE model using the `Rodas5P` method (with absolute and relative tolerance both equal `1e-8`) and use [forward automatic differentiation](https://github.com/JuliaDiff/ForwardDiff.jl) for both gradient and hessian computation. More details on available ODE solver options can be found in [the PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/API/#PEtabODEProblem).
@@ -392,7 +394,7 @@ where we simulate our ODE model using the `Rodas5P` method (with absolute and re
### [Optimisation methods and options](@id petab_optimisation_optimisers)
For our examples, we have used the `Optim.IPNewton` optimisation method. PEtab.jl supports [several additional optimisation methods](https://sebapersson.github.io/PEtab.jl/stable/pest_algs/). Furthermore, `calibrate`'s `options` argument permits the customisation of the options for any used optimiser. E.g. to designate the maximum number of iterations of the `Optim.IPNewton` method we would use:
```@example petab1
res = calibrate(petab_problem, p0, IPNewton(); options=Optim.Options(iterations = 10000))
res = calibrate(petab_problem, p0, IPNewton(); options = Optim.Options(iterations = 10000))
nothing # hide
```

@@ -402,7 +404,7 @@ Please read the [PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/

To record all the parameter sets evaluated (and their objective values) during the optimisation procedure, provide the `save_trace=true` argument to `calibrate` (or `calibrate_multistart`):
```@example petab1
res = calibrate(petab_problem, p0, IPNewton(); save_trace=true)
res = calibrate(petab_problem, p0, IPNewton(); save_trace = true)
nothing # hide
```
This is required for the various [optimisation evaluation plots](@ref petab_plotting) provided by PEtab.jl. If desired, this information can be accessed in the calibration output's `.xtrace` and `.ftrace` fields.
@@ -439,8 +441,8 @@ using Optim
using QuasiMonteCarlo
mkdir("OptimisationRuns") # hide
res_ms = calibrate_multistart(petab_problem, IPNewton(), 10; dirsave = "OptimisationRuns",
sampling_method=QuasiMonteCarlo.SobolSample())
res_ms = calibrate_multistart(petab_problem, IPNewton(), 10; dirsave = "OptimisationRuns", sampling_method=QuasiMonteCarlo.SobolSample()) # hide
sampling_method = QuasiMonteCarlo.SobolSample())
res_ms = calibrate_multistart(petab_problem, IPNewton(), 10; dirsave = "OptimisationRuns", sampling_method = QuasiMonteCarlo.SobolSample()) # hide
nothing # hide
```
The best result across all runs can still be retrieved using `get_ps(res_ms, petab_problem)`, with the results of the individual runs being stored in the `res_ms.runs` field.
@@ -452,8 +454,8 @@ nothing # hide
```
where `"OptimisationRuns"` is the name of the save directory (specified in `calibrate_multistart`). If the OptimisationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we load the second run to be saved in that folder:
```@example petab1
res_ms = PEtabMultistartResult("OptimisationRuns"; which_run=2)
rm("OptimisationRuns", recursive=true) # hide
res_ms = PEtabMultistartResult("OptimisationRuns"; which_run = 2)
rm("OptimisationRuns", recursive = true) # hide
nothing # hide
```
By default, `which_run` loads the first run saved to that directory.
@@ -479,7 +481,7 @@ Whenever we have several events or not, we bundle them together in a single vect
```@example petab1
params = [par_kB, par_kD, par_kP] # hide
events = [event1, event2]
petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0, events=events)
petab_model = PEtabModel(rn, observables, measurements, params; speciemap = u0, events = events)
nothing # hide
```

@@ -493,7 +495,7 @@ There exist various types of graphs that can be used to evaluate the parameter f

To, for a single start calibration run, plot, for each iteration of the optimization process, the best objective value achieved so far, run:
```@example petab1
res = calibrate(petab_problem, p0, IPNewton(); save_trace=true) # hide
res = calibrate(petab_problem, p0, IPNewton(); save_trace = true) # hide
plot(res)
```

0 comments on commit 4efe8ab

Please sign in to comment.