examples/kinetics_decay_chain_units_varied_params.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter variation studies of kinetic systems\n",
    "This notebook shows how one can explore the impact of a certain parameter on a kinetic model. We will also use units explicitly for our parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "from itertools import chain\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "from chempy import Substance, Reaction, ReactionSystem\n",
    "from chempy.kinetics.rates import Arrhenius, MassAction\n",
    "from chempy.kinetics.ode import get_odesys\n",
    "from chempy.printing.numbers import number_to_scientific_latex\n",
    "from chempy.units import SI_base_registry, default_units as u\n",
    "sp.init_printing()\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a generic model representing a decay-chain with two decays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, B, C = map(Substance, 'ABC')\n",
    "r1 = Reaction({'A'}, {'B'}, MassAction(Arrhenius(unique_keys=('A1', 'Ea_R_1'))))\n",
    "r2 = Reaction({'B'}, {'C'}, MassAction(Arrhenius(unique_keys=('A2', 'Ea_R_2'))))\n",
    "rsys = ReactionSystem([r1, r2])\n",
    "rsys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Arrhenius\" behaviour means that the rate of reaction depends exponentially on the inverse absolute temperature.\n",
    "\n",
    "We will use units on all our parameters in this notebook. This will prevent us from incorrect conversions or using parameters of the wrong dimensionality where they don't belong:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {'A1': 1e11/u.s, 'A2': 2e11/u.s, 'Ea_R_1': 8e3*u.K, 'Ea_R_2': 8.5e3*u.K, 'temperature': 300*u.K}\n",
    "c0 = defaultdict(lambda: 0*u.molar, {'A': 1*u.molar})\n",
    "variables = c0.copy()\n",
    "variables.update(params)\n",
    "rsys.rates(variables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "odesys, extra = get_odesys(rsys, include_params=False, lower_bounds=0)\n",
    "print(dict(zip(odesys.dep, odesys.names)))\n",
    "print(dict(zip(odesys.params, odesys.param_names)))\n",
    "odesys.exprs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the influence of ``Ea_R_2``, we will choose three temperatures: 8100, 8200 and 8300 K (this is all fictive so never mind the very high temperatures):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params2 = params.copy()\n",
    "pk = 'Ea_R_2'\n",
    "params2[pk] = [8.1e3, 8.2e3, 8.3e3]*u.K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the integartion & plotting the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res2 = odesys.integrate(7*u.s, c0, params2, integrator='cvode')\n",
    "fig, axes = plt.subplots(1, len(res2), figsize=(14, 4))\n",
    "for r, ax in zip(res2, axes):\n",
    "    r.plot(ax=ax)\n",
    "    ax.set_title('$%s = %s$' % (pk.replace('_', '\\\\_'), number_to_scientific_latex(r.named_param('Ea_R_2'))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use ``ipywidgets`` to get interactive controls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_and_plot(T_C=25):\n",
    "    res = odesys.integrate(7*u.s, c0, dict(params, temperature=(T_C+273.15)*u.K), integrator='cvode')\n",
    "    res.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(integrate_and_plot)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}