# -*- coding: utf-8 -*-
"""
Ganjar Alfian
ganjar.alfian@ugm.ac.id
"""


import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn import metrics
from sklearn.tree import DecisionTreeClassifier
from sklearn.neural_network import MLPClassifier

from sklearn.model_selection import KFold, StratifiedKFold
import numpy as np

import xgboost as xgb
from xgboost import plot_importance, plot_tree
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.feature_selection import SelectFromModel

import seaborn as sns; sns.set()
from sklearn.metrics import mean_squared_log_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

from sklearn import linear_model
from numpy import array
from sklearn.neural_network import MLPRegressor

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from scipy.stats import kurtosis, skew
import tensorflow as tf
from sklearn.metrics import mean_absolute_error

import sys
if not sys.warnoptions:
    import warnings
    warnings.simplefilter("ignore")
#from sklearn.metrics import mean_absolute_percentage_error
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble import AdaBoostRegressor
#import graphviz
import math
import tensorflow as tf    
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor

import sys
if not sys.warnoptions:
    import warnings
    warnings.simplefilter("ignore")
    

# split a univariate sequence into samples
def split_sequences_direc(sequences, n_steps_in, n_steps_out):
	X, y = list(), list()
	for i in range(len(sequences)):
		# find the end of this pattern
		end_ix = i + n_steps_in
		out_end_ix = end_ix + n_steps_out
		# check if we are beyond the dataset
		if out_end_ix > len(sequences):
			break
		# gather input and output parts of the pattern
		seq_x, seq_y = sequences[i:end_ix, :-1], sequences[end_ix:out_end_ix, -1]
		X.append(seq_x)
		y.append(seq_y)
	return array(X), array(y)

def split_sequence(sequence, n_steps_in, n_steps_out):
	X, y = list(), list()
	for i in range(len(sequence)):
		# find the end of this pattern
		end_ix = i + n_steps_in
		out_end_ix = end_ix + n_steps_out
		# check if we are beyond the dataset
		if out_end_ix > len(sequence):
			break
		# gather input and output parts of the pattern
		seq_x, seq_y = sequence[i:end_ix, :-1], sequence[out_end_ix-1:out_end_ix, -1]
		X.append(seq_x)
		y.append(seq_y)
	return array(X), array(y)


def stats_features(input_data):
    inp = list()
    for i in range(len(input_data)):
        inp2=list()
        inp2=input_data[i]
        min=float(np.min(inp2))
        max=float(np.max(inp2))
        diff=(max-min)
        std=float(np.std(inp2))
        mean=float(np.mean(inp2))
        median=float(np.median(inp2))
        kurt=float(kurtosis(inp2))
        sk=float(skew(inp2))
        inp2=np.append(inp2,min)
        inp2=np.append(inp2,max)
        inp2=np.append(inp2,diff)
        inp2=np.append(inp2,std)
        inp2=np.append(inp2,mean)
        inp2=np.append(inp2,median)
        inp2=np.append(inp2,kurt)
        inp2=np.append(inp2,sk)
        #print(list(inp2))
        inp=np.append(inp,inp2)    
    inp=inp.reshape(len(input_data),-1)
    #print(inp)
    return inp

def mean_absolute_percentage_error(y_true, y_pred): 
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

def xgboost(X_train, X_test, y_train, y_test):
   
    model = xgb.XGBRegressor(random_state=0 )
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

def dec_tree(X_train, X_test, y_train, y_test):
   
    model = DecisionTreeRegressor(random_state=0)
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

def rf(X_train, X_test, y_train, y_test):
   
    model = RandomForestRegressor(random_state=0)
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

def adaboost(X_train, X_test, y_train, y_test):
   
    model = AdaBoostRegressor(random_state=0)
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae


def lin_reg(X_train, X_test, y_train, y_test):

    model = linear_model.LinearRegression()
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

def svr(X_train, X_test, y_train, y_test):

    model = SVR(C=200.0, epsilon=0.2)
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

def knn(X_train, X_test, y_train, y_test):

    model = KNeighborsRegressor()
    model = model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    
    
    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

def mlp(X_train, X_test, y_train, y_test):
    mlp_model = MLPRegressor( activation= 'relu', early_stopping=True,  max_iter=1000, random_state=0 )
    #random state = 2 for PH 3, count only
    #random state = 42 for PH 6, count only
    mlp_model.fit(X_train, y_train)
    y_pred = mlp_model.predict(X_test)

    rmse = math.sqrt(mean_squared_error(y_test,y_pred))
    mape=mean_absolute_percentage_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae= mean_absolute_error(y_test, y_pred)
        
    return y_pred, rmse, mape, r2, mae

   
if __name__ == '__main__':
    pd.options.display.max_columns=None
    df_ori = pd.read_csv('seoul hourly/2020_hourly.csv', sep = ',', encoding= 'unicode_escape')
    #df_ori = df_ori.dropna()
    #df_ori = df_ori[df_ori['Functioning Day'] != 'No']
    
    '''
    fig3 = plt.figure(figsize=(6,6))
    dt_plot = df_ori
    g = sns.heatmap(dt_plot.corr(),cmap="BrBG",annot=False, center=0)
    dt_plot.corr()
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    plt.tight_layout()
    plt.show()
    '''
    
    df_X = df_ori[['hour','day_week','month','count','count']]
    #df_X = df_ori[['count','count']]
    in_seq = df_X.astype(float).values
        
    n_steps_in, n_steps_out = 24, 24
    X, y = split_sequence(in_seq, n_steps_in, n_steps_out)
    
    n_input = X.shape[1] * X.shape[2]
    X_new = X.reshape((X.shape[0], n_input))
    
    X_train, X_test, y_train, y_test = train_test_split(X_new, y, test_size=0.25, random_state=0, shuffle = False)
    
    #X_train, X_test, y_train, y_test = train_test_split(X_train, y_train, test_size=0.25, random_state=0, shuffle = False)
    #from sklearn.preprocessing import MinMaxScaler
    #sc = MinMaxScaler()
    #X_train = sc.fit_transform(X_train)
    #X_test = sc.transform (X_test)

    
    print('Model evaluation')
    
    print('---------------')    
    print('SVR')
    y_pred_svr, rmse, mape, r2, mae= svr(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')
    
    print('---------------')    
    print('KNN')
    y_pred_knn, rmse, mape, r2, mae= knn(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')
        
    print('---------------')    
    print('DT')
    y_pred_dt, rmse, mape, r2, mae= dec_tree(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')

    print('---------------')    
    print('Adaboost')
    y_pred_ada, rmse, mape, r2, mae= adaboost(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')

    print('---------------')    
    print('RF')
    y_pred_rf, rmse, mape, r2, mae = rf(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')      
    
    print('---------------')
    print('Linear Reg')
    y_pred_lr, rmse, mape, r2, mae = lin_reg(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')
    
    print('---------------')
    print('XGB')
    y_pred_xgb, rmse, mape, r2, mae = xgboost(X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')  
    

    print('---------------')    
    print('MLP')
    y_pred_mlp, rmse, mape, r2, mae = mlp (X_train, X_test, y_train, y_test)
    print('rmse', round(rmse, 3))
    print('mae',round (mae, 3))
    print('r2',round(r2, 3))
    print('---------------')
    
    '''
    #-----------------------------
    fig1=plt.figure(figsize=(12,4))
    plt.title('PH = '+str(n_steps_out)+' hours', fontsize=15)
    plt.plot(y_test[-100:], label='Test data')
    plt.plot(y_pred_mlp[-100:], label='Proposed model')
    plt.legend(loc='upper right', prop={'size':15})
    plt.xlabel('Time t', fontsize=15)
    plt.ylabel('Number of rented bikes', fontsize=15)
    #plt.show()    
    fig1.savefig('MLP_rented_bikes_'+str(n_steps_out)+'.png', format='png', dpi=500)
    '''
   

    

    
