ECG-RNN-Time-series-Anomaly-Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fcac29da",
   "metadata": {},
   "source": [
    "# DOWNLOADING DATASET:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5358ba1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import os\n",
    "from pathlib import Path\n",
    "import pickle\n",
    "from shutil import unpack_archive\n",
    "\n",
    "urls = dict()\n",
    "urls['ecg']=['http://www.cs.ucr.edu/~eamonn/discords/ECG_data.zip',\n",
    "             'http://www.cs.ucr.edu/~eamonn/discords/mitdbx_mitdbx_108.txt',\n",
    "             'http://www.cs.ucr.edu/~eamonn/discords/qtdbsele0606.txt',\n",
    "             'http://www.cs.ucr.edu/~eamonn/discords/chfdbchf15.txt',\n",
    "             'http://www.cs.ucr.edu/~eamonn/discords/qtdbsel102.txt']\n",
    "urls['gesture']=['http://www.cs.ucr.edu/~eamonn/discords/ann_gun_CentroidA']\n",
    "urls['space_shuttle']=['http://www.cs.ucr.edu/~eamonn/discords/TEK16.txt',\n",
    "                       'http://www.cs.ucr.edu/~eamonn/discords/TEK17.txt',\n",
    "                       'http://www.cs.ucr.edu/~eamonn/discords/TEK14.txt']\n",
    "urls['respiration']=['http://www.cs.ucr.edu/~eamonn/discords/nprs44.txt',\n",
    "                     'http://www.cs.ucr.edu/~eamonn/discords/nprs43.txt']\n",
    "urls['power_demand']=['http://www.cs.ucr.edu/~eamonn/discords/power_data.txt']\n",
    "\n",
    "for dataname in urls:\n",
    "    raw_dir = Path('dataset', dataname, 'raw')\n",
    "    raw_dir.mkdir(parents=True, exist_ok=True)\n",
    "    for url in urls[dataname]:\n",
    "        filename = raw_dir.joinpath(Path(url).name)\n",
    "        print('Downloading', url)\n",
    "        resp =requests.get(url)\n",
    "        filename.write_bytes(resp.content)\n",
    "        if filename.suffix=='':\n",
    "            filename.rename(filename.with_suffix('.txt'))\n",
    "        print('Saving to', filename.with_suffix('.txt'))\n",
    "        if filename.suffix=='.zip':\n",
    "            print('Extracting to', filename)\n",
    "            unpack_archive(str(filename), extract_dir=str(raw_dir))\n",
    "\n",
    "\n",
    "\n",
    "    for filepath in raw_dir.glob('*.txt'):\n",
    "        with open(str(filepath)) as f:\n",
    "            # Label anomaly points as 1 in the dataset\n",
    "            labeled_data=[]\n",
    "            for i, line in enumerate(f):\n",
    "                tokens = [float(token) for token in line.split()]\n",
    "                if raw_dir.parent.name== 'ecg':\n",
    "                    # Remove time-step channel\n",
    "                    tokens.pop(0)\n",
    "                if filepath.name == 'chfdbchf15.txt':\n",
    "                    tokens.append(1.0) if 2250 < i < 2400 else tokens.append(0.0)\n",
    "                elif filepath.name == 'xmitdb_x108_0.txt':\n",
    "                    tokens.append(1.0) if 4020 < i < 4400 else tokens.append(0.0)\n",
    "                elif filepath.name == 'mitdb__100_180.txt':\n",
    "                    tokens.append(1.0) if 1800 < i < 1990 else tokens.append(0.0)\n",
    "                elif filepath.name == 'chfdb_chf01_275.txt':\n",
    "                    tokens.append(1.0) if 2330 < i < 2600 else tokens.append(0.0)\n",
    "                elif filepath.name == 'ltstdb_20221_43.txt':\n",
    "                    tokens.append(1.0) if 650 < i < 780 else tokens.append(0.0)\n",
    "                elif filepath.name == 'ltstdb_20321_240.txt':\n",
    "                    tokens.append(1.0) if 710 < i < 850 else tokens.append(0.0)\n",
    "                elif filepath.name == 'chfdb_chf13_45590.txt':\n",
    "                    tokens.append(1.0) if 2800 < i < 2960 else tokens.append(0.0)\n",
    "                elif filepath.name == 'stdb_308_0.txt':\n",
    "                    tokens.append(1.0) if 2290 < i < 2550 else tokens.append(0.0)\n",
    "                elif filepath.name == 'qtdbsel102.txt':\n",
    "                    tokens.append(1.0) if 4230 < i < 4430 else tokens.append(0.0)\n",
    "                elif filepath.name == 'ann_gun_CentroidA.txt':\n",
    "                    tokens.append(1.0) if 2070 < i < 2810 else tokens.append(0.0)\n",
    "                elif filepath.name == 'TEK16.txt':\n",
    "                    tokens.append(1.0) if 4270 < i < 4370 else tokens.append(0.0)\n",
    "                elif filepath.name == 'TEK17.txt':\n",
    "                    tokens.append(1.0) if 2100 < i < 2145 else tokens.append(0.0)\n",
    "                elif filepath.name == 'TEK14.txt':\n",
    "                    tokens.append(1.0) if 1100 < i < 1200 or 1455 < i < 1955 else tokens.append(0.0)\n",
    "                elif filepath.name == 'nprs44.txt':\n",
    "                    tokens.append(1.0) if 16192 < i < 16638 or 20457 < i < 20911 else tokens.append(0.0)\n",
    "                elif filepath.name == 'nprs43.txt':\n",
    "                    tokens.append(1.0) if 12929 < i < 13432 or 14877 < i < 15086 or 15729 < i < 15924 else tokens.append(0.0)\n",
    "                elif filepath.name == 'power_data.txt':\n",
    "                    tokens.append(1.0) if 8254 < i < 8998 or 11348 < i < 12143 or 33883 < i < 34601 else tokens.append(0.0)\n",
    "                labeled_data.append(tokens)\n",
    "\n",
    "            # Fill in the point where there is no signal value.\n",
    "            if filepath.name == 'ann_gun_CentroidA.txt':\n",
    "                for i, datapoint in enumerate(labeled_data):\n",
    "                    for j,channel in enumerate(datapoint[:-1]):\n",
    "                        if channel == 0.0:\n",
    "                            labeled_data[i][j] = 0.5 * labeled_data[i - 1][j] + 0.5 * labeled_data[i + 1][j]\n",
    "\n",
    "            # Save the labeled dataset as .pkl extension\n",
    "            labeled_whole_dir = raw_dir.parent.joinpath('labeled', 'whole')\n",
    "            labeled_whole_dir.mkdir(parents=True, exist_ok=True)\n",
    "            with open(str(labeled_whole_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                pickle.dump(labeled_data, pkl)\n",
    "\n",
    "            # Divide the labeled dataset into trainset and testset, then save them\n",
    "            labeled_train_dir = raw_dir.parent.joinpath('labeled','train')\n",
    "            labeled_train_dir.mkdir(parents=True,exist_ok=True)\n",
    "            labeled_test_dir = raw_dir.parent.joinpath('labeled','test')\n",
    "            labeled_test_dir.mkdir(parents=True,exist_ok=True)\n",
    "            if filepath.name == 'chfdb_chf13_45590.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[:2439], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[2439:3726], pkl)\n",
    "            elif filepath.name == 'chfdb_chf01_275.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[:1833], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[1833:3674], pkl)\n",
    "            elif filepath.name == 'chfdbchf15.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[3381:14244], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[33:3381], pkl)\n",
    "            elif filepath.name == 'qtdbsel102.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[10093:44828], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[211:10093], pkl)\n",
    "            elif filepath.name == 'mitdb__100_180.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[2328:5271], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[73:2328], pkl)\n",
    "            elif filepath.name == 'stdb_308_0.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[2986:5359], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[265:2986], pkl)\n",
    "            elif filepath.name == 'ltstdb_20321_240.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[1520:3531], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[73:1520], pkl)\n",
    "            elif filepath.name == 'xmitdb_x108_0.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[424:3576], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[3576:5332], pkl)\n",
    "            elif filepath.name == 'ltstdb_20221_43.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[1121:3731], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[0:1121], pkl)\n",
    "            elif filepath.name == 'ann_gun_CentroidA.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[3000:], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[:3000], pkl)\n",
    "            elif filepath.name == 'nprs44.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[363:12955], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[12955:24082], pkl)\n",
    "            elif filepath.name == 'nprs43.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[4285:10498], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[10498:17909], pkl)\n",
    "            elif filepath.name == 'power_data.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[15287:33432], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[501:15287], pkl)\n",
    "            elif filepath.name == 'TEK17.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[2469:4588], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[1543:2469], pkl)\n",
    "            elif filepath.name == 'TEK16.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[521:3588], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[3588:4539], pkl)\n",
    "            elif filepath.name == 'TEK14.txt':\n",
    "                with open(str(labeled_train_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[2089:4098], pkl)\n",
    "                with open(str(labeled_test_dir.joinpath(filepath.name).with_suffix('.pkl')), 'wb') as pkl:\n",
    "                    pickle.dump(labeled_data[97:2089], pkl)\n",
    "\n",
    "nyc_taxi_raw_path = Path('dataset/nyc_taxi/raw/nyc_taxi.csv')\n",
    "labeled_data = []\n",
    "with open(str(nyc_taxi_raw_path),'r') as f:\n",
    "    for i, line in enumerate(f):\n",
    "        tokens = [float(token) for token in line.strip().split(',')[1:]]\n",
    "        tokens.append(1) if 150 < i < 250 or   \\\n",
    "                            5970 < i < 6050 or \\\n",
    "                            8500 < i < 8650 or \\\n",
    "                            8750 < i < 8890 or \\\n",
    "                            10000 < i < 10200 or \\\n",
    "                            14700 < i < 14800 \\\n",
    "                          else tokens.append(0)\n",
    "        labeled_data.append(tokens)\n",
    "nyc_taxi_train_path = nyc_taxi_raw_path.parent.parent.joinpath('labeled','train',nyc_taxi_raw_path.name).with_suffix('.pkl')\n",
    "nyc_taxi_train_path.parent.mkdir(parents=True, exist_ok=True)\n",
    "with open(str(nyc_taxi_train_path),'wb') as pkl:\n",
    "    pickle.dump(labeled_data[:13104], pkl)\n",
    "\n",
    "nyc_taxi_test_path = nyc_taxi_raw_path.parent.parent.joinpath('labeled','test',nyc_taxi_raw_path.name).with_suffix('.pkl')\n",
    "nyc_taxi_test_path.parent.mkdir(parents=True, exist_ok=True)\n",
    "with open(str(nyc_taxi_test_path),'wb') as pkl:\n",
    "    pickle.dump(labeled_data[13104:], pkl)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaa67d22",
   "metadata": {},
   "source": [
    "# PREPROCESSING:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a8eb02f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "   \n",
    "import os\n",
    "import torch\n",
    "from torch import device\n",
    "import glob\n",
    "import datetime\n",
    "import numpy as np\n",
    "import shutil\n",
    "from pathlib import Path\n",
    "import pickle\n",
    "\n",
    "def normalization(seqData,max,min):\n",
    "    return (seqData -min)/(max-min)\n",
    "\n",
    "def standardization(seqData,mean,std):\n",
    "    return (seqData-mean)/std\n",
    "\n",
    "def reconstruct(seqData,mean,std):\n",
    "    return seqData*std+mean\n",
    "\n",
    "class PickleDataLoad(object):\n",
    "    def __init__(self, data_type, filename, augment_test_data=True):\n",
    "        self.augment_test_data=augment_test_data\n",
    "        self.trainData, self.trainLabel = self.preprocessing(Path('dataset',data_type,'labeled','train',filename),train=True)\n",
    "        self.testData, self.testLabel = self.preprocessing(Path('dataset',data_type,'labeled','test',filename),train=False)\n",
    "\n",
    "    def augmentation(self,data,label,noise_ratio=0.05,noise_interval=0.0005,max_length=100000):\n",
    "        noiseSeq = torch.randn(data.size())\n",
    "        augmentedData = data.clone()\n",
    "        augmentedLabel = label.clone()\n",
    "        for i in np.arange(0, noise_ratio, noise_interval):\n",
    "            scaled_noiseSeq = noise_ratio * self.std.expand_as(data) * noiseSeq\n",
    "            augmentedData = torch.cat([augmentedData, data + scaled_noiseSeq], dim=0)\n",
    "            augmentedLabel = torch.cat([augmentedLabel, label])\n",
    "            if len(augmentedData) > max_length:\n",
    "                augmentedData = augmentedData[:max_length]\n",
    "                augmentedLabel = augmentedLabel[:max_length]\n",
    "                break\n",
    "\n",
    "        return augmentedData, augmentedLabel\n",
    "\n",
    "    def preprocessing(self, path, train=True):\n",
    "        \"\"\" Read, Standardize, Augment \"\"\"\n",
    "\n",
    "        with open(str(path), 'rb') as f:\n",
    "            data = torch.FloatTensor(pickle.load(f))\n",
    "            label = data[:,-1]\n",
    "            data = data[:,:-1]\n",
    "        if train:\n",
    "            self.mean = data.mean(dim=0)\n",
    "            self.std= data.std(dim=0)\n",
    "            self.length = len(data)\n",
    "            data,label = self.augmentation(data,label)\n",
    "        else:\n",
    "            if self.augment_test_data:\n",
    "                data, label = self.augmentation(data, label)\n",
    "\n",
    "        data = standardization(data,self.mean,self.std)\n",
    "\n",
    "        return data,label\n",
    "\n",
    "    def batchify(self,args,data, bsz):\n",
    "        nbatch = data.size(0) // bsz\n",
    "        trimmed_data = data.narrow(0,0,nbatch * bsz)\n",
    "        batched_data = trimmed_data.contiguous().view(bsz, -1, trimmed_data.size(-1)).transpose(0,1)\n",
    "        batched_data = batched_data.to(device(args.device))\n",
    "        return \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b7e207",
   "metadata": {},
   "source": [
    "# Anomaly Detection :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8046a46a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.autograd import Variable\n",
    "import torch\n",
    "import numpy as np\n",
    "\n",
    "def fit_norm_distribution_param(args, model, train_dataset, channel_idx=0):\n",
    "    predictions = []\n",
    "    organized = []\n",
    "    errors = []\n",
    "    with torch.no_grad():\n",
    "        # Turn on evaluation mode which disables dropout.\n",
    "        model.eval()\n",
    "        pasthidden = model.init_hidden(1)\n",
    "        for t in range(len(train_dataset)):\n",
    "            out, hidden = model.forward(train_dataset[t].unsqueeze(0), pasthidden)\n",
    "            predictions.append([])\n",
    "            organized.append([])\n",
    "            errors.append([])\n",
    "            predictions[t].append(out.data.cpu()[0][0][channel_idx])\n",
    "            pasthidden = model.repackage_hidden(hidden)\n",
    "            for prediction_step in range(1,args.prediction_window_size):\n",
    "                out, hidden = model.forward(out, hidden)\n",
    "                predictions[t].append(out.data.cpu()[0][0][channel_idx])\n",
    "\n",
    "            if t >= args.prediction_window_size:\n",
    "                for step in range(args.prediction_window_size):\n",
    "                    organized[t].append(predictions[step+t-args.prediction_window_size][args.prediction_window_size-1-step])\n",
    "                organized[t]= torch.FloatTensor(organized[t]).to(args.device)\n",
    "                errors[t] = organized[t] - train_dataset[t][0][channel_idx]\n",
    "                errors[t] = errors[t].unsqueeze(0)\n",
    "\n",
    "    errors_tensor = torch.cat(errors[args.prediction_window_size:],dim=0)\n",
    "    mean = errors_tensor.mean(dim=0)\n",
    "    cov = errors_tensor.t().mm(errors_tensor)/errors_tensor.size(0) - mean.unsqueeze(1).mm(mean.unsqueeze(0))\n",
    "    # cov: positive-semidefinite and symmetric.\n",
    "\n",
    "    return mean, cov\n",
    "\n",
    "\n",
    "def anomalyScore(args, model, dataset, mean, cov, channel_idx=0, score_predictor=None):\n",
    "    predictions = []\n",
    "    rearranged = []\n",
    "    errors = []\n",
    "    hiddens = []\n",
    "    predicted_scores = []\n",
    "    with torch.no_grad():\n",
    "        # Turn on evaluation mode which disables dropout.\n",
    "        model.eval()\n",
    "        pasthidden = model.init_hidden(1)\n",
    "        for t in range(len(dataset)):\n",
    "            out, hidden = model.forward(dataset[t].unsqueeze(0), pasthidden)\n",
    "            predictions.append([])\n",
    "            rearranged.append([])\n",
    "            errors.append([])\n",
    "            hiddens.append(model.extract_hidden(hidden))\n",
    "            if score_predictor is not None:\n",
    "                predicted_scores.append(score_predictor.predict(model.extract_hidden(hidden).numpy()))\n",
    "\n",
    "            predictions[t].append(out.data.cpu()[0][0][channel_idx])\n",
    "            pasthidden = model.repackage_hidden(hidden)\n",
    "            for prediction_step in range(1, args.prediction_window_size):\n",
    "                out, hidden = model.forward(out, hidden)\n",
    "                predictions[t].append(out.data.cpu()[0][0][channel_idx])\n",
    "\n",
    "            if t >= args.prediction_window_size:\n",
    "                for step in range(args.prediction_window_size):\n",
    "                    rearranged[t].append(\n",
    "                        predictions[step + t - args.prediction_window_size][args.prediction_window_size - 1 - step])\n",
    "                rearranged[t] =torch.FloatTensor(rearranged[t]).to(args.device).unsqueeze(0)\n",
    "                errors[t] = rearranged[t] - dataset[t][0][channel_idx]\n",
    "            else:\n",
    "                rearranged[t] = torch.zeros(1,args.prediction_window_size).to(args.device)\n",
    "                errors[t] = torch.zeros(1, args.prediction_window_size).to(args.device)\n",
    "\n",
    "    predicted_scores = np.array(predicted_scores)\n",
    "    scores = []\n",
    "    for error in errors:\n",
    "        mult1 = error-mean.unsqueeze(0) # [ 1 * prediction_window_size ]\n",
    "        mult2 = torch.inverse(cov) # [ prediction_window_size * prediction_window_size ]\n",
    "        mult3 = mult1.t() # [ prediction_window_size * 1 ]\n",
    "        score = torch.mm(mult1,torch.mm(mult2,mult3))\n",
    "        scores.append(score[0][0])\n",
    "\n",
    "    scores = torch.stack(scores)\n",
    "    rearranged = torch.cat(rearranged,dim=0)\n",
    "    errors = torch.cat(errors,dim=0)\n",
    "\n",
    "    return scores, rearranged, errors, hiddens, predicted_scores\n",
    "\n",
    "\n",
    "def get_precision_recall(args, score, label, num_samples, beta=1.0, sampling='log', predicted_score=None):\n",
    "    '''\n",
    "    :param args:\n",
    "    :param score: anomaly scores\n",
    "    :param label: anomaly labels\n",
    "    :param num_samples: the number of threshold samples\n",
    "    :param beta:\n",
    "    :param scale:\n",
    "    :return:\n",
    "    '''\n",
    "    if predicted_score is not None:\n",
    "        score = score - torch.FloatTensor(predicted_score).squeeze().to(args.device)\n",
    "\n",
    "    maximum = score.max()\n",
    "    if sampling=='log':\n",
    "        # Sample thresholds logarithmically\n",
    "        # The sampled thresholds are logarithmically spaced between: math:`10 ^ {start}` and: math:`10 ^ {end}`.\n",
    "        th = torch.logspace(0, torch.log10(torch.tensor(maximum)), num_samples).to(args.device)\n",
    "    else:\n",
    "        # Sample thresholds equally\n",
    "        # The sampled thresholds are equally spaced points between: attr:`start` and: attr:`end`\n",
    "        th = torch.linspace(0, maximum, num_samples).to(args.device)\n",
    "\n",
    "    precision = []\n",
    "    recall = []\n",
    "\n",
    "    for i in range(len(th)):\n",
    "        anomaly = (score > th[i]).float()\n",
    "        idx = anomaly * 2 + label\n",
    "        tn = (idx == 0.0).sum().item()  # tn\n",
    "        fn = (idx == 1.0).sum().item()  # fn\n",
    "        fp = (idx == 2.0).sum().item()  # fp\n",
    "        tp = (idx == 3.0).sum().item()  # tp\n",
    "\n",
    "        p = tp / (tp + fp + 1e-7)\n",
    "        r = tp / (tp + fn + 1e-7)\n",
    "\n",
    "        if p != 0 and r != 0:\n",
    "            precision.append(p)\n",
    "            recall.append(r)\n",
    "\n",
    "    precision = torch.FloatTensor(precision)\n",
    "    recall = torch.FloatTensor(recall)\n",
    "\n",
    "\n",
    "    f1 = (1 + beta ** 2) * (precision * recall).div(beta ** 2 * precision + recall + 1e-7)\n",
    "\n",
    "    return precision, recall, f1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61c88a5f",
   "metadata": {},
   "source": [
    "# MODEL DEFINITION:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29efd432",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.nn as nn\n",
    "import torch\n",
    "from torch.autograd import Variable\n",
    "import torch.nn.functional as F\n",
    "import shutil\n",
    "from pathlib import Path\n",
    "\n",
    "class RNNPredictor(nn.Module):\n",
    "    \"\"\"Container module with an encoder, a recurrent module, and a decoder.\"\"\"\n",
    "\n",
    "    def __init__(self, rnn_type, enc_inp_size, rnn_inp_size, rnn_hid_size, dec_out_size, nlayers, dropout=0.5,\n",
    "                 tie_weights=False,res_connection=False):\n",
    "        super(RNNPredictor, self).__init__()\n",
    "        self.enc_input_size = enc_inp_size\n",
    "\n",
    "        self.drop = nn.Dropout(dropout)\n",
    "        self.encoder = nn.Linear(enc_inp_size, rnn_inp_size)\n",
    "        if rnn_type in ['LSTM', 'GRU']:\n",
    "            self.rnn = getattr(nn, rnn_type)(rnn_inp_size, rnn_hid_size, nlayers, dropout=dropout)\n",
    "        elif rnn_type == 'SRU':\n",
    "            from cuda_functional import SRU, SRUCell\n",
    "            self.rnn = SRU(input_size=rnn_inp_size,hidden_size=rnn_hid_size,num_layers=nlayers,dropout=dropout,\n",
    "                           use_tanh=False,use_selu=True,layer_norm=True)\n",
    "        else:\n",
    "            try:\n",
    "                nonlinearity = {'RNN_TANH': 'tanh', 'RNN_RELU': 'relu'}[rnn_type]\n",
    "            except KeyError:\n",
    "                raise ValueError( \"\"\"An invalid option for `--model` was supplied,\n",
    "                                 options are ['LSTM', 'GRU', 'SRU', 'RNN_TANH' or 'RNN_RELU']\"\"\")\n",
    "            self.rnn = nn.RNN(rnn_inp_size, rnn_hid_size, nlayers, nonlinearity=nonlinearity, dropout=dropout)\n",
    "        self.decoder = nn.Linear(rnn_hid_size, dec_out_size)\n",
    "\n",
    "\n",
    "        if tie_weights:\n",
    "            if rnn_hid_size != rnn_inp_size:\n",
    "                raise ValueError('When using the tied flag, nhid must be equal to emsize')\n",
    "            self.decoder.weight = self.encoder.weight\n",
    "        self.res_connection=res_connection\n",
    "        self.init_weights()\n",
    "        self.rnn_type = rnn_type\n",
    "        self.rnn_hid_size = rnn_hid_size\n",
    "        self.nlayers = nlayers\n",
    "        #self.layerNorm1=nn.LayerNorm(normalized_shape=rnn_inp_size)\n",
    "        #self.layerNorm2=nn.LayerNorm(normalized_shape=rnn_hid_size)\n",
    "\n",
    "    def init_weights(self):\n",
    "        initrange = 0.1\n",
    "        self.encoder.weight.data.uniform_(-initrange, initrange)\n",
    "        self.decoder.bias.data.fill_(0)\n",
    "        self.decoder.weight.data.uniform_(-initrange, initrange)\n",
    "\n",
    "    def forward(self, input, hidden, return_hiddens=False, noise=False):\n",
    "        emb = self.drop(self.encoder(input.contiguous().view(-1,self.enc_input_size))) # [(seq_len x batch_size) * feature_size]\n",
    "        emb = emb.contiguous().view(-1, input.size(1), self.rnn_hid_size) # [ seq_len * batch_size * feature_size]\n",
    "        if noise:\n",
    "            # emb_noise = Variable(torch.randn(emb.size()))\n",
    "            # hidden_noise = Variable(torch.randn(hidden[0].size()))\n",
    "            # if next(self.parameters()).is_cuda:\n",
    "            #     emb_noise=emb_noise.cuda()\n",
    "            #     hidden_noise=hidden_noise.cuda()\n",
    "            # emb = emb+emb_noise\n",
    "            hidden = (F.dropout(hidden[0],training=True,p=0.9),F.dropout(hidden[1],training=True,p=0.9))\n",
    "\n",
    "        #emb = self.layerNorm1(emb)\n",
    "        output, hidden = self.rnn(emb, hidden)\n",
    "        #output = self.layerNorm2(output)\n",
    "\n",
    "        output = self.drop(output)\n",
    "        decoded = self.decoder(output.contiguous().view(output.size(0)*output.size(1), output.size(2))) # [(seq_len x batch_size) * feature_size]\n",
    "        decoded = decoded.contiguous().view(output.size(0), output.size(1), decoded.size(1)) # [ seq_len * batch_size * feature_size]\n",
    "        if self.res_connection:\n",
    "            decoded = decoded + input\n",
    "        if return_hiddens:\n",
    "            return decoded,hidden,output\n",
    "\n",
    "        return decoded, hidden\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    def init_hidden(self, bsz):\n",
    "        weight = next(self.parameters()).data\n",
    "        if self.rnn_type == 'LSTM':\n",
    "            return (Variable(weight.new(self.nlayers, bsz, self.rnn_hid_size).zero_()),\n",
    "                    Variable(weight.new(self.nlayers, bsz, self.rnn_hid_size).zero_()))\n",
    "        else:\n",
    "            return Variable(weight.new(self.nlayers, bsz, self.rnn_hid_size).zero_())\n",
    "\n",
    "    def repackage_hidden(self,h):\n",
    "        \"\"\"Wraps hidden states in new Variables, to detach them from their history.\"\"\"\n",
    "        if type(h) == tuple:\n",
    "            return tuple(self.repackage_hidden(v) for v in h)\n",
    "        else:\n",
    "            return h.detach()\n",
    "\n",
    "    def save_checkpoint(self,state, is_best):\n",
    "        print(\"=> saving checkpoint ..\")\n",
    "        args = state['args']\n",
    "        checkpoint_dir = r'/Users/immanuvelprathaps/Desktop/RNN-Time-series-Anomaly-Detection-master/save/ecg/checkpoint/chfdb_chf13_45590.pth' \n",
    "#         checkpoint_dir = Path('save',args.data,'checkpoint')\n",
    "#         checkpoint_dir.mkdir(parents=True,exist_ok=True)\n",
    "#         checkpoint = checkpoint_dir.joinpath(args.filename).with_suffix('.pth')\n",
    "\n",
    "#         torch.save(state, checkpoint)\n",
    "        if is_best:\n",
    "            model_best_dir = Path('save',args.data,'model_best')\n",
    "            model_best_dir.mkdir(parents=True,exist_ok=True)\n",
    "\n",
    "            shutil.copyfile(checkpoint, model_best_dir.joinpath(args.filename).with_suffix('.pth'))\n",
    "\n",
    "        print('=> checkpoint saved.')\n",
    "\n",
    "    def extract_hidden(self, hidden):\n",
    "        if self.rnn_type == 'LSTM':\n",
    "            return hidden[0][-1].data.cpu()  # hidden state last layer (hidden[1] is cell state)\n",
    "        else:\n",
    "            return hidden[-1].data.cpu()  # last layer\n",
    "\n",
    "    def initialize(self,args,feature_dim):\n",
    "        self.__init__(rnn_type = args.model,\n",
    "                           enc_inp_size=feature_dim,\n",
    "                           rnn_inp_size = args.emsize,\n",
    "                           rnn_hid_size = args.nhid,\n",
    "                           dec_out_size=feature_dim,\n",
    "                           nlayers = args.nlayers,\n",
    "                           dropout = args.dropout,\n",
    "                           tie_weights= args.tied,\n",
    "                           res_connection=args.res_connection)\n",
    "        self.to(args.device)\n",
    "\n",
    "    def load_checkpoint(self, args, checkpoint, feature_dim):\n",
    "        start_epoch = checkpoint['epoch'] +1\n",
    "        best_val_loss = checkpoint['best_loss']\n",
    "        args_ = checkpoint['args']\n",
    "        args_.resume = args.resume\n",
    "        args_.pretrained = args.pretrained\n",
    "        args_.epochs = args.epochs\n",
    "        args_.save_interval = args.save_interval\n",
    "        args_.prediction_window_size=args.prediction_window_size\n",
    "        self.initialize(args_, feature_dim=feature_dim)\n",
    "        self.load_state_dict(checkpoint['state_dict'])\n",
    "\n",
    "        return args_, start_epoch, best_val_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75babd56",
   "metadata": {},
   "source": [
    "# TRAINING MODEL:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0ea7d0e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import argparse\n",
    "import time\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import preprocess_data\n",
    "from model import model\n",
    "from torch import optim\n",
    "import os\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "from pathlib import Path\n",
    "from anomalyDetector import fit_norm_distribution_param\n",
    "# from torchvision.utils import save_image\n",
    "\n",
    "parser = argparse.ArgumentParser(description='PyTorch RNN Prediction Model on Time-series Dataset')\n",
    "parser.add_argument('--data', type=str, default='ecg',\n",
    "                    help='type of the dataset (ecg, gesture, power_demand, space_shuttle, respiration, nyc_taxi')\n",
    "parser.add_argument('--filename', type=str, default='chfdb_chf13_45590.pkl',\n",
    "                    help='filename of the dataset')\n",
    "parser.add_argument('--model', type=str, default='LSTM',\n",
    "                    help='type of recurrent net (RNN_TANH, RNN_RELU, LSTM, GRU, SRU)')\n",
    "parser.add_argument('--augment', type=bool, default=True,\n",
    "                    help='augment')\n",
    "parser.add_argument('--emsize', type=int, default=32,\n",
    "                    help='size of rnn input features')\n",
    "parser.add_argument('--nhid', type=int, default=32,\n",
    "                    help='number of hidden units per layer')\n",
    "parser.add_argument('--nlayers', type=int, default=2,\n",
    "                    help='number of layers')\n",
    "parser.add_argument('--res_connection', action='store_true',\n",
    "                    help='residual connection')\n",
    "parser.add_argument('--lr', type=float, default=0.0002,\n",
    "                    help='initial learning rate')\n",
    "parser.add_argument('--weight_decay', type=float, default=1e-4,\n",
    "                    help='weight decay')\n",
    "parser.add_argument('--clip', type=float, default=10,\n",
    "                    help='gradient clipping')\n",
    "parser.add_argument('--epochs', type=int, default=3500,\n",
    "                    help='upper epoch limit')\n",
    "parser.add_argument('--batch_size', type=int, default=64, metavar='N',\n",
    "                    help='batch size')\n",
    "parser.add_argument('--eval_batch_size', type=int, default=64, metavar='N',\n",
    "                    help='eval_batch size')\n",
    "parser.add_argument('--bptt', type=int, default=50,\n",
    "                    help='sequence length')\n",
    "parser.add_argument('--teacher_forcing_ratio', type=float, default=0.7,\n",
    "                    help='teacher forcing ratio (deprecated)')\n",
    "parser.add_argument('--dropout', type=float, default=0.2,\n",
    "                    help='dropout applied to layers (0 = no dropout)')\n",
    "parser.add_argument('--tied', action='store_true',\n",
    "                    help='tie the word embedding and softmax weights (deprecated)')\n",
    "parser.add_argument('--seed', type=int, default=1111,\n",
    "                    help='random seed')\n",
    "parser.add_argument('--device', type=str, default='cpu',\n",
    "                    help='cuda or cpu')\n",
    "parser.add_argument('--log_interval', type=int, default=10, metavar='N',\n",
    "                    help='report interval')\n",
    "parser.add_argument('--save_interval', type=int, default=10, metavar='N',\n",
    "                    help='save interval')\n",
    "parser.add_argument('--save_fig', action='store_true', default=True)\n",
    "#                     help='save figure')\n",
    "# parser.add_argument('--path_save', type=str, default='./result/nyc_taxi',\n",
    "#                     help='file path to save data')\n",
    "# parser.add_argument('--path_load', type=str, default='./result/nyc_taxi',\n",
    "#                     help='file path to load dataset from')\n",
    "parser.add_argument('--resume','-r',\n",
    "                    help='use checkpoint model parameters as initial parameters (default: False)',\n",
    "                    action=\"store_true\")\n",
    "parser.add_argument('--pretrained','-p',\n",
    "                    help='use checkpoint model parameters and do not train anymore (default: False)',\n",
    "                    action=\"store_true\")\n",
    "parser.add_argument('--prediction_window_size', type=int, default=10,\n",
    "                    help='prediction_window_size')\n",
    "# args = parser.parse_args()\n",
    "args, unknown = parser.parse_known_args()\n",
    "# Set the random seed manually for reproducibility.\n",
    "torch.manual_seed(args.seed)\n",
    "torch.cuda.manual_seed(args.seed)\n",
    "\n",
    "###############################################################################\n",
    "# Load data\n",
    "###############################################################################\n",
    "TimeseriesData = preprocess_data.PickleDataLoad(data_type=args.data, filename=args.filename, \n",
    "                                                augment_test_data=args.augment)\n",
    "train_dataset = TimeseriesData.batchify(args,TimeseriesData.trainData, args.batch_size)\n",
    "test_dataset = TimeseriesData.batchify(args,TimeseriesData.testData, args.eval_batch_size)\n",
    "gen_dataset = TimeseriesData.batchify(args,TimeseriesData.testData, 1)\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "# Build the model\n",
    "###############################################################################\n",
    "feature_dim = TimeseriesData.trainData.size(1)\n",
    "model = model.RNNPredictor(rnn_type = args.model,\n",
    "                           enc_inp_size=feature_dim,\n",
    "                           rnn_inp_size = args.emsize,\n",
    "                           rnn_hid_size = args.nhid,\n",
    "                           dec_out_size=feature_dim,\n",
    "                           nlayers = args.nlayers,\n",
    "                           dropout = args.dropout,\n",
    "                           tie_weights= args.tied,\n",
    "                           res_connection=args.res_connection).to(args.device)\n",
    "optimizer = optim.Adam(model.parameters(), lr= args.lr,weight_decay=args.weight_decay)\n",
    "criterion = nn.MSELoss()\n",
    "###############################################################################\n",
    "# Training code\n",
    "###############################################################################\n",
    "def get_batch(args,source, i):\n",
    "    seq_len = min(args.bptt, len(source) - 1 - i)\n",
    "    data = source[i:i+seq_len] # [ seq_len * batch_size * feature_size ]\n",
    "    target = source[i+1:i+1+seq_len] # [ (seq_len x batch_size x feature_size) ]\n",
    "    return data, target\n",
    "\n",
    "def generate_output(args,epoch, model, gen_dataset, disp_uncertainty=True,startPoint=500, endPoint=3500):\n",
    "    if args.save_fig:\n",
    "        # Turn on evaluation mode which disables dropout.\n",
    "        model.eval()\n",
    "        hidden = model.init_hidden(1)\n",
    "        outSeq = []\n",
    "        upperlim95 = []\n",
    "        lowerlim95 = []\n",
    "        with torch.no_grad():\n",
    "            for i in range(endPoint):\n",
    "                if i>=startPoint:\n",
    "#                     if disp_uncertainty and epoch > 40:\n",
    "#                         outs = []\n",
    "#                         model.train()\n",
    "#                         for i in range(20):\n",
    "#                             out_, hidden_ = model.forward(out+0.01*Variable(torch.randn(out.size())).cuda(),hidden,noise=True)\n",
    "#                             outs.append(out_)\n",
    "#                         model.eval()\n",
    "#                         outs = torch.cat(outs,dim=0)\n",
    "#                         out_mean = torch.mean(outs,dim=0) # [bsz * feature_dim]\n",
    "#                         out_std = torch.std(outs,dim=0) # [bsz * feature_dim]\n",
    "#                         upperlim95.append(out_mean + 2.58*out_std/np.sqrt(20))\n",
    "#                         lowerlim95.append(out_mean - 2.58*out_std/np.sqrt(20))\n",
    "\n",
    "                    out, hidden = model.forward(out, hidden)\n",
    "\n",
    "                    #print(out_mean,out)\n",
    "\n",
    "                else:\n",
    "                    out, hidden = model.forward(gen_dataset[i].unsqueeze(0), hidden)\n",
    "                outSeq.append(out.data.cpu()[0][0].unsqueeze(0))\n",
    "\n",
    "\n",
    "        outSeq = torch.cat(outSeq,dim=0) # [seqLength * feature_dim]\n",
    "\n",
    "        target= preprocess_data.reconstruct(gen_dataset.cpu(), TimeseriesData.mean, TimeseriesData.std)\n",
    "        outSeq = preprocess_data.reconstruct(outSeq, TimeseriesData.mean, TimeseriesData.std)\n",
    "#         if epoch>40:\n",
    "#             upperlim95 = torch.cat(upperlim95, dim=0)\n",
    "#             lowerlim95 = torch.cat(lowerlim95, dim=0)\n",
    "#             upperlim95 = preprocess_data.reconstruct(upperlim95.data.cpu().numpy(),TimeseriesData.mean,TimeseriesData.std)\n",
    "#             lowerlim95 = preprocess_data.reconstruct(lowerlim95.data.cpu().numpy(),TimeseriesData.mean,TimeseriesData.std)\n",
    "\n",
    "        plt.figure(figsize=(15,5))\n",
    "        for i in range(target.size(-1)):\n",
    "            plt.plot(target[:,:,i].numpy(), label='Target'+str(i),\n",
    "                     color='black', marker='.', linestyle='--', markersize=1, linewidth=0.5)\n",
    "            plt.plot(range(startPoint), outSeq[:startPoint,i].numpy(), label='1-step predictions for target'+str(i),\n",
    "                     color='green', marker='.', linestyle='--', markersize=1.5, linewidth=1)\n",
    "#             if epoch>40:\n",
    "#                 plt.plot(range(startPoint, endPoint), upperlim95[:,i].numpy(), label='upperlim'+str(i),\n",
    "#                          color='skyblue', marker='.', linestyle='--', markersize=1.5, linewidth=1)\n",
    "#                 plt.plot(range(startPoint, endPoint), lowerlim95[:,i].numpy(), label='lowerlim'+str(i),\n",
    "#                          color='skyblue', marker='.', linestyle='--', markersize=1.5, linewidth=1)\n",
    "            plt.plot(range(startPoint, endPoint), outSeq[startPoint:,i].numpy(), label='Recursive predictions for target'+str(i),\n",
    "                     color='blue', marker='.', linestyle='--', markersize=1.5, linewidth=1)\n",
    "\n",
    "        plt.xlim([startPoint-500, endPoint])\n",
    "        plt.xlabel('Index',fontsize=15)\n",
    "        plt.ylabel('Value',fontsize=15)\n",
    "        plt.title('Time-series Prediction on ' + args.data + ' Dataset', fontsize=18, fontweight='bold')\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "        plt.text(startPoint-500+10, target.min(), 'Epoch: '+str(epoch),fontsize=15)\n",
    "#         save_dir = Path('result',args.data,args.filename).with_suffix('').joinpath('fig_prediction')\n",
    "#         save_dir.mkdir(parents=True,exist_ok=True)\n",
    "        save_dir = r'/Users/immanuvelprathaps/Desktop/RNN-Time-series-Anomaly-Detection-master/result/ecg/chfdb_chf13_45590 fig_prediction fig_epoch'+ str(epoch)+ '.png'\n",
    "#         plt.savefig(save_dir.joinpath('fig_epoch'+str(epoch)).with_suffix('.png'))\n",
    "        plt.savefig(save_dir)\n",
    "        plt.show()\n",
    "        plt.close()\n",
    "        \n",
    "        return outSeq\n",
    "\n",
    "    else:\n",
    "        pass\n",
    "\n",
    "\n",
    "\n",
    "def evaluate_1step_pred(args, model, test_dataset):\n",
    "    # Turn on evaluation mode which disables dropout.\n",
    "    model.eval()\n",
    "    total_loss = 0\n",
    "    with torch.no_grad():\n",
    "        hidden = model.init_hidden(args.eval_batch_size)\n",
    "        for nbatch, i in enumerate(range(0, test_dataset.size(0) - 1, args.bptt)):\n",
    "\n",
    "            inputSeq, targetSeq = get_batch(args,test_dataset, i)\n",
    "            outSeq, hidden = model.forward(inputSeq, hidden)\n",
    "\n",
    "            loss = criterion(outSeq.contiguous().view(args.batch_size,-1), targetSeq.contiguous().view(args.batch_size,-1))\n",
    "            hidden = model.repackage_hidden(hidden)\n",
    "            total_loss+= loss.item()\n",
    "\n",
    "    return total_loss / nbatch\n",
    "\n",
    "def train(args, model, train_dataset,epoch):\n",
    "\n",
    "    with torch.enable_grad():\n",
    "        # Turn on training mode which enables dropout.\n",
    "        model.train()\n",
    "        total_loss = 0\n",
    "        start_time = time.time()\n",
    "        hidden = model.init_hidden(args.batch_size)\n",
    "        for batch, i in enumerate(range(0, train_dataset.size(0) - 1, args.bptt)):\n",
    "            inputSeq, targetSeq = get_batch(args,train_dataset, i)\n",
    "            # inputSeq: [ seq_len * batch_size * feature_size ]\n",
    "            # targetSeq: [ seq_len * batch_size * feature_size ]\n",
    "\n",
    "            # Starting each batch, we detach the hidden state from how it was previously produced.\n",
    "            # If we didn't, the model would try backpropagating all the way to start of the dataset.\n",
    "            hidden = model.repackage_hidden(hidden)\n",
    "            hidden_ = model.repackage_hidden(hidden)\n",
    "            optimizer.zero_grad()\n",
    "\n",
    "            '''Loss1: Free running loss'''\n",
    "            outVal = inputSeq[0].unsqueeze(0)\n",
    "            outVals=[]\n",
    "            hids1 = []\n",
    "            for i in range(inputSeq.size(0)):\n",
    "                outVal, hidden_, hid = model.forward(outVal, hidden_,return_hiddens=True)\n",
    "                outVals.append(outVal)\n",
    "                hids1.append(hid)\n",
    "            outSeq1 = torch.cat(outVals,dim=0)\n",
    "            hids1 = torch.cat(hids1,dim=0)\n",
    "            loss1 = criterion(outSeq1.contiguous().view(args.batch_size,-1), targetSeq.contiguous().view(args.batch_size,-1))\n",
    "\n",
    "            '''Loss2: Teacher forcing loss'''\n",
    "            outSeq2, hidden, hids2 = model.forward(inputSeq, hidden, return_hiddens=True)\n",
    "            loss2 = criterion(outSeq2.contiguous().view(args.batch_size, -1), targetSeq.contiguous().view(args.batch_size, -1))\n",
    "\n",
    "            '''Loss3: Simplified Professor forcing loss'''\n",
    "            loss3 = criterion(hids1.contiguous().view(args.batch_size,-1), hids2.contiguous().view(args.batch_size,-1).detach())\n",
    "\n",
    "            '''Total loss = Loss1+Loss2+Loss3'''\n",
    "            loss = loss1+loss2+loss3\n",
    "            loss.backward()\n",
    "\n",
    "            # `clip_grad_norm` helps prevent the exploding gradient problem in RNNs / LSTMs.\n",
    "            torch.nn.utils.clip_grad_norm_(model.parameters(), args.clip)\n",
    "            optimizer.step()\n",
    "\n",
    "            total_loss += loss.item()\n",
    "\n",
    "            if batch % args.log_interval == 0 and batch > 0:\n",
    "                cur_loss = total_loss / args.log_interval\n",
    "                elapsed = time.time() - start_time\n",
    "                print('| epoch {:3d} | {:5d}/{:5d} batches | ms/batch {:5.4f} | '\n",
    "                      'loss {:5.2f} '.format(\n",
    "                    epoch, batch, len(train_dataset) // args.bptt,\n",
    "                                  elapsed * 1000 / args.log_interval, cur_loss))\n",
    "                total_loss = 0\n",
    "                start_time = time.time()\n",
    "\n",
    "def evaluate(args, model, test_dataset):\n",
    "    # Turn on evaluation mode which disables dropout.\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        total_loss = 0\n",
    "        hidden = model.init_hidden(args.eval_batch_size)\n",
    "        nbatch = 1\n",
    "        for nbatch, i in enumerate(range(0, test_dataset.size(0) - 1, args.bptt)):\n",
    "            inputSeq, targetSeq = get_batch(args,test_dataset, i)\n",
    "            # inputSeq: [ seq_len * batch_size * feature_size ]\n",
    "            # targetSeq: [ seq_len * batch_size * feature_size ]\n",
    "            hidden_ = model.repackage_hidden(hidden)\n",
    "            '''Loss1: Free running loss'''\n",
    "            outVal = inputSeq[0].unsqueeze(0)\n",
    "            outVals=[]\n",
    "            hids1 = []\n",
    "            for i in range(inputSeq.size(0)):\n",
    "                outVal, hidden_, hid = model.forward(outVal, hidden_,return_hiddens=True)\n",
    "                outVals.append(outVal)\n",
    "                hids1.append(hid)\n",
    "            outSeq1 = torch.cat(outVals,dim=0)\n",
    "            hids1 = torch.cat(hids1,dim=0)\n",
    "            loss1 = criterion(outSeq1.contiguous().view(args.batch_size,-1), targetSeq.contiguous().view(args.batch_size,-1))\n",
    "\n",
    "            '''Loss2: Teacher forcing loss'''\n",
    "            outSeq2, hidden, hids2 = model.forward(inputSeq, hidden, return_hiddens=True)\n",
    "            loss2 = criterion(outSeq2.contiguous().view(args.batch_size, -1), targetSeq.contiguous().view(args.batch_size, -1))\n",
    "\n",
    "            '''Loss3: Simplified Professor forcing loss'''\n",
    "            loss3 = criterion(hids1.contiguous().view(args.batch_size,-1), hids2.contiguous().view(args.batch_size,-1).detach())\n",
    "\n",
    "            '''Total loss = Loss1+Loss2+Loss3'''\n",
    "            loss = loss1+loss2+loss3\n",
    "\n",
    "            total_loss += loss.item()\n",
    "\n",
    "    return total_loss / (nbatch+1)\n",
    "\n",
    "\n",
    "# Loop over epochs.\n",
    "if args.resume or args.pretrained:\n",
    "    print(\"=> loading checkpoint \")\n",
    "    checkpoint = torch.load(Path('save', args.data, 'checkpoint', args.filename).with_suffix('.pth'))\n",
    "    args, start_epoch, best_val_loss = model.load_checkpoint(args,checkpoint,feature_dim)\n",
    "    optimizer.load_state_dict((checkpoint['optimizer']))\n",
    "    del checkpoint\n",
    "    epoch = start_epoch\n",
    "    print(\"=> loaded checkpoint\")\n",
    "else:\n",
    "    epoch = 1\n",
    "    start_epoch = 1\n",
    "    best_val_loss = float('inf')\n",
    "    print(\"=> Start training from scratch\")\n",
    "print('-' * 89)\n",
    "print(args)\n",
    "print('-' * 89)\n",
    "\n",
    "if not args.pretrained:\n",
    "    # At any point you can hit Ctrl + C to break out of training early.\n",
    "    try:\n",
    "        for epoch in range(start_epoch, args.epochs+1):\n",
    "\n",
    "            epoch_start_time = time.time()\n",
    "            train(args,model,train_dataset,epoch)\n",
    "            val_loss = evaluate(args,model,test_dataset)\n",
    "            print('-' * 89)\n",
    "            print('| end of epoch {:3d} | time: {:5.2f}s | valid loss {:5.4f} | '.format(epoch, (time.time() - epoch_start_time),                                                                                        val_loss))\n",
    "            print('-' * 89)\n",
    "\n",
    "            generate_output(args,epoch,model,gen_dataset,startPoint=1500)\n",
    "\n",
    "            if epoch%args.save_interval==0:\n",
    "                # Save the model if the validation loss is the best we've seen so far.\n",
    "                is_best = val_loss < best_val_loss\n",
    "                best_val_loss = min(val_loss, best_val_loss)\n",
    "                model_dictionary = {'epoch': epoch,\n",
    "                                    'best_loss': best_val_loss,\n",
    "                                    'state_dict': model.state_dict(),\n",
    "                                    'optimizer': optimizer.state_dict(),\n",
    "                                    'args':args\n",
    "                                    }\n",
    "                model.save_checkpoint(model_dictionary, is_best)\n",
    "\n",
    "    except KeyboardInterrupt:\n",
    "        print('-' * 89)\n",
    "        print('Exiting from training early')\n",
    "\n",
    "\n",
    "# Calculate mean and covariance for each channel's prediction errors, and save them with the trained model\n",
    "print('=> calculating mean and covariance')\n",
    "means, covs = list(),list()\n",
    "train_dataset = TimeseriesData.batchify(args, TimeseriesData.trainData, bsz=1)\n",
    "for channel_idx in range(model.enc_input_size):\n",
    "    mean, cov = fit_norm_distribution_param(args,model,train_dataset[:TimeseriesData.length],channel_idx)\n",
    "    means.append(mean), covs.append(cov)\n",
    "model_dictionary = {'epoch': max(epoch,start_epoch),\n",
    "                    'best_loss': best_val_loss,\n",
    "                    'state_dict': model.state_dict(),\n",
    "                    'optimizer': optimizer.state_dict(),\n",
    "                    'args': args,\n",
    "                    'means': means,\n",
    "                    'covs': covs\n",
    "                    }\n",
    "model.save_checkpoint(model_dictionary, True)\n",
    "print('-' * 89)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25bc0eeb",
   "metadata": {},
   "source": [
    "# PLOTTING ANOMALY SCORES:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bee47075",
   "metadata": {},
   "outputs": [],
   "source": [
    "import argparse\n",
    "import torch\n",
    "import pickle\n",
    "import preprocess_data\n",
    "from model import model\n",
    "from torch import optim\n",
    "from pathlib import Path\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from anomalyDetector import fit_norm_distribution_param\n",
    "from anomalyDetector import anomalyScore\n",
    "from anomalyDetector import get_precision_recall\n",
    "parser = argparse.ArgumentParser(description='PyTorch RNN Anomaly Detection Model')\n",
    "parser.add_argument('--prediction_window_size', type=int, default=10,\n",
    "                    help='prediction_window_size')\n",
    "parser.add_argument('--data', type=str, default='ecg',\n",
    "                    help='type of the dataset (ecg, gesture, power_demand, space_shuttle, respiration, nyc_taxi')\n",
    "parser.add_argument('--filename', type=str, default='chfdb_chf13_45590.pkl',\n",
    "                    help='filename of the dataset')\n",
    "parser.add_argument('--save_fig', action='store_true', default=True)\n",
    "#                     help='save results as figures')\n",
    "# parser.add_argument('--path_save', type=str, default='./result/nyc_taxi',\n",
    "#                     help='file path to save data')\n",
    "# parser.add_argument('--path_load', type=str, default='./result/nyc_taxi',\n",
    "#                     help='file path to load dataset from')\n",
    "parser.add_argument('--compensate', action='store_true',\n",
    "                    help='compensate anomaly score using anomaly score esimation')\n",
    "parser.add_argument('--beta', type=float, default=1.0,\n",
    "                    help='beta value for f-beta score')\n",
    "\n",
    "\n",
    "# args_ = parser.parse_args()\n",
    "args_, unknown = parser.parse_known_args()\n",
    "\n",
    "print('-' * 89)\n",
    "print(\"=> loading checkpoint \")\n",
    "checkpoint = torch.load(str(Path('save',args_.data,'checkpoint',args_.filename).with_suffix('.pth')))\n",
    "args = checkpoint['args']\n",
    "args.prediction_window_size= args_.prediction_window_size\n",
    "args.beta = args_.beta\n",
    "args.save_fig = args_.save_fig\n",
    "args.compensate = args_.compensate\n",
    "print(\"=> loaded checkpoint\")\n",
    "\n",
    "\n",
    "# Set the random seed manually for reproducibility.\n",
    "torch.manual_seed(args.seed)\n",
    "torch.cuda.manual_seed(args.seed)\n",
    "\n",
    "###############################################################################\n",
    "# Load data\n",
    "###############################################################################\n",
    "TimeseriesData = preprocess_data.PickleDataLoad(data_type=args.data,filename=args.filename, augment_test_data=False)\n",
    "train_dataset = TimeseriesData.batchify(args,TimeseriesData.trainData[:TimeseriesData.length], bsz=1)\n",
    "test_dataset = TimeseriesData.batchify(args,TimeseriesData.testData, bsz=1)\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "# Build the model\n",
    "###############################################################################\n",
    "nfeatures = TimeseriesData.trainData.size(-1)\n",
    "model = model.RNNPredictor(rnn_type = args.model,\n",
    "                           enc_inp_size=nfeatures,\n",
    "                           rnn_inp_size = args.emsize,\n",
    "                           rnn_hid_size = args.nhid,\n",
    "                           dec_out_size=nfeatures,\n",
    "                           nlayers = args.nlayers,\n",
    "                           res_connection=args.res_connection).to(args.device)\n",
    "model.load_state_dict(checkpoint['state_dict'])\n",
    "#del checkpoint\n",
    "\n",
    "scores, predicted_scores, precisions, recalls, f_betas = list(), list(), list(), list(), list()\n",
    "targets, mean_predictions, oneStep_predictions, Nstep_predictions = list(), list(), list(), list()\n",
    "try:\n",
    "    # For each channel in the dataset\n",
    "    for channel_idx in range(nfeatures):\n",
    "        ''' 1. Load mean and covariance if they are pre-calculated, if not calculate them. '''\n",
    "        # Mean and covariance are calculated on train dataset.\n",
    "        if 'means' in checkpoint.keys() and 'covs' in checkpoint.keys():\n",
    "            print('=> loading pre-calculated mean and covariance')\n",
    "            mean, cov = checkpoint['means'][channel_idx], checkpoint['covs'][channel_idx]\n",
    "        else:\n",
    "            print('=> calculating mean and covariance')\n",
    "            mean, cov = fit_norm_distribution_param(args, model, train_dataset, channel_idx=channel_idx)\n",
    "\n",
    "        ''' 2. Train anomaly score predictor using support vector regression (SVR). (Optional) '''\n",
    "        # An anomaly score predictor is trained\n",
    "        # given hidden layer output and the corresponding anomaly score on train dataset.\n",
    "        # Predicted anomaly scores on test dataset can be used for the baseline of the adaptive threshold.\n",
    "        if args.compensate:\n",
    "            print('=> training an SVR as anomaly score predictor')\n",
    "            train_score, _, _, hiddens, _ = anomalyScore(args, model, train_dataset, mean, cov, channel_idx=channel_idx)\n",
    "            score_predictor = GridSearchCV(SVR(), cv=5,param_grid={\"C\": [1e0, 1e1, 1e2],\"gamma\": np.logspace(-1, 1, 3)})\n",
    "            score_predictor.fit(torch.cat(hiddens,dim=0).numpy(), train_score.cpu().numpy())\n",
    "        else:\n",
    "            score_predictor=None\n",
    "\n",
    "        ''' 3. Calculate anomaly scores'''\n",
    "        # Anomaly scores are calculated on the test dataset\n",
    "        # given the mean and the covariance calculated on the train dataset\n",
    "        print('=> calculating anomaly scores')\n",
    "        score, sorted_prediction, sorted_error, _, predicted_score = anomalyScore(args, model, test_dataset, mean, cov,\n",
    "                                                                                  score_predictor=score_predictor,\n",
    "                                                                                  channel_idx=channel_idx)\n",
    "\n",
    "        ''' 4. Evaluate the result '''\n",
    "        # The obtained anomaly scores are evaluated by measuring precision, recall, and f_beta scores\n",
    "        # The precision, recall, f_beta scores are are calculated repeatedly,\n",
    "        # sampling the threshold from 1 to the maximum anomaly score value, either equidistantly or logarithmically.\n",
    "        print('=> calculating precision, recall, and f_beta')\n",
    "        precision, recall, f_beta = get_precision_recall(args, score, num_samples=1000, beta=args.beta,\n",
    "                                                         label=TimeseriesData.testLabel.to(args.device))\n",
    "        print('data: ',args.data,' filename: ',args.filename,\n",
    "              ' f-beta (no compensation): ', f_beta.max().item(),' beta: ',args.beta)\n",
    "        if args.compensate:\n",
    "            precision, recall, f_beta = get_precision_recall(args, score, num_samples=1000, beta=args.beta,\n",
    "                                                             label=TimeseriesData.testLabel.to(args.device),\n",
    "                                                             predicted_score=predicted_score)\n",
    "            print('data: ',args.data,' filename: ',args.filename,\n",
    "                  ' f-beta    (compensation): ', f_beta.max().item(),' beta: ',args.beta)\n",
    "\n",
    "\n",
    "        target = preprocess_data.reconstruct(test_dataset.cpu()[:, 0, channel_idx],\n",
    "                                             TimeseriesData.mean[channel_idx],\n",
    "                                             TimeseriesData.std[channel_idx]).numpy()\n",
    "        mean_prediction = preprocess_data.reconstruct(sorted_prediction.mean(dim=1).cpu(),\n",
    "                                                      TimeseriesData.mean[channel_idx],\n",
    "                                                      TimeseriesData.std[channel_idx]).numpy()\n",
    "        oneStep_prediction = preprocess_data.reconstruct(sorted_prediction[:, -1].cpu(),\n",
    "                                                         TimeseriesData.mean[channel_idx],\n",
    "                                                         TimeseriesData.std[channel_idx]).numpy()\n",
    "        Nstep_prediction = preprocess_data.reconstruct(sorted_prediction[:, 0].cpu(),\n",
    "                                                       TimeseriesData.mean[channel_idx],\n",
    "                                                       TimeseriesData.std[channel_idx]).numpy()\n",
    "        sorted_errors_mean = sorted_error.abs().mean(dim=1).cpu()\n",
    "        sorted_errors_mean *= TimeseriesData.std[channel_idx]\n",
    "        sorted_errors_mean = sorted_errors_mean.numpy()\n",
    "        score = score.cpu()\n",
    "        scores.append(score), targets.append(target), predicted_scores.append(predicted_score)\n",
    "        mean_predictions.append(mean_prediction), oneStep_predictions.append(oneStep_prediction)\n",
    "        Nstep_predictions.append(Nstep_prediction)\n",
    "        precisions.append(precision), recalls.append(recall), f_betas.append(f_beta)\n",
    "\n",
    "\n",
    "        if args.save_fig:\n",
    "            save_dir = Path('result',args.data,args.filename).with_suffix('').joinpath('fig_detection')\n",
    "            save_dir.mkdir(parents=True,exist_ok=True)\n",
    "            plt.plot(precision.cpu().numpy(),label='precision')\n",
    "            plt.plot(recall.cpu().numpy(),label='recall')\n",
    "            plt.plot(f_beta.cpu().numpy(), label='f1')\n",
    "            plt.legend()\n",
    "            plt.xlabel('Threshold (log scale)')\n",
    "            plt.ylabel('Value')\n",
    "            plt.title('Anomaly Detection on ' + args.data + ' Dataset', fontsize=18, fontweight='bold')\n",
    "            plt.savefig(str(save_dir.joinpath('fig_f_beta_channel'+str(channel_idx)).with_suffix('.png')))\n",
    "            plt.close()\n",
    "\n",
    "\n",
    "            fig, ax1 = plt.subplots(figsize=(15,5))\n",
    "            ax1.plot(target,label='Target',\n",
    "                     color='black',  marker='.', linestyle='--', markersize=1, linewidth=0.5)\n",
    "            ax1.plot(mean_prediction, label='Mean predictions',\n",
    "                     color='purple', marker='.', linestyle='--', markersize=1, linewidth=0.5)\n",
    "            ax1.plot(oneStep_prediction, label='1-step predictions',\n",
    "                     color='green', marker='.', linestyle='--', markersize=1, linewidth=0.5)\n",
    "            ax1.plot(Nstep_prediction, label=str(args.prediction_window_size) + '-step predictions',\n",
    "                     color='blue', marker='.', linestyle='--', markersize=1, linewidth=0.5)\n",
    "            ax1.plot(sorted_errors_mean,label='Absolute mean prediction errors',\n",
    "                     color='orange', marker='.', linestyle='--', markersize=1, linewidth=1.0)\n",
    "            ax1.legend(loc='upper left')\n",
    "            ax1.set_ylabel('Value',fontsize=15)\n",
    "            ax1.set_xlabel('Index',fontsize=15)\n",
    "            ax2 = ax1.twinx()\n",
    "            ax2.plot(score.numpy().reshape(-1, 1), label='Anomaly scores from \\nmultivariate normal distribution',\n",
    "                     color='red', marker='.', linestyle='--', markersize=1, linewidth=1)\n",
    "            if args.compensate:\n",
    "                ax2.plot(predicted_score, label='Predicted anomaly scores from SVR',\n",
    "                         color='cyan', marker='.', linestyle='--', markersize=1, linewidth=1)\n",
    "            #ax2.plot(score.reshape(-1,1)/(predicted_score+1),label='Anomaly scores from \\nmultivariate normal distribution',\n",
    "            #        color='hotpink', marker='.', linestyle='--', markersize=1, linewidth=1)\n",
    "            ax2.legend(loc='upper right')\n",
    "            ax2.set_ylabel('anomaly score',fontsize=15)\n",
    "            #plt.axvspan(2830,2900 , color='yellow', alpha=0.3)\n",
    "            plt.title('Anomaly Detection on ' + args.data + ' Dataset', fontsize=18, fontweight='bold')\n",
    "            plt.tight_layout()\n",
    "            plt.xlim([0,len(test_dataset)])\n",
    "            plt.savefig(str(save_dir.joinpath('fig_scores_channel'+str(channel_idx)).with_suffix('.png')))\n",
    "            #plt.show()\n",
    "            plt.close()\n",
    "\n",
    "\n",
    "except KeyboardInterrupt:\n",
    "    print('-' * 89)\n",
    "    print('Exiting from training early')\n",
    "\n",
    "\n",
    "print('=> saving the results as pickle extensions')\n",
    "save_dir = Path('result', args.data, args.filename).with_suffix('')\n",
    "save_dir.mkdir(parents=True, exist_ok=True)\n",
    "pickle.dump(targets, open(str(save_dir.joinpath('target.pkl')),'wb'))\n",
    "pickle.dump(mean_predictions, open(str(save_dir.joinpath('mean_predictions.pkl')),'wb'))\n",
    "pickle.dump(oneStep_predictions, open(str(save_dir.joinpath('oneStep_predictions.pkl')),'wb'))\n",
    "pickle.dump(Nstep_predictions, open(str(save_dir.joinpath('Nstep_predictions.pkl')),'wb'))\n",
    "pickle.dump(scores, open(str(save_dir.joinpath('score.pkl')),'wb'))\n",
    "pickle.dump(predicted_scores, open(str(save_dir.joinpath('predicted_scores.pkl')),'wb'))\n",
    "pickle.dump(precisions, open(str(save_dir.joinpath('precision.pkl')),'wb'))\n",
    "pickle.dump(recalls, open(str(save_dir.joinpath('recall.pkl')),'wb'))\n",
    "pickle.dump(f_betas, open(str(save_dir.joinpath('f_beta.pkl')),'wb'))\n",
    "print('-' * 89)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54bdc251",
   "metadata": {},
   "source": [
    "# CONVERTING PNG TO GIF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45a27af9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio\n",
    "import glob\n",
    "import os\n",
    "import argparse\n",
    "import re\n",
    "\n",
    "def tryint(s):\n",
    "    try:\n",
    "        return int(s)\n",
    "    except:\n",
    "        return s\n",
    "\n",
    "def alphanum_key(s):\n",
    "    \"\"\" Turn a string into a list of string and number chunks.\n",
    "        \"z23a\" -> [\"z\", 23, \"a\"]\n",
    "    \"\"\"\n",
    "    return [ tryint(c) for c in re.split('([0-9]+)', s) ]\n",
    "\n",
    "parser = argparse.ArgumentParser(description='PNG to GIF')\n",
    "parser.add_argument('--path', type=str, default='/Users/immanuvelprathaps/Desktop/RNN-Time-series-Anomaly-Detection-master/result/ecg',\n",
    "                    help='file path ')\n",
    "\n",
    "# args = parser.parse_args()\n",
    "args, unknown = parser.parse_known_args()\n",
    "\n",
    "\n",
    "images = []\n",
    "filenames = glob.glob(args.path+'/*')\n",
    "\n",
    "filenames.sort(key=os.path.getmtime)\n",
    "#filenames.sort(key=alphanum_key)\n",
    "#print(filenames)\n",
    "for filename in filenames:\n",
    "    images.append(imageio.imread(filename))\n",
    "imageio.mimsave(args.path+'/fig.gif',images,duration=0.2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}