SF Bay Area Bike Data

# Notebook setup
import pandas as pd
import numpy as np
import glob
from pandas.tseries.holiday import USFederalHolidayCalendar
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import TimeSeriesSplit
from sklearn.model_selection import train_test_split


from sklearn.ensemble import RandomForestRegressor
from sklearn.multioutput import MultiOutputRegressor
from sklearn.ensemble import GradientBoostingRegressor

from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import GridSearchCV
from sklearn.externals import joblib

from rfpimp import *
from pprint import pprint

import folium

sns.set()

%matplotlib inline
# Set the option to display the max number of columns and rows 
pd.set_option("display.max_columns", 200)
pd.set_option("display.max_rows", 20000)

# Read the `station_data` 
station_df = pd.read_csv("station_data.csv")
station_df.head(2)

# Summary of station_data
station_df.info()

# Read the 'trip_data': trip_df
trip_df = pd.read_csv("trip_data.csv", 
                    parse_dates=['Start Date', 'End Date'], 
                    infer_datetime_format=True)

trip_df.head(2)

# Summary of trip_df
trip_df.info()

# Pattern of weather attributes datasets: pattern
pattern = 'kaggle_data\*.csv'

# Save all matching files with glob function: weather_files
weather_files = glob.glob(pattern)

# Aggregate all the datasets in a list
# by subsetting the 'datetime' and 'San Francisco' column of each dataset
weather_df_lst = [pd.read_csv(file, usecols=["datetime", "San Francisco"]) for file in weather_files]

# Concat all the dataframes in the weather_df_lst
weather_df = pd.concat(weather_df_lst, axis=1)
print(weather_df.head(2))

# Set the first 'datetime' column as index and
# Drop the other 'datetime' columns
weather_df = weather_df.set_index(weather_df.iloc[:, 0]).drop("datetime", axis=1)

# Convert the index to datetime
weather_df.index=pd.to_datetime(weather_df.index)

# Set the column names 
column_names = ['Humidity', 'Pressure', 'Temperature', 'Weather Description', 'Wind Direction', 'Wind Speed']
weather_df.columns = column_names
weather_df

# Subset the date interval [2014-09-01: 2015-08-31] (interval of bike trip dates)
weather_df = weather_df["2014-09-01": "2015-08-31"]
weather_df.head(3)

weather_df.info()

# Summary statistics of station_df
weather_df.describe().T

# Look at the missing values in 3 datasets
print(station_df.isna().values.any())
print(trip_df.isna().values.any())
print(weather_df.isna().values.any())

# Impute the single missing value in weather data
weather_df["Wind Direction"].fillna(method='ffill', inplace=True)

# Create dummy variables from 'Weather Description' column categories
dummies= pd.get_dummies(weather_df["Weather Description"], drop_first=True)
# Add the dummy variables to the weather_df and
# Drop the original features from the weather_df
weather_df= pd.concat([weather_df, dummies], axis=1).drop("Weather Description", axis=1)
weather_df.head(2)

fig, ax=plt.subplots(figsize=(16, 5))
sns.countplot(trip_df["Start Station"], ax=ax);

## Moved stations
moved_stations = [23, 24, 49, 69, 72]
new_stations1 = [85, 86, 87, 88, 89]
new_stations2 = [90]

# Create station list
station_lst = station_df["Id"].tolist()
station_moved = station_df[station_df["Id"].isin(moved_stations)]
station_new1 = station_df[station_df["Id"].isin(new_stations1)]
station_new2 = station_df[station_df["Id"].isin(new_stations2)]

# Create the list of coordinates
coord_list1 = list(zip(station_moved["Lat"], station_moved["Long"]))
coord_list2 = list(zip(station_new1["Lat"], station_new1["Long"]))
coord_list3 = list(zip(station_new2["Lat"], station_new2["Long"]))

# Zip the coordinates and the stations ids
moved_coords1 = list(zip(coord_list1, moved_stations))
moved_coords2 = list(zip(coord_list2, new_stations1))
moved_coords3 = list(zip(coord_list3, new_stations2))

## See the moved stations
stations_map = folium.Map(location=[37.56236, -122.150876], 
                          tiles='cartodbpositron',
                          zoom_start=10) 

# Add the pointers to the map by iterating over the coordinates
for point, station in moved_coords1: #
            marker = folium.Marker(location=point, popup='<i>Mt. Hood Meadows</i>', tooltip=station)
            marker.add_to(stations_map)    
            
for point, station in moved_coords2:
    marker2= folium.CircleMarker(location=point, popup='<i>Mt. Hood Meadows</i>', tooltip=station, radius=7)
    marker2.add_to(stations_map)
    
for point, station in moved_coords3:
    marker3= folium.CircleMarker(location=point, popup='<i>Mt. Hood Meadows</i>', tooltip=station, radius=7)
    marker3.add_to(stations_map)

stations_map

moved_stations=[23, 24, 49, 69, 72]
new_stations1=[85, 86, 87, 88, 89]
new_stations2=[90]

# Zip the moved stations and the new ones
replace_zip= list(zip(moved_stations, new_stations1))

# Replace the moved station values in 'Start Station' column with the new ones
for s1, s2 in replace_zip:
    trip_df.loc[trip_df["Start Station"]==s1, "Start Station"]=s2

# Replace the station 89 in 'Start Station' column with 90    
trip_df.loc[trip_df["Start Station"]==89, "Start Station"]=90

# Replace the moved station values in 'End Station' column with the new ones
for s1, s2 in replace_zip:
    trip_df.loc[trip_df["End Station"]==s1, "End Station"]=s2

# Replace the station 89 in 'End Station' column with 90         
trip_df.loc[trip_df["End Station"]==89, "End Station"]=90

# Drop the old stations
station_df= station_df[~station_df['Id'].isin(moved_stations + [89])]

# Check the duplicates in 'station_df' values
print(station_df["Id"].nunique()==len(station_df["Id"]))

# Check the duplicates in 'Trip ID' values in trip_df
print(trip_df["Trip ID"].nunique()==len(trip_df["Trip ID"]))

# Check the duplicates in 'Date' values in weather_df
print(weather_df.index.nunique()==len(weather_df))

# Use the 'Start Date' column to create the time features
trip_df["Month"]= trip_df["Start Date"].dt.month
trip_df["Day"]= trip_df["Start Date"].dt.dayofweek
trip_df["Hour Start"]= trip_df["Start Date"].dt.hour

# Create the trip duration column
trip_df["Duration"]= trip_df["End Date"]- trip_df["Start Date"]
# Convert the Duration into minutes
trip_df['Duration']=trip_df['Duration']/np.timedelta64(1,'m')

# Plot the boxplot
ax=sns.boxplot(x="Duration", data=trip_df, width=0.2, orient="h")

# Add the violinplot on the same figure
sns.violinplot(x="Duration", data=trip_df, bw=.2, orient="h")

# Add the label and title
ax.set(xlabel='Duration', title="Violin & Box Plot of Duration");

# Find the inter quantile of the duration
q1=trip_df['Duration'].quantile(0.25)
q3= trip_df['Duration'].quantile(0.75)
iqr = q3 - q1

print("Lower bound of outliers:", q1 - 1.5 * iqr)
print("Upper bound of outliers:", q3 + 1.5 * iqr, "\n")

# Only take the rows with duration is less then 120 min: trips
trips = trip_df[trip_df["Duration"] < 120]

# Print the percetage of removed data
print("The percentage of data removed:",\
      np.around((trip_df["Duration"] > 120).sum()/len(trip_df["Duration"])*100, decimals=2))

# Create a figure and two axes for boxplot and distplot
fig, (ax_box, ax_dist) = plt.subplots(nrows=2, 
                                      sharex=True, 
                                      figsize=(10, 6), 
                                      gridspec_kw={"height_ratios": (.15, .85)})
 
# Plot the boxplot on the ax_box axis
sns.boxplot(trips["Duration"], ax=ax_box)
# Plot the distplot on the ax_dist axis
sns.distplot(trips["Duration"], ax=ax_dist)
 # Remove xlabel for the boxplot
ax_box.set(xlabel='');

# Create the 'weekdays' and 'weekends' dataframes
weekdays=trips[trips["Start Date"].dt.weekday <5]
weekends=trips[trips["Start Date"].dt.weekday >=5]

def group_plot(df, lst, agg_func, labelx, labely, title):
    '''
    Takes a dataframe, a list of columns to groupby, an aggregation function for the groups, 
    3 strings: label of x-axis, label of y-axis and title returns a line plot
    Plots the data with labels and title.
    '''
    plt.style.use('fivethirtyeight')
    group= getattr(df.groupby(lst), agg_func)().unstack(lst[1])/1000
    ax= group.plot(figsize=(11,4), rot=70, fontsize=11)        
    ax.set_xlabel(labelx)
    ax.set_ylabel(labely)
    ax.set_title(title)
    plt.show()

labelx= "Month"
labely="Number of trips"
title= "Monthly weekdays trips by subscriber type (thousands)"
group_plot(weekdays,["Month", "Subscriber Type"], 'size', labelx, labely, title)

title= "Monthly weekend trips by subscriber type (thousands)"
group_plot(weekends,["Month", "Subscriber Type"], 'size', labelx, labely, title)

labelx= "Day"
title= "Daily trips by subscriber type (thousands)"
group_plot(trips,["Day", "Subscriber Type"], 'size', labelx, labely, title)

labelx= "Hour"
title= "Hourly weekdays trips by subscriber type (thousands)"
group_plot(weekdays,["Hour Start", "Subscriber Type"], 'size', labelx, labely, title)

title= "Hourly weekends trips by subscriber type (thousands)"
group_plot(weekends,["Hour Start", "Subscriber Type"], 'size', labelx, labely, title)

# Group the trips based on 'Hour Start' and 'Day'
# Aggregate by size of the groups
hourly=trips.groupby(["Hour Start", "Day"]).size().reset_index("Day").reset_index()
# Change the name of the columns
hourly.columns=["Hour Start", "Day", "Total"]
hourly.head(2)

# Plot the hourly rental means in the weekdays and weekends on poinplot
fig, ax=plt.subplots(figsize=(10, 5))

sns.pointplot(x="Hour Start", 
              y="Total", 
              hue="Day", 
              data=hourly, 
              errwidth=0,
              scale=.5, 
              ax=ax);

hue_labels = ["Sun", "Mon","Tue","Wed","Thu","Fri", "Sat"]
leg_handles = ax.get_legend_handles_labels()[0]
ax.legend(leg_handles, hue_labels);

ax.set_title("Hourly trips across the days");

# Resample trips_df by hour and get the size of each hour samples
hourly_total=trip_df.set_index("Start Date").resample("H").size()
# Create hourly_trips dataframe 
hourly_trips=hourly_total.to_frame(name="Total")
# Change the index name
hourly_trips.index.rename('Date', inplace=True)

# Lets add again time variables 
# Use the index to create the time features
hourly_trips["Month"]= hourly_trips.index.month
hourly_trips["Day"]= hourly_trips.index.dayofweek
hourly_trips["Hour"]= hourly_trips.index.hour

# Drop the 'Total' column
hourly_trips.drop("Total", axis=1, inplace=True)

# Add month dummy variables: m0, m1,..m12
dummy_months = pd.get_dummies(hourly_trips["Month"], drop_first=True)

# Add the dummy variables to the hourly_trips and
# Drop the original features from the hourly_trips
hourly_trips = pd.concat([hourly_trips, dummy_months], axis=1).drop("Month", axis=1)
# Rename the month columns
hourly_trips.columns = ["m"+ str(col) if type(col)==int else col for col in hourly_trips.columns]
hourly_trips.head(2)

# Boolean dictionary to map the 'True' values into 1 and "False" values into 0
bool_dct = {True:1, False:0}
# Only two flags would be enough to represent 3 partitions of the days
# hourly_trips["WD1"]=hourly_trips["Day"].isin([0,1,2]).map(bool_dct)
hourly_trips["WD2"] = hourly_trips["Day"].isin([3,4]).map(bool_dct)
hourly_trips["WKD"] = hourly_trips["Day"].isin([5,6]).map(bool_dct)
hourly_trips.head(2)

## Add day dummy variables: d0, d1,..d6
dummy_days = pd.get_dummies(hourly_trips["Day"], drop_first=True)

# Add the dummy variables to the hourly_trips and
# Drop the original features from the hourly_trips
hourly_trips= pd.concat([hourly_trips, dummy_days], axis=1).drop("Day", axis=1)
# Rename the day columns
hourly_trips.columns= ["d"+ str(col) if type(col)==int else col for col in hourly_trips.columns]
hourly_trips.head(2)

## Add hour dummy variables: h0, h1,..h23
dummy_hours= pd.get_dummies(hourly_trips["Hour"], drop_first=True)

# Add the dummy variables to the hourly_trips and
# Drop the original features from the hourly_trips
hourly_trips= pd.concat([hourly_trips, dummy_hours], axis=1).drop("Hour", axis=1)
# Rename the hour columns
hourly_trips.columns= ["h"+str(col) if type(col)==int else col for col in hourly_trips.columns]
hourly_trips.head(2)

# Use pandas US calendar
calendar = USFederalHolidayCalendar()
holidays = calendar.holidays(start=hourly_trips.index.min(), end=hourly_trips.index.max())
# Convert the calendar index into array
holidays = pd.to_datetime(holidays, format='%Y/%m/%d').date
# Create 'Holiday' column with binary values
hourly_trips["Holiday"]=[1 if day==True else 0 for day in hourly_trips.index.isin(holidays)]
hourly_trips.head(2)

features_df = pd.concat([hourly_trips, weather_df], axis=1).reset_index(drop=True)
features_df.head(2)

features_df.info()

# Create a list containing station ids
station_lst=station_df["Id"].tolist()

# Create departure columns' names
# 'd' stands for departures
departure_lst=[str(station) + "d" for station in station_lst]
# Create arrival columns' names
# 'a' stands for arrivals
arrival_lst=[str(station)+ "a" for station in station_lst]

# Create arrival and departure dictionaries
arrival_dict=dict(zip(station_lst, arrival_lst))
departure_dict=dict(zip(station_lst, departure_lst))

# Create targets dataframe: station_matrix
station_matrix=pd.DataFrame(columns=arrival_lst + departure_lst)

# Concat the trips_df and station_matrix: trips_extended
trips_extended=pd.concat([trip_df, station_matrix], sort=False).reset_index()

# Iterate over the trips_extended df and 
# assign 1 to the arrival column of the corresponding the Start Station on the same the row
# if on a row the 'Start Station' is 50 then the 50_a column value of the same row will be 1
# The same procedure is applied for the departure columns also ie
# if on a row the 'End Station' is 70 then the '70_d' column value of the same row will be 1
for i in trips_extended.index:
    trips_extended.at[i, arrival_dict[trips_extended.at[i,"Start Station"]]]=1
    trips_extended.at[i, departure_dict[trips_extended.at[i,"End Station"]]]=1

# After assigning the 1 value for each station's relevant arrival and departure columns
# fill the other NaN values with 0
trips_extended= trips_extended.fillna(0)
trips_extended.head(2)

# Sort the 'trips_extended' df by the date index
trips_extended= trips_extended.set_index("Start Date").sort_index()

# Create a dataframe with hourly resampling the 'trips_extended' df
# Sum the number of departures and arrivals in each station and each hour
trips_extended_hourly= trips_extended.resample('H').sum()

# See the departures distibution of station 70
trips_extended_hourly["70d"].plot(kind="hist");

# Create stations_hourly df by taking just the stations column
stations_hourly=trips_extended_hourly[arrival_lst+departure_lst].reset_index(drop=True)
stations_hourly.head(2)

stations_hourly.shape

# Targets dataset:y
y= stations_hourly
# Features dataset:X
X= features_df

# Find the starting indice of the last 10 percent for hold-out datasets
idx= int(len(y)* 0.90)
# Create the features train dataset: X_train
X_train= X.iloc[:idx, :]
# Create the features test dataset: X_test
X_test= X.iloc[idx:, :]
# Create the targets train dataset: y_train
y_train= y.iloc[:idx, :]
# Create the targets test dataset: y_test
y_test= y.iloc[idx: ,:]

print("X_train shape:", X_train.shape, "X_test shape:", X_test.shape)
print("y_train shape:", y_train.shape, "y_test shape:", y_test.shape)

# Add the time index to stations_hourly df
stations_hourly.index= trips_extended_hourly.index
# Take the averages of the columns
averages= stations_hourly.mean()
# Convert the averages into array
baseline= averages.values
# Stack vertically the baseline to match the y_test shape
baseline= np.tile(averages,(y_test.shape[0], 1))
# Check the shape of baseline
#print("baseline shape:", baseline.shape)

# Baseline errors are the mean of the differences between y_test and the baseline predictions
baseline_error= round (mean_squared_error(y_test.values, baseline)**0.5, 2)
print('Average baseline error: ', baseline_error)

def get_val_score_rf(X_tr, y_tr, dct):
    '''
    takes X_train, y_train and a dictionary containing hyperparameter inputs for random forest -> 
    returns the cross validation and X_train scores of the model. Both are calculated as 
    root mean squared errors 
    '''
    # Instantiate a RandomForestRegressor object by unpacking the parameters dictionary
    forest = RandomForestRegressor(**dct)
    # Instantiate TimeSeriesSplit object with 5 splits
    split=TimeSeriesSplit(n_splits=5)
    # Calculate the performance with cross-validation
    score = cross_val_score(forest, X_tr, y_tr, cv=split, scoring='neg_mean_squared_error', n_jobs=-1)
    # Take the sqrt all the splits 
    score= np.sqrt(abs(score))
    print("Cross-validation scores:", score,"\n","Cross-validation mean:", f'{np.mean(score):.2f}')
    # fit the model
    forest.fit(X_train, y_train)
    # Predict the train set to see if there is over fitting
    pred_train=forest.predict(X_train)
    print(f'RMSE Train:{np.sqrt(mean_squared_error(y_train, pred_train)):.2f}')  

# Create the params dictionary with basic hyperparameters, include the random_state
initial_params={'n_estimators':100,
        'max_depth':3,
        'min_samples_leaf':2,
        'random_state':37}

# Call the get_val_score_rf on the params dictionary
get_val_score_rf(X_train, y_train, initial_params)

rf1 = RandomForestRegressor(n_estimators=100, n_jobs=-1)
rf1.fit(X_train, y_train)
# Train the regressor
rf1.predict(X_test)

# Create a pandas Series of features importances: importances
# containing the feature names as index and their importances as values
importances_rf = pd.Series(data=rf1.feature_importances_ , index= X_train.columns)

# Get the sorted importance values: importance_sorted
importance_sorted=importances_rf.sort_values()

# Plot the sorted importance values by using horizontal bars
importance_sorted.plot(kind="barh", color='lightgreen', figsize=(10, 12));

# Sort the importances in descending order
importance_desc=importances_rf.sort_values(ascending=False)
sorted_features=importance_desc.index

# Cumulate the importances
cumulative_importances = np.cumsum(importance_desc)

# list of x locations for plotting
x_values = list(range(len(importance_desc)))

plt.figure(figsize=(14,5))
# Make a line graph
plt.plot(x_values, cumulative_importances, 'g-')
# Draw line at 99% of importance retained
plt.hlines(y = 0.99, xmin=0, xmax=len(importance_desc), color = 'r', linestyles = 'dashed')
# Format x ticks and labels
plt.xticks(x_values, sorted_features, rotation = 'vertical')
# Axis labels and title
plt.xlabel('Variable'); plt.ylabel('Cumulative Importance'); plt.title('Cumulative Importances');

# Add 1 because Python is zero-indexed
print('Number of features for 99.5% importance:', np.where(cumulative_importances > 0.995)[0][0] + 1)

# Extract the names of the most important features
important_features = [feature for feature in importance_desc[0:48].index]

# Create training and testing sets with only the important features
X_train_red_sk = X_train[important_features] # red_sk stands for reduced sklearn
X_test_red_sk = X_test[important_features]

# Sanity check
print('Important train features shape:', X_train_red_sk.shape)
print('Important test features shape:', X_test_red_sk.shape)

# Check the scores of the X_train_red_sk to compare with the first feature df X_train
get_val_score_rf(X_train_red_sk, y_train, initial_params)

# Instantiate a random forest object
rf_perm = RandomForestRegressor(n_estimators=100, n_jobs=-1)
# fit the model
rf_perm.fit(X_train, y_train)
# Apply the permutation function
imp = importances(rf_perm, X_test, y_test)
# plot the future importances
viz = plot_importances(imp)
viz.view()

# Cumulate the importances
cumulative_imp = np.cumsum(imp.values)
# Add 1 because Python is zero-indexed
print('Number of features for 99.9% importance:', np.where(cumulative_imp > 0.999)[0][0] + 1)

# Extract the names of the most important features
important_perm = [feature for feature in imp[0:20].index]
# Create training and testing sets with only the important features
X_train_perm = X_train[important_perm] # red_sk stands for reduced sklearn
X_test_perm= X_test[important_perm]

# Sanity check
print('Important train features shape:', X_train_perm.shape)
print('Important test features shape:', X_test_perm.shape)

# Check the scores of the X_train_perm to compare with 
# the first feature df(X_train) and the X_train_red_sk 
get_val_score_rf(X_train_perm, y_train, params)

# Create the base model to tune
forest_tune = RandomForestRegressor()
# print and see all the parameters of RandomForestRegressor
pprint(forest_tune.get_params())

# Number of trees in random forest
n_estimators = [int(x) for x in np.linspace(start = 100, stop =1000, num = 3)]
# Number of features to consider at every split
max_features = ['auto', 'sqrt']
# Maximum number of levels in tree
max_depth = [int(x) for x in np.linspace(3, 50, num = 2)]
max_depth.append(None)
# Minimum number of samples required to split a node
min_samples_split = [2, 10]
# Minimum number of samples required at each leaf node
min_samples_leaf = [2, 4]
# Method of selecting samples for training each tree
bootstrap = [True, False]

# Create the random grid
random_grid = {'n_estimators': n_estimators,
               'max_depth': max_depth,
               #'max_features': max_features,
               #'min_samples_split': min_samples_split,
               #'bootstrap': bootstrap
               'min_samples_leaf': min_samples_leaf}

pprint(random_grid)

# Create the splits for timeseries data 
split=TimeSeriesSplit(n_splits=2)
# Instantiate a Random Search object with parameters: n_iter = 10, cv = 2, 
forest_random = RandomizedSearchCV(estimator = forest_tune, 
                                   param_distributions = random_grid, 
                                   n_iter = 10, 
                                   cv = split,
                                   scoring='neg_mean_absolute_error',
                                   verbose=3, 
                                   random_state=37, 
                                   n_jobs = -1)

# Fit the random search model on reduced X_train_red_sk
#forest_random.fit(X_train_red_sk, y_train) #(pls uncomment to execute)

# Create the best params dictionary
best_params_random={'n_estimators':100,
                    'max_depth':50,
                    'min_samples_leaf':4,
                    'random_state':37}
# Check the scores of the X_train_red_sk
get_val_score_rf(X_train_red_sk, y_train, best_params_random)

# Create the parameter grid based on the results of random search 
param_grid = {'max_depth':[50, 60],
              'min_samples_leaf':[4, 5],
              'n_estimators':[100, 150]}
# Instantiate a RandomForestRegressor object
rf = RandomForestRegressor()
# Instantiate the grid search model
grid_search = GridSearchCV(estimator = rf, 
                           param_grid = param_grid, 
                           cv = split, 
                           scoring='neg_mean_absolute_error',
                           n_jobs = -1, 
                           verbose = 2)
pprint(param_grid)

# (pls uncomment to execute)
# Fit the grid search to the data 
#grid_search.fit(X_train_red_sk, y_train)
#grid_search.best_params_

# Create the best params dictionary
best_params_grid={'n_estimators':150,
                  'max_depth':60,
                  'min_samples_leaf':5,
                  'random_state':37}
# Check the scores 
get_val_score_rf(X_train_red_sk, y_train, best_params_grid)

# Instantiate a RandomForestRegressor object with the best parameters of random search
model_rf=RandomForestRegressor(**best_params_random)
# Fit the model 
model_rf.fit(X_train_red_sk, y_train)
# Predict the train set
pred_final=model_rf.predict(X_test_red_sk)

print(f'RMSE Test(hold-out):{np.sqrt(mean_squared_error(y_test, pred_final)):.2f}')

# Instantiate a RandomForestRegressor object with the best parameters of random search
model_rf=RandomForestRegressor(**best_params_random)
# Fit the model 
model_rf.fit(X_train_perm, y_train)
# Predict the train set
pred_final=model_rf.predict(X_test_perm)

print(f'RMSE Test(hold-out):{np.sqrt(mean_squared_error(y_test, pred_final)):.2f}')

# Save the model for future use 
filename = 'model_rf.model'
joblib.dump(model_rf, filename)

# Check if the model properly saved
model=joblib.load(filename)
pred_joblib=model.predict(X_test_perm)
print(f'RMSE Test(hold-out):{np.sqrt(mean_squared_error(y_test, pred_joblib)):.2f}')

def net_change(model, df):
    '''
    Takes our model and dataframe to predict returns->
    a dataframe with the half number of columns of the predictons made by our model
    by taking the difference of the departures and arrival columns
    '''
    # Create a dataframe from predictions
    predictions_df = pd.DataFrame(model.predict(df))
    # Create list of column names
    columns_lst = list(predictions_df.columns)
    # Create a list to aggregate the differences of departures and arrivals
    difference_cols =[]
    # Iterate over the column names list till the half point
    for idx in range(len(columns_lst)//2):
        # Collect the column differences in the list
        difference_cols.append(predictions_df[idx+70]-predictions_df[idx])
    # At the end of the loop concat all the columns at once
    net_rate_df = pd.concat(difference_cols, axis=1)
    return net_rate_df

# Call the net_change on our model with test data
net_change(model_rf, X_test_perm).head(2)

# Create an empty list to aggregate the net change in every row(hour)
columns_lst=[]

## Loop over stations_hourly df columns to find the net change (arrivals-departures)
for station in stations_hourly.columns:
    if station[-1]=="a":
        station_arr = station
        station_dep = station[:-1]+ "d"
        # Append the column difference in the list
        columns_lst.append((stations_hourly[station_arr]- stations_hourly[station_dep]))

# Create the change_df by concatenating the columns in the column_lst
net_df = pd.concat(columns_lst, axis=1) 
# Create the column names of the net_df
net_cols =[]
for station in stations_hourly.columns:
    if station[-1]=="a":
        net_cols.append(station[:-1]+ "c")
# Add the column names        
net_df.columns = net_cols       
net_df.head(2)

net_df.shape

# Get the predictions of the test 
predicted_net_change= net_change(model_rf, X_test_perm)
# Actual net changes
actual_net_change=y_test_net
# Get the Root Mean Squared Errors of differences
rmse_model_rf=round(mean_squared_error(y_test_net, predicted_net_change)**0.5, 2)
rmse_model_rf

# Targets dataset:y
y_net= net_df

# Find the starting indice of the last 10 percent for hold-out datasets
idx= int(len(y_net)* 0.90)

# Create the targets train dataset: y_train_net
y_train_net= y_net.iloc[:idx, :]
# Create the targets test dataset: y_test_net
y_test_net= y_net.iloc[idx: ,:]

# Check the shapes
print("X_train shape:", X_train.shape, "X_test shape:", X_test.shape)
print("y_train shape:", y_train_net.shape, "y_test shape:", y_test_net.shape)

# Take the averages of the columns
averages= net_df.mean()
# Convert the averages into array
baseline= averages.values
# Stack vertically the baseline to match the y_test shape
baseline= np.tile(averages,(y_test_net.shape[0], 1))
# Check the shape of baseline
#print("baseline shape:", baseline.shape)
# Baseline errors are the mean of the differences between y_test and the baseline predictions
baseline_error= round (mean_squared_error(y_test_net.values, baseline)**0.5, 2)
print('Average baseline error: ', baseline_error)

# Create the params dictionary with basic hyperparameters, include the random_state
params={'n_estimators':100,
        'max_depth':3,
        'min_samples_leaf':2,
        'random_state':37}

# Call the get_val_score_rf on X_train_net the params dictionary
get_val_score_rf(X_train, y_train_net, params)

# Instantiate a random forest object
rf_perm = RandomForestRegressor(n_estimators=100, n_jobs=-1)
# fit the model
rf_perm.fit(X_train, y_train_net)
# Apply the permutation function
imp = importances(rf_perm, X_test, y_test_net)
# Cumulate the importances
cumulative_imp = np.cumsum(imp.values)
# Add 1 because Python is zero-indexed
print('Number of features for 99.9% importance:', np.where(cumulative_imp > 0.999)[0][0] + 1)

# Extract the names of the most important features. Let's keep 20 of them
important_net = [feature for feature in imp[0:20].index]
# Create training and testing sets with only the important features
X_train_red_net = X_train[important_perm] # red_sk stands for reduced sklearn
X_test_red_net= X_test[important_perm]

# Sanity check
print('Important train features shape:', X_train_red_net.shape)
print('Important test features shape:', X_test_red_net.shape)

# Parameters from RandomizedSearchCV
params_random={'n_estimators':100,
               'max_depth':50,
               'min_samples_leaf':4,
               'random_state':37}

# Call the get_val_score_rf on X_train_net the params dictionary
get_val_score_rf(X_train_red_net, y_train_net, params_random)

params_grid={'n_estimators':100,
             'max_depth':60,
             'min_samples_leaf':5,
             'random_state':37}

# Call the get_val_score_rf on X_train, y_train_net and the params dictionary
get_val_score_rf(X_train_red_net, y_train_net, params_grid)