%load_ext watermark
%watermark -a '' -u -d -v -p numpy,pandas,matplotlib,sklearn

"""
Explanation: Copyright (c) 2015, 2016 Sebastian Raschka
https://github.com/1iyiwei/pyml
MIT License
Python Machine Learning - Code Examples
Chapter 6 - Learning Best Practices for Model Evaluation and Hyperparameter Tuning
Deciphering the black magic
How to
* train a model
* evaluate a model
* fine tune (and debug) a model 
Note that the optional watermark extension is a small IPython notebook plugin that I developed to make the code reproducible. You can just skip the following line(s).
End of explanation
"""


from IPython.display import Image
%matplotlib inline

# Added version check for recent scikit-learn 0.18 checks
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_version

"""
Explanation: The use of watermark is optional. You can install this IPython extension via "pip install watermark". For more information, please see: https://github.com/rasbt/watermark.
Overview

Streamlining workflows with pipelines
Loading the Breast Cancer Wisconsin dataset
Combining transformers and estimators in a pipeline
Using k-fold cross-validation to assess model performance
The holdout method
K-fold cross-validation
Debugging algorithms with learning and validation curves
Diagnosing bias and variance problems with learning curves
Addressing overfitting and underfitting with validation curves
Fine-tuning machine learning models via grid search
Tuning hyperparameters via grid search
Algorithm selection with nested cross-validation
Looking at different performance evaluation metrics
Reading a confusion matrix
Optimizing the precision and recall of a classification model
Plotting a receiver operating characteristic
The scoring metrics for multiclass classification
Summary
End of explanation
"""


import pandas as pd

#wdbc_source = 'https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data'
wdbc_source = '../datasets/wdbc/wdbc.data'

df = pd.read_csv(wdbc_source, header=None)

df.shape

"""
Explanation: Streamlining workflows with pipelines
Pipeline 
* useful tool for parameter tuning and model evaluation.
* fit transformers (e.g. standardization) and predictors (e.g. SVM) together for easy manipulation
Dataset example
* breast cancer
Loading the Breast Cancer Wisconsin dataset
Malignant versus benign tumor cells based on 30 features (large enough for experiments in this chapter)
End of explanation
"""


if False:
    wdbc_another_source = 'https://raw.githubusercontent.com/1iyiwei/pyml/master/code/datasets/wdbc/wdbc.data'
    df = pd.read_csv(wdbc_another_source, header=None)

df.head()

from sklearn.preprocessing import LabelEncoder

X = df.loc[:, 2:].values
y = df.loc[:, 1].values
le = LabelEncoder()
y = le.fit_transform(y)
le.transform(['M', 'B'])

import numpy as np
np.bincount(y)

if Version(sklearn_version) < '0.18':
    from sklearn.cross_validation import train_test_split
else:
    from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = \
    train_test_split(X, y, test_size=0.20, random_state=1)

"""
Explanation: <hr>

Note:
If the link to the Breast Cancer Wisconsin dataset dataset provided above does not work for you, you can find a local copy in this repository at ./../datasets/wdbc/wdbc.data.
Or you could fetch it via
End of explanation
"""


from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline

# construct a pipeline consisting of different components
# each entry contains an identifier, and the corresponding component
pipe_lr = Pipeline([('scl', StandardScaler()),
                    ('pca', PCA(n_components=2)),
                    ('clf', LogisticRegression(random_state=1))])

# now a single fit or predict call will go through all components
pipe_lr.fit(X_train, y_train)
print('Test Accuracy: %.3f' % pipe_lr.score(X_test, y_test))
y_pred = pipe_lr.predict(X_test)

"""
Explanation: Combining transformers and estimators in a pipeline
We can chain together multiple ML components into a single entity called "pipeline" for easy manipulation.
This is also a common practice for building neural networks in deep learning libraries (e.g. TensorFlow and Torch7).
<img src='./images/06_01.png' width=80%>
End of explanation
"""


import numpy as np

# compute the data indices for each fold
if Version(sklearn_version) < '0.18':
    from sklearn.cross_validation import StratifiedKFold
    from sklearn.cross_validation import KFold
else:
    from sklearn.model_selection import StratifiedKFold
    from sklearn.model_selection import KFold
    

if Version(sklearn_version) < '0.18':
    kfold = KFold(y=y_train, n_folds=10, random_state=1)
else:
    kfold = KFold(n_splits=10, random_state=1).split(X_train, y_train)

scores = []
for k, (train, test) in enumerate(kfold):
    pipe_lr.fit(X_train[train], y_train[train]) # fit, not partial_fit
    score = pipe_lr.score(X_train[test], y_train[test])
    scores.append(score)
    print('Fold: %s, Class dist.: %s, Acc: %.3f' % (k+1,
          np.bincount(y_train[train]), score))
    
print('\nCV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

"""
Explanation: Using k-fold cross validation to assess model performance
Need to partition available data into 3 main categories:
* Training: fit model parameters
* Test: evaluate model performance
* Validation: select model hyper-parameters
Need to separate them to avoid over-fit
Data is often limited, so need smart partition
Different ways for partition
Review about parameters
Parameters
* optimized by algorithms, such as weight values $\mathbf{w}$ of a linear classifier
Hyper-parameters
* picked by humans, such as weight dimension $\| \mathbf{w} \|$ of a linear classifier
We can write code to help us pick hyper-parameters.
The holdout method
Test set
* only used for evaluation in the end (to make sure the model really generalizes)
Training & validation sets
* if we have only one model and one set of hyper-parameters, we can just train
* for multiple models or sets of hyper-parameters, we use a separate validation set
* train each model/hyper-parameter-set $\rightarrow$ validate $\rightarrow$ pick best
<img src='./images/06_02.png' width=80%> 
Example
* 3 models/hyper-parameter-sets: $M_1$, $M_2$, $M_3$
* optimize each model with the training set
* evaluate their performance using the validation set and pick the best one $M_{best} \in {M_1, M_2, M_3}$
* test $M_{best}$ with test set
Disadvantage of such fixed partition
* model performance can be sensitive to the particular choice
K-fold cross-validation
Example with $K=10$:
<img src='./images/06_03.png' width=80%> 
Vary the data partition to reduce variance for model training
* Usuall applied for training-validation data partition (with separate test data)
* But the principle can be applied in various ways, e.g. for test data
* The "test fold" is for validation in the figure
Question: use 1 model or 10 separate models for above?
The iterations should be independent to avoid cumulative learning
* fit() versus partial_fit()
How to choose k?
* large k for small data sets to ensure sufficient training data
End of explanation
"""


if Version(sklearn_version) < '0.18':
    kfold = StratifiedKFold(y=y_train, 
                            n_folds=10,
                            random_state=1)
else:
    kfold = StratifiedKFold(n_splits=10,
                            random_state=1).split(X_train, y_train)

scores = []
for k, (train, test) in enumerate(kfold):
    pipe_lr.fit(X_train[train], y_train[train]) # fit, not partial_fit
    score = pipe_lr.score(X_train[test], y_train[test])
    scores.append(score)
    print('Fold: %s, Class dist.: %s, Acc: %.3f' % (k+1,
          np.bincount(y_train[train]), score))
    
print('\nCV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

# even simpler code...
# allows parallel execution via n_jobs
# same result as above
if Version(sklearn_version) < '0.18':
    from sklearn.cross_validation import cross_val_score
else:
    from sklearn.model_selection import cross_val_score

scores = cross_val_score(estimator=pipe_lr,
                         X=X_train,
                         y=y_train,
                         cv=10,
                         n_jobs=1)
print('CV accuracy scores: %s' % scores)
print('CV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

"""
Explanation: Straitified k fold
* preserves class proportions within each fold instead of purely random (notice the code below)
* helps reduce variance in data and thus the trained model
End of explanation
"""


import matplotlib.pyplot as plt

if Version(sklearn_version) < '0.18':
    from sklearn.learning_curve import learning_curve
else:
    from sklearn.model_selection import learning_curve



pipe_lr = Pipeline([('scl', StandardScaler()),
                    ('clf', LogisticRegression(penalty='l2', random_state=0))])

train_sizes, train_scores, test_scores =\
                learning_curve(estimator=pipe_lr,
                               X=X_train,
                               y=y_train,
                               train_sizes=np.linspace(0.1, 1.0, 10), # 10 intervals
                               cv=10, # k = 10 for cross validation
                               n_jobs=1)

# the returned scores are stacked in columns (n_ticks, num_folds)
train_mean = np.mean(train_scores, axis=1)
train_std = np.std(train_scores, axis=1)
test_mean = np.mean(test_scores, axis=1)
test_std = np.std(test_scores, axis=1)

plt.plot(train_sizes, train_mean,
         color='blue', marker='o',
         markersize=5, label='training accuracy')

plt.fill_between(train_sizes,
                 train_mean + train_std,
                 train_mean - train_std,
                 alpha=0.15, color='blue')

plt.plot(train_sizes, test_mean,
         color='green', linestyle='--',
         marker='s', markersize=5,
         label='validation accuracy')

plt.fill_between(train_sizes,
                 test_mean + test_std,
                 test_mean - test_std,
                 alpha=0.15, color='green')

plt.grid()
plt.xlabel('Number of training samples')
plt.ylabel('Accuracy')
plt.legend(loc='lower right')
plt.ylim([0.8, 1.0])
plt.tight_layout()
# plt.savefig('./figures/learning_curve.png', dpi=300)
plt.show()

"""
Explanation: Debugging algorithms with learning and validation curves
Two simple and useful diagnostic tools for model training and performance 
Diagnosing bias and variance problems with learning curves
<img src='./images/06_04.png' width=70% align=right>
Learning curve
* plots performance against data size
Goal
* low bias
* low variance
Observations: when training size $\uparrow$ 
* training performance $\downarrow$
* validation performance $\uparrow$
High bias
* model too simple (with respect to the data)
* increase model complexity
* decrease regularization
High variance
* model too complex
* dimensionality reduction
* decrease model complexity
* increase regularization
Plot learning curves via scikit-learn
End of explanation
"""


if Version(sklearn_version) < '0.18':
    from sklearn.learning_curve import validation_curve
else:
    from sklearn.model_selection import validation_curve



param_range = [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]
train_scores, test_scores = validation_curve(
                estimator=pipe_lr, 
                X=X_train, 
                y=y_train, 
                param_name='clf__C', # (inverse) regularization strength
                param_range=param_range,
                cv=10)

#print('all parameters: %s' % pipe_lr.get_params())

train_mean = np.mean(train_scores, axis=1)
train_std = np.std(train_scores, axis=1)
test_mean = np.mean(test_scores, axis=1)
test_std = np.std(test_scores, axis=1)

plt.plot(param_range, train_mean, 
         color='blue', marker='o', 
         markersize=5, label='training accuracy')

plt.fill_between(param_range, train_mean + train_std,
                 train_mean - train_std, alpha=0.15,
                 color='blue')

plt.plot(param_range, test_mean, 
         color='green', linestyle='--', 
         marker='s', markersize=5, 
         label='validation accuracy')

plt.fill_between(param_range, 
                 test_mean + test_std,
                 test_mean - test_std, 
                 alpha=0.15, color='green')

plt.grid()
plt.xscale('log')
plt.legend(loc='lower right')
plt.xlabel('Parameter C')
plt.ylabel('Accuracy')
plt.ylim([0.8, 1.0])
plt.tight_layout()
# plt.savefig('./figures/validation_curve.png', dpi=300)
plt.show()

"""
Explanation: Addressing over- and underfitting with validation curves
Validation curves are similar to learning curves, but plot against parameters instead of training size.
The sklearn API allows selection of parameter and range.
End of explanation
"""


from sklearn.svm import SVC
if Version(sklearn_version) < '0.18':
    from sklearn.grid_search import GridSearchCV
else:
    from sklearn.model_selection import GridSearchCV

pipe_svc = Pipeline([('scl', StandardScaler()),
            ('clf', SVC(random_state=1))])

param_range = [0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0]

# can specify multiple grids, e.g. 2 below
param_grid = [{'clf__C': param_range, 
               'clf__kernel': ['linear']},
                 {'clf__C': param_range, 
                  'clf__gamma': param_range, 
                  'clf__kernel': ['rbf']}]

gs = GridSearchCV(estimator=pipe_svc, 
                  param_grid=param_grid, 
                  scoring='accuracy', 
                  cv=10,
                  n_jobs=-1)
gs = gs.fit(X_train, y_train)
print(gs.best_score_)
print(gs.best_params_)

clf = gs.best_estimator_
clf.fit(X_train, y_train)
print('Test accuracy: %.3f' % clf.score(X_test, y_test))

"""
Explanation: Which $C$ value should we pick from the experiment above?
Hint: the goal includes:
* low bias, i.e. high accuracy
* low variance, i.e. similar training and test/validation accuracy
Fine-tuning machine learning models via search
Manually selecting hyper-parameters can be tedious
Automatic search through a set of hyper-parameters, and find the best one
* write your own code (e.g. loops)
* grid and random search in scikit-learn
Tuning hyperparameters via grid search
Given a set of parameters $p_i$ each with $N_i$ different values, search through all combinations $\prod_i N_i$.
Each such set is like a grid.
Example
Two parameters with options ${a_1, a_2}$ and ${b_1, b_2, b_3 }$
Grid search will try all combinations: 
* $(a_1, b_1)$
* $(a_1, b_2)$
* $(a_1, b_3)$
* $(a_2, b_1)$
* $(a_2, b_2)$
* $(a_3, b_3)$
End of explanation
"""


if Version(sklearn_version) < '0.18':
    from sklearn.grid_search import RandomizedSearchCV
else:
    from sklearn.model_selection import RandomizedSearchCV
    
from sklearn.svm import SVC
from scipy.stats import uniform

pipe_svc = Pipeline([('scl', StandardScaler()), ('clf', SVC(random_state=1))])

param_range = [0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0]
param_range = uniform(scale=1000)
 
param_distributions = {'clf__C': param_range,
                       'clf__kernel': ['linear', 'rbf'],
                       'clf__gamma': param_range}

rs = RandomizedSearchCV(estimator=pipe_svc,
                        param_distributions=param_distributions,
                        scoring='accuracy',
                        n_iter=100, # number of settings to try
                        cv=10,
                        n_jobs=1)

rs.fit(X_train, y_train)
print(rs.best_score_)
print(rs.best_params_)

clf = rs.best_estimator_
clf.fit(X_train, y_train)
print('Test accuracy: %.3f' % clf.score(X_test, y_test))

"""
Explanation: Tuning hyperparameters via randomized search
Grid search can suffer from curse of dimensionality
* $i$ parameters each with $N_i$ different values, the total combination is $\prod_i N_i$
* $i = 10$, $N_i = 10$ $\rightarrow$ $\prod_i N_i = 10^{10}$
Random sampling with a per-defined budget (number of samples) can help
End of explanation
"""


# cross validation for outer loop
# grid search for inner loop

outer_cv = 5
inner_cv = 2

gs = GridSearchCV(estimator=pipe_svc,
                  param_grid=param_grid,
                  scoring='accuracy',
                  cv=inner_cv)

scores = cross_val_score(gs, X_train, y_train, scoring='accuracy', cv=outer_cv)
print('CV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

# compare decision tre with SVM above
# and see whose performance is better
from sklearn.tree import DecisionTreeClassifier

gs = GridSearchCV(estimator=DecisionTreeClassifier(random_state=0),
                  param_grid=[{'max_depth': [1, 2, 3, 4, 5, 6, 7, None]}],
                  scoring='accuracy',
                  cv=inner_cv)
scores = cross_val_score(gs, X_train, y_train, scoring='accuracy', cv=outer_cv)
print('CV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

"""
Explanation: Algorithm selection with nested cross-validation
Partition data into training, validation, and test sets in nested hierarchy
<img src='./images/06_07.png' width=80%>
Another way to think about it: rotate the test set within the original set
Generalize k-fold with fixed test set
<img src='./images/06_02.png' width=80%>
<img src='./images/06_03.png' width=80%>
How to code?
Use grid search instead of classifier/pipeline with cross_val_score()
End of explanation
"""


from sklearn.metrics import confusion_matrix

pipe_svc.fit(X_train, y_train)
y_pred = pipe_svc.predict(X_test)
confmat = confusion_matrix(y_true=y_test, y_pred=y_pred)
print(confmat)

fig, ax = plt.subplots(figsize=(2.5, 2.5))
ax.matshow(confmat, cmap=plt.cm.Blues, alpha=0.3)
for i in range(confmat.shape[0]):
    for j in range(confmat.shape[1]):
        ax.text(x=j, y=i, s=confmat[i, j], va='center', ha='center')

plt.xlabel('predicted label')
plt.ylabel('true label')

plt.tight_layout()
# plt.savefig('./figures/confusion_matrix.png', dpi=300)
plt.show()

"""
Explanation: Looking at different performance evaluation metrics
<a href="https://en.wikipedia.org/wiki/Precision_and_recall"><img src="https://upload.wikimedia.org/wikipedia/commons/thumb/2/26/Precisionrecall.svg/320px-Precisionrecall.svg.png" align=right></a>
Quality
* <b>T</b>rue
* <b>F</b>alse
Prediction
* <b>P</b>ositive
* <b>N</b>egative
ACC (accuracy) = $\frac{TP+TN}{TP+TN+FP+FN}$
<br>
ERR (error) = 1 - ACC
TPR (true positive rate) = $\frac{TP}{TP+FN}$
<br>
FPR (false positive rate) = $\frac{FP}{FP+TN}$
PRE (precision) = $\frac{TP}{TP+FP}$
<br>
REC (recall) = TPR
F1 = $2\frac{PRE \times REC}{PRE + REC}$
Reading a confusion matrix
<img src='./images/06_08.png' width=40%>
End of explanation
"""


le.transform(['M', 'B'])

confmat = confusion_matrix(y_true=y_test, y_pred=y_pred)
print(confmat)

"""
Explanation: Additional Note
Remember that we previously encoded the class labels so that malignant samples are the "postive" class (1), and benign samples are the "negative" class (0):
End of explanation
"""


confmat = confusion_matrix(y_true=y_test, y_pred=y_pred)
print(confmat)

"""
Explanation: Next, we printed the confusion matrix like so:
End of explanation
"""


confmat = confusion_matrix(y_true=y_test, y_pred=y_pred, labels=[1, 0])
print(confmat)

"""
Explanation: Note that the (true) class 0 samples that are correctly predicted as class 0 (true negatives) are now in the upper left corner of the matrix (index 0, 0). In order to change the ordering so that the true negatives are in the lower right corner (index 1,1) and the true positves are in the upper left, we can use the labels argument like shown below:
End of explanation
"""


from sklearn.metrics import precision_score, recall_score, f1_score

print('Precision: %.3f' % precision_score(y_true=y_test, y_pred=y_pred))
print('Recall: %.3f' % recall_score(y_true=y_test, y_pred=y_pred))
print('F1: %.3f' % f1_score(y_true=y_test, y_pred=y_pred))

from sklearn.metrics import make_scorer

scorer = make_scorer(f1_score, pos_label=0)

c_gamma_range = [0.01, 0.1, 1.0, 10.0]

param_grid = [{'clf__C': c_gamma_range,
               'clf__kernel': ['linear']},
              {'clf__C': c_gamma_range,
               'clf__gamma': c_gamma_range,
               'clf__kernel': ['rbf']}]

gs = GridSearchCV(estimator=pipe_svc,
                  param_grid=param_grid,
                  scoring=scorer,
                  cv=10,
                  n_jobs=-1)
gs = gs.fit(X_train, y_train)
print(gs.best_score_)
print(gs.best_params_)

"""
Explanation: We conclude:
Assuming that class 1 (malignant) is the positive class in this example, our model correctly classified 71 of the samples that belong to class 0 (true negatives) and 40 samples that belong to class 1 (true positives), respectively. However, our model also incorrectly misclassified 1 sample from class 0 as class 1 (false positive), and it predicted that 2 samples are benign although it is a malignant tumor (false negatives).
Optimizing the precision and recall of a classification model
Previously we have been using accuracy as the score for optimization
* e.g. grid search
We can easily plug in other performance metrics.
End of explanation
"""


# modified from http://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html
import numpy as np
from sklearn.metrics import roc_curve

y_truth = np.array([0, 0, 1, 1])

#y_score = np.array([0, 0, 1, 1]) # completely right
y_score = np.array([0.1, 0.4, 0.35, 0.8]) # partially right
#y_score = np.array([1, 1, 0, 0]) # completely wrong

# compute fpr and tpr with thresholds
fpr, tpr, thresholds = roc_curve(y_truth, y_score, pos_label=1)

print(fpr)
print(tpr)
print(thresholds)

"""
Explanation: Plotting a receiver operating characteristic
https://en.wikipedia.org/wiki/Receiver_operating_characteristic
Plot TPR versus FPR of a model by varing the decision threshold
* roc_curve() in scikit-learn
End of explanation
"""


from sklearn.metrics import roc_curve, auc
from scipy import interp

pipe_lr = Pipeline([('scl', StandardScaler()),
                    ('pca', PCA(n_components=2)),
                    ('clf', LogisticRegression(penalty='l2', 
                                               random_state=0, 
                                               C=100.0))])

# intentionally use only 2 features to make the task harder and the curves more interesting
X_train2 = X_train[:, [4, 14]]
X_test2 = X_test[:, [4, 14]]


if Version(sklearn_version) < '0.18':
    cv = StratifiedKFold(y_train, 
                         n_folds=3, 
                         random_state=1)
    
else:
    cv = list(StratifiedKFold(n_splits=3, 
                              random_state=1).split(X_train, y_train))

fig = plt.figure(figsize=(7, 5))

mean_tpr = 0.0
mean_fpr = np.linspace(0, 1, 100)
all_tpr = []

for i, (train, test) in enumerate(cv):
    probas = pipe_lr.fit(X_train2[train],
                         y_train[train]).predict_proba(X_train2[test])

    fpr, tpr, thresholds = roc_curve(y_train[test],
                                     probas[:, 1],
                                     pos_label=1)
    mean_tpr += interp(mean_fpr, fpr, tpr)
    #mean_tpr[0] = 0.0
    roc_auc = auc(fpr, tpr)
    plt.plot(fpr,
             tpr,
             lw=1,
             label='ROC fold %d (area = %0.2f)'
                   % (i+1, roc_auc))

plt.plot([0, 1],
         [0, 1],
         linestyle='--',
         color=(0.6, 0.6, 0.6),
         label='random guessing')

mean_tpr /= len(cv)
mean_tpr[0] = 0.0
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
plt.plot(mean_fpr, mean_tpr, 'k--',
         label='mean ROC (area = %0.2f)' % mean_auc, lw=2)
plt.plot([0, 0, 1],
         [0, 1, 1],
         lw=2,
         linestyle=':',
         color='black',
         label='perfect performance')

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('false positive rate')
plt.ylabel('true positive rate')
plt.title('Receiver Operator Characteristic')
plt.legend(loc="lower right")

plt.tight_layout()
# plt.savefig('./figures/roc.png', dpi=300)
plt.show()

pipe_lr = pipe_lr.fit(X_train2, y_train)
y_pred2 = pipe_lr.predict(X_test2)

# area under curve (AUC) can give us more information than accuracy
# which is a single cut-off point (e.g. 0.5 for probabilities) of ROC
from sklearn.metrics import roc_auc_score, accuracy_score
print('ROC AUC: %.3f' % roc_auc_score(y_true=y_test, y_score=y_pred2))
print('Accuracy: %.3f' % accuracy_score(y_true=y_test, y_pred=y_pred2))

# compare with using all the original training data
# notice the better ROC
fig = plt.figure(figsize=(7, 5))

mean_tpr = 0.0
mean_fpr = np.linspace(0, 1, 100)
all_tpr = []

for i, (train, test) in enumerate(cv):
    probas = pipe_lr.fit(X_train[train],
                         y_train[train]).predict_proba(X_train[test])

    fpr, tpr, thresholds = roc_curve(y_train[test],
                                     probas[:, 1],
                                     pos_label=1)
    mean_tpr += interp(mean_fpr, fpr, tpr)
    #mean_tpr[0] = 0.0
    roc_auc = auc(fpr, tpr)
    plt.plot(fpr,
             tpr,
             lw=1,
             label='ROC fold %d (area = %0.2f)'
                   % (i+1, roc_auc))

plt.plot([0, 1],
         [0, 1],
         linestyle='--',
         color=(0.6, 0.6, 0.6),
         label='random guessing')

mean_tpr /= len(cv)
mean_tpr[0] = 0.0
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
plt.plot(mean_fpr, mean_tpr, 'k--',
         label='mean ROC (area = %0.2f)' % mean_auc, lw=2)
plt.plot([0, 0, 1],
         [0, 1, 1],
         lw=2,
         linestyle=':',
         color='black',
         label='perfect performance')

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('false positive rate')
plt.ylabel('true positive rate')
plt.title('Receiver Operator Characteristic')
plt.legend(loc="lower right")

plt.tight_layout()
# plt.savefig('./figures/roc.png', dpi=300)
plt.show()

"""
Explanation: If threshold is $0.8$
* classification = $[0, 0, 0, 1]$
* ground truth = $[0, 0, 1, 1]$
* tpr = $0.5$, fpr = $0$
If threshold is $0.35$
* classification = $[0, 1, 1, 1]$
* ground truth = $[0, 0, 1, 1]$
* tpr = $1$, fpr = $0.5$
End of explanation
"""


# make a scorer and plug in to stuff we showed above
pre_scorer = make_scorer(score_func=precision_score, 
                         pos_label=1, 
                         greater_is_better=True, 
                         average='micro')

"""
Explanation: Can plot other types of curves such as precision-recall curve
* exercise
The scoring metrics for multiclass classification
Performance metrics above are defined for binary classification
Extension for multi-class classification (one versus rest):
* micro
* macro
For $K$ classes, PRE (precision) can be computed as follows:
$$
PRE_{micro} = \frac{\sum_{i=1}^K TP_i}{\sum_{i=1}^K TP_i + \sum_{i=1}^K FP_i}
$$
$$
PRE_{macro} = \frac{\sum_{i=1}^K PRE_i}{K}
$$
Ditto other metrics (e.g. recall, accuracy, F1)
Micro
* weigh each sample equally
Macro
* weigh each class equally
* can specify weights in API
One versus one:
<a href="https://www.microsoft.com/en-us/research/people/cmbishop/">
<img src="./images/one_vs_one.svg" alt="one vs one">
</a>
One versus all/rest:
<a href="https://www.microsoft.com/en-us/research/people/cmbishop/">
<img src="./images/one_vs_rest.svg" alt="one vs rest">
</a>
End of explanation
"""


