deepmlp_classif.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.io as sio\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import os\n",
    "\n",
    "import iisignature\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from normalizers import *\n",
    "\n",
    "class Path():\n",
    "    def __init__(self, vals, dims, reps=range(1, 7), norm_vals=True):\n",
    "        self.vals = vals.reset_index(drop=True, inplace=False)\n",
    "        self.dims = dims\n",
    "        if norm_vals:\n",
    "            self.vals = better_normalize(self.vals, reps, dims=self.dims)\n",
    "\n",
    "class DB5Path(Path):\n",
    "    def __init__(self, vals, exercise, repetition, dims=None):\n",
    "        if dims is None:\n",
    "            temp_dims = vals.shape[1] - 2\n",
    "        else:\n",
    "            temp_dims = dims\n",
    "        super().__init__(vals, temp_dims)\n",
    "        self.exercise = exercise\n",
    "        self.repetition = repetition\n",
    "    \n",
    "    def naive_padding(self, length) -> np.array:\n",
    "        tempdf = pd.DataFrame(np.random.standard_normal(size = (length, self.dims)))\n",
    "        tempdf['stimulus'] = pd.Series(np.ones(length) * self.exercise)\n",
    "        tempdf['repetition'] = pd.Series(np.ones(length) * self.repetition)\n",
    "        return tempdf\n",
    "    \n",
    "    def get_windows(self, window_size, overlap, use_padding=False) -> np.array:\n",
    "        \"\"\"\n",
    "        returns a list of dataframes\n",
    "        \"\"\"\n",
    "        begin = 0\n",
    "        ans = []\n",
    "        n = len(self.vals)\n",
    "        while(begin < n):\n",
    "            if (begin + window_size < n):\n",
    "                ans.append(self.vals[begin:begin + window_size])\n",
    "            else:\n",
    "                # need to pad it out\n",
    "                if use_padding:\n",
    "                    overflow = (n - begin) % window_size\n",
    "                    ans.append(pd.concat([self.vals[begin:], self.naive_padding(window_size - overflow)], ignore_index=True))\n",
    "                else:\n",
    "                    pass  # other option is to drop the last bit\n",
    "            begin += window_size - overlap\n",
    "        return ans\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shift_stim(stims : pd.Series, shift : int):\n",
    "    # shifts all stimuli by an amount\n",
    "    stims[stims != 0] += shift\n",
    "    return stims\n",
    "\n",
    "def load_subject(data_paths, col='emg'):\n",
    "    df = pd.DataFrame([])\n",
    "    tot_stim = 0\n",
    "    for i, data_path in enumerate(data_paths):\n",
    "        mat = sio.loadmat(data_path)\n",
    "        if i == 0:\n",
    "            df = pd.DataFrame(mat[col])\n",
    "            df['stimulus'] = mat['restimulus']\n",
    "            df['repetition'] = mat['repetition']\n",
    "        else:\n",
    "            df2 = pd.DataFrame(mat[col])\n",
    "            df2['stimulus'] = shift_stim(pd.DataFrame(mat['restimulus']), tot_stim)\n",
    "            df2['repetition'] = mat['repetition']\n",
    "            df = pd.concat([df, df2], ignore_index=True)\n",
    "        tot_stim += len(np.unique(mat['restimulus'])) - 1\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>14</th>\n",
       "      <th>15</th>\n",
       "      <th>16</th>\n",
       "      <th>17</th>\n",
       "      <th>18</th>\n",
       "      <th>19</th>\n",
       "      <th>20</th>\n",
       "      <th>21</th>\n",
       "      <th>stimulus</th>\n",
       "      <th>repetition</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>81.193825</td>\n",
       "      <td>-25.613001</td>\n",
       "      <td>17.336512</td>\n",
       "      <td>-6.904053</td>\n",
       "      <td>12.986491</td>\n",
       "      <td>18.831488</td>\n",
       "      <td>1.902927</td>\n",
       "      <td>19.194145</td>\n",
       "      <td>6.613120</td>\n",
       "      <td>-22.310045</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.972004</td>\n",
       "      <td>27.264063</td>\n",
       "      <td>21.547808</td>\n",
       "      <td>-0.084706</td>\n",
       "      <td>13.291741</td>\n",
       "      <td>65.754509</td>\n",
       "      <td>2.470584</td>\n",
       "      <td>24.564764</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>81.147217</td>\n",
       "      <td>-25.613001</td>\n",
       "      <td>17.336512</td>\n",
       "      <td>-6.949311</td>\n",
       "      <td>12.898676</td>\n",
       "      <td>18.659233</td>\n",
       "      <td>1.944808</td>\n",
       "      <td>19.242105</td>\n",
       "      <td>6.309145</td>\n",
       "      <td>-22.280998</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.901077</td>\n",
       "      <td>27.264063</td>\n",
       "      <td>21.477554</td>\n",
       "      <td>-0.084706</td>\n",
       "      <td>13.291741</td>\n",
       "      <td>65.754509</td>\n",
       "      <td>2.470584</td>\n",
       "      <td>24.564764</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>81.100609</td>\n",
       "      <td>-25.613001</td>\n",
       "      <td>17.336512</td>\n",
       "      <td>-6.994570</td>\n",
       "      <td>12.810861</td>\n",
       "      <td>18.486980</td>\n",
       "      <td>1.986689</td>\n",
       "      <td>19.290066</td>\n",
       "      <td>6.005171</td>\n",
       "      <td>-22.251951</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.830149</td>\n",
       "      <td>27.264063</td>\n",
       "      <td>21.407303</td>\n",
       "      <td>-0.084706</td>\n",
       "      <td>13.291741</td>\n",
       "      <td>65.754509</td>\n",
       "      <td>2.470584</td>\n",
       "      <td>24.564764</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>81.054001</td>\n",
       "      <td>-25.613001</td>\n",
       "      <td>17.336512</td>\n",
       "      <td>-7.039829</td>\n",
       "      <td>12.723045</td>\n",
       "      <td>18.314728</td>\n",
       "      <td>2.028570</td>\n",
       "      <td>19.338028</td>\n",
       "      <td>5.701196</td>\n",
       "      <td>-22.222904</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.759222</td>\n",
       "      <td>27.264063</td>\n",
       "      <td>21.337051</td>\n",
       "      <td>-0.084706</td>\n",
       "      <td>13.291741</td>\n",
       "      <td>65.754509</td>\n",
       "      <td>2.470584</td>\n",
       "      <td>24.564764</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>81.007385</td>\n",
       "      <td>-25.613001</td>\n",
       "      <td>17.336512</td>\n",
       "      <td>-7.085087</td>\n",
       "      <td>12.635230</td>\n",
       "      <td>18.142475</td>\n",
       "      <td>2.070451</td>\n",
       "      <td>19.385988</td>\n",
       "      <td>5.397222</td>\n",
       "      <td>-22.193857</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.688294</td>\n",
       "      <td>27.264063</td>\n",
       "      <td>21.266798</td>\n",
       "      <td>-0.084706</td>\n",
       "      <td>13.291741</td>\n",
       "      <td>65.754509</td>\n",
       "      <td>2.470584</td>\n",
       "      <td>24.564764</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           0          1          2         3          4          5         6  \\\n",
       "0  81.193825 -25.613001  17.336512 -6.904053  12.986491  18.831488  1.902927   \n",
       "1  81.147217 -25.613001  17.336512 -6.949311  12.898676  18.659233  1.944808   \n",
       "2  81.100609 -25.613001  17.336512 -6.994570  12.810861  18.486980  1.986689   \n",
       "3  81.054001 -25.613001  17.336512 -7.039829  12.723045  18.314728  2.028570   \n",
       "4  81.007385 -25.613001  17.336512 -7.085087  12.635230  18.142475  2.070451   \n",
       "\n",
       "           7         8          9  ...        14         15         16  \\\n",
       "0  19.194145  6.613120 -22.310045  ... -3.972004  27.264063  21.547808   \n",
       "1  19.242105  6.309145 -22.280998  ... -3.901077  27.264063  21.477554   \n",
       "2  19.290066  6.005171 -22.251951  ... -3.830149  27.264063  21.407303   \n",
       "3  19.338028  5.701196 -22.222904  ... -3.759222  27.264063  21.337051   \n",
       "4  19.385988  5.397222 -22.193857  ... -3.688294  27.264063  21.266798   \n",
       "\n",
       "         17         18         19        20         21  stimulus  repetition  \n",
       "0 -0.084706  13.291741  65.754509  2.470584  24.564764         0           0  \n",
       "1 -0.084706  13.291741  65.754509  2.470584  24.564764         0           0  \n",
       "2 -0.084706  13.291741  65.754509  2.470584  24.564764         0           0  \n",
       "3 -0.084706  13.291741  65.754509  2.470584  24.564764         0           0  \n",
       "4 -0.084706  13.291741  65.754509  2.470584  24.564764         0           0  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subjects = range(1, 11)\n",
    "exercises = [1]\n",
    "root = './signatures_data/'\n",
    "df = pd.DataFrame([])\n",
    "for subject in subjects:\n",
    "    data_paths = [os.path.join(root, f\"s{subject}/S{subject}_E{exercise}_A1.mat\") for exercise in exercises]\n",
    "    cur_df = load_subject(data_paths, col='glove')\n",
    "    df = pd.concat([df, cur_df], ignore_index=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_log_signature(path, depth, dims=None):\n",
    "    if dims is None:\n",
    "        dims = path.shape[1] # take second, assumes (timesteps, dims)\n",
    "    s = iisignature.prepare(dims, depth)\n",
    "    return iisignature.logsig(path, s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "repetitions = [1, 2, 3, 4, 5, 6]\n",
    "stims = [1, 3]\n",
    "train_set = [1, 3, 4, 6]\n",
    "test_set = [2, 5]\n",
    "train_windows = []\n",
    "test_windows = []\n",
    "for stimulus in stims:\n",
    "    for repetition in repetitions:\n",
    "        dbpth = DB5Path(df[(df.stimulus == stimulus) & (df.repetition == repetition)], stimulus, repetition)\n",
    "        if repetition in train_set:\n",
    "            train_windows += dbpth.get_windows(window_size=200, overlap=100)\n",
    "        elif repetition in test_set:\n",
    "            test_windows += dbpth.get_windows(window_size=200, overlap=100)\n",
    "train_x, train_y = splitxy(train_windows, [5,8,10], ['stimulus'])\n",
    "test_x, test_y = splitxy(test_windows, [5,8,10], ['stimulus'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(692, 336, (692, 200, 3), (692,))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(train_windows), len(test_windows), train_x.shape, train_y.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 692/692 [00:58<00:00, 11.77it/s]\n",
      "100%|██████████| 336/336 [00:28<00:00, 11.83it/s]\n"
     ]
    }
   ],
   "source": [
    "def splitxy(dfs : list, x_cols, y_col):\n",
    "    x = []\n",
    "    y = []\n",
    "    for df in dfs:\n",
    "        x.append(df[x_cols].to_numpy())\n",
    "        y.append(df[y_col].to_numpy()[0]) # assume everything is of the same class\n",
    "    return np.array(x), np.array(y).reshape(-1)\n",
    "\n",
    "train_x, train_y = splitxy(train_windows, range(16), ['stimulus'])\n",
    "test_x, test_y = splitxy(test_windows, range(16), ['stimulus'])\n",
    "\n",
    "train_signatures = []\n",
    "for window in tqdm(train_x):\n",
    "    train_signatures.append(compute_log_signature(window, depth=4))\n",
    "test_signatures = []\n",
    "for window in tqdm(test_x):\n",
    "    test_signatures.append(compute_log_signature(window, depth=4))\n",
    "train_signatures = np.array(train_signatures)\n",
    "test_signatures = np.array(test_signatures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(pred, true):\n",
    "    assert len(pred) == len(true), \"check lengths match\"\n",
    "    return np.sum((pred == true)) / len(pred)\n",
    "\n",
    "assert accuracy(np.array([1, 1, 1, 2]), np.array([1, 2, 1, 2])) == 0.75"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(692, 17816) (692,)\n",
      "(336, 17816) (336,)\n"
     ]
    }
   ],
   "source": [
    "print(train_signatures.shape, train_y.shape)\n",
    "print(test_signatures.shape, test_y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.15287352, -0.15799403,  0.0984855 ,  0.068694  ,  0.80948655,\n",
       "         1.84939927,  0.39183229,  0.05124822, -0.09199243,  0.33404976]),\n",
       " 1)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_signatures[0, :10], train_y[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-16 {color: black;background-color: white;}#sk-container-id-16 pre{padding: 0;}#sk-container-id-16 div.sk-toggleable {background-color: white;}#sk-container-id-16 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-16 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-16 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-16 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-16 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-16 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-16 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-16 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-16 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-16 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-16 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-16 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-16 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-16 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-16 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-16 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-16 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-16 div.sk-item {position: relative;z-index: 1;}#sk-container-id-16 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-16 div.sk-item::before, #sk-container-id-16 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-16 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-16 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-16 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-16 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-16 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-16 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-16 div.sk-label-container {text-align: center;}#sk-container-id-16 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-16 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-16\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>MLPClassifier(alpha=0.001, early_stopping=True, hidden_layer_sizes=[50, 50],\n",
       "              max_iter=2000, random_state=1, solver=&#x27;lbfgs&#x27;, warm_start=True)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-16\" type=\"checkbox\" checked><label for=\"sk-estimator-id-16\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">MLPClassifier</label><div class=\"sk-toggleable__content\"><pre>MLPClassifier(alpha=0.001, early_stopping=True, hidden_layer_sizes=[50, 50],\n",
       "              max_iter=2000, random_state=1, solver=&#x27;lbfgs&#x27;, warm_start=True)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "MLPClassifier(alpha=0.001, early_stopping=True, hidden_layer_sizes=[50, 50],\n",
       "              max_iter=2000, random_state=1, solver='lbfgs', warm_start=True)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mlp = MLPClassifier(                \n",
    "                solver=\"lbfgs\",\n",
    "                random_state=1,\n",
    "                max_iter=2000,\n",
    "                early_stopping=True,\n",
    "                alpha=1e-3,\n",
    "                hidden_layer_sizes=[50, 50,], warm_start=True)\n",
    "\n",
    "mlp.fit(train_signatures, train_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training set accuracy: 0.9985549132947977\n"
     ]
    }
   ],
   "source": [
    "train_pred = mlp.predict(train_signatures)\n",
    "print(\"training set accuracy:\", accuracy(train_pred, train_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8869047619047619"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred = mlp.predict(test_signatures)\n",
    "accuracy(pred, test_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAGVCAYAAACIDTRmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzkUlEQVR4nO3de1hU5do/8O8agQGBGcSCYRQUw/M5LaOsPJCobcOt1c+iwkO6MzUPeXwLUzxQVma4Tc3Mw341s4OU7sLIs4kkKG41RFFUEsHeTYBgwDBr/f5gMzUbsBlnwThrfT/Xta6rWae5KeP2vp9nPUuQJEkCERGRQmmcHQAREVFDYqIjIiJFY6IjIiJFY6IjIiJFY6IjIiJFY6IjIiJFY6IjIiJFc3N2AERE1LDKy8tRWVkpy708PDzg6ekpy70aCxMdEZGClZeXI7SVD/Kvm2W5n8FgQE5OjkslOyY6IiIFq6ysRP51My6nt4bO17HRqpIbIlr1uoTKykomOiIiurP4+Arw8RUcuocIx653FiY6IiIVMEsizA6ubGyWRHmCaWScdUlERIrGio6ISAVESBDhWEnn6PXOwkRHRKQCIkQ42nh0/A7OwdYlEREpGis6IiIVMEsSzA6+Z9vR652FiY6ISAXUPEbH1iURESkaKzoiIhUQIcGs0oqOiY6ISAXYuiQiIlIoVnRERCrAWZdERKRo4n82R+/hiti6JCIiRWNFR0SkAmYZZl06er2zMNEREamAWYIMr+mRJ5bGxtYlEREpGis6IiIVUPNkFCY6IiIVECHADMHhe7giti6JiEjRWNEREamAKFVvjt7DFTHRERGpgFmG1qWj1zsLW5dERKRorOiIiFRAzRUdEx0RkQqIkgBRcnDWpYPXOwtbl0REpGis6IiIVICtSyIiUjQzNDA72MQzyxRLY2PrkoiIFI0VHRGRCkgyTEaRXHQyChMdEZEKqHmMjq1LIiJSNFZ0REQqYJY0MEsOTkbhWpdERHSnEiFAdLCJJ8I1M51LJzpRFJGXlwdfX18Igmv2jomI/kiSJNy4cQNGoxEaDUeX5ODSiS4vLw/BwcHODoOISHa5ublo2bKlbPdT82QUl050vr6+AIDLx1tD58O/+VDD+mu7rs4OgVSgCiYcxjeW329ykWeMjq3LRlfTrtT5aKDzZaKjhuUmuDs7BFKD/+QSDsfIx6UTHRER2aZ6MoqDby9g65KIiO5UogxrXbrqrEv2+4iIqEEcPHgQw4YNg9FohCAISExMrPfcl156CYIgYMWKFVb7CwsLER0dDZ1OBz8/P4wbNw6lpaV2xcFER0SkAjWTURzd7FFWVobu3btj1apVtzxvx44dOHr0KIxGY61j0dHROHPmDJKTk7Fr1y4cPHgQEyZMsCsOti6JiFRAhKbRHxgfMmQIhgwZcstzrl69iilTpmD37t14/PHHrY5lZmYiKSkJx44dQ+/evQEAK1euxNChQ/HOO+/UmRjrwoqOiIjsUlJSYrVVVFTc1n1EUcTzzz+PWbNmoXPnzrWOp6SkwM/Pz5LkACAiIgIajQapqak2fw8THRGRCpglQZYNAIKDg6HX6y1bfHz8bcX01ltvwc3NDa+88kqdx/Pz8xEQEGC1z83NDf7+/sjPz7f5e9i6JCJSAXneMF7duszNzYVOp7Ps12q1dt8rPT0d77//Po4fP97gzwyyoiMiIrvodDqr7XYS3aFDh3D9+nWEhITAzc0Nbm5uuHz5Ml599VW0bt0aAGAwGHD9+nWr66qqqlBYWAiDwWDzd7GiIyJSAVHSQHRwCTBRxiXAnn/+eURERFjti4yMxPPPP48xY8YAAMLDw1FUVIT09HT06tULALB3716Ioog+ffrY/F1MdEREKiBn69JWpaWlyM7OtnzOyclBRkYG/P39ERISgubNm1ud7+7uDoPBgPbt2wMAOnbsiMGDB2P8+PFYs2YNTCYTJk+ejFGjRtk84xJg65KIiBpIWloaevbsiZ49ewIAZsyYgZ49e2L+/Pk232PLli3o0KEDBg4ciKFDh6Jv37748MMP7YqDFR0RkQqIgGXWpCP3sEe/fv0g2dHuvHTpUq19/v7+2Lp1q53fbI2JjohIBeR5YNw1m4CuGTUREZGNWNEREamAPC9edc3aiImOiEgF1Pw+OtdMz0RERDZiRUdEpAJqbl26ZtREREQ2YkVHRKQC8qyM4pq1ERMdEZEKiJIA0dEHxh283llcMz0TERHZiBUdEZEKiDK0Ll11ZRQmOiIiFZDnNT2umehcM2oiIiIbsaIjIlIBMwSYHVzZxNHrnYWJjohIBdi6JCIiUihWdEREKmCG461HszyhNDomOiIiFWDrkoiISKFY0RERqYCa317AREdEpAKSDC9elVz08QLXTM9EREQ2YkVHRKQCbF0SEZGi8TU9RERECsWKjohIBfiGcSIiUjS2LomIiBSKFR0RkQqI0Dj8hnC+YZyIiO5YZkmA2cHWo6PXO4trpmciIiIbsaIjIlIBNU9GYaIjIlIBSYbX9EguujKKa0ZNRERkI1Z0REQqYIYgwxvG2bokIqI7lCg5PsYmSjIF08jYuiQiIkVjRUdEpAKiDJNRHL3eWZjoiIhUQJThDeOOXu8srpmeiYiIbMSKjohIBdS8BBgTHRGRCqh5jM41oyYiIrIRKzoiIhUQIcNal5yMQkREdyrpP7MuHdkkOxPdwYMHMWzYMBiNRgiCgMTERMsxk8mEOXPmoGvXrvD29obRaMQLL7yAvLw8q3sUFhYiOjoaOp0Ofn5+GDduHEpLS+2Kg4mOiIgaRFlZGbp3745Vq1bVOnbz5k0cP34csbGxOH78OL788ktkZWXhiSeesDovOjoaZ86cQXJyMnbt2oWDBw9iwoQJdsXB1iURkQrI+ZqekpISq/1arRZarbbW+UOGDMGQIUPqvJder0dycrLVvr///e+4//77ceXKFYSEhCAzMxNJSUk4duwYevfuDQBYuXIlhg4dinfeeQdGo9GmuFnRERGpQM2sS0c3AAgODoZer7ds8fHxssRYXFwMQRDg5+cHAEhJSYGfn58lyQFAREQENBoNUlNTbb4vKzoiIrJLbm4udDqd5XNd1Zy9ysvLMWfOHDzzzDOWe+fn5yMgIMDqPDc3N/j7+yM/P9/mezPRERGpgJytS51OZ5XoHGUymfD0009DkiSsXr1atvvWYKIjIlKBO3Wty5okd/nyZezdu9cqgRoMBly/ft3q/KqqKhQWFsJgMNj8HRyjIyIip6hJcufPn8f333+P5s2bWx0PDw9HUVER0tPTLfv27t0LURTRp08fm7+HFR0RkQrI2bq0VWlpKbKzsy2fc3JykJGRAX9/fwQFBeHJJ5/E8ePHsWvXLpjNZsu4m7+/Pzw8PNCxY0cMHjwY48ePx5o1a2AymTB58mSMGjXK5hmXABMdEZEqOCPRpaWloX///pbPM2bMAADExMRgwYIF+PrrrwEAPXr0sLpu37596NevHwBgy5YtmDx5MgYOHAiNRoORI0ciISHBrjiY6IiIqEH069cPkiTVe/xWx2r4+/tj69atDsXBREdEpALOqOjuFEx0CnHqqDc++yAA5081RWGBO95Yn4MHhxRbjr8zLQTJ2/2trunVrwRLt14EAJw84oPZT4bVee+Eb7LQvsdvDRc8ubQufUrx1Mu/oG3Xm2huqMKCsa2RkqS3Oic4rBzjXr+Gbg+UookbcPmcFovGt8YvVz2cFLX6MNE5ycGDB/H2228jPT0d165dw44dOzB8+HBnhuSyym9q0Kbzb4h8phBx40LrPKd3/xK8+t4Vy2d3j9/bBp16l+GTjNNW529aFoSMwz5o151Jjurn2VTExTOe2P2JP974+FKt40GtKrA8MRtJ2/zxj3cCcfNGE7RqX47Kctf8pUmux6mJrmbBz7Fjx2LEiBHODMXl3TfgBu4bcOOW57h7SPAPqLLpWJUJSNmtQ9TY/4PA30d0C2n7dEjbV//Dw6Pn5uPHvTqsX/z7LLlrlx1fSYPsI8Hx5+D+fETtzuTURHerBT9Jfv9K8cHTXTvDV29G976lGD37GnT+5jrPTflOjxu/umHQ/yts5ChJSQRBwv0DS/DZBwFYsvUCwrqUI/+KB7b9PaBWe5Malppbly71wHhFRQVKSkqsNrJN734lmPX+Zby1/QLGvXYNp1J88NpzbWCuO89h9yfN0avfDdxtNDVuoKQofndVoamPiP83+TrS9ukw75k2+CFJh/kfXULXB+x7pxjR7XKpySjx8fFYuHChs8NwSf2GF1n+ObRjOUI7/YbR4Z3wryM+6Pmw9S+cX/Lckb7fF/+z9lLjBkmKI/znr9Ipu3XYse5uAMDFM17o1PsmHn/h3zh11MeJ0akLKzoXMW/ePBQXF1u23NxcZ4fksoJaVULvX4W8S7XHSr771B++zaoQPqi4jiuJbFdS2ARVJuDyOU+r/bnntQhoUemkqNSpJtE5urkil6ro6nu5H9nvlzx3lPzaBP4B1q1JSapOdBFP/go3dycFR4pRZdLg3MmmaHlPhdX+Fm0qcP1nPlpAjcOlEh3V77cyDfJyfv9LQH6uBy6c9oKvXxV8m5nxv+8a0PfxIjQLqMK1Sx74aLERxtAK9OpnPVMz47AP8q9oMfjZfzf2j0AuyrOpGcbQ36szQ3Al2nT+DTeKmuCXqx747IMA/M+ayzh91Bsnj/igd/8beOCxEsx68h4nRq0+am5dOjXR3WrBz5CQECdG5nrOnWxq9cD32gUtAACPPV2IKfG5yMn0RPJnoSgraYLmgVW499ESxMzOh4fWesJw0ifN0al3KULaWv8NnKg+7br/hre/uGD5/NLCPADAd582w7vTQ3AkSY+EuS0wavJ1TFx0FT9frH5Y/MyPHJ9rTJIkQHIwUTl6vbMIki2LjTWQ/fv3Wy34WSMmJgYbN2780+tLSkqg1+vx67k20Pm61HAjuaBIYw9nh0AqUCWZsB9fobi4WJaXm9b8nnzoq8lw83Zs6KeqrAI/RP1dttgai1Mruj9b8JOIiORxp754tTFwjI6ISAXUPEbHfh8RESkaKzoiIhVQ82QUJjoiIhVg65KIiEihWNEREakAW5dERKRokgytS1dNdGxdEhGRorGiIyJSAQnVi7Y7eg9XxERHRKQCIgQIKl0Zha1LIiJSNFZ0REQqwFmXRESkaKIkQOAD40RERMrDio6ISAUkSYZZly467ZKJjohIBdQ8RsfWJRERKRorOiIiFVBzRcdER0SkApx1SUREpFCs6IiIVICzLomISNGqE52jY3QyBdPI2LokIiJFY0VHRKQCnHVJRESKJsHx98m5aOeSrUsiIlI2VnRERCrA1iURESmbinuXbF0SEZGisaIjIlIDGVqXcNHWJSs6IiIVqFkZxdHNHgcPHsSwYcNgNBohCAISExP/KyYJ8+fPR1BQELy8vBAREYHz589bnVNYWIjo6GjodDr4+flh3LhxKC0ttSsOJjoiImoQZWVl6N69O1atWlXn8WXLliEhIQFr1qxBamoqvL29ERkZifLycss50dHROHPmDJKTk7Fr1y4cPHgQEyZMsCsOti6JiFTAGbMuhwwZgiFDhtRzLwkrVqzA66+/jqioKADA5s2bERgYiMTERIwaNQqZmZlISkrCsWPH0Lt3bwDAypUrMXToULzzzjswGo02xcGKjohIDSRBng1ASUmJ1VZRUWF3ODk5OcjPz0dERIRln16vR58+fZCSkgIASElJgZ+fnyXJAUBERAQ0Gg1SU1Nt/i4mOiIisktwcDD0er1li4+Pt/se+fn5AIDAwECr/YGBgZZj+fn5CAgIsDru5uYGf39/yzm2YOuSiEgF5HxNT25uLnQ6nWW/Vqt17MYNjBUdEZEaSDJtAHQ6ndV2O4nOYDAAAAoKCqz2FxQUWI4ZDAZcv37d6nhVVRUKCwst59iCiY6IiBpdaGgoDAYD9uzZY9lXUlKC1NRUhIeHAwDCw8NRVFSE9PR0yzl79+6FKIro06ePzd/F1iURkQo4Y9ZlaWkpsrOzLZ9zcnKQkZEBf39/hISEYNq0aVi8eDHatm2L0NBQxMbGwmg0Yvjw4QCAjh07YvDgwRg/fjzWrFkDk8mEyZMnY9SoUTbPuARsTHRff/21zTd84oknbD6XiIgaUSOvVZmWlob+/ftbPs+YMQMAEBMTg40bN2L27NkoKyvDhAkTUFRUhL59+yIpKQmenp6Wa7Zs2YLJkydj4MCB0Gg0GDlyJBISEuyKQ5CkPx+e1Ghs63AKggCz2WxXAI4oKSmBXq/Hr+faQOfLLiw1rEhjD2eHQCpQJZmwH1+huLjYasLH7ar5PRny4XxovDz//IJbEH8rx5UJcbLF1lhsquhEUWzoOIiIqAGp+TU9DpVBf1ymhYiI7mAyzrp0NXYnOrPZjEWLFqFFixbw8fHBxYsXAQCxsbFYv3697AESERE5wu5Et2TJEmzcuBHLli2Dh4eHZX+XLl3w0UcfyRocERHJRZBpcz12J7rNmzfjww8/RHR0NJo0aWLZ3717d5w9e1bW4IiISCZsXdru6tWrCAsLq7VfFEWYTCZZgiIiIpKL3YmuU6dOOHToUK39n3/+OXr27ClLUEREJDMVV3R2r4wyf/58xMTE4OrVqxBFEV9++SWysrKwefNm7Nq1qyFiJCIiR/3hNTsO3cMF2V3RRUVFYefOnfj+++/h7e2N+fPnIzMzEzt37sRjjz3WEDESERHdttta6/Lhhx9GcnKy3LEQEVEDkfM1Pa7mthd1TktLQ2ZmJoDqcbtevXrJFhQREclMjjE2tSS6n3/+Gc888wx++OEH+Pn5AQCKiorw4IMPYtu2bWjZsqXcMRIREd02u8foXnzxRZhMJmRmZqKwsBCFhYXIzMyEKIp48cUXGyJGIiJyVM1kFEc3F2R3RXfgwAEcOXIE7du3t+xr3749Vq5ciYcffljW4IiISB6CVL05eg9XZHdFFxwcXOeD4Waz2a4X4RERETUGuxPd22+/jSlTpiAtLc2yLy0tDVOnTsU777wja3BERCQTPjB+a82aNYMg/N6bLSsrQ58+feDmVn15VVUV3NzcMHbsWMsr0ImI6A6i4gfGbUp0K1asaOAwiIiIGoZNiS4mJqah4yAioobE5+huT3l5OSorK6326XQ6hwIiIqIGoOJEZ/dklLKyMkyePBkBAQHw9vZGs2bNrDYiIqI7id2Jbvbs2di7dy9Wr14NrVaLjz76CAsXLoTRaMTmzZsbIkYiInIUZ13abufOndi8eTP69euHMWPG4OGHH0ZYWBhatWqFLVu2IDo6uiHiJCIiR6h41qXdFV1hYSHatGkDoHo8rrCwEADQt29fHDx4UN7oiIiIHGR3omvTpg1ycnIAAB06dMD27dsBVFd6NYs8ExHRnaVmCTBHN1dkd6IbM2YMTp48CQCYO3cuVq1aBU9PT0yfPh2zZs2SPUAiIpIBx+hsN336dMs/R0RE4OzZs0hPT0dYWBi6desma3BERESOcug5OgBo1aoVWrVqJUcsREREsrMp0SUkJNh8w1deeeW2gyEiooYhQIbX9MgSSeOzKdG99957Nt1MEASnJLoRXXvDTXBv9O8ldZl/McXZIZAKlN0QsZ+jQLKyKdHVzLIkIiIXpeLn6BweoyMiIhfAtS6JiIiUiRUdEZEaqLiiY6IjIlIBOVY2Uc3KKERERK7kthLdoUOH8NxzzyE8PBxXr14FAPzjH//A4cOHZQ2OiIhkouIlwOxOdF988QUiIyPh5eWFEydOoKKiAgBQXFyMpUuXyh4gERHJgInOdosXL8aaNWuwbt06uLv//pD2Qw89hOPHj8saHBERkaPsnoySlZWFRx55pNZ+vV6PoqIiOWIiIiKZcTKKHQwGA7Kzs2vtP3z4sOWFrEREdIepWRnF0c0F2Z3oxo8fj6lTpyI1NRWCICAvLw9btmzBzJkzMXHixIaIkYiI6LbZ3bqcO3cuRFHEwIEDcfPmTTzyyCPQarWYOXMmpkyZ0hAxEhGRo/jAuO0EQcBrr72GWbNmITs7G6WlpejUqRN8fHwaIj4iIpIBx+hug4eHBzp16oT777+fSY6IiGoxm82IjY1FaGgovLy8cM8992DRokWQpN8zpiRJmD9/PoKCguDl5YWIiAicP39e1jjsruj69+8PQah/QHLv3r0OBURERA3ACa3Lt956C6tXr8amTZvQuXNnpKWlYcyYMdDr9ZZ3ly5btgwJCQnYtGkTQkNDERsbi8jISPz000/w9PR0MOBqdie6Hj16WH02mUzIyMjA6dOnERMTI0tQREQkMxlal/YmuiNHjiAqKgqPP/44AKB169b45JNP8OOPP1bfTpKwYsUKvP7664iKigIAbN68GYGBgUhMTMSoUaMcDLia3YmuvreNL1iwAKWlpQ4HREREd7aSkhKrz1qtFlqtttZ5Dz74ID788EOcO3cO7dq1w8mTJ3H48GEsX74cQPVLvfPz8xEREWG5Rq/Xo0+fPkhJSZEt0cm2qPNzzz2Hjz/+WK7bERGRnGRcAiw4OBh6vd6yxcfH1/mVc+fOxahRo9ChQwe4u7ujZ8+emDZtGqKjowEA+fn5AIDAwECr6wIDAy3H5CDba3pSUlJk66cSEZHMZByjy83NhU6ns+yuq5oDgO3bt2PLli3YunUrOnfujIyMDEybNg1Go7FRh7rsTnQjRoyw+ixJEq5du4a0tDTExsbKFhgREd2ZdDqdVaKrz6xZsyxVHQB07doVly9fRnx8PGJiYmAwGAAABQUFCAoKslxXUFBQaz6II+xOdHq93uqzRqNB+/btERcXh0GDBskWGBERyccZz9HdvHkTGo31CFmTJk0giiIAIDQ0FAaDAXv27LEktpKSEqSmpsq60pZdic5sNmPMmDHo2rUrmjVrJlsQRESkPMOGDcOSJUsQEhKCzp0748SJE1i+fDnGjh0LoHoBkmnTpmHx4sVo27at5fECo9GI4cOHyxaHXYmuSZMmGDRoEDIzM5noiIjollauXInY2Fi8/PLLuH79OoxGI/72t79h/vz5lnNmz56NsrIyTJgwAUVFRejbty+SkpJknfNhd+uyS5cuuHjxIkJDQ2ULgoiIGpgTHhj39fXFihUrsGLFinrPEQQBcXFxiIuLcyy2W7itF6/OnDkTu3btwrVr11BSUmK1ERHRnadmjM7RzRXZXNHFxcXh1VdfxdChQwEATzzxhNVSYJIkQRAEmM1m+aMkIiK6TTYnuoULF+Kll17Cvn37GjIeIiJqKC5akTnK5kRXs9r0o48+2mDBEBFRA1Hx++jsGqO71VsLiIiI7kR2zbps167dnya7wsJChwIiIiL5qfnFq3YluoULF9ZaGYWIiFyAiluXdiW6UaNGISAgoKFiISIikp3NiY7jc0REroutSxvUzLokIiIXxNbln6tZbZqIiMiVyPbiVSIiuoOxoiMiIiVT8xid3Ys6ExERuRJWdEREasDWJRERKZqKEx1bl0REpGis6IiIVEDNk1GY6IiI1ICtSyIiImViRUdEpAJsXRIRkbKxdUlERKRMrOiIiNRAxRUdEx0RkQoI/9kcvYcrYuuSiIgUjRUdEZEasHVJRERKpubHC9i6JCIiRWNFR0SkBmxdEhGR4rloonIUW5dERKRorOiIiFRAzZNRmOiIiNRAxWN0bF0SEZGisaIjIlIBNbcuWdEREZGisaIjIlIDFY/RMdEREakAW5dEREQKxYqOiEgN2LokIiJFU3GiY+uSiIgUjYmOiEgFaiajOLrZ6+rVq3juuefQvHlzeHl5oWvXrkhLS7MclyQJ8+fPR1BQELy8vBAREYHz58/L+JMz0RERqYMk02aHX3/9FQ899BDc3d3x7bff4qeffsK7776LZs2aWc5ZtmwZEhISsGbNGqSmpsLb2xuRkZEoLy937Of9A47RERFRg3jrrbcQHByMDRs2WPaFhoZa/lmSJKxYsQKvv/46oqKiAACbN29GYGAgEhMTMWrUKFniYEVHRKQCgiTJsgFASUmJ1VZRUVHnd3799dfo3bs3nnrqKQQEBKBnz55Yt26d5XhOTg7y8/MRERFh2afX69GnTx+kpKTI9rMz0SlUl/tLsOCjc9hy9ASScn5E+GO/Wh33u8uEV9++iC1HTyDxpzQs3pgFY2v5WgWkXJd/9MEnL96D5Q90QVybe3H2O3295/7ztWDEtbkXRz++u87jVRUC1j7eAXFt7kX+T14NFTIBsrYug4ODodfrLVt8fHydX3nx4kWsXr0abdu2xe7duzFx4kS88sor2LRpEwAgPz8fABAYGGh1XWBgoOWYHNi6VChPLxE5mU3x3fa7MH9t9n8dlfDG2nOoqtJg4YS2uFnaBCPG5SP+f89iwmNdUfFbE6fETK6h8qYGgR1voudT/4ftE++p97yzu/X4OcMbvoGV9Z7z/Vst4BtgQkFmQ0RKDSU3Nxc6nc7yWavV1nmeKIro3bs3li5dCgDo2bMnTp8+jTVr1iAmJqZRYgWcXNGtXr0a3bp1g06ng06nQ3h4OL799ltnhqQYaQf8sOndljjynX+tYy1Cy9Hx3jL8/fVWOPcvH/x80QsrX28NrVZE/yf+7YRoyZW07VeCAa9eQ4fI4nrPKcl3x7cLg/HX9y5B41b3DIbz+3W4eEiHx/7nakOFSn8g56zLmt/ZNVt9iS4oKAidOnWy2texY0dcuXIFAGAwGAAABQUFVucUFBRYjsnBqYmuZcuWePPNN5Geno60tDQMGDAAUVFROHPmjDPDUjx3j+o/rZUVv//nlyQBpkoNOvcudVZYpBCSCCS+2hoPji9AQLu62+Glv7hh1/+EYPi7l+DuJTZyhCrlhFmXDz30ELKysqz2nTt3Dq1atQJQPTHFYDBgz549luMlJSVITU1FeHi4vT9hvZya6IYNG4ahQ4eibdu2aNeuHZYsWQIfHx8cPXq0zvMrKipqDYKS/XIveKLgqgfGzP4ZProquLmLeOpvebjbWAn/gPrbTES2+GFNIDRNJNw/+pc6j0sS8NXsVuj17P/B2O1mI0dHjWn69Ok4evQoli5diuzsbGzduhUffvghJk2aBAAQBAHTpk3D4sWL8fXXX+PUqVN44YUXYDQaMXz4cNniuGPG6MxmMz777DOUlZXVm8nj4+OxcOHCRo5MecxVGix6qS2mv5WDz08eh7kKOPGDHj/u00MQnB0dubK8U15I3RiACTvP1vtn6cdNd6OytAn6TpRvsgH9OWe8veC+++7Djh07MG/ePMTFxSE0NBQrVqxAdHS05ZzZs2ejrKwMEyZMQFFREfr27YukpCR4eno6FuwfOD3RnTp1CuHh4SgvL4ePjw927NhRq6dbY968eZgxY4blc0lJCYKDgxsrVEXJPu2NSY93QVPfKri7SygudMeKHWdw/pS3s0MjF3blmA/K/u2GFX27WPZJZgHJS1sidUMAph46g0spvvj5hDeWdOhpde26qA7oGlWI4e9cbuyw1cFJa13+5S9/wV/+8pd6jwuCgLi4OMTFxTkQ2K05PdG1b98eGRkZKC4uxueff46YmBgcOHCgzmSn1WrrHfSk23PzRvUfAWPrcrTtWobNy1s6OSJyZd3+Wog2D92w2rdldBi6Di9Ej6eqJzoNnp+L/jPyLMdvXHfHlpi2eDIhBy16lDVqvKQOTk90Hh4eCAsLAwD06tULx44dw/vvv4+1a9c6OTLX5tnUDGOr3ycCGIIr0KZjGW4Uu+GXPC0eHlqI4n+74XqeB1p3+A0T519GynfNcPxQ/c9EEQFAZZkGhZd//wtnUa4W+T95wUtfBX0LE5o2M1udr3GT4HO3CXe1qX6oWN/CBMBkOe7hXT0ZpVmrCuiCTKCGoeYXrzo90f03URTrfcqebNeuaxmWbTtr+fy32OrpvMmf34V3Z7WBf0AlJrx2BX53mVD4izv2fHkXtq40OitcciF5p5pi87PtLJ+/W1LdBeg+8t+IepttxzuWil/T49REN2/ePAwZMgQhISG4ceMGtm7div3792P37t3ODEsR/pWqw+DQ++s9/tVGA77aKN9zKqQerR8oxfyLx20+f+qhWz8u5Ney0q77EdnLqYnu+vXreOGFF3Dt2jXo9Xp069YNu3fvxmOPPebMsIiIFMlVW4+OcmqiW79+vTO/nohIPSSpenP0Hi6IizoTEZGi3XGTUYiISH6cdUlERMqm4lmXbF0SEZGisaIjIlIBQazeHL2HK2KiIyJSA7YuiYiIlIkVHRGRCnDWJRERKRsfGCciIlImVnRERCrA1iURESkbZ10SEREpEys6IiIVYOuSiIiUjbMuiYiIlIkVHRGRCrB1SUREysZZl0RERMrEio6ISAXYuiQiImUTperN0Xu4ILYuiYhI0VjRERGpgYonozDRERGpgAAZxuhkiaTxsXVJRESKxoqOiEgNVLwEGBMdEZEKqPnxArYuiYhI0VjRERGpAWddEhGRkgmSBMHBMTZHr3cWti6JiEjRWNEREamB+J/N0Xu4ICY6IiIVYOuSiIhIoVjRERGpAWddEhGRoql4ZRS2LomISNFY0RERqQCXACMiImWraV06ut2mN998E4IgYNq0aZZ95eXlmDRpEpo3bw4fHx+MHDkSBQUFMvyw1pjoiIioQR07dgxr165Ft27drPZPnz4dO3fuxGeffYYDBw4gLy8PI0aMkP37meiIiFRAEOXZ7FVaWoro6GisW7cOzZo1s+wvLi7G+vXrsXz5cgwYMAC9evXChg0bcOTIERw9elTGn5yJjohIHWRsXZaUlFhtFRUV9X7tpEmT8PjjjyMiIsJqf3p6Okwmk9X+Dh06ICQkBCkpKbL+6Ex0RERkl+DgYOj1essWHx9f53nbtm3D8ePH6zyen58PDw8P+Pn5We0PDAxEfn6+rPFy1iURkRrI+MB4bm4udDqdZbdWq611am5uLqZOnYrk5GR4eno6+MWOYUVHRKQCNWtdOroBgE6ns9rqSnTp6em4fv067r33Xri5ucHNzQ0HDhxAQkIC3NzcEBgYiMrKShQVFVldV1BQAIPBIOvPzoqOiIhkN3DgQJw6dcpq35gxY9ChQwfMmTMHwcHBcHd3x549ezBy5EgAQFZWFq5cuYLw8HBZY2GiIyJSg0ZeAszX1xddunSx2uft7Y3mzZtb9o8bNw4zZsyAv78/dDodpkyZgvDwcDzwwAOOxflfmOiIiNRAguPvk5N5ZZT33nsPGo0GI0eOREVFBSIjI/HBBx/I+yVgoiMiokayf/9+q8+enp5YtWoVVq1a1aDfy0RHRKQCan7xKhMdEZEaSJBhjE6WSBodHy8gIiJFY0VHRKQGKn7xKhMdEZEaiAAEGe7hgti6JCIiRWNFR0SkApx1SUREyqbiMTq2LomISNFY0RERqYGKKzomOiIiNVBxomPrkoiIFI0VHRGRGqj4OTomOiIiFVDz4wVsXRIRkaKxoiMiUgMVT0ZhoiMiUgNRAgQHE5XomomOrUsiIlI0VnRERGrA1iURESmbDInORV8x7tKJTvrPf7QqyeTkSEgNym646ENE5FLKSqv/nEkuWj3diVw60d24cQMAcKhyh5MjITXY183ZEZCa3LhxA3q9Xr4bsnXpmoxGI3Jzc+Hr6wtBcPSRf/UoKSlBcHAwcnNzodPpnB0OKRj/rNlPkiTcuHEDRqNR3huLEhxuPbrorEuXTnQajQYtW7Z0dhguS6fT8ZcPNQr+WbOPrJUcuXaiIyIiG0li9eboPVwQEx0RkRqoeIyOD4yrkFarxRtvvAGtVuvsUEjh+GeN7gSCxDmsRESKVVJSAr1ej4gWL8FN49hfOKrECnx/dQ2Ki4tdasyVrUsiIjVg65KIiEiZWNEREamBBBkqOlkiaXRMdEREasDWJRERkTKxoiMiUgNRBODgA9+iaz4wzopO5XJzczF27Fhnh0EKkZmZiQ0bNuDs2bMAgLNnz2LixIkYO3Ys9u7d6+ToVK6mdeno5oKY6FSusLAQmzZtcnYYpABJSUno0aMHZs6ciZ49eyIpKQmPPPIIsrOzcfnyZQwaNIjJjpyCrUuF+/rrr295/OLFi40UCSldXFwcZs2ahcWLF2Pbtm149tlnMXHiRCxZsgQAMG/ePLz55psYMGCAkyNVKRVPRuHKKAqn0WggCMItX+IoCALMZnMjRkVKpNfrkZ6ejrCwMIiiCK1Wix9//BE9e/YEAJw+fRoRERHIz893cqTqYlkZxX8M3DQeDt2rSqzE94UbXG5lFLYuFS4oKAhffvklRFGsczt+/LizQyQFqXkvpEajgaenp9XrZnx9fVFcXOys0EjFmOgUrlevXkhPT6/3+J9Ve0S2at26Nc6fP2/5nJKSgpCQEMvnK1euICgoyBmhEQBJEmXZXBHH6BRu1qxZKCsrq/d4WFgY9u3b14gRkVJNnDjRqgXepUsXq+Pffvstx+ecSZIcf0O4i/6lmGN0REQKVjNGN9DvBbgJDo7RSZXYU7TZ5cboWNEREamBJMHhxSpdtC5ioiMiUgNRBAQHx9hcdIyOk1GIiKhBxMfH47777oOvry8CAgIwfPhwZGVlWZ1TXl6OSZMmoXnz5vDx8cHIkSNRUFAgaxxMdEREauCEJcAOHDiASZMm4ejRo0hOTobJZMKgQYOsJshNnz4dO3fuxGeffYYDBw4gLy8PI0aMkPVHZ6IjxRg9ejSGDx9u+dyvXz9Mmzat0ePYv38/BEFAUVFRvecIgoDExESb77lgwQL06NHDobguXboEQRCQkZHh0H3INUmiKMsGVE9w+eNWUVFR53cmJSVh9OjR6Ny5M7p3746NGzfiypUrlkeeiouLsX79eixfvhwDBgxAr169sGHDBhw5cgRHjx6V7WdnoqMGNXr0aAiCAEEQ4OHhgbCwMMTFxaGqqqrBv/vLL7/EokWLbDrXluRERNWCg4Oh1+stW3x8vE3X1SwY4O/vDwBIT0+HyWRCRESE5ZwOHTogJCQEKSkpssXLySjU4AYPHowNGzagoqIC33zzDSZNmgR3d3fMmzev1rmVlZXw8HBsCnSNmv+ZiAiyzrrMzc21erxAq9X+6aWiKGLatGl46KGHLM9Y5ufnw8PDA35+flbnBgYGyrpUHCs6anBarRYGgwGtWrXCxIkTERERYVlsuqbduGTJEhiNRrRv3x5A9f9ITz/9NPz8/ODv74+oqChcunTJck+z2YwZM2bAz88PzZs3x+zZs2ut8PLfrcuKigrMmTMHwcHB0Gq1CAsLw/r163Hp0iX0798fANCsWTMIgoDRo0cDqP6fMz4+HqGhofDy8kL37t3x+eefW33PN998g3bt2sHLywv9+/e3itNWc+bMQbt27dC0aVO0adMGsbGxMJlMtc5bu3YtgoOD0bRpUzz99NO1ltT66KOP0LFjR3h6eqJDhw744IMP7I6FFEqU5NkA6HQ6q82WRDdp0iScPn0a27Zta+iftBYmOmp0Xl5eqKystHzes2cPsrKykJycjF27dsFkMiEyMhK+vr44dOgQfvjhB/j4+GDw4MGW6959911s3LgRH3/8MQ4fPozCwkLs2LHjlt/7wgsv4JNPPkFCQgIyMzOxdu1a+Pj4IDg4GF988QUAICsrC9euXcP7778PoHrW2ObNm7FmzRqcOXMG06dPx3PPPYcDBw4AqE7II0aMwLBhw5CRkYEXX3wRc+fOtfvfia+vLzZu3IiffvoJ77//PtatW4f33nvP6pzs7Gxs374dO3fuRFJSEk6cOIGXX37ZcnzLli2YP38+lixZgszMTCxduhSxsbF8DRM53eTJk7Fr1y7s27cPLVu2tOw3GAyorKysNWRQUFAAg8EgXwASUQOKiYmRoqKiJEmSJFEUpeTkZEmr1UozZ860HA8MDJQqKios1/zjH/+Q2rdvL4miaNlXUVEheXl5Sbt375YkSZKCgoKkZcuWWY6bTCapZcuWlu+SJEl69NFHpalTp0qSJElZWVkSACk5ObnOOPft2ycBkH799VfLvvLycqlp06bSkSNHrM4dN26c9Mwzz0iSJEnz5s2TOnXqZHV8zpw5te713wBIO3bsqPf422+/LfXq1cvy+Y033pCaNGki/fzzz5Z93377raTRaKRr165JkiRJ99xzj7R161ar+yxatEgKDw+XJEmScnJyJADSiRMn6v1eUp7i4mIJgDTA4ylpkPZZh7YBHk9JAKTi4mKbvlsURWnSpEmS0WiUzp07V+t4UVGR5O7uLn3++eeWfWfPnpUASCkpKbL9O+AYHTW4Xbt2wcfHByaTCaIo4tlnn8WCBQssx7t27Wo1Lnfy5ElkZ2fD19fX6j7l5eW4cOECiouLce3aNfTp08dyzM3NDb179653geqMjAw0adIEjz76qM1xZ2dn4+bNm3jssces9ldWVlpePZOZmWkVBwCEh4fb/B01Pv30UyQkJODChQsoLS1FVVVVrSWWQkJC0KJFC6vvEUURWVlZ8PX1xYULFzBu3DiMHz/eck5VVZXVGwRIvSRRgiQ4NkZX3/9f9Zk0aRK2bt2Kr776Cr6+vpZxN71eDy8vL+j1eowbNw4zZsyAv78/dDodpkyZgvDwcDzwwAMOxfpHTHTU4Pr374/Vq1fDw8MDRqMRbm7Wf+y8vb2tPpeWlqJXr17YsmVLrXvdfffdtxWDl5eX3deUlpYCAP75z39aJRjAtsF3W6WkpCA6OhoLFy5EZGQk9Ho9tm3bhnfffdfuWNetW1cr8TZp0kS2WInssXr1agDV4+V/tGHDBss4+HvvvQeNRoORI0eioqICkZGRso8tM9FRg/P29kZYWJjN599777349NNPERAQUO/CsUFBQUhNTcUjjzwCoLpySU9Px7333lvn+V27doUoijhw4IDVVOYaNRXlH1ff79SpE7RaLa5cuVJvJdixY8dab3G39/mfI0eOoFWrVnjttdcs+y5fvlzrvCtXriAvLw9Go9HyPRqNBu3bt0dgYCCMRiMuXryI6Ohou76fVEISATTuEmC2VICenp5YtWoVVq1adbtR/SlORqE7TnR0NO666y5ERUXh0KFDyMnJwf79+/HKK6/g559/BgBMnToVb775JhITE3H27Fm8/PLLt3wGrnXr1oiJicHYsWORmJhouef27dsBAK1atYIgCNi1axd++eUXlJaWwtfXFzNnzsT06dOxadMmXLhwAcePH8fKlSstEzxeeuklnD9/HrNmzUJWVha2bt2KjRs32vXztm3bFleuXMG2bdtw4cIFJCQk1DmxxtPTEzExMTh58iQOHTqEV155BU8//bRl0H7hwoWIj49HQkICzp07h1OnTmHDhg1Yvny5XfGQMkmiJMvmipjo6I7TtGlTHDx4ECEhIRgxYgQ6duyIcePGoby83FLhvfrqq3j++ecRExOD8PBw+Pr64q9//est77t69Wo8+eSTePnll9GhQweMHz/eshRRixYtsHDhQsydOxeBgYGYPHkyAGDRokWIjY1FfHw8OnbsiMGDB+Of//wnQkNDAVSPm33xxRdITExE9+7dsWbNGixdutSun/eJJ57A9OnTMXnyZPTo0QNHjhxBbGxsrfPCwsIwYsQIDB06FIMGDUK3bt2sWjwvvvgiPvroI2zYsAFdu3bFo48+io0bN1piJVIrvo+OiEjBat5H1xdD4QZ3h+5VBRMO4xu+j46IiO4cHh4eMBgMOJz/jSz3MxgMsq1e1FhY0RERKVx5ebnVIg2O8PDwgKenpyz3aixMdEREpGicjEJERIrGREdERIrGREdERIrGREdERIrGREdERIrGREdERIrGREdERIrGREdERIr2/wFe5CDOSHALLAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay\n",
    "cm = confusion_matrix(test_y, pred)\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "ax = fig.add_subplot(111)\n",
    "cm_display = ConfusionMatrixDisplay(cm,\n",
    "                                    display_labels=stims).plot(ax=ax,\n",
    "                                                                       xticks_rotation=\"vertical\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Device: cuda:0\n"
     ]
    }
   ],
   "source": [
    "import torch \n",
    "\n",
    "torch.manual_seed(0)\n",
    "if torch.cuda.is_available():\n",
    "    torch.cuda.manual_seed_all(0)\n",
    "\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "print(\"Device:\", device)\n",
    "\n",
    "class LinearClassifierNet(torch.nn.Module):\n",
    "    def __init__(self, input_dim, output_dim):\n",
    "        torch.nn.Module.__init__(self)\n",
    "        self.linear1 = torch.nn.Linear(input_dim, 2048)\n",
    "        self.dropout1 = torch.nn.Dropout(0.5)\n",
    "        self.output = torch.nn.Linear(2048, output_dim)\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = self.linear1(x)\n",
    "        x = self.dropout1(x)\n",
    "        x = self.output(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsig, tlabel = torch.Tensor(train_signatures).to(device), torch.Tensor(train_y).to(device)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Epoch 0:   0%|          | 0/692 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([17816]) tensor(1, device='cuda:0')\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                \r"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "CUDA error: CUBLAS_STATUS_ALLOC_FAILED when calling `cublasCreate(handle)`",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[16], line 20\u001b[0m\n\u001b[1;32m     18\u001b[0m output \u001b[39m=\u001b[39m model(inputs) \n\u001b[1;32m     19\u001b[0m loss \u001b[39m=\u001b[39m loss_fn(output, label) \n\u001b[0;32m---> 20\u001b[0m loss\u001b[39m.\u001b[39;49mbackward() \n\u001b[1;32m     21\u001b[0m optimizer\u001b[39m.\u001b[39mstep() \n\u001b[1;32m     22\u001b[0m c_loss \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m loss\u001b[39m.\u001b[39mitem() \n",
      "File \u001b[0;32m~/miniconda3/envs/signatures/lib/python3.11/site-packages/torch/_tensor.py:488\u001b[0m, in \u001b[0;36mTensor.backward\u001b[0;34m(self, gradient, retain_graph, create_graph, inputs)\u001b[0m\n\u001b[1;32m    478\u001b[0m \u001b[39mif\u001b[39;00m has_torch_function_unary(\u001b[39mself\u001b[39m):\n\u001b[1;32m    479\u001b[0m     \u001b[39mreturn\u001b[39;00m handle_torch_function(\n\u001b[1;32m    480\u001b[0m         Tensor\u001b[39m.\u001b[39mbackward,\n\u001b[1;32m    481\u001b[0m         (\u001b[39mself\u001b[39m,),\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    486\u001b[0m         inputs\u001b[39m=\u001b[39minputs,\n\u001b[1;32m    487\u001b[0m     )\n\u001b[0;32m--> 488\u001b[0m torch\u001b[39m.\u001b[39;49mautograd\u001b[39m.\u001b[39;49mbackward(\n\u001b[1;32m    489\u001b[0m     \u001b[39mself\u001b[39;49m, gradient, retain_graph, create_graph, inputs\u001b[39m=\u001b[39;49minputs\n\u001b[1;32m    490\u001b[0m )\n",
      "File \u001b[0;32m~/miniconda3/envs/signatures/lib/python3.11/site-packages/torch/autograd/__init__.py:197\u001b[0m, in \u001b[0;36mbackward\u001b[0;34m(tensors, grad_tensors, retain_graph, create_graph, grad_variables, inputs)\u001b[0m\n\u001b[1;32m    192\u001b[0m     retain_graph \u001b[39m=\u001b[39m create_graph\n\u001b[1;32m    194\u001b[0m \u001b[39m# The reason we repeat same the comment below is that\u001b[39;00m\n\u001b[1;32m    195\u001b[0m \u001b[39m# some Python versions print out the first line of a multi-line function\u001b[39;00m\n\u001b[1;32m    196\u001b[0m \u001b[39m# calls in the traceback and some print out the last line\u001b[39;00m\n\u001b[0;32m--> 197\u001b[0m Variable\u001b[39m.\u001b[39;49m_execution_engine\u001b[39m.\u001b[39;49mrun_backward(  \u001b[39m# Calls into the C++ engine to run the backward pass\u001b[39;49;00m\n\u001b[1;32m    198\u001b[0m     tensors, grad_tensors_, retain_graph, create_graph, inputs,\n\u001b[1;32m    199\u001b[0m     allow_unreachable\u001b[39m=\u001b[39;49m\u001b[39mTrue\u001b[39;49;00m, accumulate_grad\u001b[39m=\u001b[39;49m\u001b[39mTrue\u001b[39;49;00m)\n",
      "\u001b[0;31mRuntimeError\u001b[0m: CUDA error: CUBLAS_STATUS_ALLOC_FAILED when calling `cublasCreate(handle)`"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "../aten/src/ATen/native/cuda/Loss.cu:199: nll_loss_forward_reduce_cuda_kernel_1d: block: [0,0,0], thread: [0,0,0] Assertion `t >= 0 && t < n_classes` failed.\n"
     ]
    }
   ],
   "source": [
    "BATCH_SIZE = 8\n",
    "LR = 0.0001\n",
    "EPOCHS = 20\n",
    "\n",
    "model = LinearClassifierNet(input_dim=train_signatures.shape[1], output_dim=1).to(device) # go from signature item to corresponding label \n",
    "loss_fn = torch.nn.CrossEntropyLoss()\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=LR)\n",
    "\n",
    "for e in range(EPOCHS): \n",
    "    # train eval\n",
    "    model.train()\n",
    "    c_loss = 0.0 \n",
    "    for idx in tqdm(range(tsig.shape[0]), desc=f\"Epoch {e}\", leave=False): \n",
    "    #for i, data in tqdm(enumerate(tsig), desc=f\"Epoch {e}\", leave=False): \n",
    "        inputs, label = tsig[idx], tlabel[idx].to(dtype=torch.long)\n",
    "        print(inputs.shape, label)\n",
    "        optimizer.zero_grad() \n",
    "        output = model(inputs) \n",
    "        loss = loss_fn(output, label) \n",
    "        loss.backward() \n",
    "        optimizer.step() \n",
    "        c_loss += loss.item() \n",
    "    \n",
    "    print(f\"train loss: {c_loss}\")\n",
    "\n",
    "    # test eval \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([17816])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tsig[0].shape"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "signatures",
   "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.11.0"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}