import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cPickle as pickle
%matplotlib notebook
plt.style.use('ggplot')

from sklearn.preprocessing import StandardScaler
from sklearn.cross_validation import train_test_split, KFold
from sklearn.metrics import confusion_matrix, roc_auc_score, precision_score, recall_score, classification_report
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.feature_selection import f_classif
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.ensemble.partial_dependence import plot_partial_dependence, partial_dependence
from sklearn.grid_search import GridSearchCV

"""
Explanation: First iteration of modeling, going through several different classification algorithms. Of these, Gradient Boosted Classifier using unscaled data worked the best.
For reference/example only.
End of explanation
"""


masterdf_rrt = pd.read_csv('RRT_modeling_table_13hr_raw.csv')
masterdf_nonrrt = pd.read_csv('NonRRT_modeling_table_13hr_raw.csv')

col_use = ['on_iv', 'bu-nal',
           'age', 'sex', 'obese', 'smoker', 'prev_rrt',
           'DBP_mean', 'DBP_recent', # take the mean of all the measurements & the most recently observed point
            'SBP_mean', 'SBP_recent',
            'MAP_mean', 'MAP_recent', # mean arterial pressure
             'temp_mean', 'temp_recent',# temperature
             'SPO2_mean', 'SPO2_recent',
            'RR_mean', 'RR_recent', # respiratory rate
            'pulse_mean', 'pulse_recent',
           'anticoagulants', 'narcotics', 'narc-ans', #narcotic analgesics
            'antipsychotics', 'chemo'
]

colsfornan = ['DBP_mean', 'DBP_recent', # take the mean of all the measurements & the most recently observed point
            'SBP_mean', 'SBP_recent',
            'MAP_mean', 'MAP_recent', # mean arterial pressure
             'temp_mean', 'temp_recent',# temperature
             'SPO2_mean', 'SPO2_recent',
            'RR_mean', 'RR_recent', # respiratory rate
            'pulse_mean', 'pulse_recent']


# take out: rrt_ce_id, encntr_id, event_end_dt_tm, timestart, timeend,
# RASS score, GCS score, HR, CO2, O2
X_rrt = masterdf_rrt[col_use]

# if 'obese' is Nan, then set the patient to be not obese.
obesenanmask = np.where(pd.isnull(X_rrt['obese']))[0]
X_rrt.loc[obesenanmask, 'obese'] = 0

# Write out rows that are not all 0/NaNs across. (if all nans, remove this sample)
X_rrt = X_rrt.loc[np.where(X_rrt.ix[:, colsfornan].sum(axis=1, skipna=True)!=0)[0]]

#reset index
X_rrt = X_rrt.reset_index(drop=True)

# take out: encntr_id, not_rrt_time, timestart, timeend,
#  RASS score, GCS score, HR, CO2, O2
X_notrrt = masterdf_nonrrt[col_use]

# if 'obese' is Nan, then set the patient to be not obese.
obesenanmask = np.where(pd.isnull(X_notrrt['obese']))[0]
X_notrrt.loc[obesenanmask, 'obese'] = 0

# Write out rows that are NOT all 0/NaNs across.
X_notrrt = X_notrrt.iloc[np.where(X_notrrt.ix[:, colsfornan].sum(axis=1, skipna=True)!=0)[0]]

#reset index
X_notrrt = X_notrrt.reset_index(drop=True)


# make sure to reset index if haven't previously -- I did not run this before saving processed .p files
X_notrrt = X_notrrt.reset_index(drop=True)
X_rrt = X_rrt.reset_index(drop=True)


# additional for non-rrt: drop samples with lots of NaNs since we have plenty of samples
# DROP THE ROWS WHERE PULSE IS NAN
X_notrrt = X_notrrt.ix[np.where(pd.isnull(X_notrrt['pulse_mean'])!=True)[0]]
X_notrrt = X_notrrt.reset_index(drop=True)

# Do a similar thing for all rows with significant nans:
X_notrrt = X_notrrt.ix[np.where(pd.isnull(X_notrrt['RR_mean'])!=True)[0]]
X_notrrt = X_notrrt.reset_index(drop=True)
X_notrrt = X_notrrt.ix[np.where(pd.isnull(X_notrrt['MAP_mean'])!=True)[0]]
X_notrrt = X_notrrt.reset_index(drop=True)
X_notrrt = X_notrrt.ix[np.where(pd.isnull(X_notrrt['temp_mean'])!=True)[0]]
X_notrrt = X_notrrt.reset_index(drop=True)
X_notrrt = X_notrrt.ix[np.where(pd.isnull(X_notrrt['SPO2_mean'])!=True)[0]]
X_notrrt = X_notrrt.reset_index(drop=True)

# add labels to indicate positive or negative class
X_rrt['label'] = 1
X_notrrt['label'] = 0

# Combine the tables
XY = pd.concat([X_rrt, X_notrrt])
XY = XY.reset_index(drop=True)
y = XY.pop('label')
X = XY

# Fill nans with mean of columns
X = X.fillna(X.mean())

# map genders to 1/0
X['is_male'] = X['sex'].map({'M': 1, 'F': 0})
X.pop('sex')

X.describe().T

y.describe().T

"""
Explanation: Read in data
We did not share our modeling data, so you will have to create your own. The pipeline tool can help you do this. If you save the results to a csv, masterdf_rrt and masterdf_nonrrt are dataframes with the modeling data for each of the positive and negative classes, respectively.
End of explanation
"""


X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

len(X_train)

def score_printout(X_test, y_test, fittedModel):
    print "AUC-ROC Score of model: ", roc_auc_score(y_test, fittedModel.predict_proba(X_test)[:,1])
    print "Precision Score of model: ", precision_score(y_test, fittedModel.predict(X_test))
    print "Recall Score of model: ", recall_score(y_test, fittedModel.predict(X_test))

"""
Explanation: Modeling portion
End of explanation
"""


lrCV = LogisticRegressionCV(Cs=10, class_weight=None, cv=5, dual=False,
           fit_intercept=True, intercept_scaling=1.0, max_iter=100,
           multi_class='ovr', n_jobs=1, penalty='l2', random_state=1,
           refit=True, scoring='roc_auc', solver='lbfgs', tol=0.0001, verbose=0)
lrCV.fit(X, y)
lrCV.scores_

# Try different solver
lrCV = LogisticRegressionCV(Cs=5, class_weight=None, cv=5, dual=False,
           fit_intercept=True, intercept_scaling=1.0, max_iter=100,
           multi_class='ovr', n_jobs=1, penalty='l2', random_state=1,
           refit=True, scoring='roc_auc', solver='liblinear', tol=0.0001, verbose=0)
lrCV.fit(X_train, y_train)
lrCV.scores_

lrCV.intercept_

lrCV.n_iter_

lrCV.Cs_

score_printout(X_test, y_test, lrCV)

confusion_matrix(y_test, lrCV.predict(X_test))

"""
Explanation: Logistic Regression
End of explanation
"""


lr = LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=-1,
          penalty='l2', random_state=1, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)
lr.fit(X_train, y_train)

score_printout(X_test, y_test, lr)

print classification_report(y_test, lr.predict(X_test))

confusion_matrix(y_test, lr.predict(X_test))

lr.decision_function(X_test)

lr.coef_

lr.intercept_

"""
Explanation: compare to traditional LR...
End of explanation
"""


Xscaled = StandardScaler().fit_transform(X)
Xs_train, Xs_test, ys_train, ys_test = train_test_split(Xscaled, y, test_size=0.3)

lrs = LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=-1,
          penalty='l2', random_state=1, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)
lrs.fit(Xs_train, ys_train)

score_printout(Xs_test, ys_test, lrs)
print classification_report(ys_test, lrs.predict(Xs_test))
confusion_matrix(ys_test, lrs.predict(Xs_test))

"""
Explanation: LR & LRCV returned essentially the same results!
Let's rerun LR with scaled data. And then calculate
End of explanation
"""


rfs = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
            oob_score=False, random_state=None, verbose=0,
            warm_start=False)
rfs.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, rfs)
print classification_report(ys_test, rfs.predict(Xs_test))
confusion_matrix(ys_test, rfs.predict(Xs_test))

rf = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
            oob_score=False, random_state=None, verbose=0,
            warm_start=False)
rf.fit(X_train, y_train)
score_printout(X_test, y_test, rf)
print classification_report(y_test, rf.predict(X_test))
confusion_matrix(y_test, rf.predict(X_test))
# scaled & unscaled random forest looks very similar.

rfs = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='entropy',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=-1,
            oob_score=False, random_state=None, verbose=0,
            warm_start=False)
rfs.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, rfs)
print classification_report(ys_test, rfs.predict(Xs_test))
confusion_matrix(ys_test, rfs.predict(Xs_test))

# Increase # estimators
# Note, turning oob score on makes the result worse...
rfs = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='entropy',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, n_estimators=100, n_jobs=-1,
            oob_score=False, random_state=None, verbose=0,
            warm_start=False)
rfs.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, rfs)
print classification_report(ys_test, rfs.predict(Xs_test))
confusion_matrix(ys_test, rfs.predict(Xs_test))

"""
Explanation: Scaled version performs (slightly) better
Random Forest
End of explanation
"""


svms =  SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=3, gamma='auto', kernel='rbf',
  max_iter=-1, probability=True, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
svms.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, svms)
print classification_report(ys_test, svms.predict(Xs_test))
confusion_matrix(ys_test, svms.predict(Xs_test))

svms =  SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=3, gamma='auto', kernel='poly',
  max_iter=-1, probability=True, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
svms.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, svms)
print classification_report(ys_test, svms.predict(Xs_test))
confusion_matrix(ys_test, svms.predict(Xs_test))

svms =  SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=3, gamma='auto', kernel='linear',
  max_iter=-1, probability=True, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
svms.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, svms)
print classification_report(ys_test, svms.predict(Xs_test))
confusion_matrix(ys_test, svms.predict(Xs_test))

svms =  SVC(C=100.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=3, gamma='auto', kernel='rbf',
  max_iter=-1, probability=True, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
svms.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, svms)
print classification_report(ys_test, svms.predict(Xs_test))
confusion_matrix(ys_test, svms.predict(Xs_test))

"""
Explanation: SVM (kernel: sigmoid does not work, rbf appears best)
End of explanation
"""


ldas = LinearDiscriminantAnalysis(n_components=None, priors=None, shrinkage=None,
              solver='svd', store_covariance=False, tol=0.0001)
ldas.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, ldas)
print classification_report(ys_test, ldas.predict(Xs_test))
confusion_matrix(ys_test, ldas.predict(Xs_test))

ldas = LinearDiscriminantAnalysis(n_components=None, priors=None, shrinkage='auto',
              solver='eigen', store_covariance=False, tol=0.0001)
ldas.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, ldas)
print classification_report(ys_test, ldas.predict(Xs_test))
confusion_matrix(ys_test, ldas.predict(Xs_test))

"""
Explanation: LDA
End of explanation
"""


gbcs = GradientBoostingClassifier(init=None, learning_rate=0.1, loss='deviance',
              max_depth=3, max_features=None, max_leaf_nodes=None,
              min_samples_leaf=1, min_samples_split=2,
              min_weight_fraction_leaf=0.0, n_estimators=100,
              presort='auto', random_state=1, subsample=1.0, verbose=0,
              warm_start=False)
gbcs.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, gbcs)
print classification_report(ys_test, gbcs.predict(Xs_test))
confusion_matrix(ys_test, gbcs.predict(Xs_test))

# changed subsampling: Choosing subsample < 1.0 leads to a reduction of variance and an increase in bias.
gbcs = GradientBoostingClassifier(init=None, learning_rate=0.1, loss='deviance',
              max_depth=3, max_features=None, max_leaf_nodes=None,
              min_samples_leaf=1, min_samples_split=2,
              min_weight_fraction_leaf=0.0, n_estimators=100,
              presort='auto', random_state=1, subsample=0.75, verbose=0,
              warm_start=False)
gbcs.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, gbcs)
print classification_report(ys_test, gbcs.predict(Xs_test))
confusion_matrix(ys_test, gbcs.predict(Xs_test))

# tried changing loss to exponential --> worse.
# change max features to sqrt --> worse.
# upped min_samples_split, from 2 to 3, improved. (other values did not cause improvement)

gbcs = GradientBoostingClassifier(init=None, learning_rate=0.1, loss='deviance',
              max_depth=3, max_features=None, max_leaf_nodes=None,
              min_samples_leaf=1, min_samples_split=3,
              min_weight_fraction_leaf=0.0, n_estimators=100,
              presort='auto', random_state=1, subsample=0.75, verbose=0,
              warm_start=False)
gbcs.fit(Xs_train, ys_train)
score_printout(Xs_test, ys_test, gbcs)
print classification_report(ys_test, gbcs.predict(Xs_test))
confusion_matrix(ys_test, gbcs.predict(Xs_test))

for i, n in enumerate(X.columns.get_values()):
    print i, n

fig, axs = plot_partial_dependence(gbcs, Xs_train, range(0, 6, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)

fig, axs = plot_partial_dependence(gbcs, Xs_train, range(6, 12, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)

fig, axs = plot_partial_dependence(gbcs, Xs_train, range(12, 18, 1), feature_names=X.columns.get_values(), n_jobs=-1, grid_resolution=100)

fig, axs = plot_partial_dependence(gbcs, Xs_train, range(18, 26, 1), feature_names=X.columns.get_values(), n_jobs=-1, grid_resolution=100)
plt.subplots_adjust(top=0.9)

# GBC without scaling works better than GBC with. Hmm.
gbc = GradientBoostingClassifier(init=None, learning_rate=0.1, loss='deviance',
              max_depth=3, max_features=None, max_leaf_nodes=None,
              min_samples_leaf=1, min_samples_split=3,
              min_weight_fraction_leaf=0.0, n_estimators=100,
              presort='auto', random_state=1, subsample=0.75, verbose=0,
              warm_start=False)
gbc.fit(X_train, y_train)
score_printout(X_test, y_test, gbc)
print classification_report(y_test, gbc.predict(X_test))
confusion_matrix(y_test, gbc.predict(X_test))

fig, axs = plot_partial_dependence(gbc, X_train, range(0, 6, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbc, X_train, range(6, 12, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbc, X_train, range(12, 18, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbc, X_train, range(18, 26, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

gbcRankedFeatures = sorted(zip(X.columns, gbc.feature_importances_), 
                          key=lambda pair: pair[1], 
                          reverse=True)

def make_feature_importance_plot(featuresAndImportances, numFeatures):
    topN = featuresAndImportances[:numFeatures]
    labels = [pair[0] for pair in topN]
    values = [pair[1] for pair in topN]
    ind = np.arange(len(values))
    width = 0.35   
    plt.bar(range(numFeatures), values, width=0.8)
    ax = plt.subplot(111)
    ax.set_xticks(ind+width)
    ax.set_xticklabels(labels, rotation=60, size=12)
    plt.xlabel('Feature', size=20)
    plt.ylabel('Importance', size=20)
    plt.show()
  

plt.figure()
make_feature_importance_plot(gbcRankedFeatures, 26)
plt.tight_layout()

gbc.feature_importances_

"""
Explanation: Gradient Boosting (with partial dependence plots)
End of explanation
"""


# GBC without scaling works better than GBC with. Hmm.
gbc = GradientBoostingClassifier(init=None, learning_rate=0.1, loss='deviance',
              max_depth=3, max_features=None, max_leaf_nodes=None,
              min_samples_leaf=1, min_samples_split=3,
              min_weight_fraction_leaf=0.0, n_estimators=100,
              presort='auto', random_state=1, subsample=0.75, verbose=0,
              warm_start=False)
gbc.fit(X_train, y_train)
score_printout(X_test, y_test, gbc)
print classification_report(y_test, gbc.predict(X_test))
confusion_matrix(y_test, gbc.predict(X_test))

paramGrid = {'n_estimators': [100, 200, 300],
             'learning_rate': [0.1, 0.5, 0.75, 0.2],
             'max_depth': [3, 4, 5, 6],
             'min_samples_leaf': [1, 2],
             'subsample': [0.75],
             'loss': ['deviance'],
             'max_features': [None, 'auto']
            }

gs = GridSearchCV(GradientBoostingClassifier(), 
                  param_grid=paramGrid, 
                  scoring='roc_auc', 
                  n_jobs=-1, 
                  cv=5, 
                  verbose=10)

gs.fit(X_train, y_train)

gs.best_estimator_

"""
Explanation: Grid search for best GBC
End of explanation
"""


gbcModel = GradientBoostingClassifier(init=None, learning_rate=0.1, loss='deviance',
              max_depth=3, max_features='auto', max_leaf_nodes=None,
              min_samples_leaf=2, min_samples_split=2,
              min_weight_fraction_leaf=0.0, n_estimators=200,
              presort='auto', random_state=None, subsample=0.75, verbose=0,
              warm_start=False)
gbcModel.fit(X_train, y_train)

gbcRankedFeatures = sorted(zip(X.columns, gbcModel.feature_importances_), 
                          key=lambda pair: pair[1], 
                          reverse=True)

def make_feature_importance_plot(featuresAndImportances, numFeatures):
    topN = featuresAndImportances[:numFeatures]
    labels = [pair[0] for pair in topN]
    values = [pair[1] for pair in topN]
    ind = np.arange(len(values))
    width = 0.35   
    plt.bar(range(numFeatures), values, width=0.8)
    ax = plt.subplot(111)
    ax.set_xticks(ind+width)
    ax.set_xticklabels(labels, rotation=60, size=12)
    plt.xlabel('Feature', size=20)
    plt.ylabel('Importance', size=20)
    plt.show()
  

score_printout(X_test, y_test, gbcModel)
print classification_report(y_test, gbcModel.predict(X_test))
confusion_matrix(y_test, gbcModel.predict(X_test))

gbcRankedFeatures = sorted(zip(X.columns, gbcModel.feature_importances_), 
                          key=lambda pair: pair[1], 
                          reverse=True)

plt.figure()
make_feature_importance_plot(gbcRankedFeatures, 15)
plt.tight_layout()

fig, axs = plot_partial_dependence(gbcModel, X_train, range(1, 6, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbcModel, X_train, range(6, 12, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbcModel, X_train, range(6, 8, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbcModel, X_train, range(18, 24, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

fig, axs = plot_partial_dependence(gbcModel, X_train, range(24, 26, 1), feature_names=X.columns.get_values(), n_jobs=3, grid_resolution=50)
plt.subplots_adjust(top=0.9)

"""
Explanation: Best model so far: GBC!
End of explanation
"""


names

from mpl_toolkits.mplot3d import Axes3D

# not quite getting the results I was expecting. Needs work.
fig = plt.figure()

names = X_train.columns

target_feature = (3, 21)
pdp, (x_axis, y_axis) = partial_dependence(gbcModel, target_feature,
                                           X=X_train, grid_resolution=50)
XX, YY = np.meshgrid(x_axis, y_axis)
Z = pdp.T.reshape(XX.shape).T
ax = Axes3D(fig)
surf = ax.plot_surface(XX, YY, Z, rstride=1, cstride=1, cmap=plt.cm.BuPu)
ax.set_xlabel(names[target_feature[0]])
ax.set_ylabel(names[target_feature[1]])
ax.set_zlabel('Partial dependence')
#  pretty init view
ax.view_init(elev=22, azim=122)
plt.colorbar(surf)
plt.suptitle('')
plt.subplots_adjust(top=0.9)

plt.show()

"""
Explanation: Use 3-D plot to investigate feature interactions for weak partial dependence plots... (weak effect may be masked by stronger interaction with other features)
End of explanation
"""


