# LSTMs for Multi-Step Forecasting
#
#>>>(1) LSTM model with vector output for multi-step forecasting with univariate input data.

# Vanilla LSTM model that reads in a sequence of days of total hourly 
# power consumption and predicts a vector output of the next standard 
# day of hourly power consumption. So, output will be a set of 24 values,
# one for each hour of next day. Input is 24 values of previous day or
# 24*number-of-days. 
# A good starting point would be to use the prior day or seven days?


# univariate multi-step lstm
from math import sqrt
from numpy import split
from numpy import array
from pandas import read_csv
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Flatten
from keras.layers import LSTM

import numpy as np
 


# split a univariate dataset into train/test sets
def split_dataset(data):
    # split into standard days; each day with 24 values;
    # we have 8400 rows in the dataset, one for each hour
    # 8400 = 24 x 350 days
    # let's split it roughly into 70% and 30% for training and or testing:
    # 8400 = 5880 + 2520 = 245 + 105 days
    # let's split it roughly into 70% and 30% for training and test subsets;
    split_id = int( 0.7 * len(data) )
    # make split_id an integer number  of 24 values
    #split_id = int( split_id - (split_id % 24) )
    print('Splitting index: ', split_id)
    # split_id should be 24096 = 1004 x 24 representing about 70%
    # the remaining data represents test subset
    # split and leave the last day out, because that will only be predicted
    # tested and not used as input for a new prediction because
    # there nothing after it to do prediction for;
    train, test = data[0:split_id], data[split_id:(len(data)-1)]
    # restructure into windows of daily data
    train = array(split(train, len(train)/1))
    test = array(split(test, len(test)/1))
    # at this time train.shape = (245, 24, 7) with 7 columns, i.e., 7 variables
    # and test.shape = (104, 24, 7)
    return train, test
 
    
# evaluate one or more weekly forecasts against expected values
def evaluate_forecasts(actual, predicted):
    scores = list()
    # calculate an RMSE score for each day; 
    # actual.shape[1] should be 24 in our case;
    for i in range(actual.shape[1]):
        # calculate mse
        mse = mean_squared_error(actual[:, i], predicted[:, i])
        # calculate rmse
        rmse = sqrt(mse)
        # store
        scores.append(rmse)
    # calculate overall RMSE
    # actual.shape[0] should be 26 in our case - as the number of days in 
    # the test set;
    s = 0
    for row in range(actual.shape[0]):
        for col in range(actual.shape[1]):
            s += (actual[row, col] - predicted[row, col])**2
    score = sqrt(s / (actual.shape[0] * actual.shape[1]))
    return score, scores

# evaluate one or more forecasts against expected values
# as relative errors expressed as percentages
# 100 * (predicted - actual)/actual
# in other  words, report mean absolute percentage error (MAPE)
# which is commonly used to measure the predictive accuracy of models
def percentage_error(actual, predicted):
    res = np.empty(actual.shape)
    for j in range(actual.shape[0]):
        if actual[j] != 0:
            res[j] = (actual[j] - predicted[j]) / actual[j]
        else:
            res[j] = predicted[j] / np.mean(actual)
    return res

def mean_absolute_percentage_error(y_true, y_pred): 
    return np.mean(np.abs(percentage_error(np.asarray(y_true), np.asarray(y_pred)))) * 100

def evaluate_forecasts_MAPE(actual, predicted):
    scores = list()
    # actual.shape[1] should be 24 in our case;
    for i in range(actual.shape[1]):
        # calculate mape for this hour
        mape = mean_absolute_percentage_error(actual[:, i], predicted[:, i])
        print(mape)
        scores.append(mape)
    # calculate overall average MAPE over all 24 hours
    score = sum(scores)/len(scores)
    return score, scores

# summarize scores
def summarize_scores(name, score, scores):
    s_scores = ', '.join(['%.1f' % s for s in scores])
    print('%s: [%.3f] %s' % (name, score, s_scores))

    
# convert history into inputs and outputs
def to_supervised(train, n_input, n_out=1): # 24, 6
    # flatten data
    data = train.reshape((train.shape[0]*train.shape[1], train.shape[2]))
    X, y = list(), list()
    in_start = 0
    # step over the entire history one time step at a time
    for _ in range(len(data)):
        # define the end of the input sequence
        in_end = in_start + n_input
        out_end = in_end + n_out
        # ensure we have enough data for this instance
        if out_end <= len(data):
            # take the first 24 rows, only the 1-st column index 0 
            # (which is total Net_total) of the data that has shape 
            # (8400,6); that is "total calculated power" each hour,
            # for 24 hours, i.e., 1 days;
            x_input = data[in_start:in_end, 1] ### 0 is "Net total", 1, is "Lights", 2 is "HVAC"
            x_input = x_input.reshape((len(x_input), 1))
            X.append(x_input)
            y.append(data[in_end:out_end, 1]) ### 0 is "Net total", 1, is "Lights", 2 is "HVAC"
        # move along one time step
        # NOTE: when prediction was done for 7 days values, incrementing
        # one day at a time here was ok; it was sun-sat, mon-sun, tue-mon, 
        # wed-tue, etc.
        # this was done like that to increase amount of train data; it was
        # ok because the daily total consumption does not vary that much;
        # however, when doing 24 h for 24 values prediction, incrementing
        # hourly for the sake of increasing amount of train data makes for the
        # outputs to predict "average values" because training of model is done
        # with data 0-23h, 1-0h, 2-1h, 3-2h etc. instead of 
        # 0-23h, 0-23h, 0-23h, 0-23h, etc.
        # overlapping days I think diminishes the "daily pattern" of power 
        # and kind-of-uniformizes it;
        in_start += 1 # 24,6
    return array(X), array(y)
 
    
# train the model
def build_model(train, n_input):
    # (1) prepare data; next function uses hard coded numbers, based
    # on knowing that 1st column, is "Net_total"
    train_x, train_y = to_supervised(train, n_input)
    print('After flattening the train dataset with function to_supervised:')
    print('train_x shape:')
    print(train_x.shape)
    print('train_y shape:')
    print(train_y.shape)
    
    
    # (2) define parameters;
    # I selected these randomly for now;
    #verbose, epochs, batch_size = 1, 32, 16 # 1,70,16; 1,30,32 
    verbose, epochs, batch_size = 1, 5, 16
    n_timesteps = train_x.shape[1]
    n_features  = train_x.shape[2]
    n_outputs   = train_y.shape[1]
    print(n_timesteps)
    print(n_features)
    print(n_outputs)
    
    # (3) define model
    # Model #1: Simple Naive
    # model = Sequential()
    # model.add(LSTM(100, activation='relu', input_shape=(n_timesteps, n_features)))
    # model.add(Dense(50, activation='relu'))
    # model.add(Dense(n_outputs))
    # model.compile(loss='mse', optimizer='adam')
    # Model #2: More layers
    #model = Sequential()
    #model.add(LSTM(200, activation='relu', return_sequences=True, input_shape=(n_timesteps, n_features)))
    #model.add(LSTM(200, activation='relu', return_sequences=True))
    #model.add(LSTM(200, activation='relu'))
    #model.add(Dense(n_outputs))
    #model.compile(loss='mse', optimizer='adam') 
    model = Sequential()
    model.add(LSTM(56, activation='relu', return_sequences=True, input_shape=(n_timesteps, n_features)))
    model.add(LSTM(28, activation='relu', return_sequences=False))    
    model.add(Dense(14, activation='relu'))
    #model.add(Activation('linear')) # linear, relu, softmax, sigmoid
    model.add(Dense(n_outputs))
    model.compile(loss='mse', optimizer='adam') 

    
    # (4) fit network
    model.fit(train_x, train_y, epochs=epochs, batch_size=batch_size, verbose=verbose)
    return model
 
    
# make a forecast
def forecast(model, history, n_input):
    # flatten data
    data = array(history)
    data = data.reshape((data.shape[0]*data.shape[1], data.shape[2]))
    
    # retrieve last observations for input data
    # get the values from column 0, which is total Net_total power;
    input_x = data[-n_input:, 1] ### 0 is "Net total", 1, is "Lights", 2 is "HVAC"
    # reshape into [1, n_input, 1]
    input_x = input_x.reshape((1, len(input_x), 1))
    
    # forecast the next day, as 24 values
    yhat = model.predict(input_x, verbose=0)
    # we only want the vector forecast
    yhat = yhat[0]
    return yhat
 
  
# Walk-Forward Validation 
# evaluate a single model
def evaluate_model(train, test, n_input, my_scaler):
    # (1) fit model
    # this function does some more data manipulation; uses train
    # to generate train_x and train_y: by calling to_supervised(train, n_input);
    # those train_x and train_y take only the column with total load power
    # from the train data;
    model = build_model(train, n_input)
    
    # (2) history is a list of daily data; a "daily data point" is 24 
    # sets (for 24 h), with 6 variables in each set, where  the last one
    # in each set is the total load power for that hour;
    history = [x for x in train]
    
    # walk-forward validation over each week
    actuals = list()
    predictions = list()
    
    # create two lists to aggregate all predictions
    # we'll use them for ploting purposes only;
    plot_actuals = list()
    plot_predictions = list()
    
    
    for i in range(len(test)):
        # a) predict the day as 24 values, each values one of the 24 h
        yhat_sequence = forecast(model, history, n_input)
        
        
        # aggregate actuals and predictions for plotting
        plot_act = test[i, :]
        plot_act = plot_act[:, 1] ### 0 is "Net total", 1, is "Lights", 2 is "HVAC"
        plot_actuals.append(plot_act)
        plot_pred = yhat_sequence
        plot_predictions.append(plot_pred)
        
        
        if (i < 0):
            print ('i = %d' %(i))
            #print(plot_act)
            #print(plot_pred)
            hours = list(range(0, 1))  
            plt.plot( hours, plot_act, marker='', color='olive', 
                     linewidth=1, label='Actual')
            plt.plot( hours, plot_pred, marker='', color='red', 
                     linewidth=1, linestyle='dashed', label="Predicted")
            plt.title('LSTM Model: Input: Previous Day, Output: Next Day')
            plt.xlabel('24 Hours')
            plt.ylabel('Lights Energy [kWh]')
            #plt.ylabel('HVAC Energy [kWh]')
            #plt.ylabel('Total Energy [kWh]')
            plt.legend(loc="upper left")
            plt.show()
               
            
        # b) store the predictions
        
        #actuals.append(plot_act)
        #predictions.append(yhat_sequence)
        #print('plot_act normalized:')
        #print(plot_act)
        #print('plot_act inversed:')
        #print(my_scaler.inverse_transform(array(plot_act).reshape(-1,1)))
        #print('predicted normalized:')
        #print(yhat_sequence)
        #print('predicted inversed:')
        #print(my_scaler.inverse_transform(array(yhat_sequence).reshape(-1,1)))
        actuals.append( my_scaler.inverse_transform(array(plot_act).reshape(-1,1)) )
        predictions.append( my_scaler.inverse_transform(array(yhat_sequence).reshape(-1,1)) )
        
        # c) get real observation and add to history for predicting the next day
        history.append(test[i, :])
       
        
        
    # plot ALL actuals and ALL predictions for 105 (actually 104 because
    # for the last one we do not do it because we do not have info for
    # the next day, which would be day 106) days from test
    plot_actuals = array(plot_actuals)
    plot_predictions = array(plot_predictions)
    print('plot_actuals shape:')
    print(plot_actuals.shape)
    print('plot_predictions shape:')
    print(plot_predictions.shape)
    plot_actuals = plot_actuals.reshape(
        (plot_actuals.shape[0]*plot_actuals.shape[1], 1))
    plot_predictions = plot_predictions.reshape(
        (plot_predictions.shape[0]*plot_predictions.shape[1], 1))
    print('plot_actuals shape after reshape:')
    print(plot_actuals.shape)
    print('plot_predictions shape after reshape:')
    print(plot_predictions.shape)
    plt.plot( plot_actuals, marker='', color='olive', 
             linewidth=1, label='Actual')
    plt.plot( plot_predictions, marker='', color='red', 
             linewidth=1, linestyle='dashed', label="Predicted")
    plt.title('LSTM Model: Input: Previous Day, Output: Next Day')
    plt.xlabel('Hours of All Test Days')   
    plt.ylabel('Lights Energy [kWh]')
    #plt.ylabel('HVAC Energy [kWh]')
    #plt.ylabel('Total Energy [kWh]')
    #plt.legend(loc="upper left")
    plt.show()
    

    # (3) evaluate predictions for each day
    #print('actuals as list:')
    #print(actuals)
    actuals = array(actuals)
    predictions = array(predictions)
    
    print('Just before evaluate_forecasts_MAPE():')
    print('actuals shape:')
    print(actuals.shape)
    print('predictions shape:')
    print(predictions.shape)
    
    ###score, scores = evaluate_forecasts(test[:, :, 0], predictions)
    #score, scores = evaluate_forecasts(actuals, predictions)
    score, scores = evaluate_forecasts_MAPE(actuals, predictions)
    return score, scores





# load the new file with hourly data
dataset = read_csv('HondaHome2018_NET_POWER_Hourly.csv',
                   header=0, 
                   infer_datetime_format=True, 
                   parse_dates=['datetime'], 
                   index_col=['datetime'])



# normalize data 
from pandas import Series
from sklearn.preprocessing import MinMaxScaler
#series = Series(dataset.Lighting[:])
#series = Series(array(dataset.values[:,0]))
series = Series(dataset[dataset.columns[1]]) ### 0 is "Net total", 1, is "Lights", 2 is "HVAC"
print(series[0])
values = series.values
#print(series)
values = values.reshape((len(values), 1))
# train the normalization
min_max_scaler = MinMaxScaler(feature_range=(0, 1))
min_max_scaler = min_max_scaler.fit(values)
print('Scaler Min: %f, Max: %f' % (min_max_scaler.data_min_, min_max_scaler.data_max_))
# normalize the dataset 
normalized = min_max_scaler.transform(values)
# copy normalized column back into dataset
normalized = normalized.reshape((len(normalized), ))
dataset[dataset.columns[1]] = normalized ### 0 is "Net total", 1, is "Lights", 2 is "HVAC"



# split into train and test
train, test = split_dataset(dataset.values)

# evaluate model and get scores
n_input = 3 # 168 hours = 7 days
score, scores = evaluate_model(train, test, n_input, min_max_scaler)

# summarize scores
summarize_scores('lstm', score, scores)

# plot scores
hours = list(range(0, 1)) 
plt.plot(hours, scores, marker='o', label='lstm')
#plt.title('LSTM Model Performance Evaluation')
plt.xlabel('24 Hours of Day')
#plt.ylabel('Root Mean Squared Error, RMSE [kW]')
plt.ylabel('Mean Absolute Percentage Error, MAPE [%]')
#plt.legend(loc="upper left")
plt.show()
