PEtab Select FAMOS

In this notebook we’ll run the model selection example from PEtab Select: specification standard and supporting software for automated model selection by Pathirana et. al with basico. This assumes you installed basico with optional PEtab support using

pip install copasi-basico[petab]

We begin by downloading the model:

import os

if not os.path.exists('famos'):
    !wget -q https://github.com/copasi/basico/raw/a4ebfa0575b3a0c569ea7274bed71657207c883e/tests/famos.zip
    !unzip -qq -o famos.zip
    !rm -f famos.zip

now import basico and petab_select to evaluate the problem. The evaluate_problem function will calibrate all models suggested by petab_select returning the best overall found model.

from petab_select import Problem
from basico.petab import evaluate_problem
import basico

the main parameter to tweak is the evaluation function, that one defines what optimization methods should be run on each suggested model. The default evaluation will run a couple of steps of a global optimization method, followed by some steps of a local optimizer. That is of coure rather arbitrary and will not yield the best results in many cases. So lets define an evaluation function here just to see how it is done:

def eval_fun():
    """Evaluation function for currently loaded model

    This function does take the currently loaded model `get_current_model()` and
    perform parameter estimation on it.

    It is supposed to return the solution as returned by `run_parameter_estimation`.

    :return: Solution of parameter estimation
    :rtype: pandas.DataFrame
    """

    # run genetic algorithm, important that update_model=True to ensure
    # that the following parameter estimation is based on the result of the genetic algorithm
    basico.run_parameter_estimation(method=basico.PE.GENETIC_ALGORITHM_SR, update_model=True,
                                    settings={'method': {
                                        'Number of Generations': 50,
                                        'Population Size': 10,
                                        'Stop after # Stalled Generations': 30
                                        }})
    # refine the result with a local optimization method
    sol = basico.run_parameter_estimation(method=basico.PE.NL2SOL, update_model=True,
                                            settings={'method': {
                                                'Iteration Limit': 2000,
                                            }}
                                          )

    return sol

that being done, lets load the petab select problem and evaluate it, depending on the function chosen above, this might take a while. A parallel version is in the works for later

problem = Problem.from_yaml('./famos/petab_select_problem.yaml')
best = evaluate_problem(problem, evaluation=eval_fun, temp_dir='./tmp_selection', delete_temp_files=False)

so lets see what we got, for this we define a print function:

def print_model(model) -> None:
    """Helper method to view model attributes."""
    print(
        f"""\
Model subspace ID: {model.model_subspace_id}
PEtab YAML location: {model.model_subspace_petab_yaml}
Custom model parameters: {model.parameters}
Model hash: {model.hash}
Model ID: {model.model_id}
{problem.criterion}: {model.get_criterion(problem.criterion, compute=False)}
Model was calibrated in iteration: {model.iteration}
"""
    )

print_model(best)

Model subspace ID: M
PEtab YAML location: famos\petab\petab_problem.yaml
Custom model parameters: {'a_0ac_k05': 1, 'a_0ac_k08': 'estimate', 'a_0ac_k12': 1, 'a_0ac_k16': 1, 'a_k05_k05k08': 1, 'a_k05_k05k12': 'estimate', 'a_k05_k05k16': 1, 'a_k08_k05k08': 1, 'a_k08_k08k12': 1, 'a_k08_k08k16': 1, 'a_k12_k05k12': 'estimate', 'a_k12_k08k12': 1, 'a_k12_k12k16': 1, 'a_k16_k05k16': 1, 'a_k16_k08k16': 1, 'a_k16_k12k16': 'estimate', 'a_k05k08_k05k08k12': 1, 'a_k05k08_k05k08k16': 1, 'a_k05k12_k05k08k12': 'estimate', 'a_k05k12_k05k12k16': 1, 'a_k05k16_k05k08k16': 1, 'a_k05k16_k05k12k16': 1, 'a_k08k12_k05k08k12': 1, 'a_k08k12_k08k12k16': 1, 'a_k08k16_k05k08k16': 1, 'a_k08k16_k08k12k16': 1, 'a_k12k16_k05k12k16': 1, 'a_k12k16_k08k12k16': 'estimate', 'a_k05k08k12_4ac': 1, 'a_k05k08k16_4ac': 1, 'a_k05k12k16_4ac': 1, 'a_k08k12k16_4ac': 'estimate'}
Model hash: M-01000100001000010010000000010001
Model ID: M-01000100001000010010000000010001
Criterion.AICC: -1708.110992474285
Model was calibrated in iteration: 11

finally lets look at the path that we took:

import pandas as pd

# Read the TSV file
summary_df = pd.read_csv('output/summary.tsv', sep='\t')
summary_df

	method	# candidates	predecessor change	current model criterion	current model	candidate changes
0	forward	17	[]	inf	['a_0ac_k08', 'a_b', 'a_k05k12_k05k08k12', 'a_...	[['a_0ac_k05'], ['a_0ac_k12'], ['a_0ac_k16'], ...
1	forward	16	['a_k05_k05k12']	-1688.480722	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k05'], ['a_0ac_k12'], ['a_0ac_k16'], ...
2	backward	16	[]	-1688.480722	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
3	backward	15	['a_k12_k08k12']	-1690.781213	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
4	backward	14	['a_k08k16_k08k12k16']	-1693.062200	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
5	backward	13	['a_k05k16_k05k12k16']	-1695.323944	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
6	backward	12	['a_k05k12_k05k12k16']	-1697.564277	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
7	backward	11	['a_k08_k08k12']	-1699.774524	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
8	backward	10	['a_k08_k05k08']	-1701.947835	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
9	backward	9	['a_k08_k08k16']	-1704.098940	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
10	backward	8	['a_k08k16_k05k08k16']	-1706.240108	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
11	backward	7	['a_k12_k12k16']	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
12	forward	23	[]	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k05'], ['a_0ac_k12'], ['a_0ac_k16'], ...
13	lateral	168	[]	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k05', 'a_k05_k05k12'], ['a_0ac_k05', ...
14	lateral	0	[]	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[]

and the expected summary:

summary_df = pd.read_csv('famos/expected_summary.tsv', sep='\t')
summary_df

	method	# candidates	predecessor change	current model criterion	current model	candidate changes
0	forward	17	[]	inf	['a_0ac_k08', 'a_b', 'a_k05k12_k05k08k12', 'a_...	[['a_0ac_k05'], ['a_0ac_k12'], ['a_0ac_k16'], ...
1	forward	16	['a_k05_k05k12']	-1688.480660	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k05'], ['a_0ac_k12'], ['a_0ac_k16'], ...
2	backward	16	[]	-1688.480660	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
3	backward	15	['a_k08_k08k16']	-1690.781203	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
4	backward	14	['a_k05k16_k05k12k16']	-1693.062190	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
5	backward	13	['a_k08_k08k12']	-1695.323742	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
6	backward	12	['a_k05k12_k05k12k16']	-1697.562387	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
7	backward	11	['a_k12_k08k12']	-1699.774480	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
8	backward	10	['a_k08k16_k08k12k16']	-1701.947757	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
9	backward	9	['a_k08_k05k08']	-1704.098938	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
10	backward	8	['a_k08k16_k05k08k16']	-1706.240108	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
11	backward	7	['a_k12_k12k16']	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k08'], ['a_k05_k05k12'], ['a_k05k12_k...
12	forward	23	[]	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k05'], ['a_0ac_k12'], ['a_0ac_k16'], ...
13	lateral	168	[]	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[['a_0ac_k05', 'a_0ac_k08'], ['a_0ac_k05', 'a_...
14	lateral	0	[]	-1708.110992	['a_0ac_k08', 'a_b', 'a_k05_k05k12', 'a_k05k12...	[]