# Common imports
import numpy as np
import os

# Import and settings for plotings figures
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)

# Once the data have been downloaded and stored at ./datasets/housesales
house_sales_path = os.path.join("datasets", "housesales")

# Import panda
import pandas as pd
# Create a function to load the dataset into a panda object
def load_house_sales(path=house_sales_path):
    csv_path = os.path.join(path, "kc_house_data.csv")
    return pd.read_csv(csv_path)

# Call our function to load the dataset
housesales = load_house_sales()

# This step is done for simulating that we have same NaN values.
# In a real project, it does need to be done
shuffled_indices = np.random.permutation(len(housesales))
na_size = int(len(housesales) * 2 / 1000)
na_indices = shuffled_indices[:na_size]
housesales.loc[na_indices, 'sqft_living'] = np.nan

# View null values
housesales[housesales.isnull().any(axis=1)].head(5)

# This step is done for simulating that we have a category column
# In a real project, it does need to be done
shuffled_indices = np.random.permutation(len(housesales))
cats = ['A', 'B', 'C', 'D', 'E']
cat_size = int(len(housesales) / 5)
rem = len(housesales) % 5
for i in range(5):
    plus = 0
    if i == 4:
        plus = rem
    indices = shuffled_indices[i*cat_size:((i+1)*cat_size + plus)]
    housesales.loc[indices, "new_cat"] = cats[i]

# Explore first rows of our dataset
housesales.head()

# Show info of our dataset
housesales.info()

# Count values for discrete columns
housesales["grade"].value_counts()

# Obtain statistics for numerical columns
housesales.describe()

# Plot the histogram for each numerical value
import matplotlib.pyplot as plt
housesales.hist(bins=50, figsize=(20,15))
plt.show()

# Make this notebook's output identical at every run
np.random.seed(10)

# Split the data into train and test dataset by using our own function
# For illustration only. Sklearn has train_test_split()

import numpy as np
def split_train_test(data, test_ratio):
    shuffled_indices = np.random.permutation(len(data))
    test_set_size = int(len(data) * test_ratio)
    test_indices = shuffled_indices[:test_set_size]
    train_indices = shuffled_indices[test_set_size:]
    return data.iloc[train_indices], data.iloc[test_indices]

train_set, test_set = split_train_test(housesales, 0.2)
print(len(train_set), "train +", len(test_set), "test")

# Split the data into train and test dataset by using sklearn tools
from sklearn.model_selection import train_test_split

train_set, test_set = train_test_split(housesales, test_size=0.2, random_state=10)
print(len(train_set), "train +", len(test_set), "test")

# Plot data points by using their latitute and longitude
housesales = train_set.copy()
housesales.plot(kind="scatter", x="long", y="lat")

# Since data belongs to King County, USA, we can obtain its map
# Load and plot that map

img_path = os.path.join("datasets", "housesales", "King_County_Washington.png")
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
img=mpimg.imread(img_path)
imgplot = plt.imshow(img)
plt.show()

# Add alpha(transparency) to see density of data points
housesales.plot(kind="scatter", x="long", y="lat", alpha=0.05)

# Use color and mark size to point out price and grade features
housesales.plot(kind="scatter", x="long", y="lat", alpha=0.3,
    s=housesales["grade"]*10, label="grade", figsize=(12,8),
    c="price", cmap=plt.get_cmap("jet"), colorbar=True, 
    sharex=False
    )
plt.legend()

# Integrate its map to the previous plot

import matplotlib.pyplot as plt
import matplotlib.image as mpimg

housesales.plot(kind="scatter", x="long", y="lat", alpha=0.2,
    s=housesales["grade"]*10, label="grade", figsize=(16,8),
    c="price", cmap=plt.get_cmap("jet"), colorbar=False, 
    sharex=False
    )

plt.ylabel("Latitude", fontsize=14)
plt.xlabel("Longitude", fontsize=14)


king_img = mpimg.imread(img_path)
plt.imshow(king_img, extent=[-122.54, -120.80, 47.04, 47.85], alpha=0.4,
           cmap=plt.get_cmap("jet"))


prices = housesales["price"]
tick_values = np.linspace(prices.min(), prices.max(), 11)

cbar = plt.colorbar()
cbar.ax.set_yticklabels(["$%dk"%(round(v/1000)) for v in tick_values], fontsize=14)
cbar.set_label('Price', fontsize=16)

plt.legend()

# Generate the correlation matrix
corr_matrix = housesales.corr()

# Show features more correlated to the 'price' feature
corr_matrix['price'].sort_values(ascending=False)

# Plot scatter plots among the most correlated features
# For the main diagonal, pandas displays a histogram instead of straight line(variable against itself)
from pandas.plotting import scatter_matrix

attributes = ["price", "sqft_living", "grade","sqft_above", "bathrooms"]
fig = scatter_matrix(housesales[attributes], figsize=(12, 8))

# Plot the scatter plot between two especific features
housesales.plot(kind="scatter", x="sqft_living", y="price", alpha=0.1)
plt.axis([0, 7000, 0, 2000000])

# drop labels for training set
housesales = train_set.drop("price", axis=1) 
# Assign labels to a new variable
housesales_labels = train_set["price"].copy()

# Look at incomplete features
incomplete_rows = housesales[housesales.isnull().any(axis=1)]
incomplete_rows.head(5)

# option 1: Drop rows containing NaN values
incomplete_rows.dropna(subset=["sqft_living"])

# option 2: Drop feature containing NaN values
incomplete_rows.drop("sqft_living", axis=1).head()

# option 3: Fill out features containing NaN values with some criteria(i.e. median)
median = housesales["sqft_living"].median()
incomplete_rows["sqft_living"].fillna(median, inplace=True)
incomplete_rows.head(5)

# Option 3B: Fill out features containing NaN values by using SimpleImputer
from sklearn.impute import SimpleImputer
imputer = SimpleImputer(strategy="median")

# Remove the id, date and new_cat attribute because median can only be calculated on numerical attributes
housesales_num = housesales.drop(['id','date', 'new_cat'], axis=1)
housesales_num.head()

# Fit the Simple Imputer
imputer.fit(housesales_num)

# Query the Simple Imputer statistics
imputer.statistics_

# Show median by other mean to contrast results
housesales_num.median().values

# Apply the transformation to create our new X
X = imputer.transform(housesales_num)

# Create a dataframe from X (a numpy array object)
housesales_tr = pd.DataFrame(X, columns=housesales_num.columns,
                          index = list(housesales.index.values))

# Show that columns containing NaN values are filled out
housesales_tr.loc[incomplete_rows.index.values].head(5)

# Query the strategy used
imputer.strategy

# Create a datafrane by resetting the indexes
housesales_tr = pd.DataFrame(X, columns=housesales_num.columns)
housesales_tr.head()

# Show the categorial feature
housesales_cat = housesales['new_cat']
housesales_cat.head(10)

# (Option 1) Apply Ordinal Encoder to our categorical feature
from sklearn.preprocessing import OrdinalEncoder
ordinal_encoder = OrdinalEncoder()
housesales_cat_encoded = ordinal_encoder.fit_transform(housesales_cat.values.reshape(-1,1))
housesales_cat_encoded[:10]

# Query those categories found
ordinal_encoder.categories_

# (Option 2) Apply OneHotEnconder to our categorical feature
from sklearn.preprocessing import OneHotEncoder
cat_encoder = OneHotEncoder()
housesales_cat_1hot = cat_encoder.fit_transform(housesales_cat.values.reshape(-1,1))
housesales_cat_1hot

# By default, the OneHotEncoder class returns a sparse array, but we can convert it to a dense array
# if needed by calling the toarray() method
housesales_cat_1hot.toarray()

# (Option 2B) Alternatively, you can set sparse=False when creating the OneHotEncoder:
cat_encoder = OneHotEncoder(sparse=False)
housesales_cat_1hot = cat_encoder.fit_transform(housesales_cat.values.reshape(-1,1))
housesales_cat_1hot

# Query those categories found
cat_encoder.categories_

# Query what the columns are
housesales.columns

# Show the initial shape
housesales.shape

# Create a function to add new features a apply this one

from sklearn.preprocessing import FunctionTransformer

# get the right column indices: safer than hard-coding indices
bedrooms_ix, bathrooms_ix, floors_ix = [
    list(housesales.columns).index(col)
    for col in ("bedrooms", "bathrooms", "floors")]

def add_extra_features(X, add_bathrooms_per_floors=True):
    bedrooms_per_floor = X[:, bedrooms_ix] / X[:, floors_ix]
    if add_bathrooms_per_floors:
        bathrooms_per_floor = X[:, bathrooms_ix] / X[:, floors_ix]
        return np.c_[X, bedrooms_per_floor, bathrooms_per_floor]
    else:
        return np.c_[X, bedrooms_per_floors]

attr_adder = FunctionTransformer(add_extra_features, validate=False,
                                 kw_args={"add_bathrooms_per_floors": True})
housesales_extra_attribs = attr_adder.fit_transform(housesales.values)

# Query the current dataset's shape
housesales_extra_attribs.shape

# Create a dataframe using our dataset with extra attributes
housesales_extra_attribs = pd.DataFrame(
    housesales_extra_attribs,
    columns=list(housesales.columns)+["bedrooms_per_floor", "bathrooms_per_floor"])
housesales_extra_attribs.head()

# get the right column indices: safer than hard-coding indices
bedrooms_ix, bathrooms_ix, floors_ix = [
    list(housesales_num.columns).index(col)
    for col in ("bedrooms", "bathrooms", "floors")]

# Create and apply a pipeline
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

num_pipeline = Pipeline([
        ('imputer', SimpleImputer(strategy="median")),
        ('attribs_adder', FunctionTransformer(add_extra_features, validate=False)),
        ('std_scaler', StandardScaler()),
    ])

housesales_num_tr = num_pipeline.fit_transform(housesales_num)
housesales_num_tr

# Form a full pipeline and apply it
from sklearn.compose import ColumnTransformer

num_attribs = list(housesales_num)
cat_attribs = ["new_cat"]

full_pipeline = ColumnTransformer([
        ("num", num_pipeline, num_attribs),
        ("cat", OneHotEncoder(), cat_attribs),
    ])

housesales_prepared = full_pipeline.fit_transform(housesales)
housesales_prepared

# Query the resulting shape
housesales_prepared.shape

# Select and fit the model
from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(housesales_prepared, housesales_labels)

# Apply the full preprocessing pipeline to a few training instances and make predictions
some_data = housesales.iloc[:5]
some_labels = housesales_labels.iloc[:5]
some_data_prepared = full_pipeline.transform(some_data)

print("Predictions:", lin_reg.predict(some_data_prepared))

# Compare against the actual values:
print("Labels:", list(some_labels))

# Calculate the root mean square error
from sklearn.metrics import mean_squared_error

housesales_predictions = lin_reg.predict(housesales_prepared)
lin_mse = mean_squared_error(housesales_labels, housesales_predictions)
lin_rmse = np.sqrt(lin_mse)
lin_rmse

# Calculate the mean absolute error
from sklearn.metrics import mean_absolute_error

lin_mae = mean_absolute_error(housesales_labels, housesales_predictions)
lin_mae

# Pick up and fit the model
from sklearn.tree import DecisionTreeRegressor

tree_reg = DecisionTreeRegressor(random_state=10)
tree_reg.fit(housesales_prepared, housesales_labels)

# Calculate the root mean square error
housesales_predictions = tree_reg.predict(housesales_prepared)
tree_mse = mean_squared_error(housesales_labels, housesales_predictions)
tree_rmse = np.sqrt(tree_mse)
tree_rmse

# Pick up and fit the model
from sklearn.ensemble import RandomForestRegressor

forest_reg = RandomForestRegressor(n_estimators=10, random_state=10)
forest_reg.fit(housesales_prepared, housesales_labels)

# Calculate the root mean square error
housesales_predictions = forest_reg.predict(housesales_prepared)
forest_mse = mean_squared_error(housesales_labels, housesales_predictions)
forest_rmse = np.sqrt(forest_mse)
forest_rmse

# Show the three performance results
round(lin_rmse,2), round(tree_rmse,2), round(forest_rmse,2)

# Apply cross-validation to tree regression
from sklearn.model_selection import cross_val_score

scores = cross_val_score(tree_reg, housesales_prepared, housesales_labels,
                         scoring="neg_mean_squared_error", cv=10)
tree_rmse_scores = np.sqrt(-scores)

def display_scores(scores):
    print("Scores:", scores)
    print("Mean:", scores.mean())
    print("Standard deviation:", scores.std())

display_scores(tree_rmse_scores)

# Apply cross-validation to linear regression
lin_scores = cross_val_score(lin_reg, housesales_prepared, housesales_labels,
                             scoring="neg_mean_squared_error", cv=10)
lin_rmse_scores = np.sqrt(-lin_scores)
display_scores(lin_rmse_scores)

# Apply cross-validation Random Forest Regression
from sklearn.model_selection import cross_val_score

forest_scores = cross_val_score(forest_reg, housesales_prepared, housesales_labels,
                                scoring="neg_mean_squared_error", cv=10)
forest_rmse_scores = np.sqrt(-forest_scores)
display_scores(forest_rmse_scores)

# Search for the the best pair of hyperparameters
from sklearn.model_selection import GridSearchCV

param_grid = [
    # try 12 (3×4) combinations of hyperparameters
    {'n_estimators': [3, 10, 30], 'max_features': [4, 8, 12, 16]},
    # then try 6 (2×3) combinations with bootstrap set as False
    {'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [4, 8, 12]},
  ]

forest_reg = RandomForestRegressor(random_state=10)
# train across 5 folds, that's a total of (12+6)*5=90 rounds of training 
grid_search = GridSearchCV(forest_reg, param_grid, cv=5,
                           scoring='neg_mean_squared_error', return_train_score=True)
grid_search.fit(housesales_prepared, housesales_labels)

# Print the best parameters obtained
grid_search.best_params_

# Show the best estimator obtained
grid_search.best_estimator_

# Print all results obtained
cvres = grid_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
    print(np.sqrt(-mean_score), params)

pd.DataFrame(grid_search.cv_results_)

# Search for best hyperparameters using ranges
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint

param_distribs = {
        'n_estimators': randint(low=1, high=200),
        'max_features': randint(low=1, high=16),
    }

forest_reg = RandomForestRegressor(random_state=10)
rnd_search = RandomizedSearchCV(forest_reg, param_distributions=param_distribs,
                                n_iter=10, cv=5, scoring='neg_mean_squared_error', random_state=10)
rnd_search.fit(housesales_prepared, housesales_labels)

# Show results obtained
cvres = rnd_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
    print(np.sqrt(-mean_score), params)

# Show importance of features
feature_importances = grid_search.best_estimator_.feature_importances_
feature_importances

# Show features from most importance to least importance
extra_attribs = ["bedrooms_per_floor", "bathrooms_per_floor"]
cat_encoder = full_pipeline.named_transformers_["cat"]
cat_one_hot_attribs = list(cat_encoder.categories_[0])
attributes = num_attribs + extra_attribs + cat_one_hot_attribs
sorted(zip(feature_importances, attributes), reverse=True)

# Take the best estimator and calculate the RMSE on test data
final_model = grid_search.best_estimator_

X_test = test_set.drop("price", axis=1)
y_test = test_set["price"].copy()

X_test_prepared = full_pipeline.transform(X_test)
final_predictions = final_model.predict(X_test_prepared)

final_mse = mean_squared_error(y_test, final_predictions)
final_rmse = np.sqrt(final_mse)
final_rmse