Getting started
In this introductory tutorial we will create a PEtabODEproblem
for the small Boehm model, while simultaneously covering the main features in PEtab.jl.
To run the code you need the Boehm PEtab files which can be found here. A fully runnable example of this tutorial can be found here.
Reading a PEtab model
As a starting point we need to read the PEtab files into Julia. This is made easy by the package which given the path to the PEtab yaml-file PEtab.jl reads all PEtab files into a PEtabModel
struct. Here several things happen under the hood:
- The SBML file is translated into ModelingToolkit.jl format (e.g. allow symbolic computation of the ODE-model Jacobian).
- The observable PEtab table is translated into Julia functions for computing the observables ($h$), noise parameter ($\sigma$) and initial values ($u_0$).
- To be able to compute gradients via adjoint sensitivity analysis and/or forward sensitivity equations the derivatives of $h$ and $\sigma$ are computed symbolically with respect to the ODE-models states ($u$) and parameters.
All this happens automatically and you can find the resulting files in dirYamlFile/Julia_files/. To save time the function readPEtabModel
has the default forceBuildJlFiles=false
meaning that the Julia files are not rebuilt in case they already exist.
using PEtab
using OrdinaryDiffEq
using Printf
pathYaml = joinpath(@__DIR__, "Boehm", "Boehm_JProteomeRes2014.yaml")
petabModel = readPEtabModel(pathYaml, verbose=true)
PEtabModel for model Boehm. ODE-system has 8 states and 10 parameters.
Generated Julia files are at ...
Creating a PEtabODEProblem
Given a PEtab model we can create a PEtabODEProblem
(in the future we plan to add surrogate, SDE, etc... problems). The function setupPEtabODEProblem
accepts several (full list in API documentation), some key ones are:
odeSolverOptions
- Which ODE solver and which solver tolerances (abstol
andreltol
). Below we use the ODE solverRodas5P()
(works well for smaller models ≤ 15 states), and we use the defaultabstol, reltol .= 1e-8
.gradientMethod
- For small models like Boehm forward mode automatic differentiation (AD) is fastest, thus we choose:ForwardDiff
.hessianMethod
- For small models like Boehm with ≤20 parameters it is computationally feasible to compute the full Hessian via forward-mode AD. Thus, we choose:ForwardDiff
.
odeSolverOptions = getODESolverOptions(Rodas5P(), abstol=1e-8, reltol=1e-8)
petabProblem = setupPEtabODEProblem(petabModel, odeSolverOptions,
gradientMethod=:ForwardDiff,
hessianMethod=:ForwardDiff)
PEtabODEProblem for Boehm. ODE-states: 8. Parameters to estimate: 9 where 6 are dynamic.
---------- Problem settings ----------
Gradient method : ForwardDiff
Hessian method : ForwardDiff
--------- ODE-solver settings --------
Cost Rodas5P(). Options (abstol, reltol, maxiters) = (1.0e-08, 1.0e-08, 1.0e+04)
Gradient Rodas5P(). Options (abstol, reltol, maxiters) = (1.0e-08, 1.0e-08, 1.0e+04)
Computing the cost, gradient and hessian
The PEtabODEProblem
contains everything needed to set up an optimization problem with most available optimizers. Main fields are:
petabODEProblem.computeCost
- Given a parameter vector θ computes the cost (objective function).petabODEProblem.computeGradient!
- Given a parameter vector θ computes the gradient via chosen method.petabODEProblem.computeHessian!
- Given a parameter vector θ computes the hessian via chosen methodpetabODEProblem.lowerBounds
- A vector with the lower bounds for parameters as specified by the PEtab parameters file.petabODEProblem.upperBounds
- A vector with the upper bounds for parameters as specified by the PEtab parameters file.petabODEProblem.θ_estNames
- A vector with the names of the parameters to estimate.
- Note1 - The parameter vector θ is assumed to be on PEtab specified parameter-scale. Thus, if parameter $i$ is on the log-scale so should
θ[i]
be. - Note2 - The
computeGradient!
andcomputeHessian!
functions are in-place functions. Thus, their first argument is an already pre-allocated gradient and hessian respectively (see below).
# Parameters on log-scale
p = petabProblem.θ_nominalT
gradient = zeros(length(p))
hessian = zeros(length(p), length(p))
cost = petabProblem.computeCost(p)
petabProblem.computeGradient!(gradient, p)
petabProblem.computeHessian!(hessian, p)
@printf("Cost = %.2f\n", cost)
@printf("First element in the gradient = %.2e\n", gradient[1])
@printf("First element in the hessian = %.2f\n", hessian[1, 1])
Cost = 138.22
First element in the gradient = 2.20e-02
First element in the hessian = 2199.49
Where to go from here
The best create a PEtabODEProblem
we recommend looking at Choosing best options for a PEtab problem. We also recommend looking at Supported gradient and hessian methods.