{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "# RoadRunner in pyPESTO\n",
    "\n",
    "**After going through this notebook, you will be able to...**\n",
    "\n",
    "* ... create a pyPESTO problem using [RoadRunner](https://www.libroadrunner.org) as a simulator directly from a PEtab problem.\n",
    "* ... perform a parameter estimation using pyPESTO with RoadRunner as a simulator, setting advanced simulator features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# install pyPESTO with roadrunner support\n",
    "# %pip install pypesto[roadrunner,petab] --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# import\n",
    "from pprint import pprint\n",
    "\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import petab\n",
    "from IPython.display import Markdown, display\n",
    "\n",
    "import pypesto.objective\n",
    "import pypesto.objective.roadrunner as pypesto_rr\n",
    "import pypesto.optimize as optimize\n",
    "import pypesto.petab\n",
    "import pypesto.visualize as visualize\n",
    "\n",
    "mpl.rcParams[\"figure.dpi\"] = 100\n",
    "mpl.rcParams[\"font.size\"] = 18\n",
    "\n",
    "np.random.seed(1912)\n",
    "\n",
    "\n",
    "# name of the model that will also be the name of the python module\n",
    "model_name = \"conversion_reaction\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Creating pyPESTO problem from PEtab\n",
    "\n",
    "The [PEtab file format](https://petab.readthedocs.io/en/latest/documentation_data_format.html) stores all the necessary information to define a parameter estimation problem. This includes the model, the experimental data, the parameters to estimate, and the experimental conditions. Using the `pypesto.petab.PetabImporter` class, we can create a pyPESTO problem directly from a PEtab problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "petab_yaml = f\"./{model_name}/{model_name}.yaml\"\n",
    "\n",
    "petab_problem = petab.Problem.from_yaml(petab_yaml)\n",
    "importer = pypesto.petab.PetabImporter(\n",
    "    petab_problem, simulator_type=\"roadrunner\"\n",
    ")\n",
    "problem = importer.create_problem()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "We now have a pyPESTO problem that we can use to perform parameter estimation. We can get some information on the RoadRunnerObjective and access the RoadRunner model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "pprint(problem.objective.get_config())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# direct simulation of the model using roadrunner\n",
    "sim_res = problem.objective.roadrunner_instance.simulate(\n",
    "    times=[0, 2.5, 5, 10, 20, 50]\n",
    ")\n",
    "pprint(sim_res)\n",
    "problem.objective.roadrunner_instance.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "For more details on interacting with the roadrunner instance, we refer to the [documentation of libroadrunner](https://libroadrunner.readthedocs.io/en/latest/). However, we point out that while theoretical possible, we **strongly advice against** changing the model in that manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ret = problem.objective(\n",
    "    petab_problem.get_x_nominal(fixed=False, scaled=True),\n",
    "    mode=\"mode_fun\",\n",
    "    return_dict=True,\n",
    ")\n",
    "pprint(ret)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Optimization\n",
    "\n",
    "To optimize a problem using a RoadRunner objective, we can set additional solver options for the ODE solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "optimizer = optimize.ScipyOptimizer()\n",
    "\n",
    "solver_options = pypesto_rr.SolverOptions(\n",
    "    relative_tolerance=1e-6, absolute_tolerance=1e-12, maximum_num_steps=10000\n",
    ")\n",
    "engine = pypesto.engine.SingleCoreEngine()\n",
    "problem.objective.solver_options = solver_options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "result = optimize.minimize(\n",
    "    problem=problem,\n",
    "    optimizer=optimizer,\n",
    "    n_starts=5,  # usually a value >= 100 should be used\n",
    "    engine=engine,\n",
    ")\n",
    "display(Markdown(result.summary()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "Disclaimer: Currently there are two main things not yet fully supported with roadrunner objectives. One  is parallelization of the optimization using MultiProcessEngine. The other is explicit gradients of the objective function. While the former will be added in a near future version, we will show a workaround for the latter, until it is implemented."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Visualization Methods\n",
    "\n",
    "In order to visualize the optimization, there are a few things possible. For a more extensive explanation we refer to the \"getting started\" notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "visualize.waterfall(result);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "visualize.parameters(result);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "visualize.parameters_correlation_matrix(result);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Sensitivities via finite differences\n",
    "\n",
    "Some solvers need a way to calculate the sensitivities, which currently RoadRunner Objectives do not suport. For this scenario, we can use the FiniteDifferences objective in pypesto, which wraps a given objective into one, that computes sensitivities via finite differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# no support for sensitivities\n",
    "try:\n",
    "    ret = problem.objective(\n",
    "        petab_problem.x_nominal_free_scaled,\n",
    "        mode=\"mode_fun\",\n",
    "        return_dict=True,\n",
    "        sensi_orders=(1,),\n",
    "    )\n",
    "    pprint(ret)\n",
    "except ValueError as e:\n",
    "    pprint(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "objective_fd = pypesto.objective.FD(problem.objective)\n",
    "# support through finite differences\n",
    "try:\n",
    "    ret = objective_fd(\n",
    "        petab_problem.x_nominal_scaled,\n",
    "        mode=\"mode_fun\",\n",
    "        return_dict=True,\n",
    "        sensi_orders=(1,),\n",
    "    )\n",
    "    pprint(ret)\n",
    "except ValueError as e:\n",
    "    pprint(e)"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}