{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Anomaly Detector\n",
    "![resim.png](resim.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing the libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "plt.style.use('ggplot')\n",
    "import scipy.stats\n",
    "\n",
    "\n",
    "class simpleAnomalyDetector():\n",
    "    \"\"\"\n",
    "    USAGE\n",
    "    \n",
    "    calorieIntake = np.array([\n",
    "          [[100,20000],[300,400],[500,500]],\n",
    "          [[200,200],[200,200],[500,600]],\n",
    "          [[100,0],[100,0],[500,600]],\n",
    "          [[1000,0],[1000,0],[2000,0]],\n",
    "          [[100,100],[100,100],[500,500]],\n",
    "          [[200,200],[200,200],[500,500]],\n",
    "          [[400,300],[200,100],[500,500]]\n",
    "         ])\n",
    "    print(\"\\n\\nToy data: calorieIntake\")\n",
    "    print(calorieIntake)\n",
    "\n",
    "    # Learn Model Parameters\n",
    "    model = simpleAnomalyDetector()\n",
    "    model.fit(calorieIntake)\n",
    "\n",
    "    # Show Data\n",
    "    print(\"\\n\\n# Learn Model Parameters \\n# Show Data\")\n",
    "    print(model.df)\n",
    "\n",
    "    # Create a test case\n",
    "    print(\"\\n\\n# Create a test case\")\n",
    "    x_test = np.array([[1000,400],[200,200],[500,500]])\n",
    "    print(x_test)\n",
    "\n",
    "    # For a given time r and event c, show test value\n",
    "    r, c = 0,1\n",
    "    val = x_test[r, c]\n",
    "    print(\"\\n\\n# For a given time r = {} and event c = {}, show test value\".format(r,c))\n",
    "    print(val)\n",
    "\n",
    "    # Previously known data on that time r and event c\n",
    "    feature_size = x_test.shape[1]\n",
    "    idx = r*feature_size + c\n",
    "    print(\"\\n\\n# Previously known data on that time r = {} and event c= {}\".format(r,c))\n",
    "    print(\"# Look column {} of the dataframe\".format(idx))\n",
    "    print(model.df[idx])\n",
    "\n",
    "    ## Normal\n",
    "    x_test[r,c] = 500\n",
    "    print(\"\\n\\nPrediction for {}\".format(x_test[r,c]))\n",
    "    print(model.predict(test_matrix = x_test, row = r, col= c))\n",
    "\n",
    "\n",
    "    ## Normal\n",
    "    x_test[r,c] = 5000\n",
    "    print(\"\\n\\nPrediction for {}\".format(x_test[r,c]))\n",
    "    print(model.predict(test_matrix = x_test, row = r, col= c))\n",
    "\n",
    "\n",
    "    ## Anomaly\n",
    "    x_test[r,c] = 50000\n",
    "    print(\"\\n\\nPrediction for {}\".format(x_test[r,c]))\n",
    "    print(model.predict(test_matrix = x_test, row = r, col= c))\n",
    "    \"\"\"\n",
    "    def fit(self, data):\n",
    "        \"\"\"Read 3 dimensional array, first dimension is day, sample * features\"\"\"\n",
    "        data = data.flatten().reshape(data.shape[0], data.shape[1] * data.shape[2])\n",
    "        self.df = pd.DataFrame(data)\n",
    "        \n",
    "        # Learn parameters - central tendency and stdandard deviation\n",
    "        self.means = self.df.mean()\n",
    "        self.stds = self.df.std()\n",
    "    \n",
    "    def predict(self, test_matrix, row, col):\n",
    "        feature_size = test_matrix.shape[1]\n",
    "        idx = row*feature_size + col\n",
    "        \n",
    "        mu = self.means[idx]\n",
    "        sigma = self.stds[idx]\n",
    "        val = test_matrix[row,col]\n",
    "        \n",
    "        return np.abs(val - mu) > 3 * sigma, mu, sigma\n",
    "         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try on Toy Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Toy data: calorieIntake\n",
      "[[[  100 20000]\n",
      "  [  300   400]\n",
      "  [  500   500]]\n",
      "\n",
      " [[  200   200]\n",
      "  [  200   200]\n",
      "  [  500   600]]\n",
      "\n",
      " [[  100     0]\n",
      "  [  100     0]\n",
      "  [  500   600]]\n",
      "\n",
      " [[ 1000     0]\n",
      "  [ 1000     0]\n",
      "  [ 2000     0]]\n",
      "\n",
      " [[  100   100]\n",
      "  [  100   100]\n",
      "  [  500   500]]\n",
      "\n",
      " [[  200   200]\n",
      "  [  200   200]\n",
      "  [  500   500]]\n",
      "\n",
      " [[  400   300]\n",
      "  [  200   100]\n",
      "  [  500   500]]]\n",
      "\n",
      "\n",
      "# Learn Model Parameters \n",
      "# Show Data\n",
      "      0      1     2    3     4    5\n",
      "0   100  20000   300  400   500  500\n",
      "1   200    200   200  200   500  600\n",
      "2   100      0   100    0   500  600\n",
      "3  1000      0  1000    0  2000    0\n",
      "4   100    100   100  100   500  500\n",
      "5   200    200   200  200   500  500\n",
      "6   400    300   200  100   500  500\n",
      "\n",
      "\n",
      "# Create a test case\n",
      "[[1000  400]\n",
      " [ 200  200]\n",
      " [ 500  500]]\n",
      "\n",
      "\n",
      "# For a given time r = 0 and event c = 1, show test value\n",
      "400\n",
      "\n",
      "\n",
      "# Previously known data on that time r = 0 and event c= 1\n",
      "# Look column 1 of the dataframe\n",
      "0    20000\n",
      "1      200\n",
      "2        0\n",
      "3        0\n",
      "4      100\n",
      "5      200\n",
      "6      300\n",
      "Name: 1, dtype: int64\n",
      "\n",
      "\n",
      "Prediction for 500\n",
      "(False, 2971.4285714285716, 7509.7080026932)\n",
      "\n",
      "\n",
      "Prediction for 5000\n",
      "(False, 2971.4285714285716, 7509.7080026932)\n",
      "\n",
      "\n",
      "Prediction for 50000\n",
      "(True, 2971.4285714285716, 7509.7080026932)\n"
     ]
    }
   ],
   "source": [
    "calorieIntake = np.array([\n",
    "          [[100,20000],[300,400],[500,500]],\n",
    "          [[200,200],[200,200],[500,600]],\n",
    "          [[100,0],[100,0],[500,600]],\n",
    "          [[1000,0],[1000,0],[2000,0]],\n",
    "          [[100,100],[100,100],[500,500]],\n",
    "          [[200,200],[200,200],[500,500]],\n",
    "          [[400,300],[200,100],[500,500]]\n",
    "         ])\n",
    "print(\"\\n\\nToy data: calorieIntake\")\n",
    "print(calorieIntake)\n",
    "\n",
    "# Learn Model Parameters\n",
    "model = simpleAnomalyDetector()\n",
    "model.fit(calorieIntake)\n",
    "\n",
    "# Show Data\n",
    "print(\"\\n\\n# Learn Model Parameters \\n# Show Data\")\n",
    "print(model.df)\n",
    "\n",
    "# Create a test case\n",
    "print(\"\\n\\n# Create a test case\")\n",
    "x_test = np.array([[1000,400],[200,200],[500,500]])\n",
    "print(x_test)\n",
    "\n",
    "# For a given time r and event c, show test value\n",
    "r, c = 0,1\n",
    "val = x_test[r, c]\n",
    "print(\"\\n\\n# For a given time r = {} and event c = {}, show test value\".format(r,c))\n",
    "print(val)\n",
    "\n",
    "# Previously known data on that time r and event c\n",
    "feature_size = x_test.shape[1]\n",
    "idx = r*feature_size + c\n",
    "print(\"\\n\\n# Previously known data on that time r = {} and event c= {}\".format(r,c))\n",
    "print(\"# Look column {} of the dataframe\".format(idx))\n",
    "print(model.df[idx])\n",
    "\n",
    "## Normal\n",
    "x_test[r,c] = 500\n",
    "print(\"\\n\\nPrediction for {}\".format(x_test[r,c]))\n",
    "print(model.predict(test_matrix = x_test, row = r, col= c))\n",
    "\n",
    "\n",
    "## Normal\n",
    "x_test[r,c] = 5000\n",
    "print(\"\\n\\nPrediction for {}\".format(x_test[r,c]))\n",
    "print(model.predict(test_matrix = x_test, row = r, col= c))\n",
    "\n",
    "\n",
    "## Anomaly\n",
    "x_test[r,c] = 50000\n",
    "print(\"\\n\\nPrediction for {}\".format(x_test[r,c]))\n",
    "print(model.predict(test_matrix = x_test, row = r, col= c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1a19f8f588>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANMAAACPCAYAAABgS+5VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADRlJREFUeJzt3X9Mk3ceB/B3C27ommFbUKaTM4BmIYEwU28Z0QGD7RZlC3GeiYlb3Ey8jRkX2ZYg+8XmJLit8iPR210kxP1xyzwzm1zOxQRRSWxy64LgInGzt3mLBq9Au47uZPbHc38QHmH26dM+PKVf6vuVPIk+bXm+Jbzz/fbb5/P9GiRJkkBEs2ZMdQOI0gXDRKQTholIJwwTkU4YJiKdMExEOmGYiHTCMBHphGEi0klmsi+wNPuhZF9i3rj275OpboIwFuQUxHw8OPq9ptelUtLDRKRJOJjqFiSMYSIhSeFQqpuQMIaJxMQwEekk+GuqW5AwhomExGEekV44AUGkk9CtVLcgYQwTCUliz0SkE/ZMRDphz0SkE87mEelD4vdMRDphz0SkE4aJSCdBbbN5o6OjOHToEH766ScYDAbU1NRgw4YNOHbsGE6fPo37778fALB161asWbMGAHDixAn09vbCaDTihRdeQFlZGQBgYGAA3d3diEQiqK6uRl1dXcxrq4bp+vXrcLlc8Hq9MBgMMJvNsNlsePDBBzW9WaK4aOyZMjIy8Nxzz6GgoAA3b95EY2MjSktLAQAbN27EM888M+P5165dg9PpxMGDB+Hz+bBv3z50dHQAALq6uvDWW2/BarVi7969qn/3MSttHQ4H2tvbAQBFRUUoLCwEAHR0dMDhcGh6s0RxCYeiHyrMZjMKCiYLCBcuXIjly5fD6/UqPt/lcqG8vBwLFizAkiVLkJeXB7fbDbfbjby8PCxduhSZmZkoLy+Hy+WKee2YPdOZM2dgt9uRmTnzabW1tWhoaFDs9np6etDT0xPzwkQxKQzzfvu3VVNTg5qamqjP9Xg8+OGHH1BUVITLly/j1KlT6OvrQ0FBAZ5//nmYTCZ4vV6sWrVKfo3FYpHDZ7Va5fNWqxVXrlyJ2eSYYTIYDPD5fMjNzZ1x3ufzwWAwKL5u+hvs/jN7MNIgHI56OlZ4ppuYmIDdbsf27duxaNEiPPnkk9i8eTMA4PPPP8enn36K+vp6KO1bEe18rL95QCVM27dvx/vvv48HHnhATuno6Chu3LiBHTt2qL4hIs1C2mfzQqEQ7HY71q9fj0ceeQQAsHjxYvnx6upqHDhwAMBkjzM2NiY/5vV6YbFYAGDG+bGxMZjN5pjXjRmmsrIydHR0wO12y12fxWJBUVERjEYubERJFNJ2O5EkSfjkk0+wfPly1NbWyud9Pp8chq+++gorVqwAANhsNnR2dqK2thY+nw/Dw8MoKiqCJEkYHh6Gx+OBxWKB0+nE7t27Y15bdTbPaDRi9erVmt4YkWYae6Zvv/0WfX19yM/PxxtvvAFgchr8/PnzuHr1KgwGA3Jzc7Fz504AwIoVK/Doo4+ioaEBRqMRO3bskDuKF198Efv370ckEkFVVZUcQCWGZG92xqW+buNSX7epLdn1v7/uiXp+0c62ZDRHF/zSlsQU5F3jRPqYxQREqjBMJKYgw0SkD4XvmUTGMJGQpBDDRKQPTkAQ6YTDPCJ9cJhHpBfO5hHpQwpFUt2EhDFMJCaNwzylsvVAIIC2tjaMjIwgNzcXe/bsgclkgiRJ6O7uxoULF3Dvvfeivr5eLi48e/YsvvjiCwDApk2bUFlZGfPaDBMJSQpqC5NS2frZs2dRUlKCuro6OBwOOBwObNu2DRcuXMCNGzfQ2dmJK1eu4MiRI2hpaUEgEMDx48fR2toKAGhsbITNZoPJZFK8NsM0hxYuW5/qJggjdOu6yhO0DfPMZrNcajG9bN3lcqG5uRkAUFFRgebmZmzbtg1ff/01HnvsMRgMBqxevRq//PILfD4fLl26hNLSUjk8paWlGBgYwLp16xSvzTCRkJQ+M2ktW/f7/XLIzGYzfv75ZwCTxYA5OTnya6xWK7xeL7xe74yy9enl7EoYJhKSdCt6mLSWrSteJ4HydLWydZbLkpCkkBT1iEe0svXs7Gz4fD4Ak1W3U+vnWa1WjI6Oyq+dKk+3WCx3lLOrla0zTCQk6ZYU9VB9nULZus1mw7lz5wAA586dw9q1a+XzfX19kCQJ3333HRYtWgSz2YyysjIMDg4iEAggEAhgcHBQXpxSCStt59DYzfFUN0EYahMQYxsrop63/vNczNddvnwZ77zzDvLz8+Vh2datW7Fq1Sq0tbVhdHQUOTk5aGhokKfGu7q6MDg4iHvuuQf19fXy+pC9vb04ceIEgMmp8aqqqpjXZpjmEMN0m1qYRv8QPUw5p2KHKZU4AUFCCs+/jQMZJhKTFI49cyYihomEFAkxTES6CAfn30Qzw0RCinCYR6QPholIJxzmEekkHGGYiHTBYR6RTkKhjFQ3IWGaw3TmzBnVe5WItApHtPVMhw8fRn9/P7Kzs2G32wFgTnZaB2YRpmPHjimGiXva0mxFNIapsrISTz31FA4dOjTjfLJ3WgdUwvT6669HPS9JEvx+v+LruKctzVYwrG2YV1xcDI/HE9dzlXZaByDvtA5A3ml9VmHy+/148803cd999804L0kS3n777bgaTKRFWNJ3AiLZO60DKmFas2YNJiYmsHLlyjseKy4ujvd9ECUsJEWfGk9kDYgpc7HTOqASppdfflnxsVdffVX1hxNpFVLomeJdA2K6udhpHWDZOgkqDEPUQ4uptR+AO3dadzqdCAaD8Hg88k7rhYWF8k7roVAITqcTNptN9Tr8nomEFNQYnPb2dgwNDWF8fBwvvfQStmzZgkuXLiV9p3WAZetzimXrt6mVrf8jb2vU80/f+CwZzdEFeyYSktYhXSoxTCSkYByzZ6JhmEhI87BqnWEiMXGYR6ST4PzLEsNEYuIwj0gn87A2kGEiMQVT3QANGCYSEod5RDrRtqNtaiU9TL/PLkr2JeaN77NGUt2EeUPrbF60svW52Gkd4F3jJKgQpKiHmsrKSjQ1Nc0453A4UFJSgs7OTpSUlMDhmKz+nr7T+s6dO3HkyBEAkHdab2lpQUtLC44fP45AIKB6bYaJhBQ0RD/UFBcXyzukT3G5XKiomNzvqaKiAi6XCwAUd1ofGBiQd1o3mUzyTutq+JmJhBSOoxeK11zstA4wTCQopSGdlrJ1JXrutA4wTCSooEKYtIRnaqd1s9kc907rQ0ND8nmv1xvXmif8zERCCiscWszFTusAeyYSlNbPTNHK1uvq6tDW1obe3l55p3UAePjhh9Hf34/du3fLO60DgMlkwrPPPou9e/cCADZv3nzHpEY0SS9bfzq/Npk/fl75/ld+zzTl0n//FfPxP638Y9Tzf7n692Q0RxfsmUhIes7mzRWGiYTEMBHpJJjcTx9JwTCRkMKIpLoJCWOYSEjx3IcnGoaJhBSS2DMR6YI9E5FOQtL8Kw9UvZ3o+vXr+OabbzAxMTHjfDy3pBNpFYYU9RBZzDCdPHkSH374Ib788ku89tprch0IAHz2mbgLqNP8F5YiUQ+RxRzmnT59GgcOHEBWVhY8Hg8OHjyIkZERbNiwQXHXNYAbRNPsBefhMC9mmCKRCLKysgAAS5YsQXNzM+x2O0ZGRmKGafpt8k//jffmUeJmM6R75ZVXkJWVBaPRiIyMDLS2tmpaByJRMYd5ixcvxtWrV+X/Z2VlobGxEePj4/jxxx81XZAoHrMd5r377rv46KOP0NraCiDxdSC0iBmmXbt2zdgPFAAyMjKwa9cuvPfee5ovSqQmJIWjHlolug6EFjGHedPr4H/roYe4IyAlj1IvFG/Z+v79+wEATzzxBGpqahJeByKeDaF/i98zkZCU7s2Lp2x93759sFgs8Pv9+OCDD7Bs2TLF5yayDoQalq2TkIKRcNQjHhaLBcDk2g9r166F2+2W14EAENc6EFowTCQkrRMQExMTuHnzpvzvixcvIj8/P+F1ILTgMI+EpPULWr/fj48//njyZ4TDWLduHcrKylBYWJjQOhBacA2IOcQ1IG5TWwPid9bSqOf/M3YxGc3RBXsmEpLotw5FwzCRkMIRholIF/OxBINhIiGxZyLSSSjO75REwjCRkDgBQaQTDvOIdBKeh8O8pH9pK4qenh7Nm2KlG/4ukuOuuTePZfS38XeRHHdNmIiSjWEi0sldEyZ+RriNv4vkuGsmIIiS7a7pmYiSLe2/ZxoYGEB3dzcikQiqq6tRV1eX6ialzOHDh9Hf34/s7GzY7fZUNyftpHXPFIlE0NXVhaamJrS1teH8+fO4du1aqpuVMpWVlWhqakp1M9JWWofJ7XYjLy8PS5cuRWZmJsrLy2esl363KS4uhslkSnUz0lZah8nr9c5Y+29qTTSiZEjrMOm5JhqRmrQOk9VqxdjYmPz/2ayJRqQmrcNUWFiI4eFheDwehEIhOJ1O2Gy2VDeL0lTaf2nb39+Po0ePIhKJoKqqCps2bUp1k1Kmvb0dQ0NDGB8fR3Z2NrZs2YLHH3881c1KG2kfJqK5ktbDPKK5xDAR6YRhItIJw0SkE4aJSCcME5FOGCYinTBMRDr5PybTeQ3da0imAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a19e5fa20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "fig.set_size_inches(3,2)\n",
    "\n",
    "sns.heatmap(model.means.values.reshape(*calorieIntake[0].shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "# Previously known data on that time r = 0 and event c= 1\n",
      "# Look column 1 of the dataframe\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-27067.403439344227, 33010.26058220137)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAFDCAYAAABLHmiwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl0VPX9//HXnSQsk2FJJhBKQAtiKmAxkIAQkTW4hSpfZNGvS8FaFxBEqhXcUBEBKYLIZiVSRatWBbXLt+2PQxEVsSwGdyAiKmtIJkA2lsy9vz+QkZAQLnBnJpN5Ps6Zw8zn3jv3NW8g885dDcuyLAEAACAquMIdAAAAAKFD8wcAABBFaP4AAACiCM0fAABAFKH5AwAAiCI0fwAAAFGE5g8AACCK0PwBAABEEZo/AACAKELzBwAAEEViwx2gttu5c2fI1pWUlKSCgoKQra+uo57OoZbOop7Oop7OoZbOCnU9W7ZsaWs+tvwBAABEEZo/AACAKELzBwAAEEVo/gAAAKIIzR8AAEAUofkDAACIIjR/AAAAUYTmDwAAIIrUios8Hz58WJMmTVJFRYX8fr+6d++uYcOGKT8/X7Nnz1ZJSYnatGmjMWPGKDY2VkeOHNHcuXO1detWNWrUSOPGjVPz5s0lScuWLdOKFSvkcrk0cuRIpaWlSZJyc3O1ePFimaap/v37a9CgQeH8yAAAAGFRK7b8xcXFadKkSZoxY4aeeuop5ebmavPmzXr55ZeVnZ2tOXPmKD4+XitWrJAkrVixQvHx8Xr22WeVnZ2tV155RZK0fft2rV69Wk8//bQefPBB5eTkyDRNmaapnJwcPfDAA5o1a5Y+/PBDbd++PZwfGQAAICxqRfNnGIYaNGggSfL7/fL7/TIMQ1988YW6d+8uSerTp4/Wrl0rSVq3bp369OkjSerevbs+//xzWZaltWvXKjMzU3FxcWrevLlatGihvLw85eXlqUWLFkpOTlZsbKwyMzMD7wUAABBNasVuX0kyTVP333+/du/ercsvv1zJyclyu92KiYmRJCUmJsrn80mSfD6fvF6vJCkmJkZut1vFxcXy+Xw6//zzA+95/DLH5j/2fMuWLaH6aAAAALVGrWn+XC6XZsyYodLSUv3hD3/Qjh07TjqvZVlVxgzDqHa8pvmrs3z5ci1fvlySNG3aNCUlJdmJ74jY2NiQrq+uo57OoZbOop7Oop7OoZbOqq31rDXN3zHx8fHq0KGDtmzZorKyMvn9fsXExMjn8ykxMVHS0S13hYWF8nq98vv9Kisrk8fjCYwfc/wyx48XFhYqISGh2vVnZWUpKysr8LqgoCAYH7NaSUlJIV1fXUc9nUMtnUU9nUU9nUMtnRXqerZs2dLWfLXimL8DBw6otLRU0tEzfz/77DOlpKSoY8eOWrNmjSRp5cqVysjIkCSlp6dr5cqVkqQ1a9aoY8eOMgxDGRkZWr16tY4cOaL8/Hzt2rVL7dq103nnnaddu3YpPz9fFRUVWr16deC9AAAAokmt2PJXVFSkefPmyTRNWZalHj16KD09Xa1atdLs2bP12muvqU2bNurXr58kqV+/fpo7d67GjBkjj8ejcePGSZJat26tHj16aPz48XK5XPrNb34jl+tof3vLLbdoypQpMk1Tffv2VevWrcP2eQEAAMLFsE52oBwkSTt37gzZutjc7izq6Rxq6Szq6Szq6Rxq6Sx2+wIAACDsaP4AAACiCM0fAABAFKH5AwAAiCI0fwAAAFGE5g8AACCK0PwBAABEEZo/AACAKELzBwAAEEVo/gAAAKIIzR8AAEAUOePmzzRNvffee05mAQAAQJCdcfPn9/s1f/58J7MAAAAgyGJrmvjmm2+edFpFRYXjYQAAABBcNTZ/b731lrp06aIGDRpUmWZZVtBCAQAAIDhqbP5SUlI0YMAApaWlVZl2+PBhffjhh0ELVlukpLQM8RpDvb66jno6h1o6i3o6i3o6h1o6K3T1tLtdrsbmr2vXrjpw4EC102JiYtS7d+/TDhZpduzYGbJ1JSUlqaCgIGTrq+uop3OopbOop7Oop3OopbNCX097jWaNzd/w4cNPOi0mJkajRo06vUwAAAAIK67zBwAAEEVo/gAAAKIIzR8AAEAUofkDAACIIjR/AAAAUcR28zd9+vTA86lTpwYlDAAAAILLdvP39ddfV/scAAAAkYPdvgAAAFGE5g8AACCK0PwBAABEEZo/AACAKHJGzZ9lWU7nAAAAQAjYbv5+8YtfBJ63b98+KGEAAAAQXLabvwkTJgSeT5w4MShhAAAAEFy2mr9t27apoKCg0lhBQYG2bdsWjEwAAAAIElvN37PPPiu/319prKKiQnPnzg1KKAAAAASHreavoKBAycnJlcZatGihvXv3BiUUAAAAgsNW85eYmKitW7dWGtu6dasSEhKCEgoAAADBEWtnpuzsbM2YMUNXX321kpOTtWfPHv31r3/V4MGDg50PAAAADrLV/GVlZSk+Pl4rVqxQYWGhvF6vbr75ZnXv3t2REAUFBZo3b5727dsnwzCUlZWlq666SiUlJZo1a5b27t2rZs2a6Z577pHH45FlWVq8eLE++eQT1a9fX6NGjVLbtm0lSStXrtTSpUslSYMHD1afPn0kHd1SOW/ePB0+fFidO3fWyJEjZRiGI/nPlrl3t/TOK/KVFsuMbyRdc4NczVqEOxYAAKiDbDV/ktSjRw/16NEjKCFiYmJ00003qW3btiovL9eECRPUqVMnrVy5Ur/85S81aNAgvf3223r77bd144036pNPPtHu3bs1Z84cbdmyRYsWLdKTTz6pkpISvfnmm5o2bZqko5enycjIkMfj0fPPP6/bb79d559/vqZOnarc3Fx17tw5KJ/ndJh7d8ua9Yi0d7eOHBvcuknmPY/TAAIAAMfVitu7JSQkBLbcNWzYUCkpKfL5fFq7dq169+4tSerdu7fWrl0rSVq3bp169eolwzCUmpqq0tJSFRUVKTc3V506dZLH45HH41GnTp2Um5uroqIilZeXKzU1VYZhqFevXoH3Crt3XpH27q489uOWQAAAAKfZ3vIXKvn5+fr222/Vrl077d+/P3BSSUJCgg4cOCBJ8vl8SkpKCizj9Xrl8/nk8/nk9XoD44mJidWOH5u/OsuXL9fy5cslSdOmTau0nmDwlRb/tMXvOLGlxUoM8rrrutjY2KD//UULauks6uks6ukcaums2lrPWtX8HTx4UDNnztSIESPkdrtPOl919xY+2fF7hmGc1r2Is7KylJWVFXh94sWtnWbGN6p2/EiDhkFfd12XlJREDR1CLZ1FPZ1FPZ1DLZ0V6nq2bNnS1ny2dvt+9NFH1Y6vWbPGfqJTqKio0MyZM3XppZfq4osvliQ1adJERUVFkqSioiI1btxY0tEtd8cXs7CwUAkJCUpMTFRhYWFg3OfzKSEhQV6vt9J4YWGhEhMTHct+Vq65Qaru2D5foawjh0OfBwAA1Gm2mr+FCxdWO/7cc885EsKyLC1cuFApKSkaOHBgYDwjI0PvvfeeJOm9995T165dA+OrVq2SZVnavHmz3G63EhISlJaWpo0bN6qkpEQlJSXauHGj0tLSlJCQoIYNG2rz5s2yLEurVq1SRkaGI9nPlqtZCxn3PC7j4t6Ku7CLjIt7S0NvkbZ+LfOPf5B1wp1VAAAAzkaNu3337NkjSTJNU/n5+ZV2n+7Zs0f16tVzJMSmTZu0atUqnXPOObrvvvskSddff70GDRqkWbNmacWKFUpKStL48eMlSZ07d9aGDRs0duxY1atXT6NGjZIkeTweXXvttZo4caIkaciQIfJ4PJKkW2+9VfPnz9fhw4eVlpZWK870PcbVrIV06++UeNzmYTMmVtZrf5T14hxpxN0yXLXi3BwAABDhDKuGA+KGDx9+0gWbNm2qoUOHVjo+ri7auXNnyNZ14rEB5l9fk/Xun2X0/5WM4bfWmusSRgqOXXEOtXQW9XQW9XQOtXRWbT3mr8Ytf6+//rokadKkSXrsscfOPhVOizFwuFRWImv5u5LbI+Pq68MdCQAARDhbZ/vS+IWHYRhHj/8rK5X111dluuPlyro63LEAAEAEq3G37zH5+fl69dVXtW3bNh08eLDStAULFgQtXG3AnlYAABAJ7F7ZztaWv2eeeUbJycm6+eabVb9+/bPJBQAAgDCy1fxt375dkydPlisKzzjdsSN8J3ycyDpYLvPph6Uftso1dpKM9heFLFsk4sBl51BLZ1FPZ1FP51BLZ4W+ng5e5Ll9+/batm3b2aSBA4wGDeW6e5KUnCJz3hRZWzeFOxIAAIgwtrb8NWvWTFOmTFG3bt3UtGnTStNquhwMnGfEN5Jr3GMyn5ogc87jct33pIyUc8MdCwAARAhbW/4OHTqk9PR0+f1+FRYWVnog9IymiXLd87gUGydz1iRZe3eHOxIAAIgQtrb8HbuDBmoPo1kLue55XOaMiTJnPSLX76fJaFpL7lcMAABqLdtncGzfvl1vvvmmcnJyJB2988V3330XtGA4NSPlnKPHAB7YJ3P2JFmlxeGOBAAAajlbzd9HH32kSZMmyefzadWqVZKk8vJyvfTSS0ENh1Mz2qTKNfpBac8Omc88JutgebgjAQCAWsxW8/eXv/xFDz/8sG677bbA5V7OPfdczgCuJYz2F8l12++l7/Jkzn9S1pEj4Y4EAABqKVvN3/79+3XuuZXPKDUM4+jtx1ArGJ27y/j1WOmrjTKfnyHL7w93JAAAUAvZav7atm0b2N17zIcffqh27doFJRTOjCuzn4zrfit9skbWS3NlmWa4IwEAgFrG1tm+I0eO1BNPPKEVK1bo0KFDmjJlinbu3KmHHnoo2Plwmlz9fyWztETWX1+V3PHSsN+whRYAAATYav5SUlI0e/ZsrV+/Xunp6fJ6vUpPT1eDBg2CnQ9nwPjVdVJZiazl70rxHhkDrwt3JAAAUEvYav7++9//Kj09XZmZmcHOAwcYhiEN+41UVirrnT/LbOiRq//AcMcCAAC1gK3m74033tCCBQt08cUX69JLL1XHjh2DnQtnyXC5pF+PkVVeJuu1P8p0x8vVo2+4YwEAgDCz1fzNmDFD27dv1wcffKCFCxfqyJEjyszMVM+ePdW2bdtgZ8QZMmJi5LrtXplzHpf1p2dkNWwoI617uGMBAIAwsn2Hj1atWum6667Ts88+q/Hjx+v777/XxIkTg5kNDjDi6sk1+gHp3HYyn5sh6+tPwx0JAACEke3mT5IKCgr07rvvatGiRfrmm2/Uty+7ESOB0cAt19hHpOY/kzl3iqxvt4Q7EgAACBNbu33/9a9/6YMPPtB3332nzp07a8iQIerSpYtiY20tjlrA8DSW657HZE6fIPOZR+X6/VQZLc8JdywAABBitrq39evXa8CAAerWrRuXd4lgRlOvXPc8LvOpCTJnPSLX76fJaNYi3LEAAEAI2drt+8ADD6hXr16qV6+eioqKgp0JQWQ0/5lc9zwuHT4sc9Yjsvb5wh0JAACEkK3mr7S0VM8884xuuOEGjR07VpK0bt06vfbaa0ENh+AwUs49egzggX0yZ0+SVVoS7kgAACBEbDV/zz//vNxut+bPnx84zi81NVWrV68OajgEj3HeBXKNekDas0Pms4/LOnQw3JEAAEAI2Gr+PvvsM40cOVIJCQmBscaNG2v//v1BC4bgMzqkyfXbe6Wtm2XOf1LWkSPhjgQAAILMVvPndrtVXFxcaaygoKBSM4jIZHTJlPHrMdKXuTIXzZTl94c7EgAACCJbzV///v01c+ZMff7557IsS5s3b9a8efM0YMCAYOdDCLgu6S9j+G+kDatlLZkny7LCHQkAAASJrUu9XHPNNYqLi1NOTo78fr8WLFigrKwsXXXVVcHOhxBxZV0js7RU1t9ek9zx0tBbZBhGuGMBAACH2Wr+DMNQdna2srOzg50HYWRcfb1UViLr/70jxTeSkT0s3JEAAIDDuEUHAgzDkIbfKpWVynr7ZZnueLn60vADAFCX0PyhEsPlkn49RlZ5qaw/PyezYbxc3fuEOxYAAHCIrRM+EF2M2Fi5bv+99Itfylo8W9bG/4Y7EgAAcAjNH6plxNWT664HpXPOk/ncU7I2fR7uSAAAwAEn3e27YsUKW2/Qr18/x8KgdjEauOUaO0nmjIky506W63dPyPj5+eGOBQAAzoJhneSibo899ljguWVZ2rRpk5o2bSqv16vCwkLt27dPF1xwgSZNmuRIkPnz52vDhg1q0qSJZs6cKUkqKSnRrFmztHfvXjVr1kz33HOPPB6PLMvS4sWL9cknn6h+/foaNWqU2rZtK0lauXKlli5dKkkaPHiw+vTpI0naunWr5s2bp8OHD6tz584aOXKkrUuZcLUTAAAQCWxfpteyIScnx/rb3/5Waezvf/+7lZOTY2dxW7744gvrm2++scaPHx8YW7JkibVs2TLLsixr2bJl1pIlSyzLsqz169dbU6ZMsUzTtDZt2mRNnDjRsizLKi4utkaPHm0VFxdXem5ZljVhwgRr06ZNlmma1pQpU6wNGzbYynW0lDx48ODBgwcPHrX7YZetY/7ef/99XXnllZXGrrjiCr3//vun25SeVIcOHeTxeCqNrV27Vr1795Yk9e7dW2vXrpUkrVu3Tr169ZJhGEpNTVVpaamKioqUm5urTp06yePxyOPxqFOnTsrNzVVRUZHKy8uVmpoqwzDUq1evwHsBAABEE1uXemnatKnWrVunbt26BcbWrVunxo0bBy2YJO3fvz9w/+CEhAQdOHBAkuTz+ZSUlBSYz+v1yufzyefzyev1BsYTExOrHT82f3WWL1+u5cuXS5KmTZumQ4cOO/65TiY2NlYVFRUhW9+ZOPz1Zyp69G7F/qy1Ep6YK1d8o3BHOqlIqGekoJbOop7Oop7OoZbOCn0969may1bzN3LkSM2cOVPvvvuuvF6vCgoKtH37do0fP/6sIp4py7KqjJ3s+D3DMKqd/2SysrKUlZUVeF1QUHD6Ac9QUlJSSNd3RpJ+JteoiaqYM1l7J90t1z2Py6jfINypqhUR9YwQ1NJZ1NNZ1NM51NJZoa5ny5Ytbc1na7dvp06dNHfuXF122WVq06aNLrvsMs2dO1cXXXTRWYU8lSZNmqioqEiSVFRUFNjSeKwBPaawsFAJCQlKTExUYWFhYNzn8ykhISFwksrx8ycmJgY1e11mdOgs12/vlbZulrlgqqyKI+GOBAAAbLJ9nb9GjRqpV69eGjRokHr37q1GjYK/uy8jI0PvvfeeJOm9995T165dA+OrVq2SZVnavHmz3G63EhISlJaWpo0bN6qkpEQlJSXauHGj0tLSlJCQoIYNG2rz5s2yLEurVq1SRkZG0PPXZUZ6poybR0tffCJr0dOyTH+4IwEAABts7fbNz8/Xq6++qm3btungwYOVpi1YsMCRILNnz9aXX36p4uJi3XHHHRo2bJgGDRqkWbNmacWKFUpKSgrsZu7cubM2bNigsWPHql69eho1apQkyePx6Nprr9XEiRMlSUOGDAmcRHLrrbdq/vz5Onz4sNLS0tS5c2dHckczV88BMstKZb3xgvRyvHTTaFuXzwEAAOFz0uv8He/BBx9UcnKyLr30UtWvX7/StA4dOgQtXG2wc+fOkK0rUo+1MN9+Wdbf/yLj8v+Rce2IWtMARmo9ayNq6Szq6Szq6Rxq6azaesyfrS1/27dv1+TJk+VycTc4VGVcc4NUViLrX8skt0fGVUPDHQkAAJyErW6uffv22rZtW5CjIFIZhiHjuttkXNxb1rIlMlf+I9yRAADASdja8tesWTNNmTJF3bp1U9OmTStNGz58eFCCIbIYLpc04m5ZB8tl/fk5mQ3j5bq4d7hjAQCAE9ja8nfo0CGlp6fL7/ersLCw0gM4xoiNleu2+6TzO8paPFvWp9xFBQCA2sbWlr9jZ9MCp2LUqy/XXQ/JnPmQzIXT5Rr3qIzUC8MdCwAA/Oi0zuAoLy9Xfn6+9uzZE3gAJzIauuW6+1HJ21zm3CdkffdNuCMBAIAf2T7bd86cOfruu++qTHv99dcdD4XIZzRqLNc9j8t8aoLM2ZPk+v00GT9rFe5YAABEPVtb/hYtWqSOHTvqhRdekNvt1uLFizVgwACNHj062PkQwYzEJLnueVxyuWTOekRW4d5wRwIAIOrZav6+++473XDDDYqPj5dlWXK73brxxhvZ6odTMpJbyjXuMelQ+dEG8MC+cEcCACCq2Wr+4uLi5PcfvXdro0aNVFBQIMuyVFJSEtRwqBuM1m3kGvOIVLRX5uxJssr4dwMAQLjYav4uuOACffTRR5Kk7t2768knn9Sjjz6qjh07BjUc6g6jXXu57nxA2vmDzGefkHXoULgjAQAQlWyd8DF+/PjA8+uvv16tW7fWwYMH1atXr6AFQ91jXNhFrlvHy/zjH2QunCrX6AdlxMaFOxYAAFHFVvN3PJfLRdOHM2Zk9JRRXibrpbmyXpgt3Tpehism3LEAAIgap938AWfLdellMstKZb25WGrolm4cJcMwwh0LAICoQPOHsHBd/j8yy0pk/eMNye2Rce2vwx0JAICoQPOHsDEG3SiVlcj651sy3R65rrw23JEAAKjzbDd/JSUlOnjwoBo0aCCPxxPMTIgShmFI198ulZXKWvqiTHe8XL2vCHcsAADqtBqbv4qKCv3lL3/RypUrtX///sB406ZN1adPHw0dOlSxsWw8xJkzXC5p5DhZ5WWyXllwtAHsemm4YwEAUGfV2LktWrRIe/bs0dixY3XuuefK7XarvLxc27Zt09KlS7Vo0SLdcccdocqKOsqIjZXrjvtlPvOorJynZTVwy/hlerhjAQBQJ9V4keePP/5Y9913ny688EI1atRIMTEx8ng8uvDCCzV+/HitWbMmVDlRxxn16ss1+iEp5ecyF06VtfmLcEcCAKBOqrH5i4uLU1FRUbXT9u3bp7g4LtAL5xjueLnGPSolNpM5d7Ks778JdyQAAOqcGnf7Xn311XrsscfUr1+/Krt9//Of/2jQoEGhyokoYTRqItc9j8ucPkHm7Efl+v1UGS1ahTsWAAB1Ro3N38CBA9WqVSutWrVK69evD5zt27p1a915551KS0sLVU5EESOx2dEG8KkJMmc9Itf902UkNgt3LAAA6oRTnqqblpZGk4eQM1qkyDXuMZl/ePBoA3jfVBmNm4Y7FgAAEc+wLMsKd4jajLuOAQCASGC3o6vxhA8AAADULTR/AAAAUYTmDwAAIIrYvjfbpk2b9Itf/EKS9PXXX+uCCy4IWqjaZMeOnSFbV1JSkgoKCkK2vkhk/vMtWW+9KKP3FTJuuPPo/YFPgno6h1o6i3o6i3o6h1o6K/T1bGlrLttb/qZNmxZ4PnXq1NPPAzjAdcW1Mq68VtZ7/5S1bEm44wAAEHFsb/k7HicII5yM/7lZKi2V9X9vyoz3yHX54HBHAgAgYpxR81fTrjYg2AzDkG64XSovlfXmn2Q2jJer1+XhjgUAQEQ4o+YPCDfDFSPdMk5WeZmsl+fLcsfLyOgZ7lgAANR6nO2LiGXExsl1xwSpXXuZi56W9fn6cEcCAKDWo/lDRDPq15frroellq1lLpgqK+/LcEcCAKBW44QPRDzDHX/0PsBPTZQ5Z7L0m/Ey1q6Sr7RYZnwj6Zob5GrWItwxAaDWMvfult55hZ+bDqnt9bR9b9/jr+0Xqdf5y83N1eLFi2Wapvr3769BgwadcpmdO7nOX6SwCvfKnHqvdGC/ZJk/TWjWQsY9j9eq/3iR4tgPsNjSYlXUwh9gkYZ6Oot6OsPcu1vWrEekvbt/GuTn5hkLZz1btrR3nT/bzV+kM01Td999tx566CF5vV5NnDhRd999t1q1alXjcjR/kcX/7GTp07VVJ/zil3JdNVRyuY4+DNdPzwOvDcmIOW7MqDzP8cud+OePj7p0JjxfCM6ins4KdT0ty5KOf+iE18fGzBOn6cdfRo97fvyfsiTzx19WA9Nqer+TrLvaZY695/HLnbhuS+a/l0qbv6j6odu1l9Fv4MnXY2vdp/P5bX7OauurH9+nunUfe+/q1l11PVaN6z6hpsfGjv8Mu36QivdXKadxcW+5bv1dzf/QzpLd5i9qzvbNy8tTixYtlJycLEnKzMzU2rVrT9n8IcIcOlj9+KbPZG76LPjrN4zKTaMR82NTeZImslLjeZLpJ208Kzenhivmp7Ga3stmPuu/qyp/sUrS3t2y/jhD5iVZwa9lHWN9uLz6ej43XWb3vjV/sZ3sy/v4LzGbX2yn/FI96Ze3fvxirG7dNb2P3XXrhCbIqrLuSo/9PulgedV6PjpGfrenhnWfYf2iUd5XsvK+Cm8GwyUZx/0p4+jPLBk//bw99qhu7GTj0tGfddKPY8etJzBm8/0D035cvqKi2o9i7fMFrUynK2qaP5/PJ6/XG3jt9Xq1ZcuWKvMtX75cy5cvl3T0riZJSUkhyxgbGxvS9dVF+5N/poPVNHn1Luqm+GEjjn65mKasH/+U6T/63DKPm+Y/+kP/+OnHz+M/towlmf6f3s8yJX/1y1g/vtfR9/VXWr8s67jpVfPJNI/+Jnpsur/iuOln9hmswLp+fN/TsW2LrG1V/+/gDH33jazvvjmzZav7Avrxi9FwuVTpC+mEeY3A82Nfesc9DyyvH78gq37ZVV3+2DqOPjcqfcEaUkzl9zGO/RKi49d9knUc/4XvOrq+w59tkHVi8yfJcHtUP73Hcfkqf65jyxvH1nd8xuOaimqXP+4zGifmPSGjcexz1fD3YD9j9Q2PUd3f3XHrME6oZ+C9f8x2bPmSVxfp8CdrqtSyXnqmGv169AkZj9bMcNmtX+V/F7brF8F7UfbPelQHV/27yniD5J+pSS35jo+a5q+6vdvV/ePKyspSVtZPWzVCuRuW3b5nz7xiiPTVp1V2BR0ZfqsORNGutRP/ZZ/qx2jl5tCSLL/MPz0rrf+w6syde8h1wx1ORY0a5isLpU8+qjqhS6ZcN9913Bf/SbYwVPniD++X46m2hQV9W9mimdLH71Ud/8WFOnLdbcFee51iDr1F2r6t6s/Na0doX8NGp/FGx578uMVUknSav1zWAScWW3NRAAAZK0lEQVT7Hjp0xZCgf8c7ttvXNE3Nnz9ft99+u+Li4s46WLh4vV4VFhYGXhcWFiohISGMiRAMrmYtZN7zOAeBnybj2O7e48eu/bWs77+pekzV0JEymvB/53QZQ0fK2v5t1XoOGSEj3hO+YJHqmhukrZuq1FPX3BC+TBGKn5vOioR6nrL5c7lc+vTTT8P+W+bZOu+887Rr1y7l5+crMTFRq1ev1tixY8MdC0HgatZCuvV3SmRL6lmJhB9gkYR6Oot6Ooufm86q7fW0dbbvO++8o9LSUg0bNkyxsZG7p3jDhg168cUXZZqm+vbtq8GDB59yGc72jVzU0znU0lnU01nU0znU0lmhrqejZ/v+85//1L59+/T3v/9djRs3rjRtwYIFp58uTLp06aIuXbqEOwYAAEDY2Gr+xowZE+wcAAAACAFbzV+HDh2CnQMAAAAh4Dr1LAAAAKgraP4AAACiCM0fAABAFDmt5s80TRUVFQUrCwAAAILM1gkfpaWlWrRokdasWaPY2FgtWbJE69atU15enq677rpgZwQAAIBDbG35e/755+V2uzV//vzARZ5TU1O1evXqoIYDAACAs2xt+fvss8/03HPPVbq7R+PGjbV///6gBQMAAIDzbG35c7vdKi4urjRWUFCghARu7g4AABBJbDV//fv318yZM/X555/Lsixt3rxZ8+bN04ABA4KdDwAAAA6ytdv3mmuuUVxcnHJycuT3+7VgwQJlZWXpqquuCnY+AAAAOMhW82cYhrKzs5WdnR3sPAAAAAgiW82fJOXn5+v777/XwYMHK4337NnT8VAAAAAIDsOyLOtUMy1btkxvvvmmWrdurXr16v20sGHoscceC2rAcDOMcCcAAAA4tVN3dEfZ2vL3t7/9TdOnT1erVq3OJhMAAADCzNbZvh6PR82aNQt2FgAAAASZrd2+n3zyid5//31lZ2erSZMmlaYlJSUFLVxtsHPnzpCtKykpSQUFBSFbX11HPZ1DLZ1FPZ1FPZ1DLZ0V6nq2bNnS1ny2dvtWVFTo008/1Ycfflhl2uuvv356yQAAABA2tpq/RYsW6frrr9cll1xS6YQPAAAARBZbzZ9pmurbt69cLluHCAIAAKCWstXN/epXv9Lbb78tG4cHAgAAoBazteXv//7v/7Rv3z4tW7ZMHo+n0rQFCxYEJRgAAACcZ6v5GzNmTLBzAAAAIARsNX8dOnQIdg4AAACEwEmbv6VLl2rw4MGSar6cy/Dhw51PBQAAgKA4afNXWFhY7XMAAABErpM2f7/97W/19ddf64ILLtCoUaNCmQkAAABBUuOlXqZOnRqqHAAAAAiBGps/rusHAABQt9R4tq9lWcrPz6+xCUxOTnY8FAAAAIKjxubv8OHDp7zGX01nAgMAAKB2qbH5q1+/vl566aVQZQEAAECQ1XjMn2EYocoBAACAEOCEDwAAgChSY/P39NNPhyoHAAAAQqDGY/6SkpKCHuCjjz7SG2+8oR07dujJJ5/UeeedF5i2bNkyrVixQi6XSyNHjlRaWpokKTc3V4sXL5Zpmurfv78GDRokScrPz9fs2bNVUlKiNm3aaMyYMYqNjdWRI0c0d+5cbd26VY0aNdK4cePUvHnzoH82AACA2qbGLX+h0Lp1a917771q3759pfHt27dr9erVevrpp/Xggw8qJydHpmnKNE3l5OTogQce0KxZs/Thhx9q+/btkqSXX35Z2dnZmjNnjuLj47VixQpJ0ooVKxQfH69nn31W2dnZeuWVV0L+OQEAAGoDW82faZpBC9CqVSu1bNmyyvjatWuVmZmpuLg4NW/eXC1atFBeXp7y8vLUokULJScnKzY2VpmZmVq7dq0sy9IXX3yh7t27S5L69OmjtWvXSpLWrVunPn36SJK6d++uzz//nOMZAQBAVDpl82eapm666SYdOXIkFHkCfD6fvF5v4HViYqJ8Pl+Vca/XK5/Pp+LiYrndbsXExFSa/8T3iomJkdvtVnFxcQg/DQAAQO1Q4zF/kuRyudSyZUsVFxcrMTHxjFYyefJk7du3r8r4ddddp65du1a7zMm2zFU3fqpL0pzOMsuXL9fy5cslSdOmTQvJcY/HxMbGhnR9dR31dA61dBb1dBb1dA61dFZtrecpmz9J6tmzp6ZPn64rr7xSXq+3UuN04YUXnnL5hx9++LSDeb1eFRYWBl77fL5A83n8eGFhoRISEtSoUSOVlZXJ7/crJiam0vzH3svr9crv96usrEwej6fa9WZlZSkrKyvwuqCg4LSzn6mkpKSQrq+uo57OoZbOop7Oop7OoZbOCnU9qzuMrjq2mr9///vfkqQ33nij0rhhGJo7d+5pRrMnIyNDc+bM0cCBA1VUVKRdu3apXbt2sixLu3btUn5+vhITE7V69WqNHTtWhmGoY8eOWrNmjS655BKtXLlSGRkZkqT09HStXLlSqampWrNmjTp27MgFrAEAQFQyrDCf+fDf//5XL7zwgg4cOKD4+Hj9/Oc/14MPPihJWrp0qf7zn//I5XJpxIgR6ty5syRpw4YNevHFF2Wapvr27avBgwdLkvbs2VPlUi9xcXE6fPiw5s6dq2+//VYej0fjxo1TcnKyrXw7d+4MzgevBr9xOYt6OodaOot6Oot6OodaOqu2bvmz3fz5/X5t2rQpcPJEampq4OSKuozmL3JRT+dQS2dRT2dRT+dQS2fV1ubP1m7fHTt2aPr06Tp8+HDg+Lm4uDjdf//9atWq1VkFBQAAQOjYav4WLVqkrKws/epXvwocK/fuu+8qJydHkyZNCmpAAAAAOMfWRZ63bdumgQMHVjpJIjs7W9u2bQtWLgAAAASBrS1/iYmJ+vLLLytd1uWrr75SQkJC0ILVFikp9vafOyfU66vrqKdzqKWzqKezqKdzqKWzQldPu6fw2mr+rr/+ek2fPl3p6emBgxc3bNigMWPGnE1GAAAAhJjts3137typjz76SEVFRUpISFCPHj1sn1USyTjbN3JRT+dQS2dRT2dRT+dQS2dF3Nm+t99+u5577jlJ0vz58zVq1Chde+21zqQDAABAWJz0hI+KigoVFxdLkj7++OOQBQIAAEDwnHTL34ABA3TnnXeqUaNGOnTokO68885q51uwYEHQwgEAAMBZJ23+rrvuOg0YMEB79+7VE088wckdAAAAdUCNZ/t6vV55vV7df//96tChQ6gyAQAAIEhsXeT5l7/8ZbBzAAAAIARsNX8AAACoG2j+AAAAooit5u+jjz6qdnzNmjWOhgEAAEBw2Wr+Fi5cWO34sYtAAwAAIDLUeLbvnj17JEmmaSo/P1/H3wluz549qlevXnDTAQAAwFE1Nn9jx44NPD/xOn9NmzbV0KFDg5MKAAAAQVFj8/f6669LkiZNmqTHHnssJIEAAAAQPLaO+aPxAwAAqBtq3PJ3TH5+vl599VVt27ZNBw8erDSNe/sCAABEDlvN3zPPPKPk5GTdfPPNql+/frAzAQAAIEhsNX/bt2/X5MmT5XJxTWgAAIBIZquba9++vbZt2xbkKAAAAAg2W1v+mjVrpilTpqhbt25q2rRppWnDhw8PSjAAAAA4z1bzd+jQIaWnp8vv96uwsDDYmQAAABAktpq/UaNGBTsHAAAAQsBW83fsNm/VSU5OdiwMAAAAgstW83f8bd5OdOwuIAAAAKj9bDV/JzZ4+/bt0xtvvKH27dsHJRQAAACC44wu3Ne0aVONGDFCf/7zn53OAwAAgCA646s279y5U4cOHXIyCwAAAILM1m7fRx55RIZhBF4fOnRIP/zwg4YMGRK0YAAAAHCereavX79+lV43aNBA5557rn72s58FJRQAAACCw1bz16dPnyDHAAAAQCjYav4qKiq0dOlSrVq1SkVFRUpISFCvXr00ePBgxcbaegsAAADUArY6t5dfflnffPONfvvb36pZs2bau3ev3nrrLZWVlWnEiBFBjggAAACn2Gr+1qxZoxkzZqhRo0aSpJYtW6pNmza67777zrr5W7JkidavX6/Y2FglJydr1KhRio+PlyQtW7ZMK1askMvl0siRI5WWliZJys3N1eLFi2Wapvr3769BgwZJkvLz8zV79myVlJSoTZs2GjNmjGJjY3XkyBHNnTtXW7duVaNGjTRu3Dg1b978rHIDAABEJMuG2267zTpw4EClsf3791u33XabncVrlJuba1VUVFiWZVlLliyxlixZYlmWZf3www/Wvffeax0+fNjas2ePddddd1l+v9/y+/3WXXfdZe3evds6cuSIde+991o//PCDZVmWNXPmTOuDDz6wLMuynnvuOetf//qXZVmW9c9//tN67rnnLMuyrA8++MB6+umnbeeTePDgwYMHDx48av/DLltb/nr06KHp06dryJAhSkpKUkFBgd566y316NHjrJvPiy66KPA8NTVVa9askSStXbtWmZmZiouLU/PmzdWiRQvl5eVJklq0aBG4p3BmZqbWrl2rlJQUffHFF7r77rslHT1J5Y033tBll12mdevWaejQoZKk7t2764UXXpBlWZUuX3MyO3bsPOvPaNex2sIZ1NM51NJZ1NNZ1NM51NJZoa9nS1tz2Wr+brzxRr311lvKyclRUVGREhMTlZmZqWuvvfasIp5oxYoVyszMlCT5fD6df/75gWmJiYny+XySJK/XGxj3er3asmWLiouL5Xa7FRMTU2V+n88XWCYmJkZut1vFxcVq3LhxlQzLly/X8uXLJUnTpk1TUlKSo5+xJrGxsSFdX11HPZ1DLZ1FPZ1FPZ1DLZ1VW+tpq/mLjY3V8OHDNXz48DNayeTJk7Vv374q49ddd526du0qSVq6dKliYmJ06aWXSpIsy6r2vaobP9UWvNNZJisrS1lZWYHXoezY+Y3LWdTTOdTSWdTTWdTTOdTSWaGuZ8uWDmz5+/rrr7Vu3TrdeOONVaa98sor6tq1q1JTU0+5kocffrjG6StXrtT69esr3UnE6/WqsLAwMI/P51NiYqIkVRovLCxUQkKCGjVqpLKyMvn9fsXExFSa/9h7eb1e+f1+lZWVyePxnDI3AABAXVPjvX2XLVumDh06VDutQ4cOWrp06VkHyM3N1TvvvKP7779f9evXD4xnZGRo9erVOnLkiPLz87Vr1y61a9dO5513nnbt2qX8/HxVVFRo9erVysjIkGEY6tixY+CYwZUrVyojI0OSlJ6erpUrV0o6euZyx44dbR3vBwAAUNfUuOVv27ZtgcurnKhTp05auHDhWQfIyclRRUWFJk+eLEk6//zzddttt6l169bq0aOHxo8fL5fLpd/85jdyuY72qrfccoumTJki0zTVt29ftW7dWpJ0ww03aPbs2XrttdfUpk2bwG3p+vXrp7lz52rMmDHyeDwaN27cWecGAACIRDU2f+Xl5aqoqFC9evWqTPP7/SovLz/rAM8+++xJpw0ePFiDBw+uMt6lSxd16dKlynhycrKmTp1aZbxevXoaP3782QUFAACoA2rc7ZuSkqKNGzdWO23jxo1KSUkJSigAAAAER43NX3Z2tv74xz/q448/lmmakiTTNPXxxx/r+eefV3Z2dkhCAgAAwBk17vbt2bOn9u3bp3nz5unIkSNq3LixDhw4oHr16mno0KHq2bNnqHICAADAAae8zt/AgQPVr18/bd68WSUlJfJ4PEpNTZXb7Q5FPgAAADjI1kWe3W73Sc/6BQAAQOSo8Zg/AAAA1C00fwAAAFGE5g8AACCK0PwBAABEEZo/AACAKELzBwAAEEVo/gAAAKIIzR8AAEAUofkDAACIIjR/AAAAUYTmDwAAIIrQ/AEAAEQRmj8AAIAoQvMHAAAQRWj+AAAAogjNHwAAQBSh+QMAAIgiNH8AAABRhOYPAAAgitD8AQAARBGaPwAAgChiWJZlhTsEAAAAQoMtf7XIhAkTwh2hTqGezqGWzqKezqKezqGWzqqt9aT5AwAAiCI0fwAAAFEk5tFHH3003CHwk7Zt24Y7Qp1CPZ1DLZ1FPZ1FPZ1DLZ1VG+vJCR8AAABRhN2+AAAAUSQ23AFwVG5urhYvXizTNNW/f38NGjQo3JEi1vz587VhwwY1adJEM2fODHeciFZQUKB58+Zp3759MgxDWVlZuuqqq8IdK2IdPnxYkyZNUkVFhfx+v7p3765hw4aFO1ZEM01TEyZMUGJiYq09szJSjB49Wg0aNJDL5VJMTIymTZsW7kgRq7S0VAsXLtQPP/wgwzB05513KjU1NdyxAmj+agHTNJWTk6OHHnpIXq9XEydOVEZGhlq1ahXuaBGpT58+uuKKKzRv3rxwR4l4MTExuummm9S2bVuVl5drwoQJ6tSpE/82z1BcXJwmTZqkBg0aqKKiQo888ojS0tJq1ZdCpPnHP/6hlJQUlZeXhztKnTBp0iQ1btw43DEi3uLFi5WWlqbf/e53qqio0KFDh8IdqRJ2+9YCeXl5atGihZKTkxUbG6vMzEytXbs23LEiVocOHeTxeMIdo05ISEgIHKzcsGFDpaSkyOfzhTlV5DIMQw0aNJAk+f1++f1+GYYR5lSRq7CwUBs2bFD//v3DHQUIKCsr01dffaV+/fpJkmJjYxUfHx/mVJWx5a8W8Pl88nq9gdder1dbtmwJYyKgqvz8fH377bdq165duKNENNM0df/992v37t26/PLLdf7554c7UsT605/+pBtvvJGtfg6aMmWKJGnAgAHKysoKc5rIlJ+fr8aNG2v+/Pn67rvv1LZtW40YMSLwi19twJa/WqC6E67ZGoDa5ODBg5o5c6ZGjBght9sd7jgRzeVyacaMGVq4cKG++eYbff/99+GOFJHWr1+vJk2a1MrLaESqyZMna/r06XrggQf0r3/9S19++WW4I0Ukv9+vb7/9Vpdddpmeeuop1a9fX2+//Xa4Y1VC81cLeL1eFRYWBl4XFhYqISEhjImAn1RUVGjmzJm69NJLdfHFF4c7Tp0RHx+vDh06KDc3N9xRItKmTZu0bt06jR49WrNnz9bnn3+uOXPmhDtWREtMTJQkNWnSRF27dlVeXl6YE0Umr9crr9cb2KrfvXt3ffvtt2FOVRnNXy1w3nnnadeuXcrPz1dFRYVWr16tjIyMcMcCZFmWFi5cqJSUFA0cODDccSLegQMHVFpaKunomb+fffaZUlJSwpwqMv3v//6vFi5cqHnz5mncuHG68MILNXbs2HDHilgHDx4M7D4/ePCgPv30U51zzjlhThWZmjZtKq/Xq507d0qSPvvss1p3khzH/NUCMTExuuWWWzRlyhSZpqm+ffuqdevW4Y4VsWbPnq0vv/xSxcXFuuOOOzRs2LDAgbc4PZs2bdKqVat0zjnn6L777pMkXX/99erSpUuYk0WmoqIizZs3T6ZpyrIs9ejRQ+np6eGOBWj//v36wx/+IOnobsuePXsqLS0tzKki1y233KI5c+aooqJCzZs316hRo8IdqRLu8AEAABBF2O0LAAAQRWj+AAAAogjNHwAAQBSh+QMAAIgiNH8AAABRhOYPAAAginCdPwA4Q6NHj9a+ffsUExMjl8ulVq1aqVevXsrKypLLxe/WAGonmj8AOAv333+/OnXqpLKyMn355ZdavHix8vLyat1FXQHgGJo/AHCA2+1WRkaGmjZtqgcffFADBw5UQUGBXnvtNe3Zs0dut1t9+/bVsGHDJElTp05VWlqarrzyysB73HvvvRo2bJi6du2qF198UR988IGOHDmiZs2aaezYsdxuC4AjaP4AwEHt2rVTYmKivv76a6WkpOiuu+5Sq1at9MMPP+iJJ57Qz3/+c3Xr1k29e/fW3/72t0Dzt23bNvl8PnXp0kUbN27UV199pWeeeUZut1s7duxQfHx8mD8ZgLqCg1IAwGGJiYkqKSlRx44ddc4558jlcuncc8/VJZdcoi+//FKS1LVrV+3atUu7du2SJK1atUqZmZmKjY1VbGysDh48qB07dsiyLLVq1UoJCQnh/EgA6hC2/AGAw3w+nzwej7Zs2aI///nP+v7771VRUaGKigp1795dkhQXF6cePXro/fff15AhQ/Thhx/qd7/7nSTpwgsv1OWXX66cnBwVFBSoW7duuummm+R2u8P5sQDUEWz5AwAH5eXlyefz6YILLtCcOXOUnp6uBQsW6MUXX9SAAQNkWVZg3j59+uj999/X559/rvr16ys1NTUw7aqrrtL06dP19NNPa9euXXr33XfD8XEA1EE0fwDggLKyMq1fv17PPPOMLr30Up1zzjkqLy+Xx+NRvXr1lJeXpw8++KDSMqmpqXK5XHrppZfUq1evwHheXp62bNmiiooK1a9fX3FxcVw6BoBj2O0LAGdh+vTpiomJkWEYatWqlbKzs3XZZZdJkm699Va99NJLeuGFF9ShQwf16NFDpaWllZbv1auXXn/9dd13332BsfLycr344ovas2eP6tWrp4suukhXX311SD8XgLrLsI7fBwEACKn33ntPy5cv1+TJk8MdBUCUYD8CAITJoUOH9O9//1tZWVnhjgIgitD8AUAY5Obm6tZbb1WTJk3Us2fPcMcBEEXY7QsAABBF2PIHAAAQRWj+AAAAogjNHwAAQBSh+QMAAIgiNH8AAABRhOYPAAAgivx/TofG46Bnm6IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1a5364a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Previously known data on that time r and event c\n",
    "feature_size = x_test.shape[1]\n",
    "idx = r*feature_size + c\n",
    "print(\"\\n\\n# Previously known data on that time r = {} and event c= {}\".format(r,c))\n",
    "print(\"# Look column {} of the dataframe\".format(idx))\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches(10,5)\n",
    "\n",
    "model.df[idx].plot()\n",
    "plt.scatter(model.df.index, model.df[idx])\n",
    "for i in range(4):\n",
    "    plt.axhline(y=model.means[idx] - i * model.stds[idx], linewidth=4-i, color='b')\n",
    "for i in range(4):\n",
    "    plt.axhline(y=model.means[idx] + i * model.stds[idx], linewidth=4-i, color='b')\n",
    "plt.xlabel('Days')\n",
    "plt.ylabel('Count for Time r = {} and event c= {} '.format(r,c))\n",
    "plt.ylim(model.means[idx] - 4 * model.stds[idx], model.means[idx] + 4 * model.stds[idx])"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}