convlstm.py

# -*- coding: utf-8 -*-
"""ConvLSTM_running.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1TT_ptKSRw7tn9F5_mVK0kF9JjeLq4u0z
"""

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import datetime
from datetime import date,timedelta
from google.colab import files

import xgboost as xgb
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import KFold

import time
import datetime

import numpy as np
import pandas as pd
import tensorflow as tf

import keras
from keras.models import Sequential
from keras.layers import *
from keras.optimizers import RMSprop
from keras.callbacks import CSVLogger, EarlyStopping

from google.colab import drive
drive.mount('/content/drive')

path= '/content/drive/Shareddrives/MSc - Shiveswarran/Processed data/Nine_months_data/bus_running_times_feature_added_all.csv'
data = pd.read_csv(path)

data = data.drop(data[data['run_time_in_seconds'] > 2000].index )
data = data.loc[(data['time_of_day']>= 6) & (data['time_of_day']<19)]
data = data[data['direction']==1]

"""PrepData

"""

data['DateTime'] = pd.to_datetime(data['date'] + ' ' + data['start_time'])

ref_freq = '15min'
ix = pd.DatetimeIndex(pd.to_datetime(data['DateTime'])).floor(ref_freq)

data["DateTimeRef"] = ix

data['DayOfWeek'] = pd.to_datetime(data['date']).dt.weekday
data["DowTimeRef"] = pd.to_datetime((data['DayOfWeek'].values ) * 24 * 60 * 60 + ix.hour * 60 * 60 + ix.minute * 60, unit = 's')

data

def datetime_index_generator(start_date,end_date,start_time,end_time,freq):
  date_range = pd.date_range(start_date, end_date, freq='D')
  date_list = [i.strftime('%Y-%m-%d') for i in date_range]
  datetime_series = pd.Series()

  for date in date_list:
    begin_timestamp = pd.to_datetime(date + ' ' + start_time )
    end_timestamp = pd.to_datetime(date + ' ' + end_time )
    daily_time_range = pd.Series(index = pd.date_range(begin_timestamp,end_timestamp,freq = freq,closed='left'))

    datetime_series = pd.concat([datetime_series,daily_time_range ]).sort_index()

  return datetime_series

"""Fit_scale"""

def fit_scale(data, ref_freq = '15min'):
  means = { }
  scales = { }
  low = { }
  upr = { }

  grouping = data[data['run_time_in_seconds'].notnull()].groupby('segment', sort = False)

  for segment, data_link in grouping:
        # Fit outlier bounds using MAD
        median = data_link.groupby('DowTimeRef')['run_time_in_seconds'].median()
        error = pd.concat([data_link['DowTimeRef'], np.abs(data_link['run_time_in_seconds'] - median[data_link['DowTimeRef']].values)], axis = 1)
        mad = 1.4826 * error.groupby('DowTimeRef')['run_time_in_seconds'].median()
        
        _low = median - 3 * mad
        _upr = median + 3 * mad
        mask = (_low[data_link['DowTimeRef']].values < data_link['run_time_in_seconds']) & (data_link['run_time_in_seconds'] < _upr[data_link['DowTimeRef']].values)
        data_link_no = data_link[mask]
        
        _mean = data_link_no.groupby('DowTimeRef')["run_time_in_seconds"].mean()
        means[segment] = _mean
        low[segment] = _low
        upr[segment] = _upr
        scales[segment] = data_link_no['run_time_in_seconds'].std()

  #ix = pd.date_range('1970-01-01', '1970-01-08', freq = ref_freq, closed = 'left')
  index_series = datetime_index_generator('1970-01-01','1970-01-07','06:00:00','19:00:00',freq='15min')
  ix = index_series.index
  means_df = pd.DataFrame(data = means, index = ix).interpolate()
  low_df = pd.DataFrame(data = low, index = ix).interpolate()
  upr_df = pd.DataFrame(data = upr, index = ix).interpolate()

  means_df = means_df.fillna(method='pad').fillna(method='bfill')
  low_df = low_df.fillna(method='pad').fillna(method='bfill')
  upr_df = upr_df.fillna(method='pad').fillna(method='bfill')

  return means_df, scales, low_df, upr_df

"""Time series index generator"""

s1 = datetime_index_generator('2021-10-01','2022-02-28','06:00:00','19:00:00',freq='15min')
s2 = datetime_index_generator('2022-07-01','2022-11-01','06:00:00','19:00:00',freq='15min')

s = pd.concat([s1,s2]).sort_index()

datetime_index = s.index

datetime_index

"""Outliers"""

def remove_outliers(data, low, upr): 
    _low = low.lookup(data['DowTimeRef'], data['segment'])
    _upr = upr.lookup(data['DowTimeRef'], data['segment'])
    mask = ((_low < data['run_time_in_seconds']) & (data['run_time_in_seconds'] < _upr))
    data = data.loc[mask].copy()
    return data, (~mask).sum()

"""Transform"""

def transform(data, means_df, scales,datetime_index, freq = '15min'):
  tss = { }
  ws = { }
  removed_mean = { }
  removed_scale = { }
  ks = []
  
  for k, v in data.groupby('segment', sort = False):
        # Link Data Time Indexed
        link_time_ix = pd.DatetimeIndex(pd.to_datetime(v['DateTime']))    
        link_time_ixd = v.set_index(link_time_ix)
        
        # Link Reference Data Index
        ix_ref = link_time_ixd['DowTimeRef']  

        link_travel_time_k = link_time_ixd['run_time_in_seconds'].resample(freq).mean()
        link_travel_time_k = link_travel_time_k[link_travel_time_k.index.isin(datetime_index)]
        
        removed_mean[k] = pd.Series(data = means_df.loc[ix_ref, k].values, index = link_time_ix).resample(freq).mean()
        removed_mean[k] = removed_mean[k][removed_mean[k].index.isin(datetime_index)]

        removed_scale[k] = pd.Series(data = np.repeat(scales[k], link_travel_time_k.shape[0]), index = link_travel_time_k.index)
        tss[k] = (link_travel_time_k - removed_mean[k].values) / removed_scale[k].values
        
        ws[k] = link_time_ixd['run_time_in_seconds'].resample(freq).count()
        ws[k] = ws[k][ws[k].index.isin(datetime_index)]

        ks.append(k)

  ts = pd.DataFrame(data = tss).fillna(method='pad').fillna(0) # Link Travel Time Time Series
  df_removed_mean = pd.DataFrame(data = removed_mean, index = ts.index).fillna(method='pad').fillna(method='bfill') # Removed Mean from Link Travel Time
  df_removed_scale = pd.DataFrame(data = removed_scale, index = ts.index).fillna(method='pad').fillna(method='bfill')
  w = pd.DataFrame(data = ws).fillna(0) # Link Travel Time Weights, e.g. number of measurements

  return (ts.index, ts.values, df_removed_mean.values, df_removed_scale.values, w.values, ks)

def roll(ix, ts, removed_mean, removed_scale, w, lags, preds):
    X = np.stack([np.roll(ts, i, axis = 0) for i in range(lags, 0, -1)], axis = 1)[lags:-preds,]
    Y = np.stack([np.roll(ts, -i, axis = 0) for i in range(0, preds, 1)], axis = 1)[lags:-preds,]
    Y_ix = ix[lags:-preds]
    Y_mean = np.stack([np.roll(removed_mean, -i, axis = 0) for i in range(0, preds, 1)], axis = 1)[lags:-preds,]
    Y_scale = np.stack([np.roll(removed_scale, -i, axis = 0) for i in range(0, preds, 1)], axis = 1)[lags:-preds,]
    w_y = np.stack([np.roll(w, -i, axis = 0) for i in range(0, preds, 1)], axis = 1)[lags:-preds,]

    return X, Y, Y_ix, Y_mean, Y_scale, w_y

"""Model Definition"""

def build_model(input_timesteps, output_timesteps, num_links):
    model = Sequential()
    model.add(BatchNormalization(name = 'batch_norm_0', input_shape = (input_timesteps, num_links, 1, 1)))
    model.add(ConvLSTM2D(name ='conv_lstm_1',
                         filters = 64, kernel_size = (10, 1),                       
                         padding = 'same', 
                         return_sequences = True))
    
    model.add(Dropout(0.2, name = 'dropout_1'))
    model.add(BatchNormalization(name = 'batch_norm_1'))

    model.add(ConvLSTM2D(name ='conv_lstm_2',
                         filters = 64, kernel_size = (5, 1), 
                         padding='same',
                         return_sequences = False))
    
    model.add(Dropout(0.1, name = 'dropout_2'))
    model.add(BatchNormalization(name = 'batch_norm_2'))
    
    model.add(Flatten())
    model.add(RepeatVector(output_timesteps))
    model.add(Reshape((output_timesteps, num_links, 1, 64)))
    
    model.add(ConvLSTM2D(name ='conv_lstm_3',
                         filters = 64, kernel_size = (10, 1), 
                         padding='same',
                         return_sequences = True))
    
    model.add(Dropout(0.1, name = 'dropout_3'))
    model.add(BatchNormalization(name = 'batch_norm_3'))
    
    model.add(ConvLSTM2D(name ='conv_lstm_4',
                         filters = 64, kernel_size = (5, 1), 
                         padding='same',
                         return_sequences = True))
    
    model.add(TimeDistributed(Dense(units=1, name = 'dense_1', activation = 'relu')))
    #model.add(TimeDistributed(Dense(units=1, name = 'dense_1', activation = 'relu')))
    #model.add(Dense(units=1, name = 'dense_2'))

    optimizer = RMSprop(lr=0.0001, rho=0.9, epsilon=1e-08, decay=0.9)
    model.compile(loss = "mse", optimizer = optimizer)
    return model


def info(msg):
    print(datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') + " " + msg)

"""Load Data"""



lags = 4 * 13
preds = 3

n = len(data)

"""Train and Test"""

data_train = data[data['week_no'] < 36]
data_test = data[36 <= data['week_no']]
n_train = len(data_train)
n_test = len(data_test)
info('- Train size : {:>8} ({:.2f}%) '.format(n_train, 100. * n_train / n))
info('- Test size  : {:>8} ({:.2f}%) '.format(n_test, 100. * n_test / n))

(means, scales, low, upr) = fit_scale(data_train)

data_train_no , no_removed =  remove_outliers(data_train,low,upr)

no_removed

data_test_no , no_removed =  remove_outliers(data_test,low,upr)
no_removed

#data_train_no = data_train
#data_test_no = data_test

ix_train, ts_train, rm_mean_train, rm_scale_train, w_train, lns_train = transform(data_train_no, means, scales,datetime_index)
ix_test, ts_test, rm_mean_test, rm_scale_test, w_test, lns_test = transform(data_test_no, means, scales,datetime_index)

#Create rolling window tensor
X_train, Y_train, Y_ix_train, Y_rm_mean_train, Y_scale_train, Y_w_train = roll(ix_train, ts_train, rm_mean_train, rm_scale_train, w_train, lags, preds)
X_test, Y_test, Y_ix_test, Y_rm_mean_test, Y_scale_test, Y_w_test = roll(ix_test, ts_test, rm_mean_test, rm_scale_test, w_test, lags, preds)

X_train = X_train[:,:,:,np.newaxis,np.newaxis]
Y_train = Y_train[:,:,:,np.newaxis,np.newaxis]
X_test = X_test[:,:,:,np.newaxis,np.newaxis]
Y_test = Y_test[:,:,:,np.newaxis,np.newaxis]

info('- X_train shape : {:>20}    X_test shape : {:>20}'.format(X_train.shape, X_test.shape))
info('- Y_train shape : {:>20}    Y_test shape : {:>20}'.format(Y_train.shape, Y_test.shape))

global_start_time = time.time()
early_stopping = EarlyStopping(monitor='val_loss', patience = 3)

model = build_model(lags, preds, len(lns_train))

# Train
history = model.fit(X_train, Y_train,
                    batch_size = 32, epochs = 1,
                    shuffle = False, validation_data = (X_test, Y_test),
                    verbose = 2, callbacks = [early_stopping])
model.save('models/ConvLSTM_3x15min_10x64-5x64-10x64-5x64_Comparison.h5') 

Y_true = Y_test.squeeze() * Y_scale_test + Y_rm_mean_test
Y_naive = Y_rm_mean_test
Y_pred = model.predict(X_test).squeeze() * Y_scale_test + Y_rm_mean_test

Y_true_total = np.sum(Y_true * Y_w_test, axis = 2).squeeze()
Y_naive_total = np.sum(Y_naive * Y_w_test, axis = 2).squeeze()
Y_pred_total = np.sum(Y_pred * Y_w_test, axis = 2).squeeze()
    
for t in range(preds):
    mask = Y_true_total[:,t] > 0
    Y_true_total_t = Y_true_total[mask, t] / 60
    Y_naive_total_t = Y_naive_total[mask, t] / 60
    Y_pred_total_t = Y_pred_total[mask, t] / 60  

    error_naive_total_t = (Y_naive_total_t - Y_true_total_t)
    error_lstm_total_t = (Y_pred_total_t - Y_true_total_t)

    mae_lstm = np.mean(np.abs(error_lstm_total_t))
    rmse_lstm = np.sqrt(np.mean((error_lstm_total_t)**2))
    mape_lstm = np.mean(np.abs(error_lstm_total_t) / Y_true_total_t) * 100

    info("- t + %d - ConvLSTM - MAE: %5.2f - RMSE: %5.2f - MAPE: %5.2f" % (t + 1, mae_lstm, rmse_lstm, mape_lstm))

model.summary()

!pip install neuralplot

from neuralplot import ModelPlot

# Commented out IPython magic to ensure Python compatibility.
# %matplotlib inline

modelplot = ModelPlot(model=model, grid=False, connection=True, linewidth=0.1)
modelplot.show()

plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'], linestyle = '--')

output = pd.DataFrame(data = Y_pred[:,0,:], index = Y_ix_test, columns = lns_test)

output['Datetime']=output.index

output


"""Out of Fold Training"""

















global_start_time = time.time()
early_stopping = EarlyStopping(monitor='val_loss', patience = 3)

bootstrap_size_pct = 0.80
test_window_pct = 0.05
max_iter = 4

lags = 4 * 8
preds = 3

hist = []
for i in range(max_iter):

    info("Current window: " + str(i))
    
    # Devide into test and train
    data_train = data[:int((bootstrap_size_pct + i * test_window_pct) * n)]
    data_test = data[int((bootstrap_size_pct + i * test_window_pct) * n):int((bootstrap_size_pct + (i + 1) * test_window_pct) * n)]
    n_train = len(data_train)
    n_test = len(data_test)
    info('- Train size : {:>8} ({:.2f}%) '.format(n_train, 100. * n_train / n))
    info('- Test size  : {:>8} ({:.2f}%) '.format(n_test, 100. * n_test / n))
    
    # Mean center and scale
    (means, scales, low, upr) = fit_scale(data_train)
    assert means.shape == (4 * 24 * 7, 32)
    assert len(scales) == 32
    assert low.shape == (4 * 24 * 7, 32)
    assert upr.shape == (4 * 24 * 7, 32)
    
    data_train_no = data_train
    data_test_no = data_test
    #data_train_no, n_outliers = remove_outliers(data_train, low, upr)
    #info('- Removed {0} outliers ({1:.2f}%) from train'.format(n_outliers, 100.0 * n_outliers / len(data_train)))
    #data_test_no, n_outliers = remove_outliers(data_test, low, upr)
    #info('- Removed {0} outliers ({1:.2f}%) from test'.format(n_outliers, 100.0 * n_outliers / len(data_test)))
    
    ix_train, ts_train, rm_mean_train, rm_scale_train, w_train, lns_train = transform(data_train_no, means, scales)
    ix_test, ts_test, rm_mean_test, rm_scale_test, w_test, lns_test = transform(data_test_no, means, scales)

X_train, Y_train, Y_ix_train, Y_rm_mean_train, Y_scale_train, Y_w_train = roll(ix_train, ts_train, rm_mean_train, rm_scale_train, w_train, lags, preds)
X_test, Y_test, Y_ix_test, Y_rm_mean_test, Y_scale_test, Y_w_test = roll(ix_test, ts_test, rm_mean_test, rm_scale_test, w_test, lags, preds)
X_train = X_train[:,:,:,np.newaxis,np.newaxis]
Y_train = Y_train[:,:,:,np.newaxis,np.newaxis]
X_test = X_test[:,:,:,np.newaxis,np.newaxis]
Y_test = Y_test[:,:,:,np.newaxis,np.newaxis]
    
info('- X_train shape : {:>20}    X_test shape : {:>20}'.format(X_train.shape, X_test.shape))
info('- Y_train shape : {:>20}    Y_test shape : {:>20}'.format(Y_train.shape, Y_test.shape))
    
model = build_model(lags, preds, len(lns_train))
    
    # Train
history = model.fit(X_train, Y_train,
                        batch_size = 128, epochs = 30,
                        shuffle = False, validation_data = (X_test, Y_test),
                        verbose = 2) #, callbacks = [csv_logger, early_stopping])
hist.append(history)
model.save('models/ConvLSTM_3x15min_10x64-5x64-10x64-5x64_' + str(i) + '.h5') 

Y_true = Y_test.squeeze() * Y_scale_test + Y_rm_mean_test
Y_naive = Y_rm_mean_test
Y_pred = model.predict(X_test).squeeze() * Y_scale_test + Y_rm_mean_test
        
Y_true_total = np.sum(Y_true * Y_w_test, axis = 2).squeeze()
Y_naive_total = np.sum(Y_naive * Y_w_test, axis = 2).squeeze()
Y_pred_total = np.sum(Y_pred * Y_w_test, axis = 2).squeeze()