Chapter01/chapter1-time-series-analysis-overview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b7105c6",
   "metadata": {},
   "source": [
    "# Time series analysis on AWS\n",
    "*Chapter 1 - Time series analysis overview*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b393e4ec",
   "metadata": {},
   "source": [
    "## Initializations\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4dd4b8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --quiet tqdm kaggle tsia ruptures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9997f46",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd65af91",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.colors as mpl_colors\n",
    "import matplotlib.dates as mdates\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import ruptures as rpt\n",
    "import sys\n",
    "import tsia\n",
    "import warnings\n",
    "import zipfile\n",
    "\n",
    "from matplotlib import gridspec\n",
    "from sklearn.preprocessing import normalize\n",
    "from tqdm import tqdm\n",
    "from urllib.request import urlretrieve"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d104af8",
   "metadata": {},
   "source": [
    "### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13989034",
   "metadata": {},
   "outputs": [],
   "source": [
    "RAW_DATA = os.path.join('..', 'Data', 'raw')\n",
    "DATA = os.path.join('..', 'Data')\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "os.makedirs(RAW_DATA, exist_ok=True)\n",
    "\n",
    "%matplotlib inline\n",
    "# plt.style.use('Solarize_Light2')\n",
    "plt.style.use('fivethirtyeight')\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "plt.rcParams['figure.dpi'] = 300\n",
    "plt.rcParams['lines.linewidth'] = 0.3\n",
    "plt.rcParams['axes.titlesize'] = 6\n",
    "plt.rcParams['axes.labelsize'] = 6\n",
    "plt.rcParams['xtick.labelsize'] = 4.5\n",
    "plt.rcParams['ytick.labelsize'] = 4.5\n",
    "plt.rcParams['grid.linewidth'] = 0.2\n",
    "plt.rcParams['legend.fontsize'] = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71228c3a",
   "metadata": {},
   "source": [
    "### Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ead7b9b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def progress_report_hook(count, block_size, total_size):\n",
    "    mb = int(count * block_size // 1048576)\n",
    "    if count % 500 == 0:\n",
    "        sys.stdout.write(\"\\r{} MB downloaded\".format(mb))\n",
    "        sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "590ace42",
   "metadata": {},
   "source": [
    "### Downloading datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1a10de9",
   "metadata": {},
   "source": [
    "#### **Dataset 1:** Household energy consumption"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7821f0b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ORIGINAL_DATA = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00321/LD2011_2014.txt.zip'\n",
    "ARCHIVE_PATH  = os.path.join(RAW_DATA, 'energy-consumption.zip')\n",
    "FILE_NAME     = 'energy-consumption.csv'\n",
    "FILE_PATH     = os.path.join(DATA, 'energy', FILE_NAME)\n",
    "FILE_DIR      = os.path.dirname(FILE_PATH)\n",
    "\n",
    "if not os.path.isfile(FILE_PATH):\n",
    "    print(\"Downloading dataset (258MB), can take a few minutes depending on your connection\")\n",
    "    urlretrieve(ORIGINAL_DATA, ARCHIVE_PATH, reporthook=progress_report_hook)\n",
    "    os.makedirs(os.path.join(DATA, 'energy'), exist_ok=True)\n",
    "\n",
    "    print(\"\\nExtracting data archive\")\n",
    "    zip_ref = zipfile.ZipFile(ARCHIVE_PATH, 'r')\n",
    "    zip_ref.extractall(FILE_DIR + '/')\n",
    "    zip_ref.close()\n",
    "    \n",
    "    !rm -Rf $FILE_DIR/__MACOSX\n",
    "    !mv $FILE_DIR/LD2011_2014.txt $FILE_PATH\n",
    "    \n",
    "else:\n",
    "    print(\"File found, skipping download\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "171badbf",
   "metadata": {},
   "source": [
    "#### **Dataset 2:** Nasa Turbofan remaining useful lifetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "244b7160",
   "metadata": {},
   "outputs": [],
   "source": [
    "ok = True\n",
    "ok = ok and os.path.exists(os.path.join(DATA, 'turbofan', 'train_FD001.txt'))\n",
    "ok = ok and os.path.exists(os.path.join(DATA, 'turbofan', 'test_FD001.txt'))\n",
    "ok = ok and os.path.exists(os.path.join(DATA, 'turbofan', 'RUL_FD001.txt'))\n",
    "\n",
    "if (ok):\n",
    "    print(\"File found, skipping download\")\n",
    "\n",
    "else:\n",
    "    print('Some datasets are missing, create working directories and download original dataset from the NASA repository.')\n",
    "    \n",
    "    # Making sure the directory already exists:\n",
    "    os.makedirs(os.path.join(DATA, 'turbofan'), exist_ok=True)\n",
    "\n",
    "    # Download the dataset from the NASA repository, unzip it and set\n",
    "    # aside the first training file to work on:\n",
    "    !wget https://ti.arc.nasa.gov/c/6/ --output-document=$RAW_DATA/CMAPSSData.zip\n",
    "    !unzip $RAW_DATA/CMAPSSData.zip -d $RAW_DATA\n",
    "    !cp $RAW_DATA/train_FD001.txt $DATA/turbofan/train_FD001.txt\n",
    "    !cp $RAW_DATA/test_FD001.txt $DATA/turbofan/test_FD001.txt\n",
    "    !cp $RAW_DATA/RUL_FD001.txt $DATA/turbofan/RUL_FD001.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eedc9c2",
   "metadata": {},
   "source": [
    "#### **Dataset 3:** Human heartbeat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1314ead",
   "metadata": {},
   "outputs": [],
   "source": [
    "ECG_DATA_SOURCE = 'http://www.timeseriesclassification.com/Downloads/ECG200.zip'\n",
    "ARCHIVE_PATH  = os.path.join(RAW_DATA, 'ECG200.zip')\n",
    "FILE_NAME     = 'ecg.csv'\n",
    "FILE_PATH     = os.path.join(DATA, 'ecg', FILE_NAME)\n",
    "FILE_DIR      = os.path.dirname(FILE_PATH)\n",
    "\n",
    "if not os.path.isfile(FILE_PATH):\n",
    "    urlretrieve(ECG_DATA_SOURCE, ARCHIVE_PATH)\n",
    "    os.makedirs(os.path.join(DATA, 'ecg'), exist_ok=True)\n",
    "\n",
    "    print(\"\\nExtracting data archive\")\n",
    "    zip_ref = zipfile.ZipFile(ARCHIVE_PATH, 'r')\n",
    "    zip_ref.extractall(FILE_DIR + '/')\n",
    "    zip_ref.close()\n",
    "    \n",
    "    !mv $DATA/ecg/ECG200_TRAIN.txt $FILE_PATH\n",
    "    \n",
    "else:\n",
    "    print(\"File found, skipping download\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a297084",
   "metadata": {},
   "source": [
    "#### **Dataset 4:** Industrial pump data\n",
    "To download this dataset from Kaggle, you will need to have an account and create a token that you install on your machine. You can follow [**this link**](https://www.kaggle.com/docs/api) to get started with the Kaggle API. Once generated, make sure your Kaggle token is stored in the `~/.kaggle/kaggle.json` file, or the next cells will issue an error. In some cases, you may still have an error while using this location. Try moving your token in this location instead: `~/kaggle/kaggle.json` (not the absence of the `.` in the folder name).\n",
    "\n",
    "To get a Kaggle token, go to kaggle.com and create an account. Then navigate to **My account** and scroll down to the API section. There, click the **Create new API token** button:\n",
    "\n",
    "<img src=\"../Assets/kaggle_api.png\" />\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3617c828",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_NAME    = 'pump-sensor-data.zip'\n",
    "ARCHIVE_PATH = os.path.join(RAW_DATA, FILE_NAME)\n",
    "FILE_PATH    = os.path.join(DATA, 'pump', 'sensor.csv')\n",
    "FILE_DIR     = os.path.dirname(FILE_PATH)\n",
    "\n",
    "if not os.path.isfile(FILE_PATH):\n",
    "    if not os.path.exists('/home/ec2-user/.kaggle/kaggle.json'):\n",
    "        os.makedirs('/home/ec2-user/.kaggle/', exist_ok=True)\n",
    "        raise Exception('The kaggle.json token was not found.\\nCreating the /home/ec2-user/.kaggle/ directory: put your kaggle.json file there once you have generated it from the Kaggle website')\n",
    "    else:\n",
    "        print('The kaggle.json token file was found: making sure it is not readable by other users on this system.')\n",
    "        !chmod 600 /home/ec2-user/.kaggle/kaggle.json\n",
    "\n",
    "    os.makedirs(os.path.join(DATA, 'pump'), exist_ok=True)\n",
    "    !kaggle datasets download -d nphantawee/pump-sensor-data -p $RAW_DATA\n",
    "\n",
    "    print(\"\\nExtracting data archive\")\n",
    "    zip_ref = zipfile.ZipFile(ARCHIVE_PATH, 'r')\n",
    "    zip_ref.extractall(FILE_DIR + '/')\n",
    "    zip_ref.close()\n",
    "    \n",
    "else:\n",
    "    print(\"File found, skipping download\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c3cc24",
   "metadata": {},
   "source": [
    "#### **Dataset 5:** London household energy consumption with weather data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5b137b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_NAME    = 'smart-meters-in-london.zip'\n",
    "ARCHIVE_PATH = os.path.join(RAW_DATA, FILE_NAME)\n",
    "FILE_PATH    = os.path.join(DATA, 'energy-london', 'smart-meters-in-london.zip')\n",
    "FILE_DIR     = os.path.dirname(FILE_PATH)\n",
    "\n",
    "# Checks if the data were already downloaded:\n",
    "if os.path.exists(os.path.join(DATA, 'energy-london', 'acorn_details.csv')):\n",
    "    print(\"File found, skipping download\")\n",
    "    \n",
    "else:\n",
    "    # Downloading and unzipping datasets from Kaggle:\n",
    "    print(\"Downloading dataset (2.26G), can take a few minutes depending on your connection\")\n",
    "    os.makedirs(os.path.join(DATA, 'energy-london'), exist_ok=True)\n",
    "    !kaggle datasets download -d jeanmidev/smart-meters-in-london -p $RAW_DATA\n",
    "    \n",
    "    print('Unzipping files...')\n",
    "    zip_ref = zipfile.ZipFile(ARCHIVE_PATH, 'r')\n",
    "    zip_ref.extractall(FILE_DIR + '/')\n",
    "    zip_ref.close()\n",
    "    \n",
    "    !rm $DATA/energy-london/*zip\n",
    "    !rm $DATA/energy-london/*gz\n",
    "    !mv $DATA/energy-london/halfhourly_dataset/halfhourly_dataset/* $DATA/energy-london/halfhourly_dataset\n",
    "    !rm -Rf $DATA/energy-london/halfhourly_dataset/halfhourly_dataset\n",
    "    !mv $DATA/energy-london/daily_dataset/daily_dataset/* $DATA/energy-london/daily_dataset\n",
    "    !rm -Rf $DATA/energy-london/daily_dataset/daily_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5bd6e2e",
   "metadata": {},
   "source": [
    "## Dataset visualization\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bef12b32",
   "metadata": {},
   "source": [
    "### **1.** Household energy consumption"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9649de6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "FILE_PATH = os.path.join(DATA, 'energy', 'energy-consumption.csv')\n",
    "energy_df = pd.read_csv(FILE_PATH, sep=';', decimal=',')\n",
    "energy_df = energy_df.rename(columns={'Unnamed: 0': 'Timestamp'})\n",
    "energy_df['Timestamp'] = pd.to_datetime(energy_df['Timestamp'])\n",
    "energy_df = energy_df.set_index('Timestamp')\n",
    "energy_df.iloc[100000:, 1:5].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acc364f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5, 1.876))\n",
    "plt.plot(energy_df['MT_002'])\n",
    "plt.title('Energy consumption for household MT_002')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5ab3461",
   "metadata": {},
   "source": [
    "### **2.** NASA Turbofan data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2e44083",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_PATH = os.path.join(DATA, 'turbofan', 'train_FD001.txt')\n",
    "turbofan_df = pd.read_csv(FILE_PATH, header=None, sep=' ')\n",
    "turbofan_df.dropna(axis='columns', how='all', inplace=True)\n",
    "print('Shape:', turbofan_df.shape)\n",
    "turbofan_df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8342359",
   "metadata": {},
   "outputs": [],
   "source": [
    "columns = [\n",
    "    'unit_number',\n",
    "    'cycle',\n",
    "    'setting_1',\n",
    "    'setting_2',\n",
    "    'setting_3',\n",
    "] + ['sensor_{}'.format(s) for s in range(1,22)]\n",
    "turbofan_df.columns = columns\n",
    "turbofan_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32788dc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a RUL column and group the data by unit_number:\n",
    "turbofan_df['rul'] = 0\n",
    "grouped_data = turbofan_df.groupby(by='unit_number')\n",
    "\n",
    "# Loops through each unit number to get the lifecycle counts:\n",
    "for unit, rul in enumerate(grouped_data.count()['cycle']):\n",
    "    current_df = turbofan_df[turbofan_df['unit_number'] == (unit+1)].copy()\n",
    "    current_df['rul'] = rul - current_df['cycle']\n",
    "    turbofan_df[turbofan_df['unit_number'] == (unit+1)] = current_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b1fe985",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = turbofan_df.iloc[:, [0,1,2,3,4,5,6,25,26]].copy()\n",
    "df = df[df['unit_number'] == 1]\n",
    "\n",
    "def highlight_cols(s):\n",
    "    return f'background-color: rgba(0, 143, 213, 0.3)'\n",
    "\n",
    "df.head(10).style.applymap(highlight_cols, subset=['rul'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f313864",
   "metadata": {},
   "source": [
    "### **3.** ECG Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5fc7c8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_PATH = os.path.join(DATA, 'ecg', 'ecg.csv')\n",
    "ecg_df = pd.read_csv(FILE_PATH, header=None, sep='  ')\n",
    "print('Shape:', ecg_df.shape)\n",
    "ecg_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b2664cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['lines.linewidth'] = 0.7\n",
    "fig = plt.figure(figsize=(5,2))\n",
    "label_normal = False\n",
    "label_ischemia = False\n",
    "for i in range(0,100):\n",
    "    label = ecg_df.iloc[i, 0]\n",
    "    if (label == -1):\n",
    "        color = colors[1]\n",
    "        \n",
    "        if label_ischemia:\n",
    "            plt.plot(ecg_df.iloc[i,1:96], color=color, alpha=0.5, linestyle='--', linewidth=0.5)\n",
    "        else:\n",
    "            plt.plot(ecg_df.iloc[i,1:96], color=color, alpha=0.5, label='Ischemia', linestyle='--')\n",
    "            label_ischemia = True\n",
    "            \n",
    "    else:\n",
    "        color = colors[0]\n",
    "        \n",
    "        if label_normal:\n",
    "            plt.plot(ecg_df.iloc[i,1:96], color=color, alpha=0.5)\n",
    "        else:\n",
    "            plt.plot(ecg_df.iloc[i,1:96], color=color, alpha=0.5, label='Normal')\n",
    "            label_normal = True\n",
    "    \n",
    "plt.title('Human heartbeat activity')\n",
    "plt.legend(loc='upper right', ncol=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b904dba",
   "metadata": {},
   "source": [
    "### **4.** Industrial pump data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce940bbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_PATH = os.path.join(DATA, 'pump', 'sensor.csv')\n",
    "pump_df = pd.read_csv(FILE_PATH, sep=',')\n",
    "pump_df.drop(columns={'Unnamed: 0'}, inplace=True)\n",
    "pump_df['timestamp'] = pd.to_datetime(pump_df['timestamp'], format='%Y-%m-%d %H:%M:%S')\n",
    "pump_df = pump_df.set_index('timestamp')\n",
    "\n",
    "pump_df['machine_status'].replace(to_replace='NORMAL', value=np.nan, inplace=True)\n",
    "pump_df['machine_status'].replace(to_replace='BROKEN', value=1, inplace=True)\n",
    "pump_df['machine_status'].replace(to_replace='RECOVERING', value=1, inplace=True)\n",
    "\n",
    "print('Shape:', pump_df.shape)\n",
    "pump_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7455aa11",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_structure_df = pump_df.iloc[:, 0:10].resample('5D').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "913b88c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['hatch.linewidth'] = 0.5\n",
    "plt.rcParams['lines.linewidth'] = 0.5\n",
    "\n",
    "fig = plt.figure(figsize=(5,1))\n",
    "ax1 = fig.add_subplot(1,1,1)\n",
    "plot1 = ax1.plot(pump_df['sensor_00'], label='Healthy pump')\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "plot2 = ax2.fill_between(\n",
    "    x=pump_df.index, \n",
    "    y1=0.0, \n",
    "    y2=pump_df['machine_status'], \n",
    "    color=colors[1], \n",
    "    linewidth=0.0,\n",
    "    edgecolor='#000000',\n",
    "    alpha=0.5, \n",
    "    hatch=\"//////\", \n",
    "    label='Broken pump'\n",
    ")\n",
    "ax2.grid(False)\n",
    "ax2.set_yticks([])\n",
    "\n",
    "labels = [plot1[0].get_label(), plot2.get_label()]\n",
    "\n",
    "plt.legend(handles=[plot1[0], plot2], labels=labels, loc='lower center', ncol=2, bbox_to_anchor=(0.5, -.4))\n",
    "plt.title('Industrial pump sensor data')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3f3743d",
   "metadata": {},
   "source": [
    "### **5.** London household energy consumption with weather data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcd0191e",
   "metadata": {},
   "source": [
    "We want to filter out households that are are subject to the dToU tariff and keep only the ones with a known ACORN (i.e. not in the ACORN-U group): this will allow us to better model future analysis by adding the Acorn detail informations (which by definitions, won't be available for the ACORN-U group)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f57ceb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "household_filename = os.path.join(DATA, 'energy-london', 'informations_households.csv')\n",
    "household_df = pd.read_csv(household_filename)\n",
    "household_df = household_df[(household_df['stdorToU'] == 'Std') & (household_df['Acorn'] == 'ACORN-E')]\n",
    "print(household_df.shape)\n",
    "household_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68add26c",
   "metadata": {},
   "source": [
    "#### Associating households with they energy consumption data\n",
    "Each household (with an ID starting by `MACxxxxx` in the table above) has its consumption data stored in a block file name `block_xx`. This file is also available from the `informations_household.csv` file extracted above. We have the association between `household_id` and `block_file`: we can open each of them and keep the consumption for the households of interest. All these data will be concatenated into an `energy_df` dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fd6fd80",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "household_ids = household_df['LCLid'].tolist()\n",
    "consumption_file = os.path.join(DATA, 'energy-london', 'hourly_consumption.csv')\n",
    "min_data_points = ((pd.to_datetime('2020-12-31') - pd.to_datetime('2020-01-01')).days + 1)*24*2\n",
    "\n",
    "if os.path.exists(consumption_file):\n",
    "    print('Half-hourly consumption file already exists, loading from disk...')\n",
    "    energy_df = pd.read_csv(consumption_file)\n",
    "    energy_df['timestamp'] = pd.to_datetime(energy_df['timestamp'], format='%Y-%m-%d %H:%M:%S.%f')\n",
    "    print('Done.')\n",
    "    \n",
    "else:\n",
    "    print('Half-hourly consumption file not found. We need to generate it.')\n",
    "    \n",
    "    # We know have the block number we can use to open the right file:\n",
    "    energy_df = pd.DataFrame()\n",
    "    target_block_files = household_df['file'].unique().tolist()\n",
    "    print('- {} block files to process: '.format(len(target_block_files)), end='')\n",
    "    df_list = []\n",
    "    for block_file in tqdm(target_block_files):\n",
    "        # Reads the current block file:\n",
    "        current_filename = os.path.join(DATA, 'energy-london', 'halfhourly_dataset', '{}.csv'.format(block_file))\n",
    "        df = pd.read_csv(current_filename)\n",
    "        \n",
    "        # Set readable column names and adjust data types:\n",
    "        df.columns = ['household_id', 'timestamp', 'energy']\n",
    "        df = df.replace(to_replace='Null', value=0.0)\n",
    "        df['energy'] = df['energy'].astype(np.float64)\n",
    "        df['timestamp'] = pd.to_datetime(df['timestamp'], format='%Y-%m-%d %H:%M:%S.%f')\n",
    "        \n",
    "        # We filter on the households sampled earlier:\n",
    "        df_list.append(df[df['household_id'].isin(household_ids)].reset_index(drop=True))\n",
    "    \n",
    "    # Concatenate with the main dataframe:\n",
    "    energy_df = pd.concat(df_list, axis='index', ignore_index=True)\n",
    "    \n",
    "    datapoints = energy_df.groupby(by='household_id').count()\n",
    "    datapoints = datapoints[datapoints['timestamp'] < min_data_points]\n",
    "    hhid_to_remove = datapoints.index.tolist()\n",
    "    energy_df = energy_df[~energy_df['household_id'].isin(hhid_to_remove)]\n",
    "\n",
    "    # Let's save this dataset to disk, we will use it from now on:\n",
    "    print('Saving file to disk... ', end='')\n",
    "    energy_df.to_csv(consumption_file, index=False)\n",
    "    print('Done.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "465cdc8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "start = np.min(energy_df['timestamp'])\n",
    "end = np.max(energy_df['timestamp'])\n",
    "weather_filename = os.path.join(DATA, 'energy-london', 'weather_hourly_darksky.csv')\n",
    "\n",
    "weather_df = pd.read_csv(weather_filename)\n",
    "weather_df['time'] = pd.to_datetime(weather_df['time'], format='%Y-%m-%d %H:%M:%S')\n",
    "weather_df = weather_df.drop(columns=['precipType', 'icon', 'summary'])\n",
    "weather_df = weather_df.sort_values(by='time')\n",
    "weather_df = weather_df.set_index('time')\n",
    "weather_df = weather_df[start:end]\n",
    "\n",
    "# Let's make sure we have one datapoint per hour to match \n",
    "# the frequency used for the household energy consumption data:\n",
    "weather_df = weather_df.resample(rule='1H').mean()     # This will generate NaN values timestamp missing data\n",
    "weather_df = weather_df.interpolate(method='linear')   # This will fill the missing values with the average \n",
    "\n",
    "print(weather_df.shape)\n",
    "weather_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3427cb1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_df = energy_df.set_index(['household_id', 'timestamp'])\n",
    "energy_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f173f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "hhid = household_ids[2]\n",
    "hh_energy = energy_df.loc[hhid, :]\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "\n",
    "fig = plt.figure(figsize=(5,1))\n",
    "ax1 = fig.add_subplot(1,1,1)\n",
    "plot2 = ax1.fill_between(\n",
    "    x=weather_df.loc[start:end, 'temperature'].index, \n",
    "    y1=0.0, \n",
    "    y2=weather_df.loc[start:end, 'temperature'], \n",
    "    color=colors[1], \n",
    "    linewidth=0.0,\n",
    "    edgecolor='#000000',\n",
    "    alpha=0.25, \n",
    "    hatch=\"//////\", \n",
    "    label='Temperature'\n",
    ")\n",
    "ax1.set_ylim((0,40))\n",
    "ax1.grid(False)\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(hh_energy[start:end], label='Energy consumption', linewidth=2, color='#FFFFFF', alpha=0.5)\n",
    "plot1 = ax2.plot(hh_energy[start:end], label='Energy consumption', linewidth=0.7)\n",
    "ax2.set_title(f'Energy consumption for household {hhid}')\n",
    "\n",
    "labels = [plot1[0].get_label(), plot2.get_label()]\n",
    "plt.legend(handles=[plot1[0], plot2], labels=labels, loc='upper left', fontsize=3, ncol=2)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36684f21",
   "metadata": {},
   "outputs": [],
   "source": [
    "acorn_filename = os.path.join(DATA, 'energy-london', 'acorn_details.csv')\n",
    "acorn_df = pd.read_csv(acorn_filename, encoding='ISO-8859-1')\n",
    "acorn_df = acorn_df.sample(10).loc[:, ['MAIN CATEGORIES', 'CATEGORIES', 'REFERENCE', 'ACORN-A', 'ACORN-B', 'ACORN-E']]\n",
    "acorn_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f09fcaf8",
   "metadata": {},
   "source": [
    "## File structure exploration\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d87ca63a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display_html\n",
    "\n",
    "def display_multiple_dataframe(*args, max_rows=None, max_cols=None):\n",
    "    html_str = ''\n",
    "    for df in args:\n",
    "        html_str += df.to_html(max_cols=max_cols, max_rows=max_rows)\n",
    "        \n",
    "    display_html(html_str.replace('table','table style=\"display:inline\"'), raw=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0fb229a",
   "metadata": {},
   "outputs": [],
   "source": [
    "display_multiple_dataframe(\n",
    "    file_structure_df[['sensor_00']],\n",
    "    file_structure_df[['sensor_01']],\n",
    "    file_structure_df[['sensor_03']],\n",
    "    max_rows=10, max_cols=None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ca40f25",
   "metadata": {},
   "outputs": [],
   "source": [
    "display_multiple_dataframe(\n",
    "    file_structure_df.loc['2018-04', :].head(6),\n",
    "    file_structure_df.loc['2018-05', :].head(6),\n",
    "    file_structure_df.loc['2018-06', :].head(6),\n",
    "    max_rows=None, max_cols=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e87d94c",
   "metadata": {},
   "outputs": [],
   "source": [
    "display_multiple_dataframe(\n",
    "    file_structure_df.loc['2018-04', ['sensor_00']].head(6),\n",
    "    file_structure_df.loc['2018-05', ['sensor_00']].head(6),\n",
    "    file_structure_df.loc['2018-06', ['sensor_00']].head(6),\n",
    "    max_rows=10, max_cols=None\n",
    ")\n",
    "display_multiple_dataframe(\n",
    "    file_structure_df.loc['2018-04', ['sensor_01']].head(6),\n",
    "    file_structure_df.loc['2018-05', ['sensor_01']].head(6),\n",
    "    file_structure_df.loc['2018-06', ['sensor_01']].head(6),\n",
    "    max_rows=10, max_cols=None\n",
    ")\n",
    "print('.\\n.\\n.')\n",
    "display_multiple_dataframe(\n",
    "    file_structure_df.loc['2018-04', ['sensor_09']].head(6),\n",
    "    file_structure_df.loc['2018-05', ['sensor_09']].head(6),\n",
    "    file_structure_df.loc['2018-06', ['sensor_09']].head(6),\n",
    "    max_rows=10, max_cols=None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c8d815f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pump_df.iloc[:, [0]].resample('5D').mean()\n",
    "df2 = pump_df.iloc[:, [1]].resample('2D').mean()\n",
    "df3 = pump_df.iloc[:, [2]].resample('7D').mean()\n",
    "\n",
    "display_multiple_dataframe(\n",
    "    df1.head(10), df2.head(10), df3.head(10),\n",
    "    pd.merge(pd.merge(df1, df2, left_index=True, right_index=True, how='outer'), df3, left_index=True, right_index=True, how='outer').head(10),\n",
    "    max_rows=None, max_cols=None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7046cbb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', None)\n",
    "pd.set_option('display.max_rows', 10)\n",
    "pd.merge(pd.merge(df1, df2, left_index=True, right_index=True, how='outer'), df3, left_index=True, right_index=True, how='outer').head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ad04201",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,1))\n",
    "for i in range(len(colors)):\n",
    "    plt.plot(file_structure_df[f'sensor_0{i}'], linewidth=2, alpha=0.5, label=colors[i])\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7b85f1f",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4c78d4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5,1))\n",
    "ax1 = fig.add_subplot(1,1,1)\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "plot_sensor_0 = ax1.plot(pump_df['sensor_00'], label='Sensor 0', color=colors[0], linewidth=1, alpha=0.8)\n",
    "plot_sensor_1 = ax2.plot(pump_df['sensor_01'], label='Sensor 1', color=colors[1], linewidth=1, alpha=0.8)\n",
    "ax2.grid(False)\n",
    "plt.title('Pump sensor values (2 sensors)')\n",
    "plt.legend(handles=[plot_sensor_0[0], plot_sensor_1[0]], ncol=2, loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90daaba9",
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_pump_df = pump_df.loc[:, 'sensor_00':'sensor_14']\n",
    "reduced_pump_df = reduced_pump_df.replace([np.inf, -np.inf], np.nan)\n",
    "reduced_pump_df = reduced_pump_df.fillna(0.0)\n",
    "reduced_pump_df = reduced_pump_df.astype(np.float32)\n",
    "scaled_pump_df = pd.DataFrame(normalize(reduced_pump_df), index=reduced_pump_df.index, columns=reduced_pump_df.columns)\n",
    "scaled_pump_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbc92f9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5,1))\n",
    "\n",
    "for i in range(0,15):\n",
    "    plt.plot(scaled_pump_df.iloc[:, i], alpha=0.6)\n",
    "\n",
    "plt.title('Pump sensor values (15 sensors)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9739c28e",
   "metadata": {},
   "outputs": [],
   "source": [
    "pump_df2 = pump_df.copy()\n",
    "\n",
    "pump_df2 = pump_df2.replace([np.inf, -np.inf], np.nan)\n",
    "pump_df2 = pump_df2.fillna(0.0)\n",
    "pump_df2 = pump_df2.astype(np.float32)\n",
    "\n",
    "pump_description = pump_df2.describe().T\n",
    "constant_signals = pump_description[pump_description['min'] == pump_description['max']].index.tolist()\n",
    "pump_df2 = pump_df2.drop(columns=constant_signals)\n",
    "\n",
    "features = pump_df2.columns.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd1cf267",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hex_to_rgb(hex_color):\n",
    "    \"\"\"\n",
    "    Converts a color string in hexadecimal format to RGB format.\n",
    "    \n",
    "    PARAMS\n",
    "    ======\n",
    "        hex_color: string\n",
    "            A string describing the color to convert from hexadecimal. It can\n",
    "            include the leading # character or not\n",
    "    \n",
    "    RETURNS\n",
    "    =======\n",
    "        rgb_color: tuple\n",
    "            Each color component of the returned tuple will be a float value\n",
    "            between 0.0 and 1.0\n",
    "    \"\"\"\n",
    "    hex_color = hex_color.lstrip('#')\n",
    "    rgb_color = tuple(int(hex_color[i:i+2], base=16) / 255.0 for i in [0, 2, 4])\n",
    "    return rgb_color\n",
    "\n",
    "def plot_timeseries_strip_chart(binned_timeseries, signal_list, fig_width=12, signal_height=0.15, dates=None, day_interval=7):\n",
    "    # Build a suitable colormap:\n",
    "    colors_list = [\n",
    "        hex_to_rgb('#DC322F'), \n",
    "        hex_to_rgb('#B58900'), \n",
    "        hex_to_rgb('#2AA198')\n",
    "    ]\n",
    "    cm = mpl_colors.LinearSegmentedColormap.from_list('RdAmGr', colors_list, N=len(colors_list))\n",
    "    \n",
    "    fig = plt.figure(figsize=(fig_width, signal_height * binned_timeseries.shape[0]))\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    \n",
    "    # Devising the extent of the actual plot:\n",
    "    if dates is not None:\n",
    "        dnum = mdates.date2num(dates)\n",
    "        start = dnum[0] - (dnum[1]-dnum[0])/2.\n",
    "        stop = dnum[-1] + (dnum[1]-dnum[0])/2.\n",
    "        extent = [start, stop, 0, signal_height * (binned_timeseries.shape[0])]\n",
    "        \n",
    "    else:\n",
    "        extent = None\n",
    "        \n",
    "    # Plot the matrix:\n",
    "    im = ax.imshow(binned_timeseries, \n",
    "                   extent=extent, \n",
    "                   aspect=\"auto\", \n",
    "                   cmap=cm, \n",
    "                   origin='lower')\n",
    "    \n",
    "    # Adjusting the x-axis if we provide dates:\n",
    "    if dates is not None:\n",
    "        ax.xaxis.set_major_locator(mdates.MonthLocator())\n",
    "        ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))\n",
    "        for tick in ax.xaxis.get_major_ticks():\n",
    "            tick.label.set_fontsize(4)\n",
    "            tick.label.set_rotation(60)\n",
    "            tick.label.set_fontweight('bold')\n",
    "\n",
    "        ax.tick_params(axis='x', which='major', pad=7, labelcolor='#000000')\n",
    "        plt.xticks(ha='right')\n",
    "        \n",
    "    # Adjusting the y-axis:\n",
    "    ax.yaxis.set_major_locator(ticker.MultipleLocator(signal_height))\n",
    "    ax.set_yticklabels(signal_list, verticalalignment='bottom', fontsize=4)\n",
    "    ax.set_yticks(np.arange(len(signal_list)) * signal_height)\n",
    "\n",
    "    plt.grid()\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef704a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, Markdown, Latex\n",
    "\n",
    "# Build a list of dataframes, one per sensor:\n",
    "df_list = []\n",
    "for f in features[:1]:\n",
    "    df_list.append(pump_df2[[f]])\n",
    "\n",
    "# Discretize each signal in 3 bins:\n",
    "array = tsia.markov.discretize_multivariate(df_list)\n",
    "\n",
    "fig = plt.figure(figsize=(5.5, 0.6))\n",
    "plt.plot(pump_df2['sensor_00'], linewidth=0.7, alpha=0.6)\n",
    "plt.title('Line plot of the pump sensor 0')\n",
    "plt.show()\n",
    "\n",
    "display(Markdown('<img src=\"arrow.png\" align=\"left\" style=\"padding-left: 730px\"/>'))\n",
    "\n",
    "\n",
    "# Plot the strip chart:\n",
    "ax = plot_timeseries_strip_chart(\n",
    "    array, \n",
    "    signal_list=features[:1],\n",
    "    fig_width=5.21,\n",
    "    signal_height=0.2,\n",
    "    dates=df_list[0].index.to_pydatetime(),\n",
    "    day_interval=2\n",
    ")\n",
    "ax.set_title('Strip chart of the pump sensor 0');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59c2ff60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a list of dataframes, one per sensor:\n",
    "df_list = []\n",
    "for f in features:\n",
    "    df_list.append(pump_df2[[f]])\n",
    "\n",
    "# Discretize each signal in 3 bins:\n",
    "array = tsia.markov.discretize_multivariate(df_list)\n",
    "\n",
    "# Plot the strip chart:\n",
    "fig = plot_timeseries_strip_chart(\n",
    "    array, \n",
    "    signal_list=features,\n",
    "    fig_width=5.5,\n",
    "    signal_height=0.1,\n",
    "    dates=df_list[0].index.to_pydatetime(),\n",
    "    day_interval=2\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa605700",
   "metadata": {},
   "source": [
    "### Recurrence plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f85db253",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyts.image import RecurrencePlot\n",
    "from pyts.image import GramianAngularField\n",
    "from pyts.image import MarkovTransitionField"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74fb9a2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "hhid = household_ids[2]\n",
    "hh_energy = energy_df.loc[hhid, :]\n",
    "pump_extract_df = pump_df.iloc[:800, 0].copy()\n",
    "\n",
    "rp = RecurrencePlot(threshold='point', percentage=30)\n",
    "weather_rp = rp.fit_transform(weather_df.loc['2013-01-01':'2013-01-31']['temperature'].values.reshape(1, -1))\n",
    "energy_rp = rp.fit_transform(hh_energy['2012-07-01':'2012-07-15'].values.reshape(1, -1))\n",
    "pump_rp = rp.fit_transform(pump_extract_df.values.reshape(1, -1))\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(5.5, 2.4))\n",
    "gs = gridspec.GridSpec(nrows=3, ncols=2, width_ratios=[3,1], hspace=0.8, wspace=0.0)\n",
    "\n",
    "# Pump sensor 0:\n",
    "ax = fig.add_subplot(gs[0])\n",
    "ax.plot(pump_extract_df, label='Pump sensor 0')\n",
    "ax.set_title(f'Pump sensor 0')\n",
    "\n",
    "ax = fig.add_subplot(gs[1])\n",
    "ax.imshow(pump_rp[0], cmap='binary', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "# Energy consumption line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[2])\n",
    "plot1 = ax.plot(hh_energy['2012-07-01':'2012-07-15'], color=colors[1])\n",
    "ax.set_title(f'Energy consumption for household {hhid}')\n",
    "\n",
    "ax = fig.add_subplot(gs[3])\n",
    "ax.imshow(energy_rp[0], cmap='binary', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "# Daily temperature line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[4])\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "ax.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color=colors[2])\n",
    "ax.set_title(f'Daily temperature')\n",
    "\n",
    "ax = fig.add_subplot(gs[5])\n",
    "ax.imshow(weather_rp[0], cmap='binary', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "950062ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "hhid = household_ids[2]\n",
    "hh_energy = energy_df.loc[hhid, :]\n",
    "pump_extract_df = pump_df.iloc[:800, 0].copy()\n",
    "\n",
    "gaf = GramianAngularField(image_size=48, method='summation')\n",
    "weather_gasf = gaf.fit_transform(weather_df.loc['2013-01-01':'2013-01-31']['temperature'].values.reshape(1, -1))\n",
    "energy_gasf = gaf.fit_transform(hh_energy['2012-07-01':'2012-07-15'].values.reshape(1, -1))\n",
    "pump_gasf = gaf.fit_transform(pump_extract_df.values.reshape(1, -1))\n",
    "\n",
    "fig = plt.figure(figsize=(5.5, 2.4))\n",
    "gs = gridspec.GridSpec(nrows=3, ncols=2, width_ratios=[3,1], hspace=0.8, wspace=0.0)\n",
    "\n",
    "# Pump sensor 0:\n",
    "ax = fig.add_subplot(gs[0])\n",
    "ax.plot(pump_extract_df, label='Pump sensor 0')\n",
    "ax.set_title(f'Pump sensor 0')\n",
    "\n",
    "ax = fig.add_subplot(gs[1])\n",
    "ax.imshow(pump_gasf[0], cmap='RdBu_r', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "# Energy consumption line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[2])\n",
    "plot1 = ax.plot(hh_energy['2012-07-01':'2012-07-15'], color=colors[1])\n",
    "ax.set_title(f'Energy consumption for household {hhid}')\n",
    "\n",
    "ax = fig.add_subplot(gs[3])\n",
    "ax.imshow(energy_gasf[0], cmap='RdBu_r', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "# Daily temperature line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[4])\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "ax.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color=colors[2])\n",
    "ax.set_title(f'Daily temperature')\n",
    "\n",
    "ax = fig.add_subplot(gs[5])\n",
    "ax.imshow(weather_gasf[0], cmap='RdBu_r', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "123a169c",
   "metadata": {},
   "outputs": [],
   "source": [
    "mtf = MarkovTransitionField(image_size=48)\n",
    "\n",
    "weather_mtf = mtf.fit_transform(weather_df.loc['2013-01-01':'2013-01-31']['temperature'].values.reshape(1, -1))\n",
    "energy_mtf = mtf.fit_transform(hh_energy['2012-07-01':'2012-07-15'].values.reshape(1, -1))\n",
    "pump_mtf = mtf.fit_transform(pump_extract_df.values.reshape(1, -1))\n",
    "\n",
    "fig = plt.figure(figsize=(5.5, 2.4))\n",
    "gs = gridspec.GridSpec(nrows=3, ncols=2, width_ratios=[3,1], hspace=0.8, wspace=0.0)\n",
    "\n",
    "# Pump sensor 0:\n",
    "ax = fig.add_subplot(gs[0])\n",
    "ax.plot(pump_extract_df, label='Pump sensor 0')\n",
    "ax.set_title(f'Pump sensor 0')\n",
    "\n",
    "ax = fig.add_subplot(gs[1])\n",
    "ax.imshow(pump_mtf[0], cmap='RdBu_r', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "# Energy consumption line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[2])\n",
    "plot1 = ax.plot(hh_energy['2012-07-01':'2012-07-15'], color=colors[1])\n",
    "ax.set_title(f'Energy consumption for household {hhid}')\n",
    "\n",
    "ax = fig.add_subplot(gs[3])\n",
    "ax.imshow(energy_mtf[0], cmap='RdBu_r', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "# Daily temperature line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[4])\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "ax.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color=colors[2])\n",
    "ax.set_title(f'Daily temperature')\n",
    "\n",
    "ax = fig.add_subplot(gs[5])\n",
    "ax.imshow(weather_mtf[0], cmap='RdBu_r', origin='lower')\n",
    "ax.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f4e8ebe",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.cm as cm\n",
    "import networkx as nx\n",
    "import community\n",
    "\n",
    "def compute_network_graph(markov_field):\n",
    "    G = nx.from_numpy_matrix(markov_field[0])\n",
    "\n",
    "    # Uncover the communities in the current graph:\n",
    "    communities = community.best_partition(G)\n",
    "    nb_communities = len(pd.Series(communities).unique())\n",
    "    cmap = 'autumn'\n",
    "\n",
    "    # Compute node colors and edges colors for the modularity encoding:\n",
    "    edge_colors = [matplotlib.colors.to_hex(cm.get_cmap(cmap)(communities.get(v)/(nb_communities - 1))) for u,v in G.edges()]\n",
    "    node_colors = [communities.get(node) for node in G.nodes()]\n",
    "    node_size = [nx.average_clustering(G, [node])*90 for node in G.nodes()]\n",
    "\n",
    "    # Builds the options set to draw the network graph in the \"modularity\" configuration:\n",
    "    options = {\n",
    "        'node_size': 10,\n",
    "        'edge_color': edge_colors,\n",
    "        'node_color': node_colors,\n",
    "        'linewidths': 0,\n",
    "        'width': 0.1,\n",
    "        'alpha': 0.6,\n",
    "        'with_labels': False,\n",
    "        'cmap': cmap\n",
    "    }\n",
    "    \n",
    "    return G, options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bfa993a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5.5, 2.4))\n",
    "gs = gridspec.GridSpec(nrows=3, ncols=2, width_ratios=[3,1], hspace=0.8, wspace=0.0)\n",
    "\n",
    "# Pump sensor 0:\n",
    "ax = fig.add_subplot(gs[0])\n",
    "ax.plot(pump_extract_df, label='Pump sensor 0')\n",
    "ax.set_title(f'Pump sensor 0')\n",
    "\n",
    "ax = fig.add_subplot(gs[1])\n",
    "G, options = compute_network_graph(weather_mtf)\n",
    "nx.draw_networkx(G, **options, pos=nx.spring_layout(G), ax=ax)\n",
    "ax.axis('off')\n",
    "\n",
    "# Energy consumption line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[2])\n",
    "plot1 = ax.plot(hh_energy['2012-07-01':'2012-07-15'], color=colors[1])\n",
    "ax.set_title(f'Energy consumption for household {hhid}')\n",
    "\n",
    "ax = fig.add_subplot(gs[3])\n",
    "G, options = compute_network_graph(energy_mtf)\n",
    "nx.draw_networkx(G, **options, pos=nx.spring_layout(G), ax=ax)\n",
    "ax.axis('off')\n",
    "\n",
    "# Daily temperature line plot and recurrence plot:\n",
    "ax = fig.add_subplot(gs[4])\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "ax.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color=colors[2])\n",
    "ax.set_title(f'Daily temperature')\n",
    "\n",
    "ax = fig.add_subplot(gs[5])\n",
    "G, options = compute_network_graph(weather_mtf)\n",
    "nx.draw_networkx(G, **options, pos=nx.spring_layout(G), ax=ax)\n",
    "ax.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "979ad8ea",
   "metadata": {},
   "source": [
    "## Symbolic representation\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "102d6a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyts.bag_of_words import BagOfWords\n",
    "\n",
    "window_size, word_size = 30, 5\n",
    "bow = BagOfWords(window_size=window_size, word_size=word_size, window_step=window_size, numerosity_reduction=False)\n",
    "X = weather_df.loc['2013-01-01':'2013-01-31']['temperature'].values.reshape(1, -1)\n",
    "X_bow = bow.transform(X)\n",
    "time_index = weather_df.loc['2013-01-01':'2013-01-31']['temperature'].index\n",
    "len(X_bow[0].replace(' ', ''))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96878a1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the considered subseries\n",
    "plt.figure(figsize=(5, 2))\n",
    "splits_series = np.linspace(0, X.shape[1], 1 + X.shape[1] // window_size, dtype='int64')\n",
    "for start, end in zip(splits_series[:-1], np.clip(splits_series[1:] + 1, 0, X.shape[1])):\n",
    "    plt.plot(np.arange(start, end), X[0, start:end], 'o-', linewidth=0.5, ms=0.1)\n",
    "\n",
    "# Plot the corresponding letters\n",
    "splits_letters = np.linspace(0, X.shape[1], 1 + word_size * X.shape[1] // window_size)\n",
    "splits_letters = ((splits_letters[:-1] + splits_letters[1:]) / 2)\n",
    "splits_letters = splits_letters.astype('int64')\n",
    "\n",
    "for i, (x, text) in enumerate(zip(splits_letters, X_bow[0].replace(' ', ''))):\n",
    "    t = plt.text(x, X[0, x], text, color=\"C{}\".format(i // 5), fontsize=3.5)\n",
    "    t.set_bbox(dict(facecolor='#FFFFFF', alpha=0.5, edgecolor=\"C{}\".format(i // 5), boxstyle='round4'))\n",
    "\n",
    "plt.title('Bag-of-words representation for weather temperature')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94bb8ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyts.transformation import WEASEL\n",
    "from sklearn.preprocessing import LabelEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d501fe5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = ecg_df.iloc[:, 1:].values\n",
    "y_train = ecg_df.iloc[:, 0]\n",
    "y_train = LabelEncoder().fit_transform(y_train)\n",
    "weasel = WEASEL(word_size=3, n_bins=3, window_sizes=[10, 25], sparse=False)\n",
    "X_weasel = weasel.fit_transform(X_train, y_train)\n",
    "vocabulary_length = len(weasel.vocabulary_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eef27d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,1.5))\n",
    "width = 0.4\n",
    "x = np.arange(vocabulary_length) - width / 2\n",
    "for i in range(len(X_weasel[y_train == 0])):\n",
    "    if i == 0:\n",
    "        plt.bar(x, X_weasel[y_train == 0][i], width=width, alpha=0.25, color=colors[1], label='Time series for Ischemia')\n",
    "    else:\n",
    "        plt.bar(x, X_weasel[y_train == 0][i], width=width, alpha=0.25, color=colors[1])\n",
    "    \n",
    "for i in range(len(X_weasel[y_train == 1])):\n",
    "    if i == 0:\n",
    "        plt.bar(x+width, X_weasel[y_train == 1][i], width=width, alpha=0.25, color=colors[0], label='Time series for Normal heartbeat')\n",
    "    else:\n",
    "        plt.bar(x+width, X_weasel[y_train == 1][i], width=width, alpha=0.25, color=colors[0])\n",
    "        \n",
    "plt.xticks(\n",
    "    np.arange(vocabulary_length),\n",
    "    np.vectorize(weasel.vocabulary_.get)(np.arange(X_weasel[0].size)),\n",
    "    fontsize=2,\n",
    "    rotation=60\n",
    ")\n",
    "    \n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6677dddd",
   "metadata": {},
   "source": [
    "## Statistics\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c2ba909",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['xtick.labelsize'] = 3\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "fig = plt.figure(figsize=(5.5, 3))\n",
    "gs = gridspec.GridSpec(nrows=3, ncols=2, width_ratios=[1,1], hspace=0.8)\n",
    "\n",
    "# Pump\n",
    "ax = fig.add_subplot(gs[0])\n",
    "ax.plot(pump_extract_df, label='Pump sensor 0')\n",
    "ax.set_title(f'Pump sensor 0')\n",
    "ax.tick_params(axis='x', which='both', labelbottom=False)\n",
    "\n",
    "ax = fig.add_subplot(gs[1])\n",
    "sm.graphics.tsa.plot_acf(pump_extract_df.values.squeeze(), ax=ax, markersize=1, title='')\n",
    "ax.set_ylim(-1.2, 1.2)\n",
    "ax.tick_params(axis='x', which='major', labelsize=4)\n",
    "\n",
    "# Energy consumption\n",
    "ax = fig.add_subplot(gs[2])\n",
    "ax.plot(hh_energy['2012-07-01':'2012-07-15'], color=colors[1])\n",
    "ax.set_title(f'Energy consumption for household {hhid}')\n",
    "ax.tick_params(axis='x', which='both', labelbottom=False)\n",
    "\n",
    "ax = fig.add_subplot(gs[3])\n",
    "sm.graphics.tsa.plot_acf(hh_energy['2012-07-01':'2012-07-15'].values.squeeze(), ax=ax, markersize=1, title='')\n",
    "ax.set_ylim(-0.3, 0.3)\n",
    "ax.tick_params(axis='x', which='major', labelsize=4)\n",
    "\n",
    "# Daily temperature:\n",
    "ax = fig.add_subplot(gs[4])\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "ax.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color=colors[2])\n",
    "ax.set_title(f'Daily temperature')\n",
    "ax.tick_params(axis='x', which='both', labelbottom=False)\n",
    "\n",
    "ax = fig.add_subplot(gs[5])\n",
    "sm.graphics.tsa.plot_acf(weather_df.loc['2013-01-01':'2013-01-31']['temperature'].values.squeeze(), ax=ax, markersize=1, title='')\n",
    "ax.set_ylim(-1.2, 1.2)\n",
    "ax.tick_params(axis='x', which='major', labelsize=4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f9cd561",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import STL\n",
    "\n",
    "endog = endog.resample('30T').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2af017b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['lines.markersize'] = 1\n",
    "\n",
    "title = f'Energy consumption for household {hhid}'\n",
    "endog = hh_energy['2012-07-01':'2012-07-15']\n",
    "endog.columns = [title]\n",
    "endog = endog[title]\n",
    "stl = STL(endog, period=48)\n",
    "res = stl.fit()\n",
    "fig = res.plot()\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(5.5, 4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebb389e7",
   "metadata": {},
   "source": [
    "## Binary segmentation\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4494d2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = weather_df.loc['2013-01-01':'2013-01-31']['temperature'].values.squeeze()\n",
    "algo = rpt.Binseg(model='l2').fit(signal)\n",
    "my_bkps = algo.predict(n_bkps=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d13ccee7",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_bkps = [0] + my_bkps\n",
    "my_bkps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6379306a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5.5,1))\n",
    "start = '2012-07-01'\n",
    "end = '2012-07-15'\n",
    "plt.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color='#FFFFFF', linewidth=1.2, alpha=0.8)\n",
    "plt.plot(weather_df.loc['2013-01-01':'2013-01-31']['temperature'], color=colors[2], linewidth=0.7)\n",
    "\n",
    "plt.title(f'Daily temperature')\n",
    "plt.xticks(rotation=60, fontsize=4)\n",
    "\n",
    "weather_index = weather_df.loc['2013-01-01':'2013-01-31']['temperature'].index\n",
    "\n",
    "for index, bkps in enumerate(my_bkps[:-1]):\n",
    "    x1 = weather_index[my_bkps[index]]\n",
    "    x2 = weather_index[np.clip(my_bkps[index+1], 0, len(weather_index)-1)]\n",
    "    \n",
    "    plt.axvspan(x1, x2, color=colors[index % 5], alpha=0.2)\n",
    "\n",
    "plt.title('Daily temperature segmentation')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "conda_python3",
   "language": "python",
   "name": "conda_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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}