# IMPORTS

# Staple Inputs
import numpy as np
import pandas as pd
from pandas import DataFrame

# Visualization
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
%matplotlib inline

# Dataset
from sklearn.datasets import load_boston

# Loading the Boston dataset
data_boston = load_boston()

print(data_boston.DESCR)

boston_df = DataFrame(data_boston.data)
boston_df.columns = data_boston.feature_names

# Complete df
boston_df['PRICE'] = data_boston.target

# Separating the df into their respective df
X_fts = boston_df.drop('PRICE',axis=1)
y_price =  boston_df['PRICE']

boston_df.head()

boston_df.describe()

boston_df.dtypes

# Using the apply to function to get the total amount of unique observations
boston_df.apply(lambda x: len(x.unique()))

# Using the apply function, for every column, we find the total amount of NULL/NA values
boston_df.apply(lambda x: sum(x.isnull()))

## Indepth look at price

min_price = np.min(y_price)
max_price = np.max(y_price)
mean_price = np.mean(y_price)
median_price = np.median(y_price)
std_price = np.std(y_price)

print("Descriptive Statistics\n")
print("The maximum price for the PRICE is {0:.2f}".format(max_price*1000))
print("The minimum price for the PRICE is {0:.2f}".format(min_price*1000))
print("The mean price for the PRICE is {0:.2f}".format(mean_price*1000))
print("The median price for the PRICE is {0:.2f}".format(median_price*1000))
print("The std fo the price for the PRICE is {0:.2f}".format(std_price*1000))

all_cols = boston_df.columns

# Histogram of prices (this is the target of our dataset)
plt.hist(y_price,bins=60)

plt.xlabel('Price in $1000s')
plt.ylabel('Number of houses')
plt.show()

# creating a visual of the correlation of all the dependent variables
for i, cols in enumerate(all_cols):    
    sns.lmplot(cols, "PRICE",data = boston_df)
    plt.show()

# Compute the correlation matrix
corr = boston_df.corr()


# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(11, 9))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, annot=True, vmax=.3, center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5})

corr_matrix = boston_df.corr()
corr_matrix['PRICE'].sort_values(ascending=False)

import statsmodels.api as sm
from scipy import stats
from collections import defaultdict

global dict_adjus_R
dict_adjus_R = defaultdict(list)

def HighestPvalue(model, threshold):
    highest_pvalue = 0
    
    for index, current_pvalue in model.pvalues.items():
        if current_pvalue > highest_pvalue:
            highest_pvalue = current_pvalue
            highest_index = index
            
    if highest_pvalue > threshold: return highest_index
    else: return True

def CreateLinearReg(x, y):
    X2 = sm.add_constant(x)
    est = sm.OLS(y, X2)
    est2 = est.fit()
    return est2

def BackwardElimination(Xs, y, stats_signf):
    model_info = CreateLinearReg(Xs, y)
    p_results = HighestPvalue(model_info, stats_signf)

    dict_adjus_R[len(Xs.columns)].append([Xs.columns, model_info.rsquared_adj])
    
    if p_results is True: return model_info

    else:
        Xs.drop(p_results, axis=1, inplace=True)
        BackwardElimination(Xs, y, stats_signf)

# Statistical sigficance we would like to uses
stats_signf = 0.05
final_model = BackwardElimination(X_fts, y_price, stats_signf)
final_model.summary()

# Dictionary with some of their best models. 
# This was done to look at a group of variable and there adjusted R_Squared
dict_adjus_R

# Based on the analysis using backward elimination, I will look into using these variables
X_fts_1 = X_fts[['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 
                        'RAD', 'TAX', 'PTRATIO', 'B','LSTAT']]

X_fts_1.head()

from sklearn.metrics import r2_score

def Score_R2(y_true, y_predict):
    """ Calculates and returns the performance score between 
        true and predicted values based on the metric chosen. """

    # TODO: Calculate the performance score between 'y_true' and 'y_predict'
    score = r2_score(y_true, y_predict)

    # Return the score
    return score

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()

from sklearn.model_selection import train_test_split

# Splitting the dataset to understand how the model perform on a simple split
X_train, X_test, y_train, y_test = train_test_split(X_fts_1, y_price, test_size=0.2)
model = linreg.fit(X_train, y_train)
predictions = linreg.predict(X_test)

# Visualizing the performance of the dataset

plt.scatter(y_test, predictions)
plt.xlabel("True Values")
plt.ylabel("Predictions")
plt.show()

model_score = model.score(X_test, y_test)
print ("Model R_Square Performance:, {0:.2f}".format(model_score))

from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score

seed = 7
splits = 10
kfold = KFold(n_splits=splits, random_state=seed)

scoring = 'neg_mean_squared_error'
results = cross_val_score(linreg, X_fts_1, y_price, cv=kfold, scoring=scoring)

all_RMSE = [np.sqrt(np.abs(result)) for result in results]
RMSE = np.sqrt(np.abs(results.mean()))
RMSE_std = np.sqrt(results.std())

print("For the RMSE, the mean is {0:.3f} and the std. deviation is {1:.3f}\n".format(RMSE, RMSE_std))
print("The RMSE for all the CV are {0}".format('\n'.join(str(r) for r in all_RMSE)))

scoring = 'r2'

results = cross_val_score(linreg, X_fts_1, y_price, cv=kfold, scoring=scoring)
print("For the R Squared, the mean is {0:.3f} and the std. deviation is {1:.3f}\n".format(results.mean(), results.std()))
print("The R Squared for all the CV are \n{0}".format('\n'.join(str(r) for r in results)))

import visuals as vs

# Produce learning curves for varying training set sizes and maximum depths
vs.ModelLearning(X_fts_1, y_price)

vs.ModelComplexity(X_train, y_train)

from sklearn.metrics import make_scorer
from sklearn.model_selection import GridSearchCV
from sklearn.tree import DecisionTreeRegressor

# Using the variables from the grid search
def GridSearch_DecReg(X, y):
    """ Performs grid search over the 'max_depth' parameter for a 
        decision tree regressor trained on the input data [X, y]. """
    
    # Create cross-validation sets from the training data
    cv_sets = ShuffleSplit(X.shape[0], test_size = 0.20, random_state = 0)

    # TODO: Create a decision tree regressor object
    regressor = DecisionTreeRegressor()

    # TODO: Create a dictionary for the parameter 'max_depth' with a range from 1 to 10
    params = {'max_depth': range(1,11)}

    # TODO: Transform 'performance_metric' into a scoring function using 'make_scorer' 
    scoring_fnc = make_scorer(performance_metric)

    # TODO: Create the grid search object
    grid = GridSearchCV(regressor,params,scoring_fnc,cv=cv_sets)

    # Fit the grid search object to the data to compute the optimal model
    grid = grid.fit(X, y)

    # Return the optimal model after fitting the data
    return grid

from sklearn.externals.six import StringIO
from IPython.display import Image  
from sklearn.tree import export_graphviz
import pydotplus


def VisualizeDecisionTree(parameters, x, y):
    
    params = parameters.best_estimator_.get_params()

        
    dt_model = DecisionTreeRegressor(**params)
    dt_fit = dt_model.fit(x, y)
    ####
    
    dot_data = StringIO()
    export_graphviz(dt_fit, out_file=dot_data, special_characters=True, 
                     filled=True, rounded=True, feature_names=x.columns)

    graph = pydotplus.graph_from_dot_data(dot_data.getvalue())  
    return Image(graph.create_png())
 

grid_p = GridSearch_DecReg(X_fts_1, y_price)

# Produce the value for 'max_depth'
print("Parameter 'max_depth' is {0} for the optimal model.".format(grid_p.best_estimator_.get_params()['max_depth']))

# Checking the scores
# print(grid_p.cv_results_)

#Creating the decision tree for the specific 
X_train, X_test, y_train, y_test = train_test_split(X_fts_1, y_price, test_size=0.2)

params = grid_p.best_estimator_.get_params()

dt_model = DecisionTreeRegressor(**params)
dt_fit = dt_model.fit(X_train, y_train)

dt_scores = cross_val_score(dt_fit, X_train, y_train, cv = 10)

r2_sqr_ytest = r2_score(y_test, grid_p.best_estimator_.predict(X_test))
score_ytest = dt_fit.score(X_test, y_test)

print("""R Squared using the predicted model using the gridCV parameters
      for the y test is {0:.2f}""".format(r2_sqr_ytest))

print("The score fitting for the testing set is {0:.2f}".format(score_ytest))
print("Mean cross validation score: {0:.2f}".format(np.mean(dt_scores)))

VisualizeDecisionTree(grid_p, X_train, y_train)

# New set of variables
X_fts_2 = X_fts_1[['LSTAT','RM','DIS','CRIM','DIS']]

grid_p_2 = GridSearch_DecReg(X_fts_2, y_price)

# Produce the value for 'max_depth'
print("Parameter 'max_depth' is {0} for the optimal model.".format(grid_p_2.best_estimator_.get_params()['max_depth']))

#Creating the decision tree for the specific 
X_train, X_test, y_train, y_test = train_test_split(X_fts_2, y_price, test_size=0.2)

params = grid_p_2.best_estimator_.get_params()

dt_model = DecisionTreeRegressor(**params)
dt_fit = dt_model.fit(X_train, y_train)

dt_scores = cross_val_score(dt_fit, X_train, y_train, cv = 10)

r2_sqr_ytest = r2_score(y_test, grid_p_2.best_estimator_.predict(X_test))
score_ytest = dt_fit.score(X_test, y_test)

print("""R Squared using the predicted model using the gridCV parameters
      for the y test is {0:.2f}""".format(r2_sqr_ytest))

print("The score fitting for the testing set is {0:.2f}".format(score_ytest))
print("Mean cross validation score: {0:.2f}".format(np.mean(dt_scores)))

# Will redo the process with the top 5 variables!
vs.ModelLearning(X_fts_2, y_price)

VisualizeDecisionTree(grid_p, X_train, y_train)

from collections import defaultdict

# Let perform for all of them with differernt ranges
ranges = [0.025, 0.05, 0.1]
var_d = defaultdict(list)

from collections import defaultdict

client_data = {}

for val in boston_df.columns:
    min_val = np.min(boston_df[val])
    max_val = np.max(boston_df[val])
    sampl = np.random.uniform(low=min_val, high=max_val, size=(10,))
    client_data["{0}".format(val)] = sampl

from pandas import DataFrame
client_df = DataFrame.from_dict(client_data)
client_df = client_df[['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 
           'RAD', 'TAX', 'PTRATIO', 'B','LSTAT']]

for i, price in enumerate(grid_p.best_estimator_.predict(client_df)):
    print ("Predicted selling price for Client {0}'s home: ${1:,.2f}".format(i+1, price*1000))

import matplotlib.pyplot as plt
plt.hist(y_price, bins = 20)
for price in grid_p.best_estimator_.predict(client_df):
    plt.axvline(price, lw = 5, c = 'r')

# Using the variables from the grid search
def Temp_GridSearch(X, y):
    cv_sets = ShuffleSplit(X.shape[0], test_size = 0.20, random_state = 0)
    regressor = DecisionTreeRegressor()
    params = {'max_depth': range(1,11)}
    scoring_fnc = make_scorer(performance_metric)
    grid = GridSearchCV(regressor,params,scoring_fnc,cv=cv_sets)
    grid = grid.fit(X, y)
    return grid.best_estimator_

vs.PredictTrials(X_fts_1, y_price, Temp_GridSearch, client_df.as_matrix())