AMICI in pyPESTO
After this notebook you can…
…create a pyPESTO problem directly from an AMICI model or through PEtab.
…perform parameter estimation of your amici model and adjust advanced settings for this.
…evaluate an optimization through basic visualizations.
…inspect parameter uncertainties through profile likelihoods and MCMC sampling.
To run optimizations and/or uncertainty analysis, we turn to pyPESTO (Parameter EStimation TOolbox for python).
pyPESTO is a Python tool for parameter estimation. It provides an interface to the model simulation tool AMICI for the simulation of Ordinary Differential Equation (ODE) models specified in the SBML format. With it, we can optimize our model parameters given measurement data, we can do uncertainty analysis via profile likelihoods and/or through sampling methods. pyPESTO provides an interface to many optimizers, global and local, such as e.g. SciPy optimizers, Fides and Pyswarm. Additionally, it interfaces samplers such as pymc, emcee and some of its own samplers.
[1]:
# import
import logging
import tempfile
from pprint import pprint
import amici
import matplotlib as mpl
import numpy as np
import petab
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
import pypesto
mpl.rcParams["figure.dpi"] = 100
mpl.rcParams["font.size"] = 18
# Set seed for reproducibility
np.random.seed(1912)
# name of the model that will also be the name of the python module
model_name = "boehm_JProteomeRes2014"
# output directory
model_output_dir = "tmp/" + model_name
1. Create a pyPESTO problem
Create a pyPESTO objective from AMICI
Before we can use AMICI to simulate our model, the SBML model needs to be translated to C++ code. This is done by amici.SbmlImporter.
[2]:
sbml_file = f"./{model_name}/{model_name}.xml"
# Create an SbmlImporter instance for our SBML model
sbml_importer = amici.SbmlImporter(sbml_file)
In this example, we want to specify fixed parameters, observables and a \(\sigma\) parameter. Unfortunately, the latter two are not part of the SBML standard. However, they can be provided to amici.SbmlImporter.sbml2amici as demonstrated in the following.
Constant parameters
Constant parameters, i.e., parameters with respect to which no sensitivities are to be computed (these are often parameters specifying a certain experimental condition) are provided as a list of parameter names.
[3]:
constant_parameters = ["ratio", "specC17"]
Observables
We used SBML’s `AssignmentRule <https://sbml.org/software/libsbml/5.18.0/docs/formatted/python-api/classlibsbml_1_1_assignment_rule.html>`__ as a non-standard way to specify Model outputs within the SBML file. These rules need to be removed prior to the model import (AMICI does at this time not support these rules). This can be easily done using amici.assignmentRules2observables().
In this example, we introduced parameters named observable_* as targets of the observable AssignmentRules. Where applicable we have observable_*_sigma parameters for \(\sigma\) parameters (see below).
[4]:
# Retrieve model output names and formulae from AssignmentRules and remove the respective rules
observables = amici.assignmentRules2observables(
sbml_importer.sbml, # the libsbml model object
filter_function=lambda variable: variable.getId().startswith("observable_")
and not variable.getId().endswith("_sigma"),
)
print("Observables:")
pprint(observables)
Observables:
{'observable_pSTAT5A_rel': {'formula': '(100 * pApB + 200 * pApA * specC17) / '
'(pApB + STAT5A * specC17 + 2 * pApA * '
'specC17)',
'name': 'observable_pSTAT5A_rel'},
'observable_pSTAT5B_rel': {'formula': '-(100 * pApB - 200 * pBpB * (specC17 - '
'1)) / (STAT5B * (specC17 - 1) - pApB + '
'2 * pBpB * (specC17 - 1))',
'name': 'observable_pSTAT5B_rel'},
'observable_rSTAT5A_rel': {'formula': '(100 * pApB + 100 * STAT5A * specC17 + '
'200 * pApA * specC17) / (2 * pApB + '
'STAT5A * specC17 + 2 * pApA * specC17 '
'- STAT5B * (specC17 - 1) - 2 * pBpB * '
'(specC17 - 1))',
'name': 'observable_rSTAT5A_rel'}}
\(\sigma\) parameters
To specify measurement noise as a parameter, we simply provide a dictionary with observable names as keys and a list of (preexisting) parameter names as values to indicate which sigma parameter is to be used for which observable.
[5]:
sigma_vals = ["sd_pSTAT5A_rel", "sd_pSTAT5B_rel", "sd_rSTAT5A_rel"]
observable_names = observables.keys()
sigmas = dict(zip(list(observable_names), sigma_vals))
pprint(sigmas)
{'observable_pSTAT5A_rel': 'sd_pSTAT5A_rel',
'observable_pSTAT5B_rel': 'sd_pSTAT5B_rel',
'observable_rSTAT5A_rel': 'sd_rSTAT5A_rel'}
Generating the module
Now we can generate the python module for our model. amici.SbmlImporter.sbml2amici will symbolically derive the sensitivity equations, generate C++ code for model simulation, and assemble the python module.
[6]:
%%time
sbml_importer.sbml2amici(
model_name,
model_output_dir,
verbose=False,
observables=observables,
constant_parameters=constant_parameters,
sigmas=sigmas,
)
CPU times: user 1.53 s, sys: 12.5 ms, total: 1.55 s
Wall time: 14.4 s
Importing the module and loading the model
If everything went well, we are ready to load the newly generated model:
[7]:
model_module = amici.import_model_module(model_name, model_output_dir)
Afterwards, we can get an instance of our model from which we can retrieve information such as parameter names:
[8]:
model = model_module.getModel()
print("Model parameters:", list(model.getParameterIds()))
print("Model outputs: ", list(model.getObservableIds()))
print("Model states: ", list(model.getStateIds()))
Model parameters: ['Epo_degradation_BaF3', 'k_exp_hetero', 'k_exp_homo', 'k_imp_hetero', 'k_imp_homo', 'k_phos', 'sd_pSTAT5A_rel', 'sd_pSTAT5B_rel', 'sd_rSTAT5A_rel']
Model outputs: ['observable_pSTAT5A_rel', 'observable_pSTAT5B_rel', 'observable_rSTAT5A_rel']
Model states: ['STAT5A', 'STAT5B', 'pApB', 'pApA', 'pBpB', 'nucpApA', 'nucpApB', 'nucpBpB']
Running simulations and analyzing results
After importing the model, we can run simulations using amici.runAmiciSimulation. This requires a Model instance and a Solver instance. But, in order go gain a value of fit, we also need to provide some data.
[9]:
# we prepare our data as it is reported in the benchmark collection
# timepoints
timepoints = np.array(
[
0.0,
2.5,
5.0,
10.0,
15.0,
20.0,
30.0,
40.0,
50.0,
60.0,
80.0,
100.0,
120.0,
160.0,
200.0,
240.0,
]
)
# measurements
meas_pSTAT5A_rel = np.array(
[
7.901073,
66.363494,
81.171324,
94.730308,
95.116483,
91.441717,
91.257099,
93.672298,
88.754233,
85.269703,
81.132395,
76.135928,
65.248059,
42.599659,
25.157798,
15.430182,
]
)
meas_pSTAT5B_rel = np.array(
[
4.596533,
29.634546,
46.043806,
81.974734,
80.571609,
79.035720,
75.672380,
71.624720,
69.062863,
67.147384,
60.899476,
54.809258,
43.981290,
29.771458,
20.089017,
10.961845,
]
)
meas_rSTAT5A_rel = np.array(
[
14.723168,
33.762342,
36.799851,
49.717602,
46.928120,
47.836575,
46.928727,
40.597753,
43.783664,
44.457388,
41.327159,
41.062733,
39.235830,
36.619461,
34.893714,
32.211077,
]
)
[10]:
benchmark_parameters = np.array(
[
-1.568917588,
-4.999704894,
-2.209698782,
-1.786006548,
4.990114009,
4.197735488,
0.585755271,
0.818982819,
0.498684404,
]
)
# set timepoints for which we want to simulate the model
model.setTimepoints(timepoints)
# set fixed parameters for which we want to simulate the model
model.setFixedParameters(np.array([0.693, 0.107]))
# set parameters to optimal values found in the benchmark collection
model.setParameterScale(amici.ParameterScaling.log10)
model.setParameters(benchmark_parameters)
# Create solver instance
solver = model.getSolver()
# Run simulation using model parameters from the benchmark collection and default solver options
rdata = amici.runAmiciSimulation(model, solver)
[11]:
# Create edata instance with dimensions and timepoints
edata = amici.ExpData(
3, # number of observables
0, # number of event outputs
0, # maximum number of events
timepoints, # timepoints
)
# set observed data
edata.setObservedData(meas_pSTAT5A_rel, 0)
edata.setObservedData(meas_pSTAT5B_rel, 1)
edata.setObservedData(meas_rSTAT5A_rel, 2)
# set standard deviations to optimal values found in the benchmark collection
edata.setObservedDataStdDev(np.array(16 * [10**0.585755271]), 0)
edata.setObservedDataStdDev(np.array(16 * [10**0.818982819]), 1)
edata.setObservedDataStdDev(np.array(16 * [10**0.498684404]), 2)
[12]:
rdata = amici.runAmiciSimulation(model, solver, edata)
print("Chi2 value reported in benchmark collection: 47.9765479")
print(f"chi2 value using AMICI: {rdata['chi2']}")
Chi2 value reported in benchmark collection: 47.9765479
chi2 value using AMICI: 47.97654319310132
Creating pyPESTO objective
We are now set up to create our pyPESTO objective. This objective is a vital part of the pyPESTO infrastructure as it provides a blackbox interface to call any predefined objective function with some parameters and evaluate it. We can easily create an AmiciObjective by supplying the model, an amici solver and the data.
Keep in mind, however, that you can use ANY function you would like for this.
[13]:
# we make some more adjustments to our model and the solver
model.requireSensitivitiesForAllParameters()
solver.setSensitivityMethod(amici.SensitivityMethod.forward)
solver.setSensitivityOrder(amici.SensitivityOrder.first)
objective = pypesto.AmiciObjective(
amici_model=model, amici_solver=solver, edatas=[edata], max_sensi_order=1
)
We can now call the objective function directly for any parameter. The value that is put out is the likelihood function. If we want to interact more with the AMICI returns, we can also return this by call and e.g., retrieve the chi2 value.
[14]:
# the generic objective call
print(f"Objective value: {objective(benchmark_parameters)}")
# a call returning the AMICI data as well
obj_call_with_dict = objective(benchmark_parameters, return_dict=True)
print(
f'Chi^2 value of the same parameters: {obj_call_with_dict["rdatas"][0]["chi2"]}'
)
Objective value: 138.22199735603812
Chi^2 value of the same parameters: 47.97654319310132
Now this makes the whole process already somewhat easier, but still, getting here took us quite some coding and effort. This will only get more complicated, the more complex the model is. Therefore, in the next part, we will show you how to bypass the tedious lines of code by using PEtab.
Create a pyPESTO problem + objective from Petab
Background on PEtab
pyPESTO supports the PEtab standard. PEtab is a data format for specifying parameter estimation problems in systems biology.
A PEtab problem consist of an SBML file, defining the model topology and a set of .tsv files, defining experimental conditions, observables, measurements and parameters (and their optimization bounds, scale, priors…). All files that make up a PEtab problem can be structured in a .yaml file. The pypesto.Objective coming from a PEtab problem corresponds to the negative-log-likelihood/negative-log-posterior distribution of the parameters.
For more details on PEtab, the interested reader is referred to PEtab’s format definition, for examples the reader is referred to the PEtab benchmark collection. For demonstration purposes, a simple model of conversion-reaction will be used as the running example throughout this notebook.
[15]:
%%capture
petab_yaml = f"./{model_name}/{model_name}.yaml"
petab_problem = petab.Problem.from_yaml(petab_yaml)
importer = pypesto.petab.PetabImporter(petab_problem)
problem = importer.create_problem(verbose=False)
Compiling amici model to folder /home/docs/checkouts/readthedocs.org/user_builds/pypesto/checkouts/stable/doc/example/amici_models/0.30.0/FullModel.
[16]:
# Check the dataframes. First the parameter dataframe
petab_problem.parameter_df.head()
[16]:
| parameterName | parameterScale | lowerBound | upperBound | nominalValue | estimate | |
|---|---|---|---|---|---|---|
| parameterId | ||||||
| Epo_degradation_BaF3 | EPO_{degradation,BaF3} | log10 | 0.00001 | 100000 | 0.026983 | 1 |
| k_exp_hetero | k_{exp,hetero} | log10 | 0.00001 | 100000 | 0.000010 | 1 |
| k_exp_homo | k_{exp,homo} | log10 | 0.00001 | 100000 | 0.006170 | 1 |
| k_imp_hetero | k_{imp,hetero} | log10 | 0.00001 | 100000 | 0.016368 | 1 |
| k_imp_homo | k_{imp,homo} | log10 | 0.00001 | 100000 | 97749.379402 | 1 |
[17]:
# Check the observable dataframe
petab_problem.observable_df.head()
[17]:
| observableName | observableFormula | noiseFormula | observableTransformation | noiseDistribution | |
|---|---|---|---|---|---|
| observableId | |||||
| pSTAT5A_rel | NaN | (100 * pApB + 200 * pApA * specC17) / (pApB + ... | noiseParameter1_pSTAT5A_rel | lin | normal |
| pSTAT5B_rel | NaN | -(100 * pApB - 200 * pBpB * (specC17 - 1)) / (... | noiseParameter1_pSTAT5B_rel | lin | normal |
| rSTAT5A_rel | NaN | (100 * pApB + 100 * STAT5A * specC17 + 200 * p... | noiseParameter1_rSTAT5A_rel | lin | normal |
[18]:
# Check the measurement dataframe
petab_problem.measurement_df.head()
[18]:
| observableId | preequilibrationConditionId | simulationConditionId | measurement | time | observableParameters | noiseParameters | datasetId | |
|---|---|---|---|---|---|---|---|---|
| 0 | pSTAT5A_rel | NaN | model1_data1 | 7.901073 | 0.0 | NaN | sd_pSTAT5A_rel | model1_data1_pSTAT5A_rel |
| 1 | pSTAT5A_rel | NaN | model1_data1 | 66.363494 | 2.5 | NaN | sd_pSTAT5A_rel | model1_data1_pSTAT5A_rel |
| 2 | pSTAT5A_rel | NaN | model1_data1 | 81.171324 | 5.0 | NaN | sd_pSTAT5A_rel | model1_data1_pSTAT5A_rel |
| 3 | pSTAT5A_rel | NaN | model1_data1 | 94.730308 | 10.0 | NaN | sd_pSTAT5A_rel | model1_data1_pSTAT5A_rel |
| 4 | pSTAT5A_rel | NaN | model1_data1 | 95.116483 | 15.0 | NaN | sd_pSTAT5A_rel | model1_data1_pSTAT5A_rel |
[19]:
# check the condition dataframe
petab_problem.condition_df.head()
[19]:
| conditionName | |
|---|---|
| conditionId | |
| model1_data1 | condition1 |
This was really straightforward. With this, we are still able to do all the same things we did before and also adjust solver setting, change the model, etc.
[20]:
# call the objective function
print(f"Objective value: {problem.objective(petab_problem.x_free_indices)}")
# change things in the model
problem.objective.amici_model.requireSensitivitiesForAllParameters()
# change solver settings
print(
f"Absolute tolerance before change: {problem.objective.amici_solver.getAbsoluteTolerance()}"
)
problem.objective.amici_solver.setAbsoluteTolerance(1e-15)
print(
f"Absolute tolerance after change: {problem.objective.amici_solver.getAbsoluteTolerance()}"
)
Objective value: 928.3017253039851
Absolute tolerance before change: 1e-16
Absolute tolerance after change: 1e-15
Now we are good to go and start the first optimization.
2. Optimization
Once setup, the optimization can be done very quickly with default settings. If needed, these settings can be highly individualized and change according to the needs of our model. In this section, we shall go over some of these settings.
Optimizer
The optimizer determines the algorithm with which we optimize our model. The main disjunction is between global and local optimizers.
pyPESTO provides an interface to many optimizers, such as Fides, ScipyOptimizers, Pyswarm and many more. For a whole list of supported optimizers with settings for each optimizer you can have a look here.
[21]:
optimizer_options = {"maxiter": 1e4, "fatol": 1e-12, "frtol": 1e-12}
optimizer = optimize.FidesOptimizer(
options=optimizer_options, verbose=logging.WARN
)
History options
In some cases, it is good to trace what the optimizer did in each step, i.e., the history. There is a multitude of options on what to report here, but the most important one is trace_record which turns the history function on and off.
[22]:
# save optimizer trace
history_options = pypesto.HistoryOptions(trace_record=True)
Startpoint method
The startpoint method describes how you want to choose your startpoints, in case you do a multistart optimization. The default here is uniform meaning that each startpoint is a uniform sample from the allowed parameter space. The other two notable options are either latin_hypercube or a self-defined function. The startpoint method is an inherent attribute of the problem and can be set there.
[23]:
problem.startpoint_method = pypesto.startpoint.uniform
Optimization options
Some further possible options for the optimization. Notably allow_failed_starts, which in case of a very complicated objective function, can help get to the desired number of optimizations when turned off. As we do not need this here, we create the default options.
[24]:
opt_options = optimize.OptimizeOptions()
opt_options
[24]:
{'allow_failed_starts': True,
'report_sres': True,
'report_hess': True,
'history_beats_optimizer': True}
Running the optimization
We now only need to decide on the number of starts as well as the engine we want to use for the optimization.
[25]:
n_starts = 20 # usually a value >= 100 should be used
engine = pypesto.engine.MultiProcessEngine()
Engine will use up to 2 processes (= CPU count).
[26]:
%%time
result = optimize.minimize(
problem=problem,
optimizer=optimizer,
n_starts=n_starts,
engine=engine,
options=opt_options,
)
2025-01-13 14:50:37.357 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 140.801 and h = 3.47052e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:37.362 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 140.801: AMICI failed to integrate the forward problem
2025-01-13 14:50:37.568 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 144.836 and h = 2.68244e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:37.571 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 144.836: AMICI failed to integrate the forward problem
2025-01-13 14:50:38.133 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 168.754 and h = 2.62935e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:38.151 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 168.754: AMICI failed to integrate the forward problem
2025-01-13 14:50:38.213 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 144.339 and h = 4.63459e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:38.216 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 144.339: AMICI failed to integrate the forward problem
2025-01-13 14:50:38 fides(WARNING) Stopping as trust region radius 1.30E-16 is smaller than machine precision.
2025-01-13 14:50:38 fides(WARNING) Stopping as trust region radius 2.16E-16 is smaller than machine precision.
2025-01-13 14:50:40.705 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 145.754 and h = 4.14821e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:40.708 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 145.754: AMICI failed to integrate the forward problem
2025-01-13 14:50:41.255 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 88.6349 and h = 1.27983e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:41.258 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 88.6349: AMICI failed to integrate the forward problem
2025-01-13 14:50:41.436 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 89.0848 and h = 1.77324e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:41.439 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 89.0848: AMICI failed to integrate the forward problem
2025-01-13 14:50:41 fides(WARNING) Stopping as trust region radius 1.11E-16 is smaller than machine precision.
2025-01-13 14:50:41 fides(WARNING) Stopping as trust region radius 1.74E-16 is smaller than machine precision.
2025-01-13 14:50:42.786 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 88.997 and h = 1.18903e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:42.789 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 88.997: AMICI failed to integrate the forward problem
2025-01-13 14:50:43 fides(WARNING) Stopping as trust region radius 2.22E-16 is smaller than machine precision.
2025-01-13 14:50:44.825 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 145.423 and h = 2.39831e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:44.829 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 145.423: AMICI failed to integrate the forward problem
2025-01-13 14:50:45.126 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 144.74 and h = 4.80279e-06, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:45.129 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 144.74: AMICI failed to integrate the forward problem
2025-01-13 14:50:45 fides(WARNING) Stopping as trust region radius 1.64E-16 is smaller than machine precision.
2025-01-13 14:50:46.735 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 168.705 and h = 3.02458e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:46.738 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 168.705: AMICI failed to integrate the forward problem
2025-01-13 14:50:47.127 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 89.0118 and h = 1.6542e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:47.130 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 89.0118: AMICI failed to integrate the forward problem
2025-01-13 14:50:47 fides(WARNING) Stopping as trust region radius 1.46E-16 is smaller than machine precision.
2025-01-13 14:50:47 fides(WARNING) Stopping as trust region radius 6.29E-17 is smaller than machine precision.
2025-01-13 14:50:47.988 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 153.65 and h = 2.77433e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:47.991 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 153.65: AMICI failed to integrate the forward problem
2025-01-13 14:50:48.287 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 138.223 and h = 3.01984e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:48.290 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 138.223: AMICI failed to integrate the forward problem
2025-01-13 14:50:48.511 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 144.718 and h = 3.96807e-06, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:48.514 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 144.718: AMICI failed to integrate the forward problem
2025-01-13 14:50:48.766 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 146.864 and h = 3.19944e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:48.769 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 146.864: AMICI failed to integrate the forward problem
2025-01-13 14:50:49 fides(WARNING) Stopping as trust region radius 9.85E-17 is smaller than machine precision.
2025-01-13 14:50:49.173 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 145.475 and h = 3.88348e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:49.177 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 145.475: AMICI failed to integrate the forward problem
2025-01-13 14:50:49 fides(WARNING) Stopping as trust region radius 1.88E-16 is smaller than machine precision.
2025-01-13 14:50:50.239 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 154.613 and h = 3.19736e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:50.242 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 154.613: AMICI failed to integrate the forward problem
2025-01-13 14:50:51 fides(WARNING) Stopping as trust region radius 6.91E-17 is smaller than machine precision.
2025-01-13 14:50:51 fides(WARNING) Stopping as trust region radius 1.02E-16 is smaller than machine precision.
2025-01-13 14:50:52.286 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 89.7475 and h = 1.89272e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:52.289 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 89.7475: AMICI failed to integrate the forward problem
2025-01-13 14:50:53.423 - amici.swig_wrappers - DEBUG - [model1_data1][cvodes:cvHandleFailure:ERR_FAILURE] At t = 89.0626 and h = 1.73567e-05, the error test failed repeatedly or with |h| = hmin.
2025-01-13 14:50:53.426 - amici.swig_wrappers - ERROR - [model1_data1][FORWARD_FAILURE] AMICI forward simulation failed at t = 89.0626: AMICI failed to integrate the forward problem
2025-01-13 14:50:53 fides(WARNING) Stopping as trust region radius 5.90E-17 is smaller than machine precision.
2025-01-13 14:50:57 fides(WARNING) Stopping as trust region radius 9.09E-17 is smaller than machine precision.
2025-01-13 14:50:58 fides(WARNING) Stopping as trust region radius 6.80E-17 is smaller than machine precision.
2025-01-13 14:50:59 fides(WARNING) Stopping as trust region radius 2.11E-16 is smaller than machine precision.
2025-01-13 14:51:01 fides(WARNING) Stopping as trust region radius 2.13E-16 is smaller than machine precision.
2025-01-13 14:51:01 fides(WARNING) Stopping as trust region radius 1.07E-16 is smaller than machine precision.
CPU times: user 153 ms, sys: 37.6 ms, total: 190 ms
Wall time: 24.9 s
Now as a first step after the optimization, we can take a look at the summary of the optimizer:
[27]:
display(Markdown(result.summary()))
Optimization Result
number of starts: 20
best value: 145.75943065729737, id=17
worst value: 249.74599566316022, id=3
number of non-finite values: 0
execution time summary:
Mean execution time: 2.425s
Maximum execution time: 6.927s, id=15
Minimum execution time: 1.141s, id=3
summary of optimizer messages:
Count
Message
18
Trust Region Radius too small to proceed
2
Converged according to fval difference
best value found (approximately) 1 time(s)
number of plateaus found: 3
A summary of the best run:
Optimizer Result
optimizer used: <FidesOptimizer hessian_update=default verbose=30 options={‘maxiter’: 10000.0, ‘fatol’: 1e-12, ‘frtol’: 1e-12}>
message: Trust Region Radius too small to proceed
number of evaluations: 146
time taken to optimize: 2.051s
startpoint: [-2.59120601 0.69656874 -3.88289712 -1.74846428 -1.58175173 3.39830011 0.0917892 -0.85030875 -3.77417568]
endpoint: [-1.5233714 -4.9993203 -2.0344514 -1.83087303 -1.67387507 3.94119952 0.5077895 0.80381877 0.7960873 ]
final objective value: 145.75943065729737
final gradient value: [-0.00443302 0.02807615 -0.00045299 -0.01152386 -0.00347778 0.00157131 -0.00055203 0.00023891 0.00028335]
final hessian value: [[ 2.37307681e+03 3.59674849e-01 2.64840295e+02 2.61570330e+03 8.75216768e+02 -8.64834960e+02 6.90934989e+01 -4.51720980e+01 -2.39111935e+01] [ 3.59674849e-01 2.95943346e-04 2.94451387e-02 3.65407106e-01 1.60738911e-01 -1.70764251e-01 -1.64607049e-04 -2.70911741e-02 -3.73919541e-02] [ 2.64840295e+02 2.94451387e-02 7.03563204e+01 2.78372534e+02 1.49677494e+02 -1.35691638e+02 -3.38439756e+00 -3.48709827e-01 3.73415043e+00] [ 2.61570330e+03 3.65407106e-01 2.78372534e+02 2.91455703e+03 9.22012363e+02 -8.64555627e+02 8.16595700e+01 -4.77356288e+01 -3.38974065e+01] [ 8.75216768e+02 1.60738911e-01 1.49677494e+02 9.22012363e+02 4.25471931e+02 -4.03058210e+02 9.20521394e+00 -1.99732940e+01 1.07760880e+01] [-8.64834960e+02 -1.70764251e-01 -1.35691638e+02 -8.64555627e+02 -4.03058210e+02 1.14618499e+03 -2.12676284e+01 -1.23881600e+01 3.36521703e+01] [ 6.90934989e+01 -1.64607049e-04 -3.38439756e+00 8.16595700e+01 9.20521394e+00 -2.12676284e+01 8.48316409e+01 0.00000000e+00 0.00000000e+00] [-4.51720980e+01 -2.70911741e-02 -3.48709827e-01 -4.77356288e+01 -1.99732940e+01 -1.23881600e+01 0.00000000e+00 8.48298197e+01 0.00000000e+00] [-2.39111935e+01 -3.73919541e-02 3.73415043e+00 -3.38974065e+01 1.07760880e+01 3.36521703e+01 0.00000000e+00 0.00000000e+00 8.48297173e+01]]
We can see some informative statistics, such as the mean execution time, best and worst values, a small table on the exit messages of the optimizer as well as detailed info on the best optimizer.
As our best start is just as good as the reported benchmark value, we shall now further inspect the result thorough some useful visualisations.
3. Optimization visualization
Model fit
Probably the most useful visualization there is, is one, where we visualize the found parameter dynamics against the measurements. This way we can see whether the fit is qualitatively and/or quantitatively good.
[28]:
ax = model_fit.visualize_optimized_model_fit(
petab_problem=petab_problem, result=result, pypesto_problem=problem
)
Waterfall plot
The waterfall plot is a visualization of the final objective function values of each start. They are sorted from small to high and then plotted. Similar values will get clustered and get the same color.
This helps to determine whether the result is reproducible and whether we reliably found a local minimum that we hope to be the global one.
[29]:
visualize.waterfall(result);
Parameter plots
To visualize the parameters, there is a multitude of options:
Parameter overview
Here we plot the parameters of all starts within their bounds. This can tell us whether some bounds are always hit and might need to be questioned and if the best starts are similar or differ amongst themselves, hinting already for some non-identifiabilities.
[30]:
visualize.parameters(result);
Parameter correlation plot
To further look into possible uncertainties, we can plot the correlation of the final points. Sometimes, pairs of parameters are dependent on each other and fixing one might solve some non-identifiability.
[31]:
visualize.parameters_correlation_matrix(result);
Parameter histogram + scatter
In case we found some dependencies and for further investigation, we can also specifically look at the histograms of certain parameters and the pairwise parameter scatter plot.
[32]:
visualize.parameter_hist(result=result, parameter_name="k_exp_hetero")
visualize.parameter_hist(result=result, parameter_name="k_imp_homo");
[33]:
visualize.optimization_scatter(result, parameter_indices=[1, 4]);
However, these visualizations are only an indicator for possible uncertainties. In the next section we turn to proper uncertainty quantification.
4. Uncertainty quantification
This mainly consists of two parts:
Profile Likelihoods
MCMC sampling
Profile likelihood
The profile likelihood uses an optimization scheme to calculate the confidence intervals for each parameter. We start with the best found parameter set of the optimization. Then in each step, we increase/decrease the parameter of interest, fix it and then run one local optimization. We do this until we either hit the bounds or reach a sufficiently bad fit.
To run the profiling, we do not need a lot of setup, as we did this already for the optimization.
[34]:
%%time
result = profile.parameter_profile(
problem=problem,
result=result,
optimizer=optimizer,
engine=engine,
profile_index=[0, 1],
)
Next guess for k_exp_hetero in direction -1 is -5.0000. Step size: -0.0007.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5325. Step size: -0.0091.
2025-01-13 14:51:04 fides(WARNING) Stopping as trust region radius 7.23E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-5.0000. Start fval 145.759410, end fval 145.759413.
Next guess for k_exp_hetero in direction 1 is -4.8993. Step size: 0.1000.
2025-01-13 14:51:04 fides(WARNING) Stopping as trust region radius 1.00E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5325. Start fval 145.863527, end fval 145.763827.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5345. Step size: -0.0020.
2025-01-13 14:51:04 fides(WARNING) Stopping as trust region radius 1.96E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5345. Start fval 145.766022, end fval 145.766020.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5370. Step size: -0.0025.
2025-01-13 14:51:05 fides(WARNING) Stopping as trust region radius 1.11E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.8993. Start fval 145.770428, end fval 145.762588.
Next guess for k_exp_hetero in direction 1 is -4.8569. Step size: 0.0424.
2025-01-13 14:51:05 fides(WARNING) Stopping as trust region radius 1.65E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5370. Start fval 145.769315, end fval 145.769307.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5400. Step size: -0.0030.
2025-01-13 14:51:05 fides(WARNING) Stopping as trust region radius 9.56E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.8569. Start fval 145.764163, end fval 145.764163.
Next guess for k_exp_hetero in direction 1 is -4.8001. Step size: 0.0568.
2025-01-13 14:51:05 fides(WARNING) Stopping as trust region radius 1.97E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5400. Start fval 145.774243, end fval 145.774240.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5437. Step size: -0.0037.
2025-01-13 14:51:05 fides(WARNING) Stopping as trust region radius 1.49E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.8001. Start fval 145.766529, end fval 145.766527.
Next guess for k_exp_hetero in direction 1 is -4.7268. Step size: 0.0733.
2025-01-13 14:51:06 fides(WARNING) Stopping as trust region radius 1.14E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5437. Start fval 145.781639, end fval 145.781594.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5482. Step size: -0.0045.
2025-01-13 14:51:06 fides(WARNING) Stopping as trust region radius 1.09E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.7268. Start fval 145.770077, end fval 145.770066.
Next guess for k_exp_hetero in direction 1 is -4.6359. Step size: 0.0910.
2025-01-13 14:51:06 fides(WARNING) Stopping as trust region radius 1.08E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.6359. Start fval 145.775382, end fval 145.775373.
Next guess for k_exp_hetero in direction 1 is -4.5359. Step size: 0.1000.
2025-01-13 14:51:07 fides(WARNING) Stopping as trust region radius 7.43E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5482. Start fval 145.792671, end fval 145.792665.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5538. Step size: -0.0055.
2025-01-13 14:51:07 fides(WARNING) Stopping as trust region radius 8.59E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.5359. Start fval 145.782640, end fval 145.782638.
Next guess for k_exp_hetero in direction 1 is -4.4359. Step size: 0.1000.
2025-01-13 14:51:07 fides(WARNING) Stopping as trust region radius 1.76E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5538. Start fval 145.809271, end fval 145.809257.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5605. Step size: -0.0068.
2025-01-13 14:51:07 fides(WARNING) Stopping as trust region radius 2.05E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.4359. Start fval 145.791777, end fval 145.791767.
Next guess for k_exp_hetero in direction 1 is -4.3359. Step size: 0.1000.
2025-01-13 14:51:08 fides(WARNING) Stopping as trust region radius 1.40E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5605. Start fval 145.834165, end fval 145.834131.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5688. Step size: -0.0083.
2025-01-13 14:51:08 fides(WARNING) Stopping as trust region radius 1.71E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.3359. Start fval 145.803255, end fval 145.803251.
Next guess for k_exp_hetero in direction 1 is -4.2359. Step size: 0.1000.
2025-01-13 14:51:08 fides(WARNING) Stopping as trust region radius 1.42E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5688. Start fval 145.871464, end fval 145.871443.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5789. Step size: -0.0101.
2025-01-13 14:51:08 fides(WARNING) Stopping as trust region radius 1.42E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.2359. Start fval 145.817684, end fval 145.817641.
Next guess for k_exp_hetero in direction 1 is -4.1359. Step size: 0.1000.
2025-01-13 14:51:09 fides(WARNING) Stopping as trust region radius 8.32E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5789. Start fval 145.927429, end fval 145.927388.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.5913. Step size: -0.0124.
2025-01-13 14:51:09 fides(WARNING) Stopping as trust region radius 1.46E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.1359. Start fval 145.835717, end fval 145.835707.
Next guess for k_exp_hetero in direction 1 is -4.0359. Step size: 0.1000.
2025-01-13 14:51:09 fides(WARNING) Stopping as trust region radius 2.16E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.5913. Start fval 146.011352, end fval 146.011273.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.6065. Step size: -0.0152.
2025-01-13 14:51:10 fides(WARNING) Stopping as trust region radius 1.68E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-4.0359. Start fval 145.858367, end fval 145.858362.
Next guess for k_exp_hetero in direction 1 is -3.9359. Step size: 0.1000.
2025-01-13 14:51:10 fides(WARNING) Stopping as trust region radius 6.59E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.6065. Start fval 146.137128, end fval 146.136921.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.6251. Step size: -0.0186.
2025-01-13 14:51:10 fides(WARNING) Stopping as trust region radius 1.05E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.9359. Start fval 145.886757, end fval 145.886753.
Next guess for k_exp_hetero in direction 1 is -3.8359. Step size: 0.1000.
2025-01-13 14:51:10 fides(WARNING) Stopping as trust region radius 6.73E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.6251. Start fval 146.325513, end fval 146.325234.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.6480. Step size: -0.0229.
2025-01-13 14:51:11 fides(WARNING) Stopping as trust region radius 1.22E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.8359. Start fval 145.922321, end fval 145.922216.
Next guess for k_exp_hetero in direction 1 is -3.7359. Step size: 0.1000.
2025-01-13 14:51:11 fides(WARNING) Stopping as trust region radius 1.72E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.6480. Start fval 146.607845, end fval 146.607069.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.6762. Step size: -0.0282.
2025-01-13 14:51:11 fides(WARNING) Stopping as trust region radius 1.40E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.7359. Start fval 145.966538, end fval 145.966536.
Next guess for k_exp_hetero in direction 1 is -3.6359. Step size: 0.1000.
2025-01-13 14:51:12 fides(WARNING) Stopping as trust region radius 7.65E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.6762. Start fval 147.030368, end fval 147.028152.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.7110. Step size: -0.0348.
2025-01-13 14:51:12 fides(WARNING) Stopping as trust region radius 8.44E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.6359. Start fval 146.021790, end fval 146.021777.
Next guess for k_exp_hetero in direction 1 is -3.5359. Step size: 0.1000.
2025-01-13 14:51:12 fides(WARNING) Stopping as trust region radius 8.13E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.7110. Start fval 147.661599, end fval 147.654658.
Next guess for Epo_degradation_BaF3 in direction -1 is -1.7539. Step size: -0.0429.
2025-01-13 14:51:12 fides(WARNING) Stopping as trust region radius 2.06E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.5359. Start fval 146.090498, end fval 146.090487.
Next guess for k_exp_hetero in direction 1 is -3.4359. Step size: 0.1000.
2025-01-13 14:51:13 fides(WARNING) Stopping as trust region radius 1.51E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.7539. Start fval 148.596834, end fval 148.566670.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.4795. Step size: 0.0439.
2025-01-13 14:51:13 fides(WARNING) Stopping as trust region radius 1.57E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.4359. Start fval 146.175820, end fval 146.175628.
Next guess for k_exp_hetero in direction 1 is -3.3359. Step size: 0.1000.
2025-01-13 14:51:13 fides(WARNING) Stopping as trust region radius 6.10E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.4795. Start fval 145.871083, end fval 145.859434.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.4695. Step size: 0.0100.
2025-01-13 14:51:14 fides(WARNING) Stopping as trust region radius 1.90E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.3359. Start fval 146.280845, end fval 146.280824.
Next guess for k_exp_hetero in direction 1 is -3.2359. Step size: 0.1000.
2025-01-13 14:51:14 fides(WARNING) Stopping as trust region radius 8.34E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.4695. Start fval 145.909429, end fval 145.908646.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.4571. Step size: 0.0125.
2025-01-13 14:51:14 fides(WARNING) Stopping as trust region radius 8.76E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.2359. Start fval 146.410281, end fval 146.410211.
Next guess for k_exp_hetero in direction 1 is -3.1359. Step size: 0.1000.
2025-01-13 14:51:14 fides(WARNING) Stopping as trust region radius 2.07E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.4571. Start fval 145.983254, end fval 145.982985.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.4416. Step size: 0.0155.
2025-01-13 14:51:15 fides(WARNING) Stopping as trust region radius 6.86E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.4416. Start fval 146.094740, end fval 146.094054.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.4224. Step size: 0.0192.
2025-01-13 14:51:15 fides(WARNING) Stopping as trust region radius 8.82E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.1359. Start fval 146.568576, end fval 146.568485.
Next guess for k_exp_hetero in direction 1 is -3.0359. Step size: 0.1000.
2025-01-13 14:51:15 fides(WARNING) Stopping as trust region radius 1.21E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.4224. Start fval 146.261310, end fval 146.259605.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.3984. Step size: 0.0240.
2025-01-13 14:51:15 fides(WARNING) Stopping as trust region radius 1.86E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-3.0359. Start fval 146.760956, end fval 146.760860.
Next guess for k_exp_hetero in direction 1 is -2.9359. Step size: 0.1000.
2025-01-13 14:51:16 fides(WARNING) Stopping as trust region radius 9.05E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.3984. Start fval 146.509574, end fval 146.505446.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.3686. Step size: 0.0298.
2025-01-13 14:51:16 fides(WARNING) Stopping as trust region radius 1.04E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-2.9359. Start fval 146.992974, end fval 146.992841.
Next guess for k_exp_hetero in direction 1 is -2.8359. Step size: 0.1000.
2025-01-13 14:51:16 fides(WARNING) Stopping as trust region radius 7.17E-17 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.3686. Start fval 146.878182, end fval 146.868765.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.3319. Step size: 0.0367.
2025-01-13 14:51:17 fides(WARNING) Stopping as trust region radius 2.16E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-2.8359. Start fval 147.270068, end fval 147.269967.
Next guess for k_exp_hetero in direction 1 is -2.7359. Step size: 0.1000.
2025-01-13 14:51:17 fides(WARNING) Stopping as trust region radius 1.16E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.3319. Start fval 147.422709, end fval 147.404354.
Next guess for Epo_degradation_BaF3 in direction 1 is -1.2874. Step size: 0.0445.
2025-01-13 14:51:17 fides(WARNING) Stopping as trust region radius 6.39E-17 is smaller than machine precision.
Optimization successful for k_exp_hetero=-2.7359. Start fval 147.597493, end fval 147.597335.
Next guess for k_exp_hetero in direction 1 is -2.6359. Step size: 0.1000.
2025-01-13 14:51:18 fides(WARNING) Stopping as trust region radius 1.29E-16 is smaller than machine precision.
Optimization successful for Epo_degradation_BaF3=-1.2874. Start fval 148.225063, end fval 148.196208.
2025-01-13 14:51:18 fides(WARNING) Stopping as trust region radius 1.27E-16 is smaller than machine precision.
Optimization successful for k_exp_hetero=-2.6359. Start fval 147.979195, end fval 147.979067.
CPU times: user 101 ms, sys: 45.4 ms, total: 147 ms
Wall time: 14.3 s
We can visualize the profiles directly
[35]:
# plot profiles
pypesto.visualize.profiles(result);
Sampling
We can use MCMC sampling to get a distribution on the posterior of the parameters. Here again, we do not need a lot of setup. We only need to define a sampler, of which pyPESTO offers a multitude.
[36]:
# Sampling
sampler = sample.AdaptiveMetropolisSampler()
result = sample.sample(
problem=problem,
sampler=sampler,
n_samples=1000,
result=result,
)
Elapsed time: 0.7460423469999995
For visualization purposes, we can visualize the trace of the objective function value, as well as a scatter plot of the parameters, just like in the optimization. We do omit the scatter plot here, as it has a very large size.
[37]:
# plot objective function trace
visualize.sampling_fval_traces(result);
/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/stable/lib/python3.11/site-packages/pypesto/visualize/sampling.py:79: UserWarning: Burn in index not found in the results, the full chain will be shown.
You may want to use, e.g., `pypesto.sample.geweke_test`.
_, params_fval, _, _, _ = get_data_to_plot(
[38]:
visualize.sampling_1d_marginals(result);
/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/stable/lib/python3.11/site-packages/pypesto/visualize/sampling.py:1293: UserWarning: Burn in index not found in the results, the full chain will be shown.
You may want to use, e.g., `pypesto.sample.geweke_test`.
nr_params, params_fval, theta_lb, theta_ub, param_names = get_data_to_plot(
5. Saving results
Lastly, the whole process took quite some time, but is not necessarily finished. It is therefore very useful, to be able to save the result as is. pyPESTO uses the HDF5 format, and with two very short commands we are able to read and write a result from and to an HDF5 file.
Save result object in HDF5 File
[39]:
# create temporary file
fn = tempfile.NamedTemporaryFile(suffix=".hdf5", delete=False)
# write result with write_result function.
# Choose which parts of the result object to save with
# corresponding booleans.
store.write_result(
result=result,
filename=fn.name,
problem=True,
optimize=True,
sample=True,
profile=True,
)
Reload results
[40]:
# Read result
result2 = store.read_result(fn, problem=True)
# close file
fn.close()
As the warning already suggests, we need to assign the problem again correctly.
[41]:
result2.problem = problem
Now we are able to quickly load the results and visualize them.
Plot (reloaded) results
[42]:
# plot profiles
pypesto.visualize.profiles(result2);
