10 - Final Solution - Booking Demand Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Booking Demand Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is simple: Predict the booking demand for five 15-minute timesteps into the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment and run to install the required libraries\n",
    "# !pip3 install pygeohash simdkalman sklearn scipy==1.2 statsmodels pystan fbprophet tensorflow==2.0.0-alpha0 --upgrade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import datetime\n",
    "import multiprocessing\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pygeohash as gh\n",
    "\n",
    "import tensorflow.keras as tfk\n",
    "\n",
    "import simdkalman\n",
    "from sklearn.linear_model import HuberRegressor\n",
    "\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.tsa.holtwinters import SimpleExpSmoothing\n",
    "\n",
    "from fbprophet import Prophet\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment and run to download training.csv from the cloud\n",
    "# !wget \"https://s3-ap-southeast-1.amazonaws.com/grab-aiforsea-dataset/traffic-management.zip\"\n",
    "# !unzip \"traffic-management.zip\"\n",
    "# !mv \"Traffic Management/training.csv\" \"training.csv\"\n",
    "# !rm -r \"Traffic Management\" \"traffic-management.zip\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_location = 'training.csv'\n",
    "data_location = 'test.csv'    # <- for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>demand</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp03pn</td>\n",
       "      <td>10</td>\n",
       "      <td>14:30</td>\n",
       "      <td>0.024721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp09sw</td>\n",
       "      <td>9</td>\n",
       "      <td>6:15</td>\n",
       "      <td>0.102821</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03tu</td>\n",
       "      <td>1</td>\n",
       "      <td>12:15</td>\n",
       "      <td>0.023843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp03r2</td>\n",
       "      <td>4</td>\n",
       "      <td>22:15</td>\n",
       "      <td>0.123463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp03md</td>\n",
       "      <td>14</td>\n",
       "      <td>14:45</td>\n",
       "      <td>0.013998</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp    demand\n",
       "0   qp03pn   10     14:30  0.024721\n",
       "1   qp09sw    9      6:15  0.102821\n",
       "2   qp03tu    1     12:15  0.023843\n",
       "3   qp03r2    4     22:15  0.123463\n",
       "4   qp03md   14     14:45  0.013998"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(data_location)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Add datetime features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "timestamps = df['timestamp'].unique()\n",
    "timestamps_decoded = {}\n",
    "for timestamp in timestamps:\n",
    "    timestamps_decoded[timestamp] = list(map(int, timestamp.split(':')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 1.59 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df['hour']    = df['timestamp'].apply(lambda x : timestamps_decoded[x][0])\n",
    "df['minutes'] = df['timestamp'].apply(lambda x : timestamps_decoded[x][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 28.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df['ds'] = df['day'].astype(str) + \" \" + df['timestamp']\n",
    "df['ds'] = df['ds'].apply(lambda x : datetime.datetime.strptime(x, \"%j %H:%M\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>demand</th>\n",
       "      <th>hour</th>\n",
       "      <th>minutes</th>\n",
       "      <th>ds</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp03pn</td>\n",
       "      <td>10</td>\n",
       "      <td>14:30</td>\n",
       "      <td>0.024721</td>\n",
       "      <td>14</td>\n",
       "      <td>30</td>\n",
       "      <td>1900-01-10 14:30:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp09sw</td>\n",
       "      <td>9</td>\n",
       "      <td>6:15</td>\n",
       "      <td>0.102821</td>\n",
       "      <td>6</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-09 06:15:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03tu</td>\n",
       "      <td>1</td>\n",
       "      <td>12:15</td>\n",
       "      <td>0.023843</td>\n",
       "      <td>12</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-01 12:15:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp03r2</td>\n",
       "      <td>4</td>\n",
       "      <td>22:15</td>\n",
       "      <td>0.123463</td>\n",
       "      <td>22</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-04 22:15:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp03md</td>\n",
       "      <td>14</td>\n",
       "      <td>14:45</td>\n",
       "      <td>0.013998</td>\n",
       "      <td>14</td>\n",
       "      <td>45</td>\n",
       "      <td>1900-01-14 14:45:00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp    demand  hour  minutes                  ds\n",
       "0   qp03pn   10     14:30  0.024721    14       30 1900-01-10 14:30:00\n",
       "1   qp09sw    9      6:15  0.102821     6       15 1900-01-09 06:15:00\n",
       "2   qp03tu    1     12:15  0.023843    12       15 1900-01-01 12:15:00\n",
       "3   qp03r2    4     22:15  0.123463    22       15 1900-01-04 22:15:00\n",
       "4   qp03md   14     14:45  0.013998    14       45 1900-01-14 14:45:00"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get unique locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The dataset has 1300 unique geohashes:\n",
      " ['qp03pn' 'qp09sw' 'qp03tu' ... 'qp09j9' 'qp0917' 'qp09dz']\n"
     ]
    }
   ],
   "source": [
    "locations = df['geohash6'].unique()\n",
    "print(\"The dataset has\", len(locations), \"unique geohashes:\\n\", locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fill in the missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1900-01-01 00:00:00 --to--> 1900-01-14 23:45:00 with steps: 0:15:00\n"
     ]
    }
   ],
   "source": [
    "min_ds = df['ds'].min()\n",
    "max_ds = df['ds'].max()\n",
    "delta_ds = datetime.timedelta(minutes=15)\n",
    "print(min_ds, '--to-->', max_ds, 'with steps:', delta_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 7.95 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "present_values = df[['geohash6', 'ds']].values\n",
    "present_values = set(tuple(x) for x in present_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 26.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "missing_values = []\n",
    "for location in locations:\n",
    "    cur_time = min_ds\n",
    "    while cur_time <= max_ds:\n",
    "        if not (location, cur_time) in present_values:\n",
    "            x = {'geohash6' : location,\n",
    "                 'day' : cur_time.dayofyear,\n",
    "                 'timestamp' : str(int(cur_time.hour)) + \":\" + str(int(cur_time.minute)),\n",
    "                 'demand' : 0,\n",
    "                 'hour' : cur_time.hour,\n",
    "                 'minutes' : cur_time.minute,\n",
    "                 'ds' : cur_time}\n",
    "            missing_values.append(x)\n",
    "        cur_time += delta_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2.93 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df_missing = pd.DataFrame(missing_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>demand</th>\n",
       "      <th>hour</th>\n",
       "      <th>minutes</th>\n",
       "      <th>ds</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp03pn</td>\n",
       "      <td>10</td>\n",
       "      <td>14:30</td>\n",
       "      <td>0.024721</td>\n",
       "      <td>14</td>\n",
       "      <td>30</td>\n",
       "      <td>1900-01-10 14:30:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp09sw</td>\n",
       "      <td>9</td>\n",
       "      <td>6:15</td>\n",
       "      <td>0.102821</td>\n",
       "      <td>6</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-09 06:15:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03tu</td>\n",
       "      <td>1</td>\n",
       "      <td>12:15</td>\n",
       "      <td>0.023843</td>\n",
       "      <td>12</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-01 12:15:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp03r2</td>\n",
       "      <td>4</td>\n",
       "      <td>22:15</td>\n",
       "      <td>0.123463</td>\n",
       "      <td>22</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-04 22:15:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp03md</td>\n",
       "      <td>14</td>\n",
       "      <td>14:45</td>\n",
       "      <td>0.013998</td>\n",
       "      <td>14</td>\n",
       "      <td>45</td>\n",
       "      <td>1900-01-14 14:45:00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp    demand  hour  minutes                  ds\n",
       "0   qp03pn   10     14:30  0.024721    14       30 1900-01-10 14:30:00\n",
       "1   qp09sw    9      6:15  0.102821     6       15 1900-01-09 06:15:00\n",
       "2   qp03tu    1     12:15  0.023843    12       15 1900-01-01 12:15:00\n",
       "3   qp03r2    4     22:15  0.123463    22       15 1900-01-04 22:15:00\n",
       "4   qp03md   14     14:45  0.013998    14       45 1900-01-14 14:45:00"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.concat([df, df_missing], sort=False)\n",
    "df = df.reset_index().drop('index', axis=1)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Decode the geohashes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "locations_decoded = {}\n",
    "for location in locations:\n",
    "    locations_decoded[location] = gh.decode_exactly(location)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 1.52 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df['latitude']  = df['geohash6'].apply(lambda x : locations_decoded[x][0])\n",
    "df['longitude'] = df['geohash6'].apply(lambda x : locations_decoded[x][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>demand</th>\n",
       "      <th>hour</th>\n",
       "      <th>minutes</th>\n",
       "      <th>ds</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp03pn</td>\n",
       "      <td>10</td>\n",
       "      <td>14:30</td>\n",
       "      <td>0.024721</td>\n",
       "      <td>14</td>\n",
       "      <td>30</td>\n",
       "      <td>1900-01-10 14:30:00</td>\n",
       "      <td>-5.413513</td>\n",
       "      <td>90.664673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp09sw</td>\n",
       "      <td>9</td>\n",
       "      <td>6:15</td>\n",
       "      <td>0.102821</td>\n",
       "      <td>6</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-09 06:15:00</td>\n",
       "      <td>-5.325623</td>\n",
       "      <td>90.906372</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03tu</td>\n",
       "      <td>1</td>\n",
       "      <td>12:15</td>\n",
       "      <td>0.023843</td>\n",
       "      <td>12</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-01 12:15:00</td>\n",
       "      <td>-5.336609</td>\n",
       "      <td>90.609741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp03r2</td>\n",
       "      <td>4</td>\n",
       "      <td>22:15</td>\n",
       "      <td>0.123463</td>\n",
       "      <td>22</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-04 22:15:00</td>\n",
       "      <td>-5.402527</td>\n",
       "      <td>90.675659</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp03md</td>\n",
       "      <td>14</td>\n",
       "      <td>14:45</td>\n",
       "      <td>0.013998</td>\n",
       "      <td>14</td>\n",
       "      <td>45</td>\n",
       "      <td>1900-01-14 14:45:00</td>\n",
       "      <td>-5.391541</td>\n",
       "      <td>90.598755</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp    demand  hour  minutes                  ds  \\\n",
       "0   qp03pn   10     14:30  0.024721    14       30 1900-01-10 14:30:00   \n",
       "1   qp09sw    9      6:15  0.102821     6       15 1900-01-09 06:15:00   \n",
       "2   qp03tu    1     12:15  0.023843    12       15 1900-01-01 12:15:00   \n",
       "3   qp03r2    4     22:15  0.123463    22       15 1900-01-04 22:15:00   \n",
       "4   qp03md   14     14:45  0.013998    14       45 1900-01-14 14:45:00   \n",
       "\n",
       "   latitude  longitude  \n",
       "0 -5.413513  90.664673  \n",
       "1 -5.325623  90.906372  \n",
       "2 -5.336609  90.609741  \n",
       "3 -5.402527  90.675659  \n",
       "4 -5.391541  90.598755  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rename 'demand' to 'y' (for fbprophet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>y</th>\n",
       "      <th>hour</th>\n",
       "      <th>minutes</th>\n",
       "      <th>ds</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp03pn</td>\n",
       "      <td>10</td>\n",
       "      <td>14:30</td>\n",
       "      <td>0.024721</td>\n",
       "      <td>14</td>\n",
       "      <td>30</td>\n",
       "      <td>1900-01-10 14:30:00</td>\n",
       "      <td>-5.413513</td>\n",
       "      <td>90.664673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp09sw</td>\n",
       "      <td>9</td>\n",
       "      <td>6:15</td>\n",
       "      <td>0.102821</td>\n",
       "      <td>6</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-09 06:15:00</td>\n",
       "      <td>-5.325623</td>\n",
       "      <td>90.906372</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03tu</td>\n",
       "      <td>1</td>\n",
       "      <td>12:15</td>\n",
       "      <td>0.023843</td>\n",
       "      <td>12</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-01 12:15:00</td>\n",
       "      <td>-5.336609</td>\n",
       "      <td>90.609741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp03r2</td>\n",
       "      <td>4</td>\n",
       "      <td>22:15</td>\n",
       "      <td>0.123463</td>\n",
       "      <td>22</td>\n",
       "      <td>15</td>\n",
       "      <td>1900-01-04 22:15:00</td>\n",
       "      <td>-5.402527</td>\n",
       "      <td>90.675659</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp03md</td>\n",
       "      <td>14</td>\n",
       "      <td>14:45</td>\n",
       "      <td>0.013998</td>\n",
       "      <td>14</td>\n",
       "      <td>45</td>\n",
       "      <td>1900-01-14 14:45:00</td>\n",
       "      <td>-5.391541</td>\n",
       "      <td>90.598755</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp         y  hour  minutes                  ds  \\\n",
       "0   qp03pn   10     14:30  0.024721    14       30 1900-01-10 14:30:00   \n",
       "1   qp09sw    9      6:15  0.102821     6       15 1900-01-09 06:15:00   \n",
       "2   qp03tu    1     12:15  0.023843    12       15 1900-01-01 12:15:00   \n",
       "3   qp03r2    4     22:15  0.123463    22       15 1900-01-04 22:15:00   \n",
       "4   qp03md   14     14:45  0.013998    14       45 1900-01-14 14:45:00   \n",
       "\n",
       "   latitude  longitude  \n",
       "0 -5.413513  90.664673  \n",
       "1 -5.325623  90.906372  \n",
       "2 -5.336609  90.609741  \n",
       "3 -5.402527  90.675659  \n",
       "4 -5.391541  90.598755  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = df.rename(columns={'demand' : 'y'})\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sort by time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2.09 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df = df.sort_values(by=['ds']).reset_index().drop('index', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>y</th>\n",
       "      <th>hour</th>\n",
       "      <th>minutes</th>\n",
       "      <th>ds</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp09d7</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.091915</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.342102</td>\n",
       "      <td>90.807495</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp0d5g</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.254211</td>\n",
       "      <td>90.873413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03pw</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.016879</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.413513</td>\n",
       "      <td>90.686646</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp09jr</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.408020</td>\n",
       "      <td>90.939331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp02z1</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.484924</td>\n",
       "      <td>90.664673</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp         y  hour  minutes         ds  latitude  \\\n",
       "0   qp09d7    1       0:0  0.091915     0        0 1900-01-01 -5.342102   \n",
       "1   qp0d5g    1       0:0  0.000000     0        0 1900-01-01 -5.254211   \n",
       "2   qp03pw    1       0:0  0.016879     0        0 1900-01-01 -5.413513   \n",
       "3   qp09jr    1       0:0  0.000000     0        0 1900-01-01 -5.408020   \n",
       "4   qp02z1    1       0:0  0.000000     0        0 1900-01-01 -5.484924   \n",
       "\n",
       "   longitude  \n",
       "0  90.807495  \n",
       "1  90.873413  \n",
       "2  90.686646  \n",
       "3  90.939331  \n",
       "4  90.664673  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detect temporal anomalies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get location infos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46 36\n"
     ]
    }
   ],
   "source": [
    "n_lat  = len(df['latitude'].unique())\n",
    "n_long = len(df['longitude'].unique())\n",
    "print(n_lat, n_long)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.48492431640625 90.5877685546875\n"
     ]
    }
   ],
   "source": [
    "lat_min  = df['latitude'].min()\n",
    "long_min = df['longitude'].min()\n",
    "print(lat_min, long_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "lat_sorted  = np.sort(df['latitude'].unique())\n",
    "long_sorted = np.sort(df['longitude'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0054931640625 0.010986328125\n"
     ]
    }
   ],
   "source": [
    "lat_gap  = (lat_sorted[1:] - lat_sorted[:-1]).min()\n",
    "long_gap = (long_sorted[1:] - long_sorted[:-1]).min()\n",
    "print(lat_gap, long_gap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "latitude_i  = ((df['latitude'].values  - lat_min)  / lat_gap).astype(int)\n",
    "longitude_j = ((df['longitude'].values - long_min) / long_gap).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate density images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "timesteps = 4*24*(df['day'].values-1) + 4*df['hour'].values + df['minutes'].values//15\n",
    "demands   = df['y'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 12.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "density_images = np.zeros((4*24*61, 48, 40))\n",
    "density_mask   = np.zeros((48, 40))\n",
    "for i in range(len(timesteps)):\n",
    "    t = timesteps[i]\n",
    "    d = demands[i]\n",
    "    x = latitude_i[i]\n",
    "    y = longitude_j[i]\n",
    "    density_images[t][x+1][y+2] = d\n",
    "    density_mask[x+1][y+2] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inspect the images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAEYCAYAAABiNXvJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAU+ElEQVR4nO3dbWxe5XkH8P+FMXESCIkTxzFxEgcEaAlsCbJQpUaiKgXSbloAgVSkTZmE4MsmUanTSPtl6qRJmSZV/bIPRCuqpVWtkKhGVHVCVmjUFUUUQ9MXyCCogsTU2A7BiRPzkoRrH3wqGd//g+/nOc+LL/v/k6LHz20/z7nPca4c+5/7xdwdIhLTVe3ugIjUTwUsEpgKWCQwFbBIYCpgkcBUwCKBVSpgM9trZm+Y2VtmdqBRnRKRPFbv/wObWQeANwHcA2AUwMsAHnH318tes2HDBh8YGKjreCLL1dtvv40zZ84Y+9zVFd73TgBvufsfAMDMfgxgH4DSAh4YGMDIyEiFQ4osP4ODg6Wfq/Ij9GYAp+c8Hy3aRKRFqhQwu6UnP4+b2eNmNmJmI5OTkxUOJyLzVSngUQBb5jzvB/DH+V/k7ofcfdDdB3t6eiocTkTmq1LALwO42cy2m9k1AL4O4HBjuiUiOeoOsdz9spn9A4DnAXQAeNrdX2tYz0RkQVVSaLj7zwD8rEF9EZEaaSSWSGAqYJHAKv0ILalLly4lbZ2dnW3oiSwHugOLBKYCFglMBSwSmApYJDCFWA2mwEpaSXdgkcBUwCKBqYBFAlMBiwSmAhYJTAUsEpgKWCQwFbBIYCpgkcBUwCKBaShlC7DdL86dO5e0rV27thXdkSVEd2CRwFTAIoGpgEUCUwGLBKYQq8E++eSTpO2aa65J2hRYLSz3Wi5nugOLBKYCFglMBSwSmApYJDCFWBXMzMwkbWzU1RtvvJG07dixI2nr6Oio1J+PP/44aVuxYkWl98z16aefJm1XXVXt/qDAamG6A4sEpgIWCUwFLBKYClgkMIVY87AQCgDMLGnr6upK2lhwc/vttydtV65cSdpYEFT2nmyUUqsCK4b18fLly0nbxYsX6euvv/76hvdpOdAdWCQwFbBIYCpgkcBUwCKBKcSah4VVZVhww0Iw9p5VR11FGKV09dXpX6+ysCr3usln6Q4sEpgKWCQwFbBIYCpgkcAWLGAze9rMJszs93Paus1s2MxOFo/rmttNEWFy7sA/ALB3XtsBAEfc/WYAR4rny8758+eTP2aW/InsypUryZ9mWGrXrVUWLGB3/wWAs/Oa9wEYKj4eAnB/g/slIhnq/R24193HAKB43Fj2hWb2uJmNmNnI5ORknYcTEabpIZa7H3L3QXcf7OnpafbhRJaVegt43Mz6AKB4nGhcl0QkV71DKQ8D2A/gYPH4XMN6FMiaNWtacpzFtkNB7gJ2LPCqOoS0ndj8ZjZctJVy/hvpRwCOAbjVzEbN7FHMFu49ZnYSwD3FcxFpsQX/+XD3R0o+dXeD+yIiNdJILJHAVMAigWk+sHwuFjqVLfyX89pmyA2X2M4VQP5igO0OrBjdgUUCUwGLBKYCFglMBSwS2OL7rVwSY2NjSdu2bdsaeoyyaYIsiKoy1a+WnS8YNiqNvZadT9kxIi+opzuwSGAqYJHAVMAigamARQJblCFW2TabLGyoMtqHHYeN6gH49D32etZWNQiqshBCbkDz0Ucf0dd3dnYmbdPT00kb23GBHafsvD/88MOkbdWqVUkb+z60ajeMxUh3YJHAVMAigamARQJTAYsEFirEYtO5ctddYuEUa5uZmaHHZms+sf6wNhb6rFuXv5lF2fXIkRuWdXV10XZ2LScm0jUMr7vuuqStlnW7qqzx1YxRU1HW89IdWCQwFbBIYCpgkcBUwCKBqYBFAluUKXQti4ddvHgxaWM7JrChgmwxMzZ8D+ALop09O3/TRqCvry9pY8MMa7Fy5cpKr89RNpSSXaOBgYGkjaXI7HuzevVqehz2Pc/dAYKpuovCYkycGd2BRQJTAYsEpgIWCUwFLBLYogyxapG7xScbGscWSKtlGOfatWuTNhZ2XXvttUlbLSFL7hzjKsp2J3jvvfeStv7+/qz3ZMMz2TUHeAjGriWbn8yuW27YFd3yOEuRJUoFLBKYClgkMBWwSGDhQyyGBUQsUGGjj86cOUPfc9OmTUlb7jxUFsawRdxY2FXLcaooC31yAyuGXcve3l76tSyoy53/rRBLREJSAYsEpgIWCUwFLBJY+BCLhUFslBJrYyv6v/jii/Q4Dz74YNLGFnKbmprKOvaFCxeStrIF9dgUxSpYYFS2vWiVMIiNVGPfL4BPmWTXiI3EylV2fcumkEagO7BIYCpgkcBUwCKBLVjAZrbFzH5uZifM7DUze6Jo7zazYTM7WTzmr1QuIg2RE2JdBvBNd3/VzK4D8IqZDQP4OwBH3P2gmR0AcADAk83rKg+dalnnaD42he6BBx6gX5s7NY6tf8WCLRbGTE5O0mPfcMMNtL1eLJgqG+2Vu0MB+96cPHkyabvttttyuggAuHTpUtLGArjcNcMih1VlFrwDu/uYu79afDwN4ASAzQD2ARgqvmwIwP3N6qSIcDX9DmxmAwB2A3gJQK+7jwGzRQ5gY6M7JyKfL7uAzexaAM8C+Ia7n6/hdY+b2YiZjZT9iCgi9ckqYDPrxGzx/tDdf1I0j5tZX/H5PgDplnUA3P2Quw+6+2BPT08j+iwihQUTIJtNN74P4IS7f3fOpw4D2A/gYPH4XFN6OAcLSlio8f777ydtbIrh1q1bs74O4IuUs4CH9ZGFPuz92GsBPu2xbDvQepWNxGJYuDQ2Npa0sfW0agmx1q9fn7SdP5/9w1+i7PqyaY9RbjY5Ee4XAfwtgN+Z2fGi7duYLdxnzOxRAKcAPNycLopImQUL2N1/CaBsRvndje2OiNRCI7FEAlMBiwSmAhYJLNR8YDYEkA1nZAunVdmWEuDbYrLk9ty5c0kbG6bIhlJu2LCBHrvRKTRLh8uOzZLbF154IWnbuXNn0rZ79+46evf5quxIUTZcNErizOgOLBKYClgkMBWwSGAqYJHAQoVYDAtzqgQdZeEQC3NYiMV2V2D9YXNYWVgFVFvIjb0nW4yvlnnVbH7y66+/nrTdeOONSRsbHlkLFiZWNTo6mrRt3JhOrmNboLab7sAigamARQJTAYsEpgIWCSxUiMVCI7bSf25w8+677yZtmzdvzu4PG9nD5sqyYIyNDCubr1pl1BULp9hifmxeNcD7OT4+nrSxhfvWrctfqJSdO9uWtcq1KNuZgf3dYHPKG71DRiPoDiwSmApYJDAVsEhgKmCRwEKFWAwbpcTCLjbt8Pnnn0/aHnroIXocFgblBlasP+y1ZTsHlAVMOcdh581GFLH+APx8WEj4yiuvJG233HJL0lY2EosFgh988EHSViVIKtvatLu7O2ljO2wsRroDiwSmAhYJTAUsEpgKWCSwRRlilYU2bFoeC1nYCB4WkuzZsydpO3bsGD02+1rWTzaiiK2zlbtbQ9lx2Oghdt65u1mw0VkAv25sOuFjjz2WtNUSOLF+VgmSWHhXtiYWa2f9KXt9O+kOLBKYClgkMBWwSGAqYJHAQoVYbGobC4jYqCnWxhb0PnXqFD326dOnkza2rhU7DgtUWEjCFlwHgO3bt9P2+dgIq9z+TE9P0/dkQRKbcsnOhx3nqaeeosdhIRhbXywX+zu0Zs2a7Nc3eipjs+gOLBKYClgkMBWwSGAqYJHAVMAigS3KFLpsl4Dcua0MmyvLbNu2jbazlfrffPPNpI2ll2xxNzaUcu3atfTYubsRsOvG5tSyNLYs8WXzhNmQz61bt2b1595776XHyR26yPrD5oTXsjtH7rzuxUh3YJHAVMAigamARQJTAYsEtihDrDK5gRULkhgWDpWFH5OTk0nb+fPnk7be3t6kjQU0bMhm7pDJWuTujsACQoAP72RhW+72pDt27KDt7FrmBlYsVMsdVhqd7sAigamARQJTAYsEtmABm1mXmf3KzH5jZq+Z2XeK9m4zGzazk8Vj/lZ0ItIQOb/Vfwzgy+5+wcw6AfzSzP4HwIMAjrj7QTM7AOAAgCeb2Fcqd2E5NqqHzSVmYRXARy+xQObZZ59N2u64446s92sGtosCG2VUFvCwwKpK39mc7rJ29v1hIRabl83er2z3CfaeDPu7VnY+rbLg0X3WheJpZ/HHAewDMFS0DwG4vyk9FJFSWf98mFmHmR0HMAFg2N1fAtDr7mMAUDymg4VFpKmyCtjdr7j7LgD9AO40s9tyD2Bmj5vZiJmNlP14KiL1qekHeHefAnAUwF4A42bWBwDF40TJaw65+6C7D7I1qESkfguGWGbWA+CSu0+Z2UoAXwHwbwAOA9gP4GDx+FwzO1qGhQi5W2Ky0UfDw8P0OP39/UnbXXfdlbSxXQs2bdqUtJ07dy5pq2UKXK7c9yxbSJBdSzZqqmoox3aaOH78eNLGFtRjC++V7TTBsJFcbJQe2540d6pns+Sk0H0AhsysA7N37Gfc/admdgzAM2b2KIBTAB5uYj9FhFiwgN39twB2k/b3AdzdjE6JSB6NxBIJTAUsEtjSm19Vgk0vm5qaStrKpruxkIa9nq35xMIPFg6xEUVVsVFGbKRaLVPtcqd1MmVTPdnaXTfddFPSxr4PY2NjSRvb2pQFmQAP0C5cuJC0Vdkpoll0BxYJTAUsEpgKWCQwFbBIYEsyxGIhDRuxtWrVqqTt1ltvpe/JhoGy9xwfH0/a2JQ81sey6W6NHqHFplaWjcRiC+KzYIwFeuy8y6busdFqudavX5/1dSysAvi5NyNQbAbdgUUCUwGLBKYCFglMBSwSmApYJLAlmUKzlJVhQxzZ3FKAzxllqWZ3d3fSdvr06aSNzQcuS2jZ8M5GL6ZWy2JzLJnOHWZYNpSSHSd3Ti9L6WtZgK7dC9NVEbfnIqICFolMBSwSmApYJLAlGWIxbJgiC6GOHj1KX793796kLXcHiK1bt2b1p2y+KguNGh28zMzM0HY23DR3aCe7PmVziasMF8197TvvvEPbt23bVvex2013YJHAVMAigamARQJTAYsEtmxCrNzQZ8+ePbSdBUksPGGLw7HXnj17NmnbuJHvD5e7/WUVLKwqwwI4Ft6xYKwsxBodHU3a2BzsWvo535YtW+p+LQBMT08nbWVzjFtFd2CRwFTAIoGpgEUCUwGLBLZsQiyGjRSquvo+G+0zMjKStN13331JWzO2F83Ftlotw4IoFkKxLVnL5I6GYgEaC/lYgFYWgLFtTHfu3Jm0scCKXbcqO1fUSndgkcBUwCKBqYBFAlMBiwS2bEIsFhCtXr06aWOjpgA+0oiN7mJhDGu7fPly1vtVlbs2VC3BC5v2WEtgVQULjViIxdrKvrds7S32eraeF3tt2bpfuceuhe7AIoGpgEUCUwGLBKYCFgls2YRYDAuSWBsAdHV11X0cNhJr165ddb9fLZoRjLFrwRa+Z6Om2JajAB8Vx7DgkWEhFFvIHwD6+vqy+pMbOFUNpmqhO7BIYCpgkcBUwCKBqYBFAssuYDPrMLNfm9lPi+fdZjZsZieLx3XN66aIMLWk0E8AOAFgTfH8AIAj7n7QzA4Uz59scP+aii1Al7s1KZA/F7S3tzdpa+fc32bITYdrmXfM/keApcssaWfXt2xxQJag5x673bLuwGbWD+AvAfznnOZ9AIaKj4cA3N/YronIQnJ/hP4egH8CMHdkfK+7jwFA8UjXRDWzx81sxMxGJicnK3VWRD5rwQI2s78CMOHur9RzAHc/5O6D7j7I1vkVkfrl/A78RQB/bWZfA9AFYI2Z/ReAcTPrc/cxM+sDMNHMjopIasECdvdvAfgWAJjZlwD8o7v/jZn9O4D9AA4Wj881sZ8tUzb0kA2tYyEYw+bK1hKWVZG7o0SrlF0zFm5V2YUh9xgAsH79+qQtd2hnu1X5f+CDAO4xs5MA7imei0gL1TSZwd2PAjhafPw+gLsb3yURyaWRWCKBqYBFAlvW84GZsvCChVu5QUczAqvccGqxjfgquxZV5lvnqiUUq/I9y11IsBF0BxYJTAUsEpgKWCQwFbBIYAqx5qklbGBBB9u1IDegKQvF2HFYONXOUVdswbiVK1cmba0agcaULWqXuytF7rVsVmBFj9WyI4lIw6mARQJTAYsEpgIWCUwhVqaZmZmkjYVTuYEVW3Mpd3riYlTL9qStwHaFYKFadLoDiwSmAhYJTAUsEpgKWCSwuKlJizV6faZmjJBq59TBxTZtkS3CPjU1Rb+2bMvTCHQHFglMBSwSmApYJDAVsEhgKmCRwJRCNxhb0IzN860ltWWvb+e8WtYftuvBihUrWtGdbBpKKSKLigpYJDAVsEhgKmCRwBRiNVgzFjRrZ2DFsP4stsCKidDHWukOLBKYClgkMBWwSGAqYJHALHeLzIYczGwSwDvF0w0AzrTs4M21lM4FWFrnsxTOZZu797BPtLSAP3NgsxF3H2zLwRtsKZ0LsLTOZymdC6MfoUUCUwGLBNbOAj7UxmM32lI6F2Bpnc9SOpdE234HFpHq9CO0SGAtL2Az22tmb5jZW2Z2oNXHr8rMnjazCTP7/Zy2bjMbNrOTxeO6dvYxl5ltMbOfm9kJM3vNzJ4o2qOeT5eZ/crMflOcz3eK9pDnk6OlBWxmHQD+A8BXAewA8IiZ7WhlHxrgBwD2zms7AOCIu98M4EjxPILLAL7p7n8G4AsA/r74fkQ9n48BfNnd/wLALgB7zewLiHs+C2r1HfhOAG+5+x/c/RMAPwawr8V9qMTdfwHg7LzmfQCGio+HANzf0k7Vyd3H3P3V4uNpACcAbEbc83F3v1A87Sz+OIKeT45WF/BmAKfnPB8t2qLrdfcxYLYoAGxsc39qZmYDAHYDeAmBz8fMOszsOIAJAMPuHvp8FtLqAmYTWxWDt5mZXQvgWQDfcPfz7e5PFe5+xd13AegHcKeZ3dbuPjVTqwt4FMCWOc/7AfyxxX1ohnEz6wOA4nGizf3JZmadmC3eH7r7T4rmsOfzJ+4+BeAoZvOK8OdTptUF/DKAm81su5ldA+DrAA63uA/NcBjA/uLj/QCea2Nfstns0hrfB3DC3b8751NRz6fHzNYWH68E8BUA/4eg55Oj5QM5zOxrAL4HoAPA0+7+ry3tQEVm9iMAX8LsLJdxAP8M4L8BPANgK4BTAB529/lB16JjZnsA/C+A3wH404LW38bs78ERz+fPMRtSdWD25vSMu/+Lma1HwPPJoZFYIoFpJJZIYCpgkcBUwCKBqYBFAlMBiwSmAhYJTAUsEpgKWCSw/wdW/VuXx+SNnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(density_images[0], vmin=0.0, vmax=1.0, cmap='binary')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAEYCAYAAABiNXvJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAMU0lEQVR4nO3dX6hlZ3nH8e+v0ymR2uKkmjBkksaLUCrSRhiCUC9CMTC1pTMUArEUplCYmxYiFOrUG7FQyJV405uhBg+0KAMWZ8iNDFODFkpMTBVNx3SC0Jg6zFAG0dxUNE8vzpIcz+wze52z/z77fD+w2Xut2X/ed+358e71rHetk6pCUk+/tOoGSDo4Ayw1ZoClxgyw1JgBlhozwFJjMwU4yakkryZ5Lcn5eTVK0jg56HHgJEeA/wKeAN4AXgQ+WlX/eZfXeNBZOoCqyqT1s4zAjwGvVdX3quonwBeA0zO8n6R9miXADwDf37H8xrBO0pL88gyvnTSk3/ETOck54NwMnyNpD7ME+A3gwR3LJ4Af7H5SVV0ALoD7wNK8zRLgF4FHkrwX+B/gKeBP59KqFVrlyR3JxDqFtKcDB7iqfprkr4AvA0eAZ6vqlbm1TNJUBz6MdKAPa/AT2hFY62gRh5EkrZgBlhqbpYilOdvPz/euP7cPQx+XyRFYaswAS40ZYKkxAyw1ZhFrl70KJ+t2+d2x7Vm3QlCX7duFI7DUmAGWGjPAUmMGWGrMAEuNWYXecJOqu5MqwYuY4jhrZblrpX2ZHIGlxgyw1JgBlhozwFJjFrF2cUrfwaxyu40t1G0iR2CpMQMsNWaApcYMsNSYRawNt4hijoW+9eEILDVmgKXGDLDUmAGWGrOIJcCLzXXlCCw1ZoClxgyw1JgBlhqziLVBFjHratbrZx1Gy/wTqo7AUmMGWGrMAEuNGWCpMYtYuxyGGUn7uYZU134v6zpZq94+jsBSYwZYaswAS40ZYKmxqQFO8mySW0m+s2PdvUmuJLk+3B9bbDMlTTJmBP4ccGrXuvPA1ap6BLg6LLdTVaNv6ybJHbdZdOm3ftHUAFfVV4Hbu1afBraGx1vAmTm3S9IIBz0OfH9V3QCoqhtJ7tvriUnOAecO+DmS7mLhEzmq6gJwASCJv8mkOTpoFfpmkuMAw/2t+TVJ0lgHDfBl4Ozw+CxwaT7NWa5JhaC9bl1ZmHrbfgp18y5uLur/VaY1IsnngceBdwM3gU8CXwIuAg8BrwNPVtXuQtek92r7v2fd/uOP/fLXrd3raBkXLZjDUYKJbzA1wPNkgOfHAM9P5wA7E0tqzABLjXk+8Ehe3O1w6fJ9OwJLjRlgqTEDLDVmgKXGLGKNtG4FjGVdtO0wWLfvdj8cgaXGDLDUmAGWGjPAUmNrWcTqXFRYJbfb4eMILDVmgKXGDLDUmAGWGlvLIpa0afYqMM46e84RWGrMAEuNGWCpMQMsNWaApcasQjfV5aJr2raoc7UdgaXGDLDUmAGWGjPAUmMWsaQlcCqlpDsYYKkxAyw1ZoClxixiSROMLS6tevabI7DUmAGWGjPAUmMGWGrMIpY20iJOtxz7+mWe6ukILDVmgKXGDLDU2NQAJ3kwyVeSXEvySpKnh/X3JrmS5Ppwf2zxzZW0U6btXCc5DhyvqpeT/BrwDeAM8OfA7ap6Jsl54FhVfXzKe43ak1/17JYOuswUWoZFbItlbd99fM7EJ04dgavqRlW9PDz+MXANeAA4DWwNT9tiO9SSlmhf+8BJHgY+ALwA3F9VN2A75MB9826cpLsbfRw4yTuBLwIfq6ofjR36k5wDzh2seZLuZuo+MECSo8BzwJer6tPDuleBx6vqxrCf/HxV/daU93EfeE7cB36b+8B3/4AAnwWu/Ty8g8vA2eHxWeDSqJZoYarqjtthMKnfs26LRbznIoypQn8I+BrwbeCtYfUn2N4Pvgg8BLwOPFlVt6e8lyPwnPiXGTbDrCPwqJ/Q82KA58cAb4aF/4SWtL4MsNSYAZYaM8BSYwZYaswAS40ZYKkxAyw15kXtmnLSxtsO86QWR2CpMQMsNWaApcYMsNSYRSy1d1gKVpM4AkuNGWCpMQMsNWaApcYsYklzNvYyOfPgCCw1ZoClxgyw1JgBlhqziCWNNEtxaq/ZYrMWvByBpcYMsNSYAZYaM8BSYxax1N6yrok16T2XOetqEkdgqTEDLDVmgKXGDLDUmAGWGrMK3cBh/ssDY6xyW6z6e3AElhozwFJjBlhqzABLjVnEkpZgUVMuHYGlxgyw1JgBlhqbGuAk9yT5epJvJXklyaeG9fcmuZLk+nB/bPHNlbRTps0kyfbe969W1ZtJjgL/BjwN/Alwu6qeSXIeOFZVH5/yXqOmrax6dsu6GVsAWbft5gyyt81axKqqiW8wdQSubW8Oi0eHWwGnga1h/RZwZqYWStq3UfvASY4k+SZwC7hSVS8A91fVDYDh/r7FNVPSJKMCXFU/q6pHgRPAY0neP/YDkpxL8lKSlw7aSEmT7asKXVU/BJ4HTgE3kxwHGO5v7fGaC1V1sqpOzthWSbuMqUK/J8m7hsfvAD4MfBe4DJwdnnYWuLSoRqqnqrrjdhgkueM2aVvMY3uMmUp5HNhKcoTtwF+squeS/DtwMclfAK8DT87cGkn7MvUw0lw/zMNIB9L1MNJhtZ/DZ/v4bg92GEnS+jLAUmOtTid0Zs/duX3WwzK3uSOw1JgBlhozwFJjBlhqrFURy4LM/lnYWqxVb19HYKkxAyw1ZoClxgyw1JgBlhpbyyr0XmdoHNbq6aR+e4bSelj19nUElhozwFJjBlhqzABLja1lEWsvs1zdftXFBvWy6imSYzkCS40ZYKkxAyw1ZoClxloVscZax2LDLJx1tXxdtqUjsNSYAZYaM8BSYwZYamwji1hdZtGM1bntWixHYKkxAyw1ZoClxgyw1NhGFrEOg2UV6jatILhpHIGlxgyw1JgBlhozwFJjFrGaspAkcASWWjPAUmMGWGrMAEuNjQ5wkiNJ/iPJc8PyvUmuJLk+3B9bXDMlTbKfEfhp4NqO5fPA1ap6BLg6LK+tJHfcuujcdi3WqAAnOQH8IfCPO1afBraGx1vAmfk2TdI0Y0fgzwB/A7y1Y939VXUDYLi/b9ILk5xL8lKSl2ZqqaQ7TA1wkj8CblXVNw7yAVV1oapOVtXJg7xe0t7GzMT6PeCPk3wEuAf49ST/BNxMcryqbiQ5DtxaZEMl3WnqCFxVf1tVJ6rqYeAp4F+r6s+Ay8DZ4WlngUsLa6WkiWY5DvwM8ESS68ATw7KkJcoyJ8UnWasZ+F1OCFjlVTG8Isdi7ePvXk18ojOxpMYMsNSY5wNr3/xZvT4cgaXGDLDUmAGWGjPAUmMGWGrMAEuNGWCpMQMsNWaApcacidXAKmc5TfrsRVyTax+T+uf+2Z05AkuNGWCpMQMsNWaApcYsYmktjC1OeSrjL3IElhozwFJjBlhqzABLjRlgqTGr0Nq3VV6TeuzzDktl2hFYaswAS40ZYKkxAyw1dqiLWIs4r1Xr4bB8t47AUmMGWGrMAEuNGWCpsWUXsf4X+O/h8buH5U2wSX2BzerPJvTlN/f6h6xqylmSl6rq5Eo+fM42qS+wWf3ZpL5M4k9oqTEDLDW2ygBfWOFnz9sm9QU2qz+b1Jc7rGwfWNLs/AktNbb0ACc5leTVJK8lOb/sz59VkmeT3ErynR3r7k1yJcn14f7YKts4VpIHk3wlybUkryR5eljftT/3JPl6km8N/fnUsL5lf8ZYaoCTHAH+AfgD4H3AR5O8b5ltmIPPAad2rTsPXK2qR4Crw3IHPwX+uqp+G/gg8JfD99G1P/8H/H5V/S7wKHAqyQfp25+plj0CPwa8VlXfq6qfAF8ATi+5DTOpqq8Ct3etPg1sDY+3gDNLbdQBVdWNqnp5ePxj4BrwAH37U1X15rB4dLgVTfszxrID/ADw/R3Lbwzruru/qm7AdiiA+1bcnn1L8jDwAeAFGvcnyZEk3wRuAVeqqnV/pll2gCedpGkZfMWSvBP4IvCxqvrRqtszi6r6WVU9CpwAHkvy/lW3aZGWHeA3gAd3LJ8AfrDkNizCzSTHAYb7Wytuz2hJjrId3n+uqn8ZVrftz89V1Q+B59muV7Tvz16WHeAXgUeSvDfJrwBPAZeX3IZFuAycHR6fBS6tsC2jZfuyFZ8FrlXVp3f8U9f+vCfJu4bH7wA+DHyXpv0ZY+kTOZJ8BPgMcAR4tqr+fqkNmFGSzwOPs32Wy03gk8CXgIvAQ8DrwJNVtbvQtXaSfAj4GvBt4K1h9SfY3g/u2J/fYbtIdYTtweliVf1dkt+gYX/GcCaW1JgzsaTGDLDUmAGWGjPAUmMGWGrMAEuNGWCpMQMsNfb/yFmNtLnOy0wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(1 - density_mask, vmin=0.0, vmax=1.0, cmap='binary')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Flatten images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_images = density_images.reshape(density_images.shape[0], -1).astype('float32')\n",
    "flat_mask   = density_mask.reshape(-1).astype('float32')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load the autoencoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.02868951 -1.9072311 ]] [[0.5830229 2.27596  ]]\n"
     ]
    }
   ],
   "source": [
    "train_mu  = np.load('temporal_anomaly_detection/train_mu.npy')\n",
    "train_std = np.load('temporal_anomaly_detection/train_std.npy')\n",
    "print(train_mu, train_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "density_mask_sum = density_mask.sum()\n",
    "    \n",
    "def scaled_loss(inputs, outputs):\n",
    "  outputs *= flat_mask\n",
    "  loss = tfk.metrics.mean_squared_error(inputs, outputs)\n",
    "  return density_mask_sum * loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder = tfk.models.load_model('temporal_anomaly_detection/autoencoder.h5',\n",
    "                                    custom_objects={'scaled_loss':scaled_loss})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measure the mahalanobis distances of the encodings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "latent_vars = autoencoder.get_layer('encoder')(flat_images).numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5856,)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mahalanobis_dist = (((latent_vars - train_mu)/train_std)**2).sum(axis=1)**(0.5)\n",
    "mahalanobis_dist.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find the anomalies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data at time T is said to be anomalous if the mahalanobis distance (a multidimensional analogue of the standard deviation) of the latent encoding of the demand density map is greater than 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['is_anomaly'] = np.array(list(map(lambda x : mahalanobis_dist[x] > 4.0, timesteps)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geohash6</th>\n",
       "      <th>day</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>y</th>\n",
       "      <th>hour</th>\n",
       "      <th>minutes</th>\n",
       "      <th>ds</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>is_anomaly</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>qp09d7</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.091915</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.342102</td>\n",
       "      <td>90.807495</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>qp0d5g</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.254211</td>\n",
       "      <td>90.873413</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>qp03pw</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.016879</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.413513</td>\n",
       "      <td>90.686646</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>qp09jr</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.408020</td>\n",
       "      <td>90.939331</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>qp02z1</td>\n",
       "      <td>1</td>\n",
       "      <td>0:0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1900-01-01</td>\n",
       "      <td>-5.484924</td>\n",
       "      <td>90.664673</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  geohash6  day timestamp         y  hour  minutes         ds  latitude  \\\n",
       "0   qp09d7    1       0:0  0.091915     0        0 1900-01-01 -5.342102   \n",
       "1   qp0d5g    1       0:0  0.000000     0        0 1900-01-01 -5.254211   \n",
       "2   qp03pw    1       0:0  0.016879     0        0 1900-01-01 -5.413513   \n",
       "3   qp09jr    1       0:0  0.000000     0        0 1900-01-01 -5.408020   \n",
       "4   qp02z1    1       0:0  0.000000     0        0 1900-01-01 -5.484924   \n",
       "\n",
       "   longitude  is_anomaly  \n",
       "0  90.807495       False  \n",
       "1  90.873413       False  \n",
       "2  90.686646       False  \n",
       "3  90.939331       False  \n",
       "4  90.664673       False  "
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Group the data by geohash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 11 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# df_cut = df[df['ds'] >= datetime.datetime.strptime(str(61-14), \"%j\")]\n",
    "# df_cut = df_cut.reset_index().drop('index', axis=1)\n",
    "# gb = df_cut.groupby('geohash6')\n",
    "\n",
    "gb = df.groupby('geohash6')\n",
    "location_series = {}\n",
    "for x in gb.groups:\n",
    "    group = gb.get_group(x).reset_index().drop('index', axis=1)\n",
    "    location_series[x] = group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_inputs = [(location, location_series[location]) for location in locations]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Timestamp('1900-01-15 00:00:00'),\n",
       " Timestamp('1900-01-15 00:15:00'),\n",
       " Timestamp('1900-01-15 00:30:00'),\n",
       " Timestamp('1900-01-15 00:45:00'),\n",
       " Timestamp('1900-01-15 01:00:00')]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "future_ds = [max_ds + i*delta_ds for i in range(1, 5+1)]\n",
    "future_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stage 1 (faster but less accurate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warning: had problems running the following blocks of code on windows. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run kalman filter and the theta forecasters in parallel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stage1_worker(data_input):\n",
    "    location, df_series = data_input\n",
    "    \n",
    "    # create train-test x values\n",
    "    X_train = np.arange(len(df_series)).reshape(-1, 1)\n",
    "    X_test  = (np.arange(5) + len(df_series)).reshape(-1, 1)\n",
    "    \n",
    "    # forecast using data from previous daily cycle\n",
    "    # prev_pred = df_series['y'].values       <- deprecated\n",
    "    \n",
    "    # forecast using the kalman filter\n",
    "    kf = simdkalman.KalmanFilter(\n",
    "      state_transition = np.array([[1,1],[0,1]]),\n",
    "      process_noise = np.diag([0.1, 0.01]),\n",
    "      observation_model = np.array([[1,0]]),\n",
    "      observation_noise = 0.5)\n",
    "    kalman_pred = kf.predict(df_series['y'].values, 5).states.mean[:,0]\n",
    "    \n",
    "    # deseasonalize\n",
    "    decomposition = sm.tsa.seasonal_decompose(df_series['y'] + 1e-9,\n",
    "                                              model='multiplicative',\n",
    "                                              freq=4*24,\n",
    "                                              two_sided=False,\n",
    "                                              extrapolate_trend=1)\n",
    "    df_series['deseasonalized'] = df_series['y'].values / (decomposition.seasonal + 1e-9)\n",
    "    seasonal_test = decomposition.seasonal.values[-4*24:-4*24+5]\n",
    "    \n",
    "    # forecast using theta = 0 (with linear regression)\n",
    "    linreg = HuberRegressor().fit(X_train, df_series['deseasonalized'].values)\n",
    "    theta0_pred = linreg.predict(X_test) * seasonal_test\n",
    "    \n",
    "    # forecast using theta = 2 (with simple exponential smoothing)\n",
    "    y_theta2 = 2*df_series['deseasonalized'].values - linreg.predict(X_train)\n",
    "    ses = SimpleExpSmoothing(y_theta2).fit()\n",
    "    theta2_pred = ses.forecast(5) * seasonal_test\n",
    "    \n",
    "    # get the mean of the forecasts for theta = 0 and theta = 2\n",
    "    theta_pred = (theta0_pred + theta2_pred) / 2.0\n",
    "    \n",
    "    # clip the predictions\n",
    "    kalman_pred = np.clip(kalman_pred, a_min=0.0, a_max=1.0)\n",
    "    theta_pred  = np.clip(theta_pred , a_min=0.0, a_max=1.0)\n",
    "    \n",
    "    # return_dict[location] = {'kalman_pred' : kalman_pred,\n",
    "    #                          'theta_pred' : theta_pred}\n",
    "    return location, kalman_pred, theta_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following block takes around 30s on an 16-core machine. We can ignore the divide by zero warnings since they automatically get handled by statsmodels anyway."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of CPUs used: 4\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "print(\"# of CPUs used:\", multiprocessing.cpu_count())\n",
    "p = multiprocessing.Pool(processes=multiprocessing.cpu_count())\n",
    "stage1_results = p.map(stage1_worker, data_inputs)\n",
    "p.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Aggregate the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stage1_weights = [0.2969419 , 0.67783948]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "stage1_dfprep = []\n",
    "for result in stage1_results:\n",
    "    location, kalman_pred, theta_pred = result\n",
    "    stage1_pred = stage1_weights[0]*kalman_pred + stage1_weights[1]*theta_pred\n",
    "    for i in range(5):\n",
    "        x = {'geohash6' : location,\n",
    "             'day' : future_ds[i].timetuple().tm_yday,\n",
    "             'timestamp' : str(future_ds[i].hour) + ':' + str(future_ds[i].minute),\n",
    "             'demand' : stage1_pred[i]}\n",
    "        stage1_dfprep.append(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_stage1_results = pd.DataFrame(stage1_dfprep)\n",
    "df_stage1_results.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_stage1_results.to_csv(\"results_partial.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stage 2 (more accurate but slower)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run FBProphet in parallel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure to run stage 1 first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stage2_worker(data_input):\n",
    "    location, df_series = data_input\n",
    "    df_series['y'] += 1e-9\n",
    "    \n",
    "    prophet = Prophet(seasonality_mode='multiplicative',\n",
    "                      weekly_seasonality=True,\n",
    "                      yearly_seasonality=False)\n",
    "    prophet.add_regressor('is_anomaly')\n",
    "    prophet.fit(df_series)\n",
    "    \n",
    "    future = prophet.make_future_dataframe(periods=5, freq='15min', include_history=False)\n",
    "    future['is_anomaly'] = 0\n",
    "    prophet_pred = prophet.predict(future)['yhat'].values\n",
    "    \n",
    "    # clip\n",
    "    prophet_pred  = np.clip(prophet_pred, a_min=0.0, a_max=1.0)\n",
    "    \n",
    "    return location, prophet_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Takes around 5 mins on a 16-core machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.71 s, sys: 1.47 s, total: 3.18 s\n",
      "Wall time: 5min 7s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "p = multiprocessing.Pool(processes=multiprocessing.cpu_count())\n",
    "stage2_results = p.map(stage2_worker, data_inputs)\n",
    "p.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Aggregate the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "stage2_weights = [0.309658, 0.38132017, 0.29395973]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "stage2_preds = {}\n",
    "stage2_dfprep = []\n",
    "for result in stage1_results:\n",
    "    location, kalman_pred, theta_pred = result\n",
    "    stage2_preds[location] = stage2_weights[0]*kalman_pred + stage2_weights[1]*theta_pred\n",
    "for result in stage2_results:\n",
    "    location, prophet_pred = result\n",
    "    stage2_preds[location] += stage2_weights[2]*prophet_pred\n",
    "    for i in range(5):\n",
    "        x = {'geohash6' : location,\n",
    "             'day' : future_ds[i].timetuple().tm_yday,\n",
    "             'timestamp' : str(future_ds[i].hour) + ':' + str(future_ds[i].minute),\n",
    "             'demand' : stage2_preds[location][i]}\n",
    "        stage2_dfprep.append(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>day</th>\n",
       "      <th>demand</th>\n",
       "      <th>geohash6</th>\n",
       "      <th>timestamp</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62</td>\n",
       "      <td>0.084853</td>\n",
       "      <td>qp03wc</td>\n",
       "      <td>0:0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>62</td>\n",
       "      <td>0.090412</td>\n",
       "      <td>qp03wc</td>\n",
       "      <td>0:15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>62</td>\n",
       "      <td>0.094315</td>\n",
       "      <td>qp03wc</td>\n",
       "      <td>0:30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>62</td>\n",
       "      <td>0.097274</td>\n",
       "      <td>qp03wc</td>\n",
       "      <td>0:45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>62</td>\n",
       "      <td>0.103081</td>\n",
       "      <td>qp03wc</td>\n",
       "      <td>1:0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   day    demand geohash6 timestamp\n",
       "0   62  0.084853   qp03wc       0:0\n",
       "1   62  0.090412   qp03wc      0:15\n",
       "2   62  0.094315   qp03wc      0:30\n",
       "3   62  0.097274   qp03wc      0:45\n",
       "4   62  0.103081   qp03wc       1:0"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_stage2_results = pd.DataFrame(stage2_dfprep)\n",
    "df_stage2_results.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_stage2_results.to_csv(\"results_final.csv\", index=False)"
   ]
  },
  {
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}