Version 4 RUL multiple_models- (RandomForest Depth Analysis).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Version 03 -> Pred RUL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: texttable in /home/yash/anaconda3/lib/python3.8/site-packages (1.6.3)\n",
      "3.8.5\n"
     ]
    }
   ],
   "source": [
    "!pip install texttable\n",
    "from platform import python_version\n",
    "print(python_version())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{132: ['B0018'], 168: ['B0005', 'B0006', 'B0007']}\n"
     ]
    }
   ],
   "source": [
    "# importing required libraries\n",
    "from scipy.io import loadmat\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from pprint import pprint as pp\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "from pprint import pprint\n",
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso, BayesianRidge, ARDRegression, SGDRegressor\n",
    "from texttable import Texttable\n",
    "import math\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "# getting the battery data\n",
    "# bs_all = [\n",
    "#     'B0005', 'B0006', 'B0007', 'B0018', 'B0025', 'B0026', 'B0027', 'B0028', 'B0029', 'B0030', 'B0031', 'B0032', \n",
    "#     'B0042', 'B0043', 'B0044', 'B0045', 'B0046', \n",
    "#     'B0047', 'B0048'\n",
    "# ]\n",
    "bs_all = [\n",
    "    'B0005', 'B0006', 'B0007', 'B0018'\n",
    "]\n",
    "ds = {}\n",
    "for b in bs_all:\n",
    "    ds[b] = loadmat(f'DATA/{b}.mat')\n",
    "    \n",
    "types = {}\n",
    "times = {}\n",
    "ambient_temperatures = {}\n",
    "datas = {}\n",
    "\n",
    "for b in bs_all:\n",
    "    x = ds[b][b][\"cycle\"][0][0][0]\n",
    "    ambient_temperatures[b] = x['ambient_temperature']\n",
    "    types[b] = x['type']\n",
    "    times[b] = x['time']\n",
    "    datas[b] = x['data']\n",
    "\n",
    "# clubbing all the compatible batteries together\n",
    "# Batteries are compatible if they were recorded under similar conditions\n",
    "# And their data size match up\n",
    "bs_compt = {}\n",
    "\n",
    "for b in bs_all:\n",
    "    sz = 0\n",
    "    for j in range(datas[b].size):\n",
    "        if types[b][j] == 'discharge':\n",
    "            sz += 1\n",
    "    if bs_compt.get(sz):\n",
    "        bs_compt[sz].append(b)\n",
    "    else: \n",
    "        bs_compt[sz] = [ b ]\n",
    "pp(bs_compt)\n",
    "BSSS = bs_compt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CRITICAL TIME POINTS FOR A CYCLE\n",
    "## We will only these critical points for furthur training\n",
    "\n",
    "## TEMPERATURE_MEASURED\n",
    "## => Time at highest temperature\n",
    "\n",
    "## VOLTAGE_MEASURED\n",
    "## => Time at lowest Voltage\n",
    "\n",
    "## VOLTAGE_LOAD\n",
    "## => First time it drops below 1 volt after 1500 time\n",
    "\n",
    "def getTemperatureMeasuredCritical(tm, time):\n",
    "    high = 0\n",
    "    critical = 0\n",
    "    for i in range(len(tm)):\n",
    "        if (tm[i] > high):\n",
    "            high = tm[i]\n",
    "            critical = time[i]\n",
    "    return critical\n",
    "\n",
    "def getVoltageMeasuredCritical(vm, time):\n",
    "    low = 1e9\n",
    "    critical = 0\n",
    "    for i in range(len(vm)):\n",
    "        if (vm[i] < low):\n",
    "            low = vm[i]\n",
    "            critical = time[i]\n",
    "    return critical\n",
    "\n",
    "def getVoltageLoadCritical(vl, time):\n",
    "    for i in range(len(vl)):\n",
    "        if (time[i] > 1500 and vl[i] < 1):\n",
    "            return time[i]\n",
    "    return -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MODEL\n",
    "* Considering 1 Cycle for RUL estimation\n",
    "\n",
    "## Features\n",
    "* [CP1, CP2, CP3, Capacity] -> RUL\n",
    "\n",
    "## Remaining Useful Life\n",
    "* n = number of cycles above threshold\n",
    "* RUL of Battery after (cycle x) = (1 - (x / n)) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## X: Features\n",
    "## y: RUL\n",
    "## x: no. of cycles to merge\n",
    "def merge(X, y, x):\n",
    "    XX = []\n",
    "    yy = []\n",
    "    \n",
    "    sz = len(X)\n",
    "    for i in range(sz - x + 1):\n",
    "        curr = []\n",
    "        for j in range(x):\n",
    "            for a in X[i + j]:\n",
    "                curr.append(a)\n",
    "        XX.append(curr)\n",
    "#         val = 0\n",
    "#         for j in range(x):\n",
    "#             val += y[i + j]\n",
    "#         val /= x\n",
    "        yy.append(y[i + x - 1])\n",
    "        \n",
    "    return XX, yy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Data Structure\n",
    "    # Cycles[battery][param][cycle]\n",
    "    # Cycles[battery][Capacity][cycle]\n",
    "\n",
    "Cycles = {}\n",
    "params = ['Temperature_measured', 'Voltage_measured', 'Voltage_load', 'Time']\n",
    "\n",
    "rmses = []\n",
    "\n",
    "for bs_cmpt in bs_compt:\n",
    "    rmses.append([])\n",
    "\n",
    "# iterate over the merge hyper parameter\n",
    "for xx in range(1, 10):\n",
    "\n",
    "    results = Texttable()\n",
    "    results.add_row(['Compatible Batteries', 'Cycles', 'MAE', 'RMSE', 'R2 Score' ])\n",
    "    \n",
    "    loc = 0\n",
    "    \n",
    "    # iterate over all the battery sets\n",
    "    for bs_cmpt in bs_compt:\n",
    "        # getting data for a given set\n",
    "\n",
    "        # y contains RUL after current cycle\n",
    "        # model will train for y\n",
    "        y = []\n",
    "\n",
    "        bs = bs_compt[bs_cmpt]\n",
    "        for b in bs:\n",
    "            Cycles[b] = {}\n",
    "            for param in params:\n",
    "                Cycles[b][param] = []\n",
    "                for j in range(datas[b].size):\n",
    "                    if types[b][j] == 'discharge':\n",
    "                        Cycles[b][param].append(datas[b][j][param][0][0][0])\n",
    "\n",
    "            cap = []\n",
    "            for j in range(datas[b].size):\n",
    "                if types[b][j] == 'discharge':\n",
    "                    cap.append(datas[b][j]['Capacity'][0][0][0][0])\n",
    "            Cycles[b]['Capacity'] = np.array(cap)\n",
    "            Cycles[b]['count'] = len(Cycles[b][params[0]])\n",
    "\n",
    "            effective_cycle_count = 0\n",
    "            for x in Cycles[b]['Capacity']:\n",
    "                if (x < 1.4):\n",
    "                    break\n",
    "                effective_cycle_count += 1\n",
    "\n",
    "            for i in range(len(Cycles[b]['Capacity'])):\n",
    "                if (i < effective_cycle_count):\n",
    "                    y.append((1 - ((i + 1) / effective_cycle_count)) * 100)\n",
    "                else:\n",
    "                    y.append(0)\n",
    "\n",
    "        # preparing data for regression model\n",
    "        temperature_measured = []\n",
    "        voltage_measured = []\n",
    "        voltage_load = []\n",
    "        capacity = []\n",
    "\n",
    "        for b in bs:        \n",
    "            for c in Cycles[b]['Capacity']:\n",
    "                capacity.append(c)\n",
    "            for i in range(Cycles[b]['count']):\n",
    "                temperature_measured.append(getTemperatureMeasuredCritical(Cycles[b]['Temperature_measured'][i], Cycles[b]['Time'][i]))\n",
    "                voltage_measured.append(getVoltageMeasuredCritical(Cycles[b]['Voltage_measured'][i], Cycles[b]['Time'][i]))\n",
    "                voltage_load.append(getVoltageLoadCritical(Cycles[b]['Voltage_load'][i], Cycles[b]['Time'][i]))\n",
    "\n",
    "        # creating the model\n",
    "        X = []\n",
    "        for i in range(len(temperature_measured)):\n",
    "            X.append(np.array([temperature_measured[i], voltage_measured[i], voltage_load[i], capacity[i]]))\n",
    "    #     X.append(np.array(capacity))\n",
    "        X = np.array(X)\n",
    "        y = np.array(y)\n",
    "        \n",
    "        # merge cycles\n",
    "        X, y = merge(X, y, xx)\n",
    "\n",
    "        # creating train test split\n",
    "        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state = 0)\n",
    "\n",
    "        # fitting the model\n",
    "        regressor = LinearRegression()\n",
    "        regressor.fit(X_train, y_train)\n",
    "\n",
    "        # test\n",
    "        y_pred = regressor.predict(X_test)\n",
    "\n",
    "        # model evaluation\n",
    "        diff = 0\n",
    "        total = 0\n",
    "        rmse = 0\n",
    "        for i in range(len(y_test)):\n",
    "            diff += abs(y_test[i] - y_pred[i])\n",
    "            rmse += ((y_test[i] - y_pred[i]) * (y_test[i] - y_pred[i]))\n",
    "            total += y_test[i]\n",
    "        diff /= len(y_test)\n",
    "        total /= len(y_test)\n",
    "        rmse = math.sqrt(rmse / len(y_test))\n",
    "    #     accuracy = ((total - diff) / total) * 100\n",
    "        accuracy = r2_score(y_test, y_pred)\n",
    "\n",
    "        # Adding evaluation to result array to print in a table\n",
    "        results.add_row([ str(bs), str(Cycles[bs[0]]['count']), diff, rmse, accuracy ])\n",
    "        \n",
    "        rmses[loc].append(rmse)\n",
    "        loc += 1\n",
    "\n",
    "    # printing results\n",
    "#     print(f'Evaluation: Clubbing Compatible Batteries for cycle param: {xx}\\n{results.draw()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(rmses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs0klEQVR4nO3dd3xW9fn/8deVBSSsEMLeKyyZAQIoDkBBrdYNiCIUkLoQax3fb63tr9/autBaBwIiWJUhIq2tCmJRRAgQ9gyEkYSdMBKSkH39/kikNEQIkHOfe1zPxyMPkvuc+77fIsk753PO+XxEVTHGGBO4gtwOYIwxxl1WBMYYE+CsCIwxJsBZERhjTICzIjDGmAAX4naAi1W3bl1t0aKF2zGMMcanrF27Nl1Vo8vb5nNF0KJFCxISEtyOYYwxPkVEkn9qmw0NGWNMgLMiMMaYAGdFYIwxAc6KwBhjApwVgTHGBDgrAmOMCXBWBMYYE+CsCIzxY0XFyicJqZzMyXc7ivFiVgTG+ClV5bm/b+HX8zfx+8+3uR3HeDErAmP81FtLk/h4VQqt6kawcMMBth7McDuS8VJWBMb4oXkJqbyyeCe3d2/MZw/1p1a1UP785Q63YxkvZUVgjJ9ZmniUZxds5qq2dfnzHV2oFR7KI9e24ftd6Szfle52POOFrAiM8SMbU0/y0IfraN+gBu+M7ElYSMm3+H19m9O4djX+9OV2iottnXLz36wIjPET+9KzGTNzDXVrhPH+6F5Ur/KfyYWrhATz5A3t2Howk883HXQxpfFGVgTG+IH0rDxGvb+aYlVmje5NvRpVz9nn1q6N6diwJi8vSiSvsMiFlMZbWREY4+Oy8woZM3MNRzJzmfFAL1pFVy93v6Ag4Zmh7dl/4jQfxqd4OKXxZlYExviwgqJiHvpoHVsOZPDWiB50bxZ53v0HtIvmyjZ1efPfu8jMLfBQSuPtrAiM8VGqyrMLNvPdzjReuO0KBnaoX6HnPTO0PSdyCpjy7W6HExpfYUVgjI96dfFO5q/dz+OD2jKsd7MKP69z41rc2q0RM37Yy+GMXAcTGl9hRWCMD/pbfDJvLk1ieO+mTBzY9qKf/+T1MRQXw2tf73QgnfE1jhaBiEwSka0iskVEZotI1TLba4nI5yKysXS/0U7mMcYfLNp6mOf/voWB7evxh1s7IyIX/RpN64QzMq45n6xNZdeRUw6kNL7EsSIQkcbAY0CsqnYGgoFhZXZ7GNimql2Ba4BXRSTMqUzG+LqEfcd5bPZ6ujSpzV9HdCck+NK/hR+5rg0RYSG8+JVNPRHonB4aCgGqiUgIEA6UvZNFgRpS8itNdeA4UOhwJmN8UtLRU/xiVgKNaldjxgO9CA8LufCTzqNORBgTrmnNku1HWb33eCWlNL7IsSJQ1QPAK0AKcAjIUNXFZXZ7E+hASUFsBiaqanHZ1xKR8SKSICIJaWlpTkU2xmsdycxl1Iw1hAYH8cGY3tSJqJwD5zH9W9KgZlX+9OV2VG3qiUDl5NBQJHAr0BJoBESIyMgyu90AbCjd3g14U0Rqln0tVZ2qqrGqGhsdHe1UZGO8UmZuAaNmrOZkTj4zR/eiaZ3wSnvtamHBTBrclvUpJ1m09XClva7xLU4ODQ0C9qpqmqoWAAuAfmX2GQ0s0BJJwF6gvYOZjPEpeYVFPPjBWpKOZjHlvp50blyr0t/jjh5NaFuvOi99lUhB0TkH5CYAOFkEKUCciISXngMYCGwvZ5+BACJSH4gB9jiYyRifUVysPPnJJlbuOcZLd3bhqrbOHA2HBAfx9JD27EnPZu6aVEfew3g3J88RrALmA+soGf8PAqaKyAQRmVC62x+AfiKyGfgGeFpVbcJ0Y4A/fbmdzzce5Okh7bm9RxNH32tgh3r0ahHJ60t2kZ1n12sEGvG1E0SxsbGakJDgdgxjHDX9+z3837+280C/Fjz/s46XdK/AxVqbfII73lnBpEHtmDjo4m9SM95NRNaqamx52+zOYmO8zD82HuT//rWdG69owHM3e6YEAHo2j2RIpwZMXbab9Kw8j7yn8Q5WBMZ4kRW703ly3kZ6t6jD5Lu7ERzkmRL40a+HxJBbWMwb3+zy6Psad1kRGOMlth/K5MEP1tI8Kpxp98dSNTTY4xlaR1dnWK+mfLwqhb3p2R5/f+MOKwJjvMCBk6d54P3VRFQJYdaY3tQKD3Uty8RBbQkLCeKVRYmuZTCeZUVgjMtO5uQzasZqcvKLmDmmF41qV3M1T70aVRl7VSv+tfkQG1JPuprFeIYVgTEuyi0oYuysBFKO5TD1vljaNzjnxnpXjB/QirrVw/jTFzb1RCCwIjDGJUXFysQ561mbcoLJ93Slb+sotyOdUb1KCI8NbMuqvcdZmnjU7TjGYVYExrhAVfndP7ayaOsRnrupIzd3aeR2pHMM792MFlHhvPhlIkXFdlTgz6wIjHHB29/u5m/xyTw4oBVjrmzpdpxyhQYH8esb2pN45BSfrtvvdhzjICsCYzxs/tr9vLwokZ93a8TTQ7x7jsUbr2hA16a1ee3rneQWFLkdxzjEisAYD/o28SjPfLqJ/m2ieOnOrgR5+IaxiyUiPDu0PYcycnn/h31uxzEOsSIwxkM278/goY/W0a5+DaaM7ElYiG98+8W1iuK69vV4+9skTmTnux3HOMA3/iUaUyonv5A/fbGdWSv2seVABoU+Mn9+8rFsRs9cTWR4GDNH96JGVfduGLsUTw9pT3ZeIW8tTXI7inHA5S16aoyHTVu2l3eX/WfJioiwYLo1q03P5nWIbR5J92a1ve6H7LGsPEbNWE1hsTL3F72pV7Oq25EuWkyDGtzRowkfrExmVL8WlbpKmnGfFYHxGcez85n2/R5u6FSf3/6sEwn7jrM2+QQJ+07w5r93UawQJBDToCY9m9cmtnkdejaPpElkNY/N4FlWTn4hY2au4VBGLh+Pi6N1dHVXclSGSYPb8Y+NB5n89U5eu6eb23FMJbIiMD7jnW+TyMkv5MnrY2hcuxqNuzXm1m6NATiVW8CG1JOsTT7B2uQTLFx/kA/jUwCoX7MKsc3r0KN5JLHNI+nYqCahwc6PihYWFfPwR+vYfCCDd++LpWfzSMff00mNalfjgf4tmLpsD2OvakmnRpW/bKZxhxWB8QmHMk4za2Uyt3VvQtv6Nc7ZXqNqKFe1jT6znGNRsbLjcOaZI4a1ySf41+ZDAFQLDaZr01pnjhh6NIus9EneVJX/+WwzSxPT+ONtnRncsX6lvr5bHrq6DXNWp/LiV4l8MKa323FMJXG0CERkEjAWUEqWqxytqrll9rkGeB0IBdJV9WonMxnf9MY3u1BVHq/gylnBQUKnRrXo1KgW9/dtAZSUydnF8M53u8/cMduufvUz5xliW0TSrE74ZQ0nvbZkF/MS9vPYdW24t0/zS34db1MrPJRHrm3DH7/Yzg9J6fRvU9ftSKYSOLZUpYg0BpYDHVX1tIjMA75Q1Zln7VMbWAEMUdUUEamnqued2MSWqgw8e9KyGPzaMu6La87vbulUaa+bk19YMpy07wQJySdYl3KCU7kl6/XWrV7lP+cZWkTSuVGtCl/u+dGqZP73sy3cHduEF+/o4tr5CafkFhQx8NXviIwI5R8PX+n190KYEudbqtLpoaEQoJqIFADhwMEy20cAC1Q1BeBCJWAC0+Svd1IlJIiHr21Tqa8bHhZCv9Z16de65Lfa4mJl19EsEpKPnymHRVuPABAWEkTXJrXOHDX0bB5JZETYOa/59bYjPLdwC9fGRPPH267wuxIAqBoazK+ub8cT8zby+aaDZ87TGN/l6OL1IjIR+CNwGlisqveW2f46JUNCnYAawF9U9YNyXmc8MB6gWbNmPZOTkx3LbLzLlgMZ3PzX5TxybRuevCHG4+9/NDP3zAnohOQTbD2YQUFRyfdM6+gIejaPPHPUcDInn3unryKmfg1mj48jPMx/T8EVFys3/XU5WXkFLHniaqqEeH41NXNxzndE4OTQUCTwKXAPcBL4BJivqh+etc+bQCwwEKgGrARuUtWdP/W6NjQUWEbNWM2G1JMse+paalVz//6A3IIiNqaeJKG0HNYmnyDjdMGZ7S2iwvn0l/2Iql7FxZSe8d3ONEbNWM1vb+7otRPnmf9wa2hoELBXVdNKQywA+gEfnrXPfkpOEGcD2SKyDOgK/GQRmMCxas8xvtuZxjND23tFCUDJsEifVlH0aVWydkBxsbInPYuEfSfYeyybkX2aB0QJAAxoW5f+baL46793cWdsE2p62Y18puKcvJg6BYgTkXApGSgdCGwvs8/fgatEJEREwoE+5exjApCq8tKiROrVqMKo0qt+vFFQkNCmXg2G9W7Gs0M7BNQdtyLCM0M6cCKngHe/2+12nEqTVxh4s6w6VgSqugqYD6yj5NLRIGCqiEwQkQml+2wHvgI2AauB6aq6xalMxnf8e8dR1iaf4LGBbakWZuPP3uqKJrW4pWsj3lu+l8MZuRd+ghcrKlbeW76Xbr//msfnrKfAR+axqgyOnix2gp0j8H/FxcqNb3zP6YIiljxxtUfuAjaXLvV4Dte9+i139GjCn+/o4nacS7LzyCmemr+JDakn6dy4JlsOZHJd+3q8NaKH3/wicr5zBPYdZrzO55sOsuPwKZ4Y3M5KwAc0rRPOyLjmzEtIZdeRU27HuSj5hcX8Zckubnrje5KPZfOXYd34/JEreeG2K1iaeJRRM1aTmVtw4RfycfZdZrxKQVExk7/eSfsGNfiZF67ja8r36HVtiQgL4cWvEt2OUmEbU0/ys78u57UlOxnauSFLnriaW7s1RkQY0acZbwzrzrqUEwyfGk96Vp7bcR1lRWC8ytw1qSQfy+GpITF2x6oPqRMRxoRrWrNk+xHW7DvudpzzOp1fxB//tY3b3v6BjNMFTL8/ljeGdz/naq+fdW3EtFGx7E7L4u4pKzlw8rRLiZ1nRWC8xun8It74ZhexzSO5Nqae23HMRRrTvyX1a1bhhS+2463nHlfsTueG15cx7fu9DOvdjMVPDGDQeSYEvDamHn/7RR/SsvK4650V7E7L8mBaz7EiMF5j1sp9HD2Vx1ND2vvl1Az+rlpYMJMGtWN9ykkWbT3sdpz/kplbwLMLNjNi2ipEYPa4OF647YoK3fvQq0Ud5oyPI7+omLunrGTLgQwPJPYsKwLjFTJOF/DOt7u5Jiaa3i3ruB3HXKI7ezahdXQEL32V6DWXXy7ZdoTBk79j7poUHhzQiq8mDqBv66iLeo1OjWox78G+VA0NZvjUeFbtOeZQWndYERivMG3ZHjJOF/Dk9Z6fT8hUnpDgIJ4e0p496dnMXZPqapb0rDwenb2esR8kEBkexsKH+/PsjR0u+XLQVtHV+WRCX+rVrML9M1azdIf/zJFpRWBcl3Yqjxk/7OXmLg3p3NhWvfJ1gzvWJ7Z5JK8v2UV2XqHH319VWbj+AIMnf8dXWw7xxOB2/OORK+nSpPZlv3aj2tWY92Bf2tWvwbgPEvj7hgOXH9gLWBEY1721NIm8wmJ+ZUcDfkFEePbG9qRn5TH9+70efe+DJ08zZuYaHp+7gRZ1I/jisat4bGDbCq8lURFR1avw8bg+9GweyeNzN/C3eN+fDdmKwLgq9XgOH61K5u7YJrSsG+F2HFNJejavww2d6jN12W6PXINfXKz8LT6Z619bRvye4/z25o7Mn9Cv3GVNK0ONqqHMGtOb62Lq8dzCLby1NMlrr5SqCCsC46rXl+xCRHhsYMWWoDS+46kh7cktLOaNb3Y5+j570rIYNi2e5xZuoVvT2iyeNIAxV7Yk2OH7UKqGBjPlvp78vFsjXl6UyJ++3OGzZeC/K2cYr7fzyCkWrN/P2Ctb0rBWNbfjmErWOro69/RqyserUhjdv2WlH/EVFhUzffleXitdwe6lO7twV88mHr30ODQ4iMl3d6NmtVCmLttDRk4BL9x+heMlVNnsiMC45pVFiUSEhfDLayp3CUrjPR4f2JbQ4CBeWVy5U09sO5jJz9/+gT9/uYNrYqJZ8sTV3B3b1JX7T4KChN/f0onHrmvD3IRUHp29zuemsrYiMK5Yn3KCxduOMO6qVtQpZ+1f4x/q1azKuKta8q9Nh9iYevKyXy+3oIhXFiVyy5vLOZyRxzv39uDd+2KpV7Pq5Ye9DCLCE9fH8JubOvDF5sOMnZVATr7nr5i6VFYExhUvL0okKiKMX1xlSxz6u/FXtyYqIow/fXl5U0+sTT7OTW98z5tLk7i1W2OWPDGAoVc0rMSkl2/sVa146Y4u/JCUzsjpq8jI8Y2ZS60IjMct35XOit3HeOjaNlSvYqep/F31KiE8NrAt8XuO821i2kU/PzuvkN/9Yyt3TllJbkExs8b05tW7u1I73DuPJO/u1ZS37+3BlgOZ3DN1JUdPef+CPVYExqNUlZcX7aBRrarc26eZ23GMhwzv3YzmUeH8+csdFBVX/Khg2c40rn9tGbNW7mNU3xYsnjSAq9tFO5i0cgzp3JAZD/Qi5XgOd01ZSerxHLcjnZejRSAik0Rkq4hsEZHZIlLuQJ6I9BKRIhG508k8xn2Lth5m4/4MHh/cjqqh/rHyk7mwsJAgfn1DDIlHTrFg3f4L7n8yJ58nP9nI/TNWUzU0iE8e7MvvbulEhA8dQV7Zti4fje3DyZwC7pyygp1evGiPY0UgIo2Bx4BYVe0MBAPDytkvGHgRWORUFuMdioqVVxbvpHV0BLd3b+x2HONhN13RkK5NajH5653kFvz0VTVfbj7EoMnLWLj+AI9c24Z/PXYVsS18cyLC7s0imfdgX1Th7ndXsqESTpg7wemhoRCgmoiEAOHAwXL2eRT4FPCfGZxMuRas20/S0SyevD6GEFuCMuCICM8M7cChjFxmrth3zvajmblM+NtafvnROhrUqsLfH+nPkzfE+PyRY0yDGsyf0I8aVUO4d1o8K5LS3Y50Dse+G1X1APAKkAIcAjJUdfHZ+5QeNdwGTDnfa4nIeBFJEJGEtLSLP9lk3JdXWMTrS3bRpUkthnRu4HYc45K+raO4Niaat5cmcTInHyg5bzQvIZVBk79jaeJRnhnanoUP9adTI/+ZgLBZVDjzJ/SjSWQ4D8xcw2IvW6/ByaGhSOBWoCXQCIgQkZFldnsdeFpVz3v3hapOVdVYVY2Njvb+E0XmXB+vSuHAydP8+oYYW3QmwD01pD2n8gp5a2kSqcdzuO+91Tw1fxPtG9bky4lXMeHq1n55xFi/ZlXmPhhHx4Y1+eVH6/h07YXPlXiKk2deBgF7VTUNQEQWAP2AD8/aJxaYU/qDoS5wo4gUqupCB3MZD8vOK+TNfyfRt1UUV7ap63Yc47IODWtye/cmzFqRzIfxKQQHCf/3886M6N3M79eprh0exkdj+zD+bwn86pONZJwuYMyV7t9L42TtpgBxIhIuJT/pBwLbz95BVVuqagtVbQHMBx6yEvA/M5bv5Vh2Pr8eYkcDpsSvrm9H1dAg4lrVYfGkAYyMa+73JfCjiCohzHigFzd0qs//++c2Xvt6p+uT1Tl2RKCqq0RkPrAOKATWA1NFZELp9vOeFzD+4UR2PlOX7WFwx/r0aBbpdhzjJRrVrsba5wYT6odDQBVRJSSYt0b04JkFm/nLN7vIOF3Ab2/u6FoZOnpRrqo+Dzxf5uFyC0BVH3Ayi3HHlO92k5VfaEtQmnMEagn8KCQ4iJfu6EKtaqG8t3wvmbkFvHRHF1fOj/jO3RnG5xwuvUzwtm6NiWngzAIhxviyoCDhNzd1oHa1UF79eieZpwt5c0R3j18yG9iVbBz1l292UazKpMHt3I5ijNcSER4d2Jbf39KJJduPMPr9NWR5eK1nKwLjiL3p2cxLSGV472Y0rRPudhxjvN6ofi14/Z5urN53nBHT4jmene+x97YiMI6Y/PVOwoKDeOQ6W3TGmIr6effGvDuyJ4mHT3HPuys5nOGZmUutCEyl23owg883HmR0/xbUq+HugiHG+JpBHesza0xvDmXkcsc7K9iXnu34e1oRmEr3yqJEalYN4cEBrd2OYoxPimsVxcfj+pCTX8idU1ay7WCmo+9nRWAq1Zp9x1mamMYvr2lDrfBQt+MY47O6NKnNJxP6EhosDJu6krXJxx17LysCU2lUlZe+2kG9GlV4oF8Lt+MY4/Pa1KvBJxP6ElW9CiOnr+a7nc5MumlFYCrNt4lprNl3gkcHtqVamG9PHWyMt2gSGc68B/vSom4ESUezHHkPu6HMVIriYuWlRYk0qxPOPbFN3Y5jjF+JrlGFzx7q59iNZnZEYCrFPzcfYvuhTJ4Y3I6wEPtnZUxlc/JuY/uONZetoKiYyYsTad+gBrd0beR2HGPMRbIiMJftk4T97DuWw5PXxwTMVMLG+BMrAnNZcguK+Ms3O+nRrDYDO9RzO44x5hJYEZjL8sHKfRzJzOOpIe1t0RljfJQVgblkmbkFvP3tbga0iyauVZTbcYwxl+iCRSAiQSLSzxNhjG+ZvmwPJ3MKeOoGW3TGGF92wSJQ1WLg1Ut5cRGZJCJbRWSLiMwWkapltt8rIptKP1aISNdLeR/jeWmn8pi+fC83XdGQzo1ruR3HGHMZKjo0tFhE7pCLGAQWkcbAY0CsqnYGgoFhZXbbC1ytql2APwBTK/r6xl1vLU0ir7CYJ663RWeM8XUVvbP4CSACKBKR04AAqqo1K/D61USkAAgHDp69UVVXnPVlPNCkgnmMi/afyOHjVSnc2aMJraOrux3HGHOZKlQEqnrRC86q6gEReQVIAU4Di1V18Xme8gvgy/I2iMh4YDxAs2bNLjaKqWSvL9kFAhMHtXU7ijGmElT4qiERuUVEXin9uLkC+0cCtwItgUZAhIiM/Il9r6WkCJ4ub7uqTlXVWFWNjY6Ormhk44BdR06xYN1+7o9rTqPa1dyOY4ypBBUqAhH5MzAR2Fb6MbH0sfMZBOxV1TRVLQAWAOdcfSQiXYDpwK2qeuxiwhvPe3XxTsLDQnjoWluC0hh/UdFzBDcC3UqvIEJEZgHrgWfO85wUIE5EwikZGhoIJJy9g4g0o6Qg7lPVnReZ3XjYxtSTfLX1MI8PakudiDC34xhjKsnFTENdG/hxiZwLXi+oqqtEZD6wDiikpDimisiE0u1TgN8CUcDbpRckFapq7EVkMh708qJE6kSEMfaqVm5HMcZUoooWwQvAehFZSskVQwOAZy/0JFV9Hni+zMNTzto+FhhbwQzGRSuS0lmelM5vbupA9Sq2jIUx/uSC39EiEgQUA3FAL0qK4GlVPexwNuMlVJUXFyXSqFZVRsY1dzuOMaaSXbAIVLVYRB5R1XnAPzyQyXiZxduOsDH1JC/ecYWji2MYY9xR0ctHvxaRJ0WkqYjU+fHD0WTGKxQVK68sSqRVdAR39LD7/YzxRxUd7B1T+ufDZz2mgJ019HML1x9g19Es3hrRg5Bgm6zWGH9U0XMEz6jqXA/kMV4kv7CY15bspHPjmgzt3MDtOMYYh1R09tGHL7Sf8T+zV6ew/8Rpfn1De1uC0hg/ZucITLmy8wr567+T6NOyDgPa1nU7jjHGQXaOwJTrhS+2cyw7j6n397QlKI3xcxWdfbSl00GM9/hm+xE+WpXCgwNa0aNZpNtxjDEOO+/QkIg8ddbnd5XZ9oJToYx70k7l8dT8TXRoWNMWnTEmQFzoHMHZK4qVnVJiSCVnMS5TVZ75dBOn8gr5y7BuVAmxm8eMCQQXKgL5ic/L+9r4uI9Xp/DNjqM8O7Q97epf9FpExhgfdaEi0J/4vLyvjQ/bnZbFH/65java1mVU3xZuxzHGeNCFThZ3FZFMSn77r1b6OaVfV3U0mfGYgqJiJs3dQNXQYF65q6vdM2BMgDlvEaiqDRIHgDe+2cWm/RlMGdmD+jWt340JNDZ5TIBL2Hect5YmcVfPJgzp3NDtOMYYF1gRBLBTuQVMmreBJpHhPH9LJ7fjGGNc4mgRiMgkEdkqIltEZLaIVC2zXUTkDRFJEpFNItLDyTzmv/3+820cOHGa1+7paquOGRPAHCsCEWkMPAbEqmpnIJj/vi8BYCjQtvRjPPCOU3nMf/ti8yHmr93PI9e2oWdzmzbKmEDm9NBQCCVXG4UA4cDBMttvBT7QEvFAbRGxgWqHHc7I5dkFm+napBaPDmzrdhxjjMscKwJVPQC8AqQAh4AMVV1cZrfGQOpZX+8vfey/iMh4EUkQkYS0tDSnIgeE4mLlyU82lqw1cE83Qm2xGWMCnpNDQ5GU/MbfEmgERIjIyLK7lfPUc25UU9WpqhqrqrHR0dGVHzaAvL9iH8uT0nnu5o60iq7udhxjjBdw8tfBQcBeVU1T1QJgAdCvzD77gaZnfd2Ec4ePTCXZcTiTF7/awaAO9Rjeu+mFn2CMCQhOFkEKECci4VIyof1AYHuZff4B3F969VAcJcNHhxzMFLByC4p4fM4GalYN4c93dLE1BowxZzh2zaCqrhKR+cA6oBBYD0wVkQml26cAXwA3AklADjDaqTyB7tXFiew4fIoZD8RSt3oVt+MYY7yIoxePq+rzwPNlHp5y1nbF1kN23A9J6Uz7fi8j45pxXfv6bscxxngZu2TEz2XkFPCreRtpFR3B/97Y0e04xhgvZLeT+jFV5X8WbiY9K48F9/ejWpjNIWiMOZcdEfixhRsO8K9Nh5g0uB1dmtR2O44xxktZEfip1OM5/HbhVmKbRzLh6tZuxzHGeDErAj9UVKz8at5GFHjtnm4E20IzxpjzsHMEfujdZbtZve84r97VlaZ1wt2OY4zxcnZE4Ge2HMhg8uKd3HRFQ27vcc60TcYYcw4rAj9yOr+IiXPWU7d6Ff54W2e7e9gYUyE2NORH/vTldnanZfPR2D7UDg9zO44xxkfYEYGfWLrjKB+sTGbslS3p36au23GMMT7EisAPHMvK49fzN9G+QQ2evCHG7TjGGB9jQ0M+TlV5ZsFmMk8X8OHY3lQNtbuHjTEXx44IfNzcNal8ve0ITw2JoX2Dmm7HMcb4ICsCH7Y3PZvff76N/m2iGNO/pdtxjDE+yorARxUUFfP43A2EhQTxyl1dCbK7h40xl8jOEfioN/+dxMbUk7w1ogcNa1VzO44xxofZEYEPWpt8gjeXJnF7j8bc1KWh23GMMT7OsSIQkRgR2XDWR6aIPF5mn1oi8rmIbBSRrSJiS1VeQFZeIU/M20DDWlX5/S2d3I5jjPEDTq5ZnAh0AxCRYOAA8FmZ3R4Gtqnqz0QkGkgUkY9UNd+pXL7uD59vI/V4DnPG96VG1VC34xhj/ICnhoYGArtVNbnM4wrUkJJJcaoDxylZ6N6U46sth5mbkMovr2lN75Z13I5jjPETniqCYcDsch5/E+gAHAQ2AxNVtbjsTiIyXkQSRCQhLS3N2aRe6mhmLs8u2MQVjWsxcWA7t+MYY/yI40UgImHALcAn5Wy+AdgANKJkGOlNETnnrihVnaqqsaoaGx0d7WBa76SqPDl/E6cLinjtnm6Ehdg5fmNM5fHET5ShwDpVPVLOttHAAi2RBOwF2nsgk0/5YGUyy3am8b83daRNvepuxzHG+BlPFMFwyh8WAkih5PwBIlIfiAH2eCCTz9h15BQvfLGda2OiGdmnmdtxjDF+yNEbykQkHBgMPHjWYxMAVHUK8AdgpohsBgR4WlXTnczkS/ILi5k4ZwPVq4Tw0p1dbaEZY4wjHC0CVc0Boso8NuWszw8C1zuZwZdN/non2w5lMu3+WKJrVHE7jjHGT9lZRy8Vv+cY7y7bzfDezRjcsb7bcYwxfsyKwAtlnC7gV/M20iIqgudu7uB2HGOMn7NJ57zQ83/fwuHMXD79ZT/Cw+x/kTHGWXZE4GX+vuEACzccZOLAtnRrWtvtOMaYAGBF4EUOnDzNbxZuoUez2jx0TWu34xhjAoQVgZcoLlZ+NW8DxcXKa/d0IyTY/tcYYzzDftp4ienL9xC/5zjP39KJ5lERbscxxgQQKwIvsO1gJi8vSmRIpwbc1bOJ23GMMQHGisBluQVFPD53PZHhYbxw+xV297AxxuPs2kQXnc4vYtwHCew6msXM0b2pExHmdiRjTACyInBJTn4hY2clsHLPMV6+sytXtwu86bWNMd7BisAF2XmF/GLWGlbvPc7ku7tyW3c7L2CMcY8VgYdl5RUy5v01JCQf57V7unFrt8ZuRzLGBDgrAg86lVvA6PfXsD71JG8M787NXRq5HckYY6wIPCUzt4AHZqxm0/4M/jq8Ozde0dDtSMYYA1gReERmbgH3v7eaLQcyeHNED4Z0buB2JGOMOcOKwGEZOQXcP2MV2w5l8va9Pbi+k5WAMca7OHZDmYjEiMiGsz4yReTxcva7pnT7VhH5zqk8bjiZk8+978Wz/dAppozsaSVgjPFKjh0RqGoi0A1ARIKBA8BnZ+8jIrWBt4EhqpoiIvWcyuNpJ7LzuXf6KpLSsnj3vp5c295v/tOMMX7GU0NDA4Hdqppc5vERwAJVTQFQ1aMeyuOo46UlsDsti2n3x9rNYsYYr+apuYaGAbPLebwdECki34rIWhG5v7wni8h4EUkQkYS0tDRHg16uY1l5jJgWz560LN4bZSVgjPF+jheBiIQBtwCflLM5BOgJ3ATcADwnIu3K7qSqU1U1VlVjo6O99wdr2qk8hk+LZ9+xbGY80Iur2npvVmOM+ZEnhoaGAutU9Ug52/YD6aqaDWSLyDKgK7DTA7kq1dFTuYyYtooDJ07z/gO96ds6yu1IxhhTIZ4YGhpO+cNCAH8HrhKREBEJB/oA2z2QqVIdycxl2NR4Dp48zczRvawEjDE+xdEjgtIf7oOBB896bAKAqk5R1e0i8hWwCSgGpqvqFiczVbbDGbkMnxbP0cxcZo3pTa8WddyOZIwxF8XRIlDVHCCqzGNTynz9MvCykzmccijjNMOnxpOelc8Hv+hNz+ZWAsYY32N3Fl+iAydLSuBEdkkJ9GgW6XYkY4y5JFYElyD1eA7Dp8WTcbqAv43tQ7emtd2OZIwxl8yK4CKlHs9h2NR4TuUW8NHYPnRpUtvtSMYYc1msCC5C8rFshk+NJzu/iI/HxdG5cS23IxljzGWzIqigfenZDJ8WT25BER+P60OnRlYCxhj/YEVQAXvSshg+LZ6CIuXjcXF0aFjT7UjGGFNprAguIOloSQkUFyuzx8UR06CG25GMMaZSWRGcx64jpxg+bRUAc8bH0ba+lYAxxv94avZRn5N4+BTDp8UTJFYCxhj/ZkVQju2HMhk+LZ7gIGHO+Dja1KvudiRjjHGMFUEZWw9mMGJaPGHBQcwZ35dW0VYCxhj/ZucIzrLlQAYj31tFeGgws8fH0Twqwu1IxhjjODsiKLV5f8mRQERYCHMf7GslYIwJGFYEwMbUk4yYHk/NaqHMGR9H0zrhbkcyxhiPCfihoXUpJxj13moiI8KYPT6OxrWruR3JGGM8KqCLYG3ycUbNWENU9TBmj4ujkZWAMSYABezQ0Jp9x7n/vdVE16jC3PF9rQSMMQHLsSIQkRgR2XDWR6aIPP4T+/YSkSIRudOpPGdbtecYo2aspn6tqswZH0eDWlU98bbGGOOVHBsaUtVEoBuAiAQDB4DPyu5Xuu1FYJFTWc62cvcxxsxcQ+PIanw8rg/1algJGGMCm6eGhgYCu1U1uZxtjwKfAkedDvFDUjqjZ66maZ1qzB4XZyVgjDF4rgiGAbPLPigijYHbgCnnPOO/9xsvIgkikpCWlnZJAVYkpTNm5hpaREXw8bg4omtUuaTXMcYYf+N4EYhIGHAL8Ek5m18HnlbVovO9hqpOVdVYVY2Njo6+pBz1alalT6soPh4XR93qVgLGGPMjT1w+OhRYp6pHytkWC8wREYC6wI0iUqiqCys7RJt61flgTO/KflljjPF5niiC4ZQzLASgqi1//FxEZgL/dKIEjDHG/DRHh4ZEJBwYDCw467EJIjLByfc1xhhTcY4eEahqDhBV5rFyTwyr6gNOZjHGGFO+gL2z2BhjTAkrAmOMCXBWBMYYE+CsCIwxJsBZERhjTIATVXU7w0URkTSgvDmLKqIukF6JcSqLt+YC781muS6O5bo4/piruaqWOzWDzxXB5RCRBFWNdTtHWd6aC7w3m+W6OJbr4gRaLhsaMsaYAGdFYIwxAS7QimCq2wF+grfmAu/NZrkujuW6OAGVK6DOERhjjDlXoB0RGGOMKcOKwBhjAlxAFIGIzBCRoyKyxe0sZxORpiKyVES2i8hWEZnodiYAEakqIqtFZGNprt+7nelsIhIsIutF5J9uZ/mRiOwTkc0iskFEEtzO8yMRqS0i80VkR+m/s75ekCmm9O/px49MEXnc7VwAIjKp9N/8FhGZLSJesbC5iEwszbTVib+rgDhHICIDgCzgA1Xt7HaeH4lIQ6Chqq4TkRrAWuDnqrrN5VwCRKhqloiEAsuBiaoa72auH4nIE5SsbldTVW92Ow+UFAEQq6pedROSiMwCvlfV6aXLxoar6kmXY50hIsHAAaCPql7qjaKVlaUxJf/WO6rqaRGZB3yhqjNdztUZmAP0BvKBr4BfququynqPgDgiUNVlwHG3c5SlqodUdV3p56eA7UBjd1OBlsgq/TK09MMrfmMQkSbATcB0t7N4OxGpCQwA3gNQ1XxvKoFSA4HdbpfAWUKAaiISAoQDB13OA9ABiFfVHFUtBL4DbqvMNwiIIvAFItIC6A6scjkKcGb4ZQNwFPhaVb0iF/A68BRQ7HKOshRYLCJrRWS822FKtQLSgPdLh9Kmi0iE26HKGMZPLGXraap6AHgFSAEOARmqutjdVABsAQaISFTpqo83Ak0r8w2sCLyAiFQHPgUeV9VMt/MAqGqRqnYDmgC9Sw9PXSUiNwNHVXWt21nK0V9VewBDgYdLhyPdFgL0AN5R1e5ANvCMu5H+o3So6hbgE7ezAIhIJHAr0BJoBESIyEh3U4GqbgdeBL6mZFhoI1BYme9hReCy0jH4T4GPVHXBhfb3tNKhhG+BIe4mAaA/cEvpePwc4DoR+dDdSCVU9WDpn0eBzygZz3XbfmD/WUdz8ykpBm8xFFinqkfcDlJqELBXVdNUtYCStdb7uZwJAFV9T1V7qOoASoa5K+38AFgRuKr0pOx7wHZVnex2nh+JSLSI1C79vBol3yA7XA0FqOqzqtpEVVtQMqTwb1V1/Tc2EYkoPdlP6dDL9ZQczrtKVQ8DqSISU/rQQMDVCxHKGI6XDAuVSgHiRCS89HtzICXn7VwnIvVK/2wG3E4l/705uni9txCR2cA1QF0R2Q88r6rvuZsKKPkN9z5gc+l4PMD/qOoX7kUCoCEwq/SKjiBgnqp6zaWaXqg+8FnJzw5CgI9V9St3I53xKPBR6TDMHmC0y3kAKB3rHgw86HaWH6nqKhGZD6yjZOhlPd4z1cSnIhIFFAAPq+qJynzxgLh81BhjzE+zoSFjjAlwVgTGGBPgrAiMMSbAWREYY0yAsyIwxpgAZ0VgjDEBzorAGGMC3P8HR8ywKPikijUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwd0lEQVR4nO3dd3iV9fnH8fedTQYJkBAgg52wgxAiiqIoIohbUGmtv9payk9qHVVrW0erVesebdU6Ovy1TtyiKIiouCAghA0BAgkrYSVAyL5/f+TEYgiQkDznOeN+Xde5cpLnyXk+V64k9/mO5/sVVcUYY0zwCnE7gDHGGHdZITDGmCBnhcAYY4KcFQJjjAlyVgiMMSbIWSEwxpgg52ghEJEEEZkhIqtFZJWInNTo+A9FJM/z+FJEspzMY4wx5nBhDr/+48AsVZ0kIhFAdKPjG4HTVHWPiEwAngFOdDiTMcaYQ4hTN5SJSHtgKdBLm3EREekALFfVFEcCGWOMaZKTLYJeQAnwD0+XzyLgOlU9cITzfwp8cKwXTUxM1B49erRZSGOMCQaLFi3aqapJTR1zskWQDXwNjFLVb0TkcaBMVW9v4twxwJPAKaq6q4njU4GpAOnp6cM3bdrkSGZjjAlUIrJIVbObOubkYHERUKSq33g+nwEMayLcEOA54IKmigCAqj6jqtmqmp2U1GRBM8YYc5wcKwSquh0oFJFMz5fOBFYeeo6IpANvAD9S1bVOZTHGGHNkTs8auhb4j2fG0AbgKhGZBqCqTwN3AJ2AJ0UEoOZITRdjjDHOcLQQqOoSoPE/9qcPOX41cLWTGYwxxhyd3VlsjDFBzgqBMcYEOSsExhgT5IKmEGwrPcgf3l1BdW2d21GMMcanBE0hyCsq5R9fFPDkJ+vdjmKMMT4laArB2QO7cOHQbvx57jpWbC11O44xxviMoCkEAL8/fyAdYiK46bU8qmqsi8gYYyDICkFCdAT3XTSYVdvK+Msn+W7HMcYYnxBUhQBg7IBkLh6Wwl8/yWf5FusiMsaYoCsEAHeeO5DE2Ah+9epSKmtq3Y5jjDGuCspCEB8dzp8uHsKaHft44uN1bscxxhhXBWUhABjTrzOXZqfy1Lz1LC3c63YcY4xxTdAWAoDbzh1AcvsofvXaUiqqrYvIGBOcgroQtI8K5/5LhpBfvJ9H59h2CMaY4BTUhQBgdEYSU3LSefazDSzatMftOMYY43VBXwgAfjexP13j23GzdREZY4KQFQIgNjKMByYNYcPOAzz04Rq34xhjjFdZIfAY1SeRH43szvNfbGRhwW634xhjjNdYITjErRP6kdqhvouovKrG7TjGGOMVVggOERMZxoOTsijYVc4Ds6yLyBgTHKwQNDKyVyd+fHIP/vllAV9v2OV2HGOMcZwVgibcMj6THp2iuXnGUg5UWheRMSawWSFoQnREGA9OzqJoz0H+9MFqt+MYY4yjHC0EIpIgIjNEZLWIrBKRkxodFxF5QkTyRSRPRIY5maclRvToyE9H9eT/vt7EF/k73Y5jjDGOcbpF8DgwS1X7AVnAqkbHJwB9PY+pwFMO52mRm87OpFdiDLfMyGNfRbXbcYwxxhGOFQIRaQ+MBp4HUNUqVd3b6LQLgBe03tdAgoh0dSpTS0WFh/Lg5Cy2lR7k3veti8gYE5icbBH0AkqAf4jItyLynIjENDonBSg85PMiz9e+R0SmikiuiOSWlJQ4l7gJw7t34Gen9uKlBZv5bK13r22MMd7gZCEIA4YBT6nqCcAB4NZG50gT36eHfUH1GVXNVtXspKSktk96DDeclUHvpBhufT2PMusiMsYEGCcLQRFQpKrfeD6fQX1haHxO2iGfpwJbHcx0XKLCQ3n40qFsL6vgnvcaD3MYY4x/c6wQqOp2oFBEMj1fOhNY2ei0d4ArPbOHRgKlqrrNqUytMTQtgWmn9eaV3EI+WVPsdhxjjGkzTs8auhb4j4jkAUOBe0VkmohM8xx/H9gA5APPAtc4nKdVrhvbl4zkWG59PY/ScusiMsYEBlE9rEvep2VnZ2tubq5r119WVMqFT37BhUNTePjSLNdyGGNMS4jIIlXNbuqY3VncQoNT45l+em9eX1zEnJU73I5jjDGtZoXgOPzijL706xLHb95cxt7yKrfjGGNMq1ghOA4RYSE8fGkWew5U8ft3VrgdxxhjWsUKwXEa2C2ea8/oy1tLtjJr+Xa34xhjzHGzQtAK14zpzcBu7bntrWXsPmBdRMYY/2SFoBXCQ+u7iEoPVnPH28vdjmOMMcfFCkEr9evSnuvHZvBe3jZm5vnkvXDGGHNUVgjawM9H92JIajy3v72cnfsr3Y5jjDEtYoWgDYSFhvDQ5Cz2V9Rw+1vL8beb9Iwxwc0KQRvJSI7jhrMy+GD5dt61LiJjjB+xQtCGfnZqT4amJXDH28sp3lfhdhxjjGkWKwRtqKGLqLyqlt+9aV1Exhj/YIWgjfXpHMvN4zKZvXIHby3Z4nYcY4w5JisEDvjJKT0Z3r0Dd769gh1l1kVkjPFtVggcEBoiPDhpCFW1dfzmjWXWRWSM8WlWCBzSKymWW87ux9zVxby+2LqIjDG+ywqBg358cg9yenTkD++uYFvpQbfjGGNMk6wQOCgkRHhw8hBqapVbX7cuImOMb7JC4LDunWL4zTn9+HRtCa/mFrodxxhjDmOFwAuuOLE7J/XqxN3vrWLLXusiMsb4FisEXhASIjwwaQiqyq9n5FkXkTHGp1gh8JK0jtH8dmJ/5ufv5MUFm92OY4wx33G0EIhIgYgsE5ElIpLbxPF4EXlXRJaKyAoRucrJPG77QU46p/RJ5J6ZqyjcXe52HGOMAbzTIhijqkNVNbuJY9OBlaqaBZwOPCwiEV7I5AoR4f5JQwgR4ZYZedTVWReRMcZ9bncNKRAnIgLEAruBGncjOSsloR23TezPVxt28e9vNrkdxxhjHC8ECnwkIotEZGoTx/8C9Ae2AsuA61S1rvFJIjJVRHJFJLekpMTZxF5w2Yg0Rmckcd/7q9m064DbcYwxQc7pQjBKVYcBE4DpIjK60fGzgSVAN2Ao8BcRad/4RVT1GVXNVtXspKQkhyM7T0S4/5LBhIUKN79mXUTGGHc5WghUdavnYzHwJpDT6JSrgDe0Xj6wEejnZCZf0TW+HXecO4AFBbv555cFbscxxgQxxwqBiMSISFzDc2AcsLzRaZuBMz3nJAOZwAanMvmaScNTOaNfZx74cDWbd9ksImOMO5xsESQD80VkKbAAmKmqs0RkmohM85xzN3CyiCwDPgZ+rao7HczkU0SEey8aTIgI976/yu04xpggFebUC6vqBiCria8/fcjzrdS3FIJWl/go/ve03jw8ey1frd/FSb07uR3JGBNk3J4+aoCfje5Ft/go/jhzJbU2cGyM8TIrBD4gKjyUW8/pz4qtZby+qMjtOMaYIGOFwEecN6Qrw9ITeODDNeyvDOh76owxPsYKgY8QEe44byA791fy5Cf5bscxxgQRKwQ+ZGhaAhefkMJz8zfaonTGGK+xQuBjbh6fSagIf/pgtdtRjDFBwgqBj+ka345pp/Vm5rJtLNi42+04xpggYIXAB00d3Yuu8VHc9d4KW4fIGOM4KwQ+qF1EKLdO6MfyLWW8vtimkxpjnGWFwEedn9WNEzzTSQ/YdFJjjIOsEPgoEeH2cwdQsq+Sp+atdzuOMSaAWSHwYcPSO3Dh0G48+/kGivbYdFJjjDOsEPi4W8b3QwTun7XG7SjGmABlhcDHdUtox9TRvXl36VYWbbLppMaYtmeFwA9MO60Xye0juevdlTad1BjT5qwQ+IHoiDB+Pb4fS4tKeWvJFrfjGGMCjBUCP3Hh0BSyUuO5f9ZqyqtsOqkxpu1YIfATISHCHecNYEdZJU9/GjTbOhtjvMAKgR8Z3r0j52V142+frmfL3oNuxzHGBAgrBH7m1+MzAXhglq1OaoxpG1YI/Exqh2imju7F20u2smjTHrfjGGMCgBUCPzTttN50jovk7vdsOqkxpvUcLQQiUiAiy0RkiYjkHuGc0z3HV4jIp07mCRQxkWHcMr4fSwr38s7SrW7HMcb4OW+0CMao6lBVzW58QEQSgCeB81V1IDDZC3kCwsUnpDA4pX466cGqWrfjGGP8mNtdQz8A3lDVzQCqWuxyHr/RMJ10W2kFz3xm00mNMcfP6UKgwEciskhEpjZxPAPoICLzPOdc2dSLiMhUEckVkdySkhJHA/uTET06MnFIV57+dD3bSm06qTHm+DhdCEap6jBgAjBdREY3Oh4GDAcmAmcDt4tIRuMXUdVnVDVbVbOTkpIcjuxfbh3fj1pVHrTVSY0xx8nRQqCqWz0fi4E3gZxGpxQBs1T1gKruBD4DspzMFGjSOkZz9Sk9eePbLSwp3Ot2HGOMH3KsEIhIjIjENTwHxgHLG532NnCqiISJSDRwIrDKqUyB6poxfUiMjeSud1egatNJjTEt42SLIBmYLyJLgQXATFWdJSLTRGQagKquAmYBeZ5znlPVxsXCHENsZBi3nJ3J4s17eTdvm9txjDF+Ro71DlJEQoCRqvqldyIdXXZ2tubmNnlLQlCrrVPO/8t89hyoYu5NpxMVHup2JGOMDxGRRU1N44dmtAhUtQ54uM1TmTYVGlK/2f3W0gqetemkxpgWaG7X0EcicomIiKNpTKuM7NWJCYO68OS89ewoq3A7jjHGTzS3ENwIvAZUiUiZiOwTkTIHc5nj9JsJ/amtUx6w6aTGmGZqViFQ1ThVDVHVcFVt7/m8vdPhTMuld4rmJ6f05PXFRSwrKnU7jjF+pbKmlrvfW8nSIJuK3exZQyJyvog85Hmc62Qo0zrTx/QmMTaCu96z6aTGNJeqctuby3l+/kZ+/XpeUK3s26xCICJ/Aq4DVnoe13m+ZnxQXFQ4N43LZGHBHt5ftt3tOMb4hX98UcBri4o4sWdHVm/fF1Qr+za3RXAOcJaq/l1V/w6M93zN+KjJ2Wn079qee99fRUW1rU5qzNF8vq6EP85cybgBybz4s5EM6Nqeh2evoaqmzu1oXtGSG8oSDnke38Y5TBurn07any17D/L8/I1uxzHGZ23ceYBfvPgtfTvH8chlQwkNEW4en0nh7oO8snCz2/G8ormF4F7gWxH5p4j8C1jk+ZrxYSf3TuTsgck8+Uk+xTad1JjD7Kuo5mcv5CICz16ZTWxkGACnZySR06MjT8zNp7yqxuWUzjtmIfDcWVwHjATe8DxOUtWXHc5m2sBvz+lPVW0dD31k00mNOVRtnXLdy0vYuPMAT/5wGOmdor87JiLcMj6Tkn2V/PPLAvdCeklz7yz+hapuU9V3VPVtVbURSD/RvVMMV43qyWuLili+xaaTGtPgoY/WMHd1MXeeN4CTeycedjy7R0fO6NeZp+etp7S82oWE3tPcrqHZInKTiKSJSMeGh6PJTJv5xRl96BgdwV3vrbTppA6prVNeWbiZfRWB/Q8jULy9ZAtPzVvPlJx0fjSy+xHPu2lcJmUVNfzts/VeTOd9zS0EPwGmU79fwCLPw1Z+8xPto8K5cVwGCzbuZtZya8w54c1vt/Dr15fx7Oc2MO/r8or2csuMPHJ6dOQP5w/kaCvnDOjWnguGduPvX2wM6HG25o4R3KqqPRs9enkhn2kjl2Wn0a9LHPd+YNNJ21p1bR1/nrsOgFcXFlJTGxxTDv1R8b4Kpr6wiMTYSJ68YhgRYcd+L3zjWRnU1Cp/npvvhYTuaO4YwXQvZDEOCgsN4fZzB1C4+yD/+KLA7TgB5c3FW9i0q5zJw1PZXlbBvDW2r7YvqqypZdr/LaL0YDXPXDmcxNjIZn1f904xXDYijZcWbGbzrnKHU7rDxgiCyKg+iYztn8xfP8mnZF+l23ECQnVtHU/MXcfglHjuuWgwSXGRvBwkc8/9iaryuzeXs3jzXh6anMXAbi27FeqXZ/YlLFR4dM5ahxK6y8YIgszvJvansqaWR2bbdNK2MGNREUV7DnLDWX2JCAvh0uxU5q4uZlvpQbejmUP8/YsCZiwq4pdn9GHikK4t/v7k9lH8z8k9eGvJFlZvD7yFl5u7+mjj8QEbI/BTPRNj+J+TevDywkJWbLXppK1RVVPHX+bmk5WWwJjMzgBcPiKdOoVXFxa5nM40+GxtCfd4lo+4fmzGcb/O/57Wm9jIMB76MPBaBUctBCJyyyHPJzc6ZncW+6lrz+xLQrtw7rbppK3yam4hW/Ye5MazMr6beZLWMZpT+ybyysLN1AbR6pW+qn75iMX07RzHo5cNJSTk+PfWSoiO4OejezFn1Q4Wbdrdhindd6wWweWHPP9No2Pj2ziL8ZL4duHcOC6Trzfs5qOVO9yO45cqqmv56yf5DEtPYHTf79+MNCUnna2lFXy21gaN3VTmWT4iNER47n+yifEsH9EaV43qSWJsJA/MWhNQb6KOVQjkCM+b+tz4kSkj0shIjuXe91dRWWPTSVvqlYWFbCut4MazMg+bhz62fzKJsRG8tMAGjd1SW6dc71k+4q8/HEZax+hjf1MzxESGce0Zffhm424+W7ezTV7TFxyrEOgRnjf1ufEjDdNJN+0q519BsJZKW6qoruXJefmM6NGBUX06HXY8IiyEScPT+Hh1se0d7ZKG5SN+f4TlI1pjSk46qR3a8eCHqwNm85pjFYKshj2KgSGe5w2fDz7Wi4tIgYgsE5ElInLEWUYiMkJEakVkUgvzm1Y4tW8SZ/brzJ8/zmfnfptO2lwvfrOZHWWV3HDI2EBjl49Io7ZOeS230MvpzKHLR1xxlOUjjldEWAg3jM1g+ZYyPgiQO/WPWghUNfSQPYrDPM8bPg9v5jXGqOpQVc1u6qCIhAL3Ax+2MLtpA7+d2J+D1bU8MjvwZkI44WBVLU99up6RvToe9Z1mj8QYRvXpxEsLCgPmXaM/aMnyEa1x4QkpZCTH8vBHawLiTvKWbEzjlGuB14Fit4MEo95JsVx5Ug9eXrCZVdsCb350W/vPN5so2VfJDc2Yhnj5iHS27D3I5/mB05fsy4rLWr58xPEKDRF+NS6TDTsPMGOR/08VdroQKPCRiCwSkamND4pICnAR8LTDOcxRXHdmX9q3C+ePM2066dGUV9Xw1Lz1jOrTiRN7HT420Ni4gcl0jIngZRs0dlxlTS0//3fLl49ojXEDkhmalsDjH6/z+/W7nC4Eo1R1GDABmC4ioxsdfwz4taoe9acoIlNFJFdEcktKbEpeW4uPDueGsRl8kb+Lj1dZw+xI/u+rTew6UNWs1gBAZFgok4anMnvlDor32aCxUxqWj/j2OJePOF4Nm9dsK63g319v8so1neJoIVDVrZ6PxcCbQE6jU7KBl0WkAJgEPCkiFzbxOs+oaraqZiclJTkZOWj94MR0+nSO5Z73VwXNht0tcaCyhr99toFT+yaS3aP5y2xdPiKNmjoNiO4DX9Xa5SNa4+TeiZzaN5G/fpLv13tROFYIRCRGROIangPjgOWHnuNZqqKHqvYAZgDXqOpbTmUyRxYeGsJtE/uzcecBXviqwO04PudfXxWw+0AVN5zVsiUKeiXFMrJXR162QWNHtNXyEa1x89mZ7Cmv5jk/3ovCyRZBMjBfRJYCC4CZqjpLRKaJyDQHr2uO0+mZnTk9M4nHP17H7gNVbsfxGfsqqnnmsw2MyUxiWHqHFn//lJx0Nu8u58v1uxxIF7zacvmI1hiSmsCEQV147vMN7PLTadiOFQJV3aCqWZ7HQFW9x/P1p1X1sMFhVf2xqs5wKo9pntsm9qe8qpZHbTrpd/75RQF7y6uP+x3n2QO7kBAdzku2PHWbcWL5iNb41bgMDlbX8tdP/HNLS1+YPmp8SJ/OcfxoZHf+880m1mzf53Yc15VVVPPs5xsY278zWWkJx/UaUeGhXDIslY9WbLcb99qAU8tHtEafznFMGp7Kv7/exJa9/rcEuRUCc5jrzuxLXJRNJwX4+/yNlFXUtLr/eUpOGtW1yus2aNxqTi4f0RrXeX5HHvfDzWusEJjDdIiJ4Pqxffl83U5mLtvmdhzXlJZX8/znGxk3IJlBKa2bktincxw5PTry0oLNQV9cW8Pp5SNaIyWhHVeM7M6MRUXkF+93O06LWCEwTfrRyO4MTUvgd28uD9rdtp6fv4F9la1vDTSYcmIaBbvK+WqDDRofD28tH9Ea08f0pl14qN/tAGiFwDQpLDSExy4bSnVtHb96dWnQTX3cW17F378oYMKgLgzo1r5NXnPCoK60jwrjpQW2EF1LeXP5iNboFBvJT0/txfvLtpNXtNftOM3mmz9N4xN6JMbw+/MG8uX6XTw3f4Pbcbzq2c83cKCq7VoDUD9ofPGwVD5cvt2m57aAG8tHtMbPTu1Jh+hwHvzQf1oFVgjMUU3OTmX8wC48+OEalm8Jjj2Odx+o4h9fFDBxcFcyu8S16WtPyUmnqraONxbboHFzHLp8xMOXem/5iNaIiwrnmtP78Pm6nXy53j8WHLRCYI5KRLjv4sF0jIng+leWcLDKvxfXao6/fbaeg9W1XHdm3zZ/7cwucQzv3oEXbdC4Wb5bPuLMvpwz2LvLR7TGj07qTtf4KL/Z0tIKgTmmDjERPDQ5i/zi/dz3wSq34zhq5/5KXvhyE+dndaNvctu2BhpMyUlnQ8kBFmwMrA3Q21rD8hFnD0zmegeKspOiwkO57sy+LCncy2w/2BfcCoFpllP7JnH1KT154atNfLI6cFco/dun66msqeWXDv7jmTi4K3FRYban8VE0LB+RkRzHI5e6t3xEa0wankqvxBge+mgNtT4+2cIKgWm2m8dn0q9LHDfPWBqQd8gW76vg/77exIVDU+idFOvYddpFhHLRCSm8v3w7e8tt0Lixsopqrv7XQkJDhGevdH/5iOMVFhrCjeMyWLtjP28v2eJ2nKOyQmCaLTIslCemnEBZRQ23zMjzi77Plnh63gaqa5VrvdANcfmIdKpq6nhjsW//g/C2huUjNu0q58kfDveJ5SNa45xBXRnYrT2PzF7r08u7WyEwLZKRHMdvJvRj7upi/v1N4HRt7Cir4N/fbOKiE1LomRjj+PUGdGvP0LQEu9O4kQc/rF8+4s7zBnBS72PvAufrQkKEm8/OpGjPQZ/uCrRCYFrsxyf3YHRGEvfMXOl3t9IfyVPz1lNbp/zyDO8NSv4gJ511xftZtGmP167py95esoWnP13PD070veUjWuO0jCRO7NmRP8/Np7yqxu04TbJCYFpMRHho0hCiI8K4/pVvfbrJ2xzbSg/y4jebmTw8lfRO3uuKODerK7GRYbzow+8UveXQ5SN+f55vLh9xvOq3tOzHzv2V/OOLArfjNMkKgTkundtH8aeLB7N8SxmP+PneBX/9JJ86VaaP6ePV60ZHhHHB0G7MzNtGabn/bnPYWocuH/GUDy8f0RrDu3dgbP/OPP3pep+cIBB4P3HjNeMGdmFKTjp/+2w9X/np7ltb9h7klYWFXDoizZWBySk56VTW1PGWj88qcUpF9X+Xj3j2ymw6+fjyEa1x09mZ7K+s4elPfW+5FisEplVuP7c/PTvFcOOrS/zyXe1f5uYjiNdbAw0GpcQzJDU+KAeNVZXb3vrv8hFttbifr+rXpT0XZHXjn19uZEdZhdtxvscKgWmV6IgwHrt8KCX7KvndW8v86p9Z4e5yXsst5LIRaaQktHMtx5ScdFZv38e3hXtdy+CG5+dv9MvlI1rjhrMyqKlVnvh4ndtRvscKgWm1IakJ3HBWBu/lbePNb/2ni+Mvc/MJCRGuGdPb1RznZXUjOiKUlwJoOu6xLC3cy30frGbcAP9bPqI1uneKYUpOOq8sLGTTrgNux/mOFQLTJqad1pucHh254+0VFO4udzvOMW3adYAZi4v4QU46XePdaw0AxEbWDxq/l7eNsgr/615rqcqaWm6esZSk2EgenJzll8tHtMa1Z/QhLFR8apKFFQLTJkJDhEcuy0KAG15ZQk2tb08p/fPcfMJChGtOd7c10GBKTjoHq2t5e8lWt6M47s8f57N2x37uu3gw8e3C3Y7jdZ3bR3HVqJ68s3Qrq7aVuR0HsEJg2lBqh2j+eNEgcjft4cl5692Oc0Qbdx7gjcVFXDGyO53bR7kdB4DBKfEM7NaeF78J7EHjZUWlPPXpeiYNT2VMv85ux3HNtNG9iYsM4yEf2bzG0UIgIgUiskxElohIbhPHfygieZ7HlyKS5WQe47wLhqZwwdBuPP7xOr7d7Jt3zP7543VEhIUw7TTfaA1A/U1HU3LSWbWtjLyiwNwAqKqmjptnLKVTTAS3TxzgdhxXxUeH8/PTevPx6mJyC9xfjtwbLYIxqjpUVbObOLYROE1VhwB3A894IY9x2F0XDKJL+yhueGUJByp965b6/OL9vLVkC1ee1IOkON+as37B0G60Cw/16TVpWuMvn+Szevs+7r1oMPHRwdcl1NhVo3qQGBvpE5vXuNo1pKpfqmrD28avgVQ385i2Ed8unEcuzWLT7nLuenel23G+54mP1xEVHsrPR/dyO8ph4qLCOS+rK+8s3cp+HyugrbV8SylPfpLPxSekMHZAsttxfEJ0RBjXndmHBQW7mbe2xNUsThcCBT4SkUUiMvUY5/4U+KCpAyIyVURyRSS3pMTdH5hpnhN7deJ/T+vNK7mFzFq+ze04AKzbsY9387Zy5Uk9fPYO1ik56ZRX1fJOAA0a13cJ5dEhJoI7zgvuLqHGLhuRTlrHdjw4aw11Lm5e43QhGKWqw4AJwHQRGd3USSIyhvpC8OumjqvqM6qararZSUlJzqU1ber6sRkMSY3n1jeWsb3U/TspH/t4HdHhoUz1wdZAg6FpCfTrEhdQ3UNPzVvPqm1l3HPhIBKiI9yO41MiwkK48awMVm4rY+Yy994wOVoIVHWr52Mx8CaQ0/gcERkCPAdcoKr+uWCNaVJEWAiPXTaUyuo6bnptqavveFZvL2Nm3jauGtWTjjG++89IRPjBieks21LKsgAYNF65tYw/z13HBUO7MW5gF7fj+KTzs1LITI7jkdlrqXZp2rVjhUBEYkQkruE5MA5Y3uicdOAN4Eeq6jt3V5g20yspltvPHcD8/J38/YuNruV4fM464iLDuPrUnq5laK4LhqYQFR7CSwv9u1VQXVs/SyghOpzfnzfQ7Tg+KzREuOnsTDbuPMCMRUWuZHCyRZAMzBeRpcACYKaqzhKRaSIyzXPOHUAn4MkjTTE1/m9KThpnDUjmgVlrXLmBZsXWUj5Yvp2rTunpF10T8e3CmTi4G29/u8XnZl21xNPz1rNiaxl/vHAwHXy4FeYLxvbvzLD0BB6bs5aK6lqvX9+xQqCqG1Q1y/MYqKr3eL7+tKo+7Xl+tap28EwvPdIUU+PnRIQ/XVw/ZfD6l5d4/Rf98TnriIsK46en+H5roMEPTkzjQFUt7+X556Dxmu37eGLuOs4d0pXxg6xL6FhEhJvP7seOskpe+KrA69e3O4uNV3SKjeTBSUNYs2Mff/pgtdeuu6yolI9W7uDqU3r51XIGw9I7kJEcy4sLCt2O0mI1tfVjQu2jwvnD+dYl1Fwn9e7E6Iwknpy33utrTlkhMF5zemZnfnxyD/75ZQGfemne9GNz1hLfLpyrTunhleu1lYY7jZcW7mXFVv8aNP7bZxtYtqWUuy8c5LPTdH3VLWdnsre8muc+8+7mNVYIjFfdOqEfGcmx3PTaUnbtr3T0WksL9/Lx6mJ+dmpP2kf5T2ugwUUnpBAZFsLLftQqWLdjH4/PWcc5g7sEzR4DbWlQSjwTB3flufkb2enw38ehrBAYr4oKD+Xxy0+gtLyaW99wdiObR+esJSE6nB+P8p+xgUMlREdwzuCuvPXtFsqrfH/QuKa2jptm5BETGcpdFwxyO47funFcBpU1dfz1k3yvXdMKgfG6/l3bc8v4TGav3MHLC515t7t48x7mrSnh56N7ExsZ5sg1vGFKTjr7KmuYmecbd2cfzXPzN7K0cC93XTCIROsSOm69k2KZNCyV/3y9maI93tnbwwqBccVPRvXklD6J3PXuSjaU7G/z13909lo6xkRw5Und2/y1vWlEjw70Torx+TuN84v388jstYwf2IVzh1iXUGtdN7YvCDw2xztbWlohMK4ICREempxFZHgI17+ypE3vqMwt2M3n63Yy7bRexPhxawD+O2i8ePNeVm/3jU1MGqutU26esZToiFDuvnAQIsG145gTuiW048qR3XljcRHrduxz/HpWCIxrusRHcd9Fg8krKuWxOW13Y/mjc9aSGBvJj0b2aLPXdNMlw1KJCPXdQeO/z9/It5v38ofzB/rc0t7+7JoxfYiOCOPhj5xfdMEKgXHVhMFduTQ7lSfnrWfBxtZv0PHNhl18kb+Laaf1ol1EaBskdF+HmAjGD+rCG4uLOFjl/btOj2Z9yX4e+mgNZw1I5vysbm7HCSgdYyK4+tSezFqxnaWFex29lhUC47o7zxtIesdobnhlCaUHW3cjzaNz1pIUF8kVI/17bKCxKTnplFXU8L6LK1Q2Vlun3DIjj6jwUO6xLiFHXH1qLzrGRPCgw1taWiEwrouJDOOxy4ayvayCO95efuxvOIIv1+/k6w27mX56b6LCA6M10GBkr470TIzhZR9aiO6fXxawaNMe7jxvgM/s/RxoYiPDuOb03szP38kX+Tsdu44VAuMTTkjvwC/P6MvbS7by9pItLf5+VeXR2Wvp0j6Ky3PSHUjorvpB4zQWFuzxyuDhsWzceYAHP1zNmf06c9EJKW7HCWhXjOxOt/goHvjQuS0trRAYnzF9TG+Gd+/AbW8ub/H86S/yd7GwYA/TxwRea6DBJcNSCQ8VXnJ50LiuTvn1jDwiQkO49+LB1iXksKjwUK4fm8HSwr18uGKHI9ewQmB8RlhoCI9eOhQFbnxlKbXN3MhGVXlk9hq6xUdx6Yg0Z0O6qFNsJGcP7MLri4tcWaq4wQtfFbCgYDe3nzuAZOsS8oqLh6XQv2t7Cnc7c4OZFQLjU9I7RfOH8weyoGA3T3+6vlnf8+naEhZv3sv0M/oQGRaYrYEGU3LSKT1Yzazl2125/qZdB7h/1hpOz0xi0vBUVzIEo7DQEN79xSh+5tA2q1YIjM+5eFgKE4d05dHZa8kr2nvUc1WVR+esIyWhHZOHB25roMFJvTrRvVO0K3ca13lmCYWFCPdZl5DXhYU69+/aCoHxOSLCvRcOJikukutfXnLUBdc+WVPM0sK9XHtGHyLCAv/XOSREuHxEOt9s3M16B5bmOJp/f7OJbzbWdwl1jW/n1WsbZwX+X47xS/HR4Tx8aRYbdx3g7vdWNXlO/UyhdaR1bMclQdRNMWl4KmEhwstebBUU7i7nTx+sZnRGEpOzg+dnHSysEBifdXLvRKae2ouXFmzmoxWH94nPWVXMsi2l/PKMvoQ72Gz2NUlxkYwbmMyMRUVU1jg/aNzQJRQi1iUUqILnr8f4pRvHZTCwW3tufWMZxWUV33294b6BHp2ig3Ie++Uj0tlTXu3YdMJDvbhgM19t2MXvJvYnJcG6hAKRFQLj0yLDQnn88qEcqKzhphl51HmmlH64Ygcrt5XxyzP7OjqI5qtO6ZNIaod2jncPFe0p5773V3FKn0QuD+CpucEu+P6CjN/p0zmO2yb257O1JfzrqwLq6pTH5qylV1JM0C50FhJSvzz1l+t3sXHnAUeuoarc+voyAP50iXUJBTJHC4GIFIjIMhFZIiK5TRwXEXlCRPJFJE9EhjmZx/ivK0Z254x+nbnvg9U8MXcdq7fv47ogbQ00mDw8ldAQcWz9oZcXFjI/fye/Oac/qR2iHbmG8Q3e+Csao6pDVTW7iWMTgL6ex1TgKS/kMX5IRLj/kiHERYbx2Jx19Okcy7lDgrM10KBz+yjG9u/MjNwiqmrabmMfgC17D3LPzFWc3LsTPwjAtZvM97n9duoC4AWt9zWQICK2z51pUlJcJA9OHkJEWAg3jcskNMS6Ki7PSWfXgSpmr2y7QWNV5TdvLKNOlfsvGUKI/ZwDntOFQIGPRGSRiExt4ngKcOgKWkWer32PiEwVkVwRyS0pKXEoqvEHZ/RLZukd4xg/qIvbUXzC6L5JpCS0a9Puoddyi/hsbQm3TuhHWkfrEgoGTheCUao6jPouoOkiMrrR8abeahy20piqPqOq2aqanZSU5ERO40cCZeexthAaIlw2Io3P1+1k867WL0i2rfQgd7+3khN7duSKEwNrcx9zZI4WAlXd6vlYDLwJ5DQ6pQg4dE5aKrDVyUzGBJpLs9MIEVrdKmjoEqqpUx6YZF1CwcSxQiAiMSIS1/AcGAc03n7qHeBKz+yhkUCpqvrOXnzG+IEu8VGc0S+ZV3OLqK49/kHj1xdvYd6aEm4Zn0n3TjFtmND4OidbBMnAfBFZCiwAZqrqLBGZJiLTPOe8D2wA8oFngWsczGNMwJqSk8bO/ZV8vOr4Bo23l1bwh3dXkNOjI/9zUo+2DWd8XphTL6yqG4CsJr7+9CHPFZjuVAZjgsVpGUl0jY/ixQWFjB/Usol3qsrv3lxGdW0d91uXUFBye/qoMaYNhIWGcGl2Gp+vK2nxLlZvLdnCx6uLufnsfvRMtC6hYGSFwJgAcemINAR4Nbf5exoXl1Xw+3dWMrx7B358cg/HshnfZoXAmACRktCO0zM788rCQmqaMWisqvz2zeVUVNfywKQhdoNeELNCYEwAmZKTTvG+SuauLj7mue8s3cqcVTu4aVwmvZNivZDO+CorBMYEkDGZSSS3jzzmnsYl+yq5850VnJCewE9O6emldMZXWSEwJoA0DBp/uraELXsPNnmOqnLbW8sor6rlwUlZ1iVkrBAYE2guzU5DgVcXNj1o/F7eNj5csYMbz8qgT2frEjJWCIwJOGkdoxndN4lXcw8fNN65v5I73l5OVloCV1uXkPGwQmBMAJqSk8620go+Xfv91XrvfHsFBypreWjSkKDe1Md8n/0mGBOAzuzfmcTY7w8az8zbxsxl27hubF/6Jse5mM74GisExgSg8NAQLs1OZe7qYraXVrDL0yU0OCWen4/u5XY842OsEBgToC4fkU6d1t9p/Pt3V1JWUc1Dk7OsS8gcxrFF54wx7krvFM2pfRN5+tP1lFfV8quzMsjsYl1C5nD21sCYADYlJ53yqloGdmvPtNN7ux3H+ChrERgTwM4akMzVp/RkyonphFuXkDkCKwTGBLDw0BBuO3eA2zGMj7O3CMYYE+SsEBhjTJCzQmCMMUHOCoExxgQ5KwTGGBPkrBAYY0yQs0JgjDFBzgqBMcYEOVFVtzO0iIiUAJuO89sTgZ1tGKet+Gou8N1slqtlLFfLBGKu7qqa1NQBvysErSEiuaqa7XaOxnw1F/huNsvVMparZYItl3UNGWNMkLNCYIwxQS7YCsEzbgc4Al/NBb6bzXK1jOVqmaDKFVRjBMYYYw4XbC0CY4wxjQRFIRCRv4tIsYgsdzvLoUQkTUQ+EZFVIrJCRK5zOxOAiESJyAIRWerJ9Qe3Mx1KREJF5FsRec/tLA1EpEBElonIEhHJdTtPAxFJEJEZIrLa83t2kg9kyvT8nBoeZSJyvdu5AETkBs/v/HIReUlEotzOBCAi13kyrXDiZxUUXUMiMhrYD7ygqoPcztNARLoCXVV1sYjEAYuAC1V1pcu5BIhR1f0iEg7MB65T1a/dzNVARG4EsoH2qnqu23mgvhAA2arqU3PPReRfwOeq+pyIRADRqrrX5VjfEZFQYAtwoqoe7/1BbZUlhfrf9QGqelBEXgXeV9V/upxrEPAykANUAbOA/1XVdW11jaBoEajqZ8But3M0pqrbVHWx5/k+YBWQ4m4q0Hr7PZ+Gex4+8Y5BRFKBicBzbmfxdSLSHhgNPA+gqlW+VAQ8zgTWu10EDhEGtBORMCAa2OpyHoD+wNeqWq6qNcCnwEVteYGgKAT+QER6ACcA37gcBfiu+2UJUAzMVlWfyAU8BtwC1LmcozEFPhKRRSIy1e0wHr2AEuAfnq6050Qkxu1QjVwOvOR2CABV3QI8BGwGtgGlqvqRu6kAWA6MFpFOIhINnAOkteUFrBD4ABGJBV4HrlfVMrfzAKhqraoOBVKBHE/z1FUici5QrKqL3M7ShFGqOgyYAEz3dEe6LQwYBjylqicAB4Bb3Y30X56uqvOB19zOAiAiHYALgJ5ANyBGRK5wNxWo6irgfmA29d1CS4GatryGFQKXefrgXwf+o6pvuJ2nMU9XwjxgvLtJABgFnO/pj38ZOENE/u1upHqqutXzsRh4k/r+XLcVAUWHtOZmUF8YfMUEYLGq7nA7iMdYYKOqlqhqNfAGcLLLmQBQ1edVdZiqjqa+m7vNxgfACoGrPIOyzwOrVPURt/M0EJEkEUnwPG9H/R/IaldDAar6G1VNVdUe1HcpzFVV19+xiUiMZ7AfT9fLOOqb865S1e1AoYhker50JuDqRIRGpuAj3UIem4GRIhLt+ds8k/pxO9eJSGfPx3TgYtr45xbWli/mq0TkJeB0IFFEioA7VfV5d1MB9e9wfwQs8/THA/xWVd93LxIAXYF/eWZ0hACvqqrPTNX0QcnAm/X/OwgDXlTVWe5G+s61wH883TAbgKtczgOAp6/7LODnbmdpoKrfiMgMYDH1XS/f4jt3GL8uIp2AamC6qu5pyxcPiumjxhhjjsy6howxJshZITDGmCBnhcAYY4KcFQJjjAlyVgiMMSbIWSEwxpggZ4XAGGOCnBUCY4wJcv8POHoi9zmVueAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for rm in rmses:\n",
    "    plt.plot(range(1, len(rm) + 1), rm)\n",
    "    plt.ylabel(\"Error\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def removeFromGroup(x):\n",
    "    loc = 0\n",
    "    y = {}\n",
    "    for a in x:\n",
    "        for b in x[a]:\n",
    "            y[loc] = [ b ]\n",
    "            loc += 1\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: scikit-elm in /home/yash/anaconda3/lib/python3.8/site-packages (0.21a0)\r\n",
      "Requirement already satisfied: numpy in /home/yash/anaconda3/lib/python3.8/site-packages (from scikit-elm) (1.19.2)\r\n",
      "Requirement already satisfied: scipy in /home/yash/anaconda3/lib/python3.8/site-packages (from scikit-elm) (1.5.2)\r\n",
      "Requirement already satisfied: scikit-learn in /home/yash/anaconda3/lib/python3.8/site-packages (from scikit-elm) (0.23.2)\r\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /home/yash/anaconda3/lib/python3.8/site-packages (from scikit-learn->scikit-elm) (2.1.0)\r\n",
      "Requirement already satisfied: joblib>=0.11 in /home/yash/anaconda3/lib/python3.8/site-packages (from scikit-learn->scikit-elm) (0.17.0)\r\n"
     ]
    }
   ],
   "source": [
    "!pip install scikit-elm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 18/18 [14:10<00:00, 47.22s/it]\n"
     ]
    }
   ],
   "source": [
    "## Data Structure\n",
    "    # Cycles[battery][param][cycle]\n",
    "    # Cycles[battery][Capacity][cycle]\n",
    "\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from tqdm import tqdm    \n",
    "    \n",
    "Cycles = {}\n",
    "params = ['Temperature_measured', 'Voltage_measured', 'Voltage_load', 'Time']\n",
    "\n",
    "# remove batteries from group\n",
    "bs_compt = BSSS\n",
    "bs_compt = removeFromGroup(bs_compt)\n",
    "\n",
    "\n",
    "depths_test = []\n",
    "depths_train = []\n",
    "\n",
    "for bs_cmpt in bs_compt:\n",
    "    depths_test.append([bs_compt[bs_cmpt][0]])\n",
    "    depths_train.append([bs_compt[bs_cmpt][0]])\n",
    "    \n",
    "# iterate over the depths\n",
    "for depth in tqdm(range(3, 21)):\n",
    "\n",
    "    rmses = []\n",
    "    rmses_train = []\n",
    "\n",
    "    for bs_cmpt in bs_compt:\n",
    "        rmses.append([bs_compt[bs_cmpt][0]])\n",
    "        rmses_train.append([bs_compt[bs_cmpt][0]])\n",
    "\n",
    "    \n",
    "    # iterate over the merge hyper parameter\n",
    "    for xx in range(1, 41):\n",
    "\n",
    "        results = Texttable()\n",
    "        results.add_row([ 'Compatible Batteries', 'Cycles', 'MAE', 'RMSE', 'R2 Score' ])\n",
    "\n",
    "        loc = 0\n",
    "\n",
    "        # iterate over all the battery sets\n",
    "        for bs_cmpt in bs_compt:\n",
    "            # getting data for a given set\n",
    "\n",
    "            # y contains RUL after current cycle\n",
    "            # model will train for y\n",
    "            y = []\n",
    "\n",
    "            bs = bs_compt[bs_cmpt]\n",
    "            for b in bs:\n",
    "                Cycles[b] = {}\n",
    "                for param in params:\n",
    "                    Cycles[b][param] = []\n",
    "                    for j in range(datas[b].size):\n",
    "                        if types[b][j] == 'discharge':\n",
    "                            Cycles[b][param].append(datas[b][j][param][0][0][0])\n",
    "\n",
    "                cap = []\n",
    "                for j in range(datas[b].size):\n",
    "                    if types[b][j] == 'discharge':\n",
    "                        cap.append(datas[b][j]['Capacity'][0][0][0][0])\n",
    "                Cycles[b]['Capacity'] = np.array(cap)\n",
    "                Cycles[b]['count'] = len(Cycles[b][params[0]])\n",
    "\n",
    "                effective_cycle_count = 0\n",
    "                for x in Cycles[b]['Capacity']:\n",
    "                    if (x < 1.4):\n",
    "                        break\n",
    "                    effective_cycle_count += 1\n",
    "\n",
    "                for i in range(len(Cycles[b]['Capacity'])):\n",
    "                    if (i < effective_cycle_count):\n",
    "                        y.append((1 - ((i + 1) / effective_cycle_count)) * 100)\n",
    "                    else:\n",
    "                        y.append(0)\n",
    "\n",
    "            # preparing data for regression model\n",
    "            temperature_measured = []\n",
    "            voltage_measured = []\n",
    "            voltage_load = []\n",
    "            capacity = []\n",
    "\n",
    "            for b in bs:        \n",
    "                for c in Cycles[b]['Capacity']:\n",
    "                    capacity.append(c)\n",
    "                for i in range(Cycles[b]['count']):\n",
    "                    temperature_measured.append(getTemperatureMeasuredCritical(Cycles[b]['Temperature_measured'][i], Cycles[b]['Time'][i]))\n",
    "                    voltage_measured.append(getVoltageMeasuredCritical(Cycles[b]['Voltage_measured'][i], Cycles[b]['Time'][i]))\n",
    "                    voltage_load.append(getVoltageLoadCritical(Cycles[b]['Voltage_load'][i], Cycles[b]['Time'][i]))\n",
    "\n",
    "            # creating the model\n",
    "            X = []\n",
    "            for i in range(len(temperature_measured)):\n",
    "                X.append(np.array([temperature_measured[i], voltage_measured[i], voltage_load[i], capacity[i]]))\n",
    "        #     X.append(np.array(capacity))\n",
    "            X = np.array(X)\n",
    "            y = np.array(y)\n",
    "\n",
    "            # merge cycles\n",
    "            X, y = merge(X, y, xx)\n",
    "\n",
    "            # creating train test split\n",
    "            X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state = 0)\n",
    "\n",
    "            ############## ------------------ MODEL ------------------- ####################\n",
    "\n",
    "            ### Random Forest Regressor ### (BEST)\n",
    "            regressor = RandomForestRegressor(max_depth = depth, random_state= 0)\n",
    "            regressor.fit(X_train, y_train)\n",
    "            y_pred = regressor.predict(X_test)\n",
    "            y_pred_train = regressor.predict(X_train)\n",
    "\n",
    "            ############# ----------------- MODEL -------------------- #####################\n",
    "\n",
    "            # model evaluation\n",
    "            diff = 0\n",
    "            total = 0\n",
    "            rmse = 0\n",
    "            for i in range(len(y_test)):\n",
    "                diff += abs(y_test[i] - y_pred[i])\n",
    "                rmse += ((y_test[i] - y_pred[i]) * (y_test[i] - y_pred[i]))\n",
    "                total += y_test[i]\n",
    "            diff /= len(y_test)\n",
    "            total /= len(y_test)\n",
    "            rmse = math.sqrt(rmse / len(y_test))\n",
    "            accuracy2 = ((total - diff) / total) * 100\n",
    "            accuracy = r2_score(y_test, y_pred)\n",
    "\n",
    "            # Adding evaluation to result array to print in a table\n",
    "    #         results.add_row([ str(bs), str(Cycles[bs[0]]['count']), diff, rmse, accuracy, accuracy2 ])\n",
    "\n",
    "            rmses[loc].append(rmse)\n",
    "\n",
    "            #### adding rmses of the train\n",
    "            rmse_train = 0\n",
    "            for i in range(len(y_train)):\n",
    "                rmse_train += ((y_train[i] - y_pred_train[i]) * (y_train[i] - y_pred_train[i]))\n",
    "            rmse_train = math.sqrt(rmse_train / len(y_train))\n",
    "            rmses_train[loc].append(rmse_train)\n",
    "            loc += 1\n",
    "\n",
    "        # printing results\n",
    "    #     print(f'Evaluation: Clubbing Compatible Batteries for cycle param: {xx}\\n{results.draw()}')\n",
    "    \n",
    "    ## setting errors for depths\n",
    "    for i in range(len(rmses)):\n",
    "        \n",
    "        rm = rmses[i]\n",
    "        rm_train = rmses_train[i]\n",
    "\n",
    "        mn = 100000\n",
    "        mn_train = 100000\n",
    "        loc, loc_train = -1, -1\n",
    "        for j in range(1, len(rm)):\n",
    "            if (mn > rm[j]):\n",
    "                mn = rm[j]\n",
    "                loc = j\n",
    "\n",
    "            if (mn_train > rm_train[j]):\n",
    "                mn_train = rm_train[j]\n",
    "                loc_train = j\n",
    "        \n",
    "        depths_test[i].append(mn)\n",
    "        depths_train[i].append(mn_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n",
      "Minima Test: 0.45944155427729894, Depth: 6\n",
      "Minima Train: 0.26002986467389244, Depth: 4\n",
      "Minima Test: 0.6711554317677308, Depth: 4\n",
      "Minima Train: 0.4127800071535105, Depth: 3\n",
      "Minima Test: 0.29760281341004335, Depth: 6\n",
      "Minima Train: 0.24190013292878176, Depth: 6\n",
      "Minima Test: 0.5913611058246204, Depth: 3\n",
      "Minima Train: 0.3407219138897158, Depth: 5\n"
     ]
    }
   ],
   "source": [
    "%matplotlib\n",
    "for i in range(len(depths_test)):\n",
    "    rm = depths_test[i]\n",
    "    rm_train = depths_train[i]\n",
    "    \n",
    "    mn = 100000\n",
    "    mn_train = 100000\n",
    "    loc, loc_train = -1, -1\n",
    "    for i in range(1, len(rm)):\n",
    "        if (mn > rm[i]):\n",
    "            mn = rm[i]\n",
    "            loc = i\n",
    "        \n",
    "        if (mn_train > rm_train[i]):\n",
    "            mn_train = rm_train[i]\n",
    "            loc_train = i\n",
    "            \n",
    "    print(f\"Minima Test: {mn}, Depth: {loc}\")\n",
    "    print(f\"Minima Train: {mn_train}, Depth: {loc_train}\")\n",
    "    \n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(range(3, len(rm) + 2), rm[1:])\n",
    "    ax.plot(range(3, len(rm) + 2), rm_train[1:])\n",
    "    plt.ylabel(rm[0])\n",
    "    plt.xlabel(\"Random Forest Regressor Max-Depth (Hyper-param)\")\n",
    "    plt.show()"
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}