{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting kinetic parameters to experimental data\n",
    "Author: Björn Dahlgren.\n",
    "\n",
    "Let us consider the reaction:\n",
    "\n",
    "$$\n",
    "Fe^{3+} + SCN^- \\rightarrow FeSCN^{2+}\n",
    "$$\n",
    "\n",
    "the product is strongly coloured and we have experimental data (from a stopped-flow appartus) of the absorbance as function of time after mixing for several replicates. The experiment was performed at 7 different temperatures and for one temperature, 7 different ionic strengths. For each set of conditions the experiment was re-run 7 times (replicates). In this notebook, we will determine the activation enthalpy and entropy through regressions analysis, we will also look at the ionic strength dependence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bz2, codecs, collections, functools, itertools, json\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import chempy\n",
    "import chempy.equilibria\n",
    "from chempy.electrolytes import ionic_strength\n",
    "from chempy.kinetics.arrhenius import fit_arrhenius_equation\n",
    "from chempy.printing import number_to_scientific_latex, as_per_substance_html_table\n",
    "from chempy.properties.water_density_tanaka_2001 import water_density\n",
    "from chempy.units import rescale, to_unitless, default_units as u\n",
    "from chempy.util.regression import least_squares, irls, avg_params, plot_fit, plot_least_squares_fit, plot_avg_params\n",
    "from chempy._solution import QuantityDict\n",
    "%matplotlib inline\n",
    "print(chempy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Experimental conditions, the two solutions which were mixed in 1:1 volume ratio in a stopped flow apparatus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol1 = QuantityDict(u.mM, {'SCN-': 3*u.mM, 'K+': 3*u.mM, 'Na+': 33*u.mM, 'H+': 50*u.mM, 'ClO4-': (33+50)*u.mM})\n",
    "sol2 = QuantityDict(u.mM, {'Fe+3': 6*u.mM, 'H+': 50*u.mM, 'ClO4-': (3*6+50)*u.mM})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = (sol1 + sol2)/2  # 1:1 volume ratio at mixing\n",
    "sol.quantity_name = 'concentration'\n",
    "Ibase = ionic_strength(rescale(sol/water_density(293*u.K, units=u), u.molal))\n",
    "print(Ibase)\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ionic_strength_keys = 'abcd'\n",
    "ionic_strengths = dict(zip(ionic_strength_keys, [0, 20, 40, 60]*u.molal*1e-3))\n",
    "temperature_keys = '16.5 18.5 20.5 22.5 24.5'.split()\n",
    "T0C = 273.15*u.K\n",
    "temperatures = {k: T0C + float(k)*u.K for k in temperature_keys}\n",
    "nrep = 7\n",
    "indices = lambda k: (ionic_strength_keys.index(k[0]), temperature_keys.index(k[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will read the data from a preprocessed file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = np.array([[1e-3, 0], [0, 1e-4]])  # converts 1st col: ms -> s and 2nd col to absorbance\n",
    "_reader = codecs.getreader(\"utf-8\")\n",
    "_dat = {tuple(k): np.dot(np.array(v), transform) for k, v in json.load(_reader(bz2.BZ2File('specdata.json.bz2')))}\n",
    "data = collections.defaultdict(list)\n",
    "for (tI, tT, tR), v in _dat.items():\n",
    "    k = (tI, tT)  # tokens for ionic strength and temperatures\n",
    "    data[k].append(v)\n",
    "assert len(data) == len(ionic_strengths)*len(temperatures) and all(len(serie) == nrep for serie in data.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mk_subplots(nrows=1, subplots_adjust=True, **kwargs):\n",
    "    fig, axes = plt.subplots(nrows, len(ionic_strengths), figsize=(15,6), **kwargs)\n",
    "    if subplots_adjust:\n",
    "        plt.subplots_adjust(hspace=0.001, wspace=0.001)\n",
    "    return axes\n",
    "\n",
    "def _set_axes_titles_to_ionic_strength(axes, xlim=None, xlabel=None):\n",
    "    for tI, ax in zip(ionic_strength_keys, axes):\n",
    "        ax.set_title(r'$I\\ =\\ %s$' % number_to_scientific_latex(Ibase + ionic_strengths[tI], fmt=3))\n",
    "        if xlabel is not None:\n",
    "            ax.set_xlabel(xlabel)\n",
    "        if xlim is not None:\n",
    "            ax.set_xlim(xlim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_series(series):\n",
    "    axes = mk_subplots(sharey=True, sharex=True)\n",
    "    colors = 'rgbmk'\n",
    "    for key in itertools.product(ionic_strengths, temperatures):\n",
    "        idx_I, idx_T = indices(key)\n",
    "        for serie in series[key]:\n",
    "            axes[idx_I].plot(serie[:, 0], serie[:, 1], c=colors[idx_T], alpha=0.15)\n",
    "    _set_axes_titles_to_ionic_strength(axes, xlim=[0, 3.4], xlabel='Time / s')\n",
    "    axes[0].set_ylabel('Absorbance')\n",
    "    for c, tT in zip(reversed(colors), reversed(temperature_keys)):\n",
    "        axes[0].plot([], [], c=c, label=tT + ' °C')\n",
    "    axes[0].legend(loc='best')\n",
    "plot_series(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that one data series is off: 16.5 ℃ and 0.0862 molal. Let's ignore that for now and perform the fitting, let's start with a pseudo-first order assumption (poor but simple):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_pseudo1(serie, ax=None):\n",
    "    plateau = np.mean(serie[2*serie.shape[0]//3:, 1])\n",
    "    y = np.log(np.clip(plateau - serie[:, 1], 1e-6, 1))\n",
    "    x = serie[:, 0]\n",
    "    # irls: Iteratively reweighted least squares\n",
    "    res = irls(x, y, irls.gaussian)\n",
    "    if ax is not None:\n",
    "        plot_least_squares_fit(x, y, res)\n",
    "    return res\n",
    "\n",
    "beta, vcv, info = fit_pseudo1(data['a', '16.5'][0], ax=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_all(series, fit_cb=fit_pseudo1, plot=False):\n",
    "    if plot:\n",
    "        axes = mk_subplots(nrows=len(temperatures), sharex=True, sharey=True)#, subplots_adjust=False)\n",
    "        _set_axes_titles_to_ionic_strength(axes[0, :])\n",
    "    avg = {}\n",
    "    for key in itertools.product(ionic_strengths, temperatures):\n",
    "        idx_I, idx_T = indices(key)\n",
    "        opt_params, cov_params = [], []\n",
    "        for serie in series[key]:\n",
    "            beta, vcv, nfo = fit_cb(serie)\n",
    "            opt_params.append(beta)\n",
    "            cov_params.append(vcv)\n",
    "        ax = axes[idx_T, idx_I] if plot else None\n",
    "        avg[key] = avg_params(opt_params, cov_params)\n",
    "        plot_avg_params(opt_params, cov_params, avg[key], ax=ax, flip=True, nsigma=3)\n",
    "    for tk, ax in zip(temperature_keys, axes[:, 0]):\n",
    "        ax.set_ylabel('T = %s C' % tk)\n",
    "    return avg\n",
    "\n",
    "result_pseudo1 = fit_all(data, plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pseudo_to_k2(v):\n",
    "    unit = 1/sol['Fe+3']/u.second\n",
    "    k_val = -v[0][1]*unit\n",
    "    k_err = v[1][1]*unit\n",
    "    return k_val, k_err\n",
    "\n",
    "k_pseudo1 = {k: pseudo_to_k2(v) for k, v in result_pseudo1.items()}\n",
    "k_pseudo1['a', '16.5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k2_unit = 1/u.M/u.s\n",
    "axes = mk_subplots(sharex=True, sharey=True)\n",
    "for idxI, (tI, I) in enumerate(ionic_strengths.items()):\n",
    "    series = np.empty((len(temperatures), 3))\n",
    "    for idxT, (tT, T) in enumerate(temperatures.items()):\n",
    "        kval, kerr = [to_unitless(v, k2_unit) for v in k_pseudo1[tI, tT]]\n",
    "        lnk_err = (np.log(kval + kerr) - np.log(kval - kerr))/2\n",
    "        series[idxT, :] = to_unitless(1/T, 1/u.K), np.log(kval), 1/lnk_err**2\n",
    "    x, y, w = series.T\n",
    "    res = b, vcv, r2 = least_squares(x, y, w)\n",
    "    plot_least_squares_fit(x, y, res, w**-0.5, plot_cb=functools.partial(\n",
    "            plot_fit, ax=axes[idxI], kw_data=dict(ls='None', marker='.'), nsigma=3))\n",
    "    axes[idxI].get_lines()[-1].set_label(r'$E_{\\rm a} = %.5g\\ kJ/mol$' % (-b[1]*8.314511e-3))\n",
    "    axes[idxI].legend()\n",
    "_set_axes_titles_to_ionic_strength(axes)"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}