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 204 ms, sys: 4.24 ms, total: 208 ms
Wall time: 211 ms

Visualize

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

<seaborn.axisgrid.PairGrid at 0x752a3313f410>

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.5574. Step size: -0.0219.
Optimization successful for k1=-0.5574. Start fval -7.769848, end fval -7.807714.
Next guess for k1 in direction -1 is -0.5623. Step size: -0.0049.
Optimization successful for k1=-0.5623. Start fval -7.801301, end fval -7.801328.
Next guess for k1 in direction -1 is -0.5682. Step size: -0.0059.
Optimization successful for k1=-0.5682. Start fval -7.791723, end fval -7.791733.
Next guess for k1 in direction -1 is -0.5755. Step size: -0.0072.
Optimization successful for k1=-0.5755. Start fval -7.777342, end fval -7.777366.
Next guess for k1 in direction -1 is -0.5843. Step size: -0.0088.
Optimization successful for k1=-0.5843. Start fval -7.755805, end fval -7.755870.
Next guess for k1 in direction -1 is -0.5950. Step size: -0.0107.
Optimization successful for k1=-0.5950. Start fval -7.723579, end fval -7.723756.
Next guess for k1 in direction -1 is -0.6080. Step size: -0.0129.
Optimization successful for k1=-0.6080. Start fval -7.675533, end fval -7.676221.
Next guess for k1 in direction -1 is -0.6230. Step size: -0.0150.
Optimization successful for k1=-0.6230. Start fval -7.604138, end fval -7.604386.
Next guess for k1 in direction -1 is -0.6368. Step size: -0.0139.
Optimization successful for k1=-0.6368. Start fval -7.496577, end fval -7.496577.
Next guess for k1 in direction -1 is -0.6514. Step size: -0.0145.
Optimization successful for k1=-0.6514. Start fval -7.334793, end fval -7.334792.
Next guess for k1 in direction -1 is -0.6677. Step size: -0.0163.
Optimization successful for k1=-0.6677. Start fval -7.092278, end fval -7.092278.
Next guess for k1 in direction -1 is -0.6867. Step size: -0.0190.
Optimization successful for k1=-0.6867. Start fval -6.728188, end fval -6.728188.
Next guess for k1 in direction -1 is -0.7091. Step size: -0.0224.
Optimization successful for k1=-0.7091. Start fval -6.182617, end fval -6.182617.
Next guess for k1 in direction -1 is -0.7358. Step size: -0.0267.
Optimization successful for k1=-0.7358. Start fval -5.364786, end fval -5.364786.
Next guess for k1 in direction -1 is -0.7680. Step size: -0.0322.
Optimization successful for k1=-0.7680. Start fval -4.138435, end fval -4.138435.
Next guess for k1 in direction -1 is -0.8071. Step size: -0.0390.
Optimization successful for k1=-0.8071. Start fval -2.299035, end fval -2.299035.
Next guess for k1 in direction 1 is -0.4925. Step size: 0.0429.
Optimization successful for k1=-0.4925. Start fval -7.752011, end fval -7.772847.
Next guess for k1 in direction 1 is -0.4828. Step size: 0.0097.
Optimization successful for k1=-0.4828. Start fval -7.748999, end fval -7.749252.
Next guess for k1 in direction 1 is -0.4707. Step size: 0.0121.
Optimization successful for k1=-0.4707. Start fval -7.713611, end fval -7.713651.
Next guess for k1 in direction 1 is -0.4557. Step size: 0.0150.
Optimization successful for k1=-0.4557. Start fval -7.660221, end fval -7.660300.
Next guess for k1 in direction 1 is -0.4372. Step size: 0.0186.
Optimization successful for k1=-0.4372. Start fval -7.580216, end fval -7.580376.
Next guess for k1 in direction 1 is -0.4141. Step size: 0.0231.
Optimization successful for k1=-0.4141. Start fval -7.460355, end fval -7.460628.
Next guess for k1 in direction 1 is -0.3851. Step size: 0.0289.
Optimization successful for k1=-0.3851. Start fval -7.280756, end fval -7.281266.
Next guess for k1 in direction 1 is -0.3486. Step size: 0.0365.
Optimization successful for k1=-0.3486. Start fval -7.011718, end fval -7.012621.
Next guess for k1 in direction 1 is -0.3021. Step size: 0.0465.
Optimization successful for k1=-0.3021. Start fval -6.608664, end fval -6.610213.
Next guess for k1 in direction 1 is -0.2419. Step size: 0.0602.
Optimization successful for k1=-0.2419. Start fval -6.007086, end fval -6.007187.
Next guess for k1 in direction 1 is -0.1632. Step size: 0.0787.
Optimization successful for k1=-0.1632. Start fval -5.106022, end fval -5.116744.
Next guess for k1 in direction 1 is -0.0632. Step size: 0.1000.
Optimization successful for k1=-0.0632. Start fval -3.835290, end fval -3.895238.
Next guess for k1 in direction 1 is 0.0368. Step size: 0.1000.
Optimization successful for k1=0.0368. Start fval -2.634729, end fval -2.679969.
Next guess for k2 in direction -1 is -1.5149. Step size: -0.1000.
Optimization successful for k2=-1.5149. Start fval -7.784484, end fval -7.811323.
Next guess for k2 in direction -1 is -1.5401. Step size: -0.0252.
Optimization successful for k2=-1.5401. Start fval -7.806715, end fval -7.806740.
Next guess for k2 in direction -1 is -1.5729. Step size: -0.0328.
Optimization successful for k2=-1.5729. Start fval -7.799842, end fval -7.799850.
Next guess for k2 in direction -1 is -1.6162. Step size: -0.0433.
Optimization successful for k2=-1.6162. Start fval -7.789507, end fval -7.789530.
Next guess for k2 in direction -1 is -1.6750. Step size: -0.0587.
Optimization successful for k2=-1.6750. Start fval -7.774023, end fval -7.774078.
Next guess for k2 in direction -1 is -1.7579. Step size: -0.0829.
Optimization successful for k2=-1.7579. Start fval -7.750850, end fval -7.751026.
Next guess for k2 in direction -1 is -1.8579. Step size: -0.1000.
Optimization successful for k2=-1.8579. Start fval -7.723214, end fval -7.723524.
Next guess for k2 in direction -1 is -1.9579. Step size: -0.1000.
Optimization successful for k2=-1.9579. Start fval -7.697767, end fval -7.698027.
Next guess for k2 in direction -1 is -2.0000. Step size: -0.0421.
Optimization successful for k2=-2.0000. Start fval -7.686342, end fval -7.688118.
Next guess for k2 in direction 1 is -1.3149. Step size: 0.1000.
Optimization successful for k2=-1.3149. Start fval -7.805759, end fval -7.807037.
Next guess for k2 in direction 1 is -1.2948. Step size: 0.0202.
Optimization successful for k2=-1.2948. Start fval -7.800291, end fval -7.800315.
Next guess for k2 in direction 1 is -1.2709. Step size: 0.0239.
Optimization successful for k2=-1.2709. Start fval -7.790199, end fval -7.790205.
Next guess for k2 in direction 1 is -1.2429. Step size: 0.0280.
Optimization successful for k2=-1.2429. Start fval -7.775054, end fval -7.775064.
Next guess for k2 in direction 1 is -1.2103. Step size: 0.0326.
Optimization successful for k2=-1.2103. Start fval -7.752362, end fval -7.752382.
Next guess for k2 in direction 1 is -1.1726. Step size: 0.0378.
Optimization successful for k2=-1.1726. Start fval -7.718316, end fval -7.718358.
Next guess for k2 in direction 1 is -1.1291. Step size: 0.0434.
Optimization successful for k2=-1.1291. Start fval -7.667340, end fval -7.667422.
Next guess for k2 in direction 1 is -1.0794. Step size: 0.0498.
Optimization successful for k2=-1.0794. Start fval -7.590987, end fval -7.591141.
Next guess for k2 in direction 1 is -1.0226. Step size: 0.0568.
Optimization successful for k2=-1.0226. Start fval -7.476516, end fval -7.476815.
Next guess for k2 in direction 1 is -0.9579. Step size: 0.0647.
Optimization successful for k2=-0.9579. Start fval -7.305284, end fval -7.305802.
Next guess for k2 in direction 1 is -0.8841. Step size: 0.0738.
Optimization successful for k2=-0.8841. Start fval -7.048566, end fval -7.049633.
Next guess for k2 in direction 1 is -0.7994. Step size: 0.0847.
Optimization successful for k2=-0.7994. Start fval -6.664775, end fval -6.666711.
Next guess for k2 in direction 1 is -0.6994. Step size: 0.1000.
Optimization successful for k2=-0.6994. Start fval -6.079838, end fval -6.083569.
Next guess for k2 in direction 1 is -0.5994. Step size: 0.1000.
Optimization successful for k2=-0.5994. Start fval -5.356462, end fval -5.360704.
Next guess for k2 in direction 1 is -0.4994. Step size: 0.1000.
Optimization successful for k2=-0.4994. Start fval -4.511640, end fval -4.515579.
Next guess for k2 in direction 1 is -0.3994. Step size: 0.1000.
Optimization successful for k2=-0.3994. Start fval -3.579743, end fval -3.583415.
Next guess for k2 in direction 1 is -0.2994. Step size: 0.1000.
Optimization successful for k2=-0.2994. Start fval -2.614729, end fval -2.618341.

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/v0.5.6/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.975215 -0.682108
ensemble_std         0.43972  0.310107
ensemble_median    -0.975215 -0.682108
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       False      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 0x752a30279190>,
 'std': <pypesto.result.predict.PredictionResult at 0x752a3027bc90>,
 'median': <pypesto.result.predict.PredictionResult at 0x752a30278b90>,
 'percentile 1': <pypesto.result.predict.PredictionResult at 0x752a302790d0>,
 'percentile 5': <pypesto.result.predict.PredictionResult at 0x752a302789d0>,
 'percentile 10': <pypesto.result.predict.PredictionResult at 0x752a3027bf50>,
 'percentile 25': <pypesto.result.predict.PredictionResult at 0x752a3027a710>,
 'percentile 75': <pypesto.result.predict.PredictionResult at 0x752a30164150>,
 'percentile 90': <pypesto.result.predict.PredictionResult at 0x752a30165550>,
 'percentile 95': <pypesto.result.predict.PredictionResult at 0x752a30164610>,
 'percentile 99': <pypesto.result.predict.PredictionResult at 0x752a30164090>}