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 190 ms, sys: 6.01 ms, total: 196 ms
Wall time: 197 ms

Visualize

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

<seaborn.axisgrid.PairGrid at 0x7f2f7484f450>

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.4000. Step size: -0.0233.
Optimization successful for k1=-0.4000. Start fval -13.358559, end fval -13.392485.
Next guess for k1 in direction -1 is -0.4052. Step size: -0.0052.
Optimization successful for k1=-0.4052. Start fval -13.383910, end fval -13.383912.
Next guess for k1 in direction -1 is -0.4116. Step size: -0.0064.
Optimization successful for k1=-0.4116. Start fval -13.371049, end fval -13.371050.
Next guess for k1 in direction -1 is -0.4193. Step size: -0.0078.
Optimization successful for k1=-0.4193. Start fval -13.351754, end fval -13.351755.
Next guess for k1 in direction -1 is -0.4288. Step size: -0.0095.
Optimization successful for k1=-0.4288. Start fval -13.322827, end fval -13.322829.
Next guess for k1 in direction -1 is -0.4403. Step size: -0.0115.
Optimization successful for k1=-0.4403. Start fval -13.279462, end fval -13.279468.
Next guess for k1 in direction -1 is -0.4543. Step size: -0.0140.
Optimization successful for k1=-0.4543. Start fval -13.214456, end fval -13.214470.
Next guess for k1 in direction -1 is -0.4714. Step size: -0.0171.
Optimization successful for k1=-0.4714. Start fval -13.117011, end fval -13.117050.
Next guess for k1 in direction -1 is -0.4921. Step size: -0.0208.
Optimization successful for k1=-0.4921. Start fval -12.970940, end fval -12.971040.
Next guess for k1 in direction -1 is -0.5174. Step size: -0.0252.
Optimization successful for k1=-0.5174. Start fval -12.751976, end fval -12.752258.
Next guess for k1 in direction -1 is -0.5481. Step size: -0.0307.
Optimization successful for k1=-0.5481. Start fval -12.423795, end fval -12.424665.
Next guess for k1 in direction -1 is -0.5854. Step size: -0.0373.
Optimization successful for k1=-0.5854. Start fval -11.932182, end fval -11.936138.
Next guess for k1 in direction -1 is -0.6293. Step size: -0.0439.
Optimization successful for k1=-0.6293. Start fval -11.201310, end fval -11.233449.
Next guess for k1 in direction -1 is -0.6802. Step size: -0.0509.
Optimization successful for k1=-0.6802. Start fval -10.148921, end fval -10.239896.
Next guess for k1 in direction -1 is -0.7333. Step size: -0.0531.
Optimization successful for k1=-0.7333. Start fval -8.702119, end fval -8.702119.
Next guess for k1 in direction 1 is -0.3530. Step size: 0.0237.
Optimization successful for k1=-0.3530. Start fval -13.358410, end fval -13.392422.
Next guess for k1 in direction 1 is -0.3476. Step size: 0.0054.
Optimization successful for k1=-0.3476. Start fval -13.383819, end fval -13.383821.
Next guess for k1 in direction 1 is -0.3410. Step size: 0.0066.
Optimization successful for k1=-0.3410. Start fval -13.370920, end fval -13.370920.
Next guess for k1 in direction 1 is -0.3329. Step size: 0.0081.
Optimization successful for k1=-0.3329. Start fval -13.351572, end fval -13.351572.
Next guess for k1 in direction 1 is -0.3228. Step size: 0.0100.
Optimization successful for k1=-0.3228. Start fval -13.322551, end fval -13.322553.
Next guess for k1 in direction 1 is -0.3104. Step size: 0.0124.
Optimization successful for k1=-0.3104. Start fval -13.279020, end fval -13.279022.
Next guess for k1 in direction 1 is -0.2951. Step size: 0.0154.
Optimization successful for k1=-0.2951. Start fval -13.213719, end fval -13.213723.
Next guess for k1 in direction 1 is -0.2760. Step size: 0.0191.
Optimization successful for k1=-0.2760. Start fval -13.115763, end fval -13.115772.
Next guess for k1 in direction 1 is -0.2521. Step size: 0.0238.
Optimization successful for k1=-0.2521. Start fval -12.968835, end fval -12.968852.
Next guess for k1 in direction 1 is -0.2221. Step size: 0.0300.
Optimization successful for k1=-0.2221. Start fval -12.748483, end fval -12.748514.
Next guess for k1 in direction 1 is -0.1840. Step size: 0.0382.
Optimization successful for k1=-0.1840. Start fval -12.418069, end fval -12.418124.
Next guess for k1 in direction 1 is -0.1346. Step size: 0.0494.
Optimization successful for k1=-0.1346. Start fval -11.922645, end fval -11.922746.
Next guess for k1 in direction 1 is -0.0687. Step size: 0.0659.
Optimization successful for k1=-0.0687. Start fval -11.179575, end fval -11.179739.
Next guess for k1 in direction 1 is 0.0313. Step size: 0.1000.
Optimization successful for k1=0.0313. Start fval -9.943826, end fval -9.983117.
Next guess for k1 in direction 1 is 0.1313. Step size: 0.1000.
Optimization successful for k1=0.1313. Start fval -8.821524, end fval -8.855477.
Next guess for k1 in direction 1 is 0.2313. Step size: 0.1000.
Optimization successful for k1=0.2313. Start fval -7.903020, end fval -7.929653.
Next guess for k2 in direction -1 is -0.9786. Step size: -0.0505.
Optimization successful for k2=-0.9786. Start fval -13.358861, end fval -13.392179.
Next guess for k2 in direction -1 is -0.9903. Step size: -0.0117.
Optimization successful for k2=-0.9903. Start fval -13.383450, end fval -13.383454.
Next guess for k2 in direction -1 is -1.0049. Step size: -0.0146.
Optimization successful for k2=-1.0049. Start fval -13.370366, end fval -13.370367.
Next guess for k2 in direction -1 is -1.0231. Step size: -0.0182.
Optimization successful for k2=-1.0231. Start fval -13.350740, end fval -13.350743.
Next guess for k2 in direction -1 is -1.0461. Step size: -0.0229.
Optimization successful for k2=-1.0461. Start fval -13.321307, end fval -13.321313.
Next guess for k2 in direction -1 is -1.0751. Step size: -0.0290.
Optimization successful for k2=-1.0751. Start fval -13.277162, end fval -13.277176.
Next guess for k2 in direction -1 is -1.1123. Step size: -0.0372.
Optimization successful for k2=-1.1123. Start fval -13.210942, end fval -13.210977.
Next guess for k2 in direction -1 is -1.1608. Step size: -0.0485.
Optimization successful for k2=-1.1608. Start fval -13.111686, end fval -13.111786.
Next guess for k2 in direction -1 is -1.2255. Step size: -0.0648.
Optimization successful for k2=-1.2255. Start fval -12.962932, end fval -12.963213.
Next guess for k2 in direction -1 is -1.3156. Step size: -0.0901.
Optimization successful for k2=-1.3156. Start fval -12.740002, end fval -12.740732.
Next guess for k2 in direction -1 is -1.4156. Step size: -0.1000.
Optimization successful for k2=-1.4156. Start fval -12.490456, end fval -12.491623.
Next guess for k2 in direction -1 is -1.5156. Step size: -0.1000.
Optimization successful for k2=-1.5156. Start fval -12.254937, end fval -12.255893.
Next guess for k2 in direction -1 is -1.6156. Step size: -0.1000.
Optimization successful for k2=-1.6156. Start fval -12.041954, end fval -12.042662.
Next guess for k2 in direction -1 is -1.7156. Step size: -0.1000.
Optimization successful for k2=-1.7156. Start fval -11.855257, end fval -11.855824.
Next guess for k2 in direction -1 is -1.8156. Step size: -0.1000.
Optimization successful for k2=-1.8156. Start fval -11.695516, end fval -11.695861.
Next guess for k2 in direction -1 is -1.9156. Step size: -0.1000.
Optimization successful for k2=-1.9156. Start fval -11.560994, end fval -11.561243.
Next guess for k2 in direction -1 is -2.0000. Step size: -0.0844.
Optimization successful for k2=-2.0000. Start fval -11.465076, end fval -11.465511.
Next guess for k2 in direction 1 is -0.8810. Step size: 0.0471.
Optimization successful for k2=-0.8810. Start fval -13.358505, end fval -13.392862.
Next guess for k2 in direction 1 is -0.8707. Step size: 0.0103.
Optimization successful for k2=-0.8707. Start fval -13.384478, end fval -13.384481.
Next guess for k2 in direction 1 is -0.8582. Step size: 0.0125.
Optimization successful for k2=-0.8582. Start fval -13.371906, end fval -13.371906.
Next guess for k2 in direction 1 is -0.8431. Step size: 0.0151.
Optimization successful for k2=-0.8431. Start fval -13.353039, end fval -13.353041.
Next guess for k2 in direction 1 is -0.8248. Step size: 0.0183.
Optimization successful for k2=-0.8248. Start fval -13.324760, end fval -13.324763.
Next guess for k2 in direction 1 is -0.8027. Step size: 0.0221.
Optimization successful for k2=-0.8027. Start fval -13.282372, end fval -13.282380.
Next guess for k2 in direction 1 is -0.7761. Step size: 0.0266.
Optimization successful for k2=-0.7761. Start fval -13.218835, end fval -13.218848.
Next guess for k2 in direction 1 is -0.7441. Step size: 0.0321.
Optimization successful for k2=-0.7441. Start fval -13.123570, end fval -13.123603.
Next guess for k2 in direction 1 is -0.7054. Step size: 0.0387.
Optimization successful for k2=-0.7054. Start fval -12.980699, end fval -12.980787.
Next guess for k2 in direction 1 is -0.6587. Step size: 0.0467.
Optimization successful for k2=-0.6587. Start fval -12.766378, end fval -12.766558.
Next guess for k2 in direction 1 is -0.6019. Step size: 0.0568.
Optimization successful for k2=-0.6019. Start fval -12.445163, end fval -12.445558.
Next guess for k2 in direction 1 is -0.5321. Step size: 0.0698.
Optimization successful for k2=-0.5321. Start fval -11.963831, end fval -11.964750.
Next guess for k2 in direction 1 is -0.4443. Step size: 0.0877.
Optimization successful for k2=-0.4443. Start fval -11.242640, end fval -11.244908.
Next guess for k2 in direction 1 is -0.3443. Step size: 0.1000.
Optimization successful for k2=-0.3443. Start fval -10.313189, end fval -10.317295.
Next guess for k2 in direction 1 is -0.2443. Step size: 0.1000.
Optimization successful for k2=-0.2443. Start fval -9.349437, end fval -9.353921.
Next guess for k2 in direction 1 is -0.1443. Step size: 0.1000.
Optimization successful for k2=-0.1443. Start fval -8.444544, end fval -8.448752.

# 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".
Initializing parallel chains with a combination of the starting point and prior samples with weight: 0.9.
Elapsed time: 1.0279305050000005

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

/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/v0.5.8/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.694004 -0.937787
ensemble_std        0.273818  0.416402
ensemble_median    -0.694004 -0.937787
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     False
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 0x7f2f6bfe7d50>,
 'std': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe7610>,
 'median': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe4f50>,
 'percentile 1': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe75d0>,
 'percentile 5': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe6890>,
 'percentile 10': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe57d0>,
 'percentile 25': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe5c90>,
 'percentile 75': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe4290>,
 'percentile 90': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe46d0>,
 'percentile 95': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe51d0>,
 'percentile 99': <pypesto.result.predict.PredictionResult at 0x7f2f6bfe4c90>}