SAUBER Neckartor SARIMA Test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<function get_engine at 0x7fa420a2df28>\n"
     ]
    }
   ],
   "source": [
    "#Import database setup\n",
    "import setup_env"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# SARIMA example\n",
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from statsmodels.graphics.tsaplots import plot_acf, plot_pacf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grid search sarima hyperparameters\n",
    "from math import sqrt\n",
    "from multiprocessing import cpu_count\n",
    "from joblib import Parallel\n",
    "from joblib import delayed\n",
    "from warnings import catch_warnings\n",
    "from warnings import filterwarnings\n",
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "from sklearn.metrics import mean_squared_error\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import isnan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DB Engine läuft\n",
      "Engine(postgresql://postgres:***@192.168.10.184:5430/here)\n"
     ]
    }
   ],
   "source": [
    "engine = setup_env.get_database()\n",
    "print(engine)\n",
    "try:\n",
    "    con = engine.raw_connection()\n",
    "    con.cursor().execute(\"SET SCHEMA '{}'\".format(\"here_traffic\"))\n",
    "except:\n",
    "    print(\"Con: DB Verbindung prüfen!\") \n",
    "    exit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows = \"*\"\n",
    "sel_link_id = 563933733\n",
    "sel_min_confidence = 0\n",
    "sel_max_weekday = 8\n",
    "sel_func_classes = ('4','3')\n",
    "sel_dir_travel = 'F'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql_query = f\"\"\"\n",
    "    SELECT tr.*  \n",
    "    FROM here_traffic.stuttgart_traffic tr\n",
    "    JOIN here_streets.fc_streets_all_2018q3 st on tr.link_id = st.link_id\n",
    "    WHERE tr.link_id = {sel_link_id}\n",
    "    AND tr.confidence > {sel_min_confidence}\n",
    "    AND tr.weekday_n < {sel_max_weekday}\n",
    "    --AND st.func_class in {sel_func_classes}\n",
    "    AND tr.dir_travel = 'F'\n",
    "    LIMIT 100000\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 225 ms, sys: 47.7 ms, total: 273 ms\n",
      "Wall time: 36.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "pd_read = pd.read_sql_query(sql_query, con)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(\n",
    "    pd_read,\n",
    "    columns=[\n",
    "        \"id_pk\",\n",
    "        \"link_id\",\n",
    "        \"dir_travel\",\n",
    "        \"mean_kmh\",     #2\n",
    "        \"datum_zeit\",   #5\n",
    "        \"weekday_n\",    #6\n",
    "        \"epoch_60\",     #7\n",
    "        \"confidence\",\n",
    "        \"count_n\"\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop_duplicates(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.set_index('datum_zeit', inplace=True, drop=True)\n",
    "df.sort_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df.to_pickle('link_id_'+str(sel_link_id))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df4 = df[(df.index  > '2018-05-01 00:00:00') & (df.index < '2018-09-1 00:00:00')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "fill_index = pd.date_range('2018-05-01 00:00:00', '2019-04-30 23:00:00', freq='1H')\n",
    "df = df.reindex(fill_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['link_id']=sel_link_id\n",
    "df['dir_travel']=df['dir_travel'].fillna(method='backfill')\n",
    "df['epoch_60']=df.index.hour\n",
    "df['weekday_n']=df.index.dayofweek+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"hourweek\"] = (df[\"weekday_n\"]-1)*24+df[\"epoch_60\"]\n",
    "df[\"hourweek\"] = pd.to_numeric(df.hourweek, errors='coerce')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_mean = df['count_n'].groupby(df['hourweek']).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, row in df.iterrows():\n",
    "    m_idx = row['hourweek']\n",
    "    if isnan(row['count_n']):\n",
    "        row['count_n'] = count_mean[m_idx]\n",
    "        df['count_n'].at[index] = count_mean[m_idx].round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 5.55 s, sys: 21 ms, total: 5.57 s\n",
      "Wall time: 5.64 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for index, row in df.iterrows():\n",
    "    m_idx = row['epoch_60']\n",
    "    if isnan(row['count_n']):\n",
    "        row['count_n'] = count_mean[m_idx]\n",
    "        df['count_n'].at[index] = count_mean[m_idx].round()\n",
    "        #print(index,row['count_n'],row['epoch_60'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5hddX3v8fd378mEXJmETEKuJECARJQE0gRa9VCBU+JR6NN6LLFVtCjtOcLRai/eDrWeY0tPn15QqUeriHIUpHpao9KiRjlaC5EBApLEkBACuWdIMuSemdn7e/5Ya8+sWdl7Zs/sNXuv2evzep48mbX2mr1+sy/f3299f5dl7o6IiDS/XKMLICIi9aGALyKSEQr4IiIZoYAvIpIRCvgiIhmhgC8ikhEK+CJ1YmY7zOzaEf7u68xsS9JlkmxRwJfUMLNHzOywmY0fxu+4mV04muVqhPjf5e4/cfeLG1kmGfsU8CUVzGwh8DrAgRsaWpghmFlLNftE0kYBX9LiHcBjwL3AzaWdYav/3ZHtd5rZv4U//zjc/bSZHTOz3wr3v8fMtpnZITNba2ZzIr//KjP7fvjYfjP7SLh/vJn9nZntCf/9XelKw8yuNrNdZvYnZrYP+FK5feGxbzKzDWbWZWb/bmavKffHmtlKM3s0PG6vmX3GzFor/V2l80V+f0n42nSZ2UYzuyHy2L1mdreZfdfMjprZejO7YGRvizQTBXxJi3cAXw3//ZqZzRrqF9z99eGPl7n7ZHf/upm9AfgL4K3AbOBF4AEAM5sC/AD4V2AOcCGwLnyOjwJXAsuAy4CVwMcipzsXmA6cB9xabp+ZLQfuAX4POAf4HLC2QoqqAPwBMAO4CrgG+K+V/q7oL5rZOODbwPeAmcDtwFfNLJryuQn4M2AasA34ZNkXUTJFAV8azsxeSxA0H3T3J4DngbeN8Ol+G7jH3Z9099PAh4GrwpTRm4B97v7X7n7K3Y+6+/rI733C3Q+4eydBsHx75HmLwJ+6+2l3P1lh363A59x9vbsX3P3LwGmCimQAd3/C3R9z915330FQOfyHKv/GK4HJwJ3u3u3uPwS+A6yJHPNP7v4zd+8lqESXVfnc0sQU8CUNbga+5+4vh9tfI5LWGaY5BK16ANz9GHAQmAvMJ6hMhvy98Oc5ke1Odz8V+534vvOAD4Zpli4z6wrPOSf2e5jZRWb2HTPbZ2ZHgD8naO1XYw6w092LsfLOjWzvi/x8gqCCkIxTR5M0lJlNIEi/5MNcOMB4oM3MLgOOAxMjv3LuEE+5hyDwlp5/EkF6ZTewkyDVMdjvbQy3F4T7SsotKxvftxP4pLtXkz75LPAUsMbdj5rZ+4G3VPF7pbLON7NcJOgvAJ6r8vclo9TCl0b7dYJ89lKCtMMyYAnwE4K8/gbgN8xsYjhM8ZbY7+8Hzo9s3w+8y8yWhbnzPwfWh2mT7wCzzez9YSftFDNbFfm9j5lZu5nNAO4A/s8w/5Z/AH7fzFZZYJKZ/aew7yBuCnAEOGZmlwD/ZYi/K2o9Qav9j81snJldDbyZsK9CpBIFfGm0m4EvuftL7r6v9A/4DEFe/W+BboIA+GWCfHTUx4EvhymUt7r7D4D/DnwT2AtcQNiqd/ejwHUEwXEfsBX41fB5/ifQATwD/Bx4MtxXNXfvAN4Tlv0wQWfpOysc/ocE/RRHCSqKr8ceH/B3xc7THf4Nq4GXgb8H3uHuvxhOeSV7TDdAERHJBrXwRUQyQgFfRCQjFPBFRDJCAV9EJCNSOw5/xowZvnDhwkYXQ0RkTHniiSdedvf2co+lNuAvXLiQjo6ORhdDRGRMMbMXKz2mlI6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGJBLwzeweMztgZs9WeNzM7FPhbeeeMbPLkzhvOYWis27zfj61bivrNu+nUNRaQSIikNywzHsJVgj8SoXHVwOLw3+rCNYCX1Xh2BErFJ23f3E9G3Z2cbK7wITWPMvmt3HfLavI5yzp04mIjCmJtPDd/cfAoUEOuRH4igceI7i5xewkzh31yJYDbNjZxYnuAg6c6C6wYWcXj2w5kPSpRETGnHrl8OcS3A2oZBcDb8cGgJndamYdZtbR2dk57JNs3HOEk92FAftOdhfYtOfIsJ9LRKTZpKrT1t0/7+4r3H1Fe3vZmcGDetWcqUxozQ/YN6E1z9I5U5MqoojImFWvgL+b4GbOJfPCfYm6+uKZLJvfhhW6wYtMDHP4V188M+lTiYiMOfUK+GuBd4Sjda4EXnH3vUmfJJ8z7rtlFe1bv03brp/y6TXL1WErIhJKZJSOmd0PXA3MMLNdwJ8C4wDc/X8DDwFvJLjH5wngXUmct5x8zpjYtZ2JXdu5Zsms0TqNiMiYk0jAd/c1QzzuwHuTOJeIiIxMqjptRURk9Cjgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZEQiAd/MrjezLWa2zcw+VObxBWb2IzN7ysyeMbM3JnFeERGpXs0B38zywN3AamApsMbMlsYO+xjwoLsvB24C/r7W84qIyPAk0cJfCWxz9+3u3g08ANwYO8aBqeHPZwN7EjiviIgMQ0sCzzEX2BnZ3gWsih3zceB7ZnY7MAm4NoHziojIMNSr03YNcK+7zwPeCNxnZmec28xuNbMOM+vo7OysU9FERLIhiYC/G5gf2Z4X7ou6BXgQwN0fBc4CZsSfyN0/7+4r3H1Fe3t7AkUTEZGSJAL+48BiM1tkZq0EnbJrY8e8BFwDYGZLCAK+mvAiInVUc8B3917gNuBhYDPBaJyNZvYJM7shPOyDwHvM7GngfuCd7u61nltERKqXRKct7v4Q8FBs3x2RnzcBv5LEuUREZGQ001ZEJCMU8EVEMkIBX0QkIxTwRUQyQgFfRCQjFPBFRDJCAV9EJCMU8EVEMkIBX0QkIxTwRUQyQgFfRCQjFPBFRDJCAV9EJCMU8EVEMkIBX0QkIxTwRUQyQgFfRCQjFPBFRDJCAV9EJCMU8EVEMkIBX0QkIxTwRUQyQgFfRCQjFPBFRDJCAV9EJCMSCfhmdr2ZbTGzbWb2oQrHvNXMNpnZRjP7WhLnFRGR6rXU+gRmlgfuBq4DdgGPm9lad98UOWYx8GHgV9z9sJnNrPW8IiIyPEm08FcC29x9u7t3Aw8AN8aOeQ9wt7sfBnD3AwmcV0REhiGJgD8X2BnZ3hXui7oIuMjMfmpmj5nZ9eWeyMxuNbMOM+vo7OxMoGgiIlJSr07bFmAxcDWwBvgHM2uLH+Tun3f3Fe6+or29vU5FExHJhiQC/m5gfmR7Xrgvahew1t173P0F4DmCCkBEROokiYD/OLDYzBaZWStwE7A2dsw/E7TuMbMZBCme7QmcW0REqlRzwHf3XuA24GFgM/Cgu280s0+Y2Q3hYQ8DB81sE/Aj4I/c/WCt5xYRkerVPCwTwN0fAh6K7bsj8rMDHwj/iYhIA2imrYhIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkRCI3MReppFB0HtlygI17jvCqOVO5+uKZ5HPW6GKJZJICvoyaQtF5+xfXs2FnFye7C0xozbNsfhv33bJKQV+kAZTSkVHzyJYDbNjZxYnuAg6c6C6wYWcXj2w50OiiiWRS0wf8QtFZt3k/n1q3lXWb91MoeqOLlBkb9xzhZHdhwL6T3QU27TnSoBKJZFsiKR0zux64C8gDX3D3Oysc95vAN4BfcveOJM49GKUUGutVc6YyoTXPiUjQn9CaZ+mcqQ0slUh21dzCN7M8cDewGlgKrDGzpWWOmwK8D1hf6zmrpZRCY1198UyWzW/DCt3gRSaGFe7VF89sdNFEMimJlM5KYJu7b3f3buAB4MYyx/0P4C+BUwmcsypKKTRWPmfcd8sq2rd+m7ZdP+XTa5br6kqkgZII+HOBnZHtXeG+PmZ2OTDf3b872BOZ2a1m1mFmHZ2dnTUXrJRSiFJKob7yOWNi13badj/GNUtmKdiLNNCod9qaWQ74G+CDQx3r7p939xXuvqK9vb3mcyulICLSL4mAvxuYH9meF+4rmQJcCjxiZjuAK4G1ZrYigXMPSikFEZF+SYzSeRxYbGaLCAL9TcDbSg+6+yvAjNK2mT0C/GE9RulAf0phYtd2rlkyqx6nFBFJpZoDvrv3mtltwMMEwzLvcfeNZvYJoMPd19Z6DmkuWm5BpDESGYfv7g8BD8X23VHh2KuTOKeMTZob0XiqcLNLa+lIosoFk6jo3AgYODdCKbfRpwo32xTwJTGVgoljGMGSFoPNjVDAH32qcLOt6dfSkfqpNLP5ZNuivmM0N6KxNBkxneq15pda+JKYSsGkddJMJnZtB/rnRjz63F4818LE8eM0N6KOtL5R+tQzzaYWviSmUuu99Xj/2kWaG9FYmoyYPvVc80sBXxJTKZhM6HphwHFabqFxVOGmTz3TbAr4kphKwaTUYSvpoAo3XerZr6WAL4lSMBEZnnqm2RTwRUQaqJ5pNo3SEWlyQ02G08zbxqvXml8K+CJNbKjJcJp52xhDVcKjJXMBX60ZyZJKM2snty1iYtd2zbxtgGpmpI+WTAV8tWYka4aaDKelLupvqEp4NGWq01Y3NZesGWoynJa6qL9KlWz3pNFP62Qq4GsdEcmaoSbDaeZt/VUzI320ZCrgqzUjWTPUZDjNvK2/amekj4ZMBXy1ZiSLhpoMp8ly9dXIGemZ6rQtvdBX/cYtdE+ayV9/7A80SkdE6q5R99rOVMAH3dQ8jTRUVqQ+MhfwJVm1TiDRUFmR+lHAlxFLYgKJJv6I1E+mOm0lWdXc0nAoGiqbTvW65Z7Ul1r4MmLV3NJwKLrlXvoozda81MKXEUtiAomGyiav1ta5ZqQ3r0Ra+GZ2PXAXkAe+4O53xh7/APBuoBfoBH7X3V9M4tzSOJVuSL77x9VPINFQ2WQl0a+i9XWaV80tfDPLA3cDq4GlwBozWxo77Clghbu/BvgG8L9qPa80XlITSDTxJzlJ9KtoRnry0tInkkQLfyWwzd23A5jZA8CNwKbSAe7+o8jxjwG/k8B5JQU0ryFdkuhXqXTlpjTbyDRyOeS4JHL4c4Gdke1d4b5KbgH+pdwDZnarmXWYWUdnZ2cCRRPJliT6VbS+TrKSuOpKSl07bc3sd4AVwF+Ve9zdP+/uK9x9RXt7ez2LJtIUklqYS2m25DRyOeS4JAL+bmB+ZHteuG8AM7sW+Chwg7ufTuC8IhLTyIW5pLxGLoccl0TAfxxYbGaLzKwVuAlYGz3AzJYDnyMI9hrbJTKK1DpPl0YuhxxXc8B3917gNuBhYDPwoLtvNLNPmNkN4WF/BUwG/tHMNpjZ2gpPJ5KaEQ0iSUjTVVci4/Dd/SHgodi+OyI/X5vEeUaDVmpMF83ylGaUltFsmV5aQcElfbSYmsjoyfTSCppCnj5aTE1k9GQ64Cu4DN9o59c1yzOd1K/SHDKd0tFKjcNTjxmDmuU5fLXehKaa51fqszlkOuAruAxPpfz65LZFVU/bH4oWUxueelTC6ldpHplO6WgK+fDUa8agxpFXrx7T9pX6bB6ZbuFDeoZLjQWVUmCNmDEogSQWSxuKUp/Dk+ah3pkP+FK9JNa/l2TVoxJW6rN6g/V3pEGmUzoyPGmaMSiBekzbV+qzemkf6q2AL8Oi/Hq61KsS1vtenbT3dyjgi4xxCsbpkfZ5JAr4knqa9NNPr0W6VUqxpaW/Q522kmqa9NMvTR2CaR6J0khpn0eigC+DGu1ZnEPRpJ9+g70W9aRKeKBy35G0DvVWwI9Ry6VfGm6+PFgnWNq+TKMtLR2CqoT7peE7MhwK+BFquQxUj6UUhqJJP/0Gey2+WcdyqBLul4bvyHCo0zYi7WNo6y0NN19OeydYPaXltUj7SJR6SsN3ZDgU8CPScsmcFmm4+bIm/fRLy2uRloonDdLwHRkOBfwItVwGSsvNl+PjzIHMDk1Mw5j7tFQ8aZCW70i1lMOP0JohA1UaYnbNVxsXYNXPkg7xRQdL8wOyNthhJN+RYtFxwMNDTvcWcIdiuOP46V4mjR+d0KyAH1HpzQMy+WGG9K0mqhEitSkWnaJ7X3A5drqXQtHpLTiOc+DIKYoOBXdO9xYBeOngCYrunOoJXvPtncdw4GS4/dz+o/zJN5/hF3uPcKqnyFnjclxy7lTu/M3XVPE96Q+M7mfuLe3zcE/ftgf7otul44LHwMOd0eAaP180LA84vw+vUdO7bws5tjBp/B/x+I5DHDnVA8Bj2w/2Pe+Rk8G+9S8cAuBoeMyTL3YBcOxULwDP7HqFVYumkxuFGKOAH1Ou5aIW5ehwD/pIeotFeouOu9N59DSO010ogsPeV07iTl/weXT7wbL9LI9uP8jF504pe45qygFlgkrf4/HA0X+c0x98+vbFgk8154qeL6pQdJ586TDbO49zfvskLl8wjWOng8CwYWdXX9lK+5586TAAR8Pg0bHjEA59AagUbErB5ee7XgHgRHew/Xzn8b5znw4D+u6ukwB0h+/B/iOnAegJt7+/cT+b9hzpe49O9hTZtPcIP9i0n8vPm1bmr2pCsYrnjO2UUMAfQtpblO5eNkBFW0jlgsvA5xj8HKXWYCkoFMNzHj/di0NfDv2VEz0U3ekpBF/8va+cpLcQtAwd2LTnCD2FIj2FYl9rZ8POoHVzIgxY2w4cA+BU+HrvePkE0B98pk1opbUl1xdcAFpbckyb0MrOQycH/0PGmGLR+fN/2cy2A8fo7i3S2pLjwpmTKRTB8AEVXzF8D073BK9LqSLoKcQ+HKNgx8HjfZVBSXdvkR0HjzdlwC8WnQ07u9hx8DgLz5nEsvltjS5S1TIR8EuXd0dP9dBbCAKSO+w8dILeogeXpg5b9h3F8b7gvmnPEX60pbNsi/KRLQc49+yzyp+vzJcr3rIbuG9gK7L/sTItyVgwT1K5D3IuZ32twad3Bq3B6KUnBIEfYNPeYDTTyViwLgWDV8IgX4tl89u4cOZkNr70MuRbGD+uhQtnTh5TX7pqbdjZxbYDx/oqt9O9RbYdOMa4cy6g9eC2Bpeu38JzJpWthBeeM6mBpRodlSrhtE60iksk4JvZ9cBdQB74grvfGXt8PPAV4ArgIPBb7r4jiXMP5uipHjbvPcrRMNA8u3tgQNp1OGgRli5NDx3vBqC30B+gZk89q+yH+dypEzhysne0/4S6qfRB/sjqJY0u2gC5nPGR1Uv4vfd9kMLkWdz2+7f2VUzNplLLOTd5FqQo4KsSTl8lXEnNwzLNLA/cDawGlgJrzGxp7LBbgMPufiHwt8Bf1nreakTTDSNV+jDTGwy7Gh8Gwmb7MEc/yE7/B7mUckmTXM5oPbiNCS/+lMvPm9aUwR76W85RrS058sf2N6hE5ZUq4cmb/pkJL/yE//aGxXxk9ZKmfF8qVcKFyY1P71YjiXH4K4Ft7r7d3buBB4AbY8fcCHw5/PkbwDVmNiY+DVn5MA+Wh5XGqNTYGHfw+UYX7QyqhNNVCVdiwx1+dMYTmL0FuN7d3x1uvx1Y5e63RY55NjxmV7j9fHjMy5Wed/p5S/y6j9wzojJteHoDAK9+9WUc7+5l66ZnAVi89FKAYW9X2tdMjp7qZXfXyQH9AmYwt20C+7b/Ahj561fN65nEczYjd+e5bdsh38qcObOZPD7Pts0bgdF/PfWenMndeenQSU6c7gEMyxkTxuU5vX87RnKv59Szxo24jA/+/i8/4e4ryj2WqoBvZrcCtwJMnn3BFW/80/tqKluh6BzvHp08exq+UEk+54VLXsVLh05ysieYBGIGE8blWTB9AmPkYixVr2eSzzkWlauozCwVr2etz+nuHDtd4HRPgfHj8n1/W5LSHPCvAj7u7r8Wbn8YwN3/InLMw+Exj5pZC7APaPdBTr5ixQrv6OioqWxHTvWwcfforIPz3rfdAMDdX1s7ou0kniPp56w0SmesSNvrmdRzjjWlAQDxTtyPrF7C7b8TZHvH0ntULHrdBwnUMvHKzCoG/CRG6TwOLDazRcBu4CbgbbFj1gI3A48CbwF+OFiwl8bI5YzLz5s2JsdOF4tO9zkXUpg8iydfPNx0nepjSWkAAC2tQLoHAAylVHkdW/rrkG/hUz/c2ld5jaXGUEnNAd/de83sNuBhgmGZ97j7RjP7BNDh7muBLwL3mdk24BBBpTDqpoxvYfmCNnoKxb7x9z1Fp6c3mNlZiP5z7xvjHh0DP/BvDf4vFJ2eGYvpnTyLp3d1sXx+W98l3bh8MBq3tN03G7f039j7jKRepS/lWBkbHVWu4hprgaWZBgAMVnmNxYZRIuPw3f0h4KHYvjsiP58C/nMS5xoOM+OscXnOGpcf+uAqlZZaOLb0RjzXwqfWbWXZ/DYmnzUOw1mxcDoAU84KXtqVi4LtUk7uyvPP6XuuqRNK+6YP2F65aDruzpTwd644bxqOMzl8zuUL2nCHyWe14A6vnnc2EORNS4sulVb4nBhuL5k9BXeY2Bpsl5YhmNjaguNcNGtyVX9/fO0RJ5gV5hCu0dK/XkupUu0pFPv+7ymcudxArSp9KcfK2OiSZmlNNtNErGabRZyJmbZJKi214PkguCRxh5t4h09wVWB9VwOlYWC5cEepAsuZgcHkyMp6pSuKs8PKoyXcbpsYlLclH2xPn9RKoeh0n3MB3ZNmsWFnV10WhXN3esKrrd6C01sMKoNS5QD9FUelK654hXHg6KmyX8qWqecyrms7Mya34sC4fPA6njO5FXdoCbenTQpfq3C79No5Tj5v4P0VeOn1mTQ+jzt9gfiscbkBV3atLRZ2fg+80nP3AVd70b+lWVqTpeGk8Ul8YzHN1kyVFyjgD1ulm6S0TpqZyluaVVK6Uulc/GY818Lt9z9Vl0XhzIzWFjtjLHMtXn9RO//4xK4zbv036VQnE1rzLJ41pW8fwEWzSlc3wfYl504dsB29/8Gk8Iro0rnBVVTpCuo184LgVapsly8IAnKpYrjivPJXenDm1V7pyq7gXn5Sj3vfVWDcGYuz9e3vryz7VpD04BxF9/AqDHqLRYrFYH+h0L+SZtH710wq/X5pTaVq1mT65K+/midePMwLLx9j0YzJXHHeNPI56/tslV63aAUK/RVo6b0aqkKN7stZ8ksGNVPlBQr4w1bpvqJpvcNNJZWuVNKyKNxwlO5jEF/RdPePX2h00Yal0mfrVXOmVhz2V7lPqPEpoGULzgyKpQqzlIasVIGWAupQFWp036pYBXrl+dNxhykTxoHDLy2chsOAVCkwIFVaUtp3xcJpPPh7V/HjrZ38Yu9RLpk9hdcvbiefs6pTk/H1s0prZ/VdvZYq13B/0Uevr08Bf5iaJbg0042oS/cxeGTLATbtOcLSOVN53eJ2XvvdRXRPmsW6zfvHxE1sKn22xkLZ08gsSIsagPWn7IZKlZbbt/rS2ay+dHZ9Cj6KFPCHqVxwafRdoEaiUmtyrN7OMZ8zrlkya8A9DOLpqrSP2qn02dJ9F+qrUHROtJ0/oLHQLO+BAv4IRIMLlP+ApF0ztyZHo2O9XuKfrWYyFr4njerbqhcF/BqpNZk+Y6ljvZlbk1Fj5XvSTH1b5SQ3VCKjBnxALNf3ATnZtqjRRTtDKbh0zb2KdZuD1f2uWTKL269ZzDVLZjVNoCmlq6LS2LEeDYJd836Z2+9/ird/cX3NS3qn0Vj5ngzWt9UMFPBrVOkD0j2p/3I1Hmgb8YXOUnAppasmtuYxguGWy+a3MaErXR3rlYLgI1vSVTEloZrvSRpUaiyM1b6tOAX8Gg3VmqwUaLt7i3WtBLIUXErpqk+vWc4HrruIT69Zzr3vWsnJtkUNrXTjmr01GTVWrroqNRbS2N8wEsrh16hc5+dl885m2zNG19yruOsHz52RE3zqpcPc+Jl/q2s+s5mGYVZjLIzaabaRUoMZK8OZm7lvCxTwaxb/gFxy7hS+9O87eHnxm/BcC5/9f8/TUxgYVE72FNnWeWzUR5FEOwR7C8XMBJe4tI7aaeaRUnFpnStRqdO8WUdKKeAnIPoBWbd5/4DgEg/2EKyo2RuvBIY5iqTSELfSvu9v3MeX/n1HX6v2H36yndaWHBNb800fXOLKXd2c6C5g51zS9/q9bnF7XUbLRN+3R7Yc4N53reQnWzubsjUZN9RV12XzzuZ42wX0TJpZlwqg2YdglqOAn7BywQX6g/yE1jwLpk/kpUMnBrS2g7VCcn055sE+7JW+LGbWt++2+5+ip1Dsq3hO9gRrtLzndeczLp9r+uASVS51kjM4fs5FYHlu+9qTtLbkOFLmi5+kwQJMM7YmB1Puqmv9C4coXvRmsHxd0m7NPgSzHAX8hJXNy47LDQi0r1vczju/9LO+S/mzxuWCgDPnl6rKMZf7oD75UnBzidK+07FFuABO9RQZl89x+zWLR+NPT6146qS1JUdPoQi5YE2Vkz3FoEIs88VPUhYDTCXlGkZFp+89qUea89tP78lUvxYo4CeuUl72fddeNKA1Hc1n9hSKfOHfXqg6x1zuy1IuwMdlJWcfF88fP995jG9t2DPo74xktEw8zRZPE/189yuZCzCVlGsYxSUxWS76nsTTnP/y7L4zlqhu9u+IAn7Cqu3lj+YzP7Vu67ByzOU6YMeHC0FFA3/OgjXgu3uLmcrZlxPvZ/nepv2DBpvSF/+bVT5/PF1TLk20YPrEzHacx1W66oqOlh1OmrOcM96TWJrzdG+RnAXfnax8RxTwR8Fwe/mHm2Mu1wFbyuHHh4f+7msX8Yu9RzOVsx9KPNiUUmo9BR9wVVauI7eSeLqmXJroxfDm8C8eOpG5jvO4SqPbRprmLCf+npS7Ci46rL70XC5on5yJ74gCfgoMN8dcqQMWKHtlcd3Scxvzh6VUpSGC0dEypX6WwUaRRNosaF8AAAbySURBVCuEcvnguFM9Ra6/9FwunXt2JkblDCXeMHrDklkjTnOWU2kARdTE1jxvvmxOZlJqCvgpMJIcc6UO2GYdP5y0cldh0e348Nr4KJL4VVe5fHDchNY8l849W+9RBUOlOYeb06905ZzlNKcCfkqMNMcso2OoUSTxq654PrhSmihLwaUWI72zXLSTtlj0srPgs5zmVMBPoWpzzAoeo6eaUSRx8XxwPE2UteBSi5EsxRDvpH3f1zdw2byzueumZWcE+KymORXwU6iaHLOCx+iqZhRJXLl8sNI3I1PtUgzR70C5eQ5P73qFnFnm5p5UYl7tnXjrbMWKFd7R0dHoYkiGFYo+5CiS+FVXM0/Lb6RS6/3R5/biuRYmjh/Hsvlt3Puulbz2Le+me9Isrn7zW/nWhj0DxvAY8IHrLspUwDezJ9x9RbnH1MIXqWCwUSS66qqvcq33+KqzWZxINVw1BXwzmw58HVgI7ADe6u6HY8csAz4LTAUKwCfd/eu1nFekEYYa2SOjp+zy3rFVZ7M4kWq4am3hfwhY5+53mtmHwu0/iR1zAniHu281sznAE2b2sLt31XhuEcmIcp3o5VadzdpEquGqNeDfCFwd/vxl4BFiAd/dn4v8vMfMDgDtgAK+iFSl3KidcqvOZm0i1XDVGvBnufve8Od9wKCvspmtBFqB5ys8fitwK8CCBQtqLJqINItKo3aiq84qhTO0IUfpmNkPgHKDVj8KfNnd2yLHHnb3aRWeZzbBFcDN7v7YUAXTKB0RGUp0JJVSOIGaRum4+7WDPPF+M5vt7nvDgF52GpyZTQW+C3y0mmAvIlKNZr4d4WjI1fj7a4Gbw59vBr4VP8DMWoF/Ar7i7t+o8XwiIjJCtQb8O4HrzGwrcG24jZmtMLMvhMe8FXg98E4z2xD+W1bjeUVEZJg001ZEpIkMlsOvtYUvIiJjhAK+iEhGpDalY2adwIs1PMUM4OWEijOaVM5kqZzJUjmTVY9ynufu7eUeSG3Ar5WZdVTKY6WJypkslTNZKmeyGl1OpXRERDJCAV9EJCOaOeB/vtEFqJLKmSyVM1kqZ7IaWs6mzeGLiMhAzdzCFxGRCAV8EZGMaLqAb2bXm9kWM9sW3oUrFczsHjM7YGbPRvZNN7Pvm9nW8P+yS0vXk5nNN7MfmdkmM9toZu9LY1nN7Cwz+5mZPR2W88/C/YvMbH34/n89XLyv4cwsb2ZPmdl3wu3UldPMdpjZz8P1rjrCfal638MytZnZN8zsF2a22cyuSls5zeziyNphG8zsiJm9v9HlbKqAb2Z54G5gNbAUWGNmSxtbqj73AtfH9pVuEbkYWBduN1ov8EF3XwpcCbw3fA3TVtbTwBvc/TJgGXC9mV0J/CXwt+5+IXAYuKWBZYx6H7A5sp3Wcv6quy+LjBVP2/sOcBfwr+5+CXAZweuaqnK6+5bwdVwGXEFwq9d/otHldPem+QdcBTwc2f4w8OFGlytSnoXAs5HtLcDs8OfZwJZGl7FMmb8FXJfmsgITgSeBVQSzGFvKfR4aWL55BF/uNwDfASyl5dwBzIjtS9X7DpwNvEA44CSt5YyV7T8CP01DOZuqhQ/MBXZGtneF+9JqWLeIrDczWwgsB9aTwrKGaZINBDfe+T7BrTO73L03PCQt7//fAX8MFMPtc0hnOR34npk9Ed5uFNL3vi8COoEvhSmyL5jZJNJXzqibgPvDnxtazmYL+GOWB1V+asbImtlk4JvA+939SPSxtJTV3QseXDLPA1YClzS4SGcwszcBB9z9iUaXpQqvdffLCVKi7zWz10cfTMn73gJcDnzW3ZcDx4mlRVJSTqDvBlA3AP8Yf6wR5Wy2gL8bmB/ZnhfuS6v94a0hS/f8LXuLyHozs3EEwf6r7v5/w92pLCuAu3cBPyJIjbSZWenWnWl4/38FuMHMdgAPEKR17iJ95cTdd4f/HyDIN68kfe/7LmCXu68Pt79BUAGkrZwlq4En3X1/uN3QcjZbwH8cWByOgGgluJRa2+AyDWbIW0TWm5kZ8EVgs7v/TeShVJXVzNrNrC38eQJBP8NmgsD/lvCwhpfT3T/s7vPcfSHB5/GH7v7bpKycZjbJzKaUfibIOz9Lyt53d98H7DSzi8Nd1wCbSFk5I9bQn86BRpez0R0ao9BB8kbgOYJ87kcbXZ5Iue4H9gI9BK2UWwhyueuArcAPgOkpKOdrCS4znwE2hP/emLayAq8BngrL+SxwR7j/fOBnwDaCy+jxjX5NI2W+GvhOGssZlufp8N/G0ncnbe97WKZlQEf43v8zMC2l5ZwEHATOjuxraDm1tIKISEY0W0pHREQqUMAXEckIBXwRkYxQwBcRyQgFfBGRjFDAFxHJCAV8EZGM+P8nfWwsUYCv/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3gcd33v8fd3VxfLlh35Il9kK3YuSmInxE4QdgKhuDhpE5omOcCBuC0ETopLCzzQ0tOTEA6H0hucPhRISXsIBBJSmpDSCyakhODgEiAxURwn+IItxbEjWbYlXxTZlqzL7vf8sbPyaLWr266sjebzeh4/3pmdmd93ZnY+M/ubXa25OyIiMvXFJrsAERE5OxT4IiIRocAXEYkIBb6ISEQo8EVEIkKBLyISEQp8kbPEzPaZ2bXjnPfNZra70DVJtCjwpWiY2WYzO25m5WOYx83swomsazJkrpe7P+XuF09mTfLap8CXomBmy4A3Aw7cNKnFjMDMSkYzTqTYKPClWLwXeAa4H7gtPTK46v/90PD7zOynweOfBKNfMLOTZvbuYPwHzKzJzI6Z2UYzqwnNf6mZPRE8d9jMPhGMLzezL5pZa/Dvi+l3Gma21sxazOx/mdkh4BvZxgXT3mhm28ysw8x+bmaXZ1tZM1ttZk8H0x00sy+bWVmu9Uq3F5p/ebBtOsxsh5ndFHrufjO7x8y+b2YnzGyLmV0wvt0iU4kCX4rFe4FvBf9+08wWjDSDu/9a8HClu1e6+7fN7K3A3wDvAhYB+4GHAcxsJvAj4AdADXAhsClYxl3AVcAqYCWwGvhkqLmFwBxgKbAh2zgzuwL4OvAHwFzgK8DGHF1UCeCPgXnA1cA64I9yrVd4RjMrBb4H/BCYD3wE+JaZhbt8bgX+HJgNNAF/lXUjSqQo8GXSmdk1pELzEXd/DngJ+J1xLu53ga+7+1Z37wHuBK4OuoxuBA65++fd/bS7n3D3LaH5PuPube7eTios3xNabhL4P+7e4+7dOcZtAL7i7lvcPeHuDwA9pE4kg7j7c+7+jLv3u/s+UieHt4xyHa8CKoHPunuvuz8JPAqsD03z7+7+C3fvJ3USXTXKZcsUpsCXYnAb8EN3PxIM/zOhbp0xqiF1VQ+Au58EjgKLgVpSJ5MR5wse14SG2939dMY8meOWAh8Pulk6zKwjaLMmYz7M7CIze9TMDplZJ/DXpK72R6MGaHb3ZEa9i0PDh0KPu0idICTidKNJJpWZVZDqfokHfeEA5UCVma0ETgHTQ7MsHGGRraSCN738GaS6Vw4AzaS6Ooabb0cwfG4wLi3bn5XNHNcM/JW7j6b75B+B54H17n7CzD4GvHMU86VrrTWzWCj0zwX2jHJ+iShd4ctku4VUf/YKUt0Oq4DlwFOk+vW3AW83s+nBxxRvz5j/MHB+aPgh4P1mtiroO/9rYEvQbfIosMjMPhbcpJ1pZmtC833SzKrNbB7wKeCfxrguXwU+aGZrLGWGmf1WcO8g00ygEzhpZpcAfzjCeoVtIXXV/mdmVmpma4HfJrhXIZKLAl8m223AN9z9FXc/lP4HfJlUv/oXgF5SAfgAqf7osE8DDwRdKO9y9x8B/xv4V+AgcAHBVb27nwCuIxWOh4BG4NeD5fwl0AC8CPwS2BqMGzV3bwA+ENR+nNTN0vflmPxPSd2nOEHqRPHtjOcHrVdGO73BOtwAHAH+AXivu/9qLPVK9Jh+AEVEJBp0hS8iEhEKfBGRiFDgi4hEhAJfRCQiivZz+PPmzfNly5ZNdhkiIq8pzz333BF3r872XNEG/rJly2hoaJjsMkREXlPMbH+u59SlIyISEQp8EZGIUOCLiESEAl9EJCIU+CIiEVGQwDezr5tZm5ltz/G8mdndwc/OvWhmVxai3WwSSWfTrsPcvamRTbsOk0jqbwWJiEDhPpZ5P6m/EPjNHM/fANQF/9aQ+lvga3JMO26JpPOe+7awrbmD7t4EFWVxVtVW8eDta4jHrNDNiYi8phTkCt/dfwIcG2aSm4FvesozpH7cYlEh2g7bvLuNbc0ddPUmcKCrN8G25g42724rdFMiIq85Z6sPfzGpXwNKa2Hwz7EBYGYbzKzBzBra29vH3MiO1k66exODxnX3JtjZ2jnmZYmITDVFddPW3e9193p3r6+uzvrN4GFdWjOLirL4oHEVZXFW1MwqVIkiIq9ZZyvwD5D6Mee0JcG4glp78XxW1VZhiV7wJNODPvy1F88vdFMiIq85ZyvwNwLvDT6tcxXwqrsfLHQj8Zjx4O1rqG78HlUtP+Pv11+hG7YiIoGCfErHzB4C1gLzzKwF+D9AKYC7/z/gMeBtpH7jswt4fyHazSYeM6Z37GV6x17WLV8wUc2IiLzmFCTw3X39CM878KFCtCUiIuNTVDdtRURk4ijwRUQiQoEvIhIRCnwRkYhQ4IuIRIQCX0QkIhT4IiIRocAXEYkIBb6ISEQo8EVEIkKBLyISEQp8EZGIUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhEKPBFRCJCgS8iEhEFCXwzu97MdptZk5ndkeX5c83sx2b2vJm9aGZvK0S7IiIyenkHvpnFgXuAG4AVwHozW5Ex2SeBR9z9CuBW4B/ybVdERMamEFf4q4Emd9/r7r3Aw8DNGdM4MCt4fA7QWoB2RURkDAoR+IuB5tBwSzAu7NPA75lZC/AY8JFsCzKzDWbWYGYN7e3tBShNRETSztZN2/XA/e6+BHgb8KCZDWnb3e9193p3r6+urj5LpYmIREMhAv8AUBsaXhKMC7sdeATA3Z8GpgHzCtC2iIiMUiEC/1mgzszOM7MyUjdlN2ZM8wqwDsDMlpMKfPXZiIicRXkHvrv3Ax8GHgd2kfo0zg4z+4yZ3RRM9nHgA2b2AvAQ8D5393zbFhGR0SspxELc/TFSN2PD4z4VerwTeFMh2hIRkfHRN21FRCJCgS8iEhEKfBGRiFDgi4hEhAJfRCQiFPgiIhGhwBcRiQgFvohIRCjwRUQiQoEvIhIRCnwRkYhQ4IuIRIQCX0QkIhT4IiIRocAXEYkIBb6ISEQo8EVEIkKBLyISEQUJfDO73sx2m1mTmd2RY5p3mdlOM9thZv9ciHZFRGT08v5NWzOLA/cA1wEtwLNmtjH4Hdv0NHXAncCb3P24mc3Pt10RERmbQlzhrwaa3H2vu/cCDwM3Z0zzAeAedz8O4O5tBWhXRETGoBCBvxhoDg23BOPCLgIuMrOfmdkzZnZ9tgWZ2QYzazCzhvb29gKUJiIiaWfrpm0JUAesBdYDXzWzqsyJ3P1ed6939/rq6uqzVJqISDQUIvAPALWh4SXBuLAWYKO797n7y8AeUicAERE5SwoR+M8CdWZ2npmVAbcCGzOm+Q9SV/eY2TxSXTx7C9C2iIiMUt6B7+79wIeBx4FdwCPuvsPMPmNmNwWTPQ4cNbOdwI+B/+nuR/NtW0RERi/vj2UCuPtjwGMZ4z4VeuzAnwT/RERkEuibtiIiEaHAFxGJCAW+iEhEKPBFRCJCgS8iEhEKfBGRiFDgi4hEhAJfRCQiFPgiIhGhwBcRiQgFvohIRCjwRUQiQoEvIhIRCnwRkYhQ4IuIRIQCX0QkIhT4IiIRocAXEYmIggS+mV1vZrvNrMnM7hhmuneYmZtZfSHaFRGR0cs78M0sDtwD3ACsANab2Yos080EPgpsybdNEREZu0Jc4a8Gmtx9r7v3Ag8DN2eZ7i+AzwGnC9CmiIiMUSECfzHQHBpuCcYNMLMrgVp3/34B2hMRkXGY8Ju2ZhYD/g74+Cim3WBmDWbW0N7ePtGliYhESiEC/wBQGxpeEoxLmwlcBmw2s33AVcDGbDdu3f1ed6939/rq6uoClCYiImmFCPxngTozO8/MyoBbgY3pJ939VXef5+7L3H0Z8Axwk7s3FKBtEREZpbwD3937gQ8DjwO7gEfcfYeZfcbMbsp3+SIiUhglhViIuz8GPJYx7lM5pl1biDZFRGRs9E1bEZGIUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhEKPBFRCJCgS8iEhEKfBGRiFDgi4hEhAJfRCQiFPgiIhGhwBcRiQgFvohIRCjwRUQiQoEvIhIRCnwRkYhQ4IuIRIQCX0QkIgoS+GZ2vZntNrMmM7sjy/N/YmY7zexFM9tkZksL0a6IiIxe3oFvZnHgHuAGYAWw3sxWZEz2PFDv7pcD3wH+b77tiojI2BTiCn810OTue929F3gYuDk8gbv/2N27gsFngCUFaFdERMagEIG/GGgODbcE43K5HfjPbE+Y2QYzazCzhvb29gKUJiIiaWf1pq2Z/R5QD/xttufd/V53r3f3+urq6rNZmojIlFdSgGUcAGpDw0uCcYOY2bXAXcBb3L2nAO2KiMgYFOIK/1mgzszOM7My4FZgY3gCM7sC+Apwk7u3FaBNEREZo7wD3937gQ8DjwO7gEfcfYeZfcbMbgom+1ugEvgXM9tmZhtzLE5ERCZIIbp0cPfHgMcyxn0q9PjaQrQjIiLjp2/aiohEhAJfRCQiFPgiIhGhwBcRiQgFvohIRBTkUzrFLJF0Nu9uY0drJ5fWzGLtxfOJx2yyyxIROeumdOAnks577tvCtuYOunsTVJTFWVVbxYO3r1Hoi0jkTOkunc2729jW3EFXbwIHunoTbGvuYPNufdlXRKJnSl/h72jtpLs3MWhcV2+C773Qqi4eEYmcKR34l9bMoqIsTlco9GMG/7n9EL39SXXxiEikTOkunbUXz2dVbRWW6AVPUl6SWt2e/qS6eEQkcqZ04MdjxoO3r6G68XtUtfyMGy5biPvgabp7E+xs7ZycAkVEzqIpHfiQCv3pHXupOvAMv72yhoqy+KDnK8rirKiZNUnViYicPVO6Dz9Tuovn6T0H8VgJ08tLWVVbxdqL509ou6P5LoC+LyAiEy1SgZ/u4rn67bfTO2M+n//kH094sI7muwD6voCInA2RCnw408UzvWMv65YvmPD2wt8FgDM3ip/cdZhYzNjR2kl/Ipl1ms27285KjcWmEO929I5JZKjIBf7Zluu7AH/x6E6OnOqluzdBSdzoSwy+m5y+mTwVA3+4MM72bmflknP4H9ecx66DJ7i0ZhZvrqvmqcb2nGE+mmXoBDCxdMItTpEP/PG8MMcyT7bvApSXxDh8ooee/iTAkLAHmFYaoy+R5O5NjVPqgMnVffXA+1cTjxlP/urwkHc7W14+xtZXOujtTzKtNEZZSYy+hGedHxhxGel57g/NM9ra/2vPmf3+losmfp9ktnnNhdX8tKl9UA1AweryzI+xjbPm933jWV5oCZ9wq7j//W8Ytq7Uurazs7WTFTWzeMtF1Xlv38xlXnPhPH7adGTYNsYzT6GVxCfm8zQFCXwzux74EhAHvubun814vhz4JvB64CjwbnffV4i2h9PTn+DIyd6BYG053hWMTw3vO3KKjz/yAjsPvsrpvlSYrFh0Dp9/18ohV4zP7D1KY9tJLqiewXeea2HXwc6BeZYvnMU765fwUvsp6uZXctX5cwfmv6C6kuULZ7H15cN4rISKshJmTy/j4Kunhxac7AeLMa00Tkksxlef2puzjTcsm8Oz+47R2HZySJvZ6s41TVqu4zyRdLa8fJTGwyepW1DJmvNyL2M4jtPbn+SpxiNs3X+c08E+6OpN8Nz+49z7k71cuXQ2P9xxeMg7oqSf2WfdfUm6+5IDz2XOD4y4jPQ8Xw3Nk0w625o72Hf0FMvmzmBVbRXAwLhz50znB9sP0dR+kt7+JGUlMS6cX8knblhOLI+DP7Pdyxefw4sHXs3aZmncKI3H6E/6QA0XVM/AsAmva1Vt1bDLC0+fSDpbXzk+aHtvfeU4X3vq5YHtnW3+v/7PXTS1FW49MpeZbftltjGeeSbCmvPmTMjy8w58M4sD9wDXAS3As2a20d13hia7HTju7hea2a3A54B359v2SHr6k7xytIuevtTB33ysOzU+GP7+iwfZ3vrqoDDZ3voq33/x4KAgCL8A4jEj4T4QkN19SZ5v7uCXra/Sn/CBg/CGyxbxyrEuls2dwcfW1fGHf/wVEpUL+PAHN5B058s/bhpoF6AsbsRefhrzJDe843f4/i8PDqor3EauF+Edv3lJzrAYzwt1Ig7CXQc7B603QG9/kn1HT3Hl0tksmzuDspLYkGmGE54fGNUywvNkW8/MIM3c7z39SRoPn+Bft7YQj9moQjHTSOGS2WZvwulNnDmR9fQn2XP4JGZn3iX29CdpajvJtuaOnOE61rpyBWOuk2EsZiSSg68gMvdRpm3NHTS1nRzYZ+Ndj8wTT3iZ2bZf5j5M+tjnCZ+kx/M6OJss37dwZnY18Gl3/81g+E4Ad/+b0DSPB9M8bWYlwCGg2odpfM7S5X7dJ74+rpq2vbANgNe9biWnevtp3LkdgLoVlwEMDFfV1nHkZO+Q+asry5g3sxyAE6f7OdDRnfMKOBcDHDCDitI4PYf3YkEN7s4rx7rp6ukDDIvZoGly1TWS8pIYvYlkzlrNYHFVBTOnje48n2vd584owwymlcapLI9jdubF7e6c7Elwui+R9flsywzXlbltsJEPHDOoOWcaZsbpvgTlJTGOd/WNuIz0erjDsVO9ZG629D4cDTOYVhJj7owyTvcns647DN4+udothFnTSigrieWsYzij3UfdwTqMxkivvfYTPSMeiyMZT12ZNZbEht5PG0kseA2Fj/dz51QMu80zj5MZZTFO9SYHhtOv5/F45INvfM7d67M9V4guncVAc2i4BViTaxp37zezV4G5wJHwRGa2AdgAULnognEXtGrlKoCBq4x00Kelh0+c7h844MOOtLdxtB1qahYNHJhjlZ7FHbr7Eiw+/5KBF7uZce6cCk72lNHTl6A8fVDOHb6ukYx0VewOBw61Y8l+amoWUVkep2nXDmDoybBuxWU51/3oyZ7UesSMaSUxeo4fgngZixYtHBS0mSeyuhWXUVkep6I0PnBQpoOytXk/xMuoqVlE7expnOpNbZuyILzDB3HmATYtI+DTdS2ZPZ2e4Ao12zKOnjzNcCeEsWz+1H5OBt2GNmTb1NQsYkZZjObjp8d0MhuNbCemzu5ecu0DIOdFUK79Hn7tzJ43n+7gL9COqr4gBGeUxdjd+NKg7dH40ssQL2P2vPlZ1yN8LIanT79+zWyg9oXnXzLusE+vY1/Cx3Sih1R3YXgZ3X0J9jTtJdbXPbB93Z09TXuzHieYETNIJpOk91lndx+XLJw57tDPpahu2rr7vcC9APX19f7tP7g6r+V1nu5jx4HcfzYh11vqU/3T8XgJ7Sd6WDCrnPKM7gEDSuJGf8KHvO3OvmJw9flzefuVS0ZV90jdSNnEs7yNzmSAl5Tjsem0n+hh1rQZlDf/gmTlfG583W9w+eJz+MPH7yNRuYAbX/cbWbueUgtKvQjd4XRfEq+YDbE4hzt7SLiDxQaeTySdiuYtlB1t4lN/tnFg/TK7A3Z0V0Gwzc+ZX8nn3n5Z1u6DbG+h03WG2006rF+9NGsffSLpqS6zYTZXScwGdZWkt196v2frtkhNdKaG8LZpP9FDbFZ5ah4b/Q25cJuj6cMfeK0w/D740D/fCZB1eOv+49z9ZOOQ13z6tdPR3TdiIJbFjRsvrxnU7fHZx39FclpqP7d1nqY0HhsYPnG6j+ll8YF1yzwWM6dvP9HDORWpbsz0a3bpnCs5cLx7SC3xmJFMDt1+ufZh9cxyXu3uy9rNlnO/h7hD6cnDxLuODRxX4XXPPE4gOGmEXjv9SecP3nLBuD6l98gHcz9XiMA/ANSGhpcE47JN0xJ06ZxD6ubtpIrFjE/csHxIEFBSBqSumA+9epqF51RwuPP0oIMs3Uef2X+ZLZzLSmIsmztj3HWN5ubdglnlHO7sGfHERLx0YN12HTyBX3oLxOJ8adOe1AG24haIl3D3k41cUD2DC6srh+2bdRhYZn+WA6G3P0mscgEcbRq0flcunc2VS2ezdf9xmtpPDtrmmX234enTwsP/trWF3mHuC2QuI9v0cCYYct0MDe/3gZPGMO+qwtumpz/JgY7TWcMiVyBlux+Urb8Yztxcbu3o5ucvDT60su2D4ayqreLC+ZVDLjjS65JtHcKvtXSf/zuuXDJw0t66/zhNbWf280DfeGgYkgMnicxjMXP6dF/6Xd/dzsngNbvl5WND3u5kO/Gkt1+2fVheEuM9Vy0lZpb1AmM0+92A3upLIBbn7icbB47NdO3ZjpNME/Wx7EIE/rNAnZmdRyrYbwV+J2OajcBtwNPAO4Enh+u/L5TppXGWL5o59O1pxnSXLJoJwFef2jskCPoSzg2XLeCSRbPYc+gEFy2cyZsumDfokyrvfkMtP3vpCHsOneDC+ZU8/Itmtrd2DvTHXbZ4Fv+9fknWT7cMtxEuDuoCeNcbavn5S0fYfegEFy+cyZrz5rLl5aODhj/68PNsP3Cm3UtrZrJ+9VIa206w/2gXP9h+aGjbwUGc7YB6+UgXf3nLpcRixu5DJ+hPOt965hW6+wZ/AmY400rjVJxup6Iszrlzp9PTl6CnP0lvf5KkOy0dXVnD+kBHF2+8cO6o2qhbUEl5aYzToU/vlJfGqFtQybTSoVfTWacvMW5auZjSuLFs3gyuPDd1otj6ynH2HTk1MC78fYHGtpM0tp2gpy9JSTwVUsMdy4mkUxKzQQd8Zrsrl1TxQktH1jbT3njh3CHbJj3u2X3HeO6V44PWrazEIBbn9LJr+OWBDlYuqaJ/Xh39lQv45YEOrjh39kAwlwfb69M3Xcrzrxzn5SOnOHC8m6caB/W+Aql3QYmkU14a48LqSm5cWcP+o6c4b94MrsiouznLfs7Ul0gt652vr+VfGppHnL434bR2dA8KUgNKgxNPeUmMugUzuXX1uYNqueqCuVx1wdwz+/DwCXr6kwPTrwk+zXbVBXNHnKcsbpSETtIl6Qu+0Im+teP0qEI+bKL+xlfegR/0yX8YeJzUxzK/7u47zOwzQIO7bwTuAx40sybgGKmTwoQriceoml426unfsGwO33x6/6DPzFeUxVl93twRz7T/7Yoz3TW3XLGEzbvbBj63W6jP0IfbyDb88Iarc7a7addh/mtP+6B1G8npvgStHaf5yLo6IBVYO1s7Bz5Dn/o8fHJQyMUMSuOxQZ93P/DUPiweY3FVxZA2jp3q5bvbWods83XLF3DFuaP7dMblS6rYvLt9yGf7b7/m/KzbPdf0n33H5UOmr182J2e7//GhNw1s70sWzuQbP9837LapKE2909t/rGvYdtecP7oT3Wi2Rfp7C521q/FYCV94Yg9lJTFOrLgZj5XwxR81sqq2ihnlpRg+cKKD1PEAqddOw/7jg/dRaYwPvPl8SuOxUb3Gj5/qZWPGfs5UURbn2uULeP3S2XR09bLxheGnTwd7mAO/9bpFXFBdOaq6wvtwtMdq5jzpLwLubO3kpfaTfHdb66Dp+4N3b+HuwfBxkuu7JRPxN77y/pTORKmvr/eGhoaz2uZU/ps2meuWLZAyTS+L8/frrxh0skt/6SxbyIW/0fqrgycGDqB1b/11ADZv3jxiXePd5uG6RnPgjnX6sdaQbdukv+yVDodCXgzkqqMvkeRrP315+KAtjVG29ykM5/Of/Niovrk81n2UuYxcIZde5mimP3fOdF451jVo3bK9Zs+mTbsO85GHnh9ycsw80WceJ+GTRr6vCzPL+SkdBX6GiQiCYjFcII10AI5mmbm219q1a4HsgT/aZbwWFcN63b2pkS88sWfkT50kE2A28BdkM/d7IdYlcxkjhdxI07+5rpr3feMXRXWBluvkeDZO9GkKfMlqrAfgeI0U+DJxsl1xjmSyr5LHohhOqsVWkwJfJpUCf/KM1DWS7Q/3GfAn1100cO9GXluGC/yi+hy+iBRW+jcgcr2Ty9bHr1+Bm7oU+CJTXDxmrFu+YFAXTXo4kXQa9h8f0uc80b8CJ5NDgS8SYdneARRDP7hMDAW+SMRlewcgU9PE/JV9EREpOgp8EZGIUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhEKPBFRCJCgS8iEhEKfBGRiFDgi4hEhAJfRCQi8gp8M5tjZk+YWWPw/5BfnTazVWb2tJntMLMXzezd+bQpIiLjk+8V/h3AJnevAzYFw5m6gPe6+6XA9cAXzawqz3ZFRGSM8g38m4EHgscPALdkTuDue9y9MXjcCrQB1Xm2KyIiY5Rv4C9w94PB40PAsH9Q28xWA2XASzme32BmDWbW0N7enmdpIiISNuIPoJjZj4CFWZ66Kzzg7m5mOX8R3cwWAQ8Ct7l7Mts07n4vcC+kfsR8pNpERGT0Rgx8d78213NmdtjMFrn7wSDQ23JMNwv4PnCXuz8z7mpFRGTc8u3S2QjcFjy+Dfhu5gRmVgb8O/BNd/9Onu2JiMg45Rv4nwWuM7NG4NpgGDOrN7OvBdO8C/g14H1mti34tyrPdkVEZIzy+hFzdz8KrMsyvgH4/eDxPwH/lE87IiKSP33TVkQkIhT4IiIRocAXEYkIBb6ISEQo8EVEIkKBLyISEQp8EZGIUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhEKPBFRCJCgS8iEhEKfBGRiFDgi4hEhAJfRCQiFPgiIhGRV+Cb2Rwze8LMGoP/Zw8z7SwzazGzL+fTpoiIjE++V/h3AJvcvQ7YFAzn8hfAT/JsT0RExinfwL8ZeCB4/ABwS7aJzOz1wALgh3m2JyIi45Rv4C9w94PB40OkQn0QM4sBnwf+dKSFmdkGM2sws4b29vY8SxMRkbCSkSYwsx8BC7M8dVd4wN3dzDzLdH8EPObuLWY2bFvufi9wL0B9fX22ZYmIyDiNGPjufm2u58zssJktcveDZrYIaMsy2dXAm83sj4BKoMzMTrr7cP39IiJSYPl26WwEbgse3wZ8N3MCd/9ddz/X3ZeR6tb5psI+OhJJp6vqfDoWX82mXYdJJPXGTWSy5Bv4nwWuM7NG4NpgGDOrN7Ov5VucvLYlks577ttCe91v07HkjXzkoed5z31bFPoik2TELp3huPtRYF2W8Q3A72cZfz9wfz5tymvH5t1tbGvuwONlAHT1Jggh7msAAATrSURBVNjW3MHm3W2sWz7k/r6ITDB901YmzI7WTrp7E4PGdfcm2NnaOUkViUSbAl8mzKU1s6goiw8aV1EWZ0XNrEmqSCTaFPgyYdZePJ9VtVVML4tjwPSyOKtqq1h78fzJLk0kkvLqwxcZTjxmPHj7GjbvbmNnaycramax9uL5xGPDfx9DRCaGAl8mVDxmrFu+QDdpRYqAunRERCJCgS8iEhEKfBGRiFDgi4hEhAJfRCQizL04/66JmbUD+/NYxDzgSIHKmUiqs7BUZ2GpzsI6G3UudffqbE8UbeDny8wa3L1+susYieosLNVZWKqzsCa7TnXpiIhEhAJfRCQipnLg3zvZBYyS6iws1VlYqrOwJrXOKduHLyIig03lK3wREQlR4IuIRMSUC3wzu97MdptZk5kVzY+lm9nXzazNzLaHxs0xsyfMrDH4f/Zk1hjUVGtmPzaznWa2w8w+Woy1mtk0M/uFmb0Q1PnnwfjzzGxLsP+/bWZlk1lnmpnFzex5M3s0GC66Os1sn5n90sy2mVlDMK6o9ntQU5WZfcfMfmVmu8zs6mKr08wuDrZj+l+nmX1ssuucUoFvZnHgHuAGYAWw3sxWTG5VA+4Hrs8Ydwewyd3rgE3B8GTrBz7u7iuAq4APBduw2GrtAd7q7iuBVcD1ZnYV8DngC+5+IXAcuH0Sawz7KLArNFysdf66u68KfVa82PY7wJeAH7j7JcBKUtu1qOp0993BdlwFvB7oAv6dya7T3afMP+Bq4PHQ8J3AnZNdV6ieZcD20PBuYFHweBGwe7JrzFLzd4HrirlWYDqwFVhD6luMJdleD5NY3xJSB/dbgUcBK9I69wHzMsYV1X4HzgFeJvjASbHWmVHbbwA/K4Y6p9QVPrAYaA4NtwTjitUCdz8YPD4EFNWvhJjZMuAKYAtFWGvQTbINaAOeAF4COty9P5ikWPb/F4E/A5LB8FyKs04Hfmhmz5nZhmBcse3384B24BtBF9nXzGwGxVdn2K3AQ8HjSa1zqgX+a5anTvlF8xlZM6sE/hX4mLt3hp8rllrdPeGpt8xLgNXAJZNc0hBmdiPQ5u7PTXYto3CNu19Jqkv0Q2b2a+Eni2S/lwBXAv/o7lcAp8joFimSOgEI7s3cBPxL5nOTUedUC/wDQG1oeEkwrlgdNrNFAMH/bZNcDwBmVkoq7L/l7v8WjC7KWgHcvQP4MamukSozS/90ZzHs/zcBN5nZPuBhUt06X6L46sTdDwT/t5Hqb15N8e33FqDF3bcEw98hdQIotjrTbgC2uvvhYHhS65xqgf8sUBd8AqKM1FupjZNc03A2ArcFj28j1V8+qczMgPuAXe7+d6GniqpWM6s2s6rgcQWp+wy7SAX/O4PJJr1Od7/T3Ze4+zJSr8cn3f13KbI6zWyGmc1MPybV77ydItvv7n4IaDazi4NR64CdFFmdIes5050Dk13nZN/QmIAbJG8D9pDqz71rsusJ1fUQcBDoI3WVcjupvtxNQCPwI2BOEdR5Dam3mS8C24J/byu2WoHLgeeDOrcDnwrGnw/8Amgi9Ta6fLK3aajmtcCjxVhnUM8Lwb8d6WOn2PZ7UNMqoCHY9/8BzC7SOmcAR4FzQuMmtU79aQURkYiYal06IiKSgwJfRCQiFPgiIhGhwBcRiQgFvohIRCjwRUQiQoEvIhIR/x+wzCuTdaMumAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_acf(df['count_n'], lags=72)\n",
    "\n",
    "plot_acf(df['count_n'].diff(24)[24:],lags=72)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #df.set_index('datum_zeit', inplace=True, drop=True)\n",
    "# plt.figure(figsize=(40,20))\n",
    "# result = seasonal_decompose(df4['count_n'], model='add',freq=24)\n",
    "# result.plot()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "if df.shape[0] < 1: print(\"No data in dataframe. Exiting\"),exit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "fill_index = pd.date_range('2018-05-01 00:00:00', '2019-04-30 23:00:00', freq='1H')\n",
    "df = df.reindex(fill_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['link_id']=sel_link_id\n",
    "df['dir_travel']=df['dir_travel'].fillna(method='backfill')\n",
    "df['epoch_60']=df.index.hour\n",
    "df['weekday_n']=df.index.dayofweek+1\n",
    "\n",
    "df[\"hourweek\"] = (df[\"weekday_n\"]-1)*24+df[\"epoch_60\"]\n",
    "df[\"hourweek\"] = pd.to_numeric(df.hourweek, errors='coerce')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_mean = df['count_n'].groupby(df['hourweek']).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, row in df.iterrows():\n",
    "    m_idx = row['hourweek']\n",
    "    if isnan(row['count_n']):\n",
    "        row['count_n'] = count_mean[m_idx]\n",
    "        df['count_n'].at[index] = count_mean[m_idx].round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "remain_nan = df['count_n'].isna().sum()\n",
    "if remain_nan > 0:\n",
    "    print(\"%i Remaining cols without value. Aborting.\" % (remain_nan))\n",
    "    exit(1)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[:-456]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_24 = df['count_n'].diff(24)\n",
    "count_24 = count_24[24:]\n",
    "count_1 = df['count_n'].diff()\n",
    "count_1.fillna(0,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df5xddX3n8dd7ZpKQH4QQSEJ+kaBETFAJOA1QbZuCbEEt+Oi6FKwKPKipu+Jqa7fFH4stXa12H63WSrtSQZRVkLWrG5UuaoAVrSCDhB9JjAkxmARIhh9DCAkZZuazf5xzkzs3987MnXvn/jjn/Xw88si933PuPefMPed9vud7vuccRQRmZpZ9Hc2eATMzawwHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD36xBJG2X9KZxfvY3JG2u9zxZvjjwrWVIulvSc5KmVPGZkHTyRM5XM5QuV0TcExGnNHOerP058K0lSFoK/AYQwIVNnZlRSOoaS5lZq3HgW6t4N3AvcBNwWaEwrfX/YdH7yyX9KH39w7T4IUn7JP1+Wv4eSVslPStpraQFRZ8/VdL302G7JX0kLZ8i6bOSnkj/fbZwpCFptaSdkv5c0lPAl8qVpeO+VdJ6SX2S/k3S68otrKRVkn6SjvekpM9LmlxpuQrTK/r88vRv0ydpg6QLi4bdJOk6Sd+V9IKk+yS9cnw/i2WJA99axbuBr6b/fkfSvNE+EBG/mb48LSJmRMTXJZ0D/DVwMTAfeBy4FUDS0cAPgP8LLABOBtal3/FR4CxgJXAasAr4WNHkTgBmA0uANeXKJJ0O3Aj8EXAc8AVgbYUmqkHgj4HjgbOBc4H/VGm5ij8oaRLwbeB7wFzg/cBXJRU3+VwC/CVwLLAV+ETZP6LligPfmk7SG0lC87aIeAB4DHjHOL/uD4AbI+JnEXEQ+DBwdtpk9FbgqYj424h4KSJeiIj7ij53bUTsiYhekrB8V9H3DgEfj4iDEXGgQtka4AsRcV9EDEbEl4GDJDuSYSLigYi4NyIGImI7yc7ht8a4jGcBM4BPRUR/RNwJfAe4tGicb0bETyNigGQnunKM320Z5sC3VnAZ8L2IeDp9/zWKmnWqtICkVg9AROwDngEWAotJdiajfi59vaDofW9EvFTymdKyJcCH0maWPkl96TQXlHwOSa+S9B1JT0naC3ySpLY/FguAHRExVDK/C4veP1X0ej/JDsJyziearKkkTSVpfulM28IBpgCzJJ0GvAhMK/rICaN85RMkwVv4/ukkzSu7gB0kTR0jfW5D+v7EtKyg3G1lS8t2AJ+IiLE0n/wT8CBwaUS8IOmDwNvH8LnCvC6W1FEU+icCvxjj5y2nXMO3ZnsbSXv2CpJmh5XAcuAeknb99cDvSZqWdlO8suTzu4FXFL2/BbhC0sq07fyTwH1ps8l3gPmSPpiepD1a0plFn/uYpDmSjgeuAf5nlcvyz8B7JZ2pxHRJb0nPHZQ6GtgL7JP0auA/jrJcxe4jqbX/maRJklYDv0t6rsKsEge+NdtlwJci4lcR8VThH/B5knb1zwD9JAH4ZZL26GJ/AXw5bUK5OCJ+APxX4F+AJ4FXktbqI+IF4DyScHwK2AL8dvo9/w3oAR4GHgF+lpaNWUT0AO9J5/05kpOll1cY/U9JzlO8QLKj+HrJ8GHLVTKd/nQZLgCeBv4ReHdE/Lya+bX8kR+AYmaWD67hm5nlhAPfzCwnHPhmZjnhwDczy4mW7Yd//PHHx9KlS5s9G2ZmbeWBBx54OiLmlBvWsoG/dOlSenp6mj0bZmZtRdLjlYa5ScfMLCcc+GZmOeHANzPLCQe+mVlOOPDNzHKiLoEv6UZJeyQ9WmG4JH0ufezcw5LOqMd0yxkcCtZt2s3n1m1h3abdDA75XkFmZlC/bpk3kdwh8CsVhl8ALEv/nUlyL/AzK4w7boNDwbtuuI/1O/o40D/I1MmdrFw8i5uvPJPODtV7cmZmbaUuNfyI+CHw7AijXAR8JRL3kjzcYn49pl3s7s17WL+jj/39gwSwv3+Q9Tv6uHvznnpPysys7TSqDX8hydOACnYy/HFsAEhaI6lHUk9vb2/VE9nwxF4O9A8OKzvQP8jGJ/ZW/V1mZlnTUidtI+L6iOiOiO45c8peGTyiUxfMZOrkzmFlUyd3smLBzHrNoplZ22pU4O8ieZhzwaK0rK5WnzKXlYtnocF+iCGmpW34q0+ZW+9JmZm1nUYF/lrg3WlvnbOA5yPiyXpPpLND3HzlmczZ8m1m7fwx/3Dp6T5ha2aWqksvHUm3AKuB4yXtBD4OTAKIiP8B3A68meQZn/uBK+ox3XI6O8S0vm1M69vGucvnTdRkzMzaTl0CPyIuHWV4AO+rx7TMzGx8WuqkrZmZTRwHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczy4m6BL6k8yVtlrRV0tVlhp8o6S5JD0p6WNKb6zFdMzMbu5oDX1IncB1wAbACuFTSipLRPgbcFhGnA5cA/1jrdM3MrDr1qOGvArZGxLaI6AduBS4qGSeAmenrY4An6jBdMzOrQj0CfyGwo+j9zrSs2F8A75S0E7gdeH+5L5K0RlKPpJ7e3t46zJqZmRU06qTtpcBNEbEIeDNws6Qjph0R10dEd0R0z5kzp0GzZmaWD/UI/F3A4qL3i9KyYlcCtwFExE+Ao4Dj6zBtMzMbo3oE/v3AMkknSZpMclJ2bck4vwLOBZC0nCTw3WZjZtZANQd+RAwAVwF3AJtIeuNskHStpAvT0T4EvEfSQ8AtwOUREbVO28zMxq6rHl8SEbeTnIwtLrum6PVG4A31mJaZmY2Pr7Q1M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8uJugS+pPMlbZa0VdLVFca5WNJGSRskfa0e0zUzs7HrqvULJHUC1wHnATuB+yWtjYiNReMsAz4MvCEinpM0t9bpmplZdepRw18FbI2IbRHRD9wKXFQyznuA6yLiOYCI2FOH6ZqZWRXqEfgLgR1F73emZcVeBbxK0o8l3Svp/HJfJGmNpB5JPb29vXWYNTMzK2jUSdsuYBmwGrgU+GdJs0pHiojrI6I7IrrnzJnToFkzM8uHegT+LmBx0ftFaVmxncDaiHg5In4J/IJkB2BmZg1Sj8C/H1gm6SRJk4FLgLUl43yLpHaPpONJmni21WHaZmY2RjUHfkQMAFcBdwCbgNsiYoOkayVdmI52B/CMpI3AXcB/iYhnap22mZmNXc3dMgEi4nbg9pKya4peB/An6T8zM2sCX2lrZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTnQ1ewbM2sHgUHD35j1seGIvpy6YyepT5tLZoWbPlllVHPhmoxgcCt51w32s39HHgf5Bpk7uZOXiWdx85ZkOfWsrdWnSkXS+pM2Stkq6eoTx/r2kkNRdj+lWY3AoWLdpN59bt4V1m3YzOBSNngVrU3dv3sP6HX3s7x8kgP39g6zf0cfdm/c0e9bMqlJzDV9SJ3AdcB6wE7hf0tqI2Fgy3tHAB4D7ap1mtVxDs1pseGIvB/oHh5Ud6B9k4xN7OXf5vCbNlVn16lHDXwVsjYhtEdEP3ApcVGa8vwI+DbxUh2lWxTU0q8WpC2YydXLnsLKpkztZsWBmk+bIbHzqEfgLgR1F73emZYdIOgNYHBHfHemLJK2R1COpp7e3tw6zlhiphmY2mtWnzGXl4llosB9iiGnpEeLqU+Y2e9bMqjLh3TIldQB/B3xotHEj4vqI6I6I7jlz5tRtHlxDs1p0doibrzyTOVu+zaydP+YfLj3dzYHWluoR+LuAxUXvF6VlBUcDrwHulrQdOAtY28gTt66hWa06O8S0vm3M2nUv5y6f57C3tlSPwL8fWCbpJEmTgUuAtYWBEfF8RBwfEUsjYilwL3BhRPTUYdpj4hqamVkdeulExICkq4A7gE7gxojYIOlaoCci1o78DY1RqKFN69vmnhVWV74oy2rRyPWnLhdeRcTtwO0lZddUGHd1PaZp1grc5ddq0ej1x/fSMauBu/xaLRq9/jjwzYpUe0W2u/xaLRq9/vheOmap8RxeF7r87i/aaN3l18aq0euPa/hmqfEcXrvLr9Wi0euPA98sNZ7Da3f5tVo0ev1xk45ZaryH1+7ya7Vo5PrjGr5Zys0zlnUOfLOUm2cs69ykY1bEzTOWZQ58M7MGaIVbcDjwzcwmWKvcgiP3bfh+1q2ZTbRWuQVHrmv4rbLXNbNsa5XnIue6ht8qe10zy7ZWeepergPfN74ys0ZolWs8ch34rbLXNbNsa5VrPHId+K2y1zWz7GuF5yLnOvBbZa9rZtYIue6lA76y0iZOK1xoY1Ys94Fv+TTRYewuv9aKHPiWO40I4+IuvzC8y6+PJK1Zct2Gb/nUiOsv3OXXWpED33KnEWHsLr/Wihz4ljuNCGN3+c2nVr83l9vwLXcKYfyTXzxJdHQxbcqkuodxocvv2b93Jf3T5/K3H/tj99LJuHY4UV+XGr6k8yVtlrRV0tVlhv+JpI2SHpa0TtKSekzXbDwadf1FK1xoY43TDvfmqjnwJXUC1wEXACuASyWtKBntQaA7Il4HfAP4m1qna1YLh7HVWzucqK9HDX8VsDUitkVEP3ArcFHxCBFxV0TsT9/eCyyqw3TNzFpGO5yor0fgLwR2FL3fmZZVciXwr+UGSFojqUdST29vbx1mzcysMdrhRH1De+lIeifQDfz3csMj4vqI6I6I7jlz5jRy1szMatIO9+aqRy+dXcDioveL0rJhJL0J+CjwWxFxsA7TNTNrKa1+b656BP79wDJJJ5EE/SXAO4pHkHQ68AXg/IhonVPWZk3gm6pZs9Qc+BExIOkq4A6gE7gxIjZIuhboiYi1JE04M4D/JQngVxFxYa3TnkjeKG0itENfbcuuulx4FRG3A7eXlF1T9PpN9ZhOo3ijtInim6pZM/nWCmW0wwUU1p7aoa+2ZZcDvwxvlDZR2qGvtmWXA78Mb5TZ0Io3smqHvtqWXb55WhmNuLmWTaxWPQ/jm6plQ7t26nDgl+GNsv218snRVu+rbSNr1crEWLhJpwLfXKu9+TyMTZR27tThwLdM8nkYmyjtXJlw4Fsm+eSoTZR2rkw48C2T2uFGVtae2rky4ZO2llk+OWoToZ07dTjwzVpEu3b1y6N2rUw48M1aQDt39bP24TZ8sxbQzl39sqwVr9auhWv44+BD79aShd9jpK5+7dRkkCVZPOpy4FcpiytBO8vK71Ho6re/KPTbpatfVrXy1drj5SadKvnQu7Vk5fdo565+WdXOF1hV4sCvUhZXgnaWld/D1w20nna+wKoSB36VsrgStLMs/R6+f1NryeJRlwO/SllcCdpZXn6PrPUWaQdZPOrySdsqtfNVdlmUh98jKyem21G7XmBViQN/HCqtBFnoHtiOsrZRlspib5GISP9P35cbduh90bCiMWOEg5xyw4LxHRUNpV+2v39g2PcWjrJePDhQNK9xqHzvSy8n40cy7QgYGAyC4Ol9B4m0nID+wSEIeKLvAAEsOOYopPpnhwO/TlwLy5fBoWBgaCj9Pw5tyHv2vsRgutEPDcFLLychvXXPPkg3+qHgUHgXTi4XwujFgwMAPLrr+UPTuvPne8qemL5r8x7mzTwq+fw4sqw0AEu/o1wIF5cPH798gEeMHODtYN9LyW/y0I7nh5UXfquHd5Yv37DryI4D+/uTYVt27xtW/lL6+z7+zH4A5s88ignIewd+vUx0LWzYBjesxlNhnGHlhbLhG2W54UdusMmnDn2mwneNVEMrNw9HjHDk26qCoX9gCICnnn+pbPmTzx84Yl4PpsN2PLufoSiEcXAgDenNT72QBHca4PsOJjW8e7c9c8S8FTbkx3pfLDv93hcODisfGEzKnz/w8rDyQu3whTRkABYcM5XJXR2H5hdgclcH82dOHTaejd/QULB+Rx/bn3mRpcdNZ+XiWXRksKKW+cCPCA4OHK6JEcnGV6ihFWphQbKBFzbuCNiX7qkfePw5KAnDwob2018+S0Rwx4anytbC7tjwFDOmjP5nLpdt7VYTmmgjbZSFmvQvnx4euIXy7U/vP+L7DqbDdj53YFj5y2mwPvti/xHTh8b/LisXz+LkuTPY8KunobOLKZO6OHnuDFYuntXYGcmooaHgk/+6ia179tE/MMTkrg5OnjuDj1ywPHOhX5fAl3Q+8PdAJ/DFiPhUyfApwFeA1wPPAL8fEdvrMe3RvNg/yCPpIdf+NMCTw+vD+kfZwPuLalYFhVpsoUa2ZPb0srWwJbOn4w4VtcvTRlmqo0N85ILl/NEHPsTgjHlc9d41ma2BNsP6HX1s3bPv0LZ7cGCIrXv2sX5HH2csObbJc1dfNXfLlNQJXAdcAKwALpW0omS0K4HnIuJk4DPAp2udbqsp1MIYSLoHTkkDybWw+ijeKIPhG2UedHSIyc9sZerjP+aMJcc67Oto+zMvHlGp6x8YYvszL1b4RPuqRz/8VcDWiNgWEf3ArcBFJeNcBHw5ff0N4FxNxCnoJirUwmZs/BZTf3kP//mcZbmofTZKnjZKa6ylxyVH58Umd3Ww9LjpTZqjiaOosUFS0tuB8yPiD9P37wLOjIirisZ5NB1nZ/r+sXScpyt97+wly+O8j9w4rnla/9B6AFaetpLBiENnzbdsfBSAZSteM2z8asvH+xkbvxdeGmBX34Fh7ecSLJw1laOP6mrIb9iI7xrPNKw2EcGvnj3A/oMvA0IdYuqkTk6cPRVJTVlPZh41adzLc9t7f/2BiOguN6ylAl/SGmANwIz5r3z9mz9+c03zBgwL/GZq1bBol3ArbJQHXh4kIgn74o0yz1rxt63ndzViGhHBvoODHHx5kCmTOpkxpbOp61UrB/7ZwF9ExO+k7z8MEBF/XTTOHek4P5HUBTwFzIkRJt7d3R09PT01zRskPW0eKekn2wzve8eFAFz3tbVjHlav8lb9rmqnkZeuc9Vqxd+2nt/VqPltJWeeNHvc67akioFfj1469wPLJJ0E7AIuAd5RMs5a4DLgJ8DbgTtHCnuzcjo6xBlLjs1czwmzRqk58CNiQNJVwB0k3TJvjIgNkq4FeiJiLXADcLOkrcCzJDuFhpg2qZPXLJx5qB/+YNG/0vdDcfgim6H0IpyhoZILj8ZhaCjoP+5kBmfM42ePP+eaqdWN162x89+qTv3wI+J24PaSsmuKXr8E/Id6TKtaHR3i6Braw4oV3/+j0hWupQaHgsu/dD8vnvo2oqOL6+7eymmLZnHTFb926JYLhfnrXjq85loof/2SY4ddoTrjqORnO/3E4V0+C+UrF8864p4j09OLv1676JhhV8JOn9JFQHI74Th8Jey0ycn4p5xw9KGrbQvfVbgd8SvnTh/2hzgqLT/p+MO9G4LgqElJ+ZLjph2aLsCUtHzx7KnD5rV42KJjDw8br3L3YimUDRVdSVxYzmRHz6Gd/2C64y9UBFpF4dqEfSveBp1dfO7OLbm5NqFa/lslMn+lbT0VTuIceS6n8grz/36xm4d29hGdk4HklgsP7ezjR1ufPnTLhcL3Teoc3jWsUF7aZawjHVAI0tLy0vvDA4d2LqVX/RbKj5k6fKfY1ZmUz54++YjvKszn3KOPGlY+OS0/4ZiS8nT+F8waHt5T0vJFx047YhqFYYtnHzmsmSIOHxkODAWDg4fvqVN8D53C64jDR4uFm2UlO5fyt6eodHuMcnq2P8tjvfugK/mNDg4M8VjvPjY8+TzdSw+3AZeuP4X1OCmPssMmdWrYkEJ5Z4eGV3DaJCsL13EU/62yenHVSBz4E8wPp67e4FCwf9Yr6J8+j3WbdrfUXUclMalTTDpyn9pwP9ryNAdfHn5twsGXhzj48hBnnHjsoZ3765eUHjmWLy8e1r10dtnyVScNLy/0Jjn7lccBw4+CZ06dRAC/lh65xqHvSj5zxpJZh8aF4UeoyfhxqDwCXrfomCOOrAtHrqcuHP7Am2lpeeFBOPds6S17Hcf+/kFWzK/uYTnl7rpZ7v5WxfehOvSZtHXgcAUgHS8Of89QxITcOA0c+BPOD6euTuGuo73Lfpfo6OL9tzzou45W0IrrVulRsICuCkeuU7rGdoTaIYEON0sWK6wTpd0Yu0qOXF+/5Niyf6szlszimGn1afJtB37i1QQrPJFp2uROBJl9IlO9FO46Gp2TQR1t+1DyRvC6NXb+WyVcw59ghScy3b15Dxuf2MsKPxhlRG4CG7uR1q1WbhZrBm+HCQd+A3R2iHOXzysbWN4wh2vFZopWVm7dcrNYeSNth3nhJp0mKt4w+xb9Ou+/5UHedcN9uX5AtQ+9a5f3ZrFCJapv4dl+4HsJ1/CbaNiGSTaeVTpWlY5sfOhduzw3i/noZmSu4TfRSBtmlo12ZFM49H7/ucs4d/k8b6hVKjSLFctLs1jej25G48BvovFsmFk4XPVGObFGahbLwvozkrxWosbKTTpNVNgw1+/o40D/IFPLbJjFTR5AJg5X89zk0AiVmsUgG+vPSHzSf2QO/CaqdsO84g1LK7b5rz5lbkv29im34/JGOfHK9UhZt2l35s8ZjVSJMgd+01WzYX734SfL1owf3fU8N/zolw2puVU62VrNEclNV6zyRtkEWTuyqrQu+qR/ZQ78FlRpwwTK1owH0weDTHTNrVIPiJuuWMXlX/rpmI9I7tnS642yCUY7sqpmZ97s32q03jh5729fiU/atqBKJ3Pf8tr5ZU/GdXRoxBNV9TpRV+lk6+fv3FK2vNIRycYn9ronThOMdjK3XM+p/oGhhlwrUu066hP/4+Mafguq1A55zvJ5nLN83hE147s376lYc6tnv+RKRx49jz9X1RGJ2+qbY6TmjkrNiMN25kzM0eNI6yhQ9ugia81TjeLAb0GjtUOWHq6OdKKqnhd3VWoS6F5yLOt39B1R/pbXzmfPCwfdVt9CKjV3VLszr2ewVlpH79y0my/92/ayOwKf+B8fB36LqqYdcqQdxEg1oUo9eyq12VbasVx1zjJ6Hn9uzEckbr5pPdXuzOsZrJXW0e8+8uSIvdJ84r96DvyMqLSDqLQhv/qEo6s6AVtoAqq0Y6nmiMRaT7U780rXioxnZ15pHQVGPLrwif/qOfAzrtKGjBhXm22lHYt7RrS3anfmML6LuMrtJCqto2953Xy+t3F3xaMLr3PVc+BnXKUN9rq7tjatzdZaUzU789Eu4qr2SvFKOxU329SXAz8Hym2wzWyztfY32rmhaq8UL6yfpTsbN9vUl/vh51SlPtlXnbPM96O3UY10479KfeRHui6jEl+vUV+u4efUeE/AmsHIXYErNReCr8toNgd+jvkErI3XSBWDSs2Fvi6j+WoKfEmzga8DS4HtwMUR8VzJOCuBfwJmAoPAJyLi67VM18yar1LFoNorxX302DiKGP89MST9DfBsRHxK0tXAsRHx5yXjvAqIiNgiaQHwALA8IvpG+u7u7u7o6ekZ97yZWfMMDoWDvUkkPRAR3eWG1dqkcxGwOn39ZeBuYFjgR8Qvil4/IWkPMAcYMfDNrH25WbA11dpLZ15EPJm+fgoY8deVtAqYDDxWYfgaST2Senp7e2ucNTMzKzZqDV/SD4ATygz6aPGbiAhJFduHJM0HbgYui4ihcuNExPXA9ZA06Yw2b2ZmNnajBn5EvKnSMEm7Jc2PiCfTQC97M2pJM4HvAh+NiHvHPbdmZjZutTbprAUuS19fBvyf0hEkTQa+CXwlIr5R4/TMzGycag38TwHnSdoCvCl9j6RuSV9Mx7kY+E3gcknr038ra5yumZlVqaZumRNJUi/weA1fcTzwdJ1mp53kdbnBy+5lz5dKy70kIuaU+0DLBn6tJPVU6ouaZXldbvCye9nzZTzL7ZunmZnlhAPfzCwnshz41zd7Bpokr8sNXva8yuuyV73cmW3DNzOz4bJcwzczsyIOfDOznMhc4Es6X9JmSVvTWzZnlqQbJe2R9GhR2WxJ35e0Jf3/2GbO40SRtFjSXZI2Stog6QNpeaaXX9JRkn4q6aF0uf8yLT9J0n3pev/19Ar3TJLUKelBSd9J3+di2SVtl/RIevFqT1pW1fqeqcCX1AlcB1wArAAulbSiuXM1oW4Czi8puxpYFxHLgHXp+ywaAD4UESuAs4D3pb911pf/IHBORJwGrATOl3QW8GngMxFxMvAccGUT53GifQDYVPQ+T8v+2xGxsqj/fVXre6YCH1gFbI2IbRHRD9xKcs/+TIqIHwLPlhRfRPJsAtL/39bQmWqQiHgyIn6Wvn6BJAAWkvHlj8S+9O2k9F8A5wCFe1VlbrkLJC0C3gJ8MX0vcrLsFVS1vmct8BcCO4re70zL8qSqZxRkgaSlwOnAfeRg+dMmjfUkd6f9PsnzJfoiYiAdJcvr/WeBPwMKt1g/jvwsewDfk/SApDVpWVXrux9inmGjPaMgCyTNAP4F+GBE7E0qfImsLn9EDAIrJc0iuRPtq5s8Sw0h6a3Anoh4QNLqZs9PE7wxInZJmgt8X9LPiweOZX3PWg1/F7C46P2itCxPdqfPJig8dKbsMwqyQNIkkrD/akT877Q4N8ufPhf6LuBsYJakQgUuq+v9G4ALJW0naa49B/h78rHsRMSu9P89JDv6VVS5vmct8O8HlqVn7ScDl5Dcsz9PRn1GQRakbbc3AJsi4u+KBmV6+SXNSWv2SJoKnEdy/rH69TIAAADZSURBVOIu4O3paJlbboCI+HBELIqIpSTb9p0R8QfkYNklTZd0dOE18O+AR6lyfc/clbaS3kzSztcJ3BgRn2jyLE0YSbeQPET+eGA38HHgW8BtwIkkt5e+OCJKT+y2PUlvBO4BHuFwe+5HSNrxM7v8kl5HcnKuk6TCdltEXCvpFSS13tnAg8A7I+Jg8+Z0YqVNOn8aEW/Nw7Kny/jN9G0X8LWI+ISk46hifc9c4JuZWXlZa9IxM7MKHPhmZjnhwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5z4/4X8g7bewyAxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df5xddX3n8dd7ZpKQH4QQSEJ+kaBETFAJOA1QbZuCbEEt+Oi6FKwKPKipu+Jqa7fFH4stXa12H63WSrtSQZRVkLWrG5UuaoAVrSCDhB9JjAkxmARIhh9DCAkZZuazf5xzkzs3987MnXvn/jjn/Xw88si933PuPefMPed9vud7vuccRQRmZpZ9Hc2eATMzawwHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD36xBJG2X9KZxfvY3JG2u9zxZvjjwrWVIulvSc5KmVPGZkHTyRM5XM5QuV0TcExGnNHOerP058K0lSFoK/AYQwIVNnZlRSOoaS5lZq3HgW6t4N3AvcBNwWaEwrfX/YdH7yyX9KH39w7T4IUn7JP1+Wv4eSVslPStpraQFRZ8/VdL302G7JX0kLZ8i6bOSnkj/fbZwpCFptaSdkv5c0lPAl8qVpeO+VdJ6SX2S/k3S68otrKRVkn6SjvekpM9LmlxpuQrTK/r88vRv0ydpg6QLi4bdJOk6Sd+V9IKk+yS9cnw/i2WJA99axbuBr6b/fkfSvNE+EBG/mb48LSJmRMTXJZ0D/DVwMTAfeBy4FUDS0cAPgP8LLABOBtal3/FR4CxgJXAasAr4WNHkTgBmA0uANeXKJJ0O3Aj8EXAc8AVgbYUmqkHgj4HjgbOBc4H/VGm5ij8oaRLwbeB7wFzg/cBXJRU3+VwC/CVwLLAV+ETZP6LligPfmk7SG0lC87aIeAB4DHjHOL/uD4AbI+JnEXEQ+DBwdtpk9FbgqYj424h4KSJeiIj7ij53bUTsiYhekrB8V9H3DgEfj4iDEXGgQtka4AsRcV9EDEbEl4GDJDuSYSLigYi4NyIGImI7yc7ht8a4jGcBM4BPRUR/RNwJfAe4tGicb0bETyNigGQnunKM320Z5sC3VnAZ8L2IeDp9/zWKmnWqtICkVg9AROwDngEWAotJdiajfi59vaDofW9EvFTymdKyJcCH0maWPkl96TQXlHwOSa+S9B1JT0naC3ySpLY/FguAHRExVDK/C4veP1X0ej/JDsJyziearKkkTSVpfulM28IBpgCzJJ0GvAhMK/rICaN85RMkwVv4/ukkzSu7gB0kTR0jfW5D+v7EtKyg3G1lS8t2AJ+IiLE0n/wT8CBwaUS8IOmDwNvH8LnCvC6W1FEU+icCvxjj5y2nXMO3ZnsbSXv2CpJmh5XAcuAeknb99cDvSZqWdlO8suTzu4FXFL2/BbhC0sq07fyTwH1ps8l3gPmSPpiepD1a0plFn/uYpDmSjgeuAf5nlcvyz8B7JZ2pxHRJb0nPHZQ6GtgL7JP0auA/jrJcxe4jqbX/maRJklYDv0t6rsKsEge+NdtlwJci4lcR8VThH/B5knb1zwD9JAH4ZZL26GJ/AXw5bUK5OCJ+APxX4F+AJ4FXktbqI+IF4DyScHwK2AL8dvo9/w3oAR4GHgF+lpaNWUT0AO9J5/05kpOll1cY/U9JzlO8QLKj+HrJ8GHLVTKd/nQZLgCeBv4ReHdE/Lya+bX8kR+AYmaWD67hm5nlhAPfzCwnHPhmZjnhwDczy4mW7Yd//PHHx9KlS5s9G2ZmbeWBBx54OiLmlBvWsoG/dOlSenp6mj0bZmZtRdLjlYa5ScfMLCcc+GZmOeHANzPLCQe+mVlOOPDNzHKiLoEv6UZJeyQ9WmG4JH0ufezcw5LOqMd0yxkcCtZt2s3n1m1h3abdDA75XkFmZlC/bpk3kdwh8CsVhl8ALEv/nUlyL/AzK4w7boNDwbtuuI/1O/o40D/I1MmdrFw8i5uvPJPODtV7cmZmbaUuNfyI+CHw7AijXAR8JRL3kjzcYn49pl3s7s17WL+jj/39gwSwv3+Q9Tv6uHvznnpPysys7TSqDX8hydOACnYy/HFsAEhaI6lHUk9vb2/VE9nwxF4O9A8OKzvQP8jGJ/ZW/V1mZlnTUidtI+L6iOiOiO45c8peGTyiUxfMZOrkzmFlUyd3smLBzHrNoplZ22pU4O8ieZhzwaK0rK5WnzKXlYtnocF+iCGmpW34q0+ZW+9JmZm1nUYF/lrg3WlvnbOA5yPiyXpPpLND3HzlmczZ8m1m7fwx/3Dp6T5ha2aWqksvHUm3AKuB4yXtBD4OTAKIiP8B3A68meQZn/uBK+ox3XI6O8S0vm1M69vGucvnTdRkzMzaTl0CPyIuHWV4AO+rx7TMzGx8WuqkrZmZTRwHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczy4m6BL6k8yVtlrRV0tVlhp8o6S5JD0p6WNKb6zFdMzMbu5oDX1IncB1wAbACuFTSipLRPgbcFhGnA5cA/1jrdM3MrDr1qOGvArZGxLaI6AduBS4qGSeAmenrY4An6jBdMzOrQj0CfyGwo+j9zrSs2F8A75S0E7gdeH+5L5K0RlKPpJ7e3t46zJqZmRU06qTtpcBNEbEIeDNws6Qjph0R10dEd0R0z5kzp0GzZmaWD/UI/F3A4qL3i9KyYlcCtwFExE+Ao4Dj6zBtMzMbo3oE/v3AMkknSZpMclJ2bck4vwLOBZC0nCTw3WZjZtZANQd+RAwAVwF3AJtIeuNskHStpAvT0T4EvEfSQ8AtwOUREbVO28zMxq6rHl8SEbeTnIwtLrum6PVG4A31mJaZmY2Pr7Q1M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8uJugS+pPMlbZa0VdLVFca5WNJGSRskfa0e0zUzs7HrqvULJHUC1wHnATuB+yWtjYiNReMsAz4MvCEinpM0t9bpmplZdepRw18FbI2IbRHRD9wKXFQyznuA6yLiOYCI2FOH6ZqZWRXqEfgLgR1F73emZcVeBbxK0o8l3Svp/HJfJGmNpB5JPb29vXWYNTMzK2jUSdsuYBmwGrgU+GdJs0pHiojrI6I7IrrnzJnToFkzM8uHegT+LmBx0ftFaVmxncDaiHg5In4J/IJkB2BmZg1Sj8C/H1gm6SRJk4FLgLUl43yLpHaPpONJmni21WHaZmY2RjUHfkQMAFcBdwCbgNsiYoOkayVdmI52B/CMpI3AXcB/iYhnap22mZmNXc3dMgEi4nbg9pKya4peB/An6T8zM2sCX2lrZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTnQ1ewbM2sHgUHD35j1seGIvpy6YyepT5tLZoWbPlllVHPhmoxgcCt51w32s39HHgf5Bpk7uZOXiWdx85ZkOfWsrdWnSkXS+pM2Stkq6eoTx/r2kkNRdj+lWY3AoWLdpN59bt4V1m3YzOBSNngVrU3dv3sP6HX3s7x8kgP39g6zf0cfdm/c0e9bMqlJzDV9SJ3AdcB6wE7hf0tqI2Fgy3tHAB4D7ap1mtVxDs1pseGIvB/oHh5Ud6B9k4xN7OXf5vCbNlVn16lHDXwVsjYhtEdEP3ApcVGa8vwI+DbxUh2lWxTU0q8WpC2YydXLnsLKpkztZsWBmk+bIbHzqEfgLgR1F73emZYdIOgNYHBHfHemLJK2R1COpp7e3tw6zlhiphmY2mtWnzGXl4llosB9iiGnpEeLqU+Y2e9bMqjLh3TIldQB/B3xotHEj4vqI6I6I7jlz5tRtHlxDs1p0doibrzyTOVu+zaydP+YfLj3dzYHWluoR+LuAxUXvF6VlBUcDrwHulrQdOAtY28gTt66hWa06O8S0vm3M2nUv5y6f57C3tlSPwL8fWCbpJEmTgUuAtYWBEfF8RBwfEUsjYilwL3BhRPTUYdpj4hqamVkdeulExICkq4A7gE7gxojYIOlaoCci1o78DY1RqKFN69vmnhVWV74oy2rRyPWnLhdeRcTtwO0lZddUGHd1PaZp1grc5ddq0ej1x/fSMauBu/xaLRq9/jjwzYpUe0W2u/xaLRq9/vheOmap8RxeF7r87i/aaN3l18aq0euPa/hmqfEcXrvLr9Wi0euPA98sNZ7Da3f5tVo0ev1xk45ZaryH1+7ya7Vo5PrjGr5Zys0zlnUOfLOUm2cs69ykY1bEzTOWZQ58M7MGaIVbcDjwzcwmWKvcgiP3bfh+1q2ZTbRWuQVHrmv4rbLXNbNsa5XnIue6ht8qe10zy7ZWeepergPfN74ys0ZolWs8ch34rbLXNbNsa5VrPHId+K2y1zWz7GuF5yLnOvBbZa9rZtYIue6lA76y0iZOK1xoY1Ys94Fv+TTRYewuv9aKHPiWO40I4+IuvzC8y6+PJK1Zct2Gb/nUiOsv3OXXWpED33KnEWHsLr/Wihz4ljuNCGN3+c2nVr83l9vwLXcKYfyTXzxJdHQxbcqkuodxocvv2b93Jf3T5/K3H/tj99LJuHY4UV+XGr6k8yVtlrRV0tVlhv+JpI2SHpa0TtKSekzXbDwadf1FK1xoY43TDvfmqjnwJXUC1wEXACuASyWtKBntQaA7Il4HfAP4m1qna1YLh7HVWzucqK9HDX8VsDUitkVEP3ArcFHxCBFxV0TsT9/eCyyqw3TNzFpGO5yor0fgLwR2FL3fmZZVciXwr+UGSFojqUdST29vbx1mzcysMdrhRH1De+lIeifQDfz3csMj4vqI6I6I7jlz5jRy1szMatIO9+aqRy+dXcDioveL0rJhJL0J+CjwWxFxsA7TNTNrKa1+b656BP79wDJJJ5EE/SXAO4pHkHQ68AXg/IhonVPWZk3gm6pZs9Qc+BExIOkq4A6gE7gxIjZIuhboiYi1JE04M4D/JQngVxFxYa3TnkjeKG0itENfbcuuulx4FRG3A7eXlF1T9PpN9ZhOo3ijtInim6pZM/nWCmW0wwUU1p7aoa+2ZZcDvwxvlDZR2qGvtmWXA78Mb5TZ0Io3smqHvtqWXb55WhmNuLmWTaxWPQ/jm6plQ7t26nDgl+GNsv218snRVu+rbSNr1crEWLhJpwLfXKu9+TyMTZR27tThwLdM8nkYmyjtXJlw4Fsm+eSoTZR2rkw48C2T2uFGVtae2rky4ZO2llk+OWoToZ07dTjwzVpEu3b1y6N2rUw48M1aQDt39bP24TZ8sxbQzl39sqwVr9auhWv44+BD79aShd9jpK5+7dRkkCVZPOpy4FcpiytBO8vK71Ho6re/KPTbpatfVrXy1drj5SadKvnQu7Vk5fdo565+WdXOF1hV4sCvUhZXgnaWld/D1w20nna+wKoSB36VsrgStLMs/R6+f1NryeJRlwO/SllcCdpZXn6PrPUWaQdZPOrySdsqtfNVdlmUh98jKyem21G7XmBViQN/HCqtBFnoHtiOsrZRlspib5GISP9P35cbduh90bCiMWOEg5xyw4LxHRUNpV+2v39g2PcWjrJePDhQNK9xqHzvSy8n40cy7QgYGAyC4Ol9B4m0nID+wSEIeKLvAAEsOOYopPpnhwO/TlwLy5fBoWBgaCj9Pw5tyHv2vsRgutEPDcFLLychvXXPPkg3+qHgUHgXTi4XwujFgwMAPLrr+UPTuvPne8qemL5r8x7mzTwq+fw4sqw0AEu/o1wIF5cPH798gEeMHODtYN9LyW/y0I7nh5UXfquHd5Yv37DryI4D+/uTYVt27xtW/lL6+z7+zH4A5s88ignIewd+vUx0LWzYBjesxlNhnGHlhbLhG2W54UdusMmnDn2mwneNVEMrNw9HjHDk26qCoX9gCICnnn+pbPmTzx84Yl4PpsN2PLufoSiEcXAgDenNT72QBHca4PsOJjW8e7c9c8S8FTbkx3pfLDv93hcODisfGEzKnz/w8rDyQu3whTRkABYcM5XJXR2H5hdgclcH82dOHTaejd/QULB+Rx/bn3mRpcdNZ+XiWXRksKKW+cCPCA4OHK6JEcnGV6ihFWphQbKBFzbuCNiX7qkfePw5KAnDwob2018+S0Rwx4anytbC7tjwFDOmjP5nLpdt7VYTmmgjbZSFmvQvnx4euIXy7U/vP+L7DqbDdj53YFj5y2mwPvti/xHTh8b/LisXz+LkuTPY8KunobOLKZO6OHnuDFYuntXYGcmooaHgk/+6ia179tE/MMTkrg5OnjuDj1ywPHOhX5fAl3Q+8PdAJ/DFiPhUyfApwFeA1wPPAL8fEdvrMe3RvNg/yCPpIdf+NMCTw+vD+kfZwPuLalYFhVpsoUa2ZPb0srWwJbOn4w4VtcvTRlmqo0N85ILl/NEHPsTgjHlc9d41ma2BNsP6HX1s3bPv0LZ7cGCIrXv2sX5HH2csObbJc1dfNXfLlNQJXAdcAKwALpW0omS0K4HnIuJk4DPAp2udbqsp1MIYSLoHTkkDybWw+ijeKIPhG2UedHSIyc9sZerjP+aMJcc67Oto+zMvHlGp6x8YYvszL1b4RPuqRz/8VcDWiNgWEf3ArcBFJeNcBHw5ff0N4FxNxCnoJirUwmZs/BZTf3kP//mcZbmofTZKnjZKa6ylxyVH58Umd3Ww9LjpTZqjiaOosUFS0tuB8yPiD9P37wLOjIirisZ5NB1nZ/r+sXScpyt97+wly+O8j9w4rnla/9B6AFaetpLBiENnzbdsfBSAZSteM2z8asvH+xkbvxdeGmBX34Fh7ecSLJw1laOP6mrIb9iI7xrPNKw2EcGvnj3A/oMvA0IdYuqkTk6cPRVJTVlPZh41adzLc9t7f/2BiOguN6ylAl/SGmANwIz5r3z9mz9+c03zBgwL/GZq1bBol3ArbJQHXh4kIgn74o0yz1rxt63ndzViGhHBvoODHHx5kCmTOpkxpbOp61UrB/7ZwF9ExO+k7z8MEBF/XTTOHek4P5HUBTwFzIkRJt7d3R09PT01zRskPW0eKekn2wzve8eFAFz3tbVjHlav8lb9rmqnkZeuc9Vqxd+2nt/VqPltJWeeNHvc67akioFfj1469wPLJJ0E7AIuAd5RMs5a4DLgJ8DbgTtHCnuzcjo6xBlLjs1czwmzRqk58CNiQNJVwB0k3TJvjIgNkq4FeiJiLXADcLOkrcCzJDuFhpg2qZPXLJx5qB/+YNG/0vdDcfgim6H0IpyhoZILj8ZhaCjoP+5kBmfM42ePP+eaqdWN162x89+qTv3wI+J24PaSsmuKXr8E/Id6TKtaHR3i6Braw4oV3/+j0hWupQaHgsu/dD8vnvo2oqOL6+7eymmLZnHTFb926JYLhfnrXjq85loof/2SY4ddoTrjqORnO/3E4V0+C+UrF8864p4j09OLv1676JhhV8JOn9JFQHI74Th8Jey0ycn4p5xw9KGrbQvfVbgd8SvnTh/2hzgqLT/p+MO9G4LgqElJ+ZLjph2aLsCUtHzx7KnD5rV42KJjDw8br3L3YimUDRVdSVxYzmRHz6Gd/2C64y9UBFpF4dqEfSveBp1dfO7OLbm5NqFa/lslMn+lbT0VTuIceS6n8grz/36xm4d29hGdk4HklgsP7ezjR1ufPnTLhcL3Teoc3jWsUF7aZawjHVAI0tLy0vvDA4d2LqVX/RbKj5k6fKfY1ZmUz54++YjvKszn3KOPGlY+OS0/4ZiS8nT+F8waHt5T0vJFx047YhqFYYtnHzmsmSIOHxkODAWDg4fvqVN8D53C64jDR4uFm2UlO5fyt6eodHuMcnq2P8tjvfugK/mNDg4M8VjvPjY8+TzdSw+3AZeuP4X1OCmPssMmdWrYkEJ5Z4eGV3DaJCsL13EU/62yenHVSBz4E8wPp67e4FCwf9Yr6J8+j3WbdrfUXUclMalTTDpyn9pwP9ryNAdfHn5twsGXhzj48hBnnHjsoZ3765eUHjmWLy8e1r10dtnyVScNLy/0Jjn7lccBw4+CZ06dRAC/lh65xqHvSj5zxpJZh8aF4UeoyfhxqDwCXrfomCOOrAtHrqcuHP7Am2lpeeFBOPds6S17Hcf+/kFWzK/uYTnl7rpZ7v5WxfehOvSZtHXgcAUgHS8Of89QxITcOA0c+BPOD6euTuGuo73Lfpfo6OL9tzzou45W0IrrVulRsICuCkeuU7rGdoTaIYEON0sWK6wTpd0Yu0qOXF+/5Niyf6szlszimGn1afJtB37i1QQrPJFp2uROBJl9IlO9FO46Gp2TQR1t+1DyRvC6NXb+WyVcw59ghScy3b15Dxuf2MsKPxhlRG4CG7uR1q1WbhZrBm+HCQd+A3R2iHOXzysbWN4wh2vFZopWVm7dcrNYeSNth3nhJp0mKt4w+xb9Ou+/5UHedcN9uX5AtQ+9a5f3ZrFCJapv4dl+4HsJ1/CbaNiGSTaeVTpWlY5sfOhduzw3i/noZmSu4TfRSBtmlo12ZFM49H7/ucs4d/k8b6hVKjSLFctLs1jej25G48BvovFsmFk4XPVGObFGahbLwvozkrxWosbKTTpNVNgw1+/o40D/IFPLbJjFTR5AJg5X89zk0AiVmsUgG+vPSHzSf2QO/CaqdsO84g1LK7b5rz5lbkv29im34/JGOfHK9UhZt2l35s8ZjVSJMgd+01WzYX734SfL1owf3fU8N/zolw2puVU62VrNEclNV6zyRtkEWTuyqrQu+qR/ZQ78FlRpwwTK1owH0weDTHTNrVIPiJuuWMXlX/rpmI9I7tnS642yCUY7sqpmZ97s32q03jh5729fiU/atqBKJ3Pf8tr5ZU/GdXRoxBNV9TpRV+lk6+fv3FK2vNIRycYn9ronThOMdjK3XM+p/oGhhlwrUu066hP/4+Mafguq1A55zvJ5nLN83hE147s376lYc6tnv+RKRx49jz9X1RGJ2+qbY6TmjkrNiMN25kzM0eNI6yhQ9ugia81TjeLAb0GjtUOWHq6OdKKqnhd3VWoS6F5yLOt39B1R/pbXzmfPCwfdVt9CKjV3VLszr2ewVlpH79y0my/92/ayOwKf+B8fB36LqqYdcqQdxEg1oUo9eyq12VbasVx1zjJ6Hn9uzEckbr5pPdXuzOsZrJXW0e8+8uSIvdJ84r96DvyMqLSDqLQhv/qEo6s6AVtoAqq0Y6nmiMRaT7U780rXioxnZ15pHQVGPLrwif/qOfAzrtKGjBhXm22lHYt7RrS3anfmML6LuMrtJCqto2953Xy+t3F3xaMLr3PVc+BnXKUN9rq7tjatzdZaUzU789Eu4qr2SvFKOxU329SXAz8Hym2wzWyztfY32rmhaq8UL6yfpTsbN9vUl/vh51SlPtlXnbPM96O3UY10479KfeRHui6jEl+vUV+u4efUeE/AmsHIXYErNReCr8toNgd+jvkErI3XSBWDSs2Fvi6j+WoKfEmzga8DS4HtwMUR8VzJOCuBfwJmAoPAJyLi67VM18yar1LFoNorxX302DiKGP89MST9DfBsRHxK0tXAsRHx5yXjvAqIiNgiaQHwALA8IvpG+u7u7u7o6ekZ97yZWfMMDoWDvUkkPRAR3eWG1dqkcxGwOn39ZeBuYFjgR8Qvil4/IWkPMAcYMfDNrH25WbA11dpLZ15EPJm+fgoY8deVtAqYDDxWYfgaST2Senp7e2ucNTMzKzZqDV/SD4ATygz6aPGbiAhJFduHJM0HbgYui4ihcuNExPXA9ZA06Yw2b2ZmNnajBn5EvKnSMEm7Jc2PiCfTQC97M2pJM4HvAh+NiHvHPbdmZjZutTbprAUuS19fBvyf0hEkTQa+CXwlIr5R4/TMzGycag38TwHnSdoCvCl9j6RuSV9Mx7kY+E3gcknr038ra5yumZlVqaZumRNJUi/weA1fcTzwdJ1mp53kdbnBy+5lz5dKy70kIuaU+0DLBn6tJPVU6ouaZXldbvCye9nzZTzL7ZunmZnlhAPfzCwnshz41zd7Bpokr8sNXva8yuuyV73cmW3DNzOz4bJcwzczsyIOfDOznMhc4Es6X9JmSVvTWzZnlqQbJe2R9GhR2WxJ35e0Jf3/2GbO40SRtFjSXZI2Stog6QNpeaaXX9JRkn4q6aF0uf8yLT9J0n3pev/19Ar3TJLUKelBSd9J3+di2SVtl/RIevFqT1pW1fqeqcCX1AlcB1wArAAulbSiuXM1oW4Czi8puxpYFxHLgHXp+ywaAD4UESuAs4D3pb911pf/IHBORJwGrATOl3QW8GngMxFxMvAccGUT53GifQDYVPQ+T8v+2xGxsqj/fVXre6YCH1gFbI2IbRHRD9xKcs/+TIqIHwLPlhRfRPJsAtL/39bQmWqQiHgyIn6Wvn6BJAAWkvHlj8S+9O2k9F8A5wCFe1VlbrkLJC0C3gJ8MX0vcrLsFVS1vmct8BcCO4re70zL8qSqZxRkgaSlwOnAfeRg+dMmjfUkd6f9PsnzJfoiYiAdJcvr/WeBPwMKt1g/jvwsewDfk/SApDVpWVXrux9inmGjPaMgCyTNAP4F+GBE7E0qfImsLn9EDAIrJc0iuRPtq5s8Sw0h6a3Anoh4QNLqZs9PE7wxInZJmgt8X9LPiweOZX3PWg1/F7C46P2itCxPdqfPJig8dKbsMwqyQNIkkrD/akT877Q4N8ufPhf6LuBsYJakQgUuq+v9G4ALJW0naa49B/h78rHsRMSu9P89JDv6VVS5vmct8O8HlqVn7ScDl5Dcsz9PRn1GQRakbbc3AJsi4u+KBmV6+SXNSWv2SJoKnEdy/rH69TIAAADZSURBVOIu4O3paJlbboCI+HBELIqIpSTb9p0R8QfkYNklTZd0dOE18O+AR6lyfc/clbaS3kzSztcJ3BgRn2jyLE0YSbeQPET+eGA38HHgW8BtwIkkt5e+OCJKT+y2PUlvBO4BHuFwe+5HSNrxM7v8kl5HcnKuk6TCdltEXCvpFSS13tnAg8A7I+Jg8+Z0YqVNOn8aEW/Nw7Kny/jN9G0X8LWI+ISk46hifc9c4JuZWXlZa9IxM7MKHPhmZjnhwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5z4/4X8g7bewyAxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_acf(df['count_n'].fillna(0), lags=48)\n",
    "#plot_pacf(df['count_n'], lags=48)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# p = 24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one-step sarima forecast\n",
    "def sarima_forecast(history, config):\n",
    "\torder, sorder, trend = config\n",
    "\t# define model\n",
    "\tmodel = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)\n",
    "\t# fit model\n",
    "\tmodel_fit = model.fit(disp=False)\n",
    "\t# make one step forecast\n",
    "\tyhat = model_fit.predict(len(history), len(history))\n",
    "\treturn yhat[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_test_split(data, n_test):\n",
    "\treturn data[:-n_test], data[-n_test:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2018-05-01 00:00:00     3.0\n",
      "2018-05-01 01:00:00     4.0\n",
      "2018-05-01 02:00:00     8.0\n",
      "2018-05-01 03:00:00     2.0\n",
      "2018-05-01 04:00:00     6.0\n",
      "                       ... \n",
      "2018-05-20 19:00:00    17.0\n",
      "2018-05-20 20:00:00    10.0\n",
      "2018-05-20 21:00:00     8.0\n",
      "2018-05-20 22:00:00    16.0\n",
      "2018-05-20 23:00:00     9.0\n",
      "Freq: H, Name: count_n, Length: 480, dtype: float64\n",
      " > Model[[(1, 0, 0), (0, 0, 0, 24), 'n']] 11.828\n",
      " > Model[[(1, 0, 0), (0, 0, 0, 0), 'n']] 11.828\n",
      " > Model[[(1, 0, 0), (0, 0, 1, 0), 'n']] 10.067\n",
      " > Model[[(1, 0, 0), (0, 0, 2, 0), 'n']] 10.055\n",
      " > Model[[(1, 0, 0), (0, 0, 1, 24), 'n']] 11.916\n",
      " > Model[[(1, 0, 0), (0, 1, 0, 24), 'n']] 13.974\n",
      " > Model[[(1, 0, 0), (0, 1, 1, 24), 'n']] 11.543\n",
      " > Model[[(1, 0, 0), (0, 0, 2, 24), 'n']] 12.007\n",
      " > Model[[(1, 0, 0), (1, 0, 0, 0), 'n']] 12.687\n",
      " > Model[[(1, 0, 0), (1, 0, 0, 24), 'n']] 12.126\n",
      " > Model[[(1, 0, 0), (1, 0, 1, 0), 'n']] 10.240\n",
      " > Model[[(1, 0, 0), (1, 0, 1, 24), 'n']] 10.775\n",
      " > Model[[(1, 0, 0), (1, 0, 2, 0), 'n']] 9.769\n",
      " > Model[[(1, 0, 0), (1, 0, 2, 24), 'n']] 10.562\n",
      " > Model[[(1, 0, 0), (1, 1, 0, 24), 'n']] 13.152\n",
      " > Model[[(1, 0, 0), (0, 1, 2, 24), 'n']] 11.909\n",
      " > Model[[(1, 0, 0), (1, 1, 1, 24), 'n']] 11.545\n",
      " > Model[[(1, 0, 0), (2, 0, 0, 0), 'n']] 10.070\n",
      " > Model[[(1, 0, 0), (2, 0, 0, 24), 'n']] 11.818\n",
      " > Model[[(1, 0, 0), (2, 0, 1, 0), 'n']] 9.762\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-58-8b8d27298cc3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    109\u001b[0m         \u001b[0mcfg_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msarima_configs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseasonal\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m24\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    110\u001b[0m         \u001b[0;31m# grid search\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 111\u001b[0;31m         \u001b[0mscores\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgrid_search\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcfg_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_test\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    112\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'done'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    113\u001b[0m         \u001b[0;31m# list top 3 configs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-58-8b8d27298cc3>\u001b[0m in \u001b[0;36mgrid_search\u001b[0;34m(data, cfg_list, n_test, parallel)\u001b[0m\n\u001b[1;32m     66\u001b[0m                 \u001b[0mexecutor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mParallel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_jobs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcpu_count\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'multiprocessing'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m                 \u001b[0mtasks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdelayed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscore_model\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_test\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcfg\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mcfg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcfg_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 68\u001b[0;31m                 \u001b[0mscores\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mexecutor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtasks\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     69\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     70\u001b[0m                 \u001b[0mscores\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mscore_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_test\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcfg\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mcfg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcfg_list\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/joblib/parallel.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, iterable)\u001b[0m\n\u001b[1;32m   1015\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1016\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mretrieval_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1017\u001b[0;31m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mretrieve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1018\u001b[0m             \u001b[0;31m# Make sure that we get a last message telling us we are done\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1019\u001b[0m             \u001b[0melapsed_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_start_time\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/joblib/parallel.py\u001b[0m in \u001b[0;36mretrieve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    907\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    908\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'supports_timeout'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 909\u001b[0;31m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    910\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    911\u001b[0m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.6/multiprocessing/pool.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    636\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    637\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 638\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    639\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mready\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    640\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mTimeoutError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.6/multiprocessing/pool.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    633\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    634\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 635\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_event\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    636\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    637\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.6/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    549\u001b[0m             \u001b[0msignaled\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flag\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    550\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0msignaled\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 551\u001b[0;31m                 \u001b[0msignaled\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cond\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    552\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0msignaled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    553\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.6/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    293\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m    \u001b[0;31m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    294\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 295\u001b[0;31m                 \u001b[0mwaiter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    296\u001b[0m                 \u001b[0mgotit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    297\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# one-step sarima forecast\n",
    "def sarima_forecast(history, config):\n",
    "\torder, sorder, trend = config\n",
    "\t# define model\n",
    "\tmodel = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)\n",
    "\t# fit model\n",
    "\tmodel_fit = model.fit(disp=False)\n",
    "\t# make one step forecast\n",
    "\tyhat = model_fit.predict(len(history), len(history))\n",
    "\treturn yhat[0]\n",
    " \n",
    "# root mean squared error or rmse\n",
    "def measure_rmse(actual, predicted):\n",
    "\treturn sqrt(mean_squared_error(actual, predicted))\n",
    " \n",
    "# split a univariate dataset into train/test sets\n",
    "def train_test_split(data, n_test):\n",
    "\treturn data[:-n_test], data[-n_test:]\n",
    " \n",
    "# walk-forward validation for univariate data\n",
    "def walk_forward_validation(data, n_test, cfg):\n",
    "\tpredictions = list()\n",
    "\t# split dataset\n",
    "\ttrain, test = train_test_split(data, n_test)\n",
    "\t# seed history with training dataset\n",
    "\thistory = [x for x in train]\n",
    "\t# step over each time-step in the test set\n",
    "\tfor i in range(len(test)):\n",
    "\t\t# fit model and make forecast for history\n",
    "\t\tyhat = sarima_forecast(history, cfg)\n",
    "\t\t# store forecast in list of predictions\n",
    "\t\tpredictions.append(yhat)\n",
    "\t\t# add actual observation to history for the next loop\n",
    "\t\thistory.append(test[i])\n",
    "\t# estimate prediction error\n",
    "\terror = measure_rmse(test, predictions)\n",
    "\treturn error\n",
    " \n",
    "# score a model, return None on failure\n",
    "def score_model(data, n_test, cfg, debug=False):\n",
    "\tresult = None\n",
    "\t# convert config to a key\n",
    "\tkey = str(cfg)\n",
    "\t# show all warnings and fail on exception if debugging\n",
    "\tif debug:\n",
    "\t\tresult = walk_forward_validation(data, n_test, cfg)\n",
    "\telse:\n",
    "\t\t# one failure during model validation suggests an unstable config\n",
    "\t\ttry:\n",
    "\t\t\t# never show warnings when grid searching, too noisy\n",
    "\t\t\twith catch_warnings():\n",
    "\t\t\t\tfilterwarnings(\"ignore\")\n",
    "\t\t\t\tresult = walk_forward_validation(data, n_test, cfg)\n",
    "\t\texcept:\n",
    "\t\t\terror = None\n",
    "\t# check for an interesting result\n",
    "\tif result is not None:\n",
    "\t\tprint(' > Model[%s] %.3f' % (key, result))\n",
    "\treturn (key, result)\n",
    " \n",
    "# grid search configs\n",
    "def grid_search(data, cfg_list, n_test, parallel=True):\n",
    "\tscores = None\n",
    "\tif parallel:\n",
    "\t\t# execute configs in parallel\n",
    "\t\texecutor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')\n",
    "\t\ttasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)\n",
    "\t\tscores = executor(tasks)\n",
    "\telse:\n",
    "\t\tscores = [score_model(data, n_test, cfg) for cfg in cfg_list]\n",
    "\t# remove empty results\n",
    "\tscores = [r for r in scores if r[1] != None]\n",
    "\t# sort configs by error, asc\n",
    "\tscores.sort(key=lambda tup: tup[1])\n",
    "\treturn scores\n",
    " \n",
    "# create a set of sarima configs to try\n",
    "def sarima_configs(seasonal=[0]):\n",
    "\tmodels = list()\n",
    "\t# define config lists\n",
    "\tp_params = [1, 2]\n",
    "\td_params = [0, 1]\n",
    "\tq_params = [0, 1, 2]\n",
    "\tt_params = ['n','c','t','ct']\n",
    "\tP_params = [0, 1, 2]\n",
    "\tD_params = [0, 1]\n",
    "\tQ_params = [0, 1, 2]\n",
    "\tm_params = seasonal\n",
    "\t# create config instances\n",
    "\tfor p in p_params:\n",
    "\t\tfor d in d_params:\n",
    "\t\t\tfor q in q_params:\n",
    "\t\t\t\tfor t in t_params:\n",
    "\t\t\t\t\tfor P in P_params:\n",
    "\t\t\t\t\t\tfor D in D_params:\n",
    "\t\t\t\t\t\t\tfor Q in Q_params:\n",
    "\t\t\t\t\t\t\t\tfor m in m_params:\n",
    "\t\t\t\t\t\t\t\t\tcfg = [(p,d,q), (P,D,Q,m), t]\n",
    "\t\t\t\t\t\t\t\t\tmodels.append(cfg)\n",
    "\treturn models\n",
    " \n",
    "if __name__ == '__main__':\n",
    "\t# define dataset\n",
    "\tdata = df['count_n'][:480].fillna(0)\n",
    "\tprint(data)\n",
    "\t# data split\n",
    "\tn_test = 48\n",
    "\t# model configs\n",
    "\tcfg_list = sarima_configs(seasonal=[0,24])\n",
    "\t# grid search\n",
    "\tscores = grid_search(data, cfg_list, n_test)\n",
    "\tprint('done')\n",
    "\t# list top 3 configs\n",
    "\tfor cfg, error in scores[:3]:\n",
    "\t\tprint(cfg, error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [(1, 0, 1), (0, 1, 0, 24), 'n'] 14.301852374209908\n",
    "# [(1, 0, 0), (0, 1, 0, 24), 'n'] 14.36446029101079\n",
    "# [(2, 0, 1), (0, 1, 0, 24), 'n'] 14.591357802813214"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SARIMA(p,d,q)(P,D,Q)m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = SARIMAX(df['count_n'], exog=None, order=(1, 0, 1), seasonal_order=(0, 1, 0, 24), trend=None, measurement_error=False, time_varying_regression=False, mle_regression=True, simple_differencing=False, enforce_stationarity=True, enforce_invertibility=True, hamilton_representation=False, concentrate_scale=False, trend_offset=1, use_exact_diffuse=False, dates=None, freq='1H')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = mod.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.summary().tables[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot_diagnostics()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = res.get_prediction(start=pd.to_datetime('2018-05-01'), dynamic=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "pred_ci = pred.conf_int()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for i in range(0,)\n",
    "ax = df['count_n'][:480].plot(label='observed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7, figsize=(14, 4))\n",
    "print(pred_ci)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.fill_between(pred_ci.index,\n",
    "                pred_ci.iloc[:, 0],\n",
    "                pred_ci.iloc[:, 1], color='k', alpha=.2)\n",
    "#ax.set_xlabel('Datum')\n",
    "#ax.set_ylabel('Fahrzeuge')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one-step sarima forecast\n",
    "def sarima_forecast(history, config):\n",
    "\torder, sorder, trend = config\n",
    "\t# define model\n",
    "\tmodel = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)\n",
    "\t# fit model\n",
    "\tmodel_fit = model.fit(disp=False)\n",
    "\t# make one step forecast\n",
    "\tyhat = model_fit.predict(len(history), len(history))\n",
    "\treturn yhat[0]\n",
    " \n",
    "# root mean squared error or rmse\n",
    "def measure_rmse(actual, predicted):\n",
    "\treturn sqrt(mean_squared_error(actual, predicted))\n",
    " \n",
    "# split a univariate dataset into train/test sets\n",
    "def train_test_split(data, n_test):\n",
    "\treturn data[:-n_test], data[-n_test:]\n",
    " \n",
    "# walk-forward validation for univariate data\n",
    "def walk_forward_validation(data, n_test, cfg):\n",
    "\tpredictions = list()\n",
    "\t# split dataset\n",
    "\ttrain, test = train_test_split(data, n_test)\n",
    "\t# seed history with training dataset\n",
    "\thistory = [x for x in train]\n",
    "\t# step over each time-step in the test set\n",
    "\tfor i in range(len(test)):\n",
    "\t\t# fit model and make forecast for history\n",
    "\t\tyhat = sarima_forecast(history, cfg)\n",
    "\t\t# store forecast in list of predictions\n",
    "\t\tpredictions.append(yhat)\n",
    "\t\t# add actual observation to history for the next loop\n",
    "\t\thistory.append(test[i])\n",
    "\t# estimate prediction error\n",
    "\terror = measure_rmse(test, predictions)\n",
    "\treturn error\n",
    " \n",
    "# score a model, return None on failure\n",
    "def score_model(data, n_test, cfg, debug=False):\n",
    "\tresult = None\n",
    "\t# convert config to a key\n",
    "\tkey = str(cfg)\n",
    "\t# show all warnings and fail on exception if debugging\n",
    "\tif debug:\n",
    "\t\tresult = walk_forward_validation(data, n_test, cfg)\n",
    "\telse:\n",
    "\t\t# one failure during model validation suggests an unstable config\n",
    "\t\ttry:\n",
    "\t\t\t# never show warnings when grid searching, too noisy\n",
    "\t\t\twith catch_warnings():\n",
    "\t\t\t\tfilterwarnings(\"ignore\")\n",
    "\t\t\t\tresult = walk_forward_validation(data, n_test, cfg)\n",
    "\t\texcept:\n",
    "\t\t\terror = None\n",
    "\t# check for an interesting result\n",
    "\tif result is not None:\n",
    "\t\tprint(' > Model[%s] %.3f' % (key, result))\n",
    "\treturn (key, result)\n",
    " \n",
    "# grid search configs\n",
    "def grid_search(data, cfg_list, n_test, parallel=True):\n",
    "\tscores = None\n",
    "\tif parallel:\n",
    "\t\t# execute configs in parallel\n",
    "\t\texecutor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')\n",
    "\t\ttasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)\n",
    "\t\tscores = executor(tasks)\n",
    "\telse:\n",
    "\t\tscores = [score_model(data, n_test, cfg) for cfg in cfg_list]\n",
    "\t# remove empty results\n",
    "\tscores = [r for r in scores if r[1] != None]\n",
    "\t# sort configs by error, asc\n",
    "\tscores.sort(key=lambda tup: tup[1])\n",
    "\treturn scores\n",
    " \n",
    "# create a set of sarima configs to try\n",
    "def sarima_configs(seasonal=[0]):\n",
    "\tmodels = list()\n",
    "\t# define config lists\n",
    "\tp_params = [1, 2]\n",
    "\td_params = [0, 1]\n",
    "\tq_params = [0, 1, 2]\n",
    "\tt_params = ['n','c','t','ct']\n",
    "\tP_params = [0, 1, 2]\n",
    "\tD_params = [0, 1]\n",
    "\tQ_params = [0, 1, 2]\n",
    "\tm_params = seasonal\n",
    "\t# create config instances\n",
    "\tfor p in p_params:\n",
    "\t\tfor d in d_params:\n",
    "\t\t\tfor q in q_params:\n",
    "\t\t\t\tfor t in t_params:\n",
    "\t\t\t\t\tfor P in P_params:\n",
    "\t\t\t\t\t\tfor D in D_params:\n",
    "\t\t\t\t\t\t\tfor Q in Q_params:\n",
    "\t\t\t\t\t\t\t\tfor m in m_params:\n",
    "\t\t\t\t\t\t\t\t\tcfg = [(p,d,q), (P,D,Q,m), t]\n",
    "\t\t\t\t\t\t\t\t\tmodels.append(cfg)\n",
    "\treturn models\n",
    " \n",
    "if __name__ == '__main__':\n",
    "\t# define dataset\n",
    "\tdata = df['count_n'].\n",
    "\tprint(data)\n",
    "\t# data split\n",
    "\tn_test = 48\n",
    "\t# model configs\n",
    "\tcfg_list = sarima_configs(seasonal=[0,24])\n",
    "\t# grid search\n",
    "\tscores = grid_search(data, cfg_list, n_test)\n",
    "\tprint('done')\n",
    "\t# list top 3 configs\n",
    "\tfor cfg, error in scores[:3]:\n",
    "\t\tprint(cfg, error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}