{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Univariate Time Series Anomaly Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## boilerplate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "import copy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from pylab import rcParams\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix\n",
    "\n",
    "from torch import nn, optim\n",
    "\n",
    "import torch.nn.functional as F\n",
    "import random\n",
    "import datetime\n",
    "# from arff2pandas import a2p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<torch._C.Generator at 0x7ff27834c5b0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "sns.set(style='whitegrid', palette='muted', font_scale=1.2)\n",
    "\n",
    "HAPPY_COLORS_PALETTE = [\"#01BEFE\", \"#FFDD00\", \"#FF7D00\", \"#FF006D\", \"#ADFF02\", \"#8F00FF\"]\n",
    "\n",
    "sns.set_palette(sns.color_palette(HAPPY_COLORS_PALETTE))\n",
    "\n",
    "# rcParams['figure.figsize'] = 12, 8\n",
    "\n",
    "RANDOM_SEED = 42\n",
    "np.random.seed(RANDOM_SEED)\n",
    "torch.manual_seed(RANDOM_SEED) \n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import polars as pl\n",
    "from io import StringIO\n",
    "import math\n",
    "df = pl.read_csv('../data/battery_1.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We only need 'PACK1_CRIDATA_BATT_VOL'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualize fault and non-fault regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_condition =  df['PACK1_CRIDATA_BATT_VOL'].cast(pl.Float32) != 0\n",
    "voltage_data = (df['PACK1_CRIDATA_BATT_VOL']\n",
    "                .filter(filter_condition)\n",
    "                .cast(pl.Float32))\n",
    "\n",
    "def convert_values(values):\n",
    "    numerical_values = []\n",
    "    for value in values:\n",
    "        if value == 'False':\n",
    "            numerical_values.append(0)\n",
    "        elif value == 'True':\n",
    "            numerical_values.append(1)\n",
    "        else:\n",
    "            # numerical_values.append(np.nan)\n",
    "            numerical_values.append(-1)\n",
    "    return numerical_values\n",
    "\n",
    "\n",
    "fault_data = convert_values(df['BATT_PACK_1_FAULT']\n",
    "                            .filter(filter_condition))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'fault incidents')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 102568 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 537,
       "width": 852
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(10,3 * 2))\n",
    "\n",
    "axs[0].plot(voltage_data)\n",
    "axs[0].set_title(\"voltage\")\n",
    "axs[1].scatter(range(len(fault_data)), fault_data)\n",
    "axs[1].set_title(\"fault incidents\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train region: 145000-end\n",
    "\n",
    "Test region: 45000-60000"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "train_voltage_data = voltage_data[145000:]\n",
    "test_voltage_data = voltage_data[85000:100000]\n",
    "val_voltage_data = voltage_data[120000:135000]\n",
    "fault_voltage_data = voltage_data[45000:60000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_voltage_data = voltage_data[100000:]\n",
    "test_voltage_data = voltage_data[60000:85000]\n",
    "val_voltage_data = voltage_data[85000:100000]\n",
    "fault_voltage_data = voltage_data[0:60000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'fault incidents')"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 92264 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 537,
       "width": 852
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fault region\n",
    "fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(10,3 * 2))\n",
    "\n",
    "axs[0].plot(fault_voltage_data)\n",
    "axs[0].set_title(\"voltage\")\n",
    "fault_incidents = fault_data[45000:60000]\n",
    "axs[1].scatter(range(len(fault_incidents)), fault_incidents)\n",
    "axs[1].set_title(\"fault incidents\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff17067c890>]"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 81392 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 570
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# test region\n",
    "plt.plot(test_voltage_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff1704bb790>]"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 88300 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 427,
       "width": 564
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# test region\n",
    "plt.plot(val_voltage_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff17059fc90>]"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 45116 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 572
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# train region\n",
    "# plt.plot(voltage_data[145000:])\n",
    "plt.plot(train_voltage_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100774,)"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_voltage_data.to_numpy().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(25000,)"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_voltage_data.to_numpy().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = train_voltage_data.to_numpy()\n",
    "segment_size = 100\n",
    "segments = [ torch.tensor(data[i:i + segment_size]).unsqueeze(1).float() for i in range(0, len(data), segment_size) ]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (segments[-1].shape[0] != segment_size):\n",
    "    segments.pop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([557, 100, 1])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.stack(segments).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_dataset(series, segment_size, mean, stdev):\n",
    "    # normalize the data\n",
    "    data = (series.to_numpy() - mean)/stdev\n",
    "    segments = [ torch.tensor(data[i:i + segment_size]).unsqueeze(1).float() for i in range(0, len(data), segment_size) ]\n",
    "    # reject the last segment if it doesn't fit the shape\n",
    "    if (segments[-1].shape[0] != segment_size):\n",
    "        segments.pop()\n",
    "    n_seq, seq_len, n_features = torch.stack(segments).shape\n",
    "\n",
    "    return segments, seq_len, n_features\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = train_voltage_data.to_numpy()\n",
    "mean = np.mean(data)\n",
    "stdev = np.std(data)\n",
    "\n",
    "\n",
    "segment_size = 60\n",
    "train_dataset, seq_len, n_features = create_dataset(train_voltage_data, segment_size, mean, stdev)\n",
    "val_dataset, _, _ = create_dataset(val_voltage_data, segment_size, mean, stdev)\n",
    "test_normal_dataset, _, _ = create_dataset(test_voltage_data, segment_size, mean, stdev)\n",
    "test_anomaly_dataset, _, _ = create_dataset(fault_voltage_data, segment_size, mean, stdev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Encoder Decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cuda\n"
     ]
    }
   ],
   "source": [
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "print(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Encoder(nn.Module):\n",
    "\n",
    "  def __init__(self, seq_len, n_features, embedding_dim=64):\n",
    "    super(Encoder, self).__init__()\n",
    "\n",
    "    self.seq_len, self.n_features = seq_len, n_features\n",
    "    self.embedding_dim, self.hidden_dim = embedding_dim, 2 * embedding_dim\n",
    "\n",
    "    self.rnn1 = nn.LSTM(\n",
    "      input_size=n_features,\n",
    "      hidden_size=self.hidden_dim,\n",
    "      num_layers=1,\n",
    "      batch_first=True\n",
    "    )\n",
    "    \n",
    "    self.rnn2 = nn.LSTM(\n",
    "      input_size=self.hidden_dim,\n",
    "      hidden_size=embedding_dim,\n",
    "      num_layers=1,\n",
    "      batch_first=True\n",
    "    )\n",
    "\n",
    "  def forward(self, x):\n",
    "    x = x.reshape((1, self.seq_len, self.n_features))\n",
    "\n",
    "    x, (_, _) = self.rnn1(x)\n",
    "    x, (hidden_n, _) = self.rnn2(x)\n",
    "\n",
    "    return hidden_n.reshape((self.n_features, self.embedding_dim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Decoder(nn.Module):\n",
    "\n",
    "  def __init__(self, seq_len, input_dim=64, n_features=1):\n",
    "    super(Decoder, self).__init__()\n",
    "\n",
    "    self.seq_len, self.input_dim = seq_len, input_dim\n",
    "    self.hidden_dim, self.n_features = 2 * input_dim, n_features\n",
    "\n",
    "    self.rnn1 = nn.LSTM(\n",
    "      input_size=input_dim,\n",
    "      hidden_size=input_dim,\n",
    "      num_layers=1,\n",
    "      batch_first=True\n",
    "    )\n",
    "\n",
    "    self.rnn2 = nn.LSTM(\n",
    "      input_size=input_dim,\n",
    "      hidden_size=self.hidden_dim,\n",
    "      num_layers=1,\n",
    "      batch_first=True\n",
    "    )\n",
    "\n",
    "    self.output_layer = nn.Linear(self.hidden_dim, n_features)\n",
    "\n",
    "  def forward(self, x):\n",
    "    x = x.repeat(self.seq_len, self.n_features)\n",
    "    x = x.reshape((self.n_features, self.seq_len, self.input_dim))\n",
    "\n",
    "    x, (hidden_n, cell_n) = self.rnn1(x)\n",
    "    x, (hidden_n, cell_n) = self.rnn2(x)\n",
    "    x = x.reshape((self.seq_len, self.hidden_dim))\n",
    "\n",
    "    return self.output_layer(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "class RecurrentAutoencoder(nn.Module):\n",
    "\n",
    "  def __init__(self, seq_len, n_features, embedding_dim=64):\n",
    "    super(RecurrentAutoencoder, self).__init__()\n",
    "\n",
    "    self.encoder = Encoder(seq_len, n_features, embedding_dim).to(device)\n",
    "    self.decoder = Decoder(seq_len, embedding_dim, n_features).to(device)\n",
    "\n",
    "  def forward(self, x):\n",
    "    x = self.encoder(x)\n",
    "    x = self.decoder(x)\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = RecurrentAutoencoder(seq_len, n_features, 128)\n",
    "model = model.to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(model, train_dataset, val_dataset, n_epochs):\n",
    "  optimizer = torch.optim.Adam(model.parameters(), lr=1e-5)\n",
    "  criterion = nn.L1Loss(reduction='sum').to(device)\n",
    "  history = dict(train=[], val=[])\n",
    "\n",
    "  best_model_wts = copy.deepcopy(model.state_dict())\n",
    "  best_loss = 10000.0\n",
    "  \n",
    "  for epoch in range(1, n_epochs + 1):\n",
    "    model = model.train()\n",
    "\n",
    "    train_losses = []\n",
    "    for seq_true in train_dataset:\n",
    "      optimizer.zero_grad()\n",
    "\n",
    "      seq_true = seq_true.to(device)\n",
    "      seq_pred = model(seq_true)\n",
    "\n",
    "      loss = criterion(seq_pred, seq_true)\n",
    "\n",
    "      loss.backward()\n",
    "      optimizer.step()\n",
    "\n",
    "      train_losses.append(loss.item())\n",
    "\n",
    "    val_losses = []\n",
    "    model = model.eval()\n",
    "    with torch.no_grad():\n",
    "      for seq_true in val_dataset:\n",
    "\n",
    "        seq_true = seq_true.to(device)\n",
    "        seq_pred = model(seq_true)\n",
    "\n",
    "        loss = criterion(seq_pred, seq_true)\n",
    "        val_losses.append(loss.item())\n",
    "\n",
    "    train_loss = np.mean(train_losses)\n",
    "    val_loss = np.mean(val_losses)\n",
    "\n",
    "    history['train'].append(train_loss)\n",
    "    history['val'].append(val_loss)\n",
    "\n",
    "    if val_loss < best_loss:\n",
    "      best_loss = val_loss\n",
    "      best_model_wts = copy.deepcopy(model.state_dict())\n",
    "\n",
    "    print(f'Epoch {epoch}: train loss {train_loss} val loss {val_loss}')\n",
    "\n",
    "  model.load_state_dict(best_model_wts)\n",
    "  return model.eval(), history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1: train loss 7.188362143299209 val loss 10.97716400551796\n",
      "Epoch 2: train loss 7.061609009534848 val loss 10.861721575021743\n",
      "Epoch 3: train loss 7.050498445725213 val loss 10.759196199297905\n",
      "Epoch 4: train loss 6.953788480228296 val loss 10.529170639157295\n",
      "Epoch 5: train loss 6.819667443375704 val loss 10.165040640830993\n",
      "Epoch 6: train loss 6.738479136036151 val loss 9.992541402816773\n",
      "Epoch 7: train loss 6.614317140013353 val loss 9.791294813871383\n",
      "Epoch 8: train loss 6.452131601406606 val loss 9.845977551937104\n",
      "Epoch 9: train loss 6.408218298855107 val loss 9.518717983603478\n",
      "Epoch 10: train loss 6.238182526727458 val loss 9.489320818185806\n",
      "Epoch 11: train loss 6.234057898515033 val loss 9.289050007283688\n",
      "Epoch 12: train loss 6.202264721316195 val loss 9.33216683781147\n",
      "Epoch 13: train loss 6.1222286748257755 val loss 9.25562231683731\n",
      "Epoch 14: train loss 6.28904679132722 val loss 10.131250357687474\n",
      "Epoch 15: train loss 6.073458283509507 val loss 11.110596370041371\n",
      "Epoch 16: train loss 6.161009409070512 val loss 9.410741286337375\n",
      "Epoch 17: train loss 6.021980724913723 val loss 9.270162648320198\n",
      "Epoch 18: train loss 6.039738698001061 val loss 11.261745161175728\n",
      "Epoch 19: train loss 6.115699894411213 val loss 9.005485089600086\n",
      "Epoch 20: train loss 5.929479124247496 val loss 9.085405830562115\n",
      "Epoch 21: train loss 5.9451709336347704 val loss 10.106599816918374\n",
      "Epoch 22: train loss 5.924677521709248 val loss 9.108658552646636\n",
      "Epoch 23: train loss 5.908539707831785 val loss 8.512088710725308\n",
      "Epoch 24: train loss 5.818426139967793 val loss 9.660966635107995\n",
      "Epoch 25: train loss 5.953280821415979 val loss 8.881353223979474\n",
      "Epoch 26: train loss 5.762612802924465 val loss 8.796815033614635\n",
      "Epoch 27: train loss 6.0233771145752035 val loss 9.594733444035054\n",
      "Epoch 28: train loss 5.884460732465063 val loss 8.864278678059579\n",
      "Epoch 29: train loss 5.812006496509816 val loss 9.826941768109798\n",
      "Epoch 30: train loss 5.818772723404619 val loss 9.912757200181485\n",
      "Epoch 31: train loss 5.6362157306673275 val loss 9.99971356779337\n",
      "Epoch 32: train loss 5.646149147277506 val loss 10.41991896033287\n",
      "Epoch 33: train loss 5.785268129284996 val loss 9.180850823998451\n",
      "Epoch 34: train loss 5.6033936036066185 val loss 8.335838904500008\n",
      "Epoch 35: train loss 5.639638055348907 val loss 8.63649248892069\n",
      "Epoch 36: train loss 5.496804546314098 val loss 8.519725093007088\n",
      "Epoch 37: train loss 5.45866653725212 val loss 8.703529002845288\n",
      "Epoch 38: train loss 5.5987114275740755 val loss 7.895355988562107\n",
      "Epoch 39: train loss 5.397800835968839 val loss 8.712894497811794\n",
      "Epoch 40: train loss 5.615536484239427 val loss 9.162039269328117\n",
      "Epoch 41: train loss 5.468912153788235 val loss 8.45761411190033\n",
      "Epoch 42: train loss 5.371571399009675 val loss 7.897047997713089\n",
      "Epoch 43: train loss 5.266479222919346 val loss 7.756879223942756\n",
      "Epoch 44: train loss 5.268004156937551 val loss 7.302376929938793\n",
      "Epoch 45: train loss 5.267588082587811 val loss 9.525786200582981\n",
      "Epoch 46: train loss 5.373221474869966 val loss 7.77464767986536\n",
      "Epoch 47: train loss 5.32135811762273 val loss 8.563139500975609\n",
      "Epoch 48: train loss 5.412440207110692 val loss 7.822439022362232\n",
      "Epoch 49: train loss 5.406725116318149 val loss 7.7465962873101235\n",
      "Epoch 50: train loss 5.159115838695666 val loss 7.441654227137565\n",
      "Epoch 51: train loss 5.258179729568099 val loss 7.8189039084911345\n",
      "Epoch 52: train loss 5.251028409168931 val loss 7.743886261582374\n",
      "Epoch 53: train loss 5.242734800604854 val loss 7.365176262140274\n",
      "Epoch 54: train loss 5.16607915904282 val loss 8.90385246372223\n",
      "Epoch 55: train loss 5.2259016893338135 val loss 8.867705811679363\n",
      "Epoch 56: train loss 5.194666571919587 val loss 7.8739802262187\n",
      "Epoch 57: train loss 5.249035688970258 val loss 7.216172798454761\n",
      "Epoch 58: train loss 5.070800466896737 val loss 6.85274638992548\n",
      "Epoch 59: train loss 5.001781966748019 val loss 7.9121731922626495\n",
      "Epoch 60: train loss 5.1056873312374185 val loss 7.273398827493191\n",
      "Epoch 61: train loss 4.9760447693428445 val loss 7.714040239751339\n",
      "Epoch 62: train loss 5.0596290826797485 val loss 7.861135354578495\n",
      "Epoch 63: train loss 5.159216250112683 val loss 8.330905803859235\n",
      "Epoch 64: train loss 5.089301896732502 val loss 8.29973559910059\n",
      "Epoch 65: train loss 5.077050341283232 val loss 7.0275069693923\n",
      "Epoch 66: train loss 5.104413759262878 val loss 7.527150939106941\n",
      "Epoch 67: train loss 5.011002344894154 val loss 7.2104483036994935\n",
      "Epoch 68: train loss 5.031821377153834 val loss 7.462278498470783\n",
      "Epoch 69: train loss 4.93549394425606 val loss 7.092176900625229\n",
      "Epoch 70: train loss 4.902165090153798 val loss 7.506531097173691\n",
      "Epoch 71: train loss 5.044680127577813 val loss 8.029633438766002\n",
      "Epoch 72: train loss 5.078257953991061 val loss 7.620706166148186\n",
      "Epoch 73: train loss 5.02238962424341 val loss 7.4123868288993835\n",
      "Epoch 74: train loss 4.943703004306026 val loss 7.7856830370426175\n",
      "Epoch 75: train loss 5.058452894515122 val loss 7.487752008140087\n",
      "Epoch 76: train loss 4.983188311229369 val loss 6.923008233428002\n",
      "Epoch 77: train loss 4.91655006462474 val loss 7.989273576080799\n",
      "Epoch 78: train loss 5.02005114208878 val loss 7.763792626202107\n",
      "Epoch 79: train loss 5.119667072916329 val loss 7.779459906995297\n",
      "Epoch 80: train loss 4.815837256517206 val loss 6.943657228469848\n",
      "Epoch 81: train loss 4.94981735579097 val loss 8.068133530378342\n",
      "Epoch 82: train loss 4.966924919932637 val loss 7.317310967862606\n",
      "Epoch 83: train loss 4.9220759348634315 val loss 7.525299256563186\n",
      "Epoch 84: train loss 4.896259387266146 val loss 7.622523278236389\n",
      "Epoch 85: train loss 4.785323318562172 val loss 7.6916364338994025\n",
      "Epoch 86: train loss 4.898927788957376 val loss 7.980075287938118\n",
      "Epoch 87: train loss 4.846099264995195 val loss 7.757217339038849\n",
      "Epoch 88: train loss 4.897296259353884 val loss 7.217446682631969\n",
      "Epoch 89: train loss 4.8730314556085474 val loss 7.53286054366827\n",
      "Epoch 90: train loss 4.7620205258403665 val loss 8.212237065553666\n",
      "Epoch 91: train loss 4.7782756012951495 val loss 7.190876980245114\n",
      "Epoch 92: train loss 4.864514986617783 val loss 7.239515961110592\n",
      "Epoch 93: train loss 4.675690712807716 val loss 8.097305634975433\n",
      "Epoch 94: train loss 4.7876884221515885 val loss 6.623932673871517\n",
      "Epoch 95: train loss 4.722695211113884 val loss 7.503434510350227\n",
      "Epoch 96: train loss 4.745365573961938 val loss 7.175907820641995\n",
      "Epoch 97: train loss 4.659669167756538 val loss 7.756110722780227\n",
      "Epoch 98: train loss 4.662637394221146 val loss 6.194095297515393\n",
      "Epoch 99: train loss 4.6344790222781125 val loss 7.093573769330979\n",
      "Epoch 100: train loss 4.656862292254589 val loss 7.243564149618149\n"
     ]
    }
   ],
   "source": [
    "model, history = train_model(\n",
    "  model, \n",
    "  train_dataset, \n",
    "  val_dataset, \n",
    "  n_epochs=100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 126724 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 470,
       "width": 577
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.figure().gca()\n",
    "\n",
    "ax.plot(history['train'])\n",
    "ax.plot(history['val'])\n",
    "plt.ylabel('Loss')\n",
    "plt.xlabel('Epoch')\n",
    "plt.legend(['train', 'test'])\n",
    "plt.title('Loss over training epochs')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "date = datetime.date.today().strftime('%y-%m-%d')\n",
    "MODEL_PATH = f'model_save/model_{date}.pth'\n",
    "\n",
    "torch.save(model, MODEL_PATH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reload the model\n",
    "# model.torch.load('model.pth')\n",
    "# model = model.to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check reconstruction error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(model, dataset):\n",
    "  predictions, losses = [], []\n",
    "  criterion = nn.L1Loss(reduction='sum').to(device)\n",
    "  with torch.no_grad():\n",
    "    model = model.eval()\n",
    "    for seq_true in dataset:\n",
    "      seq_true = seq_true.to(device)\n",
    "      seq_pred = model(seq_true)\n",
    "\n",
    "      loss = criterion(seq_pred, seq_true)\n",
    "\n",
    "      predictions.append(seq_pred.cpu().numpy().flatten())\n",
    "      losses.append(loss.item())\n",
    "  return predictions, losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[6.475390911102295,\n",
       " 4.057796478271484,\n",
       " 3.6951184272766113,\n",
       " 3.5229454040527344,\n",
       " 3.1105241775512695,\n",
       " 2.7158966064453125,\n",
       " 4.410896301269531,\n",
       " 3.1737751960754395,\n",
       " 5.078624725341797,\n",
       " 564.0125732421875,\n",
       " 118.96932983398438,\n",
       " 11.067346572875977,\n",
       " 7.681413650512695,\n",
       " 2.2629685401916504,\n",
       " 2.2459702491760254,\n",
       " 5.009679794311523,\n",
       " 3.5274581909179688,\n",
       " 4.4660749435424805,\n",
       " 4.193930625915527,\n",
       " 2.288841724395752,\n",
       " 5.305790424346924,\n",
       " 506.50933837890625,\n",
       " 170.89035034179688,\n",
       " 2.9803903102874756,\n",
       " 1.8732764720916748,\n",
       " 0.4762537479400635,\n",
       " 17.63483428955078,\n",
       " 13.621052742004395,\n",
       " 1.6280626058578491,\n",
       " 1.6280626058578491,\n",
       " 1.6280626058578491,\n",
       " 9.569190979003906,\n",
       " 3.950976610183716,\n",
       " 10.132783889770508,\n",
       " 4.41414213180542,\n",
       " 3.841397762298584,\n",
       " 14.235858917236328,\n",
       " 2.4932689666748047,\n",
       " 0.3788847029209137,\n",
       " 10.606548309326172,\n",
       " 2.6478137969970703,\n",
       " 3.971092700958252,\n",
       " 1.9757204055786133,\n",
       " 2.9639644622802734,\n",
       " 5.554292678833008,\n",
       " 2.3994808197021484,\n",
       " 9.777814865112305,\n",
       " 4.184451103210449,\n",
       " 3.016322135925293,\n",
       " 7.93351411819458,\n",
       " 3.2564563751220703,\n",
       " 7.106253623962402,\n",
       " 7.106253623962402,\n",
       " 7.106253623962402,\n",
       " 25.282350540161133,\n",
       " 5.763306617736816,\n",
       " 1.9292466640472412,\n",
       " 3.0547077655792236,\n",
       " 5.411063194274902,\n",
       " 2.4595675468444824,\n",
       " 2.9209418296813965,\n",
       " 5.213022232055664,\n",
       " 2.082340717315674,\n",
       " 5.661027908325195,\n",
       " 3.1858878135681152,\n",
       " 9.611069679260254,\n",
       " 3.48431396484375,\n",
       " 3.612818717956543,\n",
       " 3.256265163421631,\n",
       " 4.243818283081055,\n",
       " 3.121319055557251,\n",
       " 7.132509231567383,\n",
       " 2.606898307800293,\n",
       " 2.889389991760254,\n",
       " 5.744939804077148,\n",
       " 2.9277758598327637,\n",
       " 2.3206820487976074,\n",
       " 9.055814743041992,\n",
       " 2.2245290279388428,\n",
       " 5.559058666229248,\n",
       " 3.029712200164795,\n",
       " 2.92405366897583,\n",
       " 4.436098575592041,\n",
       " 2.8357269763946533,\n",
       " 3.2960450649261475,\n",
       " 3.0109457969665527,\n",
       " 8.913893699645996,\n",
       " 2.2014102935791016,\n",
       " 5.028566360473633,\n",
       " 4.177780628204346,\n",
       " 3.0008766651153564,\n",
       " 5.832966327667236,\n",
       " 3.4131276607513428,\n",
       " 3.6118223667144775,\n",
       " 3.449504852294922,\n",
       " 5.172828197479248,\n",
       " 4.8995442390441895,\n",
       " 3.3583712577819824,\n",
       " 3.749042510986328,\n",
       " 11.917932510375977]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, losses = predict(model, test_normal_dataset)\n",
    "losses[:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: ylabel='Density'>"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 66208 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 611
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, losses = predict(model, train_dataset)\n",
    "sns.kdeplot(losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: ylabel='Density'>"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 71656 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 602
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# _, losses = predict(model, train_dataset)\n",
    "_, losses = predict(model, test_normal_dataset)\n",
    "sns.kdeplot(losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: ylabel='Density'>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 76276 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 630
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# _, losses = predict(model, train_dataset)\n",
    "_, losses = predict(model, test_anomaly_dataset)\n",
    "sns.kdeplot(losses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute THRESHOLD with training set data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions, losses = predict(model, train_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_array = np.array(losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "stdev = np.std(loss_array)\n",
    "mean = np.mean(loss_array)\n",
    "THRESHOLD = mean + stdev "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check on test_normal_dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of intervals exceeding 1 std dev loss: 33/416\n"
     ]
    }
   ],
   "source": [
    "_, losses = predict(model, test_normal_dataset)\n",
    "exceed_count = sum(l > THRESHOLD for l in losses)\n",
    "print(f'number of intervals exceeding 1 std dev loss: {exceed_count}/{len(test_normal_dataset)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check on test_anomaly_dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of intervals exceeding 1 std dev loss: 333/416\n"
     ]
    }
   ],
   "source": [
    "_, losses = predict(model, test_anomaly_dataset)\n",
    "exceed_count = sum(l > THRESHOLD for l in losses)\n",
    "print(f'number of intervals exceeding 1 std dev loss: {exceed_count}/{len(test_normal_dataset)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot construction error vs original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_prediction(data, model, title, ax):\n",
    "    predictions, pred_losses = predict(model, [data])\n",
    "\n",
    "    ax.plot(data, label='true')\n",
    "    ax.plot(predictions[0], label='reconstructed')\n",
    "    ax.set_title(f'{title} (loss: {np.around(pred_losses[0], 2)})')\n",
    "    ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 215592 characters]",
      "text/plain": [
       "<Figure size 2200x800 with 12 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 777,
       "width": 2174
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(\n",
    "    nrows=2,\n",
    "    ncols=6,\n",
    "    sharey=True,\n",
    "    sharex=True,\n",
    "    figsize=(22, 8)\n",
    ")\n",
    "\n",
    "sample_size = 6\n",
    "sample_indices = random.sample(range(0,len(test_normal_dataset)), sample_size)\n",
    "\n",
    "sampled_test_normal_dataset = [test_normal_dataset[i] for i in sample_indices]\n",
    "sampled_test_anomaly_dataset = [test_anomaly_dataset[i] for i in sample_indices]\n",
    "\n",
    "for i, data in enumerate(sampled_test_normal_dataset):\n",
    "    plot_prediction(data, model, title='Normal', ax=axs[0, i])\n",
    "\n",
    "for i, data in enumerate(sampled_test_anomaly_dataset):\n",
    "    plot_prediction(data, model, title='Anomaly', ax=axs[1, i])\n",
    "\n",
    "fig.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assess quality of prediction of flagged intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we want to identify the count of anomalies in each interval in the \"test_anomaly_dataset\"\n",
    "# range: 45000:60000\n",
    "start_point = 0\n",
    "end_point = 60000\n",
    "fault_segments = [ fault_data[start_point+seq_len * i:start_point + seq_len * i + seq_len] for i in range(len(test_anomaly_dataset)) ]\n",
    "# len(test_anomaly_dataset)\n",
    "anomaly_count_list = [ fault_segments[i].count(1) for i in range(len(fault_segments))]\n",
    "anomaly_flag_actual = [ count > 0 for count in anomaly_count_list ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# there are 50 - 100 minute segments where there is an anomaly\n",
    "sum(np.array(anomaly_count_list) > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, losses = predict(model, test_anomaly_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "anomaly_flag_prediction = [ l > THRESHOLD for l in losses ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[642, 289],\n",
       "       [ 25,  44]])"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "confusion_matrix(anomaly_flag_actual, anomaly_flag_prediction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "68.60000000000001"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy_score(anomaly_flag_actual, anomaly_flag_prediction) * 100"
   ]
  },
  {
   "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.11.5"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
