{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "# AMICI in pyPESTO\n", "\n", "**After this notebook you can...**\n", "\n", "* ...create a pyPESTO problem directly from an AMICI model or through PEtab.\n", "* ...perform parameter estimation of your amici model and adjust advanced settings for this.\n", "* ...evaluate an optimization through basic visualizations.\n", "* ...inspect parameter uncertainties through profile likelihoods and MCMC sampling.\n", "\n", "To run optimizations and/or uncertainty analysis, we turn to pyPESTO (**P**arameter **ES**timation **TO**olbox for python).\n", "\n", "pyPESTO is a Python tool for parameter estimation. It provides an interface to the model simulation tool [AMICI](https://github.com/AMICI-dev/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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# import\n", "import logging\n", "import tempfile\n", "from pprint import pprint\n", "\n", "import amici\n", "import amici.sim.sundials as asd\n", "import matplotlib as mpl\n", "import numpy as np\n", "import petab\n", "from IPython.display import Markdown, display\n", "\n", "import pypesto\n", "import pypesto.optimize as optimize\n", "import pypesto.petab\n", "import pypesto.profile as profile\n", "import pypesto.sample as sample\n", "import pypesto.store as store\n", "import pypesto.visualize as visualize\n", "import pypesto.visualize.model_fit as model_fit\n", "\n", "mpl.rcParams[\"figure.dpi\"] = 100\n", "mpl.rcParams[\"font.size\"] = 18\n", "\n", "# Set seed for reproducibility\n", "np.random.seed(1912)\n", "\n", "\n", "# name of the model that will also be the name of the python module\n", "model_name = \"boehm_JProteomeRes2014\"\n", "\n", "# output directory\n", "model_output_dir = \"tmp/\" + model_name" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "## 1. Create a pyPESTO problem" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Create a pyPESTO objective from AMICI\n", "\n", "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`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "sbml_file = f\"./{model_name}/{model_name}.xml\"\n", "# Create an SbmlImporter instance for our SBML model\n", "sbml_importer = amici.importers.sbml.SbmlImporter(sbml_file)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": "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](https://sbml.org/). However, they can be provided to `amici.importers.sbml.SbmlImporter.sbml2amici` as demonstrated in the following." }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Fixed parameters\n", "\n", "Fixed 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 IDs." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "fixed_parameters = [\"ratio\", \"specC17\"]" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Observation model\n", "\n", "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.assignment_rules_to_observables()`.\n", "\n", "In this example, we introduced parameters named `observable_*` as targets of the observable AssignmentRules. Where applicable we have `observable_*_sigma` parameters for $\\sigma$ parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Retrieve model output names and formulae from AssignmentRules and remove the respective rules\n", "observables = amici.importers.sbml.assignment_rules_to_observables(\n", " sbml_importer.sbml_model, # the libsbml model object\n", " filter_function=lambda variable: (\n", " variable.getId().startswith(\"observable_\")\n", " and not variable.getId().endswith(\"_sigma\")\n", " ),\n", ")\n", "for observable in observables:\n", " observable.sigma = \"sd_\" + observable.id.removeprefix(\"observable_\")\n", "print(\"Observables:\")\n", "pprint(observables)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Generating the module\n", "\n", "Now we can generate the python module for our model. `SbmlImporter.sbml2amici` will symbolically derive the sensitivity equations, generate C++ code for model simulation, and assemble the python module." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "%%time\n", "\n", "sbml_importer.sbml2amici(\n", " model_name,\n", " model_output_dir,\n", " verbose=False,\n", " observation_model=observables,\n", " fixed_parameters=fixed_parameters,\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Importing the module and loading the model\n", "\n", "If everything went well, we are ready to load the newly generated model:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "model_module = amici.import_model_module(model_name, model_output_dir)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": "Afterwards, we can get an instance of our model from which we can retrieve information such as parameter names:" }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "model = model_module.get_model()\n", "\n", "print(\"Model parameters:\", list(model.get_free_parameter_ids()))\n", "print(\"Model outputs: \", list(model.get_observable_ids()))\n", "print(\"Model states: \", list(model.get_state_ids()))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Running simulations and analyzing results\n", "\n", "After importing the model, we can run simulations using `run_simulation`. 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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# we prepare our data as it is reported in the benchmark collection\n", "\n", "# timepoints\n", "timepoints = np.array(\n", " [\n", " 0.0,\n", " 2.5,\n", " 5.0,\n", " 10.0,\n", " 15.0,\n", " 20.0,\n", " 30.0,\n", " 40.0,\n", " 50.0,\n", " 60.0,\n", " 80.0,\n", " 100.0,\n", " 120.0,\n", " 160.0,\n", " 200.0,\n", " 240.0,\n", " ]\n", ")\n", "\n", "# measurements\n", "meas_pSTAT5A_rel = np.array(\n", " [\n", " 7.901073,\n", " 66.363494,\n", " 81.171324,\n", " 94.730308,\n", " 95.116483,\n", " 91.441717,\n", " 91.257099,\n", " 93.672298,\n", " 88.754233,\n", " 85.269703,\n", " 81.132395,\n", " 76.135928,\n", " 65.248059,\n", " 42.599659,\n", " 25.157798,\n", " 15.430182,\n", " ]\n", ")\n", "meas_pSTAT5B_rel = np.array(\n", " [\n", " 4.596533,\n", " 29.634546,\n", " 46.043806,\n", " 81.974734,\n", " 80.571609,\n", " 79.035720,\n", " 75.672380,\n", " 71.624720,\n", " 69.062863,\n", " 67.147384,\n", " 60.899476,\n", " 54.809258,\n", " 43.981290,\n", " 29.771458,\n", " 20.089017,\n", " 10.961845,\n", " ]\n", ")\n", "meas_rSTAT5A_rel = np.array(\n", " [\n", " 14.723168,\n", " 33.762342,\n", " 36.799851,\n", " 49.717602,\n", " 46.928120,\n", " 47.836575,\n", " 46.928727,\n", " 40.597753,\n", " 43.783664,\n", " 44.457388,\n", " 41.327159,\n", " 41.062733,\n", " 39.235830,\n", " 36.619461,\n", " 34.893714,\n", " 32.211077,\n", " ]\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "benchmark_parameters = np.array(\n", " [\n", " -1.568917588,\n", " -4.999704894,\n", " -2.209698782,\n", " -1.786006548,\n", " 4.990114009,\n", " 4.197735488,\n", " 0.585755271,\n", " 0.818982819,\n", " 0.498684404,\n", " ]\n", ")\n", "# set timepoints for which we want to simulate the model\n", "model.set_timepoints(timepoints)\n", "\n", "# set fixed parameters for which we want to simulate the model\n", "model.set_fixed_parameters(np.array([0.693, 0.107]))\n", "\n", "# set parameters to optimal values found in the benchmark collection\n", "model.set_parameter_scale(asd.ParameterScaling.log10)\n", "model.set_free_parameters(benchmark_parameters)\n", "\n", "# Create solver instance\n", "solver = model.create_solver()\n", "\n", "# Run simulation using model parameters from the benchmark collection and default solver options\n", "rdata = asd.run_simulation(model, solver)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# Create edata instance with dimensions and timepoints\n", "edata = asd.ExpData(\n", " 3, # number of observables\n", " 0, # number of event outputs\n", " 0, # maximum number of events\n", " timepoints, # timepoints\n", ")\n", "# set observed data\n", "edata.set_measurements(meas_pSTAT5A_rel, 0)\n", "edata.set_measurements(meas_pSTAT5B_rel, 1)\n", "edata.set_measurements(meas_rSTAT5A_rel, 2)\n", "\n", "# set standard deviations to optimal values found in the benchmark collection\n", "edata.set_noise_scales(np.array(16 * [10**0.585755271]), 0)\n", "edata.set_noise_scales(np.array(16 * [10**0.818982819]), 1)\n", "edata.set_noise_scales(np.array(16 * [10**0.498684404]), 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "rdata = asd.run_simulation(model, solver, edata)\n", "\n", "print(\"Chi2 value reported in benchmark collection: 47.9765479\")\n", "print(f\"chi2 value using AMICI: {rdata['chi2']}\")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Creating pyPESTO objective\n", "\n", "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.\n", "\n", "Keep in mind, however, that you can use ANY function you would like for this." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# we make some more adjustments to our model and the solver\n", "model.require_sensitivities_for_all_parameters()\n", "\n", "solver.set_sensitivity_method(asd.SensitivityMethod.forward)\n", "solver.set_sensitivity_order(asd.SensitivityOrder.first)\n", "\n", "\n", "objective = pypesto.AmiciObjective(\n", " amici_model=model, amici_solver=solver, edatas=[edata], max_sensi_order=1\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# the generic objective call\n", "print(f\"Objective value: {objective(benchmark_parameters)}\")\n", "# a call returning the AMICI data as well\n", "obj_call_with_dict = objective(benchmark_parameters, return_dict=True)\n", "print(\n", " f\"Chi^2 value of the same parameters: {obj_call_with_dict['rdatas'][0]['chi2']}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "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." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Create a pyPESTO problem + objective from Petab" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Background on PEtab\n", "\n", "\"PEtab\n", "\n", "pyPESTO supports the [PEtab](https://github.com/PEtab-dev/PEtab) standard. PEtab is a data format for specifying parameter estimation problems in systems biology.\n", "\n", "A PEtab problem consist of an [SBML](https://sbml.org) 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.\n", "\n", "For more details on PEtab, the interested reader is referred to [PEtab's format definition](https://petab.readthedocs.io/en/latest/documentation_data_format.html), for examples the reader is referred to the [PEtab benchmark collection](https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab). For demonstration purposes, a simple model of conversion-reaction will be used as the running example throughout this notebook.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "%%capture\n", "\n", "petab_yaml = f\"./{model_name}/{model_name}.yaml\"\n", "\n", "petab_problem = petab.Problem.from_yaml(petab_yaml)\n", "importer = pypesto.petab.PetabImporter(petab_problem)\n", "problem = importer.create_problem(verbose=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# Check the dataframes. First the parameter dataframe\n", "petab_problem.parameter_df.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# Check the observable dataframe\n", "petab_problem.observable_df.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# Check the measurement dataframe\n", "petab_problem.measurement_df.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# check the condition dataframe\n", "petab_problem.condition_df.head()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# call the objective function\n", "print(f\"Objective value: {problem.objective(petab_problem.x_free_indices)}\")\n", "# change things in the model\n", "problem.objective.amici_model.require_sensitivities_for_all_parameters()\n", "# change solver settings\n", "print(\n", " f\"Absolute tolerance before change: {problem.objective.amici_solver.get_absolute_tolerance()}\"\n", ")\n", "problem.objective.amici_solver.set_absolute_tolerance(1e-15)\n", "print(\n", " f\"Absolute tolerance after change: {problem.objective.amici_solver.get_absolute_tolerance()}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "Now we are good to go and start the first optimization." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "## 2. Optimization\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Optimizer\n", "\n", "The optimizer determines the algorithm with which we optimize our model. The main disjunction is between global and local optimizers.\n", "\n", "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](https://pypesto.readthedocs.io/en/latest/api.html#optimize)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "optimizer_options = {\"maxiter\": 1e4, \"fatol\": 1e-12, \"frtol\": 1e-12}\n", "\n", "optimizer = optimize.FidesOptimizer(\n", " options=optimizer_options, verbose=logging.WARN\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### History options\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# save optimizer trace\n", "history_options = pypesto.HistoryOptions(trace_record=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Startpoint method\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "problem.startpoint_method = pypesto.startpoint.uniform" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Optimization options\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "opt_options = optimize.OptimizeOptions()\n", "opt_options" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Running the optimization\n", "\n", "We now only need to decide on the number of starts as well as the engine we want to use for the optimization." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "n_starts = 20 # usually a value >= 100 should be used\n", "engine = pypesto.engine.MultiProcessEngine()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "%%time\n", "result = optimize.minimize(\n", " problem=problem,\n", " optimizer=optimizer,\n", " n_starts=n_starts,\n", " engine=engine,\n", " options=opt_options,\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "Now as a first step after the optimization, we can take a look at the summary of the optimizer:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "display(Markdown(result.summary()))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "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.\n", "\n", "As our best start is just as good as the reported benchmark value, we shall now further inspect the result thorough some useful visualisations." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "## 3. Optimization visualization" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Model fit\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "ax = model_fit.visualize_optimized_model_fit(\n", " petab_problem=petab_problem, result=result, pypesto_problem=problem\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Waterfall plot\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "visualize.waterfall(result);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Parameter plots\n", "\n", "To visualize the parameters, there is a multitude of options:" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Parameter overview\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "visualize.parameters(result);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Parameter correlation plot\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "visualize.parameters_correlation_matrix(result);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Parameter histogram + scatter\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "visualize.parameter_hist(result=result, parameter_name=\"k_exp_hetero\")\n", "visualize.parameter_hist(result=result, parameter_name=\"k_imp_homo\");" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "visualize.optimization_scatter(result, parameter_indices=[1, 4]);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "However, these visualizations are only an indicator for possible uncertainties. In the next section we turn to proper uncertainty quantification." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "## 4. Uncertainty quantification\n", "\n", "This mainly consists of two parts:\n", "* Profile Likelihoods\n", "* MCMC sampling" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Profile likelihood\n", "\n", "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.\n", "\n", "To run the profiling, we do not need a lot of setup, as we did this already for the optimization." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "%%time\n", "\n", "result = profile.parameter_profile(\n", " problem=problem,\n", " result=result,\n", " optimizer=optimizer,\n", " engine=engine,\n", " profile_index=[0, 1],\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "We can visualize the profiles directly" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# plot profiles\n", "pypesto.visualize.profiles(result);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Sampling\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# Sampling\n", "sampler = sample.AdaptiveMetropolisSampler()\n", "result = sample.sample(\n", " problem=problem,\n", " sampler=sampler,\n", " n_samples=1000,\n", " result=result,\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# plot objective function trace\n", "visualize.sampling_fval_traces(result);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "visualize.sampling_1d_marginals(result);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "## 5. Saving results\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Save result object in HDF5 File" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# create temporary file\n", "fn = tempfile.NamedTemporaryFile(suffix=\".hdf5\", delete=False)\n", "# write result with write_result function.\n", "# Choose which parts of the result object to save with\n", "# corresponding booleans.\n", "store.write_result(\n", " result=result,\n", " filename=fn.name,\n", " problem=True,\n", " optimize=True,\n", " sample=True,\n", " profile=True,\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Reload results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# Read result\n", "result2 = store.read_result(fn, problem=True)\n", "\n", "# close file\n", "fn.close()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "As the warning already suggests, we need to assign the problem again correctly." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "result2.problem = problem" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "Now we are able to quickly load the results and visualize them." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%% md\n" } }, "source": [ "### Plot (reloaded) results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# plot profiles\n", "pypesto.visualize.profiles(result2);" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 4 }