LSTM_time-series-global_nDays6_k-day3-Copy3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "5OAyOd_Q4BGW"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings\n",
    "\n",
    "import tensorflow as tf\n",
    "from tensorflow.keras import Sequential\n",
    "from tensorflow.keras.layers import LSTM, Dense\n",
    "from tensorflow.keras.layers import Bidirectional\n",
    "from tensorflow.keras.optimizers import Adam\n",
    "from tensorflow.keras.layers import Flatten\n",
    "from tensorflow.compat.v1.keras.layers import TimeDistributed\n",
    "from tensorflow.keras.layers import Conv1D\n",
    "from tensorflow.keras.layers import MaxPooling1D\n",
    "from tensorflow.keras.layers import ConvLSTM2D\n",
    "\n",
    "# This is not recommended but I am doing this to suppress warnings from SARIMAX\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "countryName = 'Brazil'\n",
    "\n",
    "nFeatures = 1\n",
    "\n",
    "nDaysMin = 6\n",
    "k = 3\n",
    "\n",
    "nValid = 10\n",
    "nTest = 10\n",
    "\n",
    "dataDir = os.path.join('data', 'JHU', 'upto07082020_forPublication')\n",
    "\n",
    "\n",
    "# confirmedFilename = 'time_series_covid19_confirmed_global.csv'\n",
    "# deathsFilename = 'time_series_covid19_deaths_global.csv'\n",
    "# recoveredFilename = 'time_series_covid19_recovered_global.csv'\n",
    "\n",
    "confirmedFilename = 'https://raw.githubusercontent.com/arkobarman/covid-19_timeSeriesAnalysis/master/data/JHU/upto07082020_forPublication/time_series_covid19_confirmed_global.csv'\n",
    "deathsFilename = 'https://raw.githubusercontent.com/arkobarman/covid-19_timeSeriesAnalysis/master/data/JHU/upto07082020_forPublication/time_series_covid19_deaths_global.csv'\n",
    "recoveredFilename = 'https://raw.githubusercontent.com/arkobarman/covid-19_timeSeriesAnalysis/master/data/JHU/upto07082020_forPublication/time_series_covid19_recovered_global.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "JgCC9bxv4BHZ"
   },
   "outputs": [],
   "source": [
    "# split a univariate sequence into samples\n",
    "def split_sequence(sequence, n_steps, k):\n",
    "    X, y = list(), list()\n",
    "    for i in range(len(sequence)):\n",
    "        # find the end of this pattern\n",
    "        end_ix = i + n_steps\n",
    "        # check if we are beyond the sequence\n",
    "        if end_ix + k >= len(sequence):\n",
    "            break\n",
    "        # gather input and output parts of the pattern\n",
    "        seq_x, seq_y = sequence[i:end_ix], sequence[end_ix:end_ix+k]\n",
    "        X.append(seq_x)\n",
    "        y.append(seq_y)\n",
    "    return np.array(X), np.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "74VgtVGK4BIf"
   },
   "outputs": [],
   "source": [
    "def meanAbsolutePercentageError(yTrueList, yPredList):\n",
    "    absErrorList = [np.abs(yTrue - yPred) for yTrue, yPred in zip(yTrueList, yPredList)]\n",
    "    absPcErrorList = [absError/yTrue for absError, yTrue in zip(absErrorList, yTrueList)]\n",
    "    MAPE = 100*np.mean(absPcErrorList)\n",
    "    return MAPE\n",
    "\n",
    "def meanAbsolutePercentageError_kDay(yTrueListList, yPredListList):\n",
    "    # Store true and predictions for day 1 in a list, day 2 in a list and so on\n",
    "    # Keep each list of these lists in a respective dict with key as day #\n",
    "    yTrueForDayK = {}\n",
    "    yPredForDayK = {}\n",
    "    for i in range(len(yTrueListList[0])):\n",
    "        yTrueForDayK[i] = []\n",
    "        yPredForDayK[i] = []\n",
    "    for yTrueList, yPredList in zip(yTrueListList, yPredListList):\n",
    "        for i in range(len(yTrueList)):\n",
    "            yTrueForDayK[i].append(yTrueList[i])\n",
    "            yPredForDayK[i].append(yPredList[i])\n",
    "            \n",
    "    # Get MAPE for each day in a list\n",
    "    MAPEList = []\n",
    "    for i in yTrueForDayK.keys():\n",
    "        MAPEList.append(meanAbsolutePercentageError(yTrueForDayK[i], yPredForDayK[i]))\n",
    "    return np.mean(MAPEList)\n",
    "\n",
    "def meanForecastError(yTrueList, yPredList):\n",
    "    forecastErrors = [yTrue - yPred for yTrue, yPred in zip(yTrueList, yPredList)]\n",
    "    MFE = np.mean(forecastErrors)\n",
    "    return MFE\n",
    "\n",
    "def meanAbsoluteError(yTrueList, yPredList):\n",
    "    absErrorList = [np.abs(yTrue - yPred) for yTrue, yPred in zip(yTrueList, yPredList)]\n",
    "    return np.mean(absErrorList)\n",
    "\n",
    "def meanSquaredError(yTrueList, yPredList):\n",
    "    sqErrorList = [np.square(yTrue - yPred) for yTrue, yPred in zip(yTrueList, yPredList)]\n",
    "    return np.mean(sqErrorList)\n",
    "\n",
    "def rootMeanSquaredError(yTrueList, yPredList):\n",
    "    return np.sqrt(meanSquaredError(yTrueList, yPredList))\n",
    "\n",
    "def medianSymmetricAccuracy(yTrueList, yPredList):\n",
    "    '''https://helda.helsinki.fi//bitstream/handle/10138/312261/2017SW001669.pdf?sequence=1'''\n",
    "    logAccRatioList = [np.abs(np.log(yPred/yTrue)) for yTrue, yPred in zip(yTrueList, yPredList)]\n",
    "    MdSA = 100*(np.exp(np.median(logAccRatioList))-1)\n",
    "    return MdSA\n",
    "\n",
    "def medianSymmetricAccuracy_kDay(yTrueListList, yPredListList):\n",
    "    # Store true and predictions for day 1 in a list, day 2 in a list and so on\n",
    "    # Keep each list of these lists in a respective dict with key as day #\n",
    "    yTrueForDayK = {}\n",
    "    yPredForDayK = {}\n",
    "    for i in range(len(yTrueListList[0])):\n",
    "        yTrueForDayK[i] = []\n",
    "        yPredForDayK[i] = []\n",
    "    for yTrueList, yPredList in zip(yTrueListList, yPredListList):\n",
    "        for i in range(len(yTrueList)):\n",
    "            yTrueForDayK[i].append(yTrueList[i])\n",
    "            yPredForDayK[i].append(yPredList[i])\n",
    "    # Get MdSA for each day in a list\n",
    "    MdSAList = []\n",
    "    for i in yTrueForDayK.keys():\n",
    "        MdSAList.append(medianSymmetricAccuracy(yTrueForDayK[i], yPredForDayK[i]))\n",
    "    return(np.mean(MdSAList))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-qIIj_j04BJI"
   },
   "outputs": [],
   "source": [
    "# Function to get all three frames for a given country\n",
    "def getCountryCovidFrDict(countryName):\n",
    "    countryCovidFrDict = {}\n",
    "    for key in covidFrDict.keys():\n",
    "        dataFr = covidFrDict[key]\n",
    "        countryCovidFrDict[key] = dataFr[dataFr['Country/Region'] == countryName]\n",
    "    return countryCovidFrDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 90
    },
    "colab_type": "code",
    "id": "UuatFXNb4BJn",
    "outputId": "c24c03d6-888e-40a2-ce89-73727964626f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total size of data points for LSTM: 109\n",
      "Size of training set: 89\n",
      "Size of validation set: 10\n",
      "Size of test set: 10\n"
     ]
    }
   ],
   "source": [
    "# Load all 3 csv files\n",
    "covidFrDict = {}\n",
    "# covidFrDict['confirmed'] = pd.read_csv(os.path.join(dataDir, confirmedFilename))\n",
    "# covidFrDict['deaths'] = pd.read_csv(os.path.join(dataDir, deathsFilename))\n",
    "covidFrDict['confirmed'] = pd.read_csv(confirmedFilename)\n",
    "covidFrDict['deaths'] = pd.read_csv(deathsFilename)\n",
    "\n",
    "# Recovered is back again!\n",
    "covidFrDict['recovered'] = pd.read_csv(recoveredFilename)\n",
    "\n",
    "countryCovidFrDict = getCountryCovidFrDict(countryName)\n",
    "\n",
    "# Get list of dates\n",
    "colNamesList = list(countryCovidFrDict['confirmed'])\n",
    "dateList = [colName for colName in colNamesList if '/20' in colName]\n",
    "dataList = [countryCovidFrDict['confirmed'][date].iloc[0] for date in dateList]\n",
    "dataDict = dict(zip(dateList, dataList))\n",
    "\n",
    "# Get time series for cases > 100 only\n",
    "daysSince = 100\n",
    "nCasesGreaterDaysSinceList = []\n",
    "datesGreaterDaysSinceList = []\n",
    "\n",
    "for key in dataDict.keys():\n",
    "    if dataDict[key] > daysSince:\n",
    "        datesGreaterDaysSinceList.append(key)\n",
    "        nCasesGreaterDaysSinceList.append(dataDict[key])\n",
    "        \n",
    "XList, yList = split_sequence(nCasesGreaterDaysSinceList, nDaysMin, k)\n",
    "\n",
    "XTrainList = XList[0:len(XList)-(nValid + nTest)]\n",
    "XValidList = XList[len(XList)-(nValid+nTest):len(XList)-(nTest)]\n",
    "XTestList = XList[-nTest:]\n",
    "\n",
    "yTrain = yList[0:len(XList)-(nValid + nTest)]\n",
    "yValid = yList[len(XList)-(nValid+nTest):len(XList)-(nTest)]\n",
    "yTest = yList[-nTest:]\n",
    "\n",
    "print('Total size of data points for LSTM:', len(yList))\n",
    "print('Size of training set:', len(yTrain))\n",
    "print('Size of validation set:', len(yValid))\n",
    "print('Size of test set:', len(yTest))\n",
    "\n",
    "# Convert from list to matrix\n",
    "XTrain = XTrainList.reshape((XTrainList.shape[0], XTrainList.shape[1], nFeatures))\n",
    "XValid = XValidList.reshape((XValidList.shape[0], XValidList.shape[1], nFeatures))\n",
    "XTest = XTestList.reshape((XTestList.shape[0], XTestList.shape[1], nFeatures))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oxqiVdp94BKe"
   },
   "source": [
    "# Vanilla LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "JESH7DG74BLY",
    "outputId": "d080b617-51ef-41f6-de96-78896c0858f4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 2.7465771581602443\n",
      "Updating best MAPE to 2.7465771581602443...\n",
      "Updating best seed to 0...\n",
      "1 2.4005917208376797\n",
      "Updating best MAPE to 2.4005917208376797...\n",
      "Updating best seed to 1...\n",
      "2 4.870962650723183\n",
      "3 99.97409002604861\n",
      "4 2.7482014585024035\n",
      "5 99.98357815080105\n",
      "6 1.3225875552823145\n",
      "Updating best MAPE to 1.3225875552823145...\n",
      "Updating best seed to 6...\n",
      "7 99.9737411453922\n",
      "8 99.99455342341963\n",
      "9 1.6915995215819433\n",
      "10 99.97382974081661\n",
      "11 2.3219817489250882\n",
      "12 99.97374537203949\n",
      "13 1.2997722516456278\n",
      "Updating best MAPE to 1.2997722516456278...\n",
      "Updating best seed to 13...\n",
      "14 4.881652515628713\n",
      "15 5.674132819117169\n",
      "16 1.3768077043821447\n",
      "17 1.838084974102778\n",
      "18 2.6167475972208454\n",
      "19 3.2192797379328204\n",
      "20 99.97638832046702\n",
      "21 0.9878167266453056\n",
      "Updating best MAPE to 0.9878167266453056...\n",
      "Updating best seed to 21...\n",
      "22 99.97651448268873\n",
      "23 4.27264006068046\n",
      "24 99.92172407286772\n",
      "25 2.3696072621945454\n",
      "26 6.033758770778522\n",
      "27 99.97374376652567\n",
      "28 3.10448586809729\n",
      "29 5.598216826913188\n",
      "30 1.8437595712337904\n",
      "31 99.97415988676848\n",
      "32 2.4130788872426585\n",
      "33 1.8685739538712351\n",
      "34 99.97377963876848\n",
      "35 1.0410605180348185\n",
      "36 1.395944526483304\n",
      "37 1.0032097680089673\n",
      "38 99.97494174695406\n",
      "39 99.97442573982221\n",
      "40 1.8542131351713553\n",
      "41 1.490270869481048\n",
      "42 6.694908266111959\n",
      "43 99.9893619433678\n",
      "44 1.6441388711860931\n",
      "45 1.4409549313459378\n",
      "46 1.3476303723060574\n",
      "47 2.1800956442298562\n",
      "48 1.7597169021859298\n",
      "49 1.3361960701500628\n",
      "50 2.3395026568327073\n",
      "51 2.489402427081792\n",
      "52 99.97362568752344\n",
      "53 99.9742466540775\n",
      "54 1.338096294672864\n",
      "55 1.0782066864679676\n",
      "56 1.700002556479263\n",
      "57 4.1691952138443495\n",
      "58 4.496498396599427\n",
      "59 99.97345323808803\n",
      "60 99.97417290617774\n",
      "61 1.241970497488867\n",
      "62 99.97418580593866\n",
      "63 1.4585887181990589\n",
      "64 2.306390450423221\n",
      "65 1.1894433536349762\n",
      "66 3.7846800711078075\n",
      "67 8.320518011483177\n",
      "68 1.0685667911227175\n",
      "69 99.97448768705634\n",
      "70 1.450615787425965\n",
      "71 1.9751097334320065\n",
      "72 99.97368699477802\n",
      "73 99.97398048103405\n",
      "74 99.97341399343141\n",
      "75 1.0094382872060192\n",
      "76 2.268241608329294\n",
      "77 1.6620595272754795\n",
      "78 99.9780795887546\n",
      "79 1.4754370989657446\n",
      "80 1.812542428460647\n",
      "81 99.9733090673994\n",
      "82 2.9632158520958334\n",
      "83 99.97672349025807\n",
      "84 1.1604934612992006\n",
      "85 99.98032762391256\n",
      "86 1.1289490249958678\n",
      "87 99.97333856488979\n",
      "88 1.8482210233600487\n",
      "89 99.97642003284302\n",
      "90 2.798439210540081\n",
      "91 3.3801508397782953\n",
      "92 1.1708950970598968\n",
      "93 2.7498061415038375\n",
      "94 4.0372502379585455\n",
      "95 3.3825108925377614\n",
      "96 1.9712749741523463\n",
      "97 1.2151685989444891\n",
      "98 7.5058651222763375\n",
      "99 1.3265464929755566\n",
      "Training model with best seed...\n"
     ]
    }
   ],
   "source": [
    "nNeurons = 100\n",
    "nFeatures = 1\n",
    "\n",
    "bestValidMAPE = 100\n",
    "bestSeed = -1\n",
    "for seed in range(100):\n",
    "    tf.random.set_seed(seed=seed)\n",
    "    \n",
    "    # define model\n",
    "    model = Sequential()\n",
    "    model.add(LSTM(nNeurons, activation='relu', input_shape=(nDaysMin, nFeatures)))\n",
    "    model.add(Dense(1))\n",
    "    opt = Adam(learning_rate=0.1)\n",
    "    model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "    # fit model\n",
    "    model.fit(XTrain, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "    yPredListList = []\n",
    "    for day in range(nTest):\n",
    "        yPredListList.append([])\n",
    "    XValidNew = XValid.copy()\n",
    "    for day in range(k):\n",
    "        yPred = model.predict(np.float32(XValidNew), verbose=0)\n",
    "        for i in range(len(yPred)):\n",
    "            yPredListList[i].append(yPred[i][0])\n",
    "        XValidNew = np.delete(XValidNew, 0, axis=1)\n",
    "        yPred = np.expand_dims(yPred, 2)\n",
    "        XValidNew = np.append(XValidNew, yPred, axis=1)\n",
    "\n",
    "#     for yTrue, yPred in zip(yTest, yPredList):\n",
    "#         print(yTrue, yPred)\n",
    "\n",
    "    MAPE = meanAbsolutePercentageError_kDay(yValid, yPredListList)\n",
    "    print(seed, MAPE)\n",
    "    if MAPE < bestValidMAPE:\n",
    "        print('Updating best MAPE to {}...'.format(MAPE))\n",
    "        bestValidMAPE = MAPE\n",
    "        print('Updating best seed to {}...'.format(seed))\n",
    "        bestSeed = seed\n",
    "\n",
    "# define model\n",
    "print('Training model with best seed...')\n",
    "tf.random.set_seed(seed=bestSeed)\n",
    "model = Sequential()\n",
    "model.add(LSTM(nNeurons, activation='relu', input_shape=(nDaysMin, nFeatures)))\n",
    "model.add(Dense(1))\n",
    "opt = Adam(learning_rate=0.1)\n",
    "model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "# fit model\n",
    "model.fit(XTrain, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "yPredListList = []\n",
    "for day in range(nTest):\n",
    "    yPredListList.append([])\n",
    "XTestNew = XTest.copy()\n",
    "for day in range(k):\n",
    "    yPred = model.predict(np.float32(XTestNew), verbose=0)\n",
    "    for i in range(len(yPred)):\n",
    "        yPredListList[i].append(yPred[i][0])\n",
    "    XTestNew = np.delete(XTestNew, 0, axis=1)\n",
    "    yPred = np.expand_dims(yPred, 2)\n",
    "    XTestNew = np.append(XTestNew, yPred, axis=1)\n",
    "    \n",
    "MAPE = meanAbsolutePercentageError_kDay(yTest, yPredListList)\n",
    "print('Test MAPE:', MAPE)\n",
    "MdSA = medianSymmetricAccuracy_kDay(yTest, yPredListList)\n",
    "print('Test MdSA:', MdSA)\n",
    "yPredVanillaLSTM = yPredListList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9kS8Y8bT4BMT"
   },
   "source": [
    "# Stacked LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "CREiWf-b4BMZ",
    "outputId": "09d9a2e2-bba5-4f10-a132-9d8a669d2c73"
   },
   "outputs": [],
   "source": [
    "# define model\n",
    "nNeurons = 50\n",
    "nFeatures = 1\n",
    "\n",
    "bestValidMAPE = 100\n",
    "bestSeed = -1\n",
    "for seed in range(100):\n",
    "    tf.random.set_seed(seed=seed)\n",
    "    model = Sequential()\n",
    "    model.add(LSTM(nNeurons, activation='relu', return_sequences=True, input_shape=(nDaysMin, nFeatures)))\n",
    "    model.add(LSTM(nNeurons, activation='relu'))\n",
    "    model.add(Dense(1))\n",
    "    opt = Adam(learning_rate=0.1)\n",
    "    model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "    # fit model\n",
    "    model.fit(XTrain, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "    yPredListList = []\n",
    "    for day in range(nTest):\n",
    "        yPredListList.append([])\n",
    "    XValidNew = XValid.copy()\n",
    "    for day in range(k):\n",
    "        yPred = model.predict(np.float32(XValidNew), verbose=0)\n",
    "        for i in range(len(yPred)):\n",
    "            yPredListList[i].append(yPred[i][0])\n",
    "        XValidNew = np.delete(XValidNew, 0, axis=1)\n",
    "        yPred = np.expand_dims(yPred, 2)\n",
    "        XValidNew = np.append(XValidNew, yPred, axis=1)\n",
    "\n",
    "#     for yTrue, yPred in zip(yTest, yPredList):\n",
    "#         print(yTrue, yPred)\n",
    "\n",
    "    MAPE = meanAbsolutePercentageError_kDay(yValid, yPredListList)\n",
    "    print(seed, MAPE)\n",
    "    if MAPE < bestValidMAPE:\n",
    "        print('Updating best MAPE to {}...'.format(MAPE))\n",
    "        bestValidMAPE = MAPE\n",
    "        print('Updating best seed to {}...'.format(seed))\n",
    "        bestSeed = seed\n",
    "\n",
    "# define model\n",
    "print('Training model with best seed...')\n",
    "tf.random.set_seed(seed=bestSeed)\n",
    "model = Sequential()\n",
    "model.add(LSTM(nNeurons, activation='relu', return_sequences=True, input_shape=(nDaysMin, nFeatures)))\n",
    "model.add(LSTM(nNeurons, activation='relu'))\n",
    "model.add(Dense(1))\n",
    "opt = Adam(learning_rate=0.1)\n",
    "model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "# fit model\n",
    "model.fit(XTrain, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "yPredListList = []\n",
    "for day in range(nTest):\n",
    "    yPredListList.append([])\n",
    "XTestNew = XTest.copy()\n",
    "for day in range(k):\n",
    "    yPred = model.predict(np.float32(XTestNew), verbose=0)\n",
    "    for i in range(len(yPred)):\n",
    "        yPredListList[i].append(yPred[i][0])\n",
    "    XTestNew = np.delete(XTestNew, 0, axis=1)\n",
    "    yPred = np.expand_dims(yPred, 2)\n",
    "    XTestNew = np.append(XTestNew, yPred, axis=1)\n",
    "    \n",
    "MAPE = meanAbsolutePercentageError_kDay(yTest, yPredListList)\n",
    "print('Test MAPE:', MAPE)\n",
    "MdSA = medianSymmetricAccuracy_kDay(yTest, yPredListList)\n",
    "print('Test MdSA:', MdSA)\n",
    "yPredStackedLSTM = yPredListList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "H87S7mZr4BO7"
   },
   "source": [
    "# Bidirectional LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "FUj1hvXg4BSB",
    "outputId": "71323a85-520a-4a2e-d3af-bd6c2af85079"
   },
   "outputs": [],
   "source": [
    "# define model\n",
    "nNeurons = 50\n",
    "nFeatures = 1\n",
    "\n",
    "bestValidMAPE = 100\n",
    "bestSeed = -1\n",
    "for seed in range(100):\n",
    "    tf.random.set_seed(seed=seed)\n",
    "    model = Sequential()\n",
    "    model.add(Bidirectional(LSTM(nNeurons, activation='relu'), input_shape=(nDaysMin, nFeatures)))\n",
    "    model.add(Dense(1))\n",
    "    opt = Adam(learning_rate=0.1)\n",
    "    model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "    # fit model\n",
    "    model.fit(XTrain, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "    yPredListList = []\n",
    "    for day in range(nTest):\n",
    "        yPredListList.append([])\n",
    "    XValidNew = XValid.copy()\n",
    "    for day in range(k):\n",
    "        yPred = model.predict(np.float32(XValidNew), verbose=0)\n",
    "        for i in range(len(yPred)):\n",
    "            yPredListList[i].append(yPred[i][0])\n",
    "        XValidNew = np.delete(XValidNew, 0, axis=1)\n",
    "        yPred = np.expand_dims(yPred, 2)\n",
    "        XValidNew = np.append(XValidNew, yPred, axis=1)\n",
    "\n",
    "#     for yTrue, yPred in zip(yTest, yPredList):\n",
    "#         print(yTrue, yPred)\n",
    "\n",
    "    MAPE = meanAbsolutePercentageError_kDay(yValid, yPredListList)\n",
    "    print(seed, MAPE)\n",
    "    if MAPE < bestValidMAPE:\n",
    "        print('Updating best MAPE to {}...'.format(MAPE))\n",
    "        bestValidMAPE = MAPE\n",
    "        print('Updating best seed to {}...'.format(seed))\n",
    "        bestSeed = seed\n",
    "\n",
    "# define model\n",
    "print('Training model with best seed...')\n",
    "tf.random.set_seed(seed=bestSeed)\n",
    "model = Sequential()\n",
    "model.add(Bidirectional(LSTM(nNeurons, activation='relu'), input_shape=(nDaysMin, nFeatures)))\n",
    "model.add(Dense(1))\n",
    "opt = Adam(learning_rate=0.1)\n",
    "model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "# fit model\n",
    "model.fit(XTrain, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "yPredListList = []\n",
    "for day in range(nTest):\n",
    "    yPredListList.append([])\n",
    "XTestNew = XTest.copy()\n",
    "for day in range(k):\n",
    "    yPred = model.predict(np.float32(XTestNew), verbose=0)\n",
    "    for i in range(len(yPred)):\n",
    "        yPredListList[i].append(yPred[i][0])\n",
    "    XTestNew = np.delete(XTestNew, 0, axis=1)\n",
    "    yPred = np.expand_dims(yPred, 2)\n",
    "    XTestNew = np.append(XTestNew, yPred, axis=1)\n",
    "    \n",
    "MAPE = meanAbsolutePercentageError_kDay(yTest, yPredListList)\n",
    "print('Test MAPE:', MAPE)\n",
    "MdSA = medianSymmetricAccuracy_kDay(yTest, yPredListList)\n",
    "print('Test MdSA:', MdSA)\n",
    "yPredBidirectionalLSTM = yPredListList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1KtwKgmp4BSw"
   },
   "source": [
    "# CNN LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "d0MujnCm4BU-",
    "outputId": "d363481a-6d5b-4b49-9ce5-b83aca35aeb0"
   },
   "outputs": [],
   "source": [
    "# Number of subsequences to break X into (we do 15 = 5x3, 5 subsequences of size 3 each)\n",
    "nSeq = 3\n",
    "nSteps = 2\n",
    "\n",
    "# define model\n",
    "nNeurons = 50\n",
    "nFeatures = 1\n",
    "nFilters = 64\n",
    "\n",
    "bestValidMAPE = 100\n",
    "bestSeed = -1\n",
    "\n",
    "# Reshape input\n",
    "XTrainCNN = XTrainList.reshape((XTrainList.shape[0], nSeq, nSteps, nFeatures))\n",
    "\n",
    "# print(XTrainCNN.shape)\n",
    "# print(XValidCNN.shape)\n",
    "# print(XTestCNN.shape)\n",
    "\n",
    "for seed in range(100):\n",
    "    tf.random.set_seed(seed=seed)\n",
    "    model = Sequential()\n",
    "    model.add(TimeDistributed(Conv1D(filters=nFilters, kernel_size=1, activation='relu'), input_shape=(None, nSteps, nFeatures)))\n",
    "    model.add(TimeDistributed(MaxPooling1D(pool_size=2)))\n",
    "    model.add(TimeDistributed(Flatten()))\n",
    "    model.add(LSTM(nNeurons, activation='relu'))\n",
    "    model.add(Dense(1))\n",
    "    opt = Adam(learning_rate=0.1)\n",
    "    model.compile(optimizer=opt, loss='mse')\n",
    "\n",
    "    # fit model\n",
    "    model.fit(XTrainCNN, yTrain[:,0], epochs=1000, verbose=0)\n",
    "    \n",
    "    yPredListList = []\n",
    "    for day in range(nTest):\n",
    "        yPredListList.append([])\n",
    "    XValidListNew = XValidList.copy()\n",
    "    for day in range(k):\n",
    "        XValidCNN = XValidListNew.reshape((XValidListNew.shape[0], nSeq, nSteps, nFeatures))\n",
    "        yPred = model.predict(np.float32(XValidCNN), verbose=0)\n",
    "        for i in range(len(yPred)):\n",
    "            yPredListList[i].append(yPred[i][0])\n",
    "        XValidListNew = np.delete(XValidListNew, 0, axis=1)\n",
    "        XValidListNew = np.append(XValidListNew, yPred, axis=1)\n",
    "\n",
    "#     for yTrue, yPred in zip(yTest, yPredList):\n",
    "#         print(yTrue, yPred)\n",
    "\n",
    "    MAPE = meanAbsolutePercentageError_kDay(yValid, yPredListList)\n",
    "    print(seed, MAPE)\n",
    "    if MAPE < bestValidMAPE:\n",
    "        print('Updating best MAPE to {}...'.format(MAPE))\n",
    "        bestValidMAPE = MAPE\n",
    "        print('Updating best seed to {}...'.format(seed))\n",
    "        bestSeed = seed\n",
    "\n",
    "# define model\n",
    "print('Training model with best seed...')\n",
    "tf.random.set_seed(seed=bestSeed)\n",
    "model = Sequential()\n",
    "model.add(TimeDistributed(Conv1D(filters=nFilters, kernel_size=1, activation='relu'), input_shape=(None, nSteps, nFeatures)))\n",
    "model.add(TimeDistributed(MaxPooling1D(pool_size=2)))\n",
    "model.add(TimeDistributed(Flatten()))\n",
    "model.add(LSTM(nNeurons, activation='relu'))\n",
    "model.add(Dense(1))\n",
    "opt = Adam(learning_rate=0.1)\n",
    "model.compile(optimizer=opt, loss='mse')\n",
    "# fit model\n",
    "model.fit(XTrainCNN, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "yPredListList = []\n",
    "for day in range(nTest):\n",
    "    yPredListList.append([])\n",
    "XTestListNew = XTestList.copy()\n",
    "for day in range(k):\n",
    "    XTestCNN = XTestListNew.reshape((XTestListNew.shape[0], nSeq, nSteps, nFeatures))\n",
    "    yPred = model.predict(np.float32(XTestCNN), verbose=0)\n",
    "    for i in range(len(yPred)):\n",
    "        yPredListList[i].append(yPred[i][0])\n",
    "    XTestListNew = np.delete(XTestListNew, 0, axis=1)\n",
    "    XTestListNew = np.append(XTestListNew, yPred, axis=1)\n",
    "    \n",
    "MAPE = meanAbsolutePercentageError_kDay(yTest, yPredListList)\n",
    "print('Test MAPE:', MAPE)\n",
    "MdSA = medianSymmetricAccuracy_kDay(yTest, yPredListList)\n",
    "print('Test MdSA:', MdSA)\n",
    "yPredCNNLSTM = yPredListList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FFW0UaoL4BWT"
   },
   "source": [
    "# ConvLSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "_1sG2Kbj4BWe",
    "outputId": "f727a035-66c3-4999-ebbf-77e8915e5da0"
   },
   "outputs": [],
   "source": [
    "# Number of subsequences to break X into (we do 15 = 5x3, 5 subsequences of size 3 each)\n",
    "nSeq = 3\n",
    "nSteps = 2\n",
    "# Each input is rows x columns, we have rows=1 and columns=nSteps\n",
    "\n",
    "# define model\n",
    "nNeurons = 50\n",
    "nFeatures = 1\n",
    "nFilters = 64\n",
    "\n",
    "bestValidMAPE = 100\n",
    "bestSeed = -1\n",
    "\n",
    "# Reshape input\n",
    "XTrainConv = XTrainList.reshape((XTrainList.shape[0], nSeq, 1, nSteps, nFeatures))\n",
    "# XValidConv = XValidList.reshape((XValidList.shape[0], nSeq, 1, nSteps, nFeatures))\n",
    "# XTestConv = XTestList.reshape((XTestList.shape[0], nSeq, 1, nSteps, nFeatures))\n",
    "\n",
    "for seed in range(100):\n",
    "    tf.random.set_seed(seed=seed)\n",
    "    model = Sequential()\n",
    "    model.add(ConvLSTM2D(filters=64, kernel_size=(1,2), activation='relu', input_shape=(nSeq, 1, nSteps, nFeatures)))\n",
    "    model.add(Flatten())\n",
    "    model.add(Dense(1))\n",
    "    opt = Adam(learning_rate=0.1)\n",
    "    model.compile(optimizer=opt, loss='mse')\n",
    "    \n",
    "    # fit model\n",
    "    model.fit(XTrainConv, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "    yPredListList = []\n",
    "    for day in range(nTest):\n",
    "        yPredListList.append([])\n",
    "    XValidListNew = XValidList.copy()\n",
    "    for day in range(k):\n",
    "        XValidConv = XValidListNew.reshape((XValidListNew.shape[0], nSeq, 1, nSteps, nFeatures))\n",
    "        yPred = model.predict(np.float32(XValidConv), verbose=0)\n",
    "        for i in range(len(yPred)):\n",
    "            yPredListList[i].append(yPred[i][0])\n",
    "        XValidListNew = np.delete(XValidListNew, 0, axis=1)\n",
    "        XValidListNew = np.append(XValidListNew, yPred, axis=1)\n",
    "\n",
    "#     for yTrue, yPred in zip(yTest, yPredList):\n",
    "#         print(yTrue, yPred)\n",
    "\n",
    "    MAPE = meanAbsolutePercentageError_kDay(yValid, yPredListList)\n",
    "    print(seed, MAPE)\n",
    "    if MAPE < bestValidMAPE:\n",
    "        print('Updating best MAPE to {}...'.format(MAPE))\n",
    "        bestValidMAPE = MAPE\n",
    "        print('Updating best seed to {}...'.format(seed))\n",
    "        bestSeed = seed\n",
    "\n",
    "# define model\n",
    "print('Training model with best seed...')\n",
    "tf.random.set_seed(seed=bestSeed)\n",
    "model = Sequential()\n",
    "model.add(ConvLSTM2D(filters=64, kernel_size=(1,2), activation='relu', input_shape=(nSeq, 1, nSteps, nFeatures)))\n",
    "model.add(Flatten())\n",
    "model.add(Dense(1))\n",
    "opt = Adam(learning_rate=0.1)\n",
    "model.compile(optimizer=opt, loss='mse')\n",
    "# fit model\n",
    "model.fit(XTrainConv, yTrain[:,0], epochs=1000, verbose=0)\n",
    "\n",
    "yPredListList = []\n",
    "for day in range(nTest):\n",
    "    yPredListList.append([])\n",
    "XTestListNew = XTestList.copy()\n",
    "for day in range(k):\n",
    "    XTestConv = XTestListNew.reshape((XTestListNew.shape[0], nSeq, 1, nSteps, nFeatures))\n",
    "    yPred = model.predict(np.float32(XTestConv), verbose=0)\n",
    "    for i in range(len(yPred)):\n",
    "        yPredListList[i].append(yPred[i][0])\n",
    "    XTestListNew = np.delete(XTestListNew, 0, axis=1)\n",
    "    XTestListNew = np.append(XTestListNew, yPred, axis=1)\n",
    "    \n",
    "MAPE = meanAbsolutePercentageError_kDay(yTest, yPredListList)\n",
    "print('Test MAPE:', MAPE)\n",
    "MdSA = medianSymmetricAccuracy_kDay(yTest, yPredListList)\n",
    "print('Test MdSA:', MdSA)\n",
    "yPredConvLSTM = yPredListList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "nJZLQhMm4BXO",
    "outputId": "a3df1dca-56e5-4e02-ceee-f86d7cd70e10"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "\n",
    "# Format y tick labels\n",
    "def y_fmt(y, pos):\n",
    "    decades = [1e9, 1e6, 1e3, 1e0, 1e-3, 1e-6, 1e-9 ]\n",
    "    suffix  = [\"G\", \"M\", \"k\", \"\" , \"m\" , \"u\", \"n\"  ]\n",
    "    if y == 0:\n",
    "        return str(0)\n",
    "    for i, d in enumerate(decades):\n",
    "        if np.abs(y) >=d:\n",
    "            val = y/float(d)\n",
    "            signf = len(str(val).split(\".\")[1])\n",
    "            if signf == 0:\n",
    "                return '{val:d} {suffix}'.format(val=int(val), suffix=suffix[i])\n",
    "            else:\n",
    "                if signf == 1:\n",
    "                    if str(val).split(\".\")[1] == \"0\":\n",
    "                        return '{val:d} {suffix}'.format(val=int(round(val)), suffix=suffix[i]) \n",
    "                tx = \"{\"+\"val:.{signf}f\".format(signf = signf) +\"} {suffix}\"\n",
    "                return tx.format(val=val, suffix=suffix[i])\n",
    "\n",
    "                #return y\n",
    "    return y\n",
    "\n",
    "\n",
    "# ARIMA\n",
    "def getPredictions(X, nDaysInFuture=5, invertible=False, pdqTuple = (1, 2, 2)):\n",
    "    p, d, q = pdqTuple\n",
    "    predList = []\n",
    "    Xcopy = X.copy()\n",
    "    for i in range(nDaysInFuture):\n",
    "        if invertible:\n",
    "            model = SARIMAX(Xcopy, order=(p, d, q))\n",
    "\n",
    "            model_fit = model.fit(disp=False)\n",
    "\n",
    "            # make prediction\n",
    "            yhat = model_fit.predict(len(Xcopy), len(Xcopy), typ='levels')\n",
    "        else:\n",
    "            model = SARIMAX(Xcopy, order=(p, d, q), enforce_invertibility=False)\n",
    "\n",
    "            model_fit = model.fit(disp=False)\n",
    "\n",
    "            # make prediction\n",
    "            yhat = model_fit.predict(len(Xcopy), len(Xcopy), typ='levels')\n",
    "        Xcopy = np.append(Xcopy, yhat)\n",
    "        predList.append(np.around(yhat[0]))\n",
    "        \n",
    "    return predList\n",
    "\n",
    "yPredARIMA = getPredictions(nCasesGreaterDaysSinceList[-15-k:-k], nDaysInFuture=k)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "datesForPlottingList = datesGreaterDaysSinceList[-k:]\n",
    "groundTruthList = nCasesGreaterDaysSinceList[-k:]\n",
    "\n",
    "\n",
    "plt.ylabel('Number of cases', fontsize=20)\n",
    "plt.plot(datesForPlottingList, groundTruthList, '-o', linewidth=3, label='Ground Truth');\n",
    "plt.plot(datesForPlottingList, yPredVanillaLSTM[-1], '-o', linewidth=3, label='Vanilla LSTM');\n",
    "plt.plot(datesForPlottingList, yPredStackedLSTM[-1], '-o', linewidth=3, label='Stacked LSTM');\n",
    "plt.plot(datesForPlottingList, yPredBidirectionalLSTM[-1], '-o', linewidth=3, label='Bidirectional LSTM');\n",
    "plt.plot(datesForPlottingList, yPredCNNLSTM[-1], '-o', linewidth=3, label='CNN LSTM');\n",
    "plt.plot(datesForPlottingList, yPredConvLSTM[-1], '-o', linewidth=3, label='ConvLSTM');\n",
    "plt.plot(datesForPlottingList, yPredARIMA, '-o', linewidth=3, label='ARIMA');\n",
    "plt.xlabel('Date', fontsize=20);\n",
    "plt.legend(fontsize=14);\n",
    "plt.xticks(fontsize=16);\n",
    "plt.yticks(fontsize=16);\n",
    "ax = plt.gca()\n",
    "ax.yaxis.set_major_formatter(FuncFormatter(y_fmt))\n",
    "# plt.grid(axis='y')\n",
    "plt.savefig(os.path.join('plots', 'predictions_{}.png'.format(countryName)), dpi=400)\n",
    "plt.savefig(os.path.join('plots', 'predictions_{}.pdf'.format(countryName)), dpi=400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1KEpSgKs4BYE"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "LSTM_time-series-global_nDays6_k-day3.ipynb",
   "provenance": []
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}