ABCD_parmest_multiexp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "import casadi as csd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg\n",
    "from scipy.stats import chi2, f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "# experimental conditions\n",
    "expcond = [{'c0': [1.0, 2.0, 0.0, 0.0]}, {'c0': [1.0, 1.0, 0.0, 0.0]}]\n",
    "meas_vars = [['ca', 'cb', 'cc'], ['ca', 'cc', 'cd']]\n",
    "meas_vars_idx = [[0, 1, 2], [0, 2, 3]]\n",
    "datasets = ['ABCD_data.csv', 'ABCD_data_2.csv']\n",
    "expdata = []\n",
    "for data in datasets:\n",
    "    data_df = pd.read_csv(data)\n",
    "    expdata.append(data_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stoichiometry matrix for\n",
    "# A + B -> C, B + C -> D\n",
    "s = np.array([[-1.0, -1.0, 1.0, 0.0],\n",
    "              [0.0, -1.0, -1.0, 1.0]\n",
    "             ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rxnfn(kf, tf, s, tgrid = [-100]):\n",
    "    \n",
    "    nr, nc = s.shape\n",
    "    \n",
    "    c = csd.MX.sym('c', nc)\n",
    "    \n",
    "    r = []\n",
    "    for i in range(nr):\n",
    "        ri = kf[i]\n",
    "        for j in range(nc):\n",
    "            if s[i, j] < 0:\n",
    "                ri = ri * c[j] ** (-s[i, j])\n",
    "        r.append(ri) \n",
    "        \n",
    "    dc = []\n",
    "    for i in range(nc):\n",
    "        dci = 0\n",
    "        for j in range(nr):\n",
    "            dci = dci + s[j, i] * r[j]\n",
    "        dc.append(dci)\n",
    "        \n",
    "    ode = {}\n",
    "    ode['x'] = c\n",
    "    ode['p'] = kf\n",
    "    ode['ode'] = csd.vertcat(*dc)\n",
    "    \n",
    "    if tgrid[0] == -100:\n",
    "        F = csd.integrator('F','cvodes',ode,{'tf':tf})\n",
    "    else:\n",
    "        F = csd.integrator('F','cvodes',ode,{'tf':tf, 'grid': tgrid, 'output_t0' : True})\n",
    "        \n",
    "    return F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "expinfo_list = [{'data': expdata[i], 'meas_var': meas_vars[i], 'meas_var_idx': meas_vars_idx[i], 'c0': expcond[i]['c0']}\n",
    "                  for i in range(len(expdata))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_exp_ssq(kf, expinfo):\n",
    "    data = expinfo['data']\n",
    "    meas_var = expinfo['meas_var']\n",
    "    meas_var_idx = expinfo['meas_var_idx']\n",
    "    c0 = expinfo['c0']\n",
    "    tgrid = np.append(0, data['t'].values)\n",
    "    ssq = 0\n",
    "    for i in range(len(tgrid) - 1):\n",
    "        F = rxnfn(kf = kf, tf = tgrid[i + 1] - tgrid[i], s = s)\n",
    "        res = F(x0 = c0, p = kf)\n",
    "        c0 = res['xf']\n",
    "        for (j, var) in enumerate(meas_var):\n",
    "            ssq = ssq + (data.iloc[i][var] - res['xf'][meas_var_idx[j]]) ** 2\n",
    "    return ssq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_exp(kf, expinfo, tf):\n",
    "    data = expinfo['data']\n",
    "    meas_var = expinfo['meas_var']\n",
    "    meas_var_idx = expinfo['meas_var_idx']\n",
    "    c0 = expinfo['c0']\n",
    "    tgrid = list(np.linspace(0, tf))\n",
    "    F = rxnfn(kf = kf, tf = tf, s = s, tgrid = tgrid)\n",
    "    res = F(x0 = c0, p = kf)\n",
    "    res_fn = csd.Function('res_fn', [kf], [res['xf']])\n",
    "    return res_fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "kf = csd.MX.sym('kf', 2)\n",
    "exp_ssq = 0\n",
    "for i in range(len(expdata)):\n",
    "    exp_ssq = exp_ssq + get_exp_ssq(kf, expinfo_list[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ssq_fn = csd.Function('exp_ssq_fn', [kf], [exp_ssq])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calcuate hessian of sum of squares with respect to p = (k1, k2)\n",
    "ssqfn_hess_calc = csd.hessian(exp_ssq_fn(kf), kf)\n",
    "ssqfn_hess = csd.Function('ssqfn_hess', [kf], [ssqfn_hess_calc[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is Ipopt version 3.12.3, running with linear solver mumps.\n",
      "NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:        0\n",
      "Number of nonzeros in Lagrangian Hessian.............:        3\n",
      "\n",
      "Total number of variables............................:        2\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:        2\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:        0\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        0\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  5.2297514e+01 0.00e+00 4.06e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1  5.2050847e+01 0.00e+00 4.11e+00  -1.0 4.10e-02   2.0 9.98e-01 1.00e+00f  1\n",
      "   2  5.1238149e+01 0.00e+00 4.42e+00  -1.0 1.32e-01   1.5 1.00e+00 1.00e+00f  1\n",
      "   3  4.7804292e+01 0.00e+00 5.95e+00  -1.0 5.07e-01   1.0 1.00e+00 1.00e+00f  1\n",
      "   4  4.5999314e+01 0.00e+00 6.92e+00  -1.0 2.30e-01   1.5 1.00e+00 1.00e+00f  1\n",
      "   5  2.6703128e+01 0.00e+00 2.41e+01  -1.0 1.40e+00   1.0 1.00e+00 1.00e+00f  1\n",
      "   6  1.7364320e+01 0.00e+00 2.64e+01  -1.0 1.26e+01   1.4 7.47e-01 2.69e-02f  2\n",
      "   7  1.4759539e+01 0.00e+00 1.74e+01  -1.0 4.58e-01   0.9 1.00e+00 3.76e-01f  2\n",
      "   8  1.3179094e+01 0.00e+00 3.34e+00  -1.0 7.23e-01   0.5 1.00e+00 1.00e+00f  1\n",
      "   9  1.2770361e+01 0.00e+00 4.49e+00  -1.0 3.11e+00    -  1.00e+00 2.35e-01f  2\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  10  1.2471252e+01 0.00e+00 1.05e+00  -1.0 1.28e-01    -  1.00e+00 1.00e+00f  1\n",
      "  11  1.2437672e+01 0.00e+00 1.60e-01  -1.0 6.87e-02    -  1.00e+00 1.00e+00f  1\n",
      "  12  1.2437199e+01 0.00e+00 4.16e-03  -2.5 6.80e-03    -  1.00e+00 1.00e+00f  1\n",
      "  13  1.2437199e+01 0.00e+00 2.73e-06  -3.8 4.28e-05    -  1.00e+00 1.00e+00f  1\n",
      "  14  1.2437199e+01 0.00e+00 4.44e-07  -5.7 1.23e-05    -  1.00e+00 1.00e+00f  1\n",
      "  15  1.2437199e+01 0.00e+00 6.01e-09  -8.6 1.69e-07    -  1.00e+00 1.00e+00f  1\n",
      "\n",
      "Number of Iterations....: 15\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   1.2437198842545014e+01    1.2437198842545014e+01\n",
      "Dual infeasibility......:   6.0075172694499997e-09    6.0075172694499997e-09\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   2.5062313939612974e-09    2.5062313939612974e-09\n",
      "Overall NLP error.......:   6.0075172694499997e-09    6.0075172694499997e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 31\n",
      "Number of objective gradient evaluations             = 16\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 0\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 0\n",
      "Number of Lagrangian Hessian evaluations             = 15\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.067\n",
      "Total CPU secs in NLP function evaluations           =      3.587\n",
      "\n",
      "EXIT: Optimal Solution Found.\n",
      "      solver  :   t_proc      (avg)   t_wall      (avg)    n_eval\n",
      "       nlp_f  | 420.61ms ( 13.57ms) 420.22ms ( 13.56ms)        31\n",
      "  nlp_grad_f  | 939.26ms ( 55.25ms) 937.78ms ( 55.16ms)        17\n",
      "  nlp_hess_l  |   3.10 s (206.44ms)   3.09 s (206.20ms)        15\n",
      "       total  |   4.48 s (  4.48 s)   4.47 s (  4.47 s)         1\n"
     ]
    }
   ],
   "source": [
    "# NLP declaration\n",
    "nlp = {'x': kf,'f': exp_ssq};\n",
    "\n",
    "# Solve using IPOPT\n",
    "solver = csd.nlpsol('solver','ipopt',nlp)\n",
    "res = solver(x0=[3, 3], lbx = 0, ubx = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DM([0.949272, 0.212212])"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_est = res['x']\n",
    "p_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_fn = sim_exp(kf, expinfo_list[1], tf = 10)\n",
    "cf = res_fn(p_est).full().T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1a238a7978>"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3xU9Z3/8dcnkzshCYFwS8AAYhRFRVMQ6UWrLqi7YtW62p/Vbm3Z7dZetq3709/2t9vtXqR1t13burularXWemnXH7IVS1V0rVQRELyAoFwEEm7hkgC5X76/P2agMUxCkjkz58yc9/PxyCNzznxnzndy+Zzv+Xwvx5xziIhI5svyuwIiIpIaCvgiIiGhgC8iEhIK+CIiIaGALyISEtl+V6Avo0aNclVVVX5XQ0QkraxZs2a/c6483nOBDfhVVVWsXr3a72qIiKQVM9ve13NK6YiIhIQCvohISCjgi4iEhAK+iEhIKOCLiISEJwHfzB4ws31m9nYfz5uZ/cDMNpvZm2Z2nhfHFRGRgfOqhf8gMK+f5y8Hpsa+FgD/4dFxRURkgDwJ+M65l4CD/RSZD/zMRb0KlJrZOC+O3VtDczv3PPceb9c1JuPtRUTSVqpy+BXAzh7btbF9H2BmC8xstZmtrq+vH9KBIlnGPc+/y2837B1aTUVEMlSqAr7F2XfCnVecc4ucczXOuZry8rgzg09qeH4OZ44vYeXWA0N6vYhIpkrV0gq1wIQe25XArmQdbNakMn726nZaO7rIz4kAsHhtHXcv28SuhhbGlxZw+9xqrp5R0ed+EZFMk6oW/hLg5thonQuARufc7mQdbNbkkbR3drNuZwMQDfZ3PvkWdQ0tOKCuoYU7n3yLby5+K+7+xWvrklU1ERHfeDUs81HgFaDazGrN7FYz+wsz+4tYkaXAVmAz8BPgL704bl9mVpVhBiu3RvuR7162iZaOrg+Uaeno4tGVO+Puv3vZpmRWT0TEF56kdJxzN57keQd80YtjDURJYQ6njy1m5bYDwFR2NbTELdfVxw3c+yovIpLOMnam7axJZby+4xDtnd2MLy2IWyZi8fqS6bO8iEg6y9iAf8HkMlo7unmztoHb51ZTEOu8PaYgJ8KNsybE3X/73OpUVlVEJCUCewOURM2cNBKAldsO8sWLTwWIOxqn5pQyjdIRkVAw10ce2281NTUu0Ttezf3+S4wuzuPhW2d5VCsRkWAzszXOuZp4z2VsSgdg1uQy1mw/REdXt99VERHxXWYH/EkjaW7v0ro6IiJkeMCfOakMiObxRUTCLqMDfvnwPKaUD9O6OiIiZHjAh+gyC6veP0Sn8vgiEnKZH/AnlXG0rZMNuw/7XRUREV9lfMC/YHJsPP5W5fFFJNwyPuCPKc6namRhbF0dEZHwyviAD9FW/mvbDtLVHcxJZiIiqRCKgD9rchmHWzvZuEd5fBEJr3AE/EnK44uIhCLgjy8tYEJZgfL4IhJqoQj4EG3lv7btIN3K44tISIUo4JdxqLmD9/Yd9bsqIiK+CE3APz4eX2kdEQmp0AT8yhHRPP4LG/f5XRUREV+EJuCbGZefNY6XN++nsaXD7+qIiKScJwHfzOaZ2SYz22xmd8R5fqKZvWBma83sTTO7wovjDtYV08fR0eV4dsNePw4vIuKrhAO+mUWAe4HLgWnAjWY2rVexbwJPOOdmADcA/57ocYfinMoSKkoLeOat3X4cXkTEV1608GcCm51zW51z7cBjwPxeZRxQHHtcAuzy4LiDFk3rjOV37+3ncKvSOiISLl4E/ApgZ4/t2ti+nr4F3GRmtcBS4Evx3sjMFpjZajNbXV9f70HVTnTF2eNo7+rmOaV1RCRkvAj4Fmdf79lNNwIPOucqgSuAh83shGM75xY552qcczXl5eUeVO1E51aWMq4kn6Vv7UnK+4uIBJUXAb8WmNBju5ITUza3Ak8AOOdeAfKBUR4ce9CysqKjdV56r54jSuuISIh4EfBXAVPNbJKZ5RLtlF3Sq8wO4BIAMzuDaMBPTs5mAK48eyztnd08/47G5ItIeCQc8J1zncBtwDLgHaKjcdab2bfN7KpYsa8DnzezN4BHgc8453xb1GbGhBGMLc5nqUbriEiIZHvxJs65pUQ7Y3vu+9sejzcAc7w4lheysox5Z43lF6/t4GhbJ0V5A/8xLF5bx93LNrGroYXxpQXcPreaq2f07qMWEQme0My07e3Ks8fR3tnN8kEstbB4bR13PvkWdQ0tOKCuoYU7n3yLxWvrkldRERGPhDbgnz9xBKOH57H0zYGnde5etomWjq4P7Gvp6OLuZZu8rp6IiOdCG/Cjo3XG8sKmfTS1dQ7oNbsaWga1X0QkSEIb8CG6tk5bZzcvbBpYWmd8acGg9ouIBEmoA35NVRmjivIGPFrn9rnVFOREPrCvICfC7XOrk1E9ERFPhTrgR2JpneUb99HcfvK0ztUzKrjrmulUlBZgQEVpAXddM12jdEQkLXgyLDOdXTF9HA+/up0XN9VzxfRxJy1/9YwKBXgRSUuhbuEDzJxUxqiiXP77DV8W8BQRSZnQB/xIlnH1uRU8u2Evew+3+l0dEZGkCX3AB7h5dhVdzvHIq9v9roqISNIo4AMTRxZyyemj+cVrO2jr7Dr5CwZp8do65ixczqQ7nmbOwuWamSsivlDAj7nlwir2H23n6UHMvB0ILccgIkGhgB/z4VNHMaV8GA/9/n1P31fLMYhIUCjgx5gZt1xYxRu1jazdcciz99VyDCISFAr4PVxzXiVFedmetvK1HIOIBIUCfg9Fedlcd34lT7+1m31HvBmiqeUYRCQoFPB7uXn2KXR0OR5dudOT99NyDCISFKFfWqG3yeVFfOy0ch5ZuZ0vXDSF3OzEz4lajkFEgkAt/Dg+M6eKfUfa+M36PX5XRUTEM2rhx/GxqeVMGhUdonnVOeMH9Brd61ZEgk4t/DiysoxPX3AKa7Yf4q3axpOW1+QqEUkHngR8M5tnZpvMbLOZ3dFHmevNbIOZrTezX3hx3GS6rqaSwtwIDw5giKYmV4lIOkg44JtZBLgXuByYBtxoZtN6lZkK3AnMcc6dCXw10eMmW3F+DtfXTOCpdXVsP9DUb1lNrhKRdOBFC38msNk5t9U51w48BszvVebzwL3OuUMAzrmB3UTWZ3950RRyIll879l3+y2nyVUikg68CPgVQM9B67WxfT2dBpxmZivM7FUzm+fBcZNudHE+n/1wFU+t28X6XX3n8jW5SkTSgRcB3+Lsc722s4GpwEXAjcB9ZlZ6whuZLTCz1Wa2ur6+3oOqJW7BR6dQUpDDd3/Tdz5ek6tEJB14MSyzFpjQY7sS6H2/wFrgVedcB7DNzDYRPQGs6lnIObcIWARQU1PT+6Thi5KCHL548RT+eelGXtlygNlTRsYtF8TJVRoqKiI9edHCXwVMNbNJZpYL3AAs6VVmMXAxgJmNIpri2erBsVPi5tlVjCvJ5zu/2YhzgTgPnZSGiopIbwkHfOdcJ3AbsAx4B3jCObfezL5tZlfFii0DDpjZBuAF4Hbn3IFEj50q+TkRvnrpVNbtbOC3G/b6XZ0B0VBREenNk5m2zrmlwNJe+/62x2MHfC32lZauPa+SRS9t5e5lm7jk9NFkR4I9Z01DRUWkt2BHrQDJjmRx+9xqNu87ypOvBz8toqGiItKbAv4gzD1zLOdMKOX7z71La4f3Nzv3koaKikhvCviDYGb873nV7G5s5eFXtvtdnX5pqKiI9KbVMgfpwimj+MjUUdz74mauOa+CkUV5flepT0EcKioi/lELfwi+eeU0mto6+fv/3uB3VUREBkwBfwiqxw7ntounsuSNXTyXwmGai9fWMWfhcibd8TRzFi7XmHoRGRQF/CH6wkVTOH3scP5m8Vs0tnQk/XiaSCUiiVLAH6Lc7Cy+e93Z1B9p466l7wzpPQbTYtdEKhFJlAJ+As6uLOXzH53MY6t2smLz/kG9drAtdk2kEpFEKeAn6K8uPY1Jo4Zxx5Nv0tzeOeDXDbbFrolUIpIoBfwE5edE+M61Z7PzYMug0iuDbbFrIpWIJEoB3wMzJ5Vx8+xTePD377Nm+8EBvWawLXZNpBKRRFlQl/utqalxq1ev9rsaA3a0rZO533+J/Jwsnv7yR8jv1Rrv7VgOv2dapyAnwl3XTAfQOvYiMiRmtsY5VxPvObXwPVKUl83Ca6ezpb6Jby5++6Tr5vfVYgfSavil5gaIpA8treChj0wt58uXTOUHz7/HjIml/K9Zp/RbPt7SB3MWLu+zMzdorfzeVynHTk5A4OoqImrhe+4rl0zlY6eV860l61m749CgX59Owy81N0AkvSjgeyySZdxzw7mMKc7nLx95nf1H2wb1+nQafunlyUmpIZHkU8BPgtLCXP7zpvM52NTOl36xls6u7gG/Np2GX3p1ctKyESKpoYDfh0RbnGdVlPBPn5jOK1sPcPdvB57iSKfhl16dnJQaEkkNddrG4VVn5HXnV7Ju5yF+/D9bObeylMunjxvQ69JlHftjdUx0CGk69VuIpDMF/Dj6a3EONpj93z+extt1h/nGL99gcnkR1WOHe1lV33lxchpfWkBdnOAexH4LkXSmlE4cXrY487Ij/MdN5zEsL5tP37+SnQebE61exkmnfguRdOZJwDezeWa2ycw2m9kd/ZS7zsycmcWdBRYUXo+UGVdSwMO3zqKts5ub7l/JviOtiVQv46RTv4VIOkt4aQUziwDvApcBtcAq4Ebn3IZe5YYDTwO5wG3OuX7XTfBzaYX+lj1IJAi9vuMQN923kollhTz+57MpKcjxoroiIscle2mFmcBm59xW51w78BgwP065fwC+CwS+eZusFud5E0fw40+fz5b6o9z64Cpa2rtO/iIREY940WlbAezssV0LzOpZwMxmABOcc782s2/09UZmtgBYADBx4kQPqjZ0yRop85Gp5dxzwwxu+8XrfOGRNSz6dA252epKEZHk8yLSWJx9x/NEZpYFfB/4+sneyDm3yDlX45yrKS8v96BqwXTF9HH80yem8+Kmer7xyzfo7g7miqUiklm8aOHXAhN6bFcCu3psDwfOAl40M4CxwBIzu+pkefxMduPMiTQ0d/Cd32wkO2J899qzyY6opS8iyeNFwF8FTDWzSUAdcAPwqWNPOucagVHHts3sReAbYQ72x3zhoil0dHXzvWff5XBLJz/61IyTrqMvIjJUCTcpnXOdwG3AMuAd4Ann3Hoz+7aZXZXo+2e6L18ylX+YfybPb9zLLQ+8xuHWDr+rJCIZSne8Coin1tXx9SfeoHrscB767ExGFeX5XSURSUO641UamH9uBffdUsOW+qN88j9fCdyMXC1fLJL+FPAD5KLq0TzyuVkcONrGJ//zFd7de8TvKgFavlgkUyjgB8z5p5TxxF/Mpts5PnHvCpat3+N3lTxdvlhXCiL+UcAPoNPHFvPUbXM4dcxw/vzhNfzrbzf5Olbfq8XkdKUg4i8F/IAaV1LA4wsu4PqaSn64fDO3PrSKxhZ/RvB4tZicbnQi4i8F/ADLz4nwnWvP5h+vPouXN+9n/o9e9iWv79XyxbrRSWZSmi59KOAHnJlx0wWn8OjnL6CpvYtP3LuCp9/cndI6eLWYXDrdoF0GRmm69KJx+Glk7+FWvvDzNby+o4Frz6vk766aRnF++iyxnKxlp8U/cxYuj3u3sorSAlbc8XEfaiT9jcPXLQ7TyJjifB7/89n88Pn3+NELm3l16wG+d/05zJo80u+qDYhX98CV4Ahqmm7x2jr9ncWhgJ9mciJZfO2Pqrno9NH81ePruOEnr7Lgo5P52mWnkZcd/HV40uUG7TIwQbwfce8ryWNpJiD0f3vK4aep8yaOYOmXP8INH5rIj/9nK/N/tIKNew77XS0JmSDej1ijwfqmgJ/GhuVlc9c107n/lhr2H23jT374Mt/9zUbdSUtSJoj3Iw5qmikIlNLJAJecMYZlX/0o/7x0I//+4haeWreLv7/qTC6dNsbvqgWK8rrJEbQ0XRDTTEGhFn6GGFmUx79efw6PL7iAwtwIn/vZaj7/s9Vx//DDSMMHwyOIaaag0LDMDNTe2c39L2/jB8+/B0TX3P+zOVWhvrlKpgwf1FXKwIT559TfsEwF/AxWe6iZv//vDTy7YS8VpQV8/Y9OY/65FUSy4t2GOLNNuuNp4v2lG7Bt4ZWprs6QaB6DDITWww+pyhGF/OTmGh753CxGDMvha0+8wZU/+B0vbtpHUE/0yZIJs3w1+kQSpYAfAnNOHcWSL36YH9w4g+b2Lj7z01V86icrebO2we+qpUwm5HU1+kQSpYAfEllZxlXnjOe5r32Mb/3JNDbtPcJVP1rBrQ+u4vUdh/yuXtIFcfjgYGXCVYr4Szn8kDrS2sFPV7zPAyu20dDcwZxTR/LFi09l9uSRmIUvx58OlMNPD353GKvTVvrU1NbJIyu3s+ilbew/2sZ5E0v50senclF1eWACv9//QEGin0WwBeGknPSAb2bzgHuACHCfc25hr+e/BnwO6ATqgc8657b3954K+KnV2tHFE6t38uP/2UpdQwtTRxdxy4VVXHNeBYW5/s3PC8I/kMhABWH4b1JH6ZhZBLgXuByYBtxoZtN6FVsL1DjnzgZ+BXw30eOKt/JzItw8u4oXvnER//LJc8jNzuKbi9/mgn9+nruWvkPtoWZf6qWRKZJOgt6x7kXTbSaw2Tm3FcDMHgPmAxuOFXDOvdCj/KvATR4cV5IgNzuL686v5NrzKliz/RA/XfE+9728jZ/8biuXTRvDpy+o4sIpI8lK0Vj+oP8DifQU9GUdvAj4FcDOHtu1wKx+yt8KPBPvCTNbACwAmDhxogdVk8GKlyP+myvP4OevbufR13awbP1eKkcU8MnzJ3BdTSUVSf5DDvo/UNAp559at8+tjpuCDMrwXy+GZcZr6sXtGDCzm4Aa4O54zzvnFjnnapxzNeXl5R5UTQajr/VmXtt2kL+edzqv3HkJP7hxBlUjh/H9597lw99Zzs0PvMbTb+6mrTM5K3Rmwvh5v2j9oNQL+vDfhDttzWw28C3n3NzY9p0Azrm7epW7FPgh8DHn3L6Tva86bVNvMB1OOw8288s1tfxy9U52N7ZSnJ/N5WeNY/6545k1eaSnyzekUys1SHUNQgeipF6yb3G4CphqZpOAOuAG4FO9KjAD+DEwbyDBXvwxmHz5hLJCvnbZaXzlkqn87r16nlq3i1+/uYvHV+9k9PA8/vjs8Vx17njOqSxJeHhn0Jbf7YuXd1ry4sSh/g/pLeGA75zrNLPbgGVEh2U+4Jxbb2bfBlY755YQTeEUAb+M/fPvcM5dleixxVtDyZdHsoyLqkdzUfVoWtq7eH7jXpas28XPX93OAyu2MaGsgLnTxjL3rLGcN3FERi/c1t+IosEEa69OHOr/kN408UqO83LMe2NLB8vW7+GZt3azYvMB2ru6GVWUy6VnjGHumWO58NSRaXEP3sHwakXOoaRi4l0RAJrDEELJTulIhjgWBLzIQZcU5HB9zQSur5nAkdYOXtxUz7L1e/j1m7t5bNVOhuVGuPDUUVxcPZqLqsszotXpVYt6sKmYvq4I7rpmOnddMz0wfQriP7XwJaXaOrv4/ZYDPLdhLy9uqj8eIE8bU8TF1aP5WHU5558yIi1b/15dIQ22ha/OWelJLXwJjLzsCBdXj+bi6tE459hSf5QXNtbz4rv7eGDFNn780lbysrP4UFUZF546kjlTRnFWRUla5P69ukIa7FjuoHbOBmnEkkSphS+BcbStk1e3HGDFlv28suUAG/ccAaA4P5sLJo9k5qQyZk4qY9q4YrIjmb2y92CCZRBb+FoDyT9aLVPSUv2RNl7ZeoDfb97P77ccYMfB6Ho+w3IjnHfKCD5UVcaHqso4Z0KJrwu8+S2IwTVVJyFdRZxIKR1JS+XD87jqnPFcdc54APY0trLq/YOsev8gr207yPefexfnokNDq8cMZ8bEUmZMHMGMiaVMGjksZev9+M3Lzvb+DCa4piLN5OW8h7BQC1/SVmNzB6/vOMTaHYdYu7OBdTsaONLWCUTTQGdXlnJWRQlnV5YwvaKEyhEFgVnjP90M9ioiFS38IKaygkAtfMlIJYU5XHz6aC4+fTQA3d3RTuC1OxpYu/MQb9Y2cv/LW+noijZqSgtzOGt8CWeOL+aMcdGvyeXDyPGxPyBdUhKDnVSWikXEgtpZHWQK+JIxsrKMqWOGM3XMcK7/0AQgOgx0054jvFXXyNt1jbxV18hPV7xPe1c3ALmRLKaOKeKMccWcPnY4p42Jfo0pzkv61UA6pSQGG1xTkWbSTOLBU8CXjJaXHeHsylLOriw9vq+jq5ut9U28s/sw7+w+zIbdh3lxUz2/WlN7vExxfjanxU4ep40pYkp5EZPLhzG+pMCzvgGvlmJIhaEE12SvgRT0pYiDSAFfQicnkkX12OFUjx3+gYB04Ggb7+49yrt7j/Du3iO8t/coz7y9m0df6zheJj8ni8mjipgyuojJo4YxadQwqkYNo2pkIaWFuYOqRzqlJIIYXFPVWZ1JFPBFYkYW5TG7KI/ZU0Ye3+eco/5oG1vrm9ha38SW+qNsqT/Kup2H+PWbu+g55qG0MIeqkdHgP3HkMCaWFTJhRAETRxYyZnj+CVcG6ZSSCGpwTZeVVINCo3QkKdKlMzIRrR1d7DzYzLb9TWw/0My2A028H3u8q7HlAyeD3OwsKkcUUDmikIrSAipHFLCnsZXHV+083p8A/o+fh3D87jKZRulISqVTZ2Qi8nMixzuJe2vv7KauoYWdB5vZcbD5+PfaQy28XdfIwab2E14TMWNcST7LN+7jnT2HGV9SwNiSfMYW5zO2JJ9RRXlJX2IiLL+7sFLAF8/53RkZhBZqbnYWk2I5/nia2zupO9RCbUMLtYda2N3Qwu7GVnY1tLBuZwPPvN1yfDjpMZEso7wojzEl+YwtzmP08HxGD8+jfHgeo3tslw3LHfLSE37/7iS5FPDFc352RqZLC7UwN7vPqwOIzinY39TGnsZW9h5uY8/hVvY2tka/H25la30TK7cdpKG544TXmsGIwlxGFeUyqih6QhhVlMfIolxGDstl5LBjj6PfC3Mjx4egpmqGrN8n5LBSwA+hZP/D+dkZmSkt1Kwsi7XY8/st19bZRf2RNvYdaWPf4Tbqj7ax/8gfvu8/2sbaHQ3UH2k74edyTG52FmWFuYwYlktudhZtnd0nlCkblsvaHYcYUZhLaWEOw/NzhpReSpcTcqZSwA+ZVPzD+TmEL52GOnohLztC5YhCKkcUnrRsS3sXB5raOHC0nYNN7ew/2sbBpnYONrdzqKmdg00dtLZ38f6BphPu3HWgqZ1P/Pvvj2+bwfC8bEoKcygtyKWkIIeSghyKC7Ipzs+huCD2lR/dHp6fzfD8HBY+szEjTsjpSgE/ZFLRAvZzCF86DXVMtYLcCJW5Jz85HLsCrGtoYWxxPp/9cBU1VWU0NLdzqKmDxpYOGlo6ONwSe9zcTmNLB7saWzjS2snhlo64Vwn9qWto4TM/fY2ivGyG52czLDebYXnZFOVFvw/Li/zhcW42hbHtwtwIhbnZaXG/hCBQwA+ZVLWA/RofHcQJQunGi99da0cXh1s7ONzSSWNLB0daOzjS2sn/Xfw2DS0n9jvkZWdxqKmdHQebOdLaSVNbJ83t8VNQ8eTnZDEsN5uC3Mjxk0Bh7HFBbjaFOREKcqNfxx7n50QoOLY/J7adGyE/J+v4dvQri9xIVkYsvKeAHzKZ3gIO6gShsHVUHguWo3v1SXd1uwGvutnV7Whu76SprYujbR0cbeuiua2Tpvau4/ub2jo52tZJS0d0X3NbF83tXTS1R08YDc0dtHR00RJ7TUtH1wmjnwYiy6Lps/ycrBO/50TIy45u5+VkkR/7npedRe6x/dlZx79yj31FIscf52VnUZyfw7TxxUP6eQ+UJwHfzOYB9wAR4D7n3MJez+cBPwPOBw4Af+qce9+LY8vghKEFHLTZl+qo/IPBnJAjWcbw/GgH8eK1B054zadmTRxSHTq6umnp6KK1o4vW9m6aOzppae/6w76Oblo7jm13x/Z10db5h8etHd20df7h+5HWTvZ3ttPW2UVbbF9bRzdtXd20DzC9dca4Yp75ykeG9JkGKuGAb2YR4F7gMqAWWGVmS5xzG3oUuxU45Jw71cxuAL4D/Gmix5bBC2oLeCjSpdUc1JFDfv38BntC9vqEmRPJIicSbVGngnOO9q5u2jq7aevopj12Emjr7KK9s/v4V15O8pfp9qKFPxPY7JzbCmBmjwHzgZ4Bfz7wrdjjXwE/MjNzQV3XIcMFrQU8FOnUag7iyKF0+vkF9YQ5UGYWS+tEoP9RtknnxSmlAtjZY7s2ti9uGedcJ9AIjOxVBjNbYGarzWx1fX29B1WTTNVfEAiavvpH/Ow3SaefXxBPmOnKixZ+vK7r3i33gZTBObcIWATRxdMSr5pkKi+DQF+pDa9SHkHsN0mnIJrpAw1SmVrzIuDXAhN6bFcCu/ooU2tm2UAJcNCDY0tIeRUE+kptrN5+kP9aU+dJyiOI/SZBDaLxgl8QT5heSXVqzYuUzipgqplNMrNc4AZgSa8yS4BbYo+vA5Yrfy+JuH1uNQU5kQ/sG0oQ6Cu18ejKnZ6mPK6eUcGKOz7OtoVXsuKOj/uee/bq5+elY8GvrqEFxweD313XTKeitAAjepNyv5eQ9kqqU2sJt/Cdc51mdhuwjOiwzAecc+vN7NvAaufcEuB+4GEz20y0ZX9DoseVcPOq1dxXCqOrj/bIroaWtBkd1J8gXnX0F/yCcJJMhlSn1jwZh++cWwos7bXvb3s8bgU+6cWx5ESZEICGwovRRn2lNiJmcYN+SUFO2oxuOZmgjdZKp34Fr6Q6tZb8gZ+SVH1dBi9eW+d31dJCX6mNG2dNiLvfjLQZ3ZJugjiaKdlSnVpTwE9z6TS8LoiunlERNz/8j1dPj7s/3vrzkNmt0FQJYr9CsvX19xfkUTriozBeBnutr9RGvP3HVpHsLZNboakSxH6FVEhlak0BP80FdXhdpsrkIYJBELR+hUyjlE6aC+NlsJ9SfQku4iW18NOQ9WkAAAYfSURBVNNcWC+D/aRWqKQrBfwMoAAkIgOhlI6ISEiohS8iMgTpOOFRAV9EpA/9raSajjOuFfBFROLoL6in601ZlMMXEYmjv6CerhMe1cIXkdAYTN69v6CerhMe1cIXkVAY7EKD/S3mlq4THhXwRSQUBrvQYH9BPV1nXCulIyKhMNi8+8lmsafjhEcFfBEJhaHk3dMxqPdHKR0RCQW/8+6L19YxZ+FyJt3xNHMWLvflJkVq4YtIKPi50GBQJmop4ItIaPiVognKRC2ldEREkiwoE7USCvhmVmZmz5rZe7HvI+KUOdfMXjGz9Wb2ppn9aSLHFBFJN0G5QXuiLfw7gOedc1OB52PbvTUDNzvnzgTmAf9mZqUJHldEJG343WF8TKIBfz7wUOzxQ8DVvQs45951zr0Xe7wL2AeUJ3hcEZG0EZSJWol22o5xzu0GcM7tNrPR/RU2s5lALrClj+cXAAsAJk6cmGDVRESCIwhj+k8a8M3sOWBsnKf+ZjAHMrNxwMPALc657nhlnHOLgEUANTU1bjDvLyIi/TtpwHfOXdrXc2a218zGxVr344ima+KVKwaeBr7pnHt1yLUVEZEhSzSHvwS4Jfb4FuCp3gXMLBf4f8DPnHO/TPB4IiIyRIkG/IXAZWb2HnBZbBszqzGz+2Jlrgc+CnzGzNbFvs5N8LgiIjJI5lwwU+U1NTVu9erVfldDRCStmNka51xNvOc001ZEJCQU8EVEQkIBX0QkJBTwRURCQgFfRCQkFPBFREJCN0ARkdBbvLbOlzthpZoCvoiEWlBuP5gKSumISKj1d/vBTKOALyKhFpTbD6aCAr6IhFpQbj+YCgr4IhJqQbn9YCqo01ZEQu1Yx6xG6YiIhEAQbj+YCkrpiIiEhAK+iEhIKOCLiISEAr6ISEgo4IuIhIQCvohISCjgi4iEhDnn/K5DXGZWD2xP4C1GAfs9qk460OfNfGH7zPq8Q3OKc6483hOBDfiJMrPVzrkav+uRKvq8mS9sn1mf13tK6YiIhIQCvohISGRywF/kdwVSTJ8384XtM+vzeixjc/giIvJBmdzCFxGRHhTwRURCIuMCvpnNM7NNZrbZzO7wuz7JZmYTzOwFM3vHzNab2Vf8rlMqmFnEzNaa2a/9rkuymVmpmf3KzDbGfs+z/a5TMpnZX8X+lt82s0fNLN/vOnnNzB4ws31m9naPfWVm9qyZvRf7PsLr42ZUwDezCHAvcDkwDbjRzKb5W6uk6wS+7pw7A7gA+GIIPjPAV4B3/K5EitwD/MY5dzpwDhn8uc2sAvgyUOOcOwuIADf4W6ukeBCY12vfHcDzzrmpwPOxbU9lVMAHZgKbnXNbnXPtwGPAfJ/rlFTOud3Ouddjj48QDQYZfeseM6sErgTu87suyWZmxcBHgfsBnHPtzrkGf2uVdNlAgZllA4XALp/r4znn3EvAwV675wMPxR4/BFzt9XEzLeBXADt7bNeS4cGvJzOrAmYAK/2tSdL9G/DXQLffFUmByUA98NNYCus+Mxvmd6WSxTlXB/wLsAPYDTQ6537rb61SZoxzbjdEG3LAaK8PkGkB3+LsC8W4UzMrAv4L+Kpz7rDf9UkWM/tjYJ9zbo3fdUmRbOA84D+cczOAJpJwqR8Usbz1fGASMB4YZmY3+VurzJFpAb8WmNBju5IMvBzszcxyiAb7R5xzT/pdnySbA1xlZu8TTdl93Mx+7m+VkqoWqHXOHbtq+xXRE0CmuhTY5pyrd851AE8CF/pcp1TZa2bjAGLf93l9gEwL+KuAqWY2ycxyiXb2LPG5TkllZkY0v/uOc+57ftcn2ZxzdzrnKp1zVUR/v8udcxnbAnTO7QF2mll1bNclwAYfq5RsO4ALzKww9rd9CRncSd3LEuCW2ONbgKe8PkC212/oJ+dcp5ndBiwj2rv/gHNuvc/VSrY5wKeBt8xsXWzf/3HOLfWxTuKtLwGPxBoxW4E/87k+SeOcW2lmvwJeJzoCbS0ZuMSCmT0KXASMMrNa4O+AhcATZnYr0RPfJz0/rpZWEBEJh0xL6YiISB8U8EVEQkIBX0QkJBTwRURCQgFfRCQkFPBFREJCAV9EJCT+P8dgzGOKVTNOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = list(np.linspace(0, 10))\n",
    "datum = expinfo_list[1]['data']\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(t, cf[:, 1])\n",
    "ax.scatter(datum['t'], datum['cb'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean sum of squares 0.0425931\n",
      "Covariance matrix\n",
      "\n",
      "[[0.00808611, -0.000135221], \n",
      " [-0.000135221, 0.000212259]]\n"
     ]
    }
   ],
   "source": [
    "## calculate covariance matrix\n",
    "\n",
    "# number of estiamted parameters\n",
    "n_est = 2\n",
    "# number of data points\n",
    "n_data = np.sum([expdata[i].shape[0] * len(meas_vars[i]) for i in range(len(expdata))])\n",
    "\n",
    "\n",
    "# hessian\n",
    "H_sol = ssqfn_hess(p_est)\n",
    "# mean sum of squares\n",
    "msq = res['f'] / (n_data - n_est)\n",
    "print(\"mean sum of squares\", msq)\n",
    "\n",
    "print(\"Covariance matrix\")\n",
    "cov = 2 * msq * linalg.inv(H_sol)\n",
    "print(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}