{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Conversion reaction\n",
    "==================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# install if not done yet\n",
    "# !apt install libatlas-base-dev swig\n",
    "# %pip install pypesto[amici] --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import importlib\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import amici\n",
    "import amici.plotting\n",
    "import numpy as np\n",
    "\n",
    "import pypesto\n",
    "import pypesto.optimize as optimize\n",
    "import pypesto.visualize as visualize\n",
    "\n",
    "# sbml file we want to import\n",
    "sbml_file = \"conversion_reaction/model_conversion_reaction.xml\"\n",
    "# name of the model that will also be the name of the python module\n",
    "model_name = \"model_conversion_reaction\"\n",
    "# directory to which the generated model code is written\n",
    "model_output_dir = \"tmp/\" + model_name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Compile AMICI model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# import sbml model, compile and generate amici module\n",
    "sbml_importer = amici.SbmlImporter(sbml_file)\n",
    "sbml_importer.sbml2amici(model_name, model_output_dir, verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Load AMICI model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load amici module (the usual starting point later for the analysis)\n",
    "sys.path.insert(0, os.path.abspath(model_output_dir))\n",
    "model_module = importlib.import_module(model_name)\n",
    "model = model_module.getModel()\n",
    "model.requireSensitivitiesForAllParameters()\n",
    "model.setTimepoints(np.linspace(0, 10, 11))\n",
    "model.setParameterScale(amici.ParameterScaling.log10)\n",
    "model.setParameters([-0.3, -0.7])\n",
    "solver = model.getSolver()\n",
    "solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n",
    "solver.setSensitivityOrder(amici.SensitivityOrder.first)\n",
    "\n",
    "# how to run amici now:\n",
    "rdata = amici.runAmiciSimulation(model, solver, None)\n",
    "amici.plotting.plotStateTrajectories(rdata)\n",
    "edata = amici.ExpData(rdata, 0.2, 0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create objective function from amici model\n",
    "# pesto.AmiciObjective is derived from pesto.Objective,\n",
    "# the general pesto objective function class\n",
    "objective = pypesto.AmiciObjective(model, solver, [edata], 1)\n",
    "\n",
    "# create optimizer object which contains all information for doing the optimization\n",
    "optimizer = optimize.ScipyOptimizer(method=\"ls_trf\")\n",
    "\n",
    "# create problem object containing all information on the problem to be solved\n",
    "problem = pypesto.Problem(objective=objective, lb=[-2, -2], ub=[2, 2])\n",
    "\n",
    "# do the optimization\n",
    "result = optimize.minimize(\n",
    "    problem=problem, optimizer=optimizer, n_starts=100, filename=None\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "visualize.waterfall(result)\n",
    "visualize.parameters(result)\n",
    "visualize.optimization_scatter(result=result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pypesto.profile as profile\n",
    "\n",
    "profile_options = profile.ProfileOptions(\n",
    "    min_step_size=0.0005,\n",
    "    delta_ratio_max=0.05,\n",
    "    default_step_size=0.005,\n",
    "    ratio_min=0.01,\n",
    ")\n",
    "\n",
    "result = profile.parameter_profile(\n",
    "    problem=problem,\n",
    "    result=result,\n",
    "    optimizer=optimizer,\n",
    "    profile_index=np.array([1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0]),\n",
    "    result_index=0,\n",
    "    profile_options=profile_options,\n",
    "    filename=None,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# specify the parameters, for which profiles should be computed\n",
    "ax = visualize.profiles(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pypesto.sample as sample\n",
    "\n",
    "sampler = sample.AdaptiveParallelTemperingSampler(\n",
    "    internal_sampler=sample.AdaptiveMetropolisSampler(), n_chains=3\n",
    ")\n",
    "\n",
    "result = sample.sample(\n",
    "    problem, n_samples=10000, sampler=sampler, result=result, filename=None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ax = visualize.sampling_scatter(result, size=[13, 6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Predict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Let's create a function, which predicts the ratio of x_1 and x_0\n",
    "import pypesto.predict as predict\n",
    "\n",
    "\n",
    "def ratio_function(amici_output_list):\n",
    "    # This (optional) function post-processes the results from AMICI and must accept one input:\n",
    "    # a list of dicts, with the fields t, x, y[, sx, sy - if sensi_orders includes 1]\n",
    "    # We need to specify the simulation condition: here, we only have one, i.e., it's [0]\n",
    "    x = amici_output_list[0][\"x\"]\n",
    "    ratio = x[:, 1] / x[:, 0]\n",
    "    # we need to output also at least one simulation condition\n",
    "    return [ratio]\n",
    "\n",
    "\n",
    "# create pypesto prediction function\n",
    "predictor = predict.AmiciPredictor(\n",
    "    objective, post_processor=ratio_function, output_ids=[\"ratio\"]\n",
    ")\n",
    "\n",
    "# run prediction\n",
    "prediction = predictor(x=model.getUnscaledParameters())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "dict(prediction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Analyze parameter ensembles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# We want to use the sample result to create a prediction\n",
    "from pypesto.ensemble import ensemble\n",
    "\n",
    "# first collect some vectors from the sampling result\n",
    "vectors = result.sample_result.trace_x[0, ::20, :]\n",
    "\n",
    "# create the collection\n",
    "my_ensemble = ensemble.Ensemble(\n",
    "    vectors,\n",
    "    x_names=problem.x_names,\n",
    "    ensemble_type=ensemble.EnsembleType.sample,\n",
    "    lower_bound=problem.lb,\n",
    "    upper_bound=problem.ub,\n",
    ")\n",
    "\n",
    "# we can also perform an approximative identifiability analysis\n",
    "summary = my_ensemble.compute_summary()\n",
    "identifiability = my_ensemble.check_identifiability()\n",
    "print(identifiability.transpose())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# run a prediction\n",
    "ensemble_prediction = my_ensemble.predict(\n",
    "    predictor, prediction_id=\"ratio_over_time\"\n",
    ")\n",
    "\n",
    "# go for some analysis\n",
    "prediction_summary = ensemble_prediction.compute_summary(\n",
    "    percentiles_list=(1, 5, 10, 25, 75, 90, 95, 99)\n",
    ")\n",
    "dict(prediction_summary)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}