{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import datetime\n",
    "from collections import defaultdict\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import urllib\n",
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "events = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the time information from the events\n",
    "for f in [\n",
    "#    \"201508_trip_data.csv\",\n",
    "#    \"201608_trip_data.csv\",\n",
    "    \"201408_babs_open_data/201408_trip_data.csv\",\n",
    "    \"201402_babs_open_data/201402_trip_data.csv\"\n",
    "]:\n",
    "    f = open(f, 'r')\n",
    "    l = f.readline()\n",
    "    for l in f:\n",
    "        l = l.split(',')\n",
    "        tripID = l[0]\n",
    "        timeString = l[2]\n",
    "        timeUnix = time.mktime(datetime.datetime.strptime(timeString, \"%m/%d/%Y %H:%M\").timetuple())\n",
    "        events[tripID] = [timeUnix, timeString]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the earliest event\n",
    "earliest = None\n",
    "for event in events:\n",
    "  if earliest == None or events[event][0] < earliest[0]:\n",
    "    earliest = events[event]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "earliestTime = earliest[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "hourly = defaultdict(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count events by hour\n",
    "for event in events:\n",
    "    t = events[event][0]\n",
    "    hour = int(t - earliestTime) // (60*60)\n",
    "    hourly[hour] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "76832"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = open(\"hourly.json\", 'w')\n",
    "f.write(str(dict(hourly)) + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "hourly = eval(open(\"hourly.json\").read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Observations sorted by hour\n",
    "hourlySorted = []\n",
    "for h in hourly:\n",
    "    hourlySorted.append((h,hourly[h]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "hourlySorted.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [x[0] for x in hourlySorted]\n",
    "Y = [x[1] for x in hourlySorted]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm4FcWZ/78vu4osEoUJqCgDmDgqoqJxy5XgFjOSML9IdIyCM4lGjBonJpD88gCJo5iMIRpxiToEE5Utiqgoi3gFBGT3si8CIgiXfeey3Zo/qtt77rnd5/RS1VXd5/08z31un+rqt96uqq639iIhBBiGYZjSo55pBRiGYRgzsAFgGIYpUdgAMAzDlChsABiGYUoUNgAMwzAlChsAhmGYEqWoASCidkQ0lYiWEtFiIvqp4z6QiDYS0QLn74acZwYQ0WoiWk5E1+l8AYZhGCYaVGwdABG1AdBGCLGIiJoCmA+gJ4DeAPYJIf6Y5/9rAF4FcAmAdgCmAOgoeMEBwzCMVRRtAQghtgghFjnX+wEsB9DWuU0ej/QEMFIIcUwIsR7AagDd1KjLMAzDqCLUGAARtQfQBcDHjlM/IlpERC8SUXPHrS2Az3Me24Qag8EwDMNYQmAD4HT/jAXwgNMSeAZAByFEFwBbADyhR0WGYRhGBw2CeCKiBpCF/9+EEG8CgBBiW46XFwC85VxvAnB6zr12jlu+TB4TYBiGiYAQwqv7PTRBWwD/C2CZEOJJ18EZHHbpBWCJcz0ewA+IqBERnQXgnwHM8RIqhLDqb+DAgcZ1SIterBPrVAp62aiTSoq2AIjoCgD/DmAxES0EIAD8CsBtRNQFQDWA9QDudgr1ZUQ0GsAyAEcB3CtUa80wDMPEpqgBEEJ8BKC+x633CjzzGIDHYujFMAzDaIZXAudQVlZmWgVPbNSLdQoG6xQcG/WyUSeVFF0Ipi1gIu4ZYhiGCQkRQSQ8CMwwDMNkDDYADMMwJQobAIZhmBKFDQDDMEyJwgaAYRimRGEDwDAMU6KwAWAYhilR2AAwDMOUKGwAGIZhShQ2AAzDMCUKGwCGYZgShQ0AwzCMIQ4cAKqrzYXPBoBhGMYQTZsCTxg8TJcNAMMwjEHWrjUXNhsAhmGYEoUNAMMwTInCBoBhGKZEYQPAMAxjEJMHI7IBYBiGKVHYADAMw5QobAAYhmFKFDYADMMwJQobAIZhmBKFDQDDMEyJwgaAYRimRGEDwDAMU6KwAWAYhjEILwRjGIZhEocNAMMwTInCBoBhGKZEYQPAMAxTorABYBiGKVHYADAMwyTEK68Ac+ea1qIGNgAMwzAJcfvtwP33m9aihqIGgIjaEdFUIlpKRIuJ6H7HvSURTSKilUQ0kYia5zzzFBGtJqJFRNRF5wswDBOfrl2BRYtMa8EkTZAWwDEADwkhzgXwDQD9iOgcAP0BTBFCdAYwFcAAACCiGwF0EEJ0BHA3gOe0aM4wjDIWLgQ+/FBvGE2aADNm6A2DCUdRAyCE2CKEWORc7wewHEA7AD0BjHC8jXB+w/n/suP/YwDNiai1Yr0ZhkkZhw8D8+eb1sI+9u0zF3aoMQAiag+gC4DZAFoLISoBaSQAuIV8WwCf5zy2yXFjGIZh8njtNXNhNwjqkYiaAhgL4AEhxH4iyt/BIvSOFoMGDfryuqysDGVlZWFFMAzDWMOqVUDr1kDz5sX9BqW8vBzl5eXqBOYQyAAQUQPIwv9vQog3HedKImothKgkojYAtjrumwCcnvN4O8etDrkGgGGY0mTZMqBxY6BDB9OaxKdzZ+CWW4BRo2rcdu0CliwBrroqmsz8yvHgwYPjKZlD0C6g/wWwTAjxZI7beAB9nOs+AN7Mcb8DAIjoMgC73a4ihmGyQ0UFMG5cfDnnngt06xZfji3s31/796BBwNVXF35myBDgyBFtKvkSZBroFQD+HUB3IlpIRAuI6AYAjwO4lohWAugOYAgACCEmAFhHRGsAPA/gXm3aMwxjjLvvBr73PTWyqqvVyLGRIO82YACwerV+XfIp2gUkhPgIQH2f2z18nrkvjlIMw2QTItMaABdcAEyYALTVNDVlwgRv99//Xk94ceCVwAzDJEZSh58Q+YdVUSH/kuaXv0w+zGKwAWAY5ksOH7ZDhmqmTVPb+hDCjtZMXNgAMAzzJU2axKulHzggZaiiqkrOovFiy5bgcpYvl/+3bvX3c/x44ftZhA0AwzCR8DIU7kyWxYvVhNGnD3DKKXXdZ80C/umfwstrXWBPgiefLHw/i7ABYBhGOeefr0bOunXe7rt3+z/z/vvRwtq8ObjfqGHYRuCVwAzDMHFJot+8h+fcRH+OHweefz7cM9deG86/rXALgGEyyIwZwF/+YlqLurzyirf77t3AnDn+z73wAjB9enH5AwYE0yO3+2rjRqBfPxlGkhw7lmx4XrABYJgM8vDDcqGWbRQq5AcO9L/34x8DP/95cflDhtRchx3M3rMnnP+4JB2eF2wAGCalPPig3GZANY0aqZdpG0G6oq68Enj7bf26uJiYVsoGgGFSypNPAkOHqpP34IPyvw1dEyp45pl4z3/0kZq9joJy7rnJheXCBoBhSoidO01rkBwffBBfxksvFZaT1MpmXbABYJgSYf9+oFUr9XI3b/Y+1erAATnAmnbWrzetgT7YADBMiXD0qB65X/0qcOutdd0vuww4/fS67llizpzoW1/Y0HpgA8AwGSTpAcVKjxM/lizx9us3xqBa59df978XpvAdOVKuFXAZP77m+rLLatYQLFwYLKzZs4Ht24OHrxM2AAxjCfv2Ab/5jVqZ776rVp4KfvKTaFtFmNp8bdIkYMWKmt89e9a+77asunYNLvN3v7NjMzk2AAxjCTNnAo88okaWW+v89rfVyCsURlhefNF7+moxYxWny2TYsOjPxsWGgt4PNgAMUyLYXBCF4ciR8IYtv8WRlbiICxsAhkkxe/cCf/5zuGcmTQJatNCjT1D8avNEwKpVhZ89cECtLu3b13U7dCi+XKLCcngQmGGY2IRdrTprltyG4OBBdTqorFF37lxzsLqpmnqYuPHS8a9/lf+rqpSoow02AAyTQbwKpfwa50knJaNLFGw8VSwMffua1iAYbAAYhqmDe7BLUHR1Z8ydq0euSn1t6MqJChsAhskgs2bVdXNbBUG6VZ59Vq0+Uamurrl29U5bgVtovMM0bAAYhqmD7X3XSTNmTDLhTJmSTDgubAAYJuWkrUYMBNO5oiL8M7oYPNj/XqGafCGdve4lfdIYGwCGYWKjujuDSO7EqTusadPUyfLDhq4eP9gAJMDPfw7cc49pLRjbueEGvfJNFUR33WUmXEAejemFGxdJde3YChuABBg2LPyh00zpEragnjwZ+Mc/gvsvdPRiPtdfD9x5p5yWGbULZvjw8M+E7Trxw2swHEj/NFNVsAEIwObNpjVgmMIsWlTcT5QCfNIkuftlkybAyy/Hl2cL3/ue9xkGUShksG3u/gHYABRl+nS53znD2MySJfoL5LVr9covhurCNOxaBz/SbAjZABRh927TGmSXqVOzc/6sacaNA95/37QWwSlWaNpecw6DzQaCDQBjjG99y8796tNKsX7tMIXq3/8eT5divPlm9POJw7xHVEOiygDZXPgDbAAYw+Su9MwaTz8tT39SxdChwPz56uQVIvc0ryCt4ChbNvhN8wyD6gN0XGwvuFXBBqAIpZIRGPX89KfhZtwU46GHgEcfVScvbbjfYu43qeoAnVzCGu00d1exAWCYmNxzD/DYY6a1yBZeFa85c5IJ+7XXvN2j7p5qs4FgA1AE3Yn39NNAhw56wzDNzJl2fwRxef55mY46cQ8lD9IiHTWq9u9TTgFefVV/GmS9tazy/ARbKGoAiOglIqokoooct4FEtJGIFjh/N+TcG0BEq4loORFdp0vxsFRV1RwyYRPvv29+ep1uohwAztSwcyfQoEFxf24BP3Nmbfddu4CPPtJfQIdd7Jjm+fNBsf09grQAhgO43sP9j0KIrs7fewBARF8DcAuArwG4EcAzRHZEQa9ewBlnhH8u67UaRi8q8k+rVjXXbktAV1hxWLo0mXCSKFF4FpCDEGIGgF0et7yiqCeAkUKIY0KI9QBWA+gWS0NFLF8ua0IME5fp001rYAdehVuYAi/qDKkwYUQpyFUX2jYbgThjAP2IaBERvUhEzR23tgA+z/GzyXFjmMT3OtfBkSPA1Veb1sIOvvgi3vNh9i9KGpsLbZVENQDPAOgghOgCYAuAJ9SpxGSVpPc610HYgkF1d0WcgimuLr/9rTxM3sWrO8qODt8aosRXUu9gg5EJMLRUFyHEtpyfLwB4y7neBOD0nHvtHDdPBg0a9OV1WVkZysrKoqgTCBsiuxS5/HLgjjv879tWYKgmye6EKIOqTz0F3H9/8PBnzABuuim8bmGZOFGdrLRTXl6O8vJyLbKDGgBCTp8/EbURQmxxfvYC4K4bHA/gFSIaCtn1888AfGfv5hoAJpvMmgXcdpv//bQY5k2bgEsvBdas8b6/fTvwla/I/yro3bvudM44+MXzAw+EMwDf+Y6cUde4sV7jPWQIoKI+mKYBYz/yK8eDCx1PFpIg00BfBTATQCci2kBEfQH8nogqiGgRgG8C+BkACCGWARgNYBmACQDuFSItnzjD+LNsmTQCfhw5AuzYoS680aO93W1oMflNp46i2/r14fyPHRs+jLAIoTae/WTZkJZFWwBCCK/6m+8RD0KIxwAYXxe5bx9w8smmtajN0aOy9uSlV1UVUK8e0KiR/E0kC5WGDZPV0QYOHZLz3m1896BnQ/gdql5dLRcUNW0aLfysVac++8zb3e89P/wQ6N5dnz6lRiZXAm/YADRrZlqLuvTpI/V69tm69zp2lGsVcsnyRmmFOPVUoG9fffL37Im+zffZZwfz16WLt/vTT8ermEQ1AGkb2LShdlwKZNIAqDrpRzVus/7ee4HVq2vf27gR+OST5HVSxZ49wMqV4Z/79FPZdZK7y+WBA3oXFXXrBpx3Xrhn/Ao2v5O4/OJi3bpw4UYhC4WnqncIclKaV9gqW1p+smxozWXSANhK7uEnXl0ENmSIqPTrB5xzjve9Qh/zf/2XHDy9+GI9enmxdq00uIXu526HXIgLL6z9+8CB6HqpIql8lJ+uKvvOVb2D7nMN0k7JGIAkC9dRo2SBqHJQ0HYKtbqSiPsdO+SURhVceaVsIUTZfVL3pnBpriSoQncclFIcZ9IAqGwCR8kMP/gB8Mwz/jM5gpKFpnxUwsb76NFySqMK3AVOvXvXuNmSFiYXgpmWD+g3sC62pLduMmkAdHDjjbV/L1wojzQsRNyaRFZqIj/9aXC/T0RYU75rlxxXCUrQj3v9emDChPD65OKXhn/4g/ozBNz3SqrwMlVIqgo3qSmlfthgZDJpALwiPW5kv/de3d9Tp8aTmSVUZeb+/cPLCzuX3M0f+cbmww+Btnk7V40YUfsZVfTvD/zqV2plJk3aKyhes/GiEDXv2xB/mTQAXriRfddd/nO0jx2rPVALxCvYih3SDdSe6umVIaqrg8lxadJEzVmrYbEhM4fl1Vdr/54xQ25w5pXmYdIgH7/8FtYPEC+edU/RFMJMPtiypbifMIT95oulnQ01fT8yaQAKRfjw4cAJJ3jfu/JK4JprarvFydAPPVR8Glox+U88IQv1oPPWDx9O7ui8rOKVJjffHOzZ/MKISOa38eO9/R86JP/HMTAuXmfm6sC0sc8PX9X2G7nyw7xj//5ypbgXKtJVJ5k0AFH5+OO6pykVY/v2wlsEbNvmf68YixcDK1bI65Yto8tx2bQpPTOThJBTMXN3nNy/P/rpaRUVNdc6a2R/+pO3e8+e3udR/PCH4eTbtBDMS2bYcFToZdogAf7f1U9+Yod+fqTWAAgBTJ5c3I/uMLp3B9q1K+wn6odyySXF/YShXTugRw+1Ml10FDDnnQeMGVPzu1+/aOcnb9oEXHBBXfckPszcMPK7FwG1C8OSHgT2Cz9pFi6M9lwSA7QbNqiRo4vUGoAdO4DrfE4cdhMvN2NE+dh37ACGDSvsJ0gNv9DJSbn34hZI779f3I/q5nI+06bJ7o64BZubhm4XCRB9+wa/WnmxsG2kVLcHySU/fe6804weudicZwqRWgMQpLAstWmYumr3QXDj6pvflN0dv/iFWvnr19f0o7/ySrhn/+d/vN3T+tFGQWVezsIWy4VQNTsoDaTWALg0a1Z7wY4fSWWoyy/XK//WW4F58/SGoQoVM1bc/7mF/u23R5frFUY+lZVq5AcJz/QJY4XYu7fwfa+tIFRRSJbuiplXV11WSa0BcDPfvn11V9y6h01v2VKzvbIu8md9zJpV+3exDzbsBz1ypPdZqr/+dTg5qlFdMNm+MV6YhWdBsHGVbtitmovdC4PKQr7QJA1TqNq2JC6pNQCFMoi7QGvtWrkHv5f/3BWeF1wAPPqovC42L98l7gfrPh9lDGDIEDlDKBe/bo5Sx+vcWj/CpGncboL8AdugaT9pUrxww5C2LlA/oswcy8q7FyO1BqAQ7g6AO3f6+8ndfriiwrsG/fDDtX8XaxL74VWwxM1gfkcT+oWXNLZ8QA08jjyyIX5MoHsh2L59pRu3hbA5TlJrAIJEapAjhwvNqpg7t/bv5s3l/88+qz07pRBR9hZXOTOIqWHVKvn/yBHv+7rjOjfP5o9xeOXnrVu95eSfJWELL7+sTlYShWYS00ALff82fNtBD4XPLG+8Ef6Z9u3D+fdbKl5VpWfVbqGMtXGjnE7ZooX0N2MGcNVV8cNUuYWGLjp3rh03UeePRyXsB3/mmd4VjU6dkhk8tqGAUoGNedEWUtsCUIXfAdeA3oxDJNcY+NXydOrg7s65eDFw9dVqZAad4VJsYV0QGbrQEddx9A+6P5BuPZKUqYMk9IySd2wwTCVvAMKiciFOfldEUh+UW7CEGSANi9+7+C3ec7Hho0gCnVszZGE76LQYl7RT8gYgbEYL2vfv4veReE3lVEXSheiKFfH3zXfJT4+w7/Lgg2r0UIXf6mwTxjcJbCy4S6VSEYXMGYDvfMe0BsEYNKjux+IuQAoyCJzvbvLDizKOoosnnzStgT8qz8xNAlu7Thh1pNYA+GWcd97xdvfLzCozoOnFWC7FPlz3fteu+nUJS/5gso01yjjkv4/KgWjbT5/yIonv0iSF3sOGvJ1aA5BE5Pnt8e1H/gEjKvDLQPnuxT4Yv6mPOrEhg9uE6UItznTkpEl6p9Y0hxGH1BoAVZhMIF0FQq7c48fllrQ9e+oP1xS6Dwg5eFCtPF2Y3g66lElrnKfWAJg4eEI1XsZn795oC8H8/L34opxPXuxksixx6qnRnvPLI+52IjqxMX8mgcn3Nr2rqQ1pnloDEHaDp7CzS/bsCSdfFffdp1Zeoe0wsoINH1IQCg0Cq2iJuiudvWQlvegtKFnvZrKd1BqAsFsCb9wY3O+xY8B//3c4+X6ELZwOHFATzqOPBt/yQMfOhDqnuerCVJ/w+vXh5Xh1Sz3wgL//GTPChxElPmzcDdRWbHjH1BoAnScjVVaaS5y4XT4A8PzzckZS0D1jChUcphk+vLifMGmVv9f/j35UvDWpqoXh7lKbzx/+EF7Wxx/H08UWsr5Fhe2t09QaAJ3YlMGCZqAwOnvJ/P73gz8fFBXx+NFH8WXkMm1a7d8vvlizxbJXvIweHf498s+IcBk82P8ZlQVFUmcL2PSdFMJ0X7/NpMYA3Hdf4ZrqNdeEO8nnN7/xv/fSS/73dG/Clb97YJx50rnPFvM/dmxxebZi28e3fLm3u216FqPYIsS0GACT2B5HqTEAw4YBTz9d8zv/YyovDz9b4/PPvd3HjQsnRyVRM0yYwsX2TGkbumvnQmRnTrrKMYC0Gcw0khoDoAO/nTiTWlGp4sxcFbJcdu2KL8MEthk0nTN9SoksxJftRixVBkB1ZNp8UIMJWrUyrUH28cpbthcSfuj8TpLYQiGJlhdvBaEQ1REmRM3caRV4TeeLshAkSN99nEIjqRpqUhk8TFx4+U3qhC2V8W7yxCzTp3XZSFqNeFEDQEQvEVElEVXkuLUkoklEtJKIJhJR85x7TxHRaiJaRERddCmuisceq+tmU+YL2kpJawa0Aa88kDQ25blSgUidUU7jRnxAsBbAcADX57n1BzBFCNEZwFQAAwCAiG4E0EEI0RHA3QCei6vgsmVyUZOUX+Oe1FQ3nXhlmqBbK+fvv++36Mtrf5ixY7NV4GTpXWzCL17nzQv/TJSws9DSsD1vFjUAQogZAPKHB3sCGOFcj3B+u+4vO899DKA5EbWOo+CwYd7bLC9eHEdqYZKalRGHQlNVXWx/B1VcdJEaOTxfPBjuKvm0F55Jrb2wOc2jjgGcJoSoBAAhxBYAbiHfFkDu5MpNjpsSjh8HnovdpqiBB4HtYuVKs+Gr6vPWMW5jE6YPilG1Fsf0d246fEDdILDSV1m61D/RRozwdo+C6X473RkzbQVORYW3+7p1yeqRj8oP1UtWdXX65uirIolKWBTDrtLI2PwdNoj4XCURtRZCVBJRGwDujPpNAE7P8dfOcfNk0KBBX16XlZWhrKwMALBmTUStFGHTwhyVawVUyVURfhj/fgv2omDjx/j888CQIerkqXpHv/TQuQ9XLmka59NpyMrLy1FeXh5fkAdBDQA5fy7jAfQB8Ljz/80c934ARhHRZQB2u11FXuQagFxsaBr5oXsriLjhMdFR1XVj2sjoDt81AF5508bWq+m5+HHDyK0cA8DgQptKhaSoASCiVwGUAWhFRBsADAQwBMAYIroLwGcAbgEAIcQEIvo2Ea0BcABAX2WaJojtzbZSJCvpkaaB5mItAJ2FZxKtcJXy05o/ixoAIcRtPrd6+PhXfKRJbQrtl6+qy8H0LCAdYX/xhfow+vUDHnoovpxiZDE9dIXx2Wf6wzA9BsCtYnWkaiUwoHf6pwvXDIozeTLwzDPAJZeok5nVuNKBX1z16aM/HpOaBWQyPySxnsEGQ5Y6A1CIJPoNk040GzKJF716yf+26heFJMYAshBfKruAsn4msO3pbaUBSCrS0rAOIK4uWa9FhSWtC3ZsQuUgsOmuvSTGGWzOV1YagKjYVHAHIYnMlwYj5xKm3/fcc4Othg4Thm5U7j1TLBwVmJwGqrJwTtPAe9JkygCExdYBqNww9u3TH14cTBWoy5bVHOXI1Cap7ktThV4SrYwkBrptoKQNgGmCfEBxTifTOZBmS40n6XUZcUiqK0sVxeI2TXv1mMLmwh/ImAGwsclousCJci8IbsskqQyehWZ8Fgq0pEiqSyyJGVO2TCjxwmoDcPiwnGqYj6o+SNP942maX24ybNsH0vKJMpvIxjEAP5LqIs3C92FDIV8Iqw3AmDHe7m+9pTdcm8YAmOTI+vx51QvBVMiLslGbSpJKE1sxagDKy703wSqWwYYO9XZP0xS0JJgyJRuF2jvvmA0/K7Vz1bJMfT82pnlayxKjBuB3vwMGDPC/7/exHD+uR5+kWbBAr/yBA/3v2bgtgB9PPJGNMQAmOFkZW7LdMFjdBeTHsWOmNQhOoQz27rvJ6aGTNBXOpnVN4gwI3esAbKyB2xqO6e7kYlhhAP7jP2r/rnQ2kP7hD739+7UAVG4Gp+pDraoK518lUTPY8eP+ZwyrDIdRj+m0SFu3EA8CG8SNnLffru2+Y0fyuuSichViD889U4uHH4Rd+Sc1h5BT6N6Pfwy0aRNMB5U1HNMLdrLeBaS6MNJZuJlOC9sLblVY0QLIn9bpdzSgS6kkTjEOHSp8P2o8zZ9f3LjEDSOskbd1PYUXpbLBmc7vMGwlLKouPAZgEDfyt2+v2cf82DH/6Z9RWbQoml5epGmes+laVCGuu66um8kFO5t8Dy5NBhtbMrb29aucwKD7HZM6PjMqVnQBAcAVV8j/QSJs3rzi8nK58MJw/m232kGx+T22b6/rZlLfm24yF3ZaqedRetg4BpCFxZC6sKILCAg36KgKk7sdAvEHmufMiR626WmgqvbYV9U9c/So2W0BwpLETCM//25r3aYWpumCNokJJTqwxgDYRFpWIT7ySHH5aeo7t+GD8KJnT9MahEO3cd+/X438YqQp7+qWowtruoBcvv51tfJU+gfk4jVbsD1zqcTk/P0PP9QfdhrT0qYB1ChdQEkNmNvUUsrH+FYQLsePA506AZ9+mlz4fusJ1q3zf+bPf9ajiw5UFypLl9bNzLbOvgjL1q36B0+TeueshBNmjC7qYHKprwNoYFoBl5075V+SFJqjv3x5cnoEZe/ecP5VD4ytXh3+GVXo7kf1GpRWTVa2NyhEqXSdBCUJIxOHTI0BbNigTtbWrd7uJg/iaN689u9ig9VRdbVxJocXtrUibEd1F2mY+PdbVxJl8DSMnKio3FVAVRg6yJQBOP980xqEI4nmp+7DbUzOAooStkmjUchgp9H4hdGhe3d9egD2FrTV1XaklR+ZMgAqOXpUfxhJbGoXpWbitxLb65lCGdzWhUR+pGkaaBKobAF88UW4sInSl3/SCBsAHxYu9HZX+QHHXW9gS+ZNot/XxoJz4sRw/m1JLxXs21ezaWPSqBoETgJb9XJhA2AhQTNNEH82TdXLRYVepo3CzTd7u9u6wlxl+LNm6Q8jDFG7gEq91ccGIMPY3IQ2NQZgkiiFVK9e4cJI4jyAYvdM+FeNKn15L6CMYbM1D4OqD8z0IJeNYwCFDoUPG8Ybb6gL3yS6W0VR5SQxEcPWFiHABiDVFMtASbUAkpgy51WoqS7osnCSlkpMz6wxbRxKATYAiti2zbQG+lHZjFe1GVxYstJaiRJO2PDdTd+8KLRaXjdpGwS2sUXmwgZAEaedpk6WqkHgpM41CAvX7PShssC56y7/e/fe6x++CUznBR4DYKzD1q0Hpk/3do+6n0taKFQ4q1wIZrLGOXeuGjkq0zXKVvNhw58xw1+OzS0Aa/YCSguHD5vWIDiqC8cXXlAjp6ws/DNpGgPwo1B6zJ+vN+yk+Ld/MxNuobgdNy7ac2HwMzJRxlGSNBrcAkgxSdd+331Xrx5pWwgWVqdC3QH/+q/xdHGxtUUURa9hw9TIDyMnKlFadn7bjN83Q3EZAAAPqklEQVR2G9CtW3ydgsAGwELCjAGcfXa0XVRtnMURdvsN1e8QtkD3aw2anAVUKAxTtfNCRImTH/1ITdhLloTfYTcshd7v0Ue93UeOVNeVVoxYXUBEtB7AHgDVAI4KIboRUUsAowCcCWA9gFuEEHti6sn4sG6d/2wNXQXOpZfqkTttWvhn0rSJmsr0iCJr9mxvdxsNQyG8tuCIGrebN3u7f+tb0eTlY2uLzCVuC6AaQJkQ4kIhhNto6Q9gihCiM4CpAAbEDIPxwVTmyj+L2K9QbdYsnNyZM/3vhR0DEMJ7s71Dh8LplDb27w9v5F5/XY8uSRL1W6iq8nZftSqcnChdQDYYh7gGgDxk9AQwwrkeAeC7McNgfAhyaE0SO3XaOAawbx9w44113Xv3DidHJSoH9/zi6uBBO8dLwrJmTTLh6O5qsX0WUFwDIABMJKK5RPSfjltrIUQlAAghtgBQOEO+NFi7Vp0sVYXzJ5+okVMIlR/91VcDU6aEe0bVhxpmO+2o1K+vP4w0Yfq9C1W0bF4ZHnca6BVCiM1EdCqASUS0EtIo5FLgNQflXJc5f0zXruH8q6zlex37WFUFdOni/4zJc5L93jHK8ZWqDMB993m7CwFMnux9r0OHcOdh1/OputlQqHiRxJ47JvHLO8ePx18MVl5ejvLcA9QVEssACCE2O/+3EdE4AN0AVBJRayFEJRG1AeBzuCJQ2wAwOghbqHXqVPu3EMD3v1/4Gb+BNJWEfY+DB/XoEZd587zdwxZgTz0VXxdGHX75s7oaGDrU+17QNC8rK0NZzuKZwYMHh1OuAJG7gIjoRCJq6lyfBOA6AIsBjAfQx/F2J4A3Y+rIRETFAFRVFfD222r0iYPN/ahBUVlLfe01/zBM14bDoHL8aMIENbJU0rev3ekRpwXQGsAbRCQcOa8IISYR0TwAo4noLgCfAbhFgZ6MIfy6GmzA5g/LC1t3ZzVJlK46P/r2VSdLFc89Z1qDwkQ2AEKIdQDq9AwLIXYC6BFHKSYYLVoAu3frqx1fdBGgsLUZC69BaNXnxprcCiJNhXYU/BYrhh2o98P0bBu/vYAKYUOaW1y/Y4rRsKH8H2Vr3iCZb8EC4K23wsvWQf7aAwAYNUptGEuWqJWXT6HBQJPbK5skyoE3Xtx7r7kzigFg6lRzYceBDUCKcQvxjRvDPxt0QZQNtRQgmdrd55/rD0M3aesC8qNjx3D+x47Vo4dObEinzBiAJk1Ma5A8cTLQTTcF87dyZfQw0obf3HpVJPHBr1kjpx6mnbBjA6bHqtK6wjwzBsBvSXeWcQuUKLXjIKuIgWj78+ggiXnvuguRHTv0yncZMiSZcGzC9oNXbCUzBoAJhw3NzzBkYRpoUjTgUz60sX69aQ3UwgYg4/i1jNJ0sE1SmO5GUAUbS32cdZZpDdSSkSzP+JGVPvwkCrWsFJzcHZIOVO75FRU2ACkmSDdO2ENW0kaaxgCuv16vfJcoZ+AypQkbgBTjFn5+m48BwJ6MHMXDLQCGUQ8bgBQTpPb73nv69UiCE04wrQHDZA82ACkmbTN54pCEAdDdBVRK6cWkAzYATCpI4mQz7gJiSg02AClm717TGmQLvw3LVMEtAMY22AAwjMM775jWgGGShQ0AwzBMicIGgEkFSYwBMEypwQaAYRLCPb+BYWyBDQCTCrIwQ4e3aGBsgw0Ak2rS1AXEBoCxDTYADJMQbAAY22ADwKSCLHQBpam1wpQGbAAYJiG4BcDYBhsAJtWkqVbNBoCxDTYATKphA8Aw0WEDwDAJwQaAsQ02AEwq8DvDOE2Fapp0ZUoDNgBMKvA72Oa665LVIw5p6q5iSgM2AEyqWbjQtAbB4RYAYxskDFVLiEgAXCViGIbxwq9oJiIIIZSsjOEWAMMwTInCBoBhGKZEYQPAMAxTorABYBiGKVHYADAMw5QobAAYhmFKFG0GgIhuIKIVRLSKiH6pKxyGYRgmGloMABHVA/A0gOsBnAvgViI6R0dYaik3rYAP5aYV8KDctAIelJtWwINy0wp4UG5aAbRv7+VanqwSgSg3rYBWdLUAugFYLYT4TAhxFMBIAD01haWQctMK+FBuJNTGjQvdLU9IizCUm1bAg3LTCnhQnmhojzwS1Ge5Ri2iUq5U2gknKBUXG10GoC2Az3N+b3TcYtOgQXC/9QK+XevW8n/bIhp+9avFZV16qfy/Zo2/n9Gj67q1aFFzfe21xcPRwXe/C/z2tzW/9+/39te0qfxfv378MLt2jS8jKt/4hrmwozB6dG2j3KtX7fuXX56cLm5YJ55Y3O+vfy3/l5UB99wjr4N8Sw8/XNetd+9A6n1Js2bh/Kvi4ou93f/0p2DPn3yyOl0KEaI4Vc+sWcCWLUCrVsDevcDx40BVlXz5pUtlJmnWTO730rWrvD73XOD994HTTwe2bpUWVQigSRO518r+/bJ52aAB0LIlsH49cOSI3E2yqgo4ehQ4cADYt0/6b9xYZqqFC4EpU4AePeT9qiopu0kT4NAh6e+aa4B33pEF3969Ur+zzgI2bJDyAODqq4ElS4AOHYDNm4E5c4DLLgPWrZN+GjcGrrgCmDYNOPtsYPly4OBBKbuiQup3+eXABx/I+OjcGRg6VBbOF18MfPKJfN9Dh+T9I0fkM0JIQ3b4sNTv0CEZR+vWyfsXXwxUVgJ79kgddu8GFi0CbroJOOkkeeRimzbAKadIw3n++TIOP/0U2LULGDMG+PrXpb6DBgF33gkMHw4MGACccQbwla/IuFq6FOjWTerVtClw7JiMZ3dZe3U1sHgxcMklUo8LL5Tx1bEjMH26TOfVq4GdO2t0adQI2LZN/j96VD5XXi7jadw44Lzz5L25c2VaNm4sr//lX2QcNWgg9WnYUP53Ddstt8j02bVLxlOvXlJeixYyrrp0kbK2bZMFXdu2Mm+uXSvjr2FD+c7z58t33rxZ+u3SRV6ffDLw+uvA7bfLON65U4bfuLHMY26enT9f6n/RRVJm27bS37Zt8v7s2TJ/XHWVrBx88IH8Djp1knHmvuNFF8m4bNUKWLFCvoebL+vXB6ZOlYXwX/4C9Osnv58TTwSaNwc2bpT6nn8+8Le/ybg7eFCmRbduMs5atpT6NGok43ziRODKK4GxY6U+9erJb/qss6SfU06RsgFg5kyZtg0ayPho0wZYtUq+76FD0v+LLwJ33y2vGzWS380dd8j0mTRJfjfXXit1b9FCvk/HjsBpp8nv3E235s3l9ymEjLMtW+R3c+qp8n3375d5s1Ej+Vzr1jJ+tm6V+fvaa2t2n335ZeDWW6Wsxo2lriefLOOhfn2ZTvXry++qUSPp3qQJ0L07MHmydKtfX75DkyZS9jnnyDR0v4Hdu2Xc7Nkj/W3eHLzyGhctewER0WUABgkhbnB+9wcghBCP5/jhjYAYhmEioGovIF0GoD6AlQC+BWAzgDkAbhVCLFceGMMwDBMJLV1AQojjRHQfgEmQ4wwvceHPMAxjF8a2g2YYhmHMYmQlcJKLxIjoJSKqJKKKHLeWRDSJiFYS0UQiap5z7ykiWk1Ei4ioS477nY6+K4nojpg6tSOiqUS0lIgWE9H9pvUiosZE9DERLXR0Gui4tyei2U4YrxFRA8e9ERGNdHSaRURn5Mga4LgvJ6LYZ3YRUT0iWkBE423QiYjWE9EnTlzNcdyM5ilHXnMiGuO841IiutRwnurkxNEC5/8eIrrfdFwR0c+IaAkRVRDRK06+MZ2nHnC+u2TLAyFEon+QRmcNgDMBNASwCMA5GsO7EkAXABU5bo8D+IVz/UsAQ5zrGwG841xfCmC2c90SwKcAmgNo4V7H0KkNgC7OdVPI8ZJzLNDrROd/fQCznbBGAfi+4/4sgLud658AeMa57g1gpHP9dQALIbsX2ztpTTHT8GcA/g5gvPPbqE4A1gJomedmNO0cmX8F0Ne5buDINq6XI7cegC8AnG5SJwBfddKvUU5eutNknoJcLFsBoDHktzcJQIck4ilWokZ82csAvJvzuz+AX2oO80zUNgArALR2rtsAWO5cPwegd46/5QBaA/gBgGdz3J/N9adAv3EAetiiF4ATAcyDXNC3FUC9/LQD8B6AS53r+gC2eqUngHddfxF1aQdgMoAy1BiAbYZ1WgegVZ6b0bQD0AzApx7utuSp6wBMN60TpAH4DLKwbABgPIBrTeZzAP8PwAs5v/8/gIfd99cZTya6gLQtEgvBaUKISgAQQmyBjLxCuuW7b4IinYmoPWQLZTZkYhvTy+lqWQhgC2Sh+ymA3UII9zTb3LT6MmwhxHEAe4joFNU6ARgK+TEIR8dWAHYZ1kkAmEhEc4noPx03o2kH4CwA24louNPl8hciOtECvVx6A3jVuTamkxDiCwBPANjgyNkDYAHM5vMlAK5yunxOBPBtyJaS9nji3UAlfiPhSuba+kFETQGMBfCAEGK/hx6J6iWEqBZCXAhZ6+4G2S0VFOU6EdFNACqFEIvy5AcNS1f6XSGEuBjyQ+1HRFfBcNpB1ma7AhgmhOgK4ABkLdW0XiCihgBuBjDGR4fEdCKiFpDb0pwJ2Ro4CcANYUSo1kkIsQKyu2cygAmQXUvHvbyq1smEAdgE4Iyc3+0ctySpJKLWAEBEbSCbf65up3voplxnZ5BpLIC/CSHetEUvABBC7IXcBOUbAFqQ3NwvX/6XOpFc99FMCLGzgK5RuALAzUS0FsBrALoDeBJAc4M6QQix2fm/DbL7rhvMp91GAJ8LIeY5v/8BaRBM6wXIPuv5Qojtzm+TOvUAsFYIsdOp0b8Bmc9M5nMIIYYLIS4WQpQB2A05Lqg/nuL27UXo76qPmkHgRpCDwF/THGZ7AItzfj8Op/8OspbkDq58GzWDK5fBe3DFvW4RU6eXAfwxz82YXgC+AmfACMAJAKY54Y6C048I2ad4j3N9L2oGx36AuoNjjSC7JWIPAjtyv4nag8BGdIIcH2nqXJ8E4CPI/m0b8tSHADo51wMdnWzQ6zUAd1qSz7sBWAygCWTN+a8A+pnO5wBOdf6fAWAZ5JiO9niK9VHGeNkbIC3cagD9NYf1KuTsg8OQ/X59nciZ4ugwKTeSILexXgPgEwBdc9z7OPquAnBHTJ2ugGziLXIy0QInTk4xpReA8xw9FkHOSPi1434WgI8d+aMANHTcGwMY7YQ9G0D7HFkDHF2XA7hOUTrmGgBjOjlhu+m22M2/JtMuR94FAOY6+r3uFARG9YI0mNsAnJzjZlqngU4+qAAwAnI2otF8DlnhWuLkq7Kk4okXgjEMw5QoPAjMMAxTorABYBiGKVHYADAMw5QobAAYhmFKFDYADMMwJQobAIZhmBKFDQDDMEyJwgaAYRimRPk/41ErPXhVr4AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6d9af0e8d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the raw observation data\n",
    "plt.plot(X,Y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot using a sliding window\n",
    "sliding = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "wSize = 24*7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "tSum = sum([x[0] for x in hourlySorted[:wSize]])\n",
    "rSum = sum([x[1] for x in hourlySorted[:wSize]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(wSize,len(hourlySorted)-1):\n",
    "    tSum += hourlySorted[i][0] - hourlySorted[i-wSize][0]\n",
    "    rSum += hourlySorted[i][1] - hourlySorted[i-wSize][1]\n",
    "    sliding.append((tSum*1.0/wSize,rSum*1.0/wSize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [x[0] for x in sliding]\n",
    "Y = [x[1] for x in sliding]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYFNXVxt8zO8MmgoqKioJLNBpcInFBx32La9TE3URMjEgS44ZmARJNNEYSNernLhg17oqJGtdJXABRwV1cQVFAkH1gmBnmfn+cvqlbt6u6qrqrumqmz+955qnq6upbd7qr3jp17rnnkFIKgiAIQvelKu0OCIIgCMkiQi8IgtDNEaEXBEHo5ojQC4IgdHNE6AVBELo5IvSCIAjdnFBCT0SziegNIppBRK/kto0lorlE9Hru7+BkuyoIgiAUQ03I/ToBNCmllljbJyilJsTcJ0EQBCFGwrpuyGdfirEvgiAIQgKEFXoF4N9ENJ2IRhrbRxHRTCK6hYj6JtA/QRAEoUQoTAoEItpQKTWPiNYD8DSAcwDMArBIKaWI6FIAGyqlzki2u4IgCEJUQgm96wNEYwGsMH3zRLQZgMeUUjt47C/JdARBEIpAKRWLezzQdUNEjUTUK7feE8CBAN4mooHGbscAeNuvDaVUpv7Gjh2beh+6Qp+y2i/pk/SpEvoVJ2GibjYA8HDOMq8BcJdS6ikimkREw8ARObMB/CTWngmCIAixECj0SqlPAQzz2H5qIj0SBEEQYqUiZ8Y2NTWl3YU8stgnIJv9kj6FQ/oUnqz2Ky4iD8ZGPgCRSvoYgiAI3Q0igirXYKwgCILQtRGhFwRB6OaI0AuCIHRzROgFQRASYtEi4JNPgAkTgKlT0+uHDMYKgiAkxF57AS+8wOvf/S7w2GPhPyuDsYIgCF2AFSuc9TTtXRF6QRCEhKjKiMJmpBuCIAjdDxF6QRC6JdOmAV9/nXYv0mftWmD16rR7wYjQC4IQK9/5DnDuuWn3In3GjgXeeSftXjAi9IIgxMby5bxcsybdfoThssuA005Lrv0PPkiu7aiI0AuCEAurVgE75EoPUReoJn3TTcCkSeU7nkTdCILQ5Vm8GJgzh9c7O3n5wQfAG28ke9xZs/jGsmgR8OyzyR4rCvbNLk2hD1N4RBAEIZDqame9o4OXu+/OA7NJipy+uay3Hi+zMj8zS0IvFr0gCLFgCpkW+qVLkz/u559H/8xLLwGffcbry5bF2x+NeeNLGxF6QRBiQbtrAA4tNJdJMnJk9M/87nfO+qkJ1MqbMgW4+273NrHoBUHo8piiri36NGhtDd7HtLYXLoy/D2kmMPMilNAT0WwieoOIZhDRK7lt/YjoKSKaRUT/JqK+yXZVEIQssnYt+6NNcW9oSK8/++9f+P3ly4Hnn3det7fH3wevqKOuYNF3AmhSSu2olNo1t20MgGeUUlsDeA7AxUl0UBCEbNPWxsuhQ51t3/iGe59yitxLL7lfn3cesMUWwJtv8uurr3Zb/atXu91OcdBVhZ489j0SwMTc+kQAR8XVKUEQug5eFrEWf03cQhqFCROATz8F7rzT+/133gHGjYv3mF1V6BWAfxPRdCLSQx8bKKUWAIBSaj6A9ZPooCAI2cZL6O1tafnsDzgg3H5vv51sP4CuIfR7KKV2AXAogFFENAIs/iYZiV4VBKGchBH61avTEbpnnnHWC83WjaNvs2cD11zDbiGvY6U5WzjUhCml1LzcciERPQJgVwALiGgDpdQCIhoI4Cu/z48znouamprQ1NRUSp8FQUiRs84CXnkFeP11fu0n9KZ49usHjB8P/Pa3we0vXsyx8d/6Vjz91TzyCOe38cIU4SVLuL9R2XxzXu6yS3FC39zcjObm5ugHDkGg0BNRI4AqpdRKIuoJ4EAA4wFMBnA6gCsAnAbgUb82xsXtABMEITWeeop93hrbHw+w0NvumrFjg4V+2jTOfgmEt7L32gv473/zt7/1lvv1hx+6o21M9LHmzAEGDy7Nwq+pcYv6978P3HtvsNDbRvD48eOL74RFGNfNBgBeJKIZAKYCeEwp9RRY4A8golkA9gNweWy9EgQhk6xdmz8Ttb0dWN8YoTvlFGDu3OLCFp98MvpnttjC/bp/f14eeGD+voUEfO1aFvlCEAHf/KZ7229+A0yc6Ly2i4306FG4zXIQaNErpT4FMMxj+2IAARGrgiB0Jy67zLHUW1s5Xr69nfPMfJVz3g4ZwhEuxQi9OZHp2WeB/fYL/ozuz223AT/6ETBgAL+O6hMPGjDW7dk55i+91B1aah9bzylI00cvM2MFQQjNl1866z16AM89x4JeWwsccQTPCD3lFH6vvR1Yd11gs83Ct2/Org2a+AQAK1cCf/87r7e3AzNncjZLwN969wt9NIW+FNeN/dn6el5+/TXw+OM8BlBuROgFQQiNLZLXXOMI/aOPAsOHs7Ctsw777uvqnOySXrz0ElviOrFY1Nw48+fzcuBAvjEMGcKvly2LLtbmsU88sfC+hdq236vJ+U1efRU47DARekEQMo4t9ESO0Gvq6ljk7e02554L7LkncMYZwA038LaoQq9F9OKL2Vffqxew4YZs6UcRetui/8c/3O9fdJH7te2HtzN3mv9HjeUgLyaip1QkH70gCKHxyrHe1uYt9KtWAY2NwE9/yj78q65yf+6vf81vP6rQ6/1Nka6v51KGxbblxZ/+FL6djg53f+zvLI3BWbHoBUEoiZYW98BrXR2/bmkBevYErr+eUwy0tADbb8/7HHaYu43Zs3lZrDifcoqzrb6ebzRff52/v5+Vb1v0QH7OHC/eey+/XVvo/Xz25UQsekEQQuMllEcc4X5tW/SAY9XqVANPP+3+zI03sgsnaqqEtWuBrbd2qksBjkXvt78XVVX57+25Z7D7Z9ttefnJJ862IKFPI/pGLHpBEIpGqXx/tQ6RbG2NZr0uX54vtj/4QeHPrF2bX8mpkNDbIqypqoovH489WSwLpQ1F6AWhAhg8mCspJVHx6ayzgGuvdW+rq2NXjT0QqfHKZtnZmR977zXj1cRL6OvqCgu9eYxf/YqX1dXe341fygTAX8Dtm4m2+tNEhF4QKoA5c3gS09FHl9aOLW5KsUB7ie2qVd5Cv3att0gqlV8dqmfPwv2JatGb1vbEicAee/C6n0X/61/7H9sv9bIt9IMGcT/tSVXlRIReECqIxx4DJk3i9U8/5TDEKPhZ4rb7ZvlyjnH3Cq+0E55plOIslyZ+TwSaUlw3Sjn98/LRB+Hn6vGKuqmqAr797Wjtx4kIvSBUGA88wMsttuA8LVGwxU0pFkhb6AHObukl1H4CqZQTKaOjczbaqHB/ihF67brp7HQLfVQfvV+KB/3U8Oc/82s9+Bp000oSEXpBqDBMy/W++/Kt6EJ4iaGX62b4cM45E1XoOzo4rcLo0bxt440L98fraaKxkd1GJjr80rbodf+I4rfoa2o4E6e+aXndDMuFCD34JFi0KO1eCEJ5MAXtyy+B++8P/1nbZ+5n0dfU8A3ES+j9LGGl+L2aGmDkSOC664LFd82a/MieddYBli51b7v2Wh430BZ9QwNw8MGOK2rt2ugWfRihnzLFCf086CBgm22iHSMuROjBSZFeeSXtXghCMtj+cDt/vJ+bw4uGBuCkk9xte1n0tbUs9KaPXueC9xNILba1tWxh9+8fnAGzpcWJ1df07s1jBJq//Q3o2xc480xHhK+/nt1Cuv32du+8+oUIct3YN7kTTnAmWJUbEfocScS6zprlZNIThCTp7AS++ML7PdvatiNbbJEuREcHW8bmZ70sei30ptg1NXEESiGB1Ba9biPIytY3BpP6eucYZ5/N4Z8At2sfQ+/X0QE89FDhY3kd22vyk2nRZwUR+hxxC/2f/sSPaWk9qgmVxaRJLKJvvhm879y57tdVVflJvPzo6HDyvWu8kpfV1nqHV9bWhreEa2qChd7r2Hpmbm0tMGGCcyPT7Zk3h003ddqZMKHwsWw6OrxvkiL0GURXxvGLiS2W666Ltz1BKIQefAyaSQrkV4iaMoXdCmGMnRUr3BWWtN/bS+i9fPSFxFsLvW5LW+CF8BP6NWvy39NPCKZF/41vcFWr9nZgzBh27wwcWPiYGr//Y/Jk9zGyQEUL/dy5wMKF8bW3bBlw1FG8noVpz0LlEcYHbJfeu+kmXoYZjFyyhAc7Ab45AP4Wve2j19vjdN34Cf2qVWxtmy6lmhp2Wz30kFuEa2s5986SJZzi2MbvWvb7P5qb2Y0mQp8RNtnEWY9DmN95h4sv2Fx8sfjqhWiMGMFCYZet88P0FV9/feF9bR+9JozQt7Y6g5/f/CYPhtppioHCrhvzOGbEj5dFX6zQt7S4xxJ0e7/7Ha+b7eo+3ngj3xzuuMP9OfNp305e5qcbixd3UaEnoioimkFEk3Ov7yCiT3LbXieiHZLrZvLoH2zpUuDjj+NrDwAuvxy4/fbS2xQqhxdf5EHOn/0s+mdHjfJ/b5NNShN6XTWqo4Nj5desAV57jTNImvilQLDdMaax5eWjL9Z109KSv93sixlnb26vruYwSPsYms03d9YLfV9vv91FhR7AzwGY9oUCcJ5Sakel1E5KqRDDQMmw+eZs0cRhlf/oR8XnpDBLptk+/zQLAwtdl+eei6+tLbdk147fBKnHHgtuQwt9dTVHtzz/PFuvelBT4+ejt103u+zirNuum1Is+pUrCwu9GVJq7uc1uLpiBS/t2rf2Tci+2UWJZkqaUEJPRIMAHArglmI+nzS6aIGff3Lp0vwTzp5Q0dnJswQffphf64r2YVm82F1n0p7okaUfXeha2FEyXrS1AdttV3ifDz7g6katrWwU2ZOfzPj4QsfRLhFzopJXeOXChd4uHf1EceedfF3op107IqaQP19TyKK3J1KZhqAp9KY2eM1ebWnhG5l22xxzjNNfs017EDeJTKHFElao/wLgArAVb3IpEc0koquIqEB1yPh55pn86AG/E33Bgvwv3b77trc7RRGA4KnXNnb7tkX/7rsyQCsUh9e4j017u1vYtCX8yCO8POccXlZVsTAWU2oPcAu97QM30VWnWlrc2+vrHf/5ySfztpNOAg48MN+i90plYOMl9I2NPNPd3m6Ku/kkYYq7riplTqDU/7Pe78EHgb32yn/asK/5KBPRkibQi0REhwFYoJSaSURNxltjlFILcgJ/M4CLAFzq1ca4ceP+t97U1ISmpiav3SJxwAHA4YdzKFMQXndp22KfM8f9Q0WdDm2KuJ4taPLww8Czz3KlekGIQhhBtgVvwQI2VnRU2S9/6bzX0MBulSBr2Qs/i94P+ylbC70d9bLrrmzpr1zJM1sBoE8f9wxXL7xCO7fbDnj/fWDIEPd2/SRRW8s5aMz/SaOf9M1Mk+b/rNFuJfM6tw05v7EQP5qbm9Hc3BztQyEJM1ywB4AjiOhQAD0A9CaiSUqpUwFAKdVORLcDOM+vAVPoozJuHM9u0/Hud9/NWfE0YeLfw/jHL7qIixgXi+m2savAa2zrRhCU4jj23Xf33yfMOW7HbR93HPDyy447wRTlxkYW1FIt+kJcdRVP4lp3Xfd2L6EHWHznzeNUBfq9MELvZdEPHsxLO6hCW9i2IDc0OOte30l7u7fQ33abe1upFr1tBI8fPz5aAwUIdN0opS5RSm2qlNoCwA8APKeUOpWIBgIAERGAowC8XaidYhk/nic0aMaOdVeTj3KyBl0w9kl1iz0iUYBnn3XW7Tu9RgZkBZtZs7j4xdKlwIwZzvZDD3XWw8zYtAVvyhTHF9+njzs+vH9/J++MWWs1DKbQ2zlmTAYM4LGFa65xby8k9DfdxH3TNDS4reKRI9lP/sYbTuUpL6HXT/DmxC7Aact+Wt9+e+CPf+T1fffN/1+8wkc//piNThPzKQFwbjhZoJTB1LuI6A0AbwDoDx+3TRyYd1PzBNHZ7oLQohs0IGrfgS+4IFz/bNrbvYU+zTSlQvbo7OSZmQDQrx+w006OT1q7L4D8sSgv2tryo0K++opFbd993UZG375sPQPAlVdG67Mp9Jtuyla7Hpy02XhjoFcv9zY/oa+q4u/DnMBYU+M25G69FfjnP4Fhw4C992aL2kvoATYO9fiEppCF3aMHL3WcPeC4hb2eYuynBaU4/7zWo7592b2cFSJJj1LqP0qpI3Lr+ymlvqWU2kEpdapSKmDYpHhMgbTF2r47myGOGnvg1e8pwD4RirXA29u9jyEWvWDi5cPVFm3UqK/2drZMzzMcqDoIwRbVujr20Y8YARx5pPs920DZZx9n0taiRVwYxBS9U07hwcmw+Am9Dvk0r/WqqvzxrldfddanTWP3lJfQH3RQvo9+r714ef75+fsPGpS/7bvf5WVYdxUR/1+nn154LkMadAkb0xRI+wSxhT7M49Kxx3pvL3aU3Pb5+Vn0IvSCidfTqBb/5593j0XZTJnifr1wId8kdCRLTQ2XCvRKvKVnrdbWOukMNPb11NzsRKBo92SYQVg/6ut5fMC+jnVJQ6+89maf7rzTWb/pJhZ6O1Taj9NO42vV6ynmmGPy/3ci/u70d+WH7fK9/fbCRcXTINNCbz/2AcFCbzJjhnd41rRp3vtHHSXX2KIuQi+EoZDbsXfv/ElIJrvv7h7cX7KEBz6HDeO8Leedx5a4lzVqCr2N1/Wk3Rq9e3NQRKlC72XR6+vUayZt0DhclApZfmhRt6mp4b4Vsui7whyZDE3SzceOUlEq33ouJPQ77RTteMVa9CNGuF/7Cb3QNWlvZ1dK1LkVYdr1o62NJzQdeWS+AaKvAfMcMyND9t8fmD6dk+y1teULc22tdy4YwH09vfACLz/91DnG8OHB/1chGhr46cOerKUt+r593durq4NDnYNi7UuhtpYnQ/pdzyedVHhQOitk2qLX6B/66qvdj7NKeVeFKTaM0bb0OzvDTXKyH6P9wiuzNFNOCM+VV3r7cEulUEk9Helx9tlsrZssXpz/eXtQUhfInjkT+Ogj9+fDWvTap33JJU7gQyEXRhhmz+YJYBts4H3cp55ybw9j0ZvhkXFTU8O/weOPe79v5r7JMl1C6LWY25MviLyt8D594jnusmXAb39beB+vG4HfBRx1EpaQLsuWcTTMggXxtjtzJj8h+J0numJTdTVb3a+84jZCdOEPMzTST+iXLfNOT+An9H59mjrVe3JSVPTArh3Wec01bMTZ7qow+W5KmKYTiJ//f+ZMXiZ5k4mTLiX0to+byLuCU5DbJEoqAttat/EKmTTjd804YrHos8306TzYp0vyLVgQfqAvLCNGADvuCPziF/6iarpbtHvFtuoB9xOnLfR1dbxt8805FNGktpb/x7A+er09jmIaN9/MS9vdMWAAfy82Qa6b0aPz3UBJMGmS+/W3vsURQOeem/yx4yCzQj96tLMetWhvIaL6zksdQB09GrjwQl7v6OD/JejmIaTDrrsCp57qxIV/+SUvtZjaBTuK4cUXednZWVjo7Xwyfpa0Hoj0Evq2Nu+Sdq+9xoEK+royJwkVqv4Uh+tGj5uFHdANct3YE7KSwpzXoNl5567hnwcyLPR/+5uz3t7OtTBvvNG9TzEi3N4ezaIvReg/+ICX+ngdHcBddxWe7i6kjx7j2WcfXl57LS/1oORzzxUXofXhh846kb/Qt7Y675lFOLy48kq29gsJvS3Oehaodgc9+6wjZEkLvSasQM6bx24urxTK//lPPH0Jg1cEYFcis0JvYmeWLAX7RA6q7VqK0Otp519/7RxbonGyT5AhsN9+wN//Hr1dM6EdkX+0yNVXOzcbPd6krVo7k+W4cfyE6FUfta3N291y+OHsc//3v51tr73GfdLXx5lnuj/jd9MohqOOcuL9g1CKM1suWpT/nh4sLgdeFn1XossIvVcoWDGVoOzQx5/+lC/cJNB91smPOjoklj6r6EgWwDuM1+bMM506CGGxf3st5mamRMAZbAXYx7711o6Fr58qTObMyc/H8tFHXBvVT5yHD3fPGN9yS07HoIXezvNkpxAuhYcfZh93GDbdlG+QtkUdJnVzHOinBhH6MtDe7j3oadbTfPll93tnn+3d1oMPOpE6m23GF5/O1R039gVWqMakkC5nnOF+HWbgPExBkEK0tLB1a+Y+79uX4/WPP97ZNnRo4Zj7s8/ON4b23tsp3BFWnHU9Vy/3Tdyum7CMHcuRLaar5+uvgSOOKM/x9ZiCCH0ZaGsLnn22227u1zfc4L3fyJGO0OtH16OOctIg20SxwM0+7ryz81mdC0SEPrvo8RSAf6MwobClPJ0pxUKvI0YOO4z90frYthtGC73fMW0RXm89jv2P4m7RNVrtJGUnnMDtXHghu3zKic5gaRp6cYVPh0H/PvZErq5G5oX+ssv45ItzmrGONjBzZe+xh7N+333FtXvggc76YYc560cfzRajORAsvvrs0NLCFcA0YbOiRsWc86EU++i1pfrPf7KYePnVly4FPvuM17XQb7SRu21b6Hv2dKo8hbXoa2o4isce+DRvNLoCU7nQ8wE6OpwKcuVMOUDEbr1y3lySIJNCr6deA3xH93LdmOGXUdFia07aMDPw6dweQDSrzbxAzP4S8UVsDsZ2hfwYlYIZ4aVJYnLb/PnOulJsqZoTbvQAqm2FNzdz3L3ZLzMZWW0tuzO8hD6qRX/11fnbw9RuTQpt0be38zjC66+Xf5yrX7/yHi8JMin05mi6PslsC9jML2M+dofFnoGnTx6i/AkmYSkk3tr/WemuGyJOupUlvHLJhBH6Up7KlGJL1Ywnr652tvtZ4aZlOX06D0pqETYnDuki4FEt+vffd2+bOtW5BuvqgHvvDddWXJgWfW2t96QqIZhMCr1JXZ1TQEHzzDNOgqmHH+Y7fVQeeCB/23//y4/JvXsHh10G8eMfu1/X1HC/xWWTn8oibeynxbBCH9XKtYvm2EJPxOd7S4u/FW76infZxe0uND9TU8Pn2po14S16IneAA8ARQVrot922uGutFLTQpzEQ3J3ItNC/8AL/uPff7xaHd95xrGe7cEJYvEbRR4xwLCazBqzJ9dfzBeGXxW/BAr5QzdJtAP8fzz7r5MioZLL2VOPlCggj9FHdO3YhaVvoAScPjXlTuO46p56xjga65BJeNjRwVI5dxo7IPyWwH151qauqHKHXOXjKiXbdrF7ddWahZpFMC/0WWzh38WXLnO1EzglfrL8uaKabDlWz99PpC8yQOJP11/c+IfWNSRKbsb85ydSyUbGfssIOxka16M3jPPggu0Fsoa+ry084pqNhAD5/TjyR0+NqXnuNnxZtWlt5Fm+plrC2qteuLf/YUn09/x8tLSL0pZBpodfTuAFgxQpn+yGHRBdMO1VrkNBrK0hfiPrzxVo0+v/IksClSZg6qOXCFvqqqmQsepuPPsoX+sWLufaoGRNvDoZefnn+RME+fQon9oo6yckMRtCvV6/m76ncQt/QwDcZM0JJiE5o2SKiKiJ6nYgm514PJqKpRPQBEd1DRLEXMdHWDcAuE83QodFPXrt+ZFDGu7o6YOJE5wLbcksOwTSfIKK4IPSAn3nDErKBl9CHsdbjeDrzEy9T6HUmSoDdln4V0vyIatFrodc3IS30abhu9KCyCH1pRPnZfg7AiDbGFQCuUkptBWApgDM8P1UCvXrlZ67UMwZ3280d+xyVMCd/z57ui/nll935Qcz3vNLImojQu9lmm+ykbbaFfs4croIUtG9U183ee3Opv4kTnW1+rkevvDXFElWc9bG1O6ixkYU2DdeN9tGvWpX/pCGEJ9QpQESDABwKwMyAsS8AHX0+EcDRpXREp4R1da6K08ZqNtoIuOce3SfOzWFiD0gBwB138DRqE3N6eSH0I7NpuZux0KbQ20WWbbTQ65Jp3Y1Vq6JbuEnFZl9yiZNILgy2iK5e7WSutDGLXET9fxsbgUsvdYf2+h3HNGJMiz4KOlummccnDPp81zehtF03ra3svukqRT6ySNh7/V8AXABAAQAR9QewRCml7Zu5ADby+Wwgb7/trse5yy7O5KP11nPKdX35ZWHrZMoUd673UaO48rtdgWarrcL1S1e3KZS6NSw6b/ibb/IMP7uUWlfmvvv46cfML7RiBfCnPxX+XJx1Bkz++Ef3pLsgooS8/v73zvopp0SrM7xyJX9Ppl/eLkV37LG8HDnS2WZa9L17AzfdFO54Q4dyoY+oSfvs78O06NOKuunokEmGpRDo6SaiwwAsUErNJKIm862wBxlnKG1TUxOamppc75sDlC0t7AYxhTDKxWRa9b/5jfc+1uF9qa0FnniCc8h7EcWi07Nw164Fvve98BdrV0CHopqppP/7X+CiizjM9JRTeJs9ppHkbMso4yd9+3JuotdeC7f/xhs7Vajmz+fkeGFYsybf/WAL5/33Awcd5M7uaFr0m2+en+2yEOYNI4jttwfeeoufThctcvI/aYv+iy/KL7a68EhbWzyZM7NMc3Mzmr1iXGMgzFe3B4AjiOhQAD0A9AZwNYC+RFSVs+oHAfjCr4FxtkltYZ48J5zAUQWmj1L7cqOmJbZTG7/5JrDJJsFuFo3uw+TJ3u9HEao//IGtq5Urw0d1dBX0I74W17Y2fpIC2PWmhd7+/ZKy6KPQ3s75W5qa2AIeM8aZfakFzqbYAfk1a/icXL688H7mOBDgtuiTnDg0cyaXUjzkEA5n1pOjGht5Fm5nZ/mTexGxwLe2dn+ht43g8ePHx9Z24IOYUuoSpdSmSqktAPwAwHNKqZMBPA/guNxupwEoOkO0tmquvdaZGGX644ot5WaHrm2/fXiRB5wTy2/ATIv1WWcFt9XQ4EzSCqqD2dWwrdL6erePfN48Xv71r+79smDRNzezsC5aBPzjH/nx617MnQv86le8fvzx4Y+la8FGTXlbW+ukJkhS6Kuq+Aa9/vruGbBr1vBEwCFD0knXW1vLN9zuLvRJUorHbQyAXxLRBwDWBXBrsQ3pC/5nP3P8g/YswqiMHl16OJa+oPweVxct4sd2v5TINtoqq6vLTsRJHJhC7/Vb6UyL9mBakhZ91HNGR8KYT4Gdnf71fXX706eHmxuhk5jV1XGOc32jCMOAAZwG5PPP/YvwJIk2sEqYR5XNAAAgAElEQVTNv18sIvSlE0nolVL/UUodkVv/VCk1XCm1lVLq+0qpou0z84LXVXdMH30xQn/QQcX2xkELvd8A1LBh0U4+sw5od7TolSos3ldd5X6dBYte33B18RnTWt5vP+9ILrv9MGNIDzzA1aG0SI8fHz4CS1eC0gXFy53zRUe3lXsgViNCXzqZ+OpMcdCpDvr3L7698ePjKcCtLflCaRaiXHT6/6yv5xO3ra381lncjB3rtmijFM3edtvk8t5EFfof/MD9GuCEeQD//nZ75hPnihXu2gZe6Dzu+veurg6etGcydKiTqz6N5F5XXplegWzto5eom+LJRAoEU+i9hOK885wBvTD89rfx5JDWFqfXgJzGTutaCP1/6gvmxReL61eW+Ne/nHXtnsgCUYVeC7ApvtqC7ezk6komptAPHhx8HG0s2ONGYamrSzeL4/nnhxuLSoLa2soYjE2SzAn9sGH5748ZA0yaVL7+aHS/4prkdPLJPFlGW49nnskVtLoydgRKGKG3M3vGidcYTyFMdxoADBwI3Oox2nTFFd7H0QTdWPT3VOwTXCWn6/38c46YE6EvnkwI/fTpzvrMmcDdd6fXF5OBA3kZl9BPnMjT6/XF/sknwK9/HU/baWG7tQo9/Wh0FA4Qv+tGj32E9f/rm7npFjj9dCdO3kTP3vaqKxD2HClW6PWkpUoUeo0IffFkQugvvdT9utjH27jZYgsuM6cHiJMiKK46y5gDdFVVwRa9V63TOIkq9Hog1f4/7H4Cjltn773zhT4ogkjf0IoVKxF68dGXQupC7xVmmKWcFrW1ha21m28u/RhmlaC4OeccTgmQFKZFX11d2I+7/fY809gki0LvR58+/BSmqzeZBAm97k+x9RMaG/k8JEov+iUtdEx/pd7g4iD1U8YrNC1LP2hdHTB7tv/7O+1U+jGSys3+4Ydcneiaa5JpH8gX+kIpdNva+Lc1i2bEHWaqBbUU100hdG4aM/x3k02ChT5KGg8vGhu5zGWWro1yMXUqL8V1UzyZEXoz2iGtivNeBPlUo4TI+ZFUVftDDuHl/PnJTXaxXR6F0JN9brnFyaMS92+tzyed5TTs/gMGRDvOhRfyDXrx4nDZJbfaKj+BWRQaGzn0uBKFXpf3FNdN8WRG6M0JTmGTRJWDcgh9Eo/iM2e6c8vstVf8xwDyLfpCaIveTAcRt9DrMYKXXw63/5o1wMUXB8fB29TWAoMGcRivWQnNj7VrgeOOK7xPISpZ6LUlH2agX/AmM0Kvs/qddhqn8c0KQY+L2toohSQsejtsM6mJSWbfg0TIHEjUEU1JWfSnnx5u/zgmrYUpDNLRUZrroZJdN5o4rrVKJTNCr9P4Zu3xLMi3GkeSJ9Oij1IwoxD2zSOpGbjmceysizamqE6ezDf3pIQ+7Hm0Zk3pUV5hLPpSo2UaGngwtlKFvrU1nvGwSiUTQt+zJ6fxBbI34GJO77/gAvd7U6fGY43rNubNi+4r9sOOCkniqWHePM47DwA33lj4eJ9+ynUGtKtr3XV5EFML/Y9+5F8/IApa6MMmTCuX0Hd0lC70K1ZUrtBnJeS6q5IJod96a8d1kzWL3hTeP/2JM2w+8wy7QoYPj+cYWhS1f3nOnHja9TpGqcybx//76tXuUo72k82iRe7XOoeRGTqrSzXecgtw++38VypacMNGucThugkzGPuHPzjRI8VQX1/ZQi+URur2sy7GoMma0H/3u8DjjzsifPXV8bT7wgvAiBG8rl03hx/Oy8GDS/ep25+PS+g32gh46CF2MWnxBtwJr9rbnVjz5cs5QZ1S7qpJgCP0o0d797kYtMDfdx9X8QoqlFEuix4A/vnP4o/R0AA8+WTxnxcqm9SF/vzz3RdJ1oSeyAlTjBMz34vOjvjOO/G1b0fyxBmv/vXX+a4h06I3C7b06eOUgrMtZy30+ibqVSA+Kua5tGRJeYQ+zGDs4MEc3VMsWZpEKHQ9Uhd6+3E2a0KfFOZYRNylBd96i/Ofm8RZ6EQpt/Vtut7s31OXglu1Kt/toIU+Tg47zFkP853G5boJEvqtt+aEdsUiPmqhFFL30dvosLvujil6RPmCZ7pFojJrlvv13/4W3yAvwCJvuoImTHD+n+23z9+/owN46qlwQh8UuROFMEJfzsHYUgINspL+WeiaZELo9aDe++8Dv/hFun0pF7ZFbwteKYnUTGv7xz/mIiylTsG3sePn9Wtt2dtcfnn+OIGX0I8cGV8flywJ3icuoS/0ZPLll8Czz5Ym9IsXF/9ZQciE0Gs/6tZbV05UgW3R29ZnKS4NU+gvv5z9u3FahLZFb/qoCw36Nje7X3tZwvaAbbEMHMg3uKC89OWYMKVnKJci9GaEkyBEJVDoiaieiKYR0QwieouIxua2305En+S2v05EO0Q9uPYbl1rEuytiiouXRRhX4ex+/XigNI5UyAsWOOu20K+zTvT2GhrynzTi+r+18RAm2VjSrptSUxQDPPaw//7Ff16obAKFXim1BsA+SqkdAQwDcAgR6Qjy85VSOyqldlJKvRn14Poi33jjqJ/s+vTs6WR67N07X+jjsugBdqfo7/r884t34+jxE6XcUT21tfw0FpSL5K673K8bGtwT0oDibvqzZrlDD0891bHku4vQA5wPauedS2tDqExCuW6UUvpyrAdH6ugH4pKis7Xg/N//ldJK10WnfaipyRf2OCpPXX89L7VrYe1a4KqrgI8+Kq1d+0ai3VBBIYB24rD581mUAXZd3XFH+Hwmra38VDF5Mvv1dQjssGE8zqOFPuiGGeeEqfZ27xQWcQn9+ecDr75aWhtCZRJK6ImoiohmAJgP4GmllC7+dykRzSSiq4gosnd9zRpOV1uJrhuAw+1OP50F2PbRP/xw8e1qYTngAF5qi1MP8Jaaslgpd3/DhsTas2efe46XPXtyG1HGEr76ipdHHul2I3V0cFv6Oyhkabe0cC3SuMIrzzoLGDIk//24hF4QiiXUqaeU6gSwIxH1AfAwEW0LYIxSakFO4G8GcBGAS70+P27cuP+tNzU1oampCUA8j81dmepq4IwzgIsuijeeXI99aAHW0S1a6EuNwLH7GjTr9oc/5PQG5uxZ3T8zvl8XwA6Dmcr6hRecdR3G+Mkn3n01WbiQl4MHhzumH7W17IK67Tbv90XohTA0Nzej2Y5YiIlIp55SajkRNQM4WCk1IbetnYhuB3Ce3+dMoTcp9QLrDtTWsjjFKfS6LX0Tra5mn7qOzS9V6G0rOSif/h//6C30Rx/NE7t0e2GFvtB3pYub+PXVZNkyjvsv9Ymyrq5w1lF9M6u0EoBCNEwjGADGjx8fW9thom4GEFHf3HoPAAcAeJ+IBua2EYCjALwdW68qiJoa4JVX8gukh+Hee4EddwSuvNIdM67FzSx319DgJBorNbKlrc3dRpCAaYG3i7TcfTcvzRvTZ58FH7/QHAPb517oprB0aXHRQjZ1dYUHorNUMU2oTMLYGBsCeJ6IZgKYBuDfSqnHAdxFRG8AeANAf/i4bYTC6IHMRx+N/tnRo7mS1IUXAt/+NpfPW7KExW7UKLfvvFcvJytmHBa9KV5BrhttMduTqWpr3e6MWbOA994LPn4UoS90U0tK6L/8kr8T3c+4QkYFoVgCXTdKqbcA5KX8V0rtl0iPKgwtdObgZmMj1xiNwscfAyeeyKK/wQb5A4wLFgAnn8zrcQi9acUH1UIl4oFXrwRjdXXO/z50aLjjx2HRP/ggPz0EJT0Lg/bRa3RU0+23A+ec4/QhqdrAghBE6l7DX/0q7R6ki9dM4MmTwyUh80rrq1RwyGAcrhstXpdeGm5AfZ99vLeb//+IEeEiYOKw6I89loulJGHRr1zJS+2yssdMBKHcpCr0vXuzBVrJ2JEYPXpwCuMwfl0vy9xP6OfNK/y5KHz0kVOTNijFQBBmP+vr+QYX9L+HEfpRo/h1oSIus2bFJ/SrVnG4bJ8+TjSPflJpb+cbXSnZKwWhFFIV+tZWsXJsoe/sdCJxglixIn+bUvnFXAB3VlAdvx6Vvn05bn3yZHd/S8GMViFiKzgooZvf+0o5tVl1acoTT8zfz4zVf+ONaP31oq6OQzw/+4xv1EuX8nZT6IPcW4KQJKkJ/RNP5IfCVSK266azk8W/2Pz0hVw31dVcYPnDD4tru6EB2HZb97Zddy2uLY19o+jZ0/sGZuIn9FOnOpk0C82wNaN//vrXcP0shGmsNDQ4Qq9nsZZaGFwQSiU1oT/0UF5W+gCVbdHfcYd3SoQo+Al9Wxswdiyw3XbFtbtqldsXPWZM6dW3PvsMeOwx5/XQocG1Vf2EfvfdwxkO5s1liy2C9w9Cl4QEWPR1qOutt/JShF5Im9QHYysdWwBOPDG868YLbdF7ucSqqqLNPrVZs8ZdECWO+PBNNuG6vJq99wbefbfwZ556yv898yawX0BcmM6zUyr9+/NM3Ftv5e/Xriu8erV/nn5BKAepCP3Ysbzca680jp4tvKbFl+K6AQpH3RQr9PoGYk7siruYCcChpX6Tj6ZNY//6yy+Hayso5cDEidH6VojNNwd+9CPHbaMzsra08P8jNV+FNCm70H/xBfC73/H6FVeU++jZw0uQS3HddHZyNSM/oa+rK06gdQ6ZjTZyfNxJTATq0SM/dbHmO98Brr02fPSKnjBmumr0elLC+/nnvNTfzWuv8fctFr2QJmUXerN4c6VmrTSprc0veVfIddPWBkya5Lx+6CHgZz9zXs+eXTiksFiL3nxK0PH7p58evZ0gCln0AFv0Yevp6n6a36X+35NOMDZsGC9bWyW6TEifsgu9jjEGROg166wD7LKLIwaFStPNmgWcdprjiz7ySPYJ62gYnUrBTwyLFXozZFML5267RW8niEIWPcBx9lGF3nw60jeRCy4orn9h+de/eIygrY2/O3HdCGlS1sSpLS2cB0RjJ7mqZJ580hEkbdHbtVkBJwZ89mzeT6cisDNDal+xTRwW/T33JCdcQRZ9Z6f//2ajZxfbQr/hhsBvf1t8H8NQW8shnmvWiEUvpE9Zhf48K5GxWPQO/fs760RODnnb165FcOFCt3gMHcpZMDV+8e11dU6x6iiYtQOOOSb658MSxqIPW//24IOBp592u26SjoD56CPgxRd5XRckqfS6C0L6lM11s3w55xYxsSsOCQ5+dUi1Rb/PPm7xuPFGQKf9X3994MADvdvVT1FRq0z9+MfOQGOSBFn0f/hD+IikX/6SE7yZFn1ra7JulCFD2LUGOL9h0scUhCDKJvR2YYa77pJCDIXwi44xRdD8Tnv1cmrQ7rCDf7sDBvAyqvumUOx6nARZ9JqTTuLJWsuWuZ8wzKcawHky0pQzpl3/hg89JBOmhHQpm+vGvrhKnTrf3fGz6E2h33JL93tavB9/vHDb22yT3RzpQRa95u9/d9YffJDdXTfeyHn5TdIU+vp6Z0wqKH+PICRJ2Wxq++ISn2Vh6uu9xXjZMmDkSF4/9lj3e7qiVJD16Nd2FjAteqU4vcD8+cFPFC+95HwvJrbQl9ONUlfnVPUKW/RcEJIgNYtehL4wn30GHHUUMGOGe/vy5ZxF8tpr+X2TE04ADj88uO1iJk3tuy8waFC0zxSDadG3tvLA5oYbOu//5S+Oi8pk992927OFPihXf5zU1XHBFwBYd93yHFMQvCib0NtWpswUDGbmzPxtK1bwIPY55+S/RxRugNvPLVSIxsb8J4gkMC16r5vRL34RrT07nUQ5M6bW1QGLF7Or7Mgjy3NMQfCibEJvx4OL0BdHa6s7FLMYihH6trbyDCiaFn0cSdP0//rnP/PT0HbblW9gVAv9JptIllYhXQJ99ERUT0TTiGgGEb1FRGNz2wcT0VQi+oCI7iGiSDeNpKegd1fi8DEXM2mqXC6P2lqnylScQn/ZZcDvf19e101zMzB9ev5kNkEoN4FCr5RaA2AfpdSOAIYBOISIhgO4AsBVSqmtACwFcEaYAx53HDB8eAk9rnDimE5fjEVfrpzqRDz7dfr0eIVez5KdObN8Qt/UxEsReiFtQkXdKKV0ZHM92N2jAOwD4MHc9okAjg7T1gEHBBeWEIAzz8wPFQTimU5fzGBsOS3ho47iLKdxCL2OMNI3tgkTyue6+da3eCkzwIW0CeVAIaIqAK8BGALgOgAfA1iqlNIJYOcC2CionTfeKL66UaVx/PFczMImDqEvJryynFWSGhv5RhSXRb9mjfvGVq4bln7yEsNGSJtQQp8T9B2JqA+AhwFsE+Ug4/TcfACLFzehST/TCr74ZbB89VXg5JNLa7vYwdhyCaQeQyil+IrZlv2/luv/0Omn77+/PMcTujbNzc1obm5OpG1SOpdr2A8Q/QbAagAXAhiolOokou8AGKuUyqsgSkQq6jEEYMoUztUyZYp7OxFXWColRTARC2CUSTxbbMGTloYOLf64YfnpT4Htt+dCIzvv7H4v6ql0/PGcutkU+y23BD74oPR+BvHaa8Cee4ab6SsINkQEpVQs8Vphom4GEFHf3HoPAAcAeBfA8wCOy+12GoBH4+iQwPhZ9OuuC2y1VWltf+97wNGhRlQcli/ntLvlQLtb4nDdfP55/vcYx5NCGHbeWUReyAZhXDcbApiY89NXAbhXKfU4Eb0H4B9E9HsAMwDcmmA/K466Om+hW7Wq9DkIu+8eLXvl2LGcQK1cg4ra3RKH0Hv9n08+WXq7gtCVCBR6pdRbAHby2P4pAAmUTAjTot90U2DMGOCss+IJr2xoiOa20TV+y5kMLC6L3kvoZbKeUGlIouCMYlr0n3/O6Xd1xE2p6Z2jCr1GF9tOGi30cSRe+8MfnPU99uBluQZjBSEriNBnFNtHX1UVX4rdYoW+XGgfvZmKGABOPDF6WxdfzE9EALD11ryUhHpCpSFCn1FsH/2ddwIvvBC+MHYhsi702kdvCv1xx3GxmmLQ6TZuuIGXYtELlYZknMkotkXf0RE9UsaPriD05gSnffcF7ruv+Pa00NfVAZdcImX9hMpDLPqMUsykpihtP/FEMm3HgRb6/ffn2P1nny2tPbOY+GWXSQlLofIQiz6j2AUz4mTHHXkZdrbr/vsDBx+cTF+8MNMWxOFmmT8fGDiw9HYEoasiQp9RtOsmiUnF/fsD66zDdUzDCOnatcCwYfH3ww8zEVkcA6e/+Y13gjhBqBRE6DNKdTW7GHR6XU1c5fwaG3nyVb9+wft2dJQvoRngDq+Mw6LX8wAEoVIRb2WG8fLTxxUx0tjIFn0Yypm5EnCEPi7XjSBUOiL0Gaa21qmfqrnllnja7tkzv20vVq7kNLvlFHrTRy8x74JQOiL0GWbFivxi2PvsE0/b2nUTxLnn8rJcs2IBx0cf1wQxQah0ROgzTrGThIII67q57TZeJhXq6YV23cSRwE0QBBH6TLPllrw87LD42w5r0euon29+M/4++KGFfvVqKcMnCHEgQp9hbs0lfq5JIDaqVy/3RCI/tND37Bl/H/yoq+OZu21tMotVEOJAhD7DaHGtqQEOPzzetvv2DSf0aVBfz32rr+dqWIIglIYIfYbp1YuXdXXATTdxYq+4qK8PN/O2nNE2mvp6LnSS5Xw8gtCVEKHPMFroGxt5Cn8pib1swubS2WGH+I4ZlnKVLBSESkGEPsNo100SA5I6Vj2IbbcF7rgj/uMXokcPLg4uCEI8iNBnGC30P/95/G03NoYrXF3uWbGaIUPKf0xB6K4ECj0RDSKi54joHSJ6i4hG57aPJaK5RPR67q+M+Q0rAx1ts/HG8bcdNrwyLaFPItJIECqVMBZ9B4BfKqW2A7AbgHOIaJvcexOUUjvl/p5MrJcVjFLJhBh+9RVw9dXB+6Ul9OWciSsI3Z1Au0kpNR/A/Nz6SiJ6D4C2MSX4rYsSxpoHxKIXhO5AJB89EQ0GMAzAtNymUUQ0k4huIaK+MfdNSJDzzw+3X0tLeSdLadZbr/zHFITuSmi7iYh6AXgAwM9zlv31AH6nlFJEdCmACQDO8PrsuHHj/rfe1NSEpqamUvosxMCAAeFSAK9YAfTunXx/bC67DBg1qvzHFYS0aG5uRnNzcyJtkwpRwoiIagD8E8ATSqk8zy4RbQbgMaVUXtQ1EakwxxDKi1LsB29vz/eHz58PbLgh70METJsG7LprOv0UhEqFiKCUisU9HtZ1cxuAd02RJyKzCucxAN6Oo0NCeSBykofZLFrkfi0x7YLQtQl03RDRHgBOAvAWEc0AoABcAuBEIhoGoBPAbAA/SbCfQgL06FE4Q+R11/FSEosJQtcmTNTNSwC8gt0knLKLo4XeRueYOeccdutIYjFB6NrIzNgKplcvHmy1MZOJpRFaKQhCvIjQVzD9+3OWSBtT6GXikiB0fUToK5jevbn4t40p9GHKDQqCkG1E6CsYv6gbU+jXX798/REEIRlE6CuYhgbv4h4tLcCee/K6TFoShK6PZBSpYBoavC360aOdQVpd/EQQhK6LWPQVTH29t0V/wgnAhRfyuuScEYSuj1j0FYyf6wYANt8cePddYOjQ8vZJEIT4EYu+gvFz3dx0E98AvvENiaMXhO6ACH0FU18PzJ7t3rZ2LS8XLix7dwRBSAgR+grmk0+A6693b9MWfr9+5e+PIAjJECpNcUkHkDTFmeWjj4CttgI6O51tX3/Nuerb2sRtIwhpkkaaYqEbst56+UVFWluBjTYSkReE7oQIfQXTqxenQDAt+tZWSUssCN0NEfoKprqac9Gb+W5E6AWh+yFCX+H06QMsX+68FqEXhO6HCH2FYwv9F19IaKUgdDdkZmyF07s357U5+miuDXvjjcBXX6XdK0EQ4iRMzdhBACYB2ABcH/ZmpdQ1RNQPwL0ANgPXjD1eKbUswb4KCaBTFT/yCPDUU8DIkcCcOWn3ShCEOAnjuukA8Eul1HYAdgMwioi2ATAGwDNKqa0BPAfg4uS6KSSFmZN+1Sqgrg7YY490+yQIQrwECr1Sar5SamZufSWA9wAMAnAkgIm53SYCOCqpTgrJoYW+sZFfL1vG2wRB6D5EGowlosEAhgGYCmADpdQCgG8GAKQWURekvp7DK7VVf/PNwJdfptsnQRDiJbTQE1EvAA8A+HnOsrfzGkiegy5IfT0wf747t03Pnun1RxCE+AkVdUNENWCRv1Mp9Whu8wIi2kAptYCIBgLwjdUYN27c/9abmprQ1NRUdIeFeKmvB+bNA9ZdF7j0UuCss4Bf/zrtXglC5dHc3Izm5uZE2g6V1IyIJgFYpJT6pbHtCgCLlVJXENFFAPoppcZ4fFaSmmWYn/wEWL0a+PBD4OWXgY4OyXMjCFmgrEnNiGgPACcB2JeIZhDR60R0MIArABxARLMA7Afg8jg6JJSXddYB7ryT890QicgLQnck0HWjlHoJQLXP2/vH2x2h3Oy4Iy/nzk23H4IgJIfko69wlAKqqpx1QRCyQZyuG0mBUOEQcdRNe3vaPREEISnEohcEQcggUmFKEARBCI0IvSAIQjdHhF4QBKGbI0IvCILQzRGhFwRB6OaI0AuCIHRzROgFQRC6OSL0giAI3RwRekEQhG6OCL0gCEI3R4ReEAShmyNCLwiC0M0RoRcEQejmiNALgiB0c0ToBUEQujki9IIgCN2cMMXBbyWiBUT0prFtLBHNzRUK18XCBUEQhAwSxqK/HcBBHtsnKKV2yv09GXO/EqW5uTntLuSRxT4B2eyX9Ckc0qfwZLVfcREo9EqpFwEs8XgrlhJXaZDFHzWLfQKy2S/pUzikT+HJar/iohQf/SgimklEtxBR39h6JAiCIMRKsUJ/PYAhSqlhAOYDmBBflwRBEIQ4IaVU8E5EmwF4TCm1Q5T3cu8HH0AQBEHIQykVi4u8JuR+BMMnT0QDlVLzcy+PAfC23wfj6qggCIJQHIFCT0R3A2gC0J+IPgMwFsA+RDQMQCeA2QB+kmAfBUEQhBII5boRBEEQui6JzYwlooOJ6H0i+oCILkrqOMbxvCZ29SOip4hoFhH924wOIqJriOjDXOTQMGP7abk+zyKiU0vs0yAieo6I3iGit4joZ2n3i4jqiWgaEc3I9WlsbvtgIpqaO8Y9RFST215HRP/I9WkKEW1qtHVxbvt7RHRgsX0y2qvKTcCbnIU+EdFsInoj9129ktuW9jnVl4juz/1/7xDR8Az0aavcd/R6brmMiH6WgX6dS0RvE9GbRHRX7rxJ+5z6ee66K68eKKVi/wPfQD4CsBmAWgAzAWyTxLGMY+4JYBiAN41tVwC4MLd+EYDLc+uHAPhXbn04gKm59X4APgbQF8A6er2EPg0EMCy33gvALADbZKBfjbllNYCpuWPdC+C43PYbAPwkt/5TANfn1r8P4B+59W0BzAC7/wbnfm8q8Tc8F8DfAUzOvU61TwA+AdDP2pb2b3cHgB/m1mty7abaJ6t/VQC+BLBJmv0CsFHu96szzqXT0jynAGwH4E0A9eBr7ykAQ8rxPZX8w/r8Q98B8ITxegyAi5I4lnXczeAW+vcBbJBbHwjgvdz6/wH4vrHfewA2APADADcY228w94uhf48A2D8r/QLQCOBVALsC+ApAlf37AXgSwPDcejWAr7x+UwBP6P2K7MsgAE+Dx4O00C9MuU+fAuhvbUvttwPQB8DHHtszcT7l2joQwAtp9wss9HPAolgDYDKAA9I8zwEcC+Bm4/WvAVyg//8kv6ekXDcbA/jceD03t63crK+UWgAAiqOENsht9+ufvf0LxNRvIhoMfuKYCv5RU+tXzkUyAzwH4mmwRbBUKdVpHdfVJ6XUWgDLiGjduPsE4C/gk17l+tgfwJKU+6QA/JuIphPRyNy2NH+7zQEsIqLbc26Sm4ioMeU+2XwfwN259dT6pZT6EsBVAD7LtbMMwOtI9zx/G8CInKumEcCh4CefxL+nSste6TfynGgIKBH1AvAAgJ8rpVZ69KOs/VJKdSqldgRb0buC3UlhiYZysssAAAK2SURBVL1PRHQYgAVKqZlW+2GPldTvt4dSahfwBTmKiEYg3d+uBsBOAK5TSu0EoAVscaZ6Pv2vcaJaAEcAuN+nH2XrFxGtA+BI8FP+RgB6AoiSfDH2Piml3ge7aZ4G8DjYJbTWa9e4+5SU0H8BYFPj9aDctnKzgIg2ADj2H/zYhlxfNjH20/2Lvd+5wZ4HANyplHo0K/0CAKXUcgDNAHYDsA4R6fPBbP9/fSKiagB9lFKLC/S1GPYAcAQRfQLgHgD7ArgaQN8U+wSl1LzcciHY7bYr0v3t5gL4XCn1au71g2Dhz8T5BPYpv6aUWpR7nWa/9gfwiVJqcc5Cfxh8nqV5nkMpdbtSahelVBOApeBxu+S/pzj8ch6+qGo4g7F14MHYbyRxLOu4gwG8Zby+Ajn/Gtjy0YMch8IZ5PgOvAc59Po6JfZpEjjTJ7LQLwADkBu4AdADwH9zx70XOT8f2Od3Vm79bDiDVD9A/iBVHdilUPJgbK7dveEejE2lT+Dxi1659Z4AXgL7n1M9pwD8B8BWufWxuf6kfp7n2r0HwGkZOc93BfAWgAawJXwHgFFpn+cA1sstNwXwLnjcJfHvqaQfNuAfOhh8t/oQwJikjmMc727waP8asF/uh7kv4ZlcP54yvwwAf8v9aG8A2MnYfnquzx8AOLXEPu0BfjSbmTtZXs99L+um1S8A2+f6MRMcAfCr3PbNAUzLtX8vgNrc9noA9+WOPRXAYKOti3N9fQ/AgTH9jqbQp9an3LH17/aWPofT/O1ybX0LwPRc3x7KXeyp9inXXiN48Ly3sS3t72ps7jx4E8BEcARgquc52LB6O3deNZXre5IJU4IgCN2cShuMFQRBqDhE6AVBELo5IvSCIAjdHBF6QRCEbo4IvSAIQjdHhF4QBKGbI0IvCILQzRGhFwRB6Ob8P4TRv2Lj+2SzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6d9623d630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X,Y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Autoregressive features\n",
    "def feature(hour):\n",
    "    previousHours = []\n",
    "    for i in [1,2,3,4,5,24,24*7,24*7*365]:\n",
    "        previousHour = hour - i\n",
    "        previousHourExists = previousHour in hourly\n",
    "        if previousHourExists:\n",
    "            # Use the feature if it doesn't exist\n",
    "            previousHours += [hourly[previousHour]]\n",
    "        else:\n",
    "            # Otherwise add a \"missing value\" indicator\n",
    "            previousHours += [0]\n",
    "    return previousHours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feature(x) for x in hourly]\n",
    "y = [hourly[x] for x in hourly]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "theta,residuals,rank,s = numpy.linalg.lstsq(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.47105067, -0.28432565,  0.10555845,  0.01445134, -0.02131945,\n",
       "        0.17501989,  0.53969635,  0.        ])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}