PEtab.jl API
PEtabModel
A PEtabModel
for parameter estimation/inference can be created by importing a PEtab parameter estimation problem in the standard format, or it can be directly defined in Julia. For the latter, observables that link the model to measurement data are provided by PEtabObservable
, parameters to estimate are defined by PEtabParameter
, and any potential events (callbacks) are specified as PEtabEvent
.
PEtab.PEtabObservable
— TypePEtabObservable(obs_formula, noise_formula; kwargs...)
Formulas defining the likelihood that links the model output to the measurement data.
obs_formula
describes how the model output relates to the measurement data, while noise_formula
describes the standard deviation (measurement error) and can be an equation or a numerical value. Both the observable and noise formulas can be a valid Julia equation. Variables used in these formulas must be either model species, model parameters, or parameters defined as PEtabParameter
. The formulas can also include time-point-specific noise and observable parameters; for more information, see the documentation.
Keyword Argument
transformation
: The transformation applied to the observable and its corresponding measurements. Valid options are:lin
(normal measurement noise),:log
,log2
or:log10
(log-normal measurement noise). See below for more details.
Description
For a measurement y
, an observable h = obs_formula
, and a standard deviation σ = noise_formula
, the PEtabObservable
defines the likelihood that links the model output to the measurement data: $\pi(y \mid h, \sigma)$. For transformation = :lin
, the measurement noise is assumed to be normally distributed, and the likelihood is given by:
\[\pi(y|h, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2}}\mathrm{exp}\bigg( -\frac{(y - h)^2}{2\sigma^2} \bigg)\]
As a special case, if $\sigma = 1$, this likelihood reduces to the least-squares objective function. For transformation = :log
, the measurement noise is assumed to be log-normally distributed, and the likelihood is given by:
\[\pi(y|h, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2 y^2}}\mathrm{exp}\bigg( -\frac{\big(\mathrm{log}(y) - \mathrm{log}(h)\big)^2}{2\sigma^2} \bigg)\]
For transformation = :log10
, the measurement noise is assumed to be log10-normally distributed, and the likelihood is given by:
\[\pi(y|h, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2 y^2}\mathrm{log}(10) }\mathrm{exp}\bigg( -\frac{\big(\mathrm{log}_{10}(y) - \mathrm{log}_{10}(h)\big)^2}{2\sigma^2} \bigg)\]
Lastly, for transformation = :log2
, the measurement noise is assumed to be log2-normally distributed, and the likelihood is given by:
\[\pi(y|h, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2 y^2}\mathrm{log}(2) }\mathrm{exp}\bigg( -\frac{\big(\mathrm{log}_{2}(y) - \mathrm{log}_{2}(h)\big)^2}{2\sigma^2} \bigg)\]
For numerical stabillity, PEtab.jl works with the log-likelihood in practice.
PEtab.PEtabParameter
— TypePEtabParameter(x; kwargs...)
Parameter estimation information for parameter x
.
All parameters to be estimated in a PEtabODEProblem
must be declared as a PEtabParameter
, and x
must be the name of a parameter that appears in the model, observable formula, or noise formula.
Keyword Arguments
lb::Float64 = 1e-3
: The lower parameter bound for parameter estimation.ub::Float64 = 1e3
: The upper parameter bound for parameter estimation.scale::Symbol = :log10
: The scale on which to estimate the parameter. Allowed options are:log10
(default),:log2
:log
, and:lin
. Estimating pest on thelog10
scale typically improves performance and is recommended.prior = nothing
: An optional continuous univariate parameter prior distribution from Distributions.jl.prior_on_linear_scale = true
: Whether the prior is on the linear scale (default) or on the transformed scale. For example, ifscale = :log10
andprior_on_linear_scale = false
, the prior acts on the transformed value;log10(x)
.estimate = true
: Whether the parameter should be estimated (default) or treated as a constant.value = nothing
: Value to use ifestimate = false
, and value retreived by theget_x
function. Defaults to the midpoint betweenlb
andub
.
Description
If a prior $\pi(x_i)$ is provided, the parameter estimation problem becomes a maximum a posteriori problem instead of a maximum likelihood problem. Practically, instead of minimizing the negative log-likelihood,$-\ell(x)$, the negative posterior is minimized:
\[\min_{\mathbf{x}} -\ell(\mathbf{x}) - \sum_{i} \pi(x_i)\]
For all parameters $i$ with a prior.
PEtab.PEtabEvent
— TypePEtabEvent(condition, affects, targets)
A model event triggered by condition
that sets the value of targets
to that of affects
.
For a collection of examples with corresponding plots, see the documentation.
Arguments
condition
: A Boolean expression that triggers the event when it transitions fromfalse
totrue
. For example, ift == c1
, the event is triggered when the model timet
equalsc1
. ForS > 2.0
, the event triggers when specieS
passes 2.0 from below.affects
: An equation of of model species and parameters that describes the effect of the event. It can be a single expression or a vector if there are multiple targets.targets
: Model species or parameters that the event acts on. Must match the dimension ofaffects
.
Then, given a dynamic model (as ReactionSystem
or ODESystem
), measurement data as a DataFrame
, and potential simulation conditions as a Dict
(see this tutorial), a PEtabModel
can be created:
PEtab.PEtabModel
— TypePEtabModel(sys, observables::Dict{String, PEtabObservable}, measurements::DataFrame,
parameters::Vector{PEtabParameter}; kwargs...)
From a ReactionSystem
or an ODESystem
model, observables
that link the model to measurements
and parameters
to estimate, create a PEtabModel
for parameter estimation.
For examples on how to create a PEtabModel
, see the documentation.
See also PEtabObservable
, PEtabParameter
, and PEtabEvent
.
Keyword Arguments
simulation_conditions = nothing
: An optional dictionary specifying initial specie values and/or model parameters for each simulation condition. Only required if the model has multiple simulation conditions.events = nothing
: Optional model events (callbacks) provided asPEtabEvent
. Multiple events should be provided as aVector
ofPEtabEvent
.verbose::Bool = false
: Whether to print progress while building thePEtabModel
.
PEtabModel(path_yaml; kwargs...)
Import a PEtab problem in the standard format with YAML file at path_yaml
into a PEtabModel
for parameter estimation.
For examples on how to import a PEtab problem, see the documentation.
Keyword Arguments
ifelse_to_callback::Bool = true
: Whether to rewriteifelse
(SBML piecewise) expressions to callbacks. This improves simulation runtime. It is strongly recommended to set this totrue
.verbose::Bool = false
: Whether to print progress while building thePEtabModel
.write_to_file::Bool = false
: Whether to write the generated Julia functions to files in the same directory as the PEtab problem. Useful for debugging.
PEtabODEProblem
From a PEtabModel
, a PEtabODEProblem
can:
PEtab.PEtabODEProblem
— TypePEtabODEProblem(model::PEtabModel; kwargs...)
From model
create a PEtabODEProblem
for parameter estimation/inference.
If no options (kwargs
) are provided, default settings are used for the ODESolver
, gradient method, and Hessian method. For more information on the default options, see the documentation.
See also ODESolver
, SteadyStateSolver
, and PEtabModel
.
Keyword Arguments
odesolver::ODESolver
: ODE solver options for computing the likelihood (objective function).odesolver_gradient::ODESolver
: ODE solver options for computing the gradient. Defaults toodesolver
if not provided.ss_solver::SteadyStateSolver
: Steady-state solver options for computing the likelihood. Only applicable for models with steady-state simulations.ss_solver_gradient::SteadyStateSolver
: Steady-state solver options for computing the gradient. Defaults toss_solver
if not provided.gradient_method::Symbol
: Method for computing the gradient. Available options and defaults are listed in the documentation.hessian_method::Symbol
: Method for computing the Hessian. As with the gradient, available options and defaults can be found in the documentation.FIM_method=nothing
: Method for computing the empirical Fisher Information Matrix (FIM). Accepts the same options ashessian_method
.sparse_jacobian=false
: Whether to use a sparse Jacobian when solving the ODE model. This can greatly improve performance for large models.sensealg
: Sensitivity algorithm for gradient computations. Available and default options depend ongradient_method
. See the documentation for details.chunksize=nothing
: Chunk size for ForwardDiff.jl when using forward-mode automatic differentiation for gradients and Hessians. If not provided, a default value is used. Tuningchunksize
can improve performance but is non-trivial.split_over_conditions::Bool=false
: Whether to split ForwardDiff.jl calls across simulation conditions when computing the gradient and/or Hessian. This improves performance for models with many condition-specific parameters, otherwise it increases runtime.reuse_sensitivities::Bool=false
: Whether to reuse forward sensitivities from the gradient computation for the Gauss-Newton Hessian approximation. Only applies whenhessian_method=:GaussNewton
andgradient_method=:ForwardEquations
. This can greatly improve performance when the optimization algorithm computes both the gradient and Hessian simultaneously.verbose::Bool = false
: Whether to print progress while building thePEtabODEProblem
.
Returns
A PEtabODEProblem
, where the key fields are:
nllh
: Compute the negative log-likelihood function for an input vectorx
;nllh(x)
. For this function and the ones below listed below, the inputx
can be either aVector
or aComponentArray
.grad!
: Compute the in-place gradient of the nllh;grad!(g, x)
.grad
: Compute the out-of-place gradient of the nllh;g = grad(x)
.hess!
: Compute the in-place Hessian of the nllh;hess!(H, x)
.hess
: Compute the out-of-place Hessian of the nllh;H = hess(x)
.FIM
: Compute the out-of-place empirical Fisher Information Matrix (FIM) of the nllh;F = FIM(x)
.chi2
: Compute the chi-squared test statistic for the nllh (see mathematical definition below);χ² = chi2(x)
.residuals
: Computes the residuals between the measurement data and model output (see the mathematical definition below);r = residuals(x)
.simulated_values
: Computes the corresponding model values for each measurement in the measurements table, in the same order as the measurements appear.lower_bounds
: Lower parameter bounds for parameter estimation, as specified when creating themodel
.upper_bounds
: Upper parameter bounds for parameter estimation, as specified when creating themodel
.
Description
Following the PEtab standard, the objective function to be used for parameter estimation created by the PEtabODEProblem
is a likelihood function, or, if priors are provided, a posterior function. The characteristics of the objective is defined in the PEtabModel
. In practice, for numerical stability, a PEtabODEProblem
works with the negative log-likelihood:
\[-\ell(\mathbf{x}) = - \sum_{i = 1}^N \ell_i(\mathbf{x}),\]
where $\ell_i$ is the likelihood for each measurement $i$. In addition, to accommodate numerical optimization packages, the PEtabODEProblem
provides efficient functions for computing the gradient $-\nabla \ell(\mathbf{x})$ and the second derivative Hessian matrix $-\nabla^2 \ell(\mathbf{x})$, using the specified or default gradient_method
and hessian_method
.
In addition to $-\ell$ and its derivatives, the PEtabODEProblem
provides functions for diagnostics and model selection. The χ² (chi-squared) value can be used for model selection [1], and it is computed as the sum of the χ² values for each measurement. For a measurement y
, an observable h = obs_formula
, and a standard deviation σ = noise_formula
, the χ² is computed as:
\[\chi^2 = \frac{(y - h)^2}{\sigma^2}\]
The residuals $r$ can be used to assess the measurement error model and are computed as:
\[r = \frac{(y - h)}{\sigma}\]
Lastly, the empirical Fisher Information Matrix (FIM) can be used for identifiability analysis [2]. It should ideally be computed with an exact Hessian method. The inverse of the FIM provides a lower bound on the covariance matrix. While the FIM can be useful, the profile-likelihood approach generally yields better results for identifiability analysis [2].
- Cedersund et al, The FEBS journal, pp 903-922 (2009).
- Raue et al, Bioinformatics, pp 1923-1929 (2009).
A PEtabODEProblem
has numerous configurable options. Two of the most important options are the ODESolver
and, for models with steady-state simulations, the SteadyStateSolver
:
PEtab.ODESolver
— TypeODESolver(solver, kwargs..)
ODE-solver options (solver
, tolerances
, etc.) used for solving the ODE model in a PEtabODEProblem
.
Any solver
from OrdinaryDiffEq.jl and Sundials.jl is supported. For solver recommendations and default options used when an ODESolver
is not provided when creating a PEtabODEProblem
, see the documentation.
More information on the available options and solvers can also be found in the documentation for DifferentialEquations.jl.
Keyword Arguments
abstol=1e-8
: Absolute tolerance when solving the ODE model. It is not recommended to increase above 1e-6 for gradients in order to obtain accurate gradients.reltol=1e-8
: Absolute tolerance when solving the ODE model. It is not recommended to increase above 1e-6 for gradients in order to obtain accurate gradients.dtmin=nothing
: Minimum acceptable step size when solving the ODE model.maxiters=10000
: Maximum number of iterations when solving the ODE model. Increasing above the default value can cause parameter estimation to take substantial longer time.verbose::Bool=true
: Whether or not warnings are displayed if the solver exits early.true
is recommended to detect if a suboptimal choice ofsolver
.adj_solver=solver
: Solver to use when solving the adjoint ODE. Only applicable ifgradient_method=:Adjoint
when creating thePEtabODEProblem
. Defaults tosolver
.abstol_adj=abstol
: Absolute tolerance when solving the adjoint ODE model. Only applicable ifgradient_method=:Adjoint
when creating thePEtabODEProblem
. Defaults toabstol
.reltol_adj=abstol
: Relative tolerance when solving the adjoint ODE model. Only applicable ifgradient_method=:Adjoint
when creating thePEtabODEProblem
. Defaults toreltol
.
PEtab.SteadyStateSolver
— TypeSteadyStateSolver(method::Symbol; kwargs...)
Steady-state solver options (method
, tolerances
, etc.) for computing the steady state, where the ODE right-hand side f
fulfills du = f(u, p, t) ≈ 0
.
The steady state can be computed in two ways. If method=:Simulate
, by simulating the ODE model until du = f(u, p, t) ≈ 0
using ODE solver options defined in the provied ODESolver
to the PEtabODEProblem. This approach is **strongly** recommended. If
method=:Rootfinding, by directly finding the root of the RHS
f(u, p, t) = 0` using any algorithm from NonlinearSolve.jl. The root-finding approach is far less reliable than the simulation approach (see description below).
Keyword Arguments
termination_check = :wrms
: Approach to check if the model has reached steady-state formethod = :Simulate
. Two approaches are supported::wrms
: Weighted root-mean-square. Terminate when: $\sqrt{\frac{1}{N} \sum_i^N \Big( \frac{du_i}{reltol * u_i + abstol}\Big)^2} < 1$:Newton
: Terminate if the step for Newton's methodΔu
is sufficiently small: $\sqrt{\frac{1}{N} \sum_i^N \Big(\frac{\Delta u_i}{reltol * u_i + abstol}\Big)^2} < 1$ The:Newton
approach requires that the Newton stepΔu
can be computed, which is only possible if the Jacobian of the RHS of the ODE model is invertible. If this is not the case, a pseudo-inverse is used ifpseudoinverse = true
, elsewrms
is used. The:Newton
termination should perform better than:wrms
as it accounts for the scale of the ODEs, but until benchmarks confirm this, we recommend:wrms
.
pseudoinverse = true
: Whether to use a pseudo-inverse if inverting the Jacobian fails fortermination_check = :Newton
.rootfinding_alg = nothing
: Root-finding algorithm from NonlinearSolve.jl to use ifmethod = :Rootfinding
. If left empty, the default NonlinearSolve.jl algorithm is used.abstol
: Absolute tolerance for the NonlinearSolve.jl root-finding problem or thetermination_check
formula. Defaults to100 * abstol_ode
, whereabstol_ode
is the tolerance for theODESolver
in thePEtabODEProblem
.reltol
: Relative tolerance for the NonlinearSolve.jl root-finding problem or thetermination_check
formula. Defaults to100 * reltol_ode
, wherereltol_ode
is the tolerance for theODESolver
in thePEtabODEProblem
.maxiters
: Maximum number of iterations to use ifmethod = :Rootfinding
.
Description
For an ODE model of the form:
\[\frac{\mathrm{d}u}{\mathrm{d}t} = du = f(u, p, t)\]
The steady state is defined by:
\[f(u, p, t) = 0.0\]
The steady state can be computed in two ways: either by simulating the ODE model from a set of initial values u₀
until du ≈ 0
, or by using a root-finding algorithm such as Newton's method to directly solve f(u, p, t) = 0.0
. While the root-finding approach is computationally more efficient (since it does not require solving the entire ODE), it is far less reliable and can converge to the wrong root. For example, in mass-action models with positive initial values, the feasible root should be positive in u
. This is generally fulfilled when computing the steady state via simulation (a negative root typically only occurs if the ODE solver fails). However, with root-finding, there is no such guarantee, as any root that satisfies f(u, p, t) = 0.0
may be returned. Consistent with this, benchmarks have shown that simulation is the most reliable method [1].
Another alternative is to solve for the steady state symbolically. If feasible, this is the most computationally efficient approach [1].
- Fiedler et al, BMC system biology, pp 1-19 (2016)
PEtab.jl provides several functions for interacting with a PEtabODEProblem
:
PEtab.get_x
— Functionget_x(prob::PEtabODEProblem; linear_scale = false)::ComponentArray
Get the nominal parameter vector with parameters in the correct order expected by prob
for parameter estimation/inference. Nominal values can optionally be specified when creating a PEtabParameter
, or in the parameters table if the problem is provided in the PEtab standard format.
For ease of interaction (e.g., changing values), the parameter vector is returned as a ComponentArray
. For how to interact with a ComponentArray
, see the documentation and the ComponentArrays.jl documentation.
See also PEtabParameter
.
Keyword argument
linear_scale
: Whether to return parameters on the linear scale. By default, parameters are returned on the scale they are estimated, which by default islog10
as this often improves parameter estimation performance.
SciMLBase.remake
— Methodremake(prob::PEtabODEProblem, xchange::Dict)::PEtabODEProblem
Fix and consequently remove a set of parameters from being estimated without rebuilding prob
.
xchange
should be provided as a Dict
with parameter names and their new values. For example, to fix k1
to 1.0
, provide xchange = Dict(:k1 => 1.0)
. The parameters to be fixed can only be those that were originally set to be estimated.
remake
is the most efficient approach for fixing parameters, as it does not re-compile the problem.
And additionally, functions for interacting with the underlying dynamic model (ODEProblem
) within a PEtabODEProblem
:
PEtab.get_u0
— Functionget_u0(res, prob::PEtabODEProblem; kwargs...)
Retrieve the ODEProblem
initial values for simulating the ODE model in prob
. res
can be a parameter estimation result (e.g., PEtabMultistartResult
) or a Vector
with parameters in the order expected by prob
(see get_x
).
For information on keyword arguements see get_ps
.
See also get_odeproblem
and get_odesol
.
PEtab.get_ps
— Functionget_ps(res, prob::PEtabODEProblem; kwargs...)
Retrieve the ODEProblem
parameter values for simulating the ODE model in prob
. res
can be a parameter estimation result (e.g., PEtabMultistartResult
) or a Vector
with parameters in the order expected by prob
(see get_x
).
See also: get_u0
, get_odeproblem
, get_odesol
.
Keyword Arguments
retmap=true
: Whether to return the values as a map in the form[k1 => val1, ...]
. Such a map can be directly used when building anODEProblem
. Iffalse
, aVector
is returned. This keyword is only applicable forget_u0
andget_ps
.cid::Symbol
: Which simulation condition to return parameters for. If not provided, defaults to the first simulation condition. For otherget
functions, theODEProblem
,u0
, orODESolution
for the specifiedcid
is returned.preeq_id
: Which potential pre-equilibration (steady-state) simulation id to use. If a validpreeq_id
is provided, the ODE is first simulated to steady state forpreeq_id
. Then the model shifts tocid
, and the parameters forcid
are returned. For otherget
functions, theODEProblem
,u0
, orODESolution
for the specifiedcid
is returned
PEtab.get_odeproblem
— Functionget_odeproblem(res, prob::PEtabODEProblem; kwargs...) -> (sol, callbacks)
Retrieve the ODEProblem
and callbacks (CallbackSet
) for simulating the ODE model in prob
. res
can be a parameter estimation result (e.g., PEtabMultistartResult
) or a Vector
with parameters in the order expected by prob
(see get_x
).
For information on keyword arguements see get_ps
.
See also: get_u0
and get_odesol
.
PEtab.get_odesol
— Functionget_odesol(res, prob::PEtabODEProblem; kwargs...)::ODESolution
Retrieve the ODESolution
from simulating the ODE model in prob
. res
can be a parameter estimation result (e.g., PEtabMultistartResult
) or a Vector
with parameters in the order expected by prob
(see get_x
).
For information on keyword arguements see get_ps
.
See also: get_u0
and get_odeproblem
.
PEtab.solve_all_conditions
— Functionsolve_all_conditions(x, prob::PEtabODEProblem, solver; kwargs)
Solve the ODE model in prob
for all simulation conditions with the provided ODE-solver.
x
should be a Vector
or ComponentArray
with parameters in the order expected by prob
(see get_x
).
Keyword Arguments
abstol=1e-8
: Absolute tolerance for the ODE solver.reltol=1e-8
: Relative tolerance for the ODE solver.maxiters=1e4
: Maximum iterations for the ODE solver.ntimepoints_save=0
: The number of time points at which to save the ODE solution for each condition. A value of 0 means the solution is saved at the solvers default time points.save_observed_t=false
: When set to true, this option overridesntimepoints_save
and saves the ODE solution only at the time points where measurement data is available.
Returns
odesols
: A dictionary containing theODESolution
for each simulation condition.
Parameter Estimation
A PEtabODEProblem
contains all the necessary information for wrapping a suitable numerical optimization library, but for convenience, PEtab.jl provides wrappers for several available optimizers. In particular, single-start parameter estimation is provided via calibrate
:
PEtab.calibrate
— Functioncalibrate(prob::PEtabODEProblem, x0, alg; kwargs...)::PEtabOptimisationResult
From starting point x0
using optimization algorithm alg
, estimate unknown model parameters for prob
, and get results as a PEtabOptimisationResult
.
x0
can be a Vector
or a ComponentArray
, where the individual parameters must be in the order expected by prob
. To get a vector in the correct order, see get_x
.
A list of available and recommended optimization algorithms (alg
) can be found in the documentation. Briefly, supported algorithms are from:
- Optim.jl:
LBFGS()
,BFGS()
, orIPNewton()
methods. - Ipopt.jl:
IpoptOptimizer()
interior-point Newton method. - Fides.py:
Fides()
Newton trust region method.
Different ways to visualize the parameter estimation result can be found in the documentation.
See also PEtabOptimisationResult
, Fides
and IpoptOptimizer
Keyword Arguments
save_trace::Bool = false
: Whether to save the optimization trace of the objective function and parameter vector. Only applicable for some algorithms; see the documentation for details.options = DEFAULT_OPTIONS
: Configurable options foralg
. The type and available options depend on which packagealg
belongs to. For example, ifalg = IPNewton()
from Optim.jl,options
should be provided as anOptim.Options()
struct. A list of configurable options can be found in the documentation.
PEtab.PEtabOptimisationResult
— TypePEtabOptimisationResult
Parameter estimation statistics from single-start optimization with calibrate
.
See also: calibrate
Fields
xmin
: Minimizing parameter vector found by the optimization.fmin
: Minimum objective value found by the optimization.x0
: Starting parameter vector.alg
: Parameter estimation algorithm used.niterations
: Number of iterations for the optimization.runtime
: Runtime in seconds for the optimization.xtrace
: Parameter vector optimization trace. Empty ifsave_trace = false
was provided tocalibrate
.ftrace
: Objective function optimization trace. Empty ifsave_trace = false
was provided tocalibrate
.converged
: Convergence flag fromalg
.original
: Original result struct returned byalg
. For example, ifalg = IPNewton()
from Optim.jl,original
is the Optim return struct.
Multi-start (recommended method) parameter estimation, is provided via calibrate_multistart
:
PEtab.calibrate_multistart
— Functioncalibrate_multistart(prob::PEtabODEProblem, alg, nmultistarts::Integer; nprocs = 1,
dirsave=nothing, kwargs...)::PEtabMultistartResult
Perform nmultistarts
parameter estimation runs from randomly sampled starting points using the optimization algorithm alg
to estimate the unknown model parameters in prob
.
A list of available and recommended optimisation algorithms (alg
) can be found in the package documentation and in the calibrate
documentation.
If nprocs > 1
, the parameter estimation runs are performed in parallel using the pmap
function from Distributed.jl with nprocs
processes. If parameter estimation on a single process (nprocs = 1
) takes longer than 5 minutes, we strongly recommend setting nprocs > 1
, as this can greatly reduce runtime. Note that nprocs
should not be larger than the number of cores on the computer.
If dirsave
is provided, intermediate results for each run are saved in dirsave
. It is strongly recommended to provide dirsave
for larger models, as parameter estimation can take hours (or even days!),and without dirsave
, all intermediate results will be lost if something goes wrong.
Different ways to visualize the parameter estimation result can be found in the documentation.
See also PEtabMultistartResult
, get_startguesses
, and calibrate
.
Keyword Arguments
sampling_method = LatinHypercubeSample()
: Method for sampling a diverse (spread out) set of starting points. See the documentation forget_startguesses
, which is the function used for sampling.sample_prior::Bool = true
: See the documentation forget_startguesses
.seed = nothing
: Seed used when generating starting points.options = DEFAULT_OPTIONS
: Configurable options foralg
. See the documentation forcalibrate
.
PEtab.get_startguesses
— Functionget_startguesses(prob::PEtabODEProblem, n::Integer; kwargs...)
Generate n
random parameter vectors within the parameter bounds in prob
.
If n = 1
, a single random vector is returned. For n > 1
, a vector of random parameter vectors is returned. In both cases, parameter vectors are returned as a ComponentArray
. For details on how to interact with a ComponentArray
, see the documentation and the ComponentArrays.jl documentation.
See also calibrate
and calibrate_multistart
.
Keyword Arguments
sampling_method = LatinHypercubeSample()
: Method for sampling a diverse (spread out) set of parameter vectors. Any algorithm from QuasiMonteCarlo is allowed, but the defaultLatinHypercubeSample
is recommended as it usually performs well.sample_prior::Bool = true
: Whether to sample random parameter values from the prior distribution if a parameter has a prior.allow_inf::Bool = false
: Whether to return parameter vectors for which the likelihood cannot be computed (typically happens because theODEProblem
cannot be solved). Often it only makes sense to use starting points with a computable likelihood for parameter estimation, hence it typically does not make sense to change this option.
PEtab.PEtabMultistartResult
— TypePEtabMultistartResult
Parameter estimation statistics from multi-start optimization with calibrate_multistart
.
See also calibrate_multistart
and PEtabOptimisationResult
.
Fields
xmin
: Best minimizer across all runs.fmin
: Best minimum across all runs.alg
: Parameter estimation algorithm.nmultistarts
: Number of parameter estimation runs.sampling_method
: Sampling method used for generating starting points.dirsave
: Path of directory where parameter estimation run statistics are saved ifdirsave
was provided tocalibrate_multistart
.runs
: Vector ofPEtabOptimisationResult
with the parameter estimation results for each run.
PEtabMultistartResult(dirres::String; which_run::String="1")
Import multistart parameter estimation results saved at dirres
.
Each time a new optimization run is performed, results are saved with unique numerical endings. Results from a specific run can be retrieved by specifying the numerical ending with which_run
.
Lastly, model selection is provided via petab_select
:
PEtab.petab_select
— Functionpetab_select(path_yaml, alg; nmultistarts = 100, kwargs...) -> path_res
For a PEtab-select problem, perform model selection with the method specified in the PEtab select problem files. Returns the path (path_res
) to a YAML-file with model selection results.
The general model selection (e.g. to use forward-search) options are specified in the PEtab-select files. For details on how to set this up, see the PEtab-select documentation.
For each round of model selection, the candidate models are parameter estimated using multi-start parameter estimation, with nmultistarts
performed for each model. The objective values obtained from parameter estimation are then used for the next round of model evaluation.
A list of available and recommended optimization algorithms (alg
) can be found in the package documentation and calibrate
documentation.
See also calibrate_multistart
.
Keyword Arguments
kwargs
: The same keywords accepted byPEtabODEProblem
andcalibrate
.
For each case case, PEtab.jl supports the usage of optimization algorithms from Optim.jl, Ipopt.jl, and Fides.py:
PEtab.Fides
— TypeFides(hessian_method; verbose::Bool=false)
Setup the Fides box-constrained Newton-trust region optimizer for parameter estimation.
See also calibrate
and calibrate_multistart
.
Arguments
hessian_method
: Method for computing the Hessian. Allowed options are:nothing
: The Hessian computed by thePEtabODEProblem
is used.:BB
: Broyden's "bad" method.:BFGS
: Broyden-Fletcher-Goldfarb-Shanno update strategy.:BG
: Broyden's "good" method.:Broyden
: Broyden-class update scheme.:SR1
: Symmetric Rank 1 update.:SSM
: Structured Secant Method.:TSSM
: Totally Structured Secant Method.
verbose
: Whether to print progress during the parameter estimation.
Description
Fides is a Newton-trust region optimizer for box-constrained optmization problems. More information on the algorithm can be found in [1]. Fides particularly excels when the full Hessian is too computationally expensive to compute, but a Gauss-Newton Hessian approximation can be computed (for more details see the documentation). In addition to supporting user Hessians via the PEtabODEProblem
, it supports several Hessian approximation methods. Aa more extensive description than above see the Fides documentation.
- Fröhlich and Sorger, PLoS computational biology, pp e1010322 (2022)
PEtab.IpoptOptimizer
— TypeIpoptOptimizer(LBFGS::Bool)
Setup the Ipopt Interior-point Newton method optmizer for parameter estimation.
Ipopt can be configured to use either the Hessian method from the PEtabODEProblem
(LBFGS=false
) or an LBFGS scheme (LBFGS=true
). For setting other Ipopt options, see IpoptOptions
.
See also calibrate
and calibrate_multistart
.
Description
Ipopt is an Interior-point Newton method for constrained non-linear optimization problems. More information on the algorithm can be found in [1].
- Wächter and Biegler, Mathematical programming, pp 25-57 (2006)
PEtab.IpoptOptions
— TypeIpoptOptions(; kwargs...)
Options for parameter estimation with IpoptOptimizer
.
More details on the options can be found in the Ipopt documentation.
See also IpoptOptimizer
, calibrate
, and calibrate_multistart
.
Keyword Arguments
print_level = 0
: Output verbosity level (valid values are 0 ≤ print_level ≤ 12)max_iter = 1000
: Maximum number of iterationstol = 1e-8
: Relative convergence toleranceacceptable_tol = 1e-6
: Acceptable relative convergence tolerancemax_wall_time 1e20
: Maximum wall time optimization is allowed to runacceptable_obj_change_tol 1e20
: Acceptance stopping criterion based on objective function change.
Parameter estimation results can be visualized using the plot-recipes detailed in this page, and with get_obs_comparison_plots
:
PEtab.get_obs_comparison_plots
— Functionget_obs_comparison_plots(res, prob::PEtabODEProblem; kwargs...)::Dict
Plot the model fit against data for all simulation conditions and all observable ids for a PEtab.jl parameter estimation result (res
).
Each entry in the returned Dict
corresponds to plot(res, prob; obsids=[obsid], cid=cid, kwargs...)
for all possible cid
and all obsid
.
As an alternative to the PEtab.jl interface to parameter estimation, a PEtabODEProblem
can be converted to an OptimizationProblem
to access the algorithms available via Optimization.jl:
PEtab.OptimizationProblem
— FunctionOptimizationProblem(prob::PEtabODEProblem; box_constraints::Bool = true)
Create an Optimization.jl OptimizationProblem
from prob
.
To use algorithms not compatible with box constraints (e.g., Optim.jl NewtonTrustRegion
), set box_constraints = false
. Note that with this option, optimizers may move outside the bounds, which can negatively impact performance. More information on how to use an OptimizationProblem
can be found in the Optimization.jl documentation.
Bayesian Inference
PEtab.jl offers wrappers to perform Bayesian inference using state-of-the-art methods such as NUTS (the same sampler used in Turing.jl) or AdaptiveMCMC.jl. It should be noted that this part of PEtab.jl is planned to be moved to a separate package, so the syntax will change and be made more user-friendly in the future.
PEtab.PEtabLogDensity
— TypePEtabLogDensity(prob::PEtabODEProblem)
Create a LogDensityProblem
using the posterior and gradient functions from prob
.
This LogDensityProblem
interface defines everything needed to perform Bayesian inference with packages such as AdvancedHMC.jl
(which includes algorithms like NUTS, used by Turing.jl
), and AdaptiveMCMC.jl
.
PEtab.to_prior_scale
— Functionto_prior_scale(xpetab, target::PEtabLogDensity)
Transforms parameter x
from the PEtab problem scale to the prior scale.
This conversion is needed for Bayesian inference, as in PEtab.jl Bayesian inference is performed on the prior scale.
To use this function, the Bijectors, LogDensityProblems, and LogDensityProblemsAD packages must be loaded: using Bijectors, LogDensityProblems, LogDensityProblemsAD
PEtab.to_chains
— Functionto_chains(res, target::PEtabLogDensity; kwargs...)::MCMCChains
Converts Bayesian inference results obtained with PEtabLogDensity
into an MCMCChains
.
res
can be the inference results from AdvancedHMC.jl or AdaptiveMCMC.jl. The returned chain has the parameters on the prior scale.
Keyword Arguments
start_time
: Optional starting time for the inference, obtained withnow()
.end_time
: Optional ending time for the inference, obtained withnow()
.
To use this function, the MCMCChains package must be loaded: using MCMCChains