# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Predicting the Sale Price of Bulldozers Using Machine Learning
#
# For this dataset, the problem we're trying to solve, or better, the question we're trying to answer is,
#
# **How well can we predict the future sale price of a bulldozer, given its characteristics previous examples of how much similar bulldozers have been sold for?**
#
# ## Data
#
# The data is downloaded from the Kaggle Bluebook for Bulldozers competition: https://www.kaggle.com/c/bluebook-for-bulldozers/data
#
# There are 3 main datasets:
#
# * Train.csv is the training set, which contains data through the end of 2011.
# * Valid.csv is the validation set, which contains data from January 1, 2012 - April 30, 2012 You make predictions on this set   throughout the majority of the competition. Your score on this set is used to create the public leaderboard.
# * Test.csv is the test set, which won't be released until the last week of the competition. It contains data from May 1, 2012 - November 2012. Your score on the test set determines your final rank for the competition.
#
# ## Evaluation
#
# The evaluation metric for this competition is the RMSLE (root mean squared log error) between the actual and predicted auction prices.
#
# For more on the evaluation of this project check: https://www.kaggle.com/c/bluebook-for-bulldozers/overview/evaluation
#
# Note: The goal for most regression evaluation metrics is to minimize the error. For example, our goal for this project will be to build a machine learning model which minimises RMSLE.
#
# ## Features
#
# For this dataset, Kaggle provide a data dictionary which contains information about what each attribute of the dataset means. You can download this file directly from the Kaggle competition page (account required) or view it on Google Sheets.
#

# ## EDA

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

df = pd.read_csv("data/TrainAndValid.csv", low_memory=False)

df.info()

df.isna().sum()

df.columns

fig, ax = plt.subplots()
ax.scatter(df['saledate'][:1000], df['SalePrice'][:1000]);

df['saledate'][:1000]

df.SalePrice.plot(kind='hist');

# ### Parsing Dates
#
# Here, we have to convert the 'saledate' column from object type to datetime type in order to make it actual time so that we can work with it.

df = pd.read_csv("data/TrainAndValid.csv", low_memory=False, parse_dates=['saledate'])

df.saledate.dtype

df.saledate[:1000]

fig, ax = plt.subplots()
ax.scatter(df['saledate'][:1000], df['SalePrice'][:1000]);

df.head()

df.head().T

df.saledate.head(20)

# ### Sorting DataFrame by Dates in Ascending Order

df.sort_values(by=['saledate'], inplace=True)
df.saledate.head()

# ### Make a copy of original DataFrame

df_tmp = df.copy()

# ### Feature Engineering
#
# Basically, we are splitting the 'saledate; column into 5 new features called day, month, year, sale day of week (The day of the week with Monday=0, Sunday=6) & sale day of year (The ordinal day of the year) 

# +
df_tmp['saleYear'] = df_tmp.saledate.dt.year
df_tmp['saleMonth'] = df_tmp.saledate.dt.month
df_tmp['saleDay'] = df_tmp.saledate.dt.day
df_tmp['saleDayOfWeek'] = df_tmp.saledate.dt.dayofweek
df_tmp['saleDayOfYear'] = df_tmp.saledate.dt.dayofyear

df_tmp.drop('saledate', axis=1, inplace=True)
# -

df_tmp.head().T

df_tmp.state.value_counts()

df_tmp.head()

len(df_tmp)

# ## Modelling

# This will not work as we have missing values along with various object data types like strings which needs
# to be converted into categories for assigning numbers.

# +
# from sklearn.ensemble import RandomForestClassifier

# model = RandomForestClassifier(n_jobs=-1, random_state=42)

# model.fit(df_tmp.drop('SalePrice', axis=1), df_tmp['SalePrice'])
# -

df_tmp.info();

df_tmp["UsageBand"].dtype

# ## Converting String to Categories
#
# We can check the different datatypes compatible with pandas here: https://pandas.pydata.org/pandas-docs/stable/reference/general_utility_functions.html#data-types-related-functionality

pd.api.types.is_string_dtype(df_tmp["UsageBand"])

# +
#Columns with strings
#Label is column name and content is column contents

for label, content in df_tmp.items():
    if pd.api.types.is_string_dtype(content):
        print(label)

# +
# If you're wondering what df_tmp.items() does, here is an example:

random_dict = {"key1" : "Hello",
               "key2" : "World"}

for key, value in random_dict.items():
    print(f"The key is: {key}", 
          f"The value is: {value}")
# -

for label, content in df_tmp.items():
    if pd.api.types.is_string_dtype(content):
        df_tmp[label] = content.astype("category").cat.as_ordered()

df_tmp.info();

df_tmp.state.cat.categories

df_tmp.state.cat.codes

# With the help of Pandas Categories, we now have a way to access all of our data in the form of numbers.
#
# But in order to train/fit model without errors we still have to handle missing data.

# ### Save Preprocessed Data as CSV

#Export current temp dataframe
df_tmp.to_csv("data/train_tmp.csv", index=False)

# +
#Import preprocessed data
df_tmp =  pd.read_csv("data/train_tmp.csv", low_memory=False)

df_tmp.head().T
# -

# ## Fill in Missing Values

# +
#Check missing data
ratio = df_tmp.isna().sum()/len(df_tmp)

df_tmp.isna().sum()
# -

# ### Fill Numerical Missing Values First

for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        print(label)

# Check which numeric columns have null values
for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        if pd.isnull(content).sum():
            print(label)

# Fill the missing numeric rows with the median and then create a new column which says if the content was filled with 
# missing values or it has actual value
for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        if pd.isnull(content).sum():
            # Add a binary column which tells if the data was missing or not
            df_tmp[label + "_is_missing"] = pd.isnull(content)
            # Fill missing values with media since it is more robust to outliers than mean
            df_tmp[label] = content.fillna(content.median())

# Why add a binary column indicating whether the data was missing or not?
#
# We can easily fill all of the missing numeric values in our dataset with the median. However, a numeric value may be missing for a reason. In other words, absence of evidence may be evidence of absence. Adding a binary column which indicates whether the value was missing or not helps to retain this information.

# Check if there are any null values
for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        if pd.isnull(content).sum():
            print(label)

# Check howmany values were missing
df_tmp.auctioneerID_is_missing.value_counts()

# ### Filling Categorical Values After Converting Them to Numbers

# Check columns which are NOT numeric
for label, content in df_tmp.items():
    if not pd.api.types.is_numeric_dtype(content):
        print(label)

# Turn categorical variables into numbers
for label, content in df_tmp.items():
    # Check columns which are NOT numeric
    if not pd.api.types.is_numeric_dtype(content):
        # Add binary column to indicate whether sample had missing values
        df_tmp[label + "_is_missing"] = pd.isnull(content)
        # Convert the categories to codes using the .codes property of Categorical class and add +1
        # We add +1 as pandas encodes missing values of categories as -1. But we don't want -ve values so we make all missing
        # values to 0
        df_tmp[label] = pd.Categorical(content).codes + 1

df_tmp.info()

df_tmp.isna().sum()

df_tmp.head().T

# ## Time for Modelling
#
# Now all of our data is numeric and there are no missing values, we should be able to build a machine learning model!
#
# Let's reinstantiate our trusty [RandomForestRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html).
#
# This will take a few minutes which is too long for interacting with it. So what we'll do is create a subset of rows to work with.

# +
# %%time

from sklearn.ensemble import RandomForestRegressor
# Instantiating Model
model = RandomForestRegressor(n_jobs=-1, random_state=42)

#Training the Model
# model.fit(df_tmp.drop('SalePrice', axis=1), df_tmp['SalePrice'])
# -

# %%time
# Score the model
model.score(df_tmp.drop("SalePrice", axis=1), df_tmp['SalePrice'])

# Why is the metric not reliable?
#
# Because we have scored it with the training data. In order to score a model we should always use validation set and test set.

# ### Splitting Data into Train & Valid Sets 

df_tmp.head()

# According to the [Kaggle data page](https://www.kaggle.com/c/bluebook-for-bulldozers/data), the validation set and test set are split according to dates.
#
# This makes sense since we're working on a time series problem.
#
# E.g. using past events to try and predict future events.
#
# Knowing this, randomly splitting our data into train and test sets using something like `train_test_split()` wouldn't work.
#
# Instead, we split our data into training, validation and test sets using the date each sample occured.
#
# In our case:
# * Training = all samples up until 2011
# * Valid = all samples form January 1, 2012 - April 30, 2012
# * Test = all samples from May 1, 2012 - November 2012
#
# For more on making good training, validation and test sets, check out the post [How (and why) to create a good validation set](https://www.fast.ai/2017/11/13/validation-sets/) by <NAME>.

df_tmp.saleYear.value_counts()

# +
# Split data into train and valid sets
df_val = df_tmp[df_tmp.saleYear == 2012]
df_train = df_tmp[df_tmp.saleYear != 2012]

len(df_val), len(df_train)

# +
# Split data into X and Y
# df_train["SalePrice"] is 
X_train, y_train = df_train.drop("SalePrice", axis=1), df_train["SalePrice"]
X_valid, y_valid = df_val.drop("SalePrice", axis=1), df_val.SalePrice

X_train.shape, y_train.shape, X_valid.shape, y_valid.shape
# -

# ### Building An Evaluation Function

# According to Kaggle for the Bluebook for Bulldozers competition, [the evaluation function](https://www.kaggle.com/c/bluebook-for-bulldozers/overview/evaluation) they use is root mean squared log error (RMSLE).
#
# **RMSLE** = Generally you don't care if you're off by $10. But you'd care if you were off by 10%. You care more about ratios/percentages rather than differences. **MAE** (mean absolute error) is more about exact differences.
#
# It's important to understand the evaluation metric you're going for.
#
# Since Scikit-Learn doesn't have a function built-in for RMSLE, we'll create our own.
#
# We can do this by taking the square root of Scikit-Learn's [mean_squared_log_error](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_log_error.html#sklearn.metrics.mean_squared_log_error) (MSLE). MSLE is the same as taking the log of mean squared error (MSE).
#
# We'll also calculate the MAE and R^2 for fun.

# +
# Create an evaluation function (this competition uses Root Mean Square Log Error)
from sklearn.metrics import mean_squared_log_error, mean_absolute_error, r2_score

def rmsle(y_test, y_preds):
    """
    Caculates root mean squared log error between predictions and
    true labels.
    
    """
    return np.sqrt(mean_squared_log_error(y_test, y_preds))

def show_scores(model):
    train_preds = model.predict(X_train)
    valid_preds = model.predict(X_valid)
    scores = {"Training MAE": mean_absolute_error(y_train, train_preds),
              "Valid MAE": mean_absolute_error(y_valid, valid_preds),
              "Training RMSLE": rmsle(y_train, train_preds),
              "Valid RMSLE": rmsle(y_valid, valid_preds), 
              "Training R^2": model.score(X_train, y_train), #can also use r2_score(y_train, train_preds)
              "Valid R^2": model.score(X_valid, y_valid)}    #can also use r2_score(y_valid, val_preds)
    return scores


# -

# ### Testing our model on a subset (to tune the hyperparameters)
#
# Retraing an entire model would take far too long to continuing experimenting as fast as we want to.
#
# So what we'll do is take a sample of the training set and tune the hyperparameters on that before training a larger model.
#
# If you're experiments are taking longer than 10-seconds (give or take how long you have to wait), you should be trying to speed things up. You can speed things up by sampling less data or using a faster computer.

# +
# # %%time
# # This takes a long time............
# # Retrain our model on training data since the model was previously trained on both train & valid data.
# model.fit(X_train, y_train)
# show_scores(model)
# -

len(X_train)

# Depending on our computer, making calculations on ~400,000 rows may take a while...
#
# Let's alter the number of samples each `n_estimator` in the [`RandomForestRegressor`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html) see's using the `max_samples` parameter.

# Change max samples in RandomForestRegressor
model = RandomForestRegressor(n_jobs=-1, max_samples=10000)

# Setting `max_samples` to 10000 means every `n_estimator` (default 100) in our `RandomForestRegressor` will only see 10000 random samples from our DataFrame instead of the entire 400,000.
#
# In other words, we'll be looking at 40x less samples which means we'll get faster computation speeds but we should expect our results to worsen (simple the model has less samples to learn patterns from).

# %%time
# Cutting down the max no. of samples each tree can see reduces training time
model.fit(X_train, y_train)

show_scores(model)

# ### Hyperparameter tuning RandomizedSearchCV

# You can increase `n_iter` to try more combinations of hyperparameters but in our case, we'll try 20 and see where it gets us.
#
# Remember, we're trying to reduce the amount of time it takes between experiments.

# +
# %%time
from sklearn.model_selection import RandomizedSearchCV

# Different RandomForestClassifer hyperparameters
rf_grid = {
        "n_estimators": np.arange(10, 100, 10),
        "max_depth": [None, 3, 5, 10],
        "min_samples_split": np.arange(2, 20, 2),
        "min_samples_leaf": np.arange(1, 20,2),
        "max_features": [0.5, 1, "sqrt", "auto"],
        "max_samples": [10000]
          }

rs_model = RandomizedSearchCV(RandomForestRegressor(),
                              param_distributions=rf_grid,
                              n_iter=20,
                              cv=5,
                              n_jobs=-1,
                              verbose=True,
                              random_state=42)

rs_model.fit(X_train, y_train)
# -

# Find the best parameters from the RandomizedSearch 
rs_model.best_params_

# Evaluate the RandomizedSearch Model
show_scores(rs_model)

# ### Train a model with the best parameters
#
# In a model I prepared earlier, I tried 100 different combinations of hyperparameters (setting `n_iter` to 100 in `RandomizedSearchCV`) and found the best results came from the ones you see below.
#
# **Note:** This kind of search on my computer (`n_iter` = 100) took ~2-hours. So it's kind of a set and come back later experiment.
#
# We'll instantiate a new model with these discovered hyperparameters and reset the `max_samples` back to its original value.

# +
# %%time
# Most ideal hyperparameters
ideal_model = RandomForestRegressor(n_estimators=90,
                                    min_samples_leaf=1,
                                    min_samples_split=14,
                                    max_features=0.5,
                                    n_jobs=-1,
                                    max_samples=None,
                                    random_state=42) # random state to reproduce same results

ideal_model.fit(X_train, y_train)
# -

# Show scores for ideal_model (trained on all data)
show_scores(ideal_model)

# Show scores for rs_model (only trained on ~10,000 examples)
show_scores(rs_model)

# With these new hyperparameters as well as using all the samples, we can see an improvement to our models performance.
#
# You can make a faster model by altering some of the hyperparameters. Particularly by lowering `n_estimators` since each increase in `n_estimators` is basically building another small model.
#
# However, lowering of `n_estimators` or altering of other hyperparameters may lead to poorer results.

# +
# %%time
# Faster model than ideal model
fast_model = RandomForestRegressor(n_estimators=40,
                                   min_samples_leaf=3,
                                   max_features=0.5,
                                   n_jobs=-1,
                                   random_state=42) # random state to reproduce same results

fast_model.fit(X_train, y_train)
# -

show_scores(fast_model)

# ### Make predictions on test data
#
# Now we've got a trained model, it's time to make predictions on the test data.
#
# Remember what we've done.
#
# Our model is trained on data prior to 2011. However, the test data is from May 1 2012 to November 2012.
#
# So what we're doing is trying to use the patterns our model has learned in the training data to predict the sale price of a Bulldozer with characteristics it's never seen before but are assumed to be similar to that of those in the training data.

# +
df_test = pd.read_csv("data/Test.csv", parse_dates=["saledate"])

df_test.head()


# +
# Now if we try to predict, we'll get error since the model was trained on training dataset which has different no. of 
# columns/features. Moreover, the test dataset has missing values and also many features need to be convert from object to
# categorical codes

#test_preds = ideal_model.predict(df_test)
# -

# ### Preprocessing the data
#
# Our model has been trained on data formatted in the same way as the training data.
#
# This means in order to make predictions on the test data, we need to take the same steps we used to preprocess the training data to preprocess the test data.
#
# Remember: Whatever you do to the training data, you have to do to the test data.
#
# Let's create a function for doing so (by copying the preprocessing steps we used above).

def preprocess_data(df):
    # Add datetime parameters for saledate
    df["saleYear"] = df.saledate.dt.year
    df["saleMonth"] = df.saledate.dt.month
    df["saleDay"] = df.saledate.dt.day
    df["saleDayOfWeek"] = df.saledate.dt.dayofweek
    df["saleDayOfYear"] = df.saledate.dt.dayofyear
    
    # Drop orginal saledate
    df.drop("saledate", axis=1, inplace=True)
    
    # Fill numeric rows with the median
    for label, content in df.items():
        if pd.api.types.is_numeric_dtype(content):
            if pd.isnull(content).sum():
                # Add a binary column which tells us if the data was missing or not
                df[label+"_is_missing"] = pd.isnull(content)
                # Fill missing numeric values with median
                df[label] = content.fillna(content.median())
        
        # Turn categorical variables into numbers and fill missing
        if not pd.api.types.is_numeric_dtype(content):
            # Add binary column to indicate whether sample had missing value
            df[label+"_is_missing"] = pd.isnull(content)
            # Turn categories into numbers and add +1
            df[label] = pd.Categorical(content).codes+1
            
    return df


# +
df_test = preprocess_data(df_test)

df_test.head()
# -

X_train.head()

# +
# # We still cannot get y_preds from test.csv as our test dataset (after preprocessing) has 101 columns whereas, 
# # our training dataset (X_train) has 102 columns (after preprocessing).
# y_preds = ideal_model.predict(df_test)
# -

# Let's find the difference between two dataframes

set(X_train.columns) - set(df_test.columns)

# In this case, it's because the test dataset wasn't missing any `auctioneerID` fields.
#
# To fix it, we'll add a column to the test dataset called `auctioneerID_is_missing` and fill it with `False`, since none of the `auctioneerID` fields are missing in the test dataset.

# Match test dataset columns to training dataset
df_test["auctioneerID_is_missing"] = False
df_test.head()

# Make predictions on the test dataset using the best model
test_preds = ideal_model.predict(df_test)

# When looking at the [Kaggle submission requirements](https://www.kaggle.com/c/bluebook-for-bulldozers/overview/evaluation), we see that if we wanted to make a submission, the data is required to be in a certain format. Namely, a DataFrame containing the `SalesID` and the predicted `SalePrice` of the bulldozer.

# Create Dataframe compatible with Kaggle data submission requirements
df_preds = pd.DataFrame()
df_preds["SalesID"] = df_test["SalesID"]
df_preds["SalePrice"] = test_preds
df_preds

# Export prediction data
df_preds.to_csv("data/test_predictions.csv", index=False)

# ## Feature Importance
#
# Since we've built a model which is able to make predictions. The people you share these predictions with (or yourself) might be curious of what parts of the data led to these predictions.
#
# This is where **feature importance** comes in. Feature importance seeks to figure out which different attributes of the data were most important when it comes to predicting the **target variable**.
#
# In our case, after our model learned the patterns in the data, which bulldozer sale attributes were most important for predicting its overall sale price?
#
# Beware: the default feature importances for random forests can lead to non-ideal results.
#
# To find which features were most important of a machine learning model, a good idea is to search something like "\[MODEL NAME\] feature importance".
#
# Doing this for our `RandomForestRegressor` leads us to find the `feature_importances_` attribute.
#
# Let's check it out.

# Find feature importance of our best model
ideal_model.feature_importances_


# Helper functions for plotting feature importance
def plot_features(columns, importances, n=20):
    df = (pd.DataFrame({"features": columns,
                        "feature_importances": importances})
         .sort_values("feature_importances", ascending=False)
         .reset_index(drop=True))
    
    #Plot the dataframe
    fig, ax = plt.subplots()
    ax.barh(df["features"][:n], df["feature_importances"][:20])
    ax.set_xlabel("Feature Importance")
    ax.set_ylabel("Features")
    ax.invert_yaxis()


plot_features(X_train.columns, ideal_model.feature_importances_)

sum(ideal_model.feature_importances_)

df.ProductSize.isna().sum()

df.ProductSize.value_counts()

df.Turbocharged.value_counts()

df.Thumb.value_counts()

df["Enclosure"].value_counts()

# **Question to finish:** Why might knowing the feature importances of a trained machine learning model be helpful?
#
# **Final challenge/extension:** What other machine learning models could you try on our dataset? 
#  https://scikit-learn.org/stable/tutorial/machine_learning_map/index.html Check out the regression section of this map, or try to look at something like CatBoost.ai or XGBooost.ai.

StarcoderdataJupyterScriptsDedupFiltered