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 355 ms, sys: 3.77 ms, total: 359 ms
Wall time: 358 ms

Visualize

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

<seaborn.axisgrid.PairGrid at 0x7f70bba4ffd0>

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.2434. Step size: -0.0243.
Optimization successful for k1=-0.2434. Start fval -7.006209, end fval -7.044481.
Next guess for k1 in direction -1 is -0.2488. Step size: -0.0054.
Optimization successful for k1=-0.2488. Start fval -7.038100, end fval -7.038100.
Next guess for k1 in direction -1 is -0.2554. Step size: -0.0066.
Optimization successful for k1=-0.2554. Start fval -7.028530, end fval -7.028532.
Next guess for k1 in direction -1 is -0.2634. Step size: -0.0081.
Optimization successful for k1=-0.2634. Start fval -7.014178, end fval -7.014178.
Next guess for k1 in direction -1 is -0.2733. Step size: -0.0098.
Optimization successful for k1=-0.2733. Start fval -6.992641, end fval -6.992642.
Next guess for k1 in direction -1 is -0.2853. Step size: -0.0120.
Optimization successful for k1=-0.2853. Start fval -6.960344, end fval -6.960347.
Next guess for k1 in direction -1 is -0.2999. Step size: -0.0146.
Optimization successful for k1=-0.2999. Start fval -6.911928, end fval -6.911929.
Next guess for k1 in direction -1 is -0.3177. Step size: -0.0178.
Optimization successful for k1=-0.3177. Start fval -6.839343, end fval -6.839348.
Next guess for k1 in direction -1 is -0.3393. Step size: -0.0216.
Optimization successful for k1=-0.3393. Start fval -6.730532, end fval -6.730546.
Next guess for k1 in direction -1 is -0.3656. Step size: -0.0263.
Optimization successful for k1=-0.3656. Start fval -6.567404, end fval -6.567443.
Next guess for k1 in direction -1 is -0.3976. Step size: -0.0320.
Optimization successful for k1=-0.3976. Start fval -6.322828, end fval -6.322946.
Next guess for k1 in direction -1 is -0.4365. Step size: -0.0389.
Optimization successful for k1=-0.4365. Start fval -5.956124, end fval -5.956465.
Next guess for k1 in direction -1 is -0.4840. Step size: -0.0474.
Optimization successful for k1=-0.4840. Start fval -5.406319, end fval -5.407498.
Next guess for k1 in direction -1 is -0.5417. Step size: -0.0578.
Optimization successful for k1=-0.5417. Start fval -4.582765, end fval -4.587306.
Next guess for k1 in direction -1 is -0.6098. Step size: -0.0680.
Optimization successful for k1=-0.6098. Start fval -3.356962, end fval -3.412292.
Next guess for k1 in direction -1 is -0.6877. Step size: -0.0780.
Optimization successful for k1=-0.6877. Start fval -1.605543, end fval -1.799530.
Next guess for k1 in direction 1 is -0.1606. Step size: 0.0584.
Optimization successful for k1=-0.1606. Start fval -6.986863, end fval -6.986983.
Next guess for k1 in direction 1 is -0.1471. Step size: 0.0135.
Optimization successful for k1=-0.1471. Start fval -6.951887, end fval -6.951889.
Next guess for k1 in direction 1 is -0.1303. Step size: 0.0167.
Optimization successful for k1=-0.1303. Start fval -6.899263, end fval -6.899264.
Next guess for k1 in direction 1 is -0.1094. Step size: 0.0209.
Optimization successful for k1=-0.1094. Start fval -6.820351, end fval -6.820352.
Next guess for k1 in direction 1 is -0.0832. Step size: 0.0262.
Optimization successful for k1=-0.0832. Start fval -6.702017, end fval -6.702018.
Next guess for k1 in direction 1 is -0.0499. Step size: 0.0333.
Optimization successful for k1=-0.0499. Start fval -6.524556, end fval -6.524558.
Next guess for k1 in direction 1 is -0.0069. Step size: 0.0430.
Optimization successful for k1=-0.0069. Start fval -6.258374, end fval -6.258376.
Next guess for k1 in direction 1 is 0.0505. Step size: 0.0574.
Optimization successful for k1=0.0505. Start fval -5.858964, end fval -5.858966.
Next guess for k1 in direction 1 is 0.1322. Step size: 0.0818.
Optimization successful for k1=0.1322. Start fval -5.260216, end fval -5.260780.
Next guess for k1 in direction 1 is 0.2322. Step size: 0.1000.
Optimization successful for k1=0.2322. Start fval -4.576879, end fval -4.592235.
Next guess for k1 in direction 1 is 0.3322. Step size: 0.1000.
Optimization successful for k1=0.3322. Start fval -4.078678, end fval -4.089897.
Next guess for k1 in direction 1 is 0.4322. Step size: 0.1000.
Optimization successful for k1=0.4322. Start fval -3.775576, end fval -3.781408.
Next guess for k1 in direction 1 is 0.5322. Step size: 0.1000.
Optimization successful for k1=0.5322. Start fval -3.627165, end fval -3.629165.
Next guess for k1 in direction 1 is 0.6322. Step size: 0.1000.
Optimization successful for k1=0.6322. Start fval -3.570244, end fval -3.570664.
Next guess for k1 in direction 1 is 0.7322. Step size: 0.1000.
Optimization successful for k1=0.7322. Start fval -3.553962, end fval -3.554011.
Next guess for k1 in direction 1 is 0.8322. Step size: 0.1000.
Optimization successful for k1=0.8322. Start fval -3.550733, end fval -3.550736.
Next guess for k1 in direction 1 is 0.9322. Step size: 0.1000.
Optimization successful for k1=0.9322. Start fval -3.550329, end fval -3.550329.
Next guess for k1 in direction 1 is 1.0322. Step size: 0.1000.
Optimization successful for k1=1.0322. Start fval -3.550301, end fval -3.550301.
Next guess for k1 in direction 1 is 1.1322. Step size: 0.1000.
Optimization successful for k1=1.1322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.2322. Step size: 0.1000.
Optimization successful for k1=1.2322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.3322. Step size: 0.1000.
Optimization successful for k1=1.3322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.4322. Step size: 0.1000.
Optimization successful for k1=1.4322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.5322. Step size: 0.1000.
Optimization successful for k1=1.5322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.6322. Step size: 0.1000.
Optimization successful for k1=1.6322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.7322. Step size: 0.1000.
Optimization successful for k1=1.7322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.8322. Step size: 0.1000.
Optimization successful for k1=1.8322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 1.9322. Step size: 0.1000.
Optimization successful for k1=1.9322. Start fval -3.550300, end fval -3.550300.
Next guess for k1 in direction 1 is 2.0000. Step size: 0.0678.
Optimization successful for k1=2.0000. Start fval -3.477989, end fval -3.550300.
Next guess for k2 in direction -1 is -0.6436. Step size: -0.0359.
Optimization successful for k2=-0.6436. Start fval -7.006270, end fval -7.044138.
Next guess for k2 in direction -1 is -0.6517. Step size: -0.0081.
Optimization successful for k2=-0.6517. Start fval -7.037583, end fval -7.037584.
Next guess for k2 in direction -1 is -0.6617. Step size: -0.0100.
Optimization successful for k2=-0.6617. Start fval -7.027752, end fval -7.027752.
Next guess for k2 in direction -1 is -0.6740. Step size: -0.0123.
Optimization successful for k2=-0.6740. Start fval -7.013000, end fval -7.013003.
Next guess for k2 in direction -1 is -0.6891. Step size: -0.0151.
Optimization successful for k2=-0.6891. Start fval -6.990879, end fval -6.990879.
Next guess for k2 in direction -1 is -0.7079. Step size: -0.0187.
Optimization successful for k2=-0.7079. Start fval -6.957700, end fval -6.957703.
Next guess for k2 in direction -1 is -0.7311. Step size: -0.0232.
Optimization successful for k2=-0.7311. Start fval -6.907947, end fval -6.907953.
Next guess for k2 in direction -1 is -0.7600. Step size: -0.0289.
Optimization successful for k2=-0.7600. Start fval -6.833334, end fval -6.833349.
Next guess for k2 in direction -1 is -0.7964. Step size: -0.0363.
Optimization successful for k2=-0.7964. Start fval -6.721436, end fval -6.721470.
Next guess for k2 in direction -1 is -0.8425. Step size: -0.0461.
Optimization successful for k2=-0.8425. Start fval -6.553598, end fval -6.553695.
Next guess for k2 in direction -1 is -0.9020. Step size: -0.0595.
Optimization successful for k2=-0.9020. Start fval -6.301988, end fval -6.302214.
Next guess for k2 in direction -1 is -0.9809. Step size: -0.0789.
Optimization successful for k2=-0.9809. Start fval -5.924880, end fval -5.925555.
Next guess for k2 in direction -1 is -1.0809. Step size: -0.1000.
Optimization successful for k2=-1.0809. Start fval -5.409407, end fval -5.411145.
Next guess for k2 in direction -1 is -1.1809. Step size: -0.1000.
Optimization successful for k2=-1.1809. Start fval -4.891370, end fval -4.893398.
Next guess for k2 in direction -1 is -1.2809. Step size: -0.1000.
Optimization successful for k2=-1.2809. Start fval -4.400330, end fval -4.402119.
Next guess for k2 in direction -1 is -1.3809. Step size: -0.1000.
Optimization successful for k2=-1.3809. Start fval -3.954077, end fval -3.955616.
Next guess for k2 in direction -1 is -1.4809. Step size: -0.1000.
Optimization successful for k2=-1.4809. Start fval -3.561307, end fval -3.562616.
Next guess for k2 in direction -1 is -1.5809. Step size: -0.1000.
Optimization successful for k2=-1.5809. Start fval -3.223931, end fval -3.224935.
Next guess for k2 in direction -1 is -1.6809. Step size: -0.1000.
Optimization successful for k2=-1.6809. Start fval -2.939273, end fval -2.940008.
Next guess for k2 in direction -1 is -1.7809. Step size: -0.1000.
Optimization successful for k2=-1.7809. Start fval -2.702330, end fval -2.702877.
Next guess for k2 in direction -1 is -1.8809. Step size: -0.1000.
Optimization successful for k2=-1.8809. Start fval -2.507227, end fval -2.507584.
Next guess for k2 in direction -1 is -1.9809. Step size: -0.1000.
Optimization successful for k2=-1.9809. Start fval -2.347764, end fval -2.348030.
Next guess for k2 in direction 1 is -0.5265. Step size: 0.0812.
Optimization successful for k2=-0.5265. Start fval -6.987192, end fval -6.987810.
Next guess for k2 in direction 1 is -0.5084. Step size: 0.0182.
Optimization successful for k2=-0.5084. Start fval -6.953114, end fval -6.953138.
Next guess for k2 in direction 1 is -0.4861. Step size: 0.0223.
Optimization successful for k2=-0.4861. Start fval -6.901130, end fval -6.901134.
Next guess for k2 in direction 1 is -0.4587. Step size: 0.0274.
Optimization successful for k2=-0.4587. Start fval -6.823157, end fval -6.823172.
Next guess for k2 in direction 1 is -0.4249. Step size: 0.0338.
Optimization successful for k2=-0.4249. Start fval -6.706261, end fval -6.706292.
Next guess for k2 in direction 1 is -0.3828. Step size: 0.0421.
Optimization successful for k2=-0.3828. Start fval -6.530994, end fval -6.531084.
Next guess for k2 in direction 1 is -0.3298. Step size: 0.0530.
Optimization successful for k2=-0.3298. Start fval -6.268230, end fval -6.268459.
Next guess for k2 in direction 1 is -0.2611. Step size: 0.0686.
Optimization successful for k2=-0.2611. Start fval -5.874248, end fval -5.874922.
Next guess for k2 in direction 1 is -0.1611. Step size: 0.1000.
Optimization successful for k2=-0.1611. Start fval -5.244008, end fval -5.246875.
Next guess for k2 in direction 1 is -0.0611. Step size: 0.1000.
Optimization successful for k2=-0.0611. Start fval -4.638262, end fval -4.642173.
Next guess for k2 in direction 1 is 0.0389. Step size: 0.1000.
Optimization successful for k2=0.0389. Start fval -4.148412, end fval -4.151227.
Next guess for k2 in direction 1 is 0.1389. Step size: 0.1000.
Optimization successful for k2=0.1389. Start fval -3.822383, end fval -3.823444.
Next guess for k2 in direction 1 is 0.2389. Step size: 0.1000.
Optimization successful for k2=0.2389. Start fval -3.648159, end fval -3.649014.
Next guess for k2 in direction 1 is 0.3389. Step size: 0.1000.
Optimization successful for k2=0.3389. Start fval -3.576994, end fval -3.577401.
Next guess for k2 in direction 1 is 0.4389. Step size: 0.1000.
Optimization successful for k2=0.4389. Start fval -3.555533, end fval -3.555614.
Next guess for k2 in direction 1 is 0.5389. Step size: 0.1000.
Optimization successful for k2=0.5389. Start fval -3.550977, end fval -3.550984.
Next guess for k2 in direction 1 is 0.6389. Step size: 0.1000.
Optimization successful for k2=0.6389. Start fval -3.550352, end fval -3.550352.
Next guess for k2 in direction 1 is 0.7389. Step size: 0.1000.
Optimization successful for k2=0.7389. Start fval -3.550302, end fval -3.550302.
Next guess for k2 in direction 1 is 0.8389. Step size: 0.1000.
Optimization successful for k2=0.8389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 0.9389. Step size: 0.1000.
Optimization successful for k2=0.9389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.0389. Step size: 0.1000.
Optimization successful for k2=1.0389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.1389. Step size: 0.1000.
Optimization successful for k2=1.1389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.2389. Step size: 0.1000.
Optimization successful for k2=1.2389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.3389. Step size: 0.1000.
Optimization successful for k2=1.3389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.4389. Step size: 0.1000.
Optimization successful for k2=1.4389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.5389. Step size: 0.1000.
Optimization successful for k2=1.5389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.6389. Step size: 0.1000.
Optimization successful for k2=1.6389. Start fval -3.550300, end fval -3.550300.
Next guess for k2 in direction 1 is 1.7389. Step size: 0.1000.
Optimization successful for k2=1.7389. Start fval -3.538701, end fval -3.538701.
Next guess for k2 in direction 1 is 1.7597. Step size: 0.0209.
Optimization successful for k2=1.7597. Start fval -3.470398, end fval -3.470398.
Next guess for k2 in direction 1 is 1.7833. Step size: 0.0236.
Optimization successful for k2=1.7833. Start fval -3.316731, end fval -3.316731.
Next guess for k2 in direction 1 is 1.8105. Step size: 0.0272.
Optimization successful for k2=1.8105. Start fval -3.035626, end fval -3.035626.
Next guess for k2 in direction 1 is 1.8421. Step size: 0.0316.
Optimization successful for k2=1.8421. Start fval -2.563486, end fval -2.563486.
Next guess for k2 in direction 1 is 1.8790. Step size: 0.0370.
Optimization successful for k2=1.8790. Start fval -1.804144, end fval -1.804144.

# 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: 1.8660506609999992

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

/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/v0.5.4/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.413377 -0.460881
ensemble_std         0.19437  0.156101
ensemble_median    -0.413377 -0.460881
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 0x7f70b8a5e6d0>,
 'std': <pypesto.result.predict.PredictionResult at 0x7f70b8a5d950>,
 'median': <pypesto.result.predict.PredictionResult at 0x7f70b8a5d7d0>,
 'percentile 1': <pypesto.result.predict.PredictionResult at 0x7f70b8a5d490>,
 'percentile 5': <pypesto.result.predict.PredictionResult at 0x7f70b8a5f310>,
 'percentile 10': <pypesto.result.predict.PredictionResult at 0x7f70b8a5e150>,
 'percentile 25': <pypesto.result.predict.PredictionResult at 0x7f70b8a5c5d0>,
 'percentile 75': <pypesto.result.predict.PredictionResult at 0x7f70b8a5c690>,
 'percentile 90': <pypesto.result.predict.PredictionResult at 0x7f70b8a5d4d0>,
 'percentile 95': <pypesto.result.predict.PredictionResult at 0x7f70b8a5dbd0>,
 'percentile 99': <pypesto.result.predict.PredictionResult at 0x7f70b8a5ead0>}