import pandas as pd
import numpy
from numpy.random import choice
from sklearn.datasets import load_boston
from h2o.estimators.random_forest import H2ORandomForestEstimator


import h2o
h2o.init()

# transfer the boston data from pandas to H2O
boston_data = load_boston()
X = pd.DataFrame(data=boston_data.data, columns=boston_data.feature_names)
X["Median_value"] = boston_data.target
X = h2o.H2OFrame.from_python(X.to_dict("list"))

# select 10% for valdation
r = X.runif(seed=123456789)
train = X[r < 0.9,:]
valid = X[r >= 0.9,:]

h2o.export_file(train, "Boston_housing_train.csv", force=True)
h2o.export_file(valid, "Boston_housing_test.csv", force=True)

"""
Explanation: H2O Tutorial
Author: Spencer Aiello
Contact: spencer@h2oai.com
This tutorial steps through a quick introduction to H2O's Python API. The goal of this tutorial is to introduce through a complete example H2O's capabilities from Python. Also, to help those that are accustomed to Scikit Learn and Pandas, the demo will be specific call outs for differences between H2O and those packages; this is intended to help anyone that needs to do machine learning on really Big Data make the transition. It is not meant to be a tutorial on machine learning or algorithms.
Detailed documentation about H2O's and the Python API is available at http://docs.h2o.ai.
Setting up your system for this demo
The following code creates two csv files using data from the Boston Housing dataset which is built into scikit-learn and adds them to the local directory
End of explanation
"""


%matplotlib inline
import matplotlib.pyplot as plt

"""
Explanation: Enable inline plotting in the Jupyter Notebook
End of explanation
"""


fr = h2o.import_file("Boston_housing_train.csv")

"""
Explanation: Intro to H2O Data Munging
Read csv data into H2O. This loads the data into the H2O column compressed, in-memory, key-value store.
End of explanation
"""


fr.head()

"""
Explanation: View the top of the H2O frame.
End of explanation
"""


fr.tail()

"""
Explanation: View the bottom of the H2O Frame
End of explanation
"""


fr["CRIM"].head() # Tab completes

"""
Explanation: Select a column
fr["VAR_NAME"]
End of explanation
"""


columns = ["CRIM", "RM", "RAD"]
fr[columns].head()

"""
Explanation: Select a few columns
End of explanation
"""


fr[2:7,:]  # explicitly select all columns with :

"""
Explanation: Select a subset of rows
Unlike in Pandas, columns may be identified by index or column name. Therefore, when subsetting by rows, you must also pass the column selection.
End of explanation
"""


# The columns attribute is exactly like Pandas
print("Columns:", fr.columns, "\n")
print("Columns:", fr.names, "\n")
print("Columns:", fr.col_names, "\n")

# There are a number of attributes to get at the shape
print("length:", str( len(fr) ), "\n")
print("shape:", fr.shape, "\n")
print("dim:", fr.dim, "\n")
print("nrow:", fr.nrow, "\n")
print("ncol:", fr.ncol, "\n")

# Use the "types" attribute to list the column types
print("types:", fr.types, "\n")

"""
Explanation: Key attributes:
      * columns, names, col_names
      * len, shape, dim, nrow, ncol
      * types
Note: 
Since the data is not in local python memory
there is no "values" attribute. If you want to 
pull all of the data into the local python memory
then do so explicitly with h2o.export_file and
reading the data into python memory from disk.
End of explanation
"""


fr.shape

"""
Explanation: Select rows based on value
End of explanation
"""


mask = fr["CRIM"]>1
fr[mask,:].shape

"""
Explanation: Boolean masks can be used to subselect rows based on a criteria.
End of explanation
"""


fr.describe()

"""
Explanation: Get summary statistics of the data and additional data distribution information.
End of explanation
"""


x = fr.names[:]
y="Median_value"
x.remove(y)

"""
Explanation: Set up the predictor and response column names
Using H2O algorithms, it's easier to reference predictor and response columns
by name in a single frame (i.e., don't split up X and y)
End of explanation
"""


# Define and fit first 400 points
model = H2ORandomForestEstimator(seed=42)
model.train(x=x, y=y, training_frame=fr[:400,:])

model.predict(fr[400:fr.nrow,:])        # Predict the rest

"""
Explanation: Machine Learning With H2O
H2O is a machine learning library built in Java with interfaces in Python, R, Scala, and Javascript. It is open source and well-documented.
Unlike Scikit-learn, H2O allows for categorical and missing data.
The basic work flow is as follows:
* Fit the training data with a machine learning algorithm
* Predict on the testing data
Simple model
End of explanation
"""


perf = model.model_performance(fr[400:fr.nrow,:])
perf.r2()      # get the r2 on the holdout data
perf.mse()     # get the mse on the holdout data
perf           # display the performance object

"""
Explanation: The performance of the model can be checked using the holdout dataset
End of explanation
"""


r = fr.runif(seed=12345)   # build random uniform column over [0,1]
train= fr[r<0.75,:]     # perform a 75-25 split
test = fr[r>=0.75,:]

model = H2ORandomForestEstimator(seed=42)
model.train(x=x, y=y, training_frame=train, validation_frame=test)

perf = model.model_performance(test)
perf.r2()

"""
Explanation: Train-Test Split
Instead of taking the first 400 observations for training, we can use H2O to create a random test train split of the data.
End of explanation
"""


model = H2ORandomForestEstimator(nfolds=10) # build a 10-fold cross-validated model
model.train(x=x, y=y, training_frame=fr)

scores = numpy.array([m.r2() for m in model.xvals]) # iterate over the xval models using the xvals attribute
print("Expected R^2: %.2f +/- %.2f \n" % (scores.mean(), scores.std()*1.96))
print("Scores:", scores.round(2))

"""
Explanation: There was a massive jump in the R^2 value. This is because the original data is not shuffled.
Cross validation
H2O's machine learning algorithms take an optional parameter nfolds to specify the number of cross-validation folds to build. H2O's cross-validation uses an internal weight vector to build the folds in an efficient manner (instead of physically building the splits).
In conjunction with the nfolds parameter, a user may specify the way in which observations are assigned to each fold with the fold_assignment parameter, which can be set to either:
        * AUTO:  Perform random assignment
        * Random: Each row has a equal (1/nfolds) chance of being in any fold.
        * Modulo: Observations are in/out of the fold based by modding on nfolds
End of explanation
"""


from sklearn.cross_validation import cross_val_score
from h2o.cross_validation import H2OKFold
from h2o.model.regression import h2o_r2_score
from sklearn.metrics.scorer import make_scorer

"""
Explanation: However, you can still make use of the cross_val_score from Scikit-Learn
Cross validation: H2O and Scikit-Learn
End of explanation
"""


model = H2ORandomForestEstimator(seed=42)

scorer = make_scorer(h2o_r2_score)   # make h2o_r2_score into a scikit_learn scorer
custom_cv = H2OKFold(fr, n_folds=10, seed=42) # make a cv 
scores = cross_val_score(model, fr[x], fr[y], scoring=scorer, cv=custom_cv)

print("Expected R^2: %.2f +/- %.2f \n" % (scores.mean(), scores.std()*1.96))
print("Scores:", scores.round(2))

"""
Explanation: You still must use H2O to make the folds. Currently, there is no H2OStratifiedKFold. Additionally, the H2ORandomForestEstimator is similar to the scikit-learn RandomForestRegressor object with its own train method.
End of explanation
"""


h2o.__PROGRESS_BAR__=False
h2o.no_progress()

"""
Explanation: There isn't much difference in the R^2 value since the fold strategy is exactly the same. However, there was a major difference in terms of computation time and memory usage.
Since the progress bar print out gets annoying let's disable that
End of explanation
"""


from sklearn import __version__
sklearn_version = __version__
print(sklearn_version)

"""
Explanation: Grid Search
Grid search in H2O is still under active development and it will be available very soon. However, it is possible to make use of Scikit's grid search infrastructure (with some performance penalties)
Randomized grid search: H2O and Scikit-Learn
End of explanation
"""


%%time
from sklearn.grid_search import RandomizedSearchCV  # Import grid search
from scipy.stats import randint, uniform

model = H2ORandomForestEstimator(seed=42)        # Define model

params = {"ntrees": randint(20,30),
          "max_depth": randint(1,10),
          "min_rows": randint(1,10),    # scikit's  min_samples_leaf
          "mtries": randint(2,fr[x].shape[1]),} # Specify parameters to test

scorer = make_scorer(h2o_r2_score)   # make h2o_r2_score into a scikit_learn scorer
custom_cv = H2OKFold(fr, n_folds=5, seed=42) # make a cv 
random_search = RandomizedSearchCV(model, params, 
                                   n_iter=10, 
                                   scoring=scorer, 
                                   cv=custom_cv, 
                                   random_state=42,
                                   n_jobs=1)       # Define grid search object

random_search.fit(fr[x], fr[y])

print("Best R^2:", random_search.best_score_, "\n")
print("Best params:", random_search.best_params_)

"""
Explanation: If you have 0.16.1, then your system can't handle complex randomized grid searches (it works in every other version of sklearn, including the soon to be released 0.16.2 and the older versions).
The steps to perform a randomized grid search:
1. Import model and RandomizedSearchCV
2. Define model
3. Specify parameters to test
4. Define grid search object
5. Fit data to grid search object
6. Collect scores
All the steps will be repeated from above.
Because 0.16.1 is installed, we use scipy to define specific distributions
ADVANCED TIP:
Turn off reference counting for spawning jobs in parallel (n_jobs=-1, or n_jobs > 1).
We'll turn it back on again in the aftermath of a Parallel job.
If you don't want to run jobs in parallel, don't turn off the reference counting.
Pattern is:
         >>> h2o.turn_off_ref_cnts()
         >>> .... parallel job ....
         >>> h2o.turn_on_ref_cnts()
End of explanation
"""


def report_grid_score_detail(random_search, charts=True):
    """Input fit grid search estimator. Returns df of scores with details"""
    df_list = []

    for line in random_search.grid_scores_:
        results_dict = dict(line.parameters)
        results_dict["score"] = line.mean_validation_score
        results_dict["std"] = line.cv_validation_scores.std()*1.96
        df_list.append(results_dict)

    result_df = pd.DataFrame(df_list)
    result_df = result_df.sort("score", ascending=False)
    
    if charts:
        for col in get_numeric(result_df):
            if col not in ["score", "std"]:
                plt.scatter(result_df[col], result_df.score)
                plt.title(col)
                plt.show()

        for col in list(result_df.columns[result_df.dtypes == "object"]):
            cat_plot = result_df.score.groupby(result_df[col]).mean()[0]
            cat_plot.sort()
            cat_plot.plot(kind="barh", xlim=(.5, None), figsize=(7, cat_plot.shape[0]/2))
            plt.show()
    return result_df

def get_numeric(X):
    """Return list of numeric dtypes variables"""
    return X.dtypes[X.dtypes.apply(lambda x: str(x).startswith(("float", "int", "bool")))].index.tolist()

report_grid_score_detail(random_search).head()

"""
Explanation: We might be tempted to think that we just had a large improvement; however we must be cautious. The function below creates a more detailed report.
End of explanation
"""


%%time

params = {"ntrees": randint(30,35),
          "max_depth": randint(5,8),
          "mtries": randint(4,6),}

custom_cv = H2OKFold(fr, n_folds=5, seed=42)           # In small datasets, the fold size can have a big
                                                       # impact on the std of the resulting scores. More
random_search = RandomizedSearchCV(model, params,      # folds --> Less examples per fold --> higher 
                                   n_iter=5,           # variation per sample
                                   scoring=scorer, 
                                   cv=custom_cv, 
                                   random_state=43, 
                                   n_jobs=1)       

random_search.fit(fr[x], fr[y])

print("Best R^2:", random_search.best_score_, "\n")
print("Best params:", random_search.best_params_)

report_grid_score_detail(random_search)

"""
Explanation: Based on the grid search report, we can narrow the parameters to search and rerun the analysis. The parameters below were chosen after a few runs:
End of explanation
"""


from h2o.transforms.preprocessing import H2OScaler
from h2o.estimators.pca import H2OPrincipalComponentAnalysisEstimator as H2OPCA

"""
Explanation: Transformations
Rule of machine learning: Don't use your testing data to inform your training data. Unfortunately, this happens all the time when preparing a dataset for the final model. But on smaller datasets, you must be especially careful.
At the moment, there are no classes for managing data transformations. On the one hand, this requires the user to tote around some extra state, but on the other, it allows the user to be more explicit about transforming H2OFrames.
Basic steps:

Remove the response variable from transformations.
Import transformer
Define transformer
Fit train data to transformer
Transform test and train data
Re-attach the response variable.

First let's normalize the data using the means and standard deviations of the training data.
Then let's perform a principal component analysis on the training data and select the top 5 components.
Using these components, let's use them to reduce the train and test design matrices.
End of explanation
"""


y_train = train.pop("Median_value")
y_test  = test.pop("Median_value")

norm = H2OScaler()
norm.fit(train)
X_train_norm = norm.transform(train)
X_test_norm  = norm.transform(test)

print(X_test_norm.shape)
X_test_norm

"""
Explanation: Normalize Data: Use the means and standard deviations from the training data.
End of explanation
"""


pca = H2OPCA(k=5)
pca.fit(X_train_norm)
X_train_norm_pca = pca.transform(X_train_norm)
X_test_norm_pca  = pca.transform(X_test_norm)

# prop of variance explained by top 5 components?

print(X_test_norm_pca.shape)
X_test_norm_pca[:5]

model = H2ORandomForestEstimator(seed=42)
model.train(x=X_train_norm_pca.names, y=y_train.names, training_frame=X_train_norm_pca.cbind(y_train))
y_hat  = model.predict(X_test_norm_pca)

h2o_r2_score(y_test,y_hat)

"""
Explanation: Then, we can apply PCA and keep the top 5 components. A user warning is expected here.
End of explanation
"""


from h2o.transforms.preprocessing import H2OScaler
from h2o.estimators.pca import H2OPrincipalComponentAnalysisEstimator as H2OPCA

from sklearn.pipeline import Pipeline                # Import Pipeline <other imports not shown>
model = H2ORandomForestEstimator(seed=42)
pipe = Pipeline([("standardize", H2OScaler()),       # Define pipeline as a series of steps
                 ("pca", H2OPCA(k=5)),
                 ("rf", model)])                     # Notice the last step is an estimator

pipe.fit(train, y_train)                             # Fit training data
y_hat = pipe.predict(test)                           # Predict testing data (due to last step being an estimator)
h2o_r2_score(y_test, y_hat)                          # Notice the final score is identical to before

"""
Explanation: Although this is MUCH simpler than keeping track of all of these transformations manually, it gets to be somewhat of a burden when you want to chain together multiple transformers.
Pipelines
"Tranformers unite!"
If your raw data is a mess and you have to perform several transformations before using it, use a pipeline to keep things simple.
Steps:

Import Pipeline, transformers, and model
Define pipeline. The first and only argument is a list of tuples where the first element of each tuple is a name you give the step and the second element is a defined transformer. The last step is optionally an estimator class (like a RandomForest).
Fit the training data to pipeline
Either transform or predict the testing data
End of explanation
"""


pipe = Pipeline([("standardize", H2OScaler()),
                 ("pca", H2OPCA()),
                 ("rf", H2ORandomForestEstimator(seed=42))])

params = {"standardize__center":    [True, False],           # Parameters to test
          "standardize__scale":     [True, False],
          "pca__k":                 randint(2, 6),
          "rf__ntrees":             randint(10,20),
          "rf__max_depth":          randint(4,10),
          "rf__min_rows":           randint(5,10), }
#           "rf__mtries":             randint(1,4),}           # gridding over mtries is 
                                                               # problematic with pca grid over 
                                                               # k above 

from sklearn.grid_search import RandomizedSearchCV
from h2o.cross_validation import H2OKFold
from h2o.model.regression import h2o_r2_score
from sklearn.metrics.scorer import make_scorer

custom_cv = H2OKFold(fr, n_folds=5, seed=42)
random_search = RandomizedSearchCV(pipe, params,
                                   n_iter=5,
                                   scoring=make_scorer(h2o_r2_score),
                                   cv=custom_cv,
                                   random_state=42,
                                   n_jobs=1)


random_search.fit(fr[x],fr[y])
results = report_grid_score_detail(random_search)
results.head()

"""
Explanation: This is so much easier!!!
But, wait a second, we did worse after applying these transformations! We might wonder how different hyperparameters for the transformations impact the final score.
Combining randomized grid search and pipelines
"Yo dawg, I heard you like models, so I put models in your models to model models."
Steps:

Import Pipeline, grid search, transformers, and estimators <Not shown below>
Define pipeline
Define parameters to test in the form: "(Step name)__(argument name)" A double underscore separates the two words.
Define grid search
Fit to grid search
End of explanation
"""


best_estimator = random_search.best_estimator_                        # fetch the pipeline from the grid search
h2o_model      = h2o.get_model(best_estimator._final_estimator._id)    # fetch the model from the pipeline

save_path = h2o.save_model(h2o_model, path=".", force=True)
print(save_path)

# assumes new session
my_model = h2o.load_model(path=save_path)

my_model.predict(X_test_norm_pca)

"""
Explanation: Currently Under Development (drop-in scikit-learn pieces):
    * Richer set of transforms (only PCA and Scale are implemented)
    * Richer set of estimators (only RandomForest is available)
    * Full H2O Grid Search
Other Tips: Model Save/Load
It is useful to save constructed models to disk and reload them between H2O sessions. Here's how:
End of explanation
"""


