From 492370f40f40f0f572e7b531227bd29124800f2f Mon Sep 17 00:00:00 2001 From: sebapersson Date: Fri, 15 Nov 2024 14:11:34 +0100 Subject: [PATCH] Update PEtab tutorial --- docs/Project.toml | 2 + docs/pages.jl | 2 +- .../petab_ode_param_fitting.md | 138 +++++++++--------- 3 files changed, 68 insertions(+), 74 deletions(-) rename docs/{unpublished => src/inverse_problems}/petab_ode_param_fitting.md (80%) diff --git a/docs/Project.toml b/docs/Project.toml index a18db09ea5..c0bee6b155 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -33,6 +33,7 @@ OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2" +PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" @@ -77,6 +78,7 @@ OrdinaryDiffEqRosenbrock = "1" OrdinaryDiffEqSDIRK = "1" OrdinaryDiffEqTsit5 = "1" OrdinaryDiffEqVerner = "1" +PEtab = "3.5" Plots = "1.40" QuasiMonteCarlo = "0.3" SciMLBase = "2.46" diff --git a/docs/pages.jl b/docs/pages.jl index c72f578bbb..8dcc001a41 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -46,7 +46,7 @@ pages = Any[ ], "Inverse problems" => Any[ "inverse_problems/optimization_ode_param_fitting.md", - # "inverse_problems/petab_ode_param_fitting.md", + "inverse_problems/petab_ode_param_fitting.md", "inverse_problems/behaviour_optimisation.md", # "inverse_problems/structural_identifiability.md", "inverse_problems/global_sensitivity_analysis.md", diff --git a/docs/unpublished/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md similarity index 80% rename from docs/unpublished/petab_ode_param_fitting.md rename to docs/src/inverse_problems/petab_ode_param_fitting.md index a57e431912..91833e0155 100644 --- a/docs/unpublished/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -1,7 +1,6 @@ # [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). - ## Introductory example Let us consider a simple catalysis network, where an enzyme ($E$) turns a substrate ($S$) into a product ($P$): ```@example petab1 @@ -30,16 +29,16 @@ 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=8) +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): 1. **Model**: The model which we want to fit to the data (a `ReactionSystem`). 2. **Observables**: The possible observables that can be measured (a `Dict{String,PEtabObservable}`). -3. **Estimation parameters**: The parameters which we want to fit to the data (a `Vector{PEtabParameter}`). +3. **Estimation parameters**: The parameters which we want to fit to the data (a `Vector{PEtabParameter}`). 4. **Experimental (or simulation) conditions**: The simulations (each corresponding to a potential experiment) carried out during each step of the optimisation process (a `Dict{String,Dict{Symbol,Float64}}`). -5. **Measurements**: The measurements to which the model is fitted (a `DataFrame`). +5. **Measurements**: The measurements to which the model is fitted (a `DataFrame`). ### Observables The observables define the quantities that we may measure in our experiments. Typically, each corresponds to a single species, however, [more complicated observables are possible](@ref petab_observables_observables). For each observable, we also need a noise formula, defining the uncertainty in its measurements. By default, PEtab assumes normally distributed noise, with a mean equal to the true value and a standard deviation which we have to define. It is also possible to use [more advanced noise formulas](@ref petab_observables_noise_formula). @@ -56,7 +55,7 @@ nothing # hide ### Parameters Each parameter of the system can either be 1. Known (described [here](@ref petab_parameters_known)). -2. Depend on experimental/simulation conditions (described [here](@ref petab_simulation_conditions)). +2. Depend on experimental/simulation conditions (described [here](@ref petab_simulation_conditions)). 3. Be an unknown that we wish to fit to data. In our case, we wish to fit all three system parameters ($kB$, $kD$, and $kP$). For each, we create a single `PEtabParameter`, and then gather these into a single vector. @@ -72,13 +71,12 @@ For each parameter, it is also possible to set [a lower and/or upper bound](@ref ### Simulation conditions Sometimes, several different experiments are performed on a system (each potentially generating several measurements). An experiment could e.g. be the time development of a system from a specific initial condition. Since each experimental condition (during the optimisation procedure, for a guess of the unknown parameters) generates a distinct simulation, these are also called simulation conditions. In our example, all data comes from a single experiment, and the simulation condition input is not required. How to define and use different experimental conditions is described [here](@ref petab_simulation_conditions). - ### Measurements Finally, we need to define the system measurements to which the parameters will be fitted. Each measurement combines: 1. The observable which is observed (here we use the id tag defined in the `observables` dictionary). 2. The time point of the measurement. 3. The measured value. - + For cases where several simulation conditions are given, we also need to provide: 4. The simulation condition which generates the measurement ([here](@ref petab_simulation_conditions) is an example where this is used). @@ -94,43 +92,39 @@ Since, in our example, all measurements are of the same observable, we can set ` ### 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 `state_map` 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). +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; state_map=u0) +petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0) nothing # hide ``` ### Fitting parameters We are now able to fit our model to the data. First, we create a `PEtabODEProblem`. Here, we use `petab_model` as the only input, but it is also possible to set various [numeric solver and automatic differentiation options](@ref petab_simulation_options) (such as method or tolerance). ```@example petab1 -petab_problem = PEtabODEProblem(petab_model; verbose=false) # hide -nothing # hide -``` -```julia petab_problem = PEtabODEProblem(petab_model) ``` Since no additional input is given, default options are selected by PEtab.jl (and generally, its choices are good). -To fit a parameter set we use the `calibrate_model` function. In addition to our `PEtabODEProblem`, we must also provide an initial guess (which can be generated with the `generate_startguesses` function) and an optimisation algorithm (which needs to be imported specifically). PEtab.jl supports various [optimisation methods and options](@ref petab_optimisation_optimisers). +To fit a parameter set we use the `calibrate` function. In addition to our `PEtabODEProblem`, we must also provide an initial guess (which can be generated with the `generate_startguesses` function) and an optimisation algorithm (which needs to be imported specifically). PEtab.jl supports various [optimisation methods and options](@ref petab_optimisation_optimisers). ```@example petab1 using Optim -p0 = generate_startguesses(petab_problem, 1) +p0 = get_startguesses(petab_problem, 1) p0 = [0.0, 0.0, 0.0] # hide -res = calibrate_model(petab_problem, p0, IPNewton()) # hide -res = calibrate_model(petab_problem, p0, IPNewton()) +res = calibrate(petab_problem, p0, IPNewton()) # hide +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=6, color=:lightblue) +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/dev/API_choosen/#PEtab.get_odeproblem), 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*](@ref structural_identifiability), which is very important for parameter fitting. +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. ### Final notes -PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options that are not specified). +PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/petab_preeq_simulations/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options that are not specified). ## [Additional features: Observables](@id petab_observables) @@ -145,7 +139,7 @@ end obs_X = PEtabObservable(X1 + X2, 0.5) ``` -A common application for this is to define an [*offset* and a *scale* for each observable](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/#Defining-the-Observable). +A common application for this is to define an [*offset* and a *scale* for each observable](https://sebapersson.github.io/PEtab.jl/stable/petab_obs_noise/). ### [Advanced observables noise formulas](@id petab_observables_noise_formula) In our basic example we assumed that the normally distributed noise had a standard deviation of `0.5`. However, this value may be a parameter (or indeed any algebraic expression). E.g, we could set @@ -176,7 +170,7 @@ 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; state_map=u0, parameter_map=[:kB => 1.0]) +petab_model_known_param = PEtabModel(rn, observables, measurements, params; speciemap=u0, parametermap=[:kB => 1.0]) nothing # hide ``` @@ -209,7 +203,7 @@ 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)`). ## [Simulation conditions](@id petab_simulation_conditions) -Sometimes, we have data from several different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define several different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. +Sometimes, we have data from different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. If, for our previous catalysis example, we had measured the system for two different initial values of $S$ ($S(0)=1.0$ and $S(0)=\tfrac{1}{2}$), these would correspond to two different simulation conditions. For each condition we define a `Dict` mapping the species to their initial condition (here, $S$ is the only species in each `Dict`): ```@example petab1 @@ -239,12 +233,12 @@ p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] # Simulate data. using OrdinaryDiffEqTsit5 -t1, d1 = let +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) data_sol.t[2:end], (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] end -t2, d2 = let +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.t[2:end], (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] @@ -259,8 +253,8 @@ par_kD = PEtabParameter(:kD) par_kP = PEtabParameter(:kP) params = [par_kB, par_kD, par_kP] -c1 = Dict(:kB => 1.0) -c2 = Dict(:kB => 0.5) +c1 = Dict(:S => 1.0) +c2 = Dict(:S => 0.5) simulation_conditions = Dict("c1" => c1, "c2" => c2) using DataFrames @@ -268,10 +262,11 @@ 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, simulation_conditions, observables, measurements, params; state_map=u0) +petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0, + simulation_conditions = simulation_conditions) nothing # hide ``` -Note that the `u0` we pass into `PEtabModel` through the `state_map` argument no longer contains the value of $S$ (as it is provided by the conditions). +Note that the `u0` we pass into `PEtabModel` through the `speciemap` argument no longer contains the value of $S$ (as it is provided by the conditions). ## [Additional features: Measurements](@id petab_simulation_measurements) @@ -300,12 +295,12 @@ measurements = vcat(m1,m2) which then can be used as input to `PEtabModel`. ### Varying parameters between different simulation conditions -Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates ($kP_1$ and $kP_2$) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters $kP_1$ and $kP_2$ are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_condition_specific/) provided instructions of how to handle them. +Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates ($kP_1$ and $kP_2$) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters $kP_1$ and $kP_2$ are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/stable/petab_cond_specific/) provided instructions of how to handle them. ## [Additional features: Initial conditions](@id petab_simulation_initial_conditions) ### [Fitting initial conditions](@id petab_simulation_initial_conditions_fitted) -Sometimes, initial conditions are uncertain quantities which we wish to fit to the data. This is possible by [defining an initial condition as a parameter](@ref dsl_description_parametric_initial_conditions): +Sometimes, initial conditions are uncertain quantities which we wish to fit to the data. This is possible by [defining an initial condition as a parameter](@ref dsl_description_parametric_initial_conditions): ```@example petab4 using Catalyst, PEtab # hide rn = @reaction_network begin @@ -334,11 +329,10 @@ nothing # hide and we can use our updated `rn`, `u0`, and `params` as input to our `PEtabModel`. ### [Uncertain initial conditions](@id petab_simulation_initial_conditions_uncertainty) -Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attaching a prior to it corresponding to our certainty about its value. +Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attaching a prior to it corresponding to our certainty about its value. Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values $S(0) = E(0) = 1.0$. ```@example petab5 - using Catalyst, PEtab rn = @reaction_network begin @@ -352,7 +346,7 @@ 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 OrdinaryDiffEq +using OrdinaryDiffEqTsit5 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) @@ -366,51 +360,49 @@ observables = Dict("obs_P" => obs_P) par_kB = PEtabParameter(:kB) par_kD = PEtabParameter(:kD) par_kP = PEtabParameter(:kP) -par_S0 = PEtabParameter(:S0) -par_E0 = PEtabParameter(:E0) +par_S0 = PEtabParameter(:S0; prior = Normal(1.0, 0.1)) +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) -petab_model = PEtabModel(rn, observables, measurements, params; state_map=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. ## [Additional features: Simulation options](@id petab_simulation_options) -While in our basic example, we do not provide any additional information to our `PEtabODEProblem`, this is an opportunity to specify how the model should be simulated, and what automatic differentiation techniques to use for the optimisation procedure (if none are provided, appropriate defaults are selected). +While in our basic example, we do not provide any additional information to our `PEtabODEProblem`, this is an opportunity to specify how the model should be simulated, and what automatic differentiation techniques to use for the optimisation procedure (if none are provided, appropriate defaults are selected). -Here is an example, taken from the [more detailed PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/Boehm/#Creating-a-PEtabODEProblem) +Here is an example, adapted from the [more detailed PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/default_options/) ```@example petab1 using OrdinaryDiffEqRosenbrock -PEtabODEProblem(petab_model, ode_solver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8), gradient_method=:ForwardDiff, hessian_method=:ForwardDiff, verbose=false); nothing # hide -``` -```julia PEtabODEProblem(petab_model, - ode_solver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8), + 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/dev/API_choosen/#PEtab.ODESolver). +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). ## [Additional features: Optimisation](@id petab_optimisation) ### [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/dev/Avaible_optimisers/#options_optimisers). Furthermore, `calibrate_model`'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: +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_model(petab_problem, p0, IPNewton(); options=Optim.Options(iterations = 10000)) +res = calibrate(petab_problem, p0, IPNewton(); options=Optim.Options(iterations = 10000)) nothing # hide ``` -Please read the [PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/Avaible_optimisers/#options_optimisers) to learn how to customise the various optimisers' properties. +Please read the [PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/pest_algs/) to learn how to customize the various optimisers' properties. ### [Optimisation path recording](@id petab_optimisation_path_recording) -To record all the parameter sets evaluated (and their objective values) during the optimisation procedure, provide the `save_trace=true` argument to `calibrate_model` (or `calibrate_model_multistart`): +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_model(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. @@ -418,54 +410,54 @@ This is required for the various [optimisation evaluation plots](@ref petab_plot ## Objective function extraction While PEtab.jl provides various tools for analysing the objective function generated by `PEtabODEProblem`, it is also possible to extract this function for customised analysis. Given a `PEtabODEProblem` ```@example petab1 -petab_problem = PEtabODEProblem(petab_model; verbose=false); nothing # hide +petab_problem = PEtabODEProblem(petab_model) +nothing # hide ``` ```julia petab_problem = PEtabODEProblem(petab_model) ``` We can find the: -1. Objective function as the `petab_problem.compute_cost`. It takes a single argument (`p`) and returns the objective value. -2. Gradient as the `petab_problem.compute_gradient!` field. It takes two arguments (`g` and `p`) with the updated gradient values being written to `g`. -3. Hessian as the `petab_problem.compute_hessian` field. It takes two arguments (`H` and `p`) with the updated hessian values being written to `H`. - +1. Objective function (negative log-likelihood) as the `petab_problem.nllh`. It takes a single argument (`p`) and returns the objective value. +2. Gradient as the `petab_problem.grad!` field. It takes two arguments (`g` and `p`) with the updated gradient values being written to `g`. +3. Hessian as the `petab_problem.hess!` field. It takes two arguments (`H` and `p`) with the updated hessian values being written to `H`. ## [Multi-start optimisation](@id petab_multistart_optimisation) -To avoid the optimisation process returning a local minimum, it is often advised to run it multiple times, using different initial guesses. PEtab.jl supports this through the `calibrate_model_multistart` function. This is identical to the `calibrate_model` function, but takes two additional arguments: +To avoid the optimisation process returning a local minimum, it is often advised to run it multiple times, using different initial guesses. PEtab.jl supports this through the `calibrate_multistart` function. This is identical to the `calibrate` function, but takes one additional arguments: -1. `n_multistart`: The number of runs to perform. -2. `dir_save`: A location to which the output is automatically saved. If `dir_save=nothing`, no output is saved. +1. `nmultistarts`: The number of runs to perform. -And one additional optional argument: +And two additional optional argument: -3. `sampling_method`: Selects the sampling method with which to select the initial guesses (`QuasiMonteCarlo.LatinHypercubeSample()` used by default). - -Because `calibrate_model_multistart` handles initial guess sampling, unlike for `calibrate_model`, no initial guess has to be provided. +2. `dirsave`: A location to which the output is automatically saved. If `dirsave=nothing`, no output is saved. It is recommended to save intermediate results for parameter estimation runs that take a long time, to not lose intermediate results if something goes wrong. +3. `sampling_method`: Selects the sampling method with which to select the initial guesses (`QuasiMonteCarlo.LatinHypercubeSample()` used by default). + +Because `calibrate_multistart` handles initial guess sampling, unlike for `calibrate`, no initial guess has to be provided. Here, we fit parameters through 10 independent optimisation runs, using QuasiMonteCarlo's `SobolSample` method, and save the result to the OptimisationRuns folder: ```@example petab1 using Optim using QuasiMonteCarlo mkdir("OptimisationRuns") # hide -res_ms = calibrate_model_multistart(petab_problem, IPNewton(), 10, "OptimisationRuns"; sampling_method=QuasiMonteCarlo.SobolSample()) -res_ms = calibrate_model_multistart(petab_problem, IPNewton(), 10, "OptimisationRuns"; sampling_method=QuasiMonteCarlo.SobolSample()) # 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 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. +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. To load the result in a later session, we can call: ```@example petab1 -res_ms = PEtabMultistartOptimisationResult("OptimisationRuns") +res_ms = PEtabMultistartResult("OptimisationRuns") nothing # hide ``` -where `"OptimisationRuns"` is the name of the save directory (specified in `calibrate_model_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: +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 = PEtabMultistartOptimisationResult("OptimisationRuns"; which_run="2") +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. - ## [Events](@id petab_events) So far, we have assumed that all experiments, after initiation, run without interference. Experiments where conditions change, or where species are added/removed during the time course, can be represented through events (related to [callbacks](@ref advanced_simulations_callbacks)). In PEtab, an event is represented through the `PEtabEvent` structure. It takes three arguments: 1. The condition for triggering the event. This can either indicate a point in time, or a boolean condition. @@ -487,27 +479,27 @@ 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; state_map=u0, events=events) +petab_model = PEtabModel(rn, observables, measurements, params; speciemap=u0, events=events) nothing # hide ``` -More details on how to use events, including how to create events with multiple targets, can be found in [PEtab.jl's documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_event/). +More details on how to use events, including how to create events with multiple targets, can be found in [PEtab.jl's documentation](https://sebapersson.github.io/PEtab.jl/stable/petab_event/). !!! note PEtab currently ignores events [created as a part of a Catalyst `ReactionSystem` model](@ref constraint_equations_events), and does not support SciML-style events [implemented through `callbacks` to `solve`](@ref advanced_simulations_callbacks). Instead, events have to use the preceding interface. ## [Plot recipes](@id petab_plotting) -There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@ref petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. +There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate` or a `calibrate_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@ref petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. 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_model(petab_problem, p0, IPNewton(); save_trace=true) # hide +res = calibrate(petab_problem, p0, IPNewton(); save_trace=true) # hide plot(res) ``` For a multi-start calibration run, the default output is instead a so-called waterfall plot: ```@example petab1 -res_ms = PEtabMultistartOptimisationResult("../assets/boehm___for_petab_tutorial") # hide +res_ms = PEtabMultistartResult("../assets/boehm___for_petab_tutorial") # hide plot(res_ms) ```