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 183 ms, sys: 6.95 ms, total: 189 ms
Wall time: 190 ms

Visualize

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

<seaborn.axisgrid.PairGrid at 0x781e1b103710>

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.5201. Step size: -0.0221.
Optimization successful for k1=-0.5201. Start fval -4.663366, end fval -4.698129.
Next guess for k1 in direction -1 is -0.5250. Step size: -0.0049.
Optimization successful for k1=-0.5250. Start fval -4.690146, end fval -4.690156.
Next guess for k1 in direction -1 is -0.5310. Step size: -0.0060.
Optimization successful for k1=-0.5310. Start fval -4.678189, end fval -4.678192.
Next guess for k1 in direction -1 is -0.5383. Step size: -0.0073.
Optimization successful for k1=-0.5383. Start fval -4.660257, end fval -4.660264.
Next guess for k1 in direction -1 is -0.5472. Step size: -0.0089.
Optimization successful for k1=-0.5472. Start fval -4.633380, end fval -4.633397.
Next guess for k1 in direction -1 is -0.5580. Step size: -0.0108.
Optimization successful for k1=-0.5580. Start fval -4.593100, end fval -4.593144.
Next guess for k1 in direction -1 is -0.5712. Step size: -0.0132.
Optimization successful for k1=-0.5712. Start fval -4.532744, end fval -4.532860.
Next guess for k1 in direction -1 is -0.5871. Step size: -0.0160.
Optimization successful for k1=-0.5871. Start fval -4.442332, end fval -4.442657.
Next guess for k1 in direction -1 is -0.6064. Step size: -0.0193.
Optimization successful for k1=-0.6064. Start fval -4.307046, end fval -4.308398.
Next guess for k1 in direction -1 is -0.6290. Step size: -0.0225.
Optimization successful for k1=-0.6290. Start fval -4.106070, end fval -4.115722.
Next guess for k1 in direction -1 is -0.6539. Step size: -0.0249.
Optimization successful for k1=-0.6539. Start fval -3.817148, end fval -3.817148.
Next guess for k1 in direction -1 is -0.6788. Step size: -0.0250.
Optimization successful for k1=-0.6788. Start fval -3.369341, end fval -3.369341.
Next guess for k1 in direction -1 is -0.7064. Step size: -0.0276.
Optimization successful for k1=-0.7064. Start fval -2.697999, end fval -2.697999.
Next guess for k1 in direction -1 is -0.7381. Step size: -0.0317.
Optimization successful for k1=-0.7381. Start fval -1.691280, end fval -1.691280.
Next guess for k1 in direction -1 is -0.7756. Step size: -0.0374.
Optimization successful for k1=-0.7756. Start fval -0.180090, end fval -0.180090.
Next guess for k1 in direction -1 is -0.8205. Step size: -0.0449.
Optimization successful for k1=-0.8205. Start fval 2.085400, end fval 2.085400.
Next guess for k1 in direction 1 is -0.4755. Step size: 0.0225.
Optimization successful for k1=-0.4755. Start fval -4.663260, end fval -4.698178.
Next guess for k1 in direction 1 is -0.4704. Step size: 0.0051.
Optimization successful for k1=-0.4704. Start fval -4.690220, end fval -4.690228.
Next guess for k1 in direction 1 is -0.4641. Step size: 0.0063.
Optimization successful for k1=-0.4641. Start fval -4.678294, end fval -4.678296.
Next guess for k1 in direction 1 is -0.4563. Step size: 0.0078.
Optimization successful for k1=-0.4563. Start fval -4.660396, end fval -4.660399.
Next guess for k1 in direction 1 is -0.4467. Step size: 0.0096.
Optimization successful for k1=-0.4467. Start fval -4.633558, end fval -4.633565.
Next guess for k1 in direction 1 is -0.4348. Step size: 0.0119.
Optimization successful for k1=-0.4348. Start fval -4.593317, end fval -4.593332.
Next guess for k1 in direction 1 is -0.4201. Step size: 0.0147.
Optimization successful for k1=-0.4201. Start fval -4.532981, end fval -4.533013.
Next guess for k1 in direction 1 is -0.4017. Step size: 0.0184.
Optimization successful for k1=-0.4017. Start fval -4.442518, end fval -4.442583.
Next guess for k1 in direction 1 is -0.3786. Step size: 0.0231.
Optimization successful for k1=-0.3786. Start fval -4.306885, end fval -4.307010.
Next guess for k1 in direction 1 is -0.3492. Step size: 0.0294.
Optimization successful for k1=-0.3492. Start fval -4.103506, end fval -4.103738.
Next guess for k1 in direction 1 is -0.3113. Step size: 0.0379.
Optimization successful for k1=-0.3113. Start fval -3.798629, end fval -3.799093.
Next guess for k1 in direction 1 is -0.2613. Step size: 0.0500.
Optimization successful for k1=-0.2613. Start fval -3.341812, end fval -3.342439.
Next guess for k1 in direction 1 is -0.1929. Step size: 0.0684.
Optimization successful for k1=-0.1929. Start fval -2.659091, end fval -2.660779.
Next guess for k1 in direction 1 is -0.0929. Step size: 0.1000.
Optimization successful for k1=-0.0929. Start fval -1.581614, end fval -1.655512.
Next guess for k1 in direction 1 is 0.0071. Step size: 0.1000.
Optimization successful for k1=0.0071. Start fval -0.698783, end fval -0.750008.
Next guess for k1 in direction 1 is 0.1071. Step size: 0.1000.
Optimization successful for k1=0.1071. Start fval 0.022981, end fval -0.010441.
Next guess for k2 in direction -1 is -1.3362. Step size: -0.0858.
Optimization successful for k2=-1.3362. Start fval -4.663138, end fval -4.697805.
Next guess for k2 in direction -1 is -1.3572. Step size: -0.0211.
Optimization successful for k2=-1.3572. Start fval -4.689665, end fval -4.689686.
Next guess for k2 in direction -1 is -1.3842. Step size: -0.0269.
Optimization successful for k2=-1.3842. Start fval -4.677489, end fval -4.677496.
Next guess for k2 in direction -1 is -1.4190. Step size: -0.0349.
Optimization successful for k2=-1.4190. Start fval -4.659206, end fval -4.659221.
Next guess for k2 in direction -1 is -1.4650. Step size: -0.0460.
Optimization successful for k2=-1.4650. Start fval -4.631784, end fval -4.631824.
Next guess for k2 in direction -1 is -1.5273. Step size: -0.0623.
Optimization successful for k2=-1.5273. Start fval -4.590702, end fval -4.590813.
Next guess for k2 in direction -1 is -1.6153. Step size: -0.0880.
Optimization successful for k2=-1.6153. Start fval -4.529174, end fval -4.529517.
Next guess for k2 in direction -1 is -1.7153. Step size: -0.1000.
Optimization successful for k2=-1.7153. Start fval -4.460174, end fval -4.460683.
Next guess for k2 in direction -1 is -1.8153. Step size: -0.1000.
Optimization successful for k2=-1.8153. Start fval -4.396362, end fval -4.396768.
Next guess for k2 in direction -1 is -1.9153. Step size: -0.1000.
Optimization successful for k2=-1.9153. Start fval -4.339589, end fval -4.339863.
Next guess for k2 in direction -1 is -2.0000. Step size: -0.0847.
Optimization successful for k2=-2.0000. Start fval -4.297226, end fval -4.297675.
Next guess for k2 in direction 1 is -1.1762. Step size: 0.0741.
Optimization successful for k2=-1.1762. Start fval -4.663090, end fval -4.698386.
Next guess for k2 in direction 1 is -1.1606. Step size: 0.0156.
Optimization successful for k2=-1.1606. Start fval -4.690534, end fval -4.690545.
Next guess for k2 in direction 1 is -1.1420. Step size: 0.0186.
Optimization successful for k2=-1.1420. Start fval -4.678774, end fval -4.678777.
Next guess for k2 in direction 1 is -1.1198. Step size: 0.0222.
Optimization successful for k2=-1.1198. Start fval -4.661139, end fval -4.661144.
Next guess for k2 in direction 1 is -1.0936. Step size: 0.0263.
Optimization successful for k2=-1.0936. Start fval -4.634706, end fval -4.634719.
Next guess for k2 in direction 1 is -1.0625. Step size: 0.0311.
Optimization successful for k2=-1.0625. Start fval -4.595069, end fval -4.595094.
Next guess for k2 in direction 1 is -1.0259. Step size: 0.0366.
Optimization successful for k2=-1.0259. Start fval -4.535618, end fval -4.535667.
Next guess for k2 in direction 1 is -0.9830. Step size: 0.0429.
Optimization successful for k2=-0.9830. Start fval -4.446582, end fval -4.446687.
Next guess for k2 in direction 1 is -0.9327. Step size: 0.0503.
Optimization successful for k2=-0.9327. Start fval -4.313188, end fval -4.313410.
Next guess for k2 in direction 1 is -0.8736. Step size: 0.0591.
Optimization successful for k2=-0.8736. Start fval -4.113156, end fval -4.113669.
Next guess for k2 in direction 1 is -0.8038. Step size: 0.0698.
Optimization successful for k2=-0.8038. Start fval -3.813716, end fval -3.814802.
Next guess for k2 in direction 1 is -0.7200. Step size: 0.0838.
Optimization successful for k2=-0.7200. Start fval -3.365726, end fval -3.368146.
Next guess for k2 in direction 1 is -0.6200. Step size: 0.1000.
Optimization successful for k2=-0.6200. Start fval -2.721878, end fval -2.726970.
Next guess for k2 in direction 1 is -0.5200. Step size: 0.1000.
Optimization successful for k2=-0.5200. Start fval -2.003354, end fval -2.009104.
Next guess for k2 in direction 1 is -0.4200. Step size: 0.1000.
Optimization successful for k2=-0.4200. Start fval -1.266921, end fval -1.271651.
Next guess for k2 in direction 1 is -0.3200. Step size: 0.1000.
Optimization successful for k2=-0.3200. Start fval -0.569721, end fval -0.573188.
Next guess for k2 in direction 1 is -0.2200. Step size: 0.1000.
Optimization successful for k2=-0.2200. Start fval 0.040669, end fval 0.038342.

# 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.0212118230000016

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

/home/docs/checkouts/readthedocs.org/user_builds/pypesto/envs/v0.5.9/lib/python3.11/site-packages/pypesto/visualize/sampling.py:1221: 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.808626 -0.926512
ensemble_std        0.453656  0.453671
ensemble_median    -0.808626 -0.926512
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     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 0x781e1b9b9110>,
 'std': <pypesto.result.predict.PredictionResult at 0x781e1b9bb4d0>,
 'median': <pypesto.result.predict.PredictionResult at 0x781e1b9b9c50>,
 'percentile 1': <pypesto.result.predict.PredictionResult at 0x781e1b9ba4d0>,
 'percentile 5': <pypesto.result.predict.PredictionResult at 0x781e1b9b91d0>,
 'percentile 10': <pypesto.result.predict.PredictionResult at 0x781e1b9ba990>,
 'percentile 25': <pypesto.result.predict.PredictionResult at 0x781e1b9b9f10>,
 'percentile 75': <pypesto.result.predict.PredictionResult at 0x781e1b9b81d0>,
 'percentile 90': <pypesto.result.predict.PredictionResult at 0x781e1b9b9890>,
 'percentile 95': <pypesto.result.predict.PredictionResult at 0x781e1b9ba950>,
 'percentile 99': <pypesto.result.predict.PredictionResult at 0x781e1b9ba910>}