Step 3: Clustering + Model Exploration

import os
import gc
import datetime
import cloudpickle

import scipy
import numpy as np

import seaborn as sns
import plotly

import dask.array as da
import dask.dataframe as dd
from dask.distributed import Client
from joblib import parallel_backend

from umap import UMAP
from sklearn.decomposition import PCA
from dask_ml.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

from sklearn.neighbors import DistanceMetric
from sklearn.cluster import KMeans, SpectralClustering, DBSCAN

from sklearn.svm import SVC
from rgf.sklearn import RGFClassifier
from dask_ml.xgboost import XGBClassifier
from sklearn.ensemble import AdaBoostClassifier

from sklearn.metrics import (
    confusion_matrix,
    homogeneity_completeness_v_measure,
    precision_recall_fscore_support,
)

# Display all cell outputs
from IPython.core.interactiveshell import InteractiveShell

InteractiveShell.ast_node_interactivity = "all"

# start a local Dask client
# use the daskboard at http://localhost:8787/status
# NOTE: ProgressBar does not show up when client is initialized
client = Client()

# inferred int64 types cause a type mismatch (int vs float) error when dask sees a null value
# null values cannot be interpreted as ints
custom_dtypes = {
    "date": "object",
    "serial_number": "object",
    "capacity_bytes": "float32",
    "failure": "float32",
    "smart_1_normalized": "float32",
    "smart_1_raw": "float32",
    "smart_5_normalized": "float32",
    "smart_5_raw": "float32",
    "smart_7_normalized": "float32",
    "smart_7_raw": "float32",
    "smart_10_normalized": "float32",
    "smart_10_raw": "float32",
    "smart_184_normalized": "float32",
    "smart_184_raw": "float32",
    "smart_187_normalized": "float32",
    "smart_187_raw": "float32",
    "smart_188_normalized": "float32",
    "smart_188_raw": "float32",
    "smart_189_normalized": "float32",
    "smart_189_raw": "float32",
    "smart_190_normalized": "float32",
    "smart_190_raw": "float32",
    "smart_193_normalized": "float32",
    "smart_193_raw": "float32",
    "smart_194_normalized": "float32",
    "smart_194_raw": "float32",
    "smart_197_normalized": "float32",
    "smart_197_raw": "float32",
    "smart_198_normalized": "float32",
    "smart_198_raw": "float32",
    "smart_240_normalized": "float32",
    "smart_240_raw": "float32",
    "smart_241_normalized": "float32",
    "smart_241_raw": "float32",
    "smart_242_normalized": "float32",
    "smart_242_raw": "float32",
}

# read all the cleaned seagate data into one dataframe
DATA_DIR = "/home/kachauha/Downloads/data_Q4_2018_clean"
seagate_df = dd.read_csv(os.path.join(DATA_DIR, "*.csv"), dtype=custom_dtypes)

# get the serial numbers for all the failed hard drives
failed_serials = seagate_df[seagate_df["failure"] == 1]["serial_number"]

# collapse time dim, drop cols that are unncessary for train data
seagate_df_grouped = seagate_df.drop(["date", "failure"], axis=1).groupby(
    "serial_number"
)

# characterize a time series by simple stats - mean,std,max,min,size,kurtosis
means = seagate_df_grouped.mean()
means = means.rename(columns={col: "mean_" + col for col in means.columns})

stds = seagate_df_grouped.std(ddof=0)
stds = stds.rename(columns={col: "std_" + col for col in stds.columns})
stds = stds.fillna(0)  # FIXME: std returns nans even for ddof=0

# maxs = seagate_df_grouped.max()
# maxs = maxs.rename(columns={col: 'max_' + col for col in maxs.columns})

# mins = seagate_df_grouped.min()
# mins = mins.rename(columns={col: 'min_' + col for col in mins.columns})

days = seagate_df_grouped.size().to_frame("days")

# join stats into one df
X_df = means.merge(right=stds, left_index=True, right_index=True)
# X_df = X_df.merge(right=maxs, left_index=True, right_index=True)
# X_df = X_df.merge(right=mins, left_index=True, right_index=True)
X_df = X_df.merge(right=days, left_index=True, right_index=True)

# convert to dask array - some algos dont like dataframes
X_arr = X_df.to_dask_array(lengths=True)

# sanity check - ensure no nans
X_df.isna().any().any().compute()

# clean up memory
del means
del stds
del days
gc.collect()

# get ground truth labels
# NOTE: dask-xgboost needs this to be dask obejct
# NOTE: if not converted to dask array, it is of nan shape
Y_arr = da.from_array(X_df.index.isin(failed_serials.values.compute()).compute())
Y_arr.shape

# true if drive had failed. used for indexing into Xtrain dask array
failed_idx = (Y_arr == 1).compute()
failed_idx.sum()

std_scaler = StandardScaler()
# rob_scaler = RobustScaler()

X_arr = std_scaler.fit_transform(X_arr)
# X = rob_scaler.fit_transform(X)

da.to_npy_stack("preprocessed/X_arr", X_arr)
da.to_npy_stack("preprocessed/Y_arr", Y_arr)

# vanilla pca for visualization in 3d
with parallel_backend("dask"):
    pcs = PCA(n_components=3)
    X_pca = pcs.fit_transform(X_arr)

# plot
working_pca = dict(
    mode="markers",
    name="Working",
    type="scatter3d",
    x=X_pca[~failed_idx, 0],
    y=X_pca[~failed_idx, 1],
    z=X_pca[~failed_idx, 2],
    marker=dict(size=2, color="rgb(0, 0, 200)"),
)

failed_pca = dict(
    mode="markers",
    name="Failed",
    type="scatter3d",
    x=X_pca[failed_idx, 0],
    y=X_pca[failed_idx, 1],
    z=X_pca[failed_idx, 2],
    marker=dict(size=2, color="rgb(200, 0, 0)"),
)

layout_pca = dict(
    title="Simple Reduction Euclidean PCA",
    scene=dict(
        xaxis=dict(zeroline=False),
        yaxis=dict(zeroline=False),
        zaxis=dict(zeroline=False),
    ),
)

fig = dict(data=[working_pca, failed_pca], layout=layout_pca)
plotly.offline.plot(
    fig,
    filename="../../../reports/figures/pca/seagate_simple_euclidean.html",
    auto_open=True,
)

#################################################################
#   KERNEL PCA RUNS OUT OF MEMORY for rbf, sigmoid, cosine      #
#################################################################
# # euclidean is not good enough for the dims we have. go nonlinear
# with parallel_backend('dask'):
#     kpca = KernelPCA(n_components=3, kernel='cosine', random_state=42)
#     X_kpca = kpca.fit_transform(X_arr)

# # plot
# working_kpca = dict(
#     mode = "markers",
#     name = "Working",
#     type = "scatter3d",
#     x = X_kpca[~failed_idx, 0],
#     y = X_kpca[~failed_idx, 1],
#     z = X_kpca[~failed_idx, 2],
#     marker = dict( size=2, color="rgb(0, 0, 200)" )
# )

# failed_kpca = dict(
#     mode = "markers",
#     name = "Failed",
#     type = "scatter3d",
#     x = X_kpca[failed_idx, 0],
#     y = X_kpca[failed_idx, 1],
#     z = X_kpca[failed_idx, 2],
#     marker = dict( size=2, color="rgb(200, 0, 0)" )
# )

# layout_kpca = dict(
#     title = 'Simple Reduction RBF PCA',
#     scene = dict(
#         xaxis = dict( zeroline=False ),
#         yaxis = dict( zeroline=False ),
#         zaxis = dict( zeroline=False ),
#     )
# )

# fig = dict( data=[working_kpca, failed_kpca], layout=layout_kpca )
# plotly.offline.plot(fig, filename='seagate_kpca_rbf.html', auto_open=True)

# visualize using umap (tsne may be too expensive)
with parallel_backend("dask"):
    ucs = UMAP(n_components=3, metric="canberra", n_neighbors=32, random_state=42)
    X_umap = ucs.fit_transform(X_arr, y=da.where(Y_arr == 0, -1, Y_arr))

# plot
working_umap = dict(
    mode="markers",
    name="Working",
    type="scatter3d",
    x=X_umap[~failed_idx, 0],
    y=X_umap[~failed_idx, 1],
    z=X_umap[~failed_idx, 2],
    marker=dict(size=2, color="rgb(0, 0, 200)"),
)

failed_umap = dict(
    mode="markers",
    name="Failed",
    type="scatter3d",
    x=X_umap[failed_idx, 0],
    y=X_umap[failed_idx, 1],
    z=X_umap[failed_idx, 2],
    marker=dict(size=2, color="rgb(200, 0, 0)"),
)

layout_umap = dict(
    title="Simple Reduction Canberra UMAP, Semisupervised",
    scene=dict(
        xaxis=dict(zeroline=False),
        yaxis=dict(zeroline=False),
        zaxis=dict(zeroline=False),
    ),
)

fig = dict(data=[working_umap, failed_umap], layout=layout_umap)
plotly.offline.plot(
    fig,
    filename="../../../reports/figures/umap/seagate_n32_canberra_ss.html",
    auto_open=True,
)

# keep 20 out of 67 dimensions
feat_reducer = UMAP(
    n_components=20, metric="chebyshev", n_neighbors=15, random_state=42
)
X_umap = feat_reducer.fit_transform(X_arr)

def save_model(model, fname, suffix=None):
    """Serialize and save a dask or sklearn model

    Arguments:
        model {dask_ml or sklearn object} -- trained model to save
        fname {str} -- name of file to write to

    Keyword Arguments:
        suffix {str} -- suffix in the filename (default: {None})
                        This serves as the identifier of the current
                        state of the model
                        If not provided, the current timestamp is used
    """
    # generate default suffix if needed
    if suffix is None:
        suffix = datetime.datetime.now().strftime("%b_%d_%Y_%H_%M_%S")

    # serialize and write
    with open(fname + "_" + suffix + ".cpkl", "wb") as file:
        cloudpickle.dump(model, file)

# split the original dataset
with parallel_backend("dask"):
    X_train, X_test, Y_train, Y_test = train_test_split(X_arr, Y_arr, stratify=Y_arr)

# split the reduced dataset
with parallel_backend("dask"):
    X_train_umap, X_test_umap, Y_train_umap, Y_test_umap = train_test_split(
        X_umap, Y_arr, stratify=Y_arr
    )

with parallel_backend("dask"):
    for k in range(2, 10):
        # find kmeans clusters for current k
        km = KMeans(n_clusters=k, random_state=42)
        km.fit(X_arr)

        # get metrics for current k
        h, c, v = homogeneity_completeness_v_measure(Y_train, km.labels_)
        print("========================================")
        print("k={}: H={}, C={}, V={}".format(k, h, c, v))
        print("conf mat:")
        print(confusion_matrix(Y_train, km.labels_))

with parallel_backend("dask"):
    for k in range(2, 6):
        # find kmeans clusters for current k
        sc = SpectralClustering(n_clusters=k, random_state=42)
        sc.fit(X_arr)

        # get metrics for current k
        h, c, v = homogeneity_completeness_v_measure(Y_train, sc.labels_)
        print("========================================")
        print("k={}: H={}, C={}, V={}".format(k, h, c, v))
        print("conf mat:")
        print(confusion_matrix(Y_train, sc.labels_))

man = DistanceMetric.get_metric("manhattan")
X_df.index.head()

foo = X_df.loc["6VDHZAXD", :].compute()
bar = X_df.loc["6VDHZATY", :].compute()

scipy.spatial.distance.cityblock(foo, bar)

with parallel_backend("dask"):
    for epsilon in range(1, 6):
        # find kmeans clusters for current k
        db = DBSCAN(metric="manhattan", eps=epsilon)
        db.fit(X_arr)

        # get metrics for current k
        h, c, v = homogeneity_completeness_v_measure(Y_train, db.labels_)
        print("========================================")
        print("e={}: H={}, C={}, V={}".format(epsilon, h, c, v))
        print("conf mat:")
        print(confusion_matrix(Y_train, db.labels_))

# vanilla classifier
xgbc = XGBClassifier(max_depth=6, random_state=42)
xgbc.fit(X_train, Y_train)

preds = xgbc.predict(X_train)
confusion_matrix(Y_train.compute(), preds.compute())

# reduced data classifier
xgbc_umap = XGBClassifier(max_depth=6, random_state=42)
# xgbc_umap = XGBClassifier(learning_rate=0.1,
#                                max_depth=3,
#                                min_child_weight=0.5,
#                                reg_lambda=0.01,
#                                reg_alpha=0.01,
#                                random_state=42)
xgbc_umap.fit(da.from_array(X_train_umap), Y_train)

preds = xgbc_umap.predict(da.from_array(X_train_umap))
confusion_matrix(Y_train.compute(), preds.compute())

gc.collect()

with parallel_backend("dask"):
    svc = SVC(C=5, class_weight="balanced", probability=True, random_state=42)
    svc.fit(X_train, Y_train)

with parallel_backend("dask"):
    prob_preds_train = svc.predict_proba(X_train)
    prob_preds_test = svc.predict_proba(X_test)

# threshold probability for model train on all data = ~0.3
TRAIN_THRESHOLD_PROBA = 0.5

# how well did it fit on the train data
confmat_train = confusion_matrix(
    Y_train, prob_preds_train[:, 1] > TRAIN_THRESHOLD_PROBA
)
sns.heatmap(
    confmat_train / confmat_train.sum(axis=1, keepdims=True),
    vmin=0,
    vmax=1,
    annot=True,
    cmap="Blues",
)

# how well does it generalize
TEST_THRESHOLD_PROBA = 0.005
confmat_test = confusion_matrix(Y_test, prob_preds_test[:, 1] > TEST_THRESHOLD_PROBA)
sns.heatmap(
    confmat_test / confmat_test.sum(axis=1, keepdims=True),
    vmin=0,
    vmax=1,
    annot=True,
    cmap="Blues",
)

# grid search for best threshold proba
for prob in np.arange(0.33, 0.36, 0.005):
    curr_preds = prob_preds_test[:, 1] > prob
    p, r, f, _ = precision_recall_fscore_support(
        y_true=Y_test, y_pred=curr_preds, average="binary"
    )
    print("p={:2f}, prec={:2f}, rec={:2f}, f={:2f}".format(prob, p, r, f))

with parallel_backend("dask"):
    # init adaboosted ensemble of svm classifiers
    ada_base_learner = SVC(
        C=2.5, class_weight="balanced", probability=True, random_state=42
    )
    ada_clf = AdaBoostClassifier(
        base_estimator=ada_base_learner,
        n_estimators=100,
        learning_rate=0.5,
        random_state=42,
    )
    # training time
    ada_clf.fit(X_train, Y_train)

# how well did it do
with parallel_backend("dask"):
    ada_probs_train = ada_clf.predict_proba(X_train)
    ada_probs_test = ada_clf.predict_proba(X_test)

TRAIN_THRESHOLD_PROBA = 0.5
ada_confmat_train = confusion_matrix(
    Y_train, ada_probs_train[:, 1] > TRAIN_THRESHOLD_PROBA
)
sns.heatmap(
    ada_confmat_train / ada_confmat_train.sum(axis=1, keepdims=True),
    vmin=0,
    vmax=1,
    annot=True,
    cmap="Blues",
)

TEST_THRESHOLD_PROBA = 0.005
ada_confmat_test = confusion_matrix(Y_test, ada_probs_test[:, 1] > TEST_THRESHOLD_PROBA)
sns.heatmap(
    ada_confmat_test / ada_confmat_test.sum(axis=1, keepdims=True),
    vmin=0,
    vmax=1,
    annot=True,
    cmap="Blues",
)

rgfc = RGFClassifier(max_leaf=5000, l2=1e-10, learning_rate=0.01)
rgfc.fit(X_train, Y_train)

%%time
with parallel_backend("dask"):
    rgfc = RGFClassifier(max_leaf=5000, l2=1e-10, learning_rate=0.01)
    rgfc.fit(X_train, Y_train)

with parallel_backend("dask"):
    rgf_prob_preds_train = rgfc.predict_proba(X_train)
    rgf_prob_preds_test = rgfc.predict_proba(X_test)

# threshold probability for model test on all data = ~0.3
TEST_THRESHOLD_PROBA = 0.5

# how well did it fit on the test data
rgf_confmat_test = confusion_matrix(
    Y_test, rgf_prob_preds_test[:, 1] > TEST_THRESHOLD_PROBA
)
sns.heatmap(
    rgf_confmat_test / rgf_confmat_test.sum(axis=1, keepdims=True),
    vmin=0,
    vmax=1,
    annot=True,
    cmap="Blues",
)