pyPESTO vs no pyPESTO
The objectives of this notebook are twofold:
General Workflow: We walk step-by-step through a process to estimate parameters of dynamical models. By following this workflow, you will gain a clear understanding of the essential steps involved and how they contribute to the overall outcome.
Benefits of pyPESTO: Throughout the notebook, we highlight the key advantages of using pyPESTO in each step of the workflow, compared to “doing things manually”. By leveraging its capabilities, you can significantly increase efficiency and effectiveness when solving your parameter optimization tasks.
This notebook is divided into several sections, each focusing on a specific aspect of the parameter estimation workflow. Here’s an overview of what you find in each section:
Contents
Objective Function: We discuss the creation of an objective function that quantifies the goodness-of-fit between a model and observed data. We will demonstrate how pyPESTO simplifies this potentially cumbersome process and provides various options for objective function definition.
Optimization: We show how to find optimal model parameters. We illustrate the general workflow and how pyPESTO allows to flexibly use different optimizers and to analyze and interpret results.
Profiling: After a successful parameter optimization, we show how pyPESTO provides profile likelihood functionality to assess uncertainty and identifiability of selected parameters.
Sampling: In addition to profiles, we use MCMC to sample from the Bayesian posterior distribution. We show how pyPESTO facilitates the use of different sampling methods.
Result Storage: This section focuses on storing and organizing the results obtained from the parameter optimization workflow, which is necessary to keep results for later processing and sharing.
By the end of this notebook, you’ll have gained valuable insights into the parameter estimation workflow for dynamical models. Moreover, you’ll have a clear understanding of the benefits pyPESTO brings to each step of this workflow. This tutorial will equip you with the knowledge and tools necessary to streamline your parameter estimation tasks and obtain accurate and reliable results.
[1]:
# install dependencies
#!pip install pypesto[amici,petab]
[1]:
# imports
import logging
import os
import random
from pprint import pprint
import amici
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import petab
import scipy.optimize
from IPython.display import Markdown, display
import pypesto.optimize as optimize
import pypesto.petab
import pypesto.profile as profile
import pypesto.sample as sample
import pypesto.store as store
import pypesto.visualize as visualize
import pypesto.visualize.model_fit as model_fit
mpl.rcParams['figure.dpi'] = 100
mpl.rcParams['font.size'] = 18
# for reproducibility
random.seed(1912)
np.random.seed(1912)
# name of the model
model_name = "boehm_JProteomeRes2014"
# output directory
model_output_dir = "tmp/" + model_name
1. Create an objective function
As application problem, we consider the model by Böhm et al., JProteomRes 2014, which describes, trained on quantitative mass spectronomy data, the process of dimerization of phosphorylated STAT5A and STAT5B, important transductors of activation signals of cytokine receptors to the nucleus. The model is available via the PEtab benchmark collection. For simulation, we use AMICI, an efficient ODE simulation and sensitivity calculation routine. PEtab is a data format specification that standardises parameter estimation problems in systems biology..
Without pyPESTO
To fit an (ODE) model to data, the model needs to be implemented in a simulation program as a function mapping parameters to simulated data. Simulations must then be mapped to experimentally observed data via formulation of a single-value cost function (e.g. squared or absolute differences, corresponding to a normal or Laplace measurement noise model). Loading the model via PEtab and AMICI already simplifies these stepss substantially compared to encoding the model and the objective function manually:
[3]:
%%capture
# PEtab problem loading
petab_yaml = f"./{model_name}/{model_name}.yaml"
petab_problem = petab.Problem.from_yaml(petab_yaml)
# AMICI model complilation
amici_model = amici.petab_import.import_petab_problem(
petab_problem, force_compile=True
)
2023-08-25 15:26:01.472 - amici.petab_import - INFO - Importing model ...
2023-08-25 15:26:01.473 - amici.petab_import - INFO - Validating PEtab problem ...
2023-08-25 15:26:01.558 - amici.petab_import - INFO - Model name is 'FullModel'.
Writing model code to '/Users/pauljonasjost/Documents/GitHub_Folders/pyPESTO/doc/example/amici_models/FullModel'.
2023-08-25 15:26:01.560 - amici.petab_import - INFO - Species: 8
2023-08-25 15:26:01.562 - amici.petab_import - INFO - Global parameters: 15
2023-08-25 15:26:01.563 - amici.petab_import - INFO - Reactions: 9
2023-08-25 15:26:01.618 - amici.petab_import - INFO - Observables: 3
2023-08-25 15:26:01.620 - amici.petab_import - INFO - Sigmas: 3
2023-08-25 15:26:01.632 - amici.petab_import - DEBUG - Adding output parameters to model: ['noiseParameter1_pSTAT5A_rel', 'noiseParameter1_pSTAT5B_rel', 'noiseParameter1_rSTAT5A_rel']
2023-08-25 15:26:01.634 - amici.petab_import - DEBUG - Adding initial assignments for dict_keys([])
2023-08-25 15:26:01.657 - amici.petab_import - DEBUG - Condition table: (1, 1)
2023-08-25 15:26:01.658 - amici.petab_import - DEBUG - Fixed parameters are ['ratio', 'specC17']
2023-08-25 15:26:01.660 - amici.petab_import - INFO - Overall fixed parameters: 2
2023-08-25 15:26:01.661 - amici.petab_import - INFO - Variable parameters: 16
2023-08-25 15:26:01.684 - amici.sbml_import - DEBUG - Finished processing SBML annotations ++ (1.34E-04s)
2023-08-25 15:26:01.717 - amici.sbml_import - DEBUG - Finished gathering local SBML symbols ++ (2.14E-02s)
2023-08-25 15:26:01.752 - amici.sbml_import - DEBUG - Finished processing SBML parameters ++ (2.24E-02s)
2023-08-25 15:26:01.765 - amici.sbml_import - DEBUG - Finished processing SBML compartments ++ (4.04E-04s)
2023-08-25 15:26:01.788 - amici.sbml_import - DEBUG - Finished processing SBML species initials +++ (6.57E-03s)
2023-08-25 15:26:01.799 - amici.sbml_import - DEBUG - Finished processing SBML rate rules +++ (7.80E-05s)
2023-08-25 15:26:01.801 - amici.sbml_import - DEBUG - Finished processing SBML species ++ (2.70E-02s)
2023-08-25 15:26:01.818 - amici.sbml_import - DEBUG - Finished processing SBML reactions ++ (5.84E-03s)
2023-08-25 15:26:01.838 - amici.sbml_import - DEBUG - Finished processing SBML rules ++ (9.97E-03s)
2023-08-25 15:26:01.849 - amici.sbml_import - DEBUG - Finished processing SBML events ++ (8.28E-05s)
2023-08-25 15:26:01.859 - amici.sbml_import - DEBUG - Finished processing SBML initial assignments++ (1.03E-04s)
2023-08-25 15:26:01.870 - amici.sbml_import - DEBUG - Finished processing SBML species references ++ (5.32E-04s)
2023-08-25 15:26:01.871 - amici.sbml_import - DEBUG - Finished importing SBML + (1.95E-01s)
2023-08-25 15:26:01.926 - amici.sbml_import - DEBUG - Finished processing SBML observables + (4.47E-02s)
2023-08-25 15:26:01.938 - amici.sbml_import - DEBUG - Finished processing SBML event observables + (2.83E-06s)
2023-08-25 15:26:02.017 - amici.de_export - DEBUG - Finished running smart_multiply ++ (2.33E-03s)
2023-08-25 15:26:02.122 - amici.de_export - DEBUG - Finished simplifying xdot +++ (8.01E-03s)
2023-08-25 15:26:02.123 - amici.de_export - DEBUG - Finished computing xdot ++ (1.85E-02s)
2023-08-25 15:26:02.147 - amici.de_export - DEBUG - Finished simplifying x0 +++ (1.72E-03s)
2023-08-25 15:26:02.149 - amici.de_export - DEBUG - Finished computing x0 ++ (1.35E-02s)
2023-08-25 15:26:02.152 - amici.de_export - DEBUG - Finished importing SbmlImporter + (1.45E-01s)
2023-08-25 15:26:02.334 - amici.de_export - DEBUG - Finished simplifying Jy ++++ (1.42E-01s)
2023-08-25 15:26:02.335 - amici.de_export - DEBUG - Finished computing Jy +++ (1.51E-01s)
2023-08-25 15:26:02.404 - amici.de_export - DEBUG - Finished simplifying y ++++ (4.67E-02s)
2023-08-25 15:26:02.405 - amici.de_export - DEBUG - Finished computing y +++ (5.86E-02s)
2023-08-25 15:26:02.425 - amici.de_export - DEBUG - Finished simplifying sigmay ++++ (1.38E-04s)
2023-08-25 15:26:02.426 - amici.de_export - DEBUG - Finished computing sigmay +++ (8.85E-03s)
2023-08-25 15:26:02.458 - amici.de_export - DEBUG - Finished writing Jy.cpp ++ (2.84E-01s)
2023-08-25 15:26:02.523 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (3.62E-02s)
2023-08-25 15:26:02.547 - amici.de_export - DEBUG - Finished simplifying dJydsigma ++++ (1.39E-02s)
2023-08-25 15:26:02.548 - amici.de_export - DEBUG - Finished computing dJydsigma +++ (7.16E-02s)
2023-08-25 15:26:02.558 - amici.de_export - DEBUG - Finished writing dJydsigma.cpp ++ (8.93E-02s)
2023-08-25 15:26:02.601 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (1.81E-02s)
2023-08-25 15:26:02.629 - amici.de_export - DEBUG - Finished simplifying dJydy ++++ (1.72E-02s)
2023-08-25 15:26:02.630 - amici.de_export - DEBUG - Finished computing dJydy +++ (5.57E-02s)
2023-08-25 15:26:02.641 - amici.de_export - DEBUG - Finished writing dJydy.cpp ++ (7.41E-02s)
2023-08-25 15:26:02.669 - amici.de_export - DEBUG - Finished simplifying Jz ++++ (9.60E-05s)
2023-08-25 15:26:02.670 - amici.de_export - DEBUG - Finished computing Jz +++ (8.37E-03s)
2023-08-25 15:26:02.682 - amici.de_export - DEBUG - Finished computing z +++ (1.74E-04s)
2023-08-25 15:26:02.701 - amici.de_export - DEBUG - Finished simplifying sigmaz ++++ (1.33E-04s)
2023-08-25 15:26:02.701 - amici.de_export - DEBUG - Finished computing sigmaz +++ (7.93E-03s)
2023-08-25 15:26:02.702 - amici.de_export - DEBUG - Finished writing Jz.cpp ++ (4.87E-02s)
2023-08-25 15:26:02.729 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (7.59E-05s)
2023-08-25 15:26:02.739 - amici.de_export - DEBUG - Finished simplifying dJzdsigma ++++ (2.09E-04s)
2023-08-25 15:26:02.740 - amici.de_export - DEBUG - Finished computing dJzdsigma +++ (1.98E-02s)
2023-08-25 15:26:02.742 - amici.de_export - DEBUG - Finished writing dJzdsigma.cpp ++ (2.80E-02s)
2023-08-25 15:26:02.769 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (6.99E-05s)
2023-08-25 15:26:02.778 - amici.de_export - DEBUG - Finished simplifying dJzdz ++++ (9.28E-05s)
2023-08-25 15:26:02.779 - amici.de_export - DEBUG - Finished computing dJzdz +++ (1.78E-02s)
2023-08-25 15:26:02.780 - amici.de_export - DEBUG - Finished writing dJzdz.cpp ++ (2.73E-02s)
2023-08-25 15:26:02.807 - amici.de_export - DEBUG - Finished simplifying Jrz ++++ (1.01E-04s)
2023-08-25 15:26:02.808 - amici.de_export - DEBUG - Finished computing Jrz +++ (7.65E-03s)
2023-08-25 15:26:02.818 - amici.de_export - DEBUG - Finished computing rz +++ (2.55E-04s)
2023-08-25 15:26:02.819 - amici.de_export - DEBUG - Finished writing Jrz.cpp ++ (2.59E-02s)
2023-08-25 15:26:02.845 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (9.20E-05s)
2023-08-25 15:26:02.856 - amici.de_export - DEBUG - Finished simplifying dJrzdsigma ++++ (9.63E-05s)
2023-08-25 15:26:02.857 - amici.de_export - DEBUG - Finished computing dJrzdsigma +++ (2.05E-02s)
2023-08-25 15:26:02.858 - amici.de_export - DEBUG - Finished writing dJrzdsigma.cpp ++ (2.91E-02s)
2023-08-25 15:26:02.886 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (6.81E-05s)
2023-08-25 15:26:02.899 - amici.de_export - DEBUG - Finished simplifying dJrzdz ++++ (1.23E-04s)
2023-08-25 15:26:02.900 - amici.de_export - DEBUG - Finished computing dJrzdz +++ (2.28E-02s)
2023-08-25 15:26:02.901 - amici.de_export - DEBUG - Finished writing dJrzdz.cpp ++ (3.12E-02s)
2023-08-25 15:26:02.928 - amici.de_export - DEBUG - Finished simplifying root ++++ (1.70E-04s)
2023-08-25 15:26:02.929 - amici.de_export - DEBUG - Finished computing root +++ (9.40E-03s)
2023-08-25 15:26:02.931 - amici.de_export - DEBUG - Finished writing root.cpp ++ (1.89E-02s)
2023-08-25 15:26:02.999 - amici.de_export - DEBUG - Finished simplifying w +++++ (3.15E-02s)
2023-08-25 15:26:03.000 - amici.de_export - DEBUG - Finished computing w ++++ (4.06E-02s)
2023-08-25 15:26:03.035 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (2.51E-02s)
2023-08-25 15:26:03.055 - amici.de_export - DEBUG - Finished simplifying dwdp ++++ (1.14E-02s)
2023-08-25 15:26:03.056 - amici.de_export - DEBUG - Finished computing dwdp +++ (1.05E-01s)
2023-08-25 15:26:03.078 - amici.de_export - DEBUG - Finished simplifying spl ++++ (1.36E-04s)
2023-08-25 15:26:03.079 - amici.de_export - DEBUG - Finished computing spl +++ (8.40E-03s)
2023-08-25 15:26:03.101 - amici.de_export - DEBUG - Finished simplifying sspl ++++ (1.03E-04s)
2023-08-25 15:26:03.103 - amici.de_export - DEBUG - Finished computing sspl +++ (1.16E-02s)
2023-08-25 15:26:03.110 - amici.de_export - DEBUG - Finished writing dwdp.cpp ++ (1.66E-01s)
2023-08-25 15:26:03.219 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (8.22E-02s)
2023-08-25 15:26:03.318 - amici.de_export - DEBUG - Finished simplifying dwdx ++++ (8.91E-02s)
2023-08-25 15:26:03.319 - amici.de_export - DEBUG - Finished computing dwdx +++ (1.91E-01s)
2023-08-25 15:26:03.359 - amici.de_export - DEBUG - Finished writing dwdx.cpp ++ (2.39E-01s)
2023-08-25 15:26:03.369 - amici.de_export - DEBUG - Finished writing create_splines.cpp ++ (3.98E-04s)
2023-08-25 15:26:03.404 - amici.de_export - DEBUG - Finished simplifying spline_values +++++ (1.11E-04s)
2023-08-25 15:26:03.405 - amici.de_export - DEBUG - Finished computing spline_values ++++ (9.25E-03s)
2023-08-25 15:26:03.417 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (1.12E-04s)
2023-08-25 15:26:03.427 - amici.de_export - DEBUG - Finished simplifying dspline_valuesdp ++++ (9.45E-05s)
2023-08-25 15:26:03.428 - amici.de_export - DEBUG - Finished computing dspline_valuesdp +++ (3.99E-02s)
2023-08-25 15:26:03.429 - amici.de_export - DEBUG - Finished writing dspline_valuesdp.cpp ++ (4.87E-02s)
2023-08-25 15:26:03.464 - amici.de_export - DEBUG - Finished simplifying spline_slopes +++++ (1.12E-04s)
2023-08-25 15:26:03.465 - amici.de_export - DEBUG - Finished computing spline_slopes ++++ (9.20E-03s)
2023-08-25 15:26:03.476 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (7.68E-05s)
2023-08-25 15:26:03.485 - amici.de_export - DEBUG - Finished simplifying dspline_slopesdp ++++ (9.28E-05s)
2023-08-25 15:26:03.486 - amici.de_export - DEBUG - Finished computing dspline_slopesdp +++ (3.93E-02s)
2023-08-25 15:26:03.487 - amici.de_export - DEBUG - Finished writing dspline_slopesdp.cpp ++ (4.71E-02s)
2023-08-25 15:26:03.523 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (7.05E-03s)
2023-08-25 15:26:03.538 - amici.de_export - DEBUG - Finished simplifying dwdw ++++ (3.78E-03s)
2023-08-25 15:26:03.539 - amici.de_export - DEBUG - Finished computing dwdw +++ (3.13E-02s)
2023-08-25 15:26:03.543 - amici.de_export - DEBUG - Finished writing dwdw.cpp ++ (4.30E-02s)
2023-08-25 15:26:03.588 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (1.64E-02s)
2023-08-25 15:26:03.599 - amici.de_export - DEBUG - Finished simplifying dxdotdw ++++ (4.12E-04s)
2023-08-25 15:26:03.600 - amici.de_export - DEBUG - Finished computing dxdotdw +++ (3.57E-02s)
2023-08-25 15:26:03.610 - amici.de_export - DEBUG - Finished writing dxdotdw.cpp ++ (5.61E-02s)
2023-08-25 15:26:03.642 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (1.69E-03s)
2023-08-25 15:26:03.656 - amici.de_export - DEBUG - Finished simplifying dxdotdx_explicit ++++ (1.36E-04s)
2023-08-25 15:26:03.657 - amici.de_export - DEBUG - Finished computing dxdotdx_explicit +++ (2.76E-02s)
2023-08-25 15:26:03.660 - amici.de_export - DEBUG - Finished writing dxdotdx_explicit.cpp ++ (3.91E-02s)
2023-08-25 15:26:03.691 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (1.63E-03s)
2023-08-25 15:26:03.702 - amici.de_export - DEBUG - Finished simplifying dxdotdp_explicit ++++ (2.59E-04s)
2023-08-25 15:26:03.703 - amici.de_export - DEBUG - Finished computing dxdotdp_explicit +++ (2.19E-02s)
2023-08-25 15:26:03.705 - amici.de_export - DEBUG - Finished writing dxdotdp_explicit.cpp ++ (3.14E-02s)
2023-08-25 15:26:03.747 - amici.de_export - DEBUG - Finished running smart_jacobian +++++ (2.99E-03s)
2023-08-25 15:26:03.834 - amici.de_export - DEBUG - Finished simplifying dydx +++++ (7.46E-02s)
2023-08-25 15:26:03.834 - amici.de_export - DEBUG - Finished computing dydx ++++ (9.83E-02s)
2023-08-25 15:26:03.853 - amici.de_export - DEBUG - Finished running smart_jacobian +++++ (2.80E-04s)
2023-08-25 15:26:03.865 - amici.de_export - DEBUG - Finished simplifying dydw +++++ (1.05E-04s)
2023-08-25 15:26:03.866 - amici.de_export - DEBUG - Finished computing dydw ++++ (2.15E-02s)
2023-08-25 15:26:03.951 - amici.de_export - DEBUG - Finished simplifying dydx ++++ (7.01E-02s)
2023-08-25 15:26:03.952 - amici.de_export - DEBUG - Finished computing dydx +++ (2.25E-01s)
2023-08-25 15:26:03.981 - amici.de_export - DEBUG - Finished writing dydx.cpp ++ (2.62E-01s)
2023-08-25 15:26:04.018 - amici.de_export - DEBUG - Finished running smart_jacobian +++++ (2.60E-04s)
2023-08-25 15:26:04.028 - amici.de_export - DEBUG - Finished simplifying dydp +++++ (9.56E-05s)
2023-08-25 15:26:04.028 - amici.de_export - DEBUG - Finished computing dydp ++++ (1.88E-02s)
2023-08-25 15:26:04.039 - amici.de_export - DEBUG - Finished simplifying dydp ++++ (9.50E-05s)
2023-08-25 15:26:04.040 - amici.de_export - DEBUG - Finished computing dydp +++ (4.04E-02s)
2023-08-25 15:26:04.042 - amici.de_export - DEBUG - Finished writing dydp.cpp ++ (5.07E-02s)
2023-08-25 15:26:04.061 - amici.de_export - DEBUG - Finished computing dzdx +++ (1.50E-04s)
2023-08-25 15:26:04.061 - amici.de_export - DEBUG - Finished writing dzdx.cpp ++ (8.54E-03s)
2023-08-25 15:26:04.079 - amici.de_export - DEBUG - Finished computing dzdp +++ (1.53E-04s)
2023-08-25 15:26:04.081 - amici.de_export - DEBUG - Finished writing dzdp.cpp ++ (9.44E-03s)
2023-08-25 15:26:04.098 - amici.de_export - DEBUG - Finished computing drzdx +++ (1.45E-04s)
2023-08-25 15:26:04.099 - amici.de_export - DEBUG - Finished writing drzdx.cpp ++ (9.13E-03s)
2023-08-25 15:26:04.117 - amici.de_export - DEBUG - Finished computing drzdp +++ (1.53E-04s)
2023-08-25 15:26:04.118 - amici.de_export - DEBUG - Finished writing drzdp.cpp ++ (8.88E-03s)
2023-08-25 15:26:04.144 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (3.40E-04s)
2023-08-25 15:26:04.154 - amici.de_export - DEBUG - Finished simplifying dsigmaydy ++++ (8.96E-05s)
2023-08-25 15:26:04.154 - amici.de_export - DEBUG - Finished computing dsigmaydy +++ (2.02E-02s)
2023-08-25 15:26:04.155 - amici.de_export - DEBUG - Finished writing dsigmaydy.cpp ++ (2.86E-02s)
2023-08-25 15:26:04.183 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (1.20E-03s)
2023-08-25 15:26:04.193 - amici.de_export - DEBUG - Finished simplifying dsigmaydp ++++ (1.58E-04s)
2023-08-25 15:26:04.194 - amici.de_export - DEBUG - Finished computing dsigmaydp +++ (1.88E-02s)
2023-08-25 15:26:04.197 - amici.de_export - DEBUG - Finished writing dsigmaydp.cpp ++ (2.99E-02s)
2023-08-25 15:26:04.208 - amici.de_export - DEBUG - Finished writing sigmay.cpp ++ (1.52E-03s)
2023-08-25 15:26:04.232 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (7.05E-05s)
2023-08-25 15:26:04.243 - amici.de_export - DEBUG - Finished simplifying dsigmazdp ++++ (1.14E-04s)
2023-08-25 15:26:04.244 - amici.de_export - DEBUG - Finished computing dsigmazdp +++ (1.93E-02s)
2023-08-25 15:26:04.244 - amici.de_export - DEBUG - Finished writing dsigmazdp.cpp ++ (2.79E-02s)
2023-08-25 15:26:04.256 - amici.de_export - DEBUG - Finished writing sigmaz.cpp ++ (5.24E-05s)
2023-08-25 15:26:04.272 - amici.de_export - DEBUG - Finished computing stau +++ (1.40E-04s)
2023-08-25 15:26:04.273 - amici.de_export - DEBUG - Finished writing stau.cpp ++ (8.13E-03s)
2023-08-25 15:26:04.292 - amici.de_export - DEBUG - Finished computing deltax +++ (1.35E-04s)
2023-08-25 15:26:04.293 - amici.de_export - DEBUG - Finished writing deltax.cpp ++ (8.38E-03s)
2023-08-25 15:26:04.313 - amici.de_export - DEBUG - Finished computing deltasx +++ (2.33E-04s)
2023-08-25 15:26:04.314 - amici.de_export - DEBUG - Finished writing deltasx.cpp ++ (1.03E-02s)
2023-08-25 15:26:04.332 - amici.de_export - DEBUG - Finished writing w.cpp ++ (9.05E-03s)
2023-08-25 15:26:04.343 - amici.de_export - DEBUG - Finished writing x0.cpp ++ (1.96E-03s)
2023-08-25 15:26:04.375 - amici.de_export - DEBUG - Finished simplifying x0_fixedParameters ++++ (2.13E-03s)
2023-08-25 15:26:04.376 - amici.de_export - DEBUG - Finished computing x0_fixedParameters +++ (1.20E-02s)
2023-08-25 15:26:04.380 - amici.de_export - DEBUG - Finished writing x0_fixedParameters.cpp ++ (2.49E-02s)
2023-08-25 15:26:04.409 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (2.16E-03s)
2023-08-25 15:26:04.420 - amici.de_export - DEBUG - Finished simplifying sx0 ++++ (9.81E-05s)
2023-08-25 15:26:04.421 - amici.de_export - DEBUG - Finished computing sx0 +++ (2.15E-02s)
2023-08-25 15:26:04.422 - amici.de_export - DEBUG - Finished writing sx0.cpp ++ (3.13E-02s)
2023-08-25 15:26:04.450 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (3.50E-04s)
2023-08-25 15:26:04.460 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (3.11E-04s)
2023-08-25 15:26:04.471 - amici.de_export - DEBUG - Finished simplifying sx0_fixedParameters ++++ (1.77E-04s)
2023-08-25 15:26:04.472 - amici.de_export - DEBUG - Finished computing sx0_fixedParameters +++ (2.92E-02s)
2023-08-25 15:26:04.475 - amici.de_export - DEBUG - Finished writing sx0_fixedParameters.cpp ++ (3.96E-02s)
2023-08-25 15:26:04.499 - amici.de_export - DEBUG - Finished writing xdot.cpp ++ (1.48E-02s)
2023-08-25 15:26:04.513 - amici.de_export - DEBUG - Finished writing y.cpp ++ (4.81E-03s)
2023-08-25 15:26:04.537 - amici.de_export - DEBUG - Finished simplifying x_rdata ++++ (1.97E-04s)
2023-08-25 15:26:04.538 - amici.de_export - DEBUG - Finished computing x_rdata +++ (8.23E-03s)
2023-08-25 15:26:04.540 - amici.de_export - DEBUG - Finished writing x_rdata.cpp ++ (1.77E-02s)
2023-08-25 15:26:04.566 - amici.de_export - DEBUG - Finished simplifying total_cl ++++ (1.07E-04s)
2023-08-25 15:26:04.566 - amici.de_export - DEBUG - Finished computing total_cl +++ (8.30E-03s)
2023-08-25 15:26:04.567 - amici.de_export - DEBUG - Finished writing total_cl.cpp ++ (1.67E-02s)
2023-08-25 15:26:04.592 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (8.56E-05s)
2023-08-25 15:26:04.601 - amici.de_export - DEBUG - Finished simplifying dtotal_cldp ++++ (9.78E-05s)
2023-08-25 15:26:04.602 - amici.de_export - DEBUG - Finished computing dtotal_cldp +++ (1.75E-02s)
2023-08-25 15:26:04.603 - amici.de_export - DEBUG - Finished writing dtotal_cldp.cpp ++ (2.57E-02s)
2023-08-25 15:26:04.630 - amici.de_export - DEBUG - Finished simplifying dtotal_cldx_rdata ++++ (1.07E-04s)
2023-08-25 15:26:04.631 - amici.de_export - DEBUG - Finished computing dtotal_cldx_rdata +++ (8.89E-03s)
2023-08-25 15:26:04.632 - amici.de_export - DEBUG - Finished writing dtotal_cldx_rdata.cpp ++ (1.72E-02s)
2023-08-25 15:26:04.662 - amici.de_export - DEBUG - Finished simplifying x_solver ++++ (1.60E-04s)
2023-08-25 15:26:04.663 - amici.de_export - DEBUG - Finished computing x_solver +++ (9.63E-03s)
2023-08-25 15:26:04.666 - amici.de_export - DEBUG - Finished writing x_solver.cpp ++ (2.06E-02s)
2023-08-25 15:26:04.692 - amici.de_export - DEBUG - Finished simplifying dx_rdatadx_solver ++++ (6.79E-04s)
2023-08-25 15:26:04.693 - amici.de_export - DEBUG - Finished computing dx_rdatadx_solver +++ (9.45E-03s)
2023-08-25 15:26:04.694 - amici.de_export - DEBUG - Finished writing dx_rdatadx_solver.cpp ++ (1.80E-02s)
2023-08-25 15:26:04.722 - amici.de_export - DEBUG - Finished simplifying dx_rdatadp ++++ (8.74E-04s)
2023-08-25 15:26:04.723 - amici.de_export - DEBUG - Finished computing dx_rdatadp +++ (1.08E-02s)
2023-08-25 15:26:04.725 - amici.de_export - DEBUG - Finished writing dx_rdatadp.cpp ++ (2.03E-02s)
2023-08-25 15:26:04.749 - amici.de_export - DEBUG - Finished running smart_jacobian ++++ (6.97E-05s)
2023-08-25 15:26:04.758 - amici.de_export - DEBUG - Finished simplifying dx_rdatadtcl ++++ (8.82E-05s)
2023-08-25 15:26:04.759 - amici.de_export - DEBUG - Finished computing dx_rdatadtcl +++ (1.66E-02s)
2023-08-25 15:26:04.760 - amici.de_export - DEBUG - Finished writing dx_rdatadtcl.cpp ++ (2.47E-02s)
2023-08-25 15:26:04.770 - amici.de_export - DEBUG - Finished writing z.cpp ++ (4.95E-05s)
2023-08-25 15:26:04.779 - amici.de_export - DEBUG - Finished writing rz.cpp ++ (5.57E-05s)
2023-08-25 15:26:04.812 - amici.de_export - DEBUG - Finished generating cpp code + (2.65E+00s)
2023-08-25 15:26:59.498 - amici.de_export - DEBUG - Finished compiling cpp code + (5.47E+01s)
2023-08-25 15:26:59.945 - amici.petab_import - INFO - Finished Importing PEtab model (5.85E+01s)
2023-08-25 15:26:59.954 - amici.petab_import - INFO - Successfully loaded model FullModel from pyPESTO/doc/example/amici_models/FullModel.
AMICI allows us to construct an objective function from the PEtab problem, already considering the noise distribution assumed for this model. We can also simulate the problem for a parameter with this simple setup.
[4]:
# Simulation with PEtab nominal parameter values
print("PEtab benchmark parameters")
pprint(amici.petab_objective.simulate_petab(petab_problem, amici_model))
# Simulation with specified parameter values
parameters = np.array([-1.5, -5.0, -2.2, -1.7, 5.0, 4.2, 0.5, 0.8, 0.5])
ids = list(amici_model.getParameterIds())
ids[6:] = ["sd_pSTAT5A_rel", "sd_pSTAT5B_rel", "sd_rSTAT5A_rel"]
print("Individualized parameters")
pprint(
amici.petab_objective.simulate_petab(
petab_problem,
amici_model,
problem_parameters={x_id: x_i for x_id, x_i in zip(ids, parameters)},
scaled_parameters=True,
)
)
PEtab benchmark parameters
{'edatas': [<Swig Object of type 'std::vector< amici::ExpData * >::value_type' at 0x12c917d80
condition 'model1_data1' starting at t=0.0 with custom parameter scales, constants, parameters
16x3 time-resolved datapoints
(48/48 measurements & 0/48 sigmas set)
10x0 event-resolved datapoints
(0/0 measurements & 0/0 sigmas set)
>],
'llh': -138.22199656856435,
'rdatas': [<ReturnDataView(<amici.amici.ReturnData; proxy of <Swig Object of type 'std::unique_ptr< amici::ReturnData >::pointer' at 0x12c917c60> >)>],
'sllh': None}
Individualized parameters
{'edatas': [<Swig Object of type 'std::vector< amici::ExpData * >::value_type' at 0x12c7dadf0
condition 'model1_data1' starting at t=0.0 with custom parameter scales, constants, parameters
16x3 time-resolved datapoints
(48/48 measurements & 0/48 sigmas set)
10x0 event-resolved datapoints
(0/0 measurements & 0/0 sigmas set)
>],
'llh': -185.54291970899519,
'rdatas': [<ReturnDataView(<amici.amici.ReturnData; proxy of <Swig Object of type 'std::unique_ptr< amici::ReturnData >::pointer' at 0x12c917060> >)>],
'sllh': None}
We see that to call the objective function, we need to supply the parameters in a dictionary format. This is not really suitable for parameter estimation, as e.g. optimization packages usually work with (numpy) arrays. Therefore we need to create some kind of parameter mapping.
[5]:
class Objective:
"""
A very basic implementation to an objective function for AMICI, that can call the objective function just based on the parameters.
"""
def __init__(self, petab_problem: petab.Problem, model: amici.Model):
"""Constructor for objective."""
self.petab_problem = petab_problem
self.model = model
self.x_ids = list(self.model.getParameterIds())
# nned to change the names for the last ones
self.x_ids[6:] = ["sd_pSTAT5A_rel", "sd_pSTAT5B_rel", "sd_rSTAT5A_rel"]
def x_dct(self, x: np.ndarray):
"""
Turn array of parameters to dictionary usable for objective call.
"""
return {x_id: x_i for x_id, x_i in zip(self.x_ids, x)}
def __call__(self, x: np.ndarray):
"""Call the objective function"""
return -amici.petab_objective.simulate_petab(
petab_problem,
amici_model,
problem_parameters=self.x_dct(x),
scaled_parameters=True,
)["llh"]
# Test it out
obj = Objective(petab_problem, amici_model)
pprint(obj(parameters))
185.54291970899519
Summary
We have constructed a very basic functioning objective function for this specific problem. AMICI and PEtab already simplify the workflow compared to coding the objective from scratch, however there are still some shortcomings. Some things that we have not yet considered and that would require additional code are:
What if we have multiple simulation conditions? We would have to adjust the parameter mapping, to use the correct parameters and constants for each condition, and aggregate results.
What if our system starts in a steady state? In this case we would have to preequilibrate, something that AMICI can do, but we would need to account for that additionally (and it would most likely also include an additional simulation condition).
For later analysis we would like to be able to not only get the objective function but also the residuals, something that we can change in AMICI but we would have to account for this flexibility additionally.
If we fix a parameter (i.e. keeping it constant while optimizing the remaining parameters), we would have to create a different parameter mapping (same for unfixing a parameter).
What if we want to include prior knowledge of parameters?
AMICI can also calculate sensitivities (
sllhin the above output). During parameter estimation, the inference (optimization/sampling/…) algorithm typically calls the objective function many times both with and without sensitivities. Thus, we need to implement the ability to call e.g. function value, gradient and Hessian matrix (or an approximation thereof), as well as combinations of these for efficiency.
This is most likely not the complete list but can already can get yield quite substantial coding effort. While each problem can be tackled, it is a lot of code lines, and we would need to rewrite it each time if we want to change something (or invest even more work and make the design of the objective function flexible).
In short: There is a need for a tool that can account for all these variables in the objective function formulation.
With pyPESTO
All the above is easily addressed by using pyPESTO’s objective function implementation. We support a multitude of objective functions (JAX, Aesara, AMICI, Julia, self-written). For PEtab models with AMICI, we take care of the parameter mapping, multiple simulation conditions (including preequilibration), changing between residuals and objective function, fixing parameters, and sensitivity calculation.
While there is a lot of possibility for individualization, in its most basic form creating an objective from a petab file accounting for all of the above boils down to just a few lines in pyPESTO:
[6]:
petab_yaml = f"./{model_name}/{model_name}.yaml"
petab_problem = petab.Problem.from_yaml(petab_yaml)
# import the petab problem and generate a pypesto problem
importer = pypesto.petab.PetabImporter(petab_problem)
problem = importer.create_problem()
# call the objective to get the objective function value and (additionally) the gradient
problem.objective(parameters, sensi_orders=(0, 1))
[6]:
(185.5429188951038,
array([ 4.96886729e+02, 1.50715517e-01, 4.44258325e+01, 7.14778242e+02,
-5.19647592e-05, -1.66953531e+02, -1.50846999e+02, -6.86643591e+01,
-1.59022641e+01]))
2. Optimization
After creating our objective function, we can now set up an optimization problem to find model parameters that best describe the data. For this we will need * parameter bounds (to restrict the search area; these are usually based on domain knowledge) * startpoints for the multistart local optimization (as dynamical system constrained optimization problems are in most cases non-convex) * an optimizer
Without pyPESTO
[7]:
# bounds
ub = 5 * np.ones(len(parameters))
lb = -5 * np.ones(len(parameters))
# number of starts
n_starts = 25
# draw uniformly distributed parameters within these bounds
x_guesses = np.random.random((n_starts, len(lb))) * (ub - lb) + lb
# optimize
results = []
for x0 in x_guesses:
results.append(
scipy.optimize.minimize(
obj,
x0,
bounds=zip(lb, ub),
tol=1e-12,
options={"maxfun": 500},
method="L-BFGS-B",
)
)
pprint(results)
[ message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 483.97739572289964
x: [-3.352e+00 -2.216e+00 -1.774e+00 -1.776e+00 6.234e-01
3.661e+00 1.261e+00 6.150e-01 2.480e-01]
nit: 18
jac: [ 2.192e+02 3.002e+02 7.278e+02 -2.483e+02 3.634e+02
-3.849e+01 -3.175e+01 -1.073e+03 -4.504e+02]
nfev: 520
njev: 52
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 242.1633515170673
x: [-1.690e+00 1.599e+00 4.108e+00 -2.908e+00 4.344e+00
2.980e+00 2.152e+00 1.482e+00 1.401e+00]
nit: 19
jac: [-1.918e+01 -8.125e+00 8.223e+00 -3.257e+00 -1.422e+01
-1.719e+01 3.437e+01 -1.304e+01 3.139e+01]
nfev: 520
njev: 52
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 214.20136458346468
x: [-6.979e-02 -3.372e+00 -1.727e+00 4.216e+00 -4.263e+00
4.777e+00 1.342e+00 1.486e+00 1.120e+00]
nit: 24
jac: [ 1.194e+01 -6.911e-01 -5.608e-01 -4.737e-01 -1.027e+00
-1.151e+01 -5.561e+00 1.893e+01 -1.621e+01]
nfev: 530
njev: 53
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.74599720689707
x: [ 1.320e+00 -2.696e+00 -5.719e-02 2.128e+00 -9.272e-01
-1.710e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 22
jac: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
-5.684e-06 5.684e-06 5.684e-06 0.000e+00]
nfev: 270
njev: 27
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 211.80059302429987
x: [-5.000e+00 4.982e+00 -4.710e+00 -4.963e+00 -4.994e+00
4.055e+00 1.533e+00 1.645e+00 7.582e-01]
nit: 15
jac: [ 1.543e+01 -5.858e-01 -1.036e+00 6.517e+00 -7.668e+00
-2.808e+00 1.728e-01 4.503e+00 1.388e+00]
nfev: 520
njev: 52
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 223.69767170122753
x: [-2.092e+00 -1.184e+00 4.852e+00 -3.410e+00 -4.006e-01
3.576e+00 1.544e+00 1.486e+00 1.247e+00]
nit: 35
jac: [-2.892e+00 1.587e+01 5.931e+00 -7.305e+00 -5.950e+00
4.433e+00 1.067e+01 -2.547e+01 2.397e+01]
nfev: 530
njev: 53
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 237.2482338229555
x: [-3.714e-01 4.988e+00 -5.000e+00 4.999e+00 -4.999e+00
5.000e+00 1.853e+00 1.778e+00 1.323e+00]
nit: 24
jac: [-4.974e-04 1.434e+00 1.573e-01 -1.395e-01 -2.983e-02
3.767e+00 3.086e+01 2.687e+01 3.920e+00]
nfev: 510
njev: 51
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.74599744332093
x: [ 2.623e+00 -3.868e+00 4.241e+00 6.765e-01 4.940e+00
-4.530e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 26
jac: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 5.684e-06 -2.842e-06 -1.137e-05]
nfev: 450
njev: 45
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.74599618803322
x: [-5.000e+00 -1.669e+00 4.782e+00 3.631e+00 -4.844e+00
-4.694e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 20
jac: [ 5.684e-06 0.000e+00 0.000e+00 5.684e-06 -2.842e-06
0.000e+00 -2.842e-06 2.842e-06 5.684e-06]
nfev: 430
njev: 43
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.74599529669632
x: [-5.000e+00 -5.000e+00 5.000e+00 -5.000e+00 -5.000e+00
-5.000e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 15
jac: [ 0.000e+00 0.000e+00 -0.000e+00 0.000e+00 0.000e+00
0.000e+00 -2.842e-06 5.684e-06 2.842e-06]
nfev: 390
njev: 39
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 213.3930555660625
x: [-4.976e+00 4.990e+00 4.959e+00 -4.994e+00 -5.000e+00
4.930e+00 1.550e+00 1.674e+00 7.217e-01]
nit: 14
jac: [ 1.565e+00 7.910e-02 1.781e+00 1.024e+00 1.101e+00
2.874e+00 3.555e-01 4.905e-01 -4.516e-01]
nfev: 510
njev: 51
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.7459952943729
x: [-5.000e+00 -5.000e+00 -5.000e+00 -5.000e+00 -5.000e+00
-5.000e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 15
jac: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
-2.842e-06 -2.842e-06 -2.842e-06 -2.842e-06]
nfev: 360
njev: 36
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 249.74612304678104
x: [ 5.000e+00 -5.000e+00 -5.000e+00 5.000e+00 5.000e+00
4.999e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 24
jac: [-2.842e-04 8.527e-06 -2.842e-06 0.000e+00 -8.527e-06
3.411e-04 -3.351e-03 4.059e-03 -4.007e-04]
nfev: 510
njev: 51
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.74599744228198
x: [ 2.788e+00 -3.974e+00 3.981e+00 4.062e+00 -4.665e+00
-3.281e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 14
jac: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 -5.684e-06]
nfev: 480
njev: 48
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 249.5523160565121
x: [-2.706e+00 -3.267e+00 1.257e+00 2.064e+00 1.969e-01
1.412e+00 1.874e+00 1.758e+00 1.297e+00]
nit: 15
jac: [-2.522e-01 -1.616e-01 6.545e-01 1.178e+00 4.326e-01
-2.745e-01 2.402e-01 -5.134e-02 4.700e-01]
nfev: 530
njev: 53
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 249.74599444226348
x: [-4.686e+00 2.111e+00 2.352e+00 3.564e+00 3.211e+00
3.715e-01 1.873e+00 1.759e+00 1.299e+00]
nit: 36
jac: [ 5.684e-06 5.684e-06 -5.684e-06 5.684e-06 5.684e-06
-8.527e-06 -9.237e-04 2.177e-03 -2.240e-03]
nfev: 520
njev: 52
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 235.88920140413381
x: [-3.579e+00 3.028e+00 4.795e+00 4.039e+00 -1.795e+00
3.737e+00 1.784e+00 1.808e+00 1.216e+00]
nit: 20
jac: [ 3.365e+00 3.241e+00 1.683e+00 2.220e+00 9.867e-01
-3.602e-01 2.902e+01 3.113e+01 -1.724e+01]
nfev: 610
njev: 61
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 220.02415658728512
x: [-4.987e+00 -4.995e+00 4.993e+00 -4.992e+00 4.825e+00
3.875e+00 1.449e+00 1.627e+00 1.031e+00]
nit: 15
jac: [ 2.576e+00 -6.263e-01 1.389e+00 2.389e-01 7.329e-01
9.629e-01 -4.125e-01 -6.170e-01 -1.864e+00]
nfev: 510
njev: 51
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 213.09054793943582
x: [-4.864e+00 -4.913e+00 3.790e+00 -4.928e+00 -4.981e+00
4.780e+00 1.559e+00 1.668e+00 7.252e-01]
nit: 21
jac: [ 3.923e+00 1.645e+00 1.015e+01 6.265e+00 4.520e+00
1.086e+01 2.107e+00 4.310e-01 3.321e-01]
nfev: 550
njev: 55
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 233.42558246297477
x: [-4.957e+00 4.967e+00 4.933e+00 4.973e+00 -3.787e+00
4.875e+00 1.551e+00 1.644e+00 1.318e+00]
nit: 27
jac: [-2.601e-01 -1.059e+00 6.951e-01 -9.444e-01 3.608e-01
3.734e+00 2.608e+00 -1.406e+00 3.095e+00]
nfev: 530
njev: 53
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.7459534303374
x: [-6.783e-01 -5.000e+00 5.000e+00 -5.000e+00 -6.724e-01
-2.462e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 13
jac: [-3.411e-05 2.785e-04 -1.421e-05 -1.222e-04 -3.979e-05
-8.527e-05 1.705e-05 2.558e-05 -7.390e-05]
nfev: 290
njev: 29
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
success: True
status: 0
fun: 249.7459974433216
x: [ 2.690e+00 -1.853e+00 2.859e+00 4.703e+00 4.025e+00
-4.232e+00 1.873e+00 1.759e+00 1.299e+00]
nit: 20
jac: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 5.684e-06 5.684e-06]
nfev: 250
njev: 25
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 249.4079576472867
x: [-1.928e+00 -6.320e-01 2.840e+00 2.231e+00 -2.678e+00
1.019e+00 1.870e+00 1.756e+00 1.299e+00]
nit: 10
jac: [ 1.235e+00 9.780e-01 6.752e-01 1.105e-01 2.556e+00
1.004e-01 4.811e-01 5.480e-02 -4.183e-02]
nfev: 580
njev: 58
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 200.88846803388742
x: [-8.127e-01 3.724e+00 1.937e+00 -2.592e+00 4.676e+00
3.893e+00 1.392e+00 1.065e+00 1.110e+00]
nit: 20
jac: [ 6.924e+01 4.523e+00 5.096e+00 1.620e+01 7.100e+00
-9.445e+01 9.871e+00 -4.982e+01 3.366e+01]
nfev: 690
njev: 69
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>,
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: False
status: 1
fun: 216.78983787278423
x: [ 9.620e-02 1.241e+00 3.601e+00 1.566e+00 -4.929e+00
5.000e+00 1.412e+00 1.423e+00 1.287e+00]
nit: 28
jac: [ 8.744e-01 -2.230e-02 -1.398e-01 -4.002e-03 -2.683e-01
-6.436e-01 4.931e+00 1.177e+01 -2.101e+00]
nfev: 510
njev: 51
hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>]
We might want to change the optimizer, like e.g. NLopt
[8]:
import nlopt
opt = nlopt.opt(
nlopt.LD_LBFGS, len(parameters)
) # only one of many possible options
opt.set_lower_bounds(lb)
opt.set_upper_bounds(ub)
def nlopt_objective(x, grad):
"""We need a wrapper function of the kind f(x,grad) for nlopt."""
r = obj(x)
return r
opt.set_min_objective(nlopt_objective)
results = []
for x0 in x_guesses:
try:
result = opt.optimize(x0)
except (
nlopt.RoundoffLimited,
nlopt.ForcedStop,
ValueError,
RuntimeError,
MemoryError,
) as e:
result = None
results.append(result)
pprint(results)
[array([-4.72537584, -3.8267521 , 0.17032373, 0.5309411 , 1.82691101,
-0.13738387, -2.19222726, 2.40671846, -2.17865902]),
array([-4.85782185, 0.88722054, 3.46323867, -0.42272524, 2.97501061,
-0.44685985, -1.9617645 , 4.19526924, -0.68191772]),
array([ 2.75539723, -3.45625301, -1.74304845, 4.19486719, -4.36152045,
3.20102079, 4.1666204 , 4.5071745 , 1.64830203]),
array([-0.13272644, -1.78655792, -2.05292081, 4.94102789, 0.68000657,
-0.41145952, 4.43118647, 2.86292729, -2.27641822]),
array([ 1.5071923 , 3.23408298, 4.40175342, -4.93504248, -3.68651524,
4.89047865, -3.50203955, -3.98810331, -0.60343463]),
array([ 3.59031468, -0.64508741, 4.57795683, -4.55808472, 1.45169025,
0.16191615, -0.9214029 , 1.8818166 , -2.04635126]),
array([-3.69691906, -1.11149925, -2.07266599, -1.6551983 , -1.05891694,
0.25590375, 0.87136513, -1.83339326, -2.29220816]),
array([ 2.62294851, -3.86768587, 4.24057642, 0.67648706, 4.94028248,
-4.53014752, -4.41436998, 0.48069498, 2.08662195]),
array([ 1.11811741, -1.66877199, 4.78163474, 3.63123695, -4.84414353,
-4.69389636, -4.22521978, 1.05436896, 1.66464083]),
array([ 2.10705352, -4.17386158, -4.95145244, -0.4940422 , 2.44773506,
-0.72754709, -3.38821849, 4.3015123 , -4.03270095]),
array([-1.21935294, 4.99254589, -3.5227032 , -4.57026229, -4.27577682,
1.65134668, -2.57941689, 3.3876373 , -3.08581727]),
array([-2.69338089, 1.3336723 , 4.00935726, 4.23436455, -4.97880599,
0.66011236, -0.92734049, -0.72506365, -1.95148656]),
array([ 4.59360518, -1.45672536, 2.53472283, 1.59953602, 4.74752881,
2.97708352, -1.75879731, 1.52861569, -4.47452224]),
array([ 2.78765335, -3.97396464, 3.98103304, 4.06162031, -4.66533684,
-3.28137522, 3.15208208, -2.66502967, -4.85197795]),
array([-2.01253459, -2.8480651 , 3.12268386, 1.19351138, -0.60901754,
0.29935873, 3.43245553, 4.09645236, -0.05881582]),
array([-4.51736894, 2.39365053, -0.85230688, 2.93845516, 3.92387757,
3.35653866, 4.52675063, 1.98365382, 3.80369101]),
array([-1.17029265, 3.16644483, 4.85261058, 4.35300099, 1.26174191,
0.82811007, 4.66370112, 3.96059639, 3.24314499]),
array([-2.59120601, 0.69656874, -3.88289712, -1.74846428, -1.58175173,
3.39830011, 0.0917892 , -0.85030875, -3.77417568]),
array([ 2.45781935, 1.53870162, -0.24553228, 1.49870916, -3.42788561,
4.98603203, 3.19947195, -4.22036418, 0.83316028]),
array([ 2.20258998, 4.3092804 , -2.2015135 , -1.86005028, 4.82608847,
2.24886943, -1.09100022, -1.53563431, 0.22579574]),
array([-0.53682615, -4.81637178, 4.95364701, -4.57752447, -0.60577667,
-2.88604054, 0.85270085, 0.07054205, -1.27549967]),
array([ 2.69036786, -1.85327759, 2.85858288, 4.70288006, 4.02501682,
-4.23172439, -0.72188414, 3.55067703, 4.87046828]),
array([-1.88273814, -0.60977486, 2.87775688, 2.25140401, -2.65489216,
1.01047951, 3.69127283, 4.44893301, 2.65440911]),
array([0.66867122, 3.97870907, 2.13991535, 0.2612146 , 4.9597103 ,
3.23568699, 4.22366861, 0.40266923, 3.15201217]),
array([ 4.03440768, 0.83760837, 3.55179682, 1.57898613, -4.87401594,
4.33049155, 3.79628273, 1.90990052, 1.02667127])]
We can already see that the NLopt library takes different arguments and has a different result output than scipy. In order to be able to compare them, we need to modify the code again. We would at the very least like the end objective function value, our starting value and some kind of exit message.
[9]:
results = []
for x0 in x_guesses:
try:
result = opt.optimize(x0)
msg = 'Finished Successfully.'
except (
nlopt.RoundoffLimited,
nlopt.ForcedStop,
ValueError,
RuntimeError,
MemoryError,
) as e:
result = None
msg = str(e)
res_complete = {
"x": result,
"x0": x0,
"fun": opt.last_optimum_value(),
"message": msg,
"exitflag": opt.last_optimize_result(),
}
results.append(res_complete)
pprint(results)
[{'exitflag': 1,
'fun': 1150653011.8393009,
'message': 'Finished Successfully.',
'x': array([-4.72537584, -3.8267521 , 0.17032373, 0.5309411 , 1.82691101,
-0.13738387, -2.19222726, 2.40671846, -2.17865902]),
'x0': array([-4.72537584, -3.8267521 , 0.17032373, 0.5309411 , 1.82691101,
-0.13738387, -2.19222726, 2.40671846, -2.17865902])},
{'exitflag': 1,
'fun': 373210706.087737,
'message': 'Finished Successfully.',
'x': array([-4.85782185, 0.88722054, 3.46323867, -0.42272524, 2.97501061,
-0.44685985, -1.9617645 , 4.19526924, -0.68191772]),
'x0': array([-4.85782185, 0.88722054, 3.46323867, -0.42272524, 2.97501061,
-0.44685985, -1.9617645 , 4.19526924, -0.68191772])},
{'exitflag': 1,
'fun': 425.9896608503201,
'message': 'Finished Successfully.',
'x': array([ 2.75539723, -3.45625301, -1.74304845, 4.19486719, -4.36152045,
3.20102079, 4.1666204 , 4.5071745 , 1.64830203]),
'x0': array([ 2.75539723, -3.45625301, -1.74304845, 4.19486719, -4.36152045,
3.20102079, 4.1666204 , 4.5071745 , 1.64830203])},
{'exitflag': 1,
'fun': 113150729.31030881,
'message': 'Finished Successfully.',
'x': array([-0.13272644, -1.78655792, -2.05292081, 4.94102789, 0.68000657,
-0.41145952, 4.43118647, 2.86292729, -2.27641822]),
'x0': array([-0.13272644, -1.78655792, -2.05292081, 4.94102789, 0.68000657,
-0.41145952, 4.43118647, 2.86292729, -2.27641822])},
{'exitflag': 1,
'fun': 2170773400755.13,
'message': 'Finished Successfully.',
'x': array([ 1.5071923 , 3.23408298, 4.40175342, -4.93504248, -3.68651524,
4.89047865, -3.50203955, -3.98810331, -0.60343463]),
'x0': array([ 1.5071923 , 3.23408298, 4.40175342, -4.93504248, -3.68651524,
4.89047865, -3.50203955, -3.98810331, -0.60343463])},
{'exitflag': 1,
'fun': 42320078.18806582,
'message': 'Finished Successfully.',
'x': array([ 3.59031468, -0.64508741, 4.57795683, -4.55808472, 1.45169025,
0.16191615, -0.9214029 , 1.8818166 , -2.04635126]),
'x0': array([ 3.59031468, -0.64508741, 4.57795683, -4.55808472, 1.45169025,
0.16191615, -0.9214029 , 1.8818166 , -2.04635126])},
{'exitflag': 1,
'fun': 243163763.80728397,
'message': 'Finished Successfully.',
'x': array([-3.69691906, -1.11149925, -2.07266599, -1.6551983 , -1.05891694,
0.25590375, 0.87136513, -1.83339326, -2.29220816]),
'x0': array([-3.69691906, -1.11149925, -2.07266599, -1.6551983 , -1.05891694,
0.25590375, 0.87136513, -1.83339326, -2.29220816])},
{'exitflag': 1,
'fun': 30002126964787.39,
'message': 'Finished Successfully.',
'x': array([ 2.62294851, -3.86768587, 4.24057642, 0.67648706, 4.94028248,
-4.53014752, -4.41436998, 0.48069498, 2.08662195]),
'x0': array([ 2.62294851, -3.86768587, 4.24057642, 0.67648706, 4.94028248,
-4.53014752, -4.41436998, 0.48069498, 2.08662195])},
{'exitflag': 1,
'fun': 12556009991670.39,
'message': 'Finished Successfully.',
'x': array([ 1.11811741, -1.66877199, 4.78163474, 3.63123695, -4.84414353,
-4.69389636, -4.22521978, 1.05436896, 1.66464083]),
'x0': array([ 1.11811741, -1.66877199, 4.78163474, 3.63123695, -4.84414353,
-4.69389636, -4.22521978, 1.05436896, 1.66464083])},
{'exitflag': 1,
'fun': 634294830118.3749,
'message': 'Finished Successfully.',
'x': array([ 2.10705352, -4.17386158, -4.95145244, -0.4940422 , 2.44773506,
-0.72754709, -3.38821849, 4.3015123 , -4.03270095]),
'x0': array([ 2.10705352, -4.17386158, -4.95145244, -0.4940422 , 2.44773506,
-0.72754709, -3.38821849, 4.3015123 , -4.03270095])},
{'exitflag': 1,
'fun': 9986849980.33512,
'message': 'Finished Successfully.',
'x': array([-1.21935294, 4.99254589, -3.5227032 , -4.57026229, -4.27577682,
1.65134668, -2.57941689, 3.3876373 , -3.08581727]),
'x0': array([-1.21935294, 4.99254589, -3.5227032 , -4.57026229, -4.27577682,
1.65134668, -2.57941689, 3.3876373 , -3.08581727])},
{'exitflag': 1,
'fun': 29198157.935426034,
'message': 'Finished Successfully.',
'x': array([-2.69338089, 1.3336723 , 4.00935726, 4.23436455, -4.97880599,
0.66011236, -0.92734049, -0.72506365, -1.95148656]),
'x0': array([-2.69338089, 1.3336723 , 4.00935726, 4.23436455, -4.97880599,
0.66011236, -0.92734049, -0.72506365, -1.95148656])},
{'exitflag': 1,
'fun': 2817631865279.2437,
'message': 'Finished Successfully.',
'x': array([ 4.59360518, -1.45672536, 2.53472283, 1.59953602, 4.74752881,
2.97708352, -1.75879731, 1.52861569, -4.47452224]),
'x0': array([ 4.59360518, -1.45672536, 2.53472283, 1.59953602, 4.74752881,
2.97708352, -1.75879731, 1.52861569, -4.47452224])},
{'exitflag': 1,
'fun': 16029712077044.059,
'message': 'Finished Successfully.',
'x': array([ 2.78765335, -3.97396464, 3.98103304, 4.06162031, -4.66533684,
-3.28137522, 3.15208208, -2.66502967, -4.85197795]),
'x0': array([ 2.78765335, -3.97396464, 3.98103304, 4.06162031, -4.66533684,
-3.28137522, 3.15208208, -2.66502967, -4.85197795])},
{'exitflag': 1,
'fun': 4468.484751390617,
'message': 'Finished Successfully.',
'x': array([-2.01253459, -2.8480651 , 3.12268386, 1.19351138, -0.60901754,
0.29935873, 3.43245553, 4.09645236, -0.05881582]),
'x0': array([-2.01253459, -2.8480651 , 3.12268386, 1.19351138, -0.60901754,
0.29935873, 3.43245553, 4.09645236, -0.05881582])},
{'exitflag': 1,
'fun': 426.9335166325137,
'message': 'Finished Successfully.',
'x': array([-4.51736894, 2.39365053, -0.85230688, 2.93845516, 3.92387757,
3.35653866, 4.52675063, 1.98365382, 3.80369101]),
'x0': array([-4.51736894, 2.39365053, -0.85230688, 2.93845516, 3.92387757,
3.35653866, 4.52675063, 1.98365382, 3.80369101])},
{'exitflag': 1,
'fun': 481.3231591295339,
'message': 'Finished Successfully.',
'x': array([-1.17029265, 3.16644483, 4.85261058, 4.35300099, 1.26174191,
0.82811007, 4.66370112, 3.96059639, 3.24314499]),
'x0': array([-1.17029265, 3.16644483, 4.85261058, 4.35300099, 1.26174191,
0.82811007, 4.66370112, 3.96059639, 3.24314499])},
{'exitflag': 1,
'fun': 36260050961.21613,
'message': 'Finished Successfully.',
'x': array([-2.59120601, 0.69656874, -3.88289712, -1.74846428, -1.58175173,
3.39830011, 0.0917892 , -0.85030875, -3.77417568]),
'x0': array([-2.59120601, 0.69656874, -3.88289712, -1.74846428, -1.58175173,
3.39830011, 0.0917892 , -0.85030875, -3.77417568])},
{'exitflag': 1,
'fun': 7147839056555.6,
'message': 'Finished Successfully.',
'x': array([ 2.45781935, 1.53870162, -0.24553228, 1.49870916, -3.42788561,
4.98603203, 3.19947195, -4.22036418, 0.83316028]),
'x0': array([ 2.45781935, 1.53870162, -0.24553228, 1.49870916, -3.42788561,
4.98603203, 3.19947195, -4.22036418, 0.83316028])},
{'exitflag': 1,
'fun': 37797579.29678398,
'message': 'Finished Successfully.',
'x': array([ 2.20258998, 4.3092804 , -2.2015135 , -1.86005028, 4.82608847,
2.24886943, -1.09100022, -1.53563431, 0.22579574]),
'x0': array([ 2.20258998, 4.3092804 , -2.2015135 , -1.86005028, 4.82608847,
2.24886943, -1.09100022, -1.53563431, 0.22579574])},
{'exitflag': 1,
'fun': 1146659.4928225973,
'message': 'Finished Successfully.',
'x': array([-0.53682615, -4.81637178, 4.95364701, -4.57752447, -0.60577667,
-2.88604054, 0.85270085, 0.07054205, -1.27549967]),
'x0': array([-0.53682615, -4.81637178, 4.95364701, -4.57752447, -0.60577667,
-2.88604054, 0.85270085, 0.07054205, -1.27549967])},
{'exitflag': 1,
'fun': 1236788.617249787,
'message': 'Finished Successfully.',
'x': array([ 2.69036786, -1.85327759, 2.85858288, 4.70288006, 4.02501682,
-4.23172439, -0.72188414, 3.55067703, 4.87046828]),
'x0': array([ 2.69036786, -1.85327759, 2.85858288, 4.70288006, 4.02501682,
-4.23172439, -0.72188414, 3.55067703, 4.87046828])},
{'exitflag': 1,
'fun': 441.8147467190984,
'message': 'Finished Successfully.',
'x': array([-1.88273814, -0.60977486, 2.87775688, 2.25140401, -2.65489216,
1.01047951, 3.69127283, 4.44893301, 2.65440911]),
'x0': array([-1.88273814, -0.60977486, 2.87775688, 2.25140401, -2.65489216,
1.01047951, 3.69127283, 4.44893301, 2.65440911])},
{'exitflag': 1,
'fun': 4453.420140298416,
'message': 'Finished Successfully.',
'x': array([0.66867122, 3.97870907, 2.13991535, 0.2612146 , 4.9597103 ,
3.23568699, 4.22366861, 0.40266923, 3.15201217]),
'x0': array([0.66867122, 3.97870907, 2.13991535, 0.2612146 , 4.9597103 ,
3.23568699, 4.22366861, 0.40266923, 3.15201217])},
{'exitflag': 1,
'fun': 324.16551565091083,
'message': 'Finished Successfully.',
'x': array([ 4.03440768, 0.83760837, 3.55179682, 1.57898613, -4.87401594,
4.33049155, 3.79628273, 1.90990052, 1.02667127]),
'x0': array([ 4.03440768, 0.83760837, 3.55179682, 1.57898613, -4.87401594,
4.33049155, 3.79628273, 1.90990052, 1.02667127])}]
This is smoothly running and does not take too much code. But again, we did not consider quite a few things regarding this problem:
The Optimizer internally uses finite differences, which is firstly inefficient and secondly inaccurate, especially for very stiff models. Constructing sensitivities and integrating them into the optimization can be quite tedious.
There is no tracking of the history, we only get the end points. If we want to analyze this in more detail we need to implement this into the objective function.
Many times, especcially for larger models, we might want to change the optimizer depending on the performance. For example, for some systems other optimizers might perform better, e.g. a global vs a local one. A detailed analysis on this would require some setup, and each optimizer takes arguments in a different form.
For bigger models and more starts, parallelization becomes a key component to ensure efficency.
Especially when considering multiple optimizers, the lack of a quasi standardised result format becomes apparent und thus one would either have to write a proper result class or individualize all downstream analysis for each optimizer.
With pyPESTO
Using pyPESTO, all the above is easily possible. A pypesto.engine.MultiProcessEngine allows to use parallelization, and optimize.ScipyOptimizer specifies to use a scipy based optimizer. Alternatively, e.g. try optimize.FidesOptimizer or optimize.NLoptOptimizer, all with consistent calls and output formats. The results of the single optimizer runs are filled into a unified pyPESTO result object.
[10]:
results_pypesto = optimize.minimize(
problem=problem,
optimizer=optimize.ScipyOptimizer(),
n_starts=n_starts,
engine=pypesto.engine.MultiProcessEngine(),
)
# a summary of the results
display(Markdown(results_pypesto.summary()))
Engine will use up to 8 processes (= CPU count).
0%| | 0/25 [00:00<?, ?it/s]2023-08-25 15:41:14.938 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 110.219 and h = 1.19944e-05, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:14.939 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 110.219: AMICI failed to integrate the forward problem
2023-08-25 15:41:18.777 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 77.9061 and h = 1.79864e-05, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:18.777 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 77.9061: AMICI failed to integrate the forward problem
2023-08-25 15:41:21.496 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 138.718 and h = 1.90748e-05, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:21.496 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 138.718: AMICI failed to integrate the forward problem
2023-08-25 15:41:21.570 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 138.718 and h = 1.90748e-05, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:21.570 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 138.718: AMICI failed to integrate the forward problem
2023-08-25 15:41:21.606 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 138.718 and h = 1.90748e-05, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:21.606 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 138.718: AMICI failed to integrate the forward problem
4%|█████▏ | 1/25 [00:12<05:03, 12.66s/it]2023-08-25 15:41:23.091 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 121.411 and h = 7.31397e-06, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:23.091 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 121.411: AMICI failed to integrate the forward problem
2023-08-25 15:41:23.133 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 120 and h = 4.38244e-05, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:23.133 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 120: AMICI failed to integrate the forward problem
2023-08-25 15:41:24.005 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:ERR_FAILURE] AMICI ERROR: in module CVODES in function CVode : At t = 38.28 and h = 1.71773e-06, the error test failed repeatedly or with |h| = hmin.
2023-08-25 15:41:24.006 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 38.28: AMICI failed to integrate the forward problem
44%|█████████████████████████████████████████████████████████▏ | 11/25 [00:16<00:12, 1.10it/s]2023-08-25 15:41:27.757 - amici.swig_wrappers - DEBUG - [model1_data1][CVODES:CVode:TOO_MUCH_WORK] AMICI ERROR: in module CVODES in function CVode : At t = 11.1928, mxstep steps taken before reaching tout.
2023-08-25 15:41:27.757 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 11.1928: AMICI failed to integrate the forward problem
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 25/25 [00:22<00:00, 1.11it/s]
Optimization Result
number of starts: 25
best value: 138.2224842858494, id=2
worst value: 254.03736660792256, id=6
number of non-finite values: 0
execution time summary:
Mean execution time: 4.508s
Maximum execution time: 10.909s, id=1
Minimum execution time: 0.390s, id=16
summary of optimizer messages:
Count
Message
17
CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
8
ABNORMAL_TERMINATION_IN_LNSRCH
best value found (approximately) 2 time(s)
number of plateaus found: 5
A summary of the best run:
Optimizer Result
optimizer used: <ScipyOptimizer method=L-BFGS-B options={‘maxfun’: 1000, ‘disp’: False}>
message: ABNORMAL_TERMINATION_IN_LNSRCH
number of evaluations: 141
time taken to optimize: 6.583s
startpoint: [-4.32816355 -4.74870318 -1.74009129 -1.26420992 3.53977881 0.54755365 2.64804722 1.62058431 1.57747828]
endpoint: [-1.56907929 -5. -2.2098163 -1.78589671 3.55917603 4.19771074 0.58569077 0.81885971 0.49858833]
final objective value: 138.2224842858494
final gradient value: [-0.00783036 0.05534759 0.00129469 -0.00675505 -0.00121895 0.00394696 -0.00021472 0.00294705 0.00089969]
Beyond the optimization itself, pyPESTO naturally provides various analysis and visualization routines:
[11]:
_, axes = plt.subplots(ncols=2, figsize=(12, 6), constrained_layout=True)
visualize.waterfall(results_pypesto, ax=axes[0])
visualize.parameters(results_pypesto, ax=axes[1]);
[11]:
<Axes: title={'center': 'Estimated parameters'}, xlabel='Parameter value', ylabel='Parameter'>
3. Profiling
Profile likelihood analysis allows to systematically asess parameter uncertainty and identifiability, tracing a maximum profile in the likelihood landscape starting from the optimal parameter value.
Without pyPESTO
When it comes to profiling, we have the main apparatus already prepared with a working optimizer and our objective function. We still need a wrapper around the objective function as well as the geneal setup for the profiling, which includes selecting startpoints and cutoffs. For the sake of computation time, we will limit the maximum number of steps the scipy optimizer takes to 50.
[12]:
# sort the results
results_sorted = sorted(results, key=lambda a: a["fun"])
results_sorted
[12]:
[{'x': array([ 4.03440768, 0.83760837, 3.55179682, 1.57898613, -4.87401594,
4.33049155, 3.79628273, 1.90990052, 1.02667127]),
'x0': array([ 4.03440768, 0.83760837, 3.55179682, 1.57898613, -4.87401594,
4.33049155, 3.79628273, 1.90990052, 1.02667127]),
'fun': 324.16551565091083,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.75539723, -3.45625301, -1.74304845, 4.19486719, -4.36152045,
3.20102079, 4.1666204 , 4.5071745 , 1.64830203]),
'x0': array([ 2.75539723, -3.45625301, -1.74304845, 4.19486719, -4.36152045,
3.20102079, 4.1666204 , 4.5071745 , 1.64830203]),
'fun': 425.9896608503201,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-4.51736894, 2.39365053, -0.85230688, 2.93845516, 3.92387757,
3.35653866, 4.52675063, 1.98365382, 3.80369101]),
'x0': array([-4.51736894, 2.39365053, -0.85230688, 2.93845516, 3.92387757,
3.35653866, 4.52675063, 1.98365382, 3.80369101]),
'fun': 426.9335166325137,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-1.88273814, -0.60977486, 2.87775688, 2.25140401, -2.65489216,
1.01047951, 3.69127283, 4.44893301, 2.65440911]),
'x0': array([-1.88273814, -0.60977486, 2.87775688, 2.25140401, -2.65489216,
1.01047951, 3.69127283, 4.44893301, 2.65440911]),
'fun': 441.8147467190984,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-1.17029265, 3.16644483, 4.85261058, 4.35300099, 1.26174191,
0.82811007, 4.66370112, 3.96059639, 3.24314499]),
'x0': array([-1.17029265, 3.16644483, 4.85261058, 4.35300099, 1.26174191,
0.82811007, 4.66370112, 3.96059639, 3.24314499]),
'fun': 481.3231591295339,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([0.66867122, 3.97870907, 2.13991535, 0.2612146 , 4.9597103 ,
3.23568699, 4.22366861, 0.40266923, 3.15201217]),
'x0': array([0.66867122, 3.97870907, 2.13991535, 0.2612146 , 4.9597103 ,
3.23568699, 4.22366861, 0.40266923, 3.15201217]),
'fun': 4453.420140298416,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-2.01253459, -2.8480651 , 3.12268386, 1.19351138, -0.60901754,
0.29935873, 3.43245553, 4.09645236, -0.05881582]),
'x0': array([-2.01253459, -2.8480651 , 3.12268386, 1.19351138, -0.60901754,
0.29935873, 3.43245553, 4.09645236, -0.05881582]),
'fun': 4468.484751390617,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-0.53682615, -4.81637178, 4.95364701, -4.57752447, -0.60577667,
-2.88604054, 0.85270085, 0.07054205, -1.27549967]),
'x0': array([-0.53682615, -4.81637178, 4.95364701, -4.57752447, -0.60577667,
-2.88604054, 0.85270085, 0.07054205, -1.27549967]),
'fun': 1146659.4928225973,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.69036786, -1.85327759, 2.85858288, 4.70288006, 4.02501682,
-4.23172439, -0.72188414, 3.55067703, 4.87046828]),
'x0': array([ 2.69036786, -1.85327759, 2.85858288, 4.70288006, 4.02501682,
-4.23172439, -0.72188414, 3.55067703, 4.87046828]),
'fun': 1236788.617249787,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-2.69338089, 1.3336723 , 4.00935726, 4.23436455, -4.97880599,
0.66011236, -0.92734049, -0.72506365, -1.95148656]),
'x0': array([-2.69338089, 1.3336723 , 4.00935726, 4.23436455, -4.97880599,
0.66011236, -0.92734049, -0.72506365, -1.95148656]),
'fun': 29198157.935426034,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.20258998, 4.3092804 , -2.2015135 , -1.86005028, 4.82608847,
2.24886943, -1.09100022, -1.53563431, 0.22579574]),
'x0': array([ 2.20258998, 4.3092804 , -2.2015135 , -1.86005028, 4.82608847,
2.24886943, -1.09100022, -1.53563431, 0.22579574]),
'fun': 37797579.29678398,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 3.59031468, -0.64508741, 4.57795683, -4.55808472, 1.45169025,
0.16191615, -0.9214029 , 1.8818166 , -2.04635126]),
'x0': array([ 3.59031468, -0.64508741, 4.57795683, -4.55808472, 1.45169025,
0.16191615, -0.9214029 , 1.8818166 , -2.04635126]),
'fun': 42320078.18806582,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-0.13272644, -1.78655792, -2.05292081, 4.94102789, 0.68000657,
-0.41145952, 4.43118647, 2.86292729, -2.27641822]),
'x0': array([-0.13272644, -1.78655792, -2.05292081, 4.94102789, 0.68000657,
-0.41145952, 4.43118647, 2.86292729, -2.27641822]),
'fun': 113150729.31030881,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-3.69691906, -1.11149925, -2.07266599, -1.6551983 , -1.05891694,
0.25590375, 0.87136513, -1.83339326, -2.29220816]),
'x0': array([-3.69691906, -1.11149925, -2.07266599, -1.6551983 , -1.05891694,
0.25590375, 0.87136513, -1.83339326, -2.29220816]),
'fun': 243163763.80728397,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-4.85782185, 0.88722054, 3.46323867, -0.42272524, 2.97501061,
-0.44685985, -1.9617645 , 4.19526924, -0.68191772]),
'x0': array([-4.85782185, 0.88722054, 3.46323867, -0.42272524, 2.97501061,
-0.44685985, -1.9617645 , 4.19526924, -0.68191772]),
'fun': 373210706.087737,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-4.72537584, -3.8267521 , 0.17032373, 0.5309411 , 1.82691101,
-0.13738387, -2.19222726, 2.40671846, -2.17865902]),
'x0': array([-4.72537584, -3.8267521 , 0.17032373, 0.5309411 , 1.82691101,
-0.13738387, -2.19222726, 2.40671846, -2.17865902]),
'fun': 1150653011.8393009,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-1.21935294, 4.99254589, -3.5227032 , -4.57026229, -4.27577682,
1.65134668, -2.57941689, 3.3876373 , -3.08581727]),
'x0': array([-1.21935294, 4.99254589, -3.5227032 , -4.57026229, -4.27577682,
1.65134668, -2.57941689, 3.3876373 , -3.08581727]),
'fun': 9986849980.33512,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([-2.59120601, 0.69656874, -3.88289712, -1.74846428, -1.58175173,
3.39830011, 0.0917892 , -0.85030875, -3.77417568]),
'x0': array([-2.59120601, 0.69656874, -3.88289712, -1.74846428, -1.58175173,
3.39830011, 0.0917892 , -0.85030875, -3.77417568]),
'fun': 36260050961.21613,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.10705352, -4.17386158, -4.95145244, -0.4940422 , 2.44773506,
-0.72754709, -3.38821849, 4.3015123 , -4.03270095]),
'x0': array([ 2.10705352, -4.17386158, -4.95145244, -0.4940422 , 2.44773506,
-0.72754709, -3.38821849, 4.3015123 , -4.03270095]),
'fun': 634294830118.3749,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 1.5071923 , 3.23408298, 4.40175342, -4.93504248, -3.68651524,
4.89047865, -3.50203955, -3.98810331, -0.60343463]),
'x0': array([ 1.5071923 , 3.23408298, 4.40175342, -4.93504248, -3.68651524,
4.89047865, -3.50203955, -3.98810331, -0.60343463]),
'fun': 2170773400755.13,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 4.59360518, -1.45672536, 2.53472283, 1.59953602, 4.74752881,
2.97708352, -1.75879731, 1.52861569, -4.47452224]),
'x0': array([ 4.59360518, -1.45672536, 2.53472283, 1.59953602, 4.74752881,
2.97708352, -1.75879731, 1.52861569, -4.47452224]),
'fun': 2817631865279.2437,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.45781935, 1.53870162, -0.24553228, 1.49870916, -3.42788561,
4.98603203, 3.19947195, -4.22036418, 0.83316028]),
'x0': array([ 2.45781935, 1.53870162, -0.24553228, 1.49870916, -3.42788561,
4.98603203, 3.19947195, -4.22036418, 0.83316028]),
'fun': 7147839056555.6,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 1.11811741, -1.66877199, 4.78163474, 3.63123695, -4.84414353,
-4.69389636, -4.22521978, 1.05436896, 1.66464083]),
'x0': array([ 1.11811741, -1.66877199, 4.78163474, 3.63123695, -4.84414353,
-4.69389636, -4.22521978, 1.05436896, 1.66464083]),
'fun': 12556009991670.39,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.78765335, -3.97396464, 3.98103304, 4.06162031, -4.66533684,
-3.28137522, 3.15208208, -2.66502967, -4.85197795]),
'x0': array([ 2.78765335, -3.97396464, 3.98103304, 4.06162031, -4.66533684,
-3.28137522, 3.15208208, -2.66502967, -4.85197795]),
'fun': 16029712077044.059,
'message': 'Finished Successfully.',
'exitflag': 1},
{'x': array([ 2.62294851, -3.86768587, 4.24057642, 0.67648706, 4.94028248,
-4.53014752, -4.41436998, 0.48069498, 2.08662195]),
'x0': array([ 2.62294851, -3.86768587, 4.24057642, 0.67648706, 4.94028248,
-4.53014752, -4.41436998, 0.48069498, 2.08662195]),
'fun': 30002126964787.39,
'message': 'Finished Successfully.',
'exitflag': 1}]
[14]:
results_pypesto.optimize_result[0]["x"][problem.x_free_indices][1:]
[14]:
array([-5. , -2.2098163 , -1.78589671, 3.55917603, 4.19771074,
0.58569077, 0.81885971, 0.49858833])
[16]:
# we optimimize the first parameter
# x_start = results_sorted[0]["x"][1:]
# x_fixed = results_sorted[0]["x"][0]
x_start = results_pypesto.optimize_result[0]["x"][problem.x_free_indices][1:]
x_fixed = results_pypesto.optimize_result[0]["x"][problem.x_free_indices][0]
fval_min = results_pypesto.optimize_result[0]["fval"]
# determine stepsize, ratios
stepsize = 0.05
ratio_min = 0.145
x_profile = [results_pypesto.optimize_result[0]["x"][problem.x_free_indices]]
fval_profile = [results_pypesto.optimize_result[0]["fval"]]
# set up for nlopt optimizer
opt = nlopt.opt(
nlopt.LD_LBFGS, len(parameters) - 1
) # only one of many possible options
opt.set_lower_bounds(lb[:-1])
opt.set_upper_bounds(ub[:-1])
for direction, bound in zip([-1, 1], (-5, 3)): # profile in both directions
print(f"direction: {direction}")
x0_curr = x_fixed
x_rest = x_start
run = True
while direction * (x0_curr - bound) < 0 and run:
x0_curr += stepsize * direction
# define objective for fixed parameter
def fix_obj(x: np.ndarray):
x = np.insert(x, 0, x0_curr)
return obj(x)
# define nlopt objective
def nlopt_objective(x, grad):
"""We need a wrapper function of the kind f(x,grad) for nlopt."""
r = fix_obj(x)
return r
opt.set_min_objective(nlopt_objective)
result = opt.optimize(x_rest)
# update profiles
if direction == 1:
x_profile.append(np.insert(result, 0, x0_curr))
fval_profile.append(opt.last_optimum_value())
if np.exp(fval_min - fval_profile[-1]) <= ratio_min:
run = False
if direction == -1:
x_profile.insert(0, np.insert(result, 0, x0_curr))
fval_profile.insert(0, opt.last_optimum_value())
if np.exp(fval_min - fval_profile[0]) <= ratio_min:
run = False
x_rest = result
direction: -1
direction: 1
[17]:
plt.plot(
[x[0] for x in x_profile], np.exp(np.min(fval_profile) - fval_profile)
);
This is a very basic implementation that still lacks a few things: * If we want to profile all parameters, we will want to parallelize this to save time. * We chose a very unflexible stepsize, in general we would want to be able to automatically adjust the stepsize during each profile calculation. * As this approach requires (multiple) optimizations under the hood, the things discussed in the last step also apply here mostly.
With pyPESTO
pyPESTO takes care of those things and integrates the profiling directly into the Result object
[18]:
result_pypesto = profile.parameter_profile(
problem=problem,
result=results_pypesto,
optimizer=optimize.ScipyOptimizer(),
engine=pypesto.engine.MultiProcessEngine(),
profile_index=[0],
)
visualize.profiles(result_pypesto);
Engine will use up to 8 processes (= CPU count).
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:13<00:00, 13.82s/it]
4. Sampling
pyPESTO also supports Bayesian sampling methods. These are used to retrieve posterior distributions and measure uncertainty globally.
Without pyPESTO
While there are many available sampling methods, setting them up for a more complex objective function can be time intensive, and comparing different ones even more so.
[19]:
import emcee
n_samples = 500
# set up the sampler
# rewrite nll to llh
def log_prob(x):
"""Log-probability density function."""
# check if parameter lies within bounds
if any(x < lb) or any(x > ub):
return -np.inf
# invert sign
return -1.0 * obj(x)
# def a function to get multiple startpoints for walkers
def get_epsilon_ball_initial_state(
center: np.ndarray,
lb: np.ndarray,
ub: np.ndarray,
nwalkers: int = 20,
epsilon: float = 1e-3,
):
"""Get walker initial positions as samples from an epsilon ball.
The ball is scaled in each direction according to the magnitude of the
center in that direction.
It is assumed that, because vectors are generated near a good point,
all generated vectors are evaluable, so evaluability is not checked.
Points that are generated outside the problem bounds will get shifted
to lie on the edge of the problem bounds.
Parameters
----------
center:
The center of the epsilon ball. The dimension should match the full
dimension of the pyPESTO problem. This will be returned as the
first position.
lb, ub:
Upper and lower bounds of the objective.
nwalkers:
Number of emcee walkers.
epsilon:
The relative radius of the ball. e.g., if `epsilon=0.5`
and the center of the first dimension is at 100, then the upper
and lower bounds of the epsilon ball in the first dimension will
be 150 and 50, respectively.
"""
# Epsilon ball
lb = center * (1 - epsilon)
ub = center * (1 + epsilon)
# Sample initial positions
dim = lb.size
lb = lb.reshape((1, -1))
ub = ub.reshape((1, -1))
# create uniform points in [0, 1]
xs = np.random.random((nwalkers - 1, dim))
# re-scale
xs = xs * (ub - lb) + lb
initial_state_after_first = xs
# Include `center` in initial positions
initial_state = np.row_stack(
(
center,
initial_state_after_first,
)
)
return initial_state
sampler = emcee.EnsembleSampler(
nwalkers=18, ndim=len(ub), log_prob_fn=log_prob
)
sampler.run_mcmc(
initial_state=get_epsilon_ball_initial_state(
results_sorted[0]["x"], lb, ub, 18
),
nsteps=n_samples,
skip_initial_state_check=True,
)
[19]:
State([[-1.24027462 1.78932702 -3.93716867 -1.83911125 -1.99677746 4.14486841
1.47913514 2.0021859 1.35932731]
[-1.47022978 1.7611698 4.98482712 0.85855303 -4.1989386 -4.3008923
1.8615283 1.95714177 1.35676185]
[-0.41107448 -0.14952555 -1.77644606 4.16117419 -4.81759753 4.76263479
1.60195281 1.74819962 1.4403907 ]
[-1.43602926 2.50778511 -2.30100058 -2.81499186 -4.86002175 3.36625301
1.68878632 1.88056349 1.14940453]
[-1.22391434 1.81872122 -4.926788 -1.83840725 -1.59900807 4.96243857
1.52924675 2.06366179 1.22987441]
[ 3.17131352 2.57330524 -1.1846534 -1.70751361 -2.87054037 0.42192922
1.73901292 1.86776607 1.410648 ]
[-1.48461283 2.07931539 -3.32203888 -2.5926305 -2.01782331 4.26739815
2.0839329 1.8225061 1.03565239]
[-1.93277149 -0.53931518 -1.76408703 4.33144966 4.79263617 -0.9410984
1.80687188 1.61994219 1.35782666]
[ 0.10231878 3.61432235 -3.84512502 -4.9434848 -1.90631217 4.65431699
2.07729333 1.65028872 1.19654252]
[ 4.50592151 -2.55538837 -1.16047637 4.24362302 4.53497182 1.87264848
1.88624933 1.70845149 1.22235004]
[ 4.72615409 3.13268638 -1.56100893 -4.8662477 2.02282208 -3.87082935
1.71348793 1.81644395 1.27623322]
[ 2.78834613 0.85239735 -0.21509618 2.03024593 -3.91778162 4.8823026
1.7798872 1.89429546 1.29492976]
[-0.32634656 3.31840234 -1.24790645 -4.29790084 -4.71308262 3.9882119
1.67219851 1.8025746 1.33922103]
[-0.08441014 1.99504729 -4.3086613 -2.44371181 -1.08546383 4.95857931
1.58357273 2.03714516 1.29240578]
[-0.10478905 2.40772042 -4.44534855 -3.06426882 -0.89430395 4.15788078
1.71021755 2.11709698 1.23181781]
[ 0.61026717 3.16617924 -3.2045833 -3.67833471 -2.67609702 4.98107667
1.64134768 2.04945557 1.06515929]
[ 4.80721281 -0.14817726 -3.47387807 0.65699343 2.30248275 2.93320564
1.94145041 1.85902189 1.20024436]
[-0.30164889 0.26109268 -1.84307512 3.18671824 -3.29807383 4.68070785
1.74777087 1.80071269 1.29463877]], log_prob=[-225.64207758 -252.53559047 -229.04464792 -225.0066885 -226.23100939
-253.38487017 -229.64580756 -252.46891095 -229.74162106 -250.5537262
-252.83686794 -251.71454896 -226.72542441 -228.79079296 -237.22532707
-227.92871341 -251.80959409 -232.78825374], blobs=None, random_state=('MT19937', array([2206932849, 687533236, 392309260, 3170464034, 53645069,
3010884295, 1924462243, 1739011224, 1215225621, 290578729,
3346691071, 1848570829, 23027121, 456591643, 3025351839,
44139322, 3859461820, 3384285855, 1545011441, 2880274270,
1612523433, 348209045, 2395282107, 139706992, 2541325984,
361020130, 1683022293, 3472867620, 989676495, 1333052438,
261248819, 846013908, 363225567, 1078525269, 3382521778,
1987817078, 1431689355, 919377321, 640858636, 1080089014,
3234408472, 2099893506, 3873028967, 1835169171, 806641627,
3825290061, 2135782189, 2804364627, 1288904372, 532697971,
1285750807, 3181725207, 1937910098, 3735350617, 877929555,
794118818, 531193134, 2968996371, 2235534554, 1078546710,
1699481864, 16632259, 2038009533, 4124018018, 1654549904,
1839175806, 281104275, 3001893995, 3549514596, 572512883,
775895305, 2476554611, 1078562900, 477044261, 3332147477,
1790764712, 1220166955, 1835496428, 2754893033, 1269592747,
1030059335, 2361857228, 3976443209, 3069245420, 2891322212,
777908704, 1732733343, 3104821860, 846811797, 2485970223,
717890732, 3822556252, 4038352219, 1021866056, 782933989,
3607286638, 2876106162, 1844124260, 1289090079, 771261560,
1552270256, 1354994831, 3061800544, 2727263367, 3030113580,
2186079388, 539503901, 877058179, 3425099351, 2714112648,
584347502, 448943255, 481046113, 2494146037, 1959281397,
2997223436, 580854431, 901139350, 4073689258, 2403752855,
1273639913, 17097930, 1189258404, 1129946182, 3861197036,
1187616964, 3950619282, 2894123197, 3052892285, 1794601679,
3107229605, 1154736540, 1445112066, 1281647315, 3823808737,
2464923304, 3066806796, 911645021, 3321406851, 2506397230,
3224207588, 34403862, 4121992940, 125096971, 3733411609,
2433840407, 1211748718, 692955217, 3920121066, 3170374543,
963071047, 2240583049, 2557131029, 2215007747, 1682863338,
1829007553, 188935160, 4233449025, 1142368962, 4126532027,
1540531607, 3427751919, 1553010111, 2479983119, 3408252102,
2263816213, 331359825, 3633921403, 3759892034, 292106085,
1864810289, 1140673266, 2800793353, 2838103537, 396634619,
2380262092, 558090601, 3954852938, 2356468210, 854842063,
3987873003, 1413040425, 1717097406, 2845933124, 200449670,
697004378, 2330358332, 913572043, 727824675, 2521505152,
3756628260, 1304545993, 237809106, 2921467337, 3517022909,
2809328755, 1400146847, 2513699124, 366244197, 2865045532,
185705230, 2728436123, 1264754284, 377298617, 2139695975,
2167647175, 223358529, 3465282111, 1175303169, 3186216422,
3649327174, 41779725, 1271572271, 1509599366, 3834341205,
776192713, 2664384316, 2403609316, 3263681045, 3055346811,
119641578, 1236369036, 1658776216, 2518401352, 4226029546,
3148558757, 2569699277, 2866355296, 2156478906, 1404501902,
2259574338, 2099399259, 1361291934, 3002098967, 1676689722,
802343793, 2988447027, 4257587183, 1160559483, 4259810484,
26038768, 3634335801, 3081765329, 2625613137, 3151957490,
925383249, 525896746, 2564842755, 2264351719, 1664592786,
4270323838, 3033360425, 754685161, 2610981497, 4055010380,
939595199, 551357476, 3155657354, 1972748719, 197478011,
2898800626, 1689855652, 953799410, 585253348, 375694973,
1377335697, 2538595639, 2825497566, 1340999129, 831526576,
3017026296, 1486493792, 3366584623, 57393291, 2269395590,
851853425, 1288518763, 249497874, 326769358, 1621412413,
478423386, 4228785772, 3199093009, 2834245505, 3430966499,
3276897556, 17435474, 3402869961, 2647167094, 1896074115,
3830180145, 1079813803, 1492462393, 1934793483, 2199874291,
3105650711, 2135627634, 2313133474, 1975487203, 1890372153,
4112771771, 1009532521, 4071594554, 3150015758, 4198705016,
3926942927, 1307590463, 2199556149, 1191234777, 3507715113,
2175050552, 3877421719, 1129190928, 2107289827, 3479211066,
2448609618, 804432187, 1598435854, 3338802337, 1787761744,
1428721688, 3471720360, 2655347578, 3314264648, 3027267759,
2007712732, 3733317522, 4012993888, 3517787824, 551121758,
2049597321, 3456036022, 3415694232, 3759659216, 2509150560,
2767078802, 171594234, 3992175113, 283686696, 4132055111,
1994172934, 3077263724, 2389273218, 1682293509, 1448618303,
3795182571, 3684132545, 1622325522, 3459644093, 2428584405,
415654718, 421558721, 1903663875, 3716389580, 3419812698,
3617346627, 1591072231, 2762520964, 116836745, 3639259734,
1005442451, 1461831630, 867361387, 1942784541, 1142795005,
1525588494, 1321625262, 162610824, 4008904733, 1776666739,
873008342, 3840442180, 2973938450, 4265481404, 4283339674,
2273252972, 71877482, 1390256942, 3544503825, 425620956,
3851338020, 2957518941, 445243979, 1074579722, 2688962277,
4273255105, 1546547539, 4024051829, 3945648095, 229231550,
595803490, 3758182796, 2169358100, 3500261562, 4192015134,
2183314072, 1545238201, 3103643224, 3841556466, 3855483966,
1662567278, 3143839091, 808076356, 480190800, 2688847279,
3994938844, 925302366, 2500422343, 610881158, 1984695872,
3101566415, 3452810700, 4264390600, 1896509376, 2705432340,
737630594, 843491200, 3532758010, 1025149261, 1657901107,
3198420133, 3883637990, 2870068863, 2458990462, 3855620477,
4085561001, 2402086898, 3598591303, 3550267891, 3130649350,
811095721, 3994393403, 4237031623, 4083059107, 3051463399,
3574114492, 3489500082, 1078191029, 1011531782, 3665502319,
2506534754, 3377378812, 4091943684, 3385579500, 873609207,
2952279524, 1124109539, 2561046657, 1209401355, 652418891,
146960807, 2284822124, 70957741, 218064618, 353348997,
193324864, 346234800, 2222422197, 907424622, 3028157175,
3359071299, 326033693, 1308837373, 3853624073, 941872757,
1348026446, 401040482, 1878332630, 2032502345, 3465082472,
620100896, 3561419166, 494354990, 238926942, 3590224542,
3575718072, 2671530629, 2301328592, 3229986077, 292475316,
1970818708, 3723688063, 3273180879, 1219909701, 3669876766,
3726886119, 4035180072, 3342544030, 4229704504, 2954320999,
3660720816, 3963744058, 4088207964, 787636590, 1028989741,
3551773942, 3067705925, 1879440107, 2690101453, 1476966661,
1164988387, 567866675, 4223115538, 2801780003, 784163621,
3001146061, 47857172, 3826349248, 591270366, 1038637042,
2849851035, 2179802647, 2327748806, 803249147, 1437242643,
2668896084, 887003105, 131613121, 1216052268, 1414385990,
2639415044, 2951259651, 744354232, 2078830196, 2862706838,
3251688536, 3902545329, 3578883028, 843511480, 2008248639,
3610132004, 622281062, 3765494681, 593697613, 1024899973,
2150321665, 3572264842, 3718275156, 3339033624, 789397804,
455982697, 195867210, 832452258, 1590638004, 2841209280,
1250620031, 4231398546, 2538639652, 1651308686, 4233459872,
3251288337, 1530737085, 2508960905, 819142661, 2454195021,
1499019860, 316344890, 1411618432, 1346866985, 2082162230,
1861144179, 3200584504, 1713787377, 180706102, 1331333666,
1253441295, 685235807, 1697835523, 3989857807, 2558228675,
828902009, 1580370495, 2751730402, 2538134001, 1555804373,
231859026, 818685043, 1092546692, 3623429586, 3779756715,
4050788987, 796440633, 1710608815, 2296686361, 3037349092,
1169055388, 3595308497, 268610246, 3144126922, 305091101,
3004394692, 4235572670, 141994113, 1728717716, 1992324897,
3387776119, 519323380, 4203830862, 2836686724, 1390785037,
4054831231, 3030165607, 916606003, 3053193754, 4131727760,
1575646449, 878167720, 38027722, 1743581095, 2239841900,
3572764997, 55813195, 3787178673, 3949825982, 2088303512,
3672572846, 2002937565, 1152259001, 2024262702, 3512380730,
1978640799, 689801872, 1484426853, 2228701662], dtype=uint32), 379, 0, 0.0))
[20]:
trace_x = np.array([sampler.get_chain(flat=True)])
trace_neglogpost = np.array([-sampler.get_log_prob(flat=True)])
[21]:
plt.plot(
trace_neglogpost.reshape(
9000,
),
"o",
alpha=0.05,
)
plt.ylim([240, 300]);
With pyPESTO
pyPESTO supports a number of samplers and unifies their usage, making a change of sampler comparatively easy. Instead of the below sample.AdaptiveMetropolisSampler, try e.g. also a sample.EmceeSampler (like above) or a sample.AdaptiveParallelTemperingSampler. It also unifies the result object to a certain extent to allow visualizations across samplers.
[22]:
# Sampling
sampler = sample.AdaptiveMetropolisSampler()
result_pypesto = sample.sample(
problem=problem,
sampler=sampler,
n_samples=1000,
result=result_pypesto,
)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1000/1000 [00:03<00:00, 311.04it/s]
Elapsed time: 3.9380469999998695
[23]:
# plot objective function trace
visualize.sampling_fval_traces(result_pypesto);
[24]:
visualize.sampling_1d_marginals(result_pypesto);
5. Storage
As the analysis itself is time consuming, it is neccesary to save the results for later usage. In this case it becomes even more apparent why one needs a unified result object, as otherwise saving will have to be adjusted each time one changes optimizer/sampler/profile startopint or other commonly changed things, or use an unsafe format such as pickling.
pyPESTO offers a unified result object and a very easy to use saving/loading-scheme:
[25]:
import tempfile
# create temporary file
fn = tempfile.mktemp(".h5")
# write result with write_result function.
# Choose which parts of the result object to save with
# corresponding booleans.
store.write_result(
result=result_pypesto,
filename=fn,
problem=True,
optimize=True,
sample=True,
profile=True,
)
[26]:
# Read result
result2 = store.read_result(fn, problem=True)
WARNING: You are loading a problem.
This problem is not to be used without a separately created objective.
[27]:
# plot profiles
visualize.sampling_1d_marginals(result2);
This concludes our brief rundown of a typical pyPESTO workflow and manual alternatives. In addition to what was shown here, pyPESTO provides a lot more functionality, including but not limited to visualization routines, diagnostics, model selection and hierarchical optimization. For further information, see the other example notebooks and the API documentation.