# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/DonnaVakalis/forecasting-MRT-gap/blob/main/model/model_BiLSTM_MSLE.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="PotS_WG4C0bC"
# # PREAMBLE (ABOUT THIS NOTEBOOK)

# + [markdown] id="_gBDZLPtC7VV"
# This notebook applies data prediction with Keras LSTM model in Python, on time series (hourly) weather data to predict daily thermal comfort exceedances indoors (e.g., when it gets too hot, and by how much). While this notebook contains only the first steps, ultimately it can help determine how much cooling is needed (in terms of mean radiant temperature or MRT), and ultimately help decide on which cooling technologies are feasible and most carbon efficient.

# + [markdown] id="7TblImSlCdcM"
# # SETUP AND IMPORTS

# + colab={"base_uri": "https://localhost:8080/"} id="h4uWABiYwGMt" outputId="ba57bef9-b2d1-4222-9562-8a53dfa9945d"
# first install requests, for importing libraries from github
# !pip install requests

# + [markdown] id="0fHyFFZ2svZc"
# Import pre-packaged libraries

# + id="GqrUpMCWCcCi"
import requests
import os
import sys
from google.colab import drive
import pickle
import glob
import numpy as np
import pandas as pd

from tensorflow.keras import layers
from tensorflow import keras 
import tensorflow as tf\

import torch # for custom loss function

from keras.models import Sequential
from keras.layers import LSTM, Dense, Bidirectional 

from numpy import array

import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error



# + [markdown] id="R8N14VdVtYgr"
# Mount drive and set paths

# + colab={"base_uri": "https://localhost:8080/"} id="ea4B0rT7tcvn" outputId="01745b48-9241-4180-a32a-c5b33d987470"
# Point to the folder for this project (forecasting-MRT-gap/data)
drive.mount('/content/gdrive')
DATA_PATH = "/content/gdrive/My Drive/Colab Notebooks/forecasting-MRT-gap/data"
FIGURES_PATH = "/content/gdrive/My Drive/Colab Notebooks/forecasting-MRT-gap/figures"
#LOSS_PATH = "/content/gdrive/My Drive/Colab Notebooks/forecasting-MRT-gap/loss"

# + [markdown] id="AAFXnOugCxQy"
# # LOAD + FILTER THE DATA
#

# + [markdown] id="VoRTvWQ-JxmY"
# The training and testing datasets were pre-processed in a separate notebook saved as separate pickles for 2015,2016 and 2017 (and separated by training and test and output variables based on a prior project) in this public repository:
#
# https://github.com/eta-lab/forecasting-comfort
#
# Details about methods from this project can be found here: <NAME>., <NAME>., <NAME>, (2021). Forecasting the impact of climate change on thermal comfort using a weighted ensemble of supervised learning models. Building & Environment.

# + [markdown] id="CJS4khlBBodr"
# ## GET OUTDOOR AIR TEMPERATURES FOR SUMMER 2015, 2016 and 2017

# + colab={"base_uri": "https://localhost:8080/", "height": 282} id="doivvZSBJv8T" outputId="71e77771-b54d-4cd8-a348-ce5585b958e0"
# 2015

dat_train = pickle.load(open(DATA_PATH +'/trainX_2015.pkl','rb'))
dat_test = pickle.load(open(DATA_PATH +'/testX_2015.pkl','rb'))

# check if train and test set are the same, and visualize temp values
df_diff = pd.concat([dat_train,dat_test]).drop_duplicates(keep=False)
print(df_diff.shape, dat_train.shape, dat_test.shape)
df_train = dat_train[['temp']]
df_test = dat_test[['temp']]
plt.plot(df_train.index, df_train.temp, df_test .index, df_test .temp)
plt.show() 

# + colab={"base_uri": "https://localhost:8080/", "height": 282} id="4AktS9jKsSzD" outputId="8ee2d401-f323-436c-c7c9-3fc4c5328750"
# 2016

dat_train = pickle.load(open(DATA_PATH +'/trainX_2016.pkl','rb'))
dat_test = pickle.load(open(DATA_PATH +'/testX_2016.pkl','rb'))

# check if 2015 train and test set are the same
df_diff = pd.concat([dat_train,dat_test]).drop_duplicates(keep=False)
print(df_diff.shape, dat_train.shape, dat_test.shape)
df_train = dat_train[['temp']]
df_test = dat_test[['temp']]
plt.plot(df_train.index, df_train.temp, df_test .index, df_test .temp)
plt.show() 

# + colab={"base_uri": "https://localhost:8080/", "height": 282} id="6oSRZIdgpvkw" outputId="6267a552-3e23-4f32-a4ef-2681f56eebc7"
# 2017

dat_train = pickle.load(open(DATA_PATH +'/trainX_2017.pkl','rb'))
dat_test = pickle.load(open(DATA_PATH +'/testX_2017.pkl','rb'))

# check if train and test set are the same, and visualize temp values
df_diff = pd.concat([dat_train,dat_test]).drop_duplicates(keep=False)
print(df_diff.shape, dat_train.shape, dat_test.shape)
df_train = dat_train[['temp']]
df_test = dat_test[['temp']]
plt.plot(df_train.index, df_train.temp, df_test .index, df_test .temp)
plt.show() 

# + [markdown] id="9Gi4O9iC5fmn"
# ## PRE-PROCESS THE X (INPUT VARIABLE) DATA

# + [markdown] id="Xto7qFE1Xy9X"
# For each year 2015,2016,2017, examine temperature data (outdoor) and prepare to remerge training/test sets (for dynamically selecting training/testing sets from whole)

# + [markdown] id="Jg-0q2wswNkq"
# ### Restrict to time frames we want 
#
# We want temperature data from all three summers, minus the smoky periods (smoke-affected days = when recorded incidences of regional forest fires)
#
# These are respectively the first two weeks of July 2015 and the first two weeks of August 2017

# + id="Ijkgnpjn61zq"
# 2015

# get training set from another year (so complete for 2015)
df_1 = pickle.load(open(DATA_PATH +'/trainX_2016.pkl','rb'))

# keep only summer dates, e.g., June, July and August:
remove_before = pd.to_datetime('2015-06-01')
remove_after = pd.to_datetime('2015-08-31')
df_2 = df_1.query('index >= @remove_before and index <= @remove_after')

# drop smoke weeks (first two weeks July 2015) 
start_smoke = pd.to_datetime('2015-07-01')
end_smoke = pd.to_datetime('2015-07-14')
df_3 = df_2.query('index < @start_smoke or index > @end_smoke')

# save 2015 outdoor temps (x variable)
df_2015_X = df_3[['temp']]

# + id="z5v_-qcrYCqn"
# 2016 

# get training set from another year (so complete for 2016)
df_1 = pickle.load(open(DATA_PATH +'/trainX_2017.pkl','rb'))

# keep only summer dates, e.g., June, July and August:
remove_before = pd.to_datetime('2016-06-01')
remove_after = pd.to_datetime('2016-08-31')
df_2 = df_1.query('index >= @remove_before and index <= @remove_after')

# save 2016 outdoor temps (x variable)
df_2016_X = df_2[['temp']]

# + id="bAfyy04v_R7X"
# 2017

# get training set from another year (so complete for 2017)
df_1 = pickle.load(open(DATA_PATH +'/trainX_2015.pkl','rb'))

# keep only summer dates, e.g., June, July and August:
remove_before = pd.to_datetime('2017-06-01')
remove_after = pd.to_datetime('2017-08-31')
df_2 = df_1.query('index >= @remove_before and index <= @remove_after')

# drop smoke weeks (first two weeks July 2015) 
start_smoke = pd.to_datetime('2017-08-01')
end_smoke = pd.to_datetime('2017-07-14')
df_3 = df_2.query('index < @start_smoke or index > @end_smoke')

# save 2017 outdoor temps (x variable)
df_2017_X = df_3[['temp']]

# + colab={"base_uri": "https://localhost:8080/", "height": 264} id="m0E6kQnz_pN7" outputId="90464eea-0caf-43b8-d3d4-b8a7b74e957b"
# plot all outdoor temperatures in cleaned dataset
plt.plot(df_2015_X.index, df_2015_X.temp, df_2016_X.index, df_2016_X.temp, df_2017_X.index, df_2017_X.temp) 
plt.show() 

# save a single dataset with all days
ls_X = [df_2015_X, df_2016_X, df_2017_X]  
df_all_X = pd.concat(ls_X)

# + [markdown] id="Oaovum2mB39U"
# ## GET INDOOR AIR TEMPERATURES (OUTPUT VARIABLE) FOR SAME TIME SPANS

# + id="0wMf2PSmB9by"
#Get raw data of indoor air temps (what we want to predict) 

df1 = pickle.load(open(DATA_PATH +'/trainY_2015.pkl','rb'))
df2 = pickle.load(open(DATA_PATH +'/trainY_2016.pkl','rb'))
df3 = pickle.load(open(DATA_PATH +'/trainY_2017.pkl','rb'))


# Put in all in one place 
ls_Y = [df1, df2, df3]  
df_Y = pd.concat(ls_Y)

# Drop duplicated date-times

df_Y = df_Y[~df_Y.index.duplicated(keep='first')]

# Drop date-times not in X set
df_all_Y = df_Y[df_Y.index.isin(df_all_X.index)]

# + [markdown] id="0ICisOYGKIzA"
# ## Plot X and Y variables

# + colab={"base_uri": "https://localhost:8080/", "height": 264} id="M53g3iJhJczi" outputId="5d1aa1fb-4fdd-4340-8b6f-fe83965137fb"
# plot Outdoor (X values) and Indoor (to be predicted)
plt.plot(df_all_X.index, df_all_X.temp, df_all_Y.index, df_all_Y) 
plt.show() # Note LSTM assumes the book ends of a gap are adjacent, so it's good to see that there are no sharp changes over the gaps

# + [markdown] id="nSl76bU4f0l1"
# ## Add feature for time of day and transform variables

# + [markdown] id="GcgO__MkZWFO"
# This step adds a variable to help the model understand the cyclicity of the 24 hours cycle
# It borrows from the tutorial provided here: https://towardsdatascience.com/single-and-multi-step-temperature-time-series-forecasting-for-vilnius-using-lstm-deep-learning-b9719a0009de
#
# We want the machine to know that hours 23 and 0 are closer to each other than hours 0 and 4. We know that the period of the cycle is 24 hours. We can use the cos(x) and sin(x) functions. The x in the functions is the hour of the day.

# + id="seFBtHHEY8c3"
# Extract the hour of day
df_all_X["hour"] = df_all_X.index.hour
# Creating the cyclical daily feature 
df_all_X["day_cos"] = [np.cos(x * (2 * np.pi / 24)) for x in df_all_X["hour"]]
df_all_X["day_sin"] = [np.sin(x * (2 * np.pi / 24)) for x in df_all_X["hour"]]

# + [markdown] id="vQia6BLyVI5n"
# # BiLSTM MODEL with Logarithmic Mean Squared Error
#
# The idea of Bidirectional LSTMs (BiSTM) is to use input information from both the past and future of a specific time step.
#
# Some of the steps for making and fitting this model follow the tutorial found here: 
# https://towardsdatascience.com/predictive-analysis-rnn-lstm-and-gru-to-predict-water-consumption-
#

# + [markdown] id="DX64sYj2gCiH"
# ### Pre-processing to setup model: Split training/test sets, scale the data, shape the 3D matrix for feeding into the model
#
# Note: for this type of time-step depending analysis, make sure training/test split shuffle is set to False (i.e., keep continguous chunks for training and testing)

# + colab={"base_uri": "https://localhost:8080/"} id="eRcvJAhyVNKH" outputId="888e38f7-6cfa-473c-b0ba-a6c17048c748"
x = df_all_X.copy().to_numpy()
y = df_all_Y.copy().to_numpy()

# need to reshape Y which contains 1 features (currently it's 0) 
y = y.reshape(-1, 1)
 
# split the data into the train and test parts.
# NB split should not be a random selection... (shuffle=False)
X_train, X_test, Y_train, Y_test = train_test_split(x, y, shuffle = False, test_size= 0.15)
 
train_dim = (X_train.shape, Y_train.shape)
test_dim = (X_test.shape, Y_test.shape)

# check
print("train_dim Xs and Y:", train_dim)
print("test_dim Xs and Y", test_dim)

#Scale and center the data:
#1) variables are better to be small values, probably in the range of 0-1.
#2) want the scale of the output variable is in the range 0-1 to match the scale of the activation function (tanh) on the output layer 

# Fit the scaler with training data (used MinMaxScaler from sklearn)
scaler_x = MinMaxScaler(feature_range = (0,1)) 
scaler_y = MinMaxScaler(feature_range = (0,1))

# Fit the scaler using available training data
input_scaler = scaler_x.fit(X_train)
output_scaler = scaler_y.fit(Y_train)

# Apply the scaler to training data
train_y_norm = output_scaler.transform(Y_train)
train_x_norm = input_scaler.transform(X_train)

# Apply the scaler to test data
#NB: we should scale the unseen tes data with the scaler fitted on the training data.
test_y_norm = output_scaler.transform(Y_test)
test_x_norm = input_scaler.transform(X_test)


# + [markdown] id="oqZqxqfdoScR"
# Create a 3D Input Dataset
#
# This involves setting the 'time steps', which in this case will be the previous 48 hours

# + colab={"base_uri": "https://localhost:8080/"} id="c9eL9ZbkoUeu" outputId="c966baa4-abf3-4989-c08b-c76e99af9e70"
#Default is one day back and one day forward
def create_dataset (X, y, time_steps = 24): 
    Xs, ys = [], []
    
    for i in range(len(X)-time_steps):
        v = X[i:i+time_steps, :]
        Xs.append(v)
        ys.append(y[i+time_steps])
        
    return np.array(Xs), np.array(ys)


TIME_STEPS = 48 # previous and next two days' weather

X_test, Y_test = create_dataset(test_x_norm, test_y_norm, TIME_STEPS)
X_train, Y_train = create_dataset(train_x_norm, train_y_norm, TIME_STEPS)
print('X_train.shape: ', X_train.shape)
print('Y_train.shape: ', Y_train.shape)
print('X_test.shape: ', X_test.shape) 
print('Y_test.shape: ', Y_test.shape)


# + [markdown] id="1PTx0Gw-ryR6"
# ## Create the BiLSTM model with special loss functiom
#
# The 'Bidirectional' wrapper will create two copies of the hidden layer, one fit in the input sequences as-is and one on a reversed copy of the input sequence. By default, the output values from these LSTMs will be concatenated.
#
# That means that instead of the TimeDistributed layer receiving 48 timesteps of X HIDDEN_LAYER units, it will now receive 48 timesteps of 2*X HIDDEN_LAYER units (X units + X units) outputs.
#
#
# For the loss function, given the sharpness of our data, we want to relax the punishing effect of large differences in large predicted values. Therefore, instead of the default MSE, we will use Mean Squared Logarithmic Error loss, or MSLE: this calculates the natural logarithm of each of the predicted values, then calculate the mean squared error.  
#

# + id="35TxUo7en1Bn"
# Create the BiLSTM model

def create_model_bilstm(units):
    model = Sequential()
    # First layer of BiLSTM
    model.add(Bidirectional(LSTM(units = units, return_sequences=True), input_shape=(X_train.shape[1], X_train.shape[2])))
    # Second layer of BiLSTM
    model.add(Bidirectional(LSTM(units = units)))
    model.add(Dense(1))
    #Compile model
    model.compile(loss='mean_squared_logarithmic_error', optimizer= 'adam', metrics=['mse']) #change of loss  
    return model

# Choose the number of units (neurons) in hidden layers
HIDDEN_LAYERS = 64
model_bilstm = create_model_bilstm(HIDDEN_LAYERS)  


# + [markdown] id="jIGCfV8otYOT"
# ## Fit the model
#

# + colab={"base_uri": "https://localhost:8080/", "height": 743} id="LAdxTA_3dAxZ" outputId="d1340f2a-097f-4130-a39f-3527c38a448a"
# 10 epoch and batch_size = 64
def fit_model(model):
    early_stop = keras.callbacks.EarlyStopping(monitor = 'val_loss',
                                               patience = 10)
    # shuffle = False because the order of the data matters
    # use 20% of train data as validation data 
    history = model.fit(X_train, Y_train, epochs = 10, validation_split = 0.2,
                    batch_size = 64, shuffle = False, callbacks = [early_stop])
    return history

history= fit_model(model_bilstm)

# Visualize the training loss versus validation loss
plt.figure(figsize = (10, 6))
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])

plt.ylabel('Loss')
plt.xlabel('epoch')
plt.legend(['Train loss', 'Validation loss'], loc='upper right')
plt.savefig(FIGURES_PATH + '/BiLSTM_loss_plot with training and test sets continguous' +'.jpg', format='jpg', dpi=1000)
plt.show()

# + [markdown] id="nYuVPF94yUHC"
# Last step! # Very important to transform the target variable back to original data space :)
#

# + [markdown] id="vhUeRu4Uy4QA"
# ## Check Predictions of BiLSTM Model

# + colab={"base_uri": "https://localhost:8080/", "height": 404} id="hnZyEh7Oy8lt" outputId="0bf46427-48e6-4787-ba4a-4a54c82f8cab"
prediction = model_bilstm.predict(X_test)
prediction = scaler_y.inverse_transform(prediction)
true_future = scaler_y.inverse_transform(np.array(Y_test))
    
plt.figure(figsize=(10, 6))
range_future = len(prediction)

plt.plot(np.arange(range_future), true_future, label='True Future')
plt.plot(np.arange(range_future), np.array(prediction),label='Prediction')

plt.title('True future vs prediction for BiLSTM model')
plt.legend(loc='upper left')
plt.xlabel('Date')
plt.ylabel('Hourly Indoor Air Temperature (degrees C)')

plt.savefig(FIGURES_PATH + '/Model Check on Test Set Keeping training test batches continguous' +'.jpg', format='jpg', dpi=1000)
plt.show()

# + colab={"base_uri": "https://localhost:8080/"} id="uImyY81S2Jfe" outputId="cb447f41-47db-4a1d-85a6-10409498faae"
# Calculate MAE and RMSE

errors = prediction - Y_test
mse = np.square(errors).mean()
rmse = np.sqrt(mse)
mae = np.abs(errors).mean()

print('Mean Absolute Error: {:.4f}'.format(mae))
print('Root Mean Square Error: {:.4f}'.format(rmse))

# + [markdown] id="qwrGhlEWvuYl"
# Conclusion: this approach as-is performs badly!

# + [markdown] id="1ZR4boOnYpPl"
# # IDEAS for IMPROVEMENT: 
#
# BiLSTM model for sharpness of data: DILATE loss function
#
# In the literature this type of input data might be called "sharp", and the default scaling and loss functions lead to "too smooth" predictions...
#
# There are some special loss functions such as  DIstortion Loss including shApe and TimE (DILATE [link text](https://papers.nips.cc/paper/2019/file/466accbac9a66b805ba50e42ad715740-Paper.pdf)) or soft-dynamic time warping (soft-DTW) [link text](http://proceedings.mlr.press/v70/cuturi17a/cuturi17a.pdf)  which are specially designed for time series like this. 

# + id="c3N8m6L2Ymes"


# + [markdown] id="6DPPLKmYgOba"
# # ARCHIVED

# + id="215cBSe61Pk3"

# soft_dtw
url = "https://raw.githubusercontent.com/vincent-leguen/DILATE/master/loss/soft_dtw.py"
r = requests.get(url) 
with open('soft_dtw.py', 'w') as f:
    f.write(r.text)
import soft_dtw

# dilate_loss
url = "https://raw.githubusercontent.com/vincent-leguen/DILATE/master/loss/dilate_loss.py"
r = requests.get(url) 
with open('dilate_loss.py', 'w') as f:
    f.write(r.text)

# + id="aG2QG1IbeCSj"
# # !cp drive/MyDrive/mylib.py 
# # !cp "$inputFile" /destination/
# import mylib

InputFile = LOSS_PATH + '/dilate_loss' + '.py'
# !cp "InputFile"

InputFile2 = LOSS_PATH + '/soft_dtw' + '.py'
# !cp "InputFile2"


import dilate_loss.py
import soft_dtw.py

# + id="6VU9vPId56Sb"
# !git clone https://github.com/vincent-leguen/DILATE.git

# #%load dilate_loss.py
# %load path_soft_dtw.py
# #%load soft_dtw.py

# + id="kWcj15cmx5yc"
# Saved the files to drive (from repo here: https://github.com/vincent-leguen/DILATE.git)
# Because the files internally refer back to a path structure that is not preserved when uploading to google drive

#py_file_location = "/content/drive/My Drive"
#sys.path.append(os.path.abspath(py_file_location))

 
#sys.path.append(os.path.abspath(LOSS_PATH))
import LOSS_PATH + dilate_loss

StarcoderdataJupyterScriptsDedupFiltered