Conversion reaction

# install if not done yet
# !apt install libatlas-base-dev swig
# %pip install pypesto[amici] --quiet

import importlib
import os
import sys

import amici
import amici.plotting
import numpy as np

import pypesto
import pypesto.optimize as optimize
import pypesto.visualize as visualize

# sbml file we want to import
sbml_file = "conversion_reaction/model_conversion_reaction.xml"
# name of the model that will also be the name of the python module
model_name = "model_conversion_reaction"
# directory to which the generated model code is written
model_output_dir = "tmp/" + model_name

Compile AMICI model

# import sbml model, compile and generate amici module
sbml_importer = amici.SbmlImporter(sbml_file)
sbml_importer.sbml2amici(model_name, model_output_dir, verbose=False)

Load AMICI model

# load amici module (the usual starting point later for the analysis)
sys.path.insert(0, os.path.abspath(model_output_dir))
model_module = importlib.import_module(model_name)
model = model_module.getModel()
model.requireSensitivitiesForAllParameters()
model.setTimepoints(np.linspace(0, 10, 11))
model.setParameterScale(amici.ParameterScaling.log10)
model.setParameters([-0.3, -0.7])
solver = model.getSolver()
solver.setSensitivityMethod(amici.SensitivityMethod.forward)
solver.setSensitivityOrder(amici.SensitivityOrder.first)

# how to run amici now:
rdata = amici.runAmiciSimulation(model, solver, None)
amici.plotting.plotStateTrajectories(rdata)
edata = amici.ExpData(rdata, 0.2, 0.0)

Optimize

%%time
# create objective function from amici model
# pesto.AmiciObjective is derived from pesto.Objective,
# the general pesto objective function class
objective = pypesto.AmiciObjective(model, solver, [edata], 1)

# create optimizer object which contains all information for doing the optimization
optimizer = optimize.ScipyOptimizer(method="ls_trf")

# create problem object containing all information on the problem to be solved
problem = pypesto.Problem(objective=objective, lb=[-2, -2], ub=[2, 2])

# do the optimization
result = optimize.minimize(
    problem=problem, optimizer=optimizer, n_starts=10, filename=None
)

CPU times: user 338 ms, sys: 7.82 ms, total: 346 ms
Wall time: 346 ms

Visualize

visualize.waterfall(result)
visualize.parameters(result)
visualize.optimization_scatter(result=result)

<seaborn.axisgrid.PairGrid at 0x7f7be571fb50>

Profiles

import pypesto.profile as profile

profile_options = profile.ProfileOptions(
    min_step_size=0.0005,
    delta_ratio_max=0.05,
    default_step_size=0.005,
    ratio_min=0.01,
)

result = profile.parameter_profile(
    problem=problem,
    result=result,
    optimizer=optimizer,
    profile_index=np.array([0,1]),
    result_index=0,
    profile_options=profile_options,
    filename=None,
)

Next guess for k1 in direction -1 is -0.3051. Step size: -0.0245.
Optimization successful for k1=-0.3051. Start fval -8.902813, end fval -8.938359.
Next guess for k1 in direction -1 is -0.3106. Step size: -0.0055.
Optimization successful for k1=-0.3106. Start fval -8.930656, end fval -8.930657.
Next guess for k1 in direction -1 is -0.3172. Step size: -0.0067.
Optimization successful for k1=-0.3172. Start fval -8.919107, end fval -8.919107.
Next guess for k1 in direction -1 is -0.3253. Step size: -0.0081.
Optimization successful for k1=-0.3253. Start fval -8.901782, end fval -8.901783.
Next guess for k1 in direction -1 is -0.3352. Step size: -0.0099.
Optimization successful for k1=-0.3352. Start fval -8.875792, end fval -8.875793.
Next guess for k1 in direction -1 is -0.3473. Step size: -0.0121.
Optimization successful for k1=-0.3473. Start fval -8.836807, end fval -8.836809.
Next guess for k1 in direction -1 is -0.3620. Step size: -0.0147.
Optimization successful for k1=-0.3620. Start fval -8.778364, end fval -8.778371.
Next guess for k1 in direction -1 is -0.3798. Step size: -0.0178.
Optimization successful for k1=-0.3798. Start fval -8.690759, end fval -8.690777.
Next guess for k1 in direction -1 is -0.4014. Step size: -0.0216.
Optimization successful for k1=-0.4014. Start fval -8.559444, end fval -8.559492.
Next guess for k1 in direction -1 is -0.4277. Step size: -0.0262.
Optimization successful for k1=-0.4277. Start fval -8.362612, end fval -8.362747.
Next guess for k1 in direction -1 is -0.4595. Step size: -0.0318.
Optimization successful for k1=-0.4595. Start fval -8.067584, end fval -8.067954.
Next guess for k1 in direction -1 is -0.4980. Step size: -0.0385.
Optimization successful for k1=-0.4980. Start fval -7.625392, end fval -7.626571.
Next guess for k1 in direction -1 is -0.5444. Step size: -0.0464.
Optimization successful for k1=-0.5444. Start fval -6.963036, end fval -6.970023.
Next guess for k1 in direction -1 is -0.5986. Step size: -0.0542.
Optimization successful for k1=-0.5986. Start fval -5.981285, end fval -6.029407.
Next guess for k1 in direction -1 is -0.6610. Step size: -0.0624.
Optimization successful for k1=-0.6610. Start fval -4.571924, end fval -4.717577.
Next guess for k1 in direction -1 is -0.7259. Step size: -0.0649.
Optimization successful for k1=-0.7259. Start fval -2.704482, end fval -2.704482.
Next guess for k1 in direction 1 is -0.2258. Step size: 0.0548.
Optimization successful for k1=-0.2258. Start fval -8.880002, end fval -8.880425.
Next guess for k1 in direction 1 is -0.2131. Step size: 0.0126.
Optimization successful for k1=-0.2131. Start fval -8.843769, end fval -8.843781.
Next guess for k1 in direction 1 is -0.1975. Step size: 0.0156.
Optimization successful for k1=-0.1975. Start fval -8.788806, end fval -8.788807.
Next guess for k1 in direction 1 is -0.1781. Step size: 0.0194.
Optimization successful for k1=-0.1781. Start fval -8.706349, end fval -8.706353.
Next guess for k1 in direction 1 is -0.1538. Step size: 0.0243.
Optimization successful for k1=-0.1538. Start fval -8.582678, end fval -8.582685.
Next guess for k1 in direction 1 is -0.1233. Step size: 0.0305.
Optimization successful for k1=-0.1233. Start fval -8.397211, end fval -8.397224.
Next guess for k1 in direction 1 is -0.0845. Step size: 0.0388.
Optimization successful for k1=-0.0845. Start fval -8.119099, end fval -8.119121.
Next guess for k1 in direction 1 is -0.0342. Step size: 0.0503.
Optimization successful for k1=-0.0342. Start fval -7.702046, end fval -7.702079.
Next guess for k1 in direction 1 is 0.0334. Step size: 0.0676.
Optimization successful for k1=0.0334. Start fval -7.076370, end fval -7.076425.
Next guess for k1 in direction 1 is 0.1334. Step size: 0.1000.
Optimization successful for k1=0.1334. Start fval -6.094134, end fval -6.119372.
Next guess for k1 in direction 1 is 0.2334. Step size: 0.1000.
Optimization successful for k1=0.2334. Start fval -5.244782, end fval -5.266632.
Next guess for k1 in direction 1 is 0.3334. Step size: 0.1000.
Optimization successful for k1=0.3334. Start fval -4.610359, end fval -4.625943.
Next guess for k1 in direction 1 is 0.4334. Step size: 0.1000.
Optimization successful for k1=0.4334. Start fval -4.216160, end fval -4.224484.
Next guess for k2 in direction -1 is -0.8494. Step size: -0.0448.
Optimization successful for k2=-0.8494. Start fval -8.902793, end fval -8.937852.
Next guess for k2 in direction -1 is -0.8597. Step size: -0.0103.
Optimization successful for k2=-0.8597. Start fval -8.929892, end fval -8.929895.
Next guess for k2 in direction -1 is -0.8724. Step size: -0.0127.
Optimization successful for k2=-0.8724. Start fval -8.917959, end fval -8.917960.
Next guess for k2 in direction -1 is -0.8881. Step size: -0.0157.
Optimization successful for k2=-0.8881. Start fval -8.900059, end fval -8.900060.
Next guess for k2 in direction -1 is -0.9077. Step size: -0.0196.
Optimization successful for k2=-0.9077. Start fval -8.873212, end fval -8.873216.
Next guess for k2 in direction -1 is -0.9322. Step size: -0.0245.
Optimization successful for k2=-0.9322. Start fval -8.832952, end fval -8.832961.
Next guess for k2 in direction -1 is -0.9632. Step size: -0.0309.
Optimization successful for k2=-0.9632. Start fval -8.772584, end fval -8.772606.
Next guess for k2 in direction -1 is -1.0025. Step size: -0.0394.
Optimization successful for k2=-1.0025. Start fval -8.682080, end fval -8.682127.
Next guess for k2 in direction -1 is -1.0534. Step size: -0.0509.
Optimization successful for k2=-1.0534. Start fval -8.546401, end fval -8.546525.
Next guess for k2 in direction -1 is -1.1209. Step size: -0.0674.
Optimization successful for k2=-1.1209. Start fval -8.342997, end fval -8.343338.
Next guess for k2 in direction -1 is -1.2209. Step size: -0.1000.
Optimization successful for k2=-1.2209. Start fval -8.013581, end fval -8.014852.
Next guess for k2 in direction -1 is -1.3209. Step size: -0.1000.
Optimization successful for k2=-1.3209. Start fval -7.678982, end fval -7.680471.
Next guess for k2 in direction -1 is -1.4209. Step size: -0.1000.
Optimization successful for k2=-1.4209. Start fval -7.360217, end fval -7.361384.
Next guess for k2 in direction -1 is -1.5209. Step size: -0.1000.
Optimization successful for k2=-1.5209. Start fval -7.069435, end fval -7.070324.
Next guess for k2 in direction -1 is -1.6209. Step size: -0.1000.
Optimization successful for k2=-1.6209. Start fval -6.812727, end fval -6.813393.
Next guess for k2 in direction -1 is -1.7209. Step size: -0.1000.
Optimization successful for k2=-1.7209. Start fval -6.591549, end fval -6.592046.
Next guess for k2 in direction -1 is -1.8209. Step size: -0.1000.
Optimization successful for k2=-1.8209. Start fval -6.404485, end fval -6.404823.
Next guess for k2 in direction -1 is -1.9209. Step size: -0.1000.
Optimization successful for k2=-1.9209. Start fval -6.248413, end fval -6.248660.
Next guess for k2 in direction -1 is -2.0000. Step size: -0.0791.
Optimization successful for k2=-2.0000. Start fval -6.144116, end fval -6.144637.
Next guess for k2 in direction 1 is -0.7047. Step size: 0.1000.
Optimization successful for k2=-0.7047. Start fval -8.864364, end fval -8.866058.
Next guess for k2 in direction 1 is -0.6829. Step size: 0.0218.
Optimization successful for k2=-0.6829. Start fval -8.822251, end fval -8.822311.
Next guess for k2 in direction 1 is -0.6564. Step size: 0.0265.
Optimization successful for k2=-0.6564. Start fval -8.756661, end fval -8.756671.
Next guess for k2 in direction 1 is -0.6242. Step size: 0.0322.
Optimization successful for k2=-0.6242. Start fval -8.658235, end fval -8.658272.
Next guess for k2 in direction 1 is -0.5850. Step size: 0.0392.
Optimization successful for k2=-0.5850. Start fval -8.510658, end fval -8.510732.
Next guess for k2 in direction 1 is -0.5371. Step size: 0.0480.
Optimization successful for k2=-0.5371. Start fval -8.289342, end fval -8.289510.
Next guess for k2 in direction 1 is -0.4778. Step size: 0.0593.
Optimization successful for k2=-0.4778. Start fval -7.957494, end fval -7.957887.
Next guess for k2 in direction 1 is -0.4033. Step size: 0.0745.
Optimization successful for k2=-0.4033. Start fval -7.460130, end fval -7.461154.
Next guess for k2 in direction 1 is -0.3033. Step size: 0.1000.
Optimization successful for k2=-0.3033. Start fval -6.694137, end fval -6.697487.
Next guess for k2 in direction 1 is -0.2033. Step size: 0.1000.
Optimization successful for k2=-0.2033. Start fval -5.895272, end fval -5.899927.
Next guess for k2 in direction 1 is -0.1033. Step size: 0.1000.
Optimization successful for k2=-0.1033. Start fval -5.163258, end fval -5.167716.
Next guess for k2 in direction 1 is -0.0033. Step size: 0.1000.
Optimization successful for k2=-0.0033. Start fval -4.587096, end fval -4.590256.
Next guess for k2 in direction 1 is 0.0967. Step size: 0.1000.
Optimization successful for k2=0.0967. Start fval -4.210830, end fval -4.212306.

# specify the parameters, for which profiles should be computed
ax = visualize.profiles(result)

Sampling

import pypesto.sample as sample

sampler = sample.AdaptiveParallelTemperingSampler(
    internal_sampler=sample.AdaptiveMetropolisSampler(), n_chains=3
)

result = sample.sample(
    problem, n_samples=1000, sampler=sampler, result=result, filename=None
)

Initializing betas with "near-exponential decay".
Elapsed time: 2.0482147180000005

ax = visualize.sampling_scatter(result, size=[13, 6])

/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/v0.5.5/lib/python3.11/site-packages/pypesto/visualize/sampling.py:1223: 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, _ = get_data_to_plot(

Predict

# Let's create a function, which predicts the ratio of x_1 and x_0
import pypesto.predict as predict


def ratio_function(amici_output_list):
    # This (optional) function post-processes the results from AMICI and must accept one input:
    # a list of dicts, with the fields t, x, y[, sx, sy - if sensi_orders includes 1]
    # We need to specify the simulation condition: here, we only have one, i.e., it's [0]
    x = amici_output_list[0]["x"]
    ratio = x[:, 1] / x[:, 0]
    # we need to output also at least one simulation condition
    return [ratio]


# create pypesto prediction function
predictor = predict.AmiciPredictor(
    objective, post_processor=ratio_function, output_ids=["ratio"]
)

# run prediction
prediction = predictor(x=model.getUnscaledParameters())

dict(prediction)

{'conditions': [{'timepoints': array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.]),
   'output_ids': ['ratio'],
   'x_names': ['k1', 'k2'],
   'output': array([0.        , 1.95196396, 2.00246152, 2.00290412, 2.00290796,
          2.00290801, 2.00290801, 2.00290799, 2.002908  , 2.00290801,
          2.002908  ]),
   'output_sensi': None,
   'output_weight': None,
   'output_sigmay': None}],
 'condition_ids': ['condition_0'],
 'comment': None,
 'parameter_ids': ['k1', 'k2']}

Analyze parameter ensembles

# We want to use the sample result to create a prediction
from pypesto.ensemble import ensemble

# first collect some vectors from the sampling result
vectors = result.sample_result.trace_x[0, ::20, :]

# create the collection
my_ensemble = ensemble.Ensemble(
    vectors,
    x_names=problem.x_names,
    ensemble_type=ensemble.EnsembleType.sample,
    lower_bound=problem.lb,
    upper_bound=problem.ub,
)

# we can also perform an approximative identifiability analysis
summary = my_ensemble.compute_summary()
identifiability = my_ensemble.check_identifiability()
print(identifiability.transpose())

parameterId               k1        k2
parameterId               k1        k2
lowerBound                -2        -2
upperBound                 2         2
ensemble_mean      -0.542604  -0.50028
ensemble_std        0.262051  0.190921
ensemble_median    -0.542604  -0.50028
within lb: 1 std        True      True
within ub: 1 std        True      True
within lb: 2 std        True      True
within ub: 2 std        True      True
within lb: 3 std        True      True
within ub: 3 std        True      True
within lb: perc 5       True      True
within lb: perc 20      True      True
within ub: perc 80      True      True
within ub: perc 95      True      True

# run a prediction
ensemble_prediction = my_ensemble.predict(
    predictor, prediction_id="ratio_over_time"
)

# go for some analysis
prediction_summary = ensemble_prediction.compute_summary(
    percentiles_list=(1, 5, 10, 25, 75, 90, 95, 99)
)
dict(prediction_summary)

{'mean': <pypesto.result.predict.PredictionResult at 0x7f7be40c5b50>,
 'std': <pypesto.result.predict.PredictionResult at 0x7f7be40c6d50>,
 'median': <pypesto.result.predict.PredictionResult at 0x7f7be40c4f50>,
 'percentile 1': <pypesto.result.predict.PredictionResult at 0x7f7be4166f50>,
 'percentile 5': <pypesto.result.predict.PredictionResult at 0x7f7be4166dd0>,
 'percentile 10': <pypesto.result.predict.PredictionResult at 0x7f7be40b5c90>,
 'percentile 25': <pypesto.result.predict.PredictionResult at 0x7f7be40b6010>,
 'percentile 75': <pypesto.result.predict.PredictionResult at 0x7f7be40b40d0>,
 'percentile 90': <pypesto.result.predict.PredictionResult at 0x7f7be40b5e90>,
 'percentile 95': <pypesto.result.predict.PredictionResult at 0x7f7be40b7010>,
 'percentile 99': <pypesto.result.predict.PredictionResult at 0x7f7be40b5110>}