import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import json
import os
import itertools
from urllib.request import urlopen
from pandas.io.json import json_normalize
from datetime import datetime, timedelta
from IPython.core.interactiveshell import InteractiveShell

InteractiveShell.ast_node_interactivity = "all"
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (10, 5)
%matplotlib inline

response = urlopen("https://datacases.s3.us-east-2.amazonaws.com/datathon-2020/Ernst+and+Young/Dubai+Weather_20180101_20200316.txt")
rawdata = json.loads(response.read().decode('utf-8', 'replace'))

df = json_normalize(rawdata)
#weather = pd.DataFrame([y for x in df['weather'].values.tolist() for y in x])

for i in range(len(df.weather)):
    if len(df.weather[i]) > 1:
        print(f"{i}, {len(df.weather[i])}")

for i in range(len(df.weather)):
    df.weather[i] = df.weather[i][0]
condition = json_normalize(df.weather)

df = pd.concat([df, condition], axis=1)
data = df.copy()

data.head(3)

data = data.reset_index()

for i in range(len(data.dt)):
    data.loc[data.index==i, 'timestamp'] = datetime.fromtimestamp(data.dt[i])
    
data = data.drop(['index', 'dt', 'dt_iso', 'weather', 'id', 'icon', 'city_name', 
                  'lat', 'lon', 'timezone', 'description'], axis=1)

data['rain.1h'].fillna(0, inplace=True)
data['rain.3h'].fillna(0, inplace=True)

data.info()

data.rename(columns={'main.temp':'temp', 'main.temp_min':'temp_min', 'main.temp_max':'temp_max', 'main.feels_like': 'temp_feel', 
                     'main.pressure':'pressure', 'main.humidity':'humidity', 'clouds.all':'cloudiness',
                     'wind.speed':'wind_speed', 'wind.deg':'wind_deg', 'rain.1h':'rain_1h', 'rain.3h':'rain_3h', 
                     'main':'condition'}, inplace=True)

data.head()

data.describe()

f = plt.figure(figsize=(20,25));
gs = f.add_gridspec(4, 2)
ax1 = f.add_subplot(gs[0, 0])
sns.boxplot(data=data, x='condition', y='temp')
plt.title('Temperature Boxplots by Weather Condition')
ax2 = f.add_subplot(gs[0, 1])
sns.boxplot(data=data, x='condition', y='pressure')
plt.title('Atm. Pressure (hPa) by Weather Condition')
ax3 = f.add_subplot(gs[1, 0])
sns.boxplot(data=data, x='condition', y='humidity')
plt.title('Humidity (%) by Weather Condition')
ax4 = f.add_subplot(gs[1, 1])
sns.boxplot(data=data, x='condition', y='cloudiness')
plt.title('Cloudiness (%) by Weather Condition')
ax5 = f.add_subplot(gs[2, 0])
sns.boxplot(data=data, x='condition', y='wind_speed')
plt.title('Wind Speed (m/s) by Weather Condition')
ax6 = f.add_subplot(gs[2, 1])
sns.boxplot(data=data, x='condition', y='wind_deg')
plt.title('Wind Direction (degrees) by Weather Condition')
ax7 = f.add_subplot(gs[3, 0])
sns.boxplot(data=data, x='condition', y='rain_1h')
plt.title('Rain volume in mm (last 1 hour) by Weather Condition')
ax8 = f.add_subplot(gs[3, 1])
sns.boxplot(data=data, x='condition', y='rain_3h')
plt.title('Rain volume in mm (last 3 hours) by Weather Condition')
plt.show();

plt.rcParams['figure.figsize'] = (10, 5)
plt.title('Temperature throughout the period');
sns.lineplot(data=data, x='timestamp', y='temp_min', color='b');
sns.lineplot(data=data, x='timestamp', y='temp', color='yellow');
sns.lineplot(data=data, x='timestamp', y='temp_max', color='r');
sns.lineplot(data=data, x='timestamp', y='temp_feel', color='g');
plt.legend(['min','avg','max','feels']);

sns.lineplot(data=data, x='timestamp', y='pressure', color='c')
plt.title('Atmospheric pressure');

sns.lineplot(data=data, x='timestamp', y='humidity', color='orange')
plt.title('Humidity');

sns.lineplot(data=data, x='timestamp', y='cloudiness', color='c')
plt.title('Cloudiness');

f=plt.figure(figsize=(20,5))
gs = f.add_gridspec(1,2)
f.add_subplot(gs[0,0])
sns.lineplot(data=data, x='timestamp', y='wind_speed', color='darkcyan')
plt.title('Wind speed');
f.add_subplot(gs[0,1])
sns.lineplot(data=data, x='timestamp', y='wind_deg', color='darkolivegreen')
plt.title('Wind direction (High=N,W; Low=E,S)');

sns.lineplot(data=data, x='timestamp', y='rain_1h', color='y')
sns.lineplot(data=data, x='timestamp', y='rain_3h', color='r')
plt.title('Rain volume (mm) by date')
plt.legend(['last 1 hour', 'last 3 hours']);

sns.lineplot(data=data, x='condition', y='rain_1h', color='y');
sns.lineplot(data=data, x='condition', y='rain_3h', color='r');
plt.title('Rain volume by weather condition');
plt.legend(['last 1 hour','last 3 hours']);

corrmatrix = data.iloc[:,:9].corr()

plt.figure(figsize=(10,6))
ax = sns.heatmap(corrmatrix, vmin=-1, vmax=1, annot=True,);
bottom, top = ax.get_ylim()
ax.set_ylim(bottom + 0.5, top - 0.5)
plt.title('Correlation matrix')
plt.show();

f = plt.figure(figsize=(20,6))
gs = f.add_gridspec(1,2)
f.add_subplot(gs[0,0])
ax1 = sns.heatmap(data.loc[data.condition=='Clear','temp':'cloudiness'].corr(), vmin=-1, vmax=1, annot=True, cmap='YlGnBu');
bottom, top = ax1.get_ylim()
ax1.set_ylim(bottom + 0.5, top - 0.5)
plt.title('Correlation matrix - Clear condition')

f.add_subplot(gs[0,1])
ax2 = sns.heatmap(data.loc[data.condition=='Thunderstorm','temp':'cloudiness'].corr(), vmin=-1, vmax=1, annot=True, cmap='YlGnBu');
bottom, top = ax2.get_ylim()
ax2.set_ylim(bottom + 0.5, top - 0.5)
plt.title('Correlation matrix - Thunderstorms')
plt.show();

f = plt.figure(figsize=(20,6))
gs = f.add_gridspec(1,2)
f.add_subplot(gs[0,0])
sns.lineplot(data=data, x='condition', y='pressure')
plt.title('Atmospheric pressure (mm) by weather condition');
gs = f.add_gridspec(1,2)
f.add_subplot(gs[0,1])
sns.lineplot(data=data, x='condition', y='humidity');
plt.title('Humidity (%) by weather condition');
plt.show();

data.tail()

from selenium import webdriver
import time

url = 'https://www.dubaipulse.gov.ae/data/allresource/rta-bus/rta_bus_ridership-open?organisation=rta&service=rta-bus&count=1212&result_per_page=1300&as_sfid=AAAAAAVQxFA0BFeFROVV-_FUrwIfaEqwRWpoZA-y-UptqSEqxmERCKYLhWrwqWh3AfDCDdi1moQM5yS3Qjy2NzBMeMFf3DsQYwQOBarG4FRgrDCOeBE9L_Tq7J9m8CMoTDSCXIY%3D&as_fid=ff49229e06fa994326e53390b91e89d1dc5e2954'

driver = webdriver.Firefox(executable_path=r'D:\Library\env\geckodriver.exe')

driver.get(url)

results = driver.find_elements_by_xpath('//div[@class="file-action col-sm-12 col-xs-5"]')
print('Number of results', len(results))

# the code below is commented out in the submitted version of this file
#i=571
#for result in results[i:]:
    #result.click()
    #time.sleep(10)

csv = pd.read_csv(r"D:\Library\Datathon\2020\Bus\bus_ridership_2018-01-16_00-00-00.csv")

def timestamping(date, time):
    datesplit = date.split('-')
    timesplit = time.split(':')
    timestamp = datetime(int(datesplit[0]),int(datesplit[1]),int(datesplit[2]),int(timesplit[0]),int(timesplit[1]),int(timesplit[2]))
    return timestamp

start = datetime.now()
csv['timestamp'] = csv.apply(lambda x: timestamping(x['txn_date'],x['txn_time']), axis=1)
print(f"time: {datetime.now() - start}")

csv = csv.sort_values(by=['route_name','end_location','timestamp']).reset_index()

start = datetime.now()
csv['diff_prev'] = csv.timestamp.diff()
print(f"time: {datetime.now() - start}")

csv.head()

csv2 = csv.loc[(csv.diff_prev > np.timedelta64(30,'s')) | (csv.diff_prev < np.timedelta64(0,'s')),:]

len(csv)
len(csv2)

csv2['duration'] = (abs(csv2.timestamp.diff(periods=-1)) / np.timedelta64(1,'s'))

csv2.iloc[-1,-1] = (csv.iloc[-1,-2] - csv2.iloc[-1,-3]) / np.timedelta64(1,'s')

csv2.head()

csv2.duration.describe()

sns.distplot(csv2.duration)
plt.title('Ride duration histogram (28-01-2018)');

ride_duration = pd.DataFrame(columns=['date','count', 'mean', 'std', 'min', '25%', '50%', '75%', 'max'])

# this cell was rerun as blank at the end, before the submission, as the output was too long and irrelevant to the article
files = []
src = "D:/Library/Datathon/2020/Bus/"
i=1

with os.scandir(src) as files:
    for file in files:
        print(f"Processing file: {i}, the time is: {datetime.now()}")
        csv = pd.read_csv(src + file.name)
        
        # sort values, create timestamp, calculate time difference (in seconds) between rides of each bus line
        csv['timestamp'] = csv.apply(lambda x: timestamping(x['txn_date'],x['txn_time']), axis=1)
        csv = csv.sort_values(by=['route_name','end_location','timestamp']).reset_index()
        csv['diff_prev'] = csv.timestamp.diff()
        
        # get the first record in each ride, on each stop, for each bus line
        csv2 = csv.loc[(csv.diff_prev > np.timedelta64(30,'s')) | (csv.diff_prev < np.timedelta64(0,'s')),:].copy()
        
        # calculate ride durations
        csv2['duration'] = (abs(csv2.timestamp.diff(periods=-1)) / np.timedelta64(1,'s'))
        
        # get date and duration descriptive statistics from each provided file
        filedate = file.name.split('_')[2]
        datetime_file = datetime(int(filedate.split('-')[0]), int(filedate.split('-')[1]), int(filedate.split('-')[2]))
        a = []
        b = []
        b.append(csv2.duration.describe().values)
        b = list(itertools.chain.from_iterable(b))
        a.append([datetime_file,b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7]])
        
        # save the date and duration descriptive statistics from each file to ride_duration data frame
        ride_duration = pd.concat([ride_duration, pd.DataFrame(a,
                                                               columns=['date','count', 'mean', 'std', 'min', '25%', '50%', 
                                                                        '75%', 'max'])])
        i+=1

ride_duration_backup = ride_duration.copy()

ride_duration.head()

sns.lineplot(data=ride_duration, x='date', y='50%');
plt.ylabel('median');
plt.title('Median bus ride duration by date');
plt.show();

data['date'] = pd.to_datetime(data.timestamp).apply(lambda x: x.date())
ride_duration['date'] = pd.to_datetime(ride_duration.date).apply(lambda x: x.date())

combined_data = data.join(ride_duration[['date','50%','std']].rename(columns={'50%':'median'}).set_index('date'), on='date')

combined_data.info()

combined_data.loc[combined_data['median'].isna()].date.unique()

combined_data['month'] = combined_data.date.apply(lambda x: str(x)[:7])

monthly_agg = combined_data.groupby('month')['median','std'].median()

combined_data = combined_data.join(monthly_agg.rename(columns={'median':'monthly_median', 'std':'monthly_std'}), on='month')

combined_data.set_index('timestamp', inplace=True)

for i in combined_data.loc[combined_data['median'].isna()].index:
    combined_data.loc[combined_data.index == i, 'median'] = combined_data.loc[combined_data.index == i, 'monthly_median']
    combined_data.loc[combined_data.index == i, 'std'] = combined_data.loc[combined_data.index == i, 'monthly_std']

combined_data.info()

combined_data.tail(10)

combined_data.condition.fillna('Clear', inplace=True)

final_data = combined_data.drop(['date', 'description', 'month', 'monthly_median', 'monthly_std'], axis=1)

final_data.head()

X = final_data.iloc[:,:-2]
y = final_data.iloc[:,-2]

from sklearn.preprocessing import LabelEncoder, OneHotEncoder

cond_labelenc = LabelEncoder()
X.iloc[:,-1] = cond_labelenc.fit_transform(X.iloc[:,-1])

X.head(2)

from sklearn.compose import ColumnTransformer

ct = ColumnTransformer([('one_hot_encoder', OneHotEncoder(categories='auto'), [11])], remainder='passthrough')
X = ct.fit_transform(X)

X[0]

y = np.array(y)

sns.distplot(y);

y = np.log1p(y)

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

from sklearn.metrics import mean_absolute_error
from sklearn.model_selection import cross_val_score

from sklearn.linear_model import LinearRegression

regressor1 = LinearRegression()
regressor1.fit(X_train, y_train)

y_pred1 = np.expm1(regressor1.predict(X_test))

MAE1 = mean_absolute_error(np.expm1(y_test), y_pred1)
print(f"Mean absolute error: {MAE1}")

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV

parameters = {'n_estimators': [2, 3, 5, 10, 15, 20, 50, 100, 200, 500, 750, 1000],
              'max_leaf_nodes': [3, 5, 10, 20, 35, 50],
              'random_state': [0]}
grid_search = GridSearchCV(estimator = RandomForestRegressor(),
                           param_grid = parameters,
                           scoring = 'neg_mean_absolute_error',
                           cv = 5,
                           n_jobs = -1)
grid_search = grid_search.fit(X_train, y_train)
print(f"Best MAE: {grid_search.best_score_ * (-1)}")
print(f"Best parameters: {grid_search.best_params_}")

regressor2 = RandomForestRegressor(n_estimators=200, max_leaf_nodes=5, random_state=0)
regressor2.fit(X_train, y_train)

y_pred2 = np.expm1(regressor2.predict(X_test))

MAE2 = mean_absolute_error(np.expm1(y_test), y_pred2)
print(f"Mean absolute error: {MAE2}")

from sklearn.svm import SVR

from sklearn.preprocessing import StandardScaler

sc_X = StandardScaler()
sc_y = StandardScaler()
sc_X_train = sc_X.fit_transform(X_train)
sc_y_train = sc_y.fit_transform(y_train.reshape(-1,1))
sc_X_test = sc_X.fit_transform(X_test)

parameters = {'C': [0.5, 1, 5, 10, 20, 50, 100],
              'kernel': ['rbf']}
grid_search = GridSearchCV(estimator = SVR(),
                           param_grid = parameters,
                           scoring = 'neg_mean_absolute_error',
                           cv = 5,
                           n_jobs = -1)
grid_search = grid_search.fit(sc_X_train, sc_y_train)
print(f"Best MAE: {grid_search.best_score_ * (-1)}")
print(f"Best parameters: {grid_search.best_params_}")

regressor3 = SVR(kernel = 'rbf', C = 20)
regressor3.fit(sc_X_train, sc_y_train)

y_pred3 = regressor3.predict(sc_X_test)
y_pred3 = np.expm1(sc_y.inverse_transform(y_pred3))

MAE3 = mean_absolute_error(np.expm1(y_test), y_pred3)
print(f"Mean absolute error: {MAE3}")

from xgboost import XGBRegressor

parameters = {'base_score': [0.1, 0.5, 0.7, 1, 2],
              'learning_rate': [0.01, 0.03, 0.1, 0.3],
              'n_estimators': [20, 50, 100, 150, 200, 300, 400, 1000],
              'max_depth': [3, 5, 7, 9]}
grid_search = GridSearchCV(estimator = XGBRegressor(),
                           param_grid = parameters,
                           scoring = 'neg_mean_absolute_error',
                           cv = 5,
                           n_jobs = -1)
#grid_search = grid_search.fit(X_train, y_train)
#print(f"Best MAE: {grid_search.best_score_ * (-1)}")
#print(f"Best parameters: {grid_search.best_params_}")

regressor4 = XGBRegressor(learning_rate=0.02, n_estimators=300, max_depth=9)
regressor4.fit(X_train, y_train)
y_pred4 = np.expm1(regressor4.predict(X_test))

MAE4 = mean_absolute_error(np.expm1(y_test), y_pred4)
print(f"Mean absolute error: {MAE4}")

models_summary = {'Multiple Linear': MAE1, 'Random Forest': MAE2, 'SVR': MAE3, 'XGB': MAE4}

sns.barplot(x=list(models_summary.keys()), y=list(models_summary.values()), color='r')
plt.title("Mean absolute error by model (the lower, the better)");

plt.figure(figsize=(20,10))
sns.lineplot(data=y_test, color='r').set(yscale='log');
sns.lineplot(data=np.log1p(y_pred4), color='b').set(yscale='log');
sns.lineplot(data=np.log1p(y_pred3), color='y').set(yscale='log');
sns.lineplot(data=np.log1p(y_pred2), color='brown').set(yscale='log');
sns.lineplot(data=np.log1p(y_pred1), color='cyan').set(yscale='log');
plt.title('Models summary: predictions vs the test set')
plt.legend(['Test set','XGB prediction','SVM prediction','RF prediction','Linear prediciton']);