from IPython.display import Image
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import time
%matplotlib inline

"""
Explanation: Goals of this Lesson


Understand the Bias-Variance Tradeoff


Ensembling Methods

Bagging
Voting 
Stacking



References

Chapter 8 of Elements of Statistical Learning by Hastie, Tibshirani, Friedman
A Few Useful Things to Know about Machine Learning
SciKit-Learn's documentation on ensemble methods

0.  Preliminaries
First we need to import Numpy, Pandas, MatPlotLib...
End of explanation
"""


### function for shuffling the data and labels
def shuffle_in_unison(features, labels):
    rng_state = np.random.get_state()
    np.random.shuffle(features)
    np.random.set_state(rng_state)
    np.random.shuffle(labels)
    
### calculate classification errors
# return a percentage: (number misclassified)/(total number of datapoints)
def calc_classification_error(predictions, class_labels):
    n = predictions.size
    num_of_errors = 0.
    for idx in xrange(n):
        if (predictions[idx] >= 0.5 and class_labels[idx]==0) or (predictions[idx] < 0.5 and class_labels[idx]==1):
            num_of_errors += 1
    return num_of_errors/n

"""
Explanation: Again we need functions for shuffling the data and calculating classification errrors.
End of explanation
"""


from sklearn.decomposition import PCA

# load the 403 x 360,000 matrix
br_paintings = np.load(open('../data/bob_ross/bob_ross_paintings.npy','rb'))

# perform PCA again
pca = PCA(n_components=400)
start_time = time.time()
pca_paintings = pca.fit_transform(br_paintings)
end_time = time.time()

# remove the br_paintings from memory
br_paintings = None

print "Training took a total of %.2f seconds." %(end_time-start_time)
print "Preserved percentage of original variance: %.2f%%" %(pca.explained_variance_ratio_.sum() * 100) 
print "Dataset is now of size: %d x %d"%(pca_paintings.shape)

"""
Explanation: 0.1  Load the dataset of paintings
We are going to use the Bob Ross paintings dataset throughout this session.  Let's again load the data and run PCA...
End of explanation
"""


br_labels_data = pd.read_csv('../data/bob_ross/elements-by-episode.csv')
br_labels_data.head()

"""
Explanation: We want to formulate a binary classification problem.  In the data folder there's a file that has labels denoting what is in each painting (tree, mountain, etc.).  Let's load it...
End of explanation
"""


labels = (br_labels_data['HILLS'] + br_labels_data['MOUNTAIN'] + br_labels_data['MOUNTAINS'] \
          + br_labels_data['SNOWY_MOUNTAIN'] > 0).astype('int8').as_matrix()
print "Contains mountain?: "+str(bool(labels[5]))
recon_img = pca.inverse_transform(pca_paintings[5,:])
plt.imshow(np.reshape(recon_img, (300, 400, 3)))
plt.show()

"""
Explanation: Let's make two classes: 1 = 'painting contains hill or mountain', 0 = 'doesn't contain hill/mountain':
End of explanation
"""


# set the random number generator for reproducability
np.random.seed(182)

# shuffle data
N = pca_paintings.shape[0]
shuffle_in_unison(pca_paintings, labels)

# split into train and test sets
train_features = pca_paintings[:int(.8*N), :]
test_features = pca_paintings[int(.8*N):, :]
train_labels = labels[:int(.8*N)]
test_labels = labels[int(.8*N):]

"""
Explanation: Make training and test split...
End of explanation
"""


np.random.seed(182)

from sklearn.linear_model import LogisticRegression

# initialize and train a logistic regression model
lr_model = LogisticRegression()
lr_model.fit(train_features, train_labels)

# compute error on test data
lr_predictions = lr_model.predict(test_features)
one_model_test_error_rate = calc_classification_error(lr_predictions, test_labels)

print "Classification error on test set: %.2f%%" %(one_model_test_error_rate*100)
# compute the baseline error since the classes are imbalanced
print "Baseline Error: %.2f%%" %((sum(test_labels)*100.)/len(test_labels))

"""
Explanation: 0.2  Run a baseline classifier
In order to see the improvements that ensembling provides, let's train a baseline logistic regression classifier for later comparison.
End of explanation
"""


### function for bootstrap resampling
def bootstrap_resample(features, labels, n_resamples):
    N = features.shape[0]
    idxs = np.arange(N)
    # numpy's choice() handles the sampling w/ replacement
    resampled_idxs = np.random.choice(idxs, size=(N,n_resamples), replace=True)
    boot_samps_x = []
    boot_samps_y = []
    for i in xrange(n_resamples):
        boot_samps_x.append(features[resampled_idxs[:,i],:])
        boot_samps_y.append(labels[resampled_idxs[:,i]])
    return boot_samps_x, boot_samps_y

"""
Explanation: 1.  The Bias-Variance Tradeoff
When faced with an important decision, its common to ask multiple people for their advice.  Why should a classification decision be any different?  If computer power is not a limiting factor--which is usually the case--why not train multiple classifiers and combine their predictions?  This is exactly what ensembling classifiers does.  In this section we'll cover three methods for combining classifiers: bagging, averaging, and stacking.  But first, let's examine why one classifier is usually not enough.  It can be formalized as a tradeoff between bias and variance.  
Recall the squared loss function: $$\mathcal{L} = \sum_{i}^{N} (y_{i} - f(\mathbf{x}{i}))^{2}. $$  This loss is over a particular training set {$\mathbf{X}, \mathbf{y}$} but we are really interested in the loss over all possible datasets we could have observed, ${\mathbf{X}, \mathbf{y}} \sim p(\mathcal{D})$: $$\mathbb{E}{p(\mathcal{D})}[\mathcal{L}] = \mathbb{E}{p(\mathcal{D})}[(y{i} - f(\mathbf{x}{i}))^{2}]. $$  After some algebraic manipulations, we can re-write the expected loss as $$\mathbb{E}{p(\mathcal{D})}[\mathcal{L}] = (f^{}(\mathbf{x}{i}) - \mathbb{E}[f(\mathbf{x}{i})])^{2} + \text{Var}[f(\mathbf{x}_{i})] + \text{constant (error)}. $$  The first term, $(f^{}(\mathbf{x}{i}) - \mathbb{E}[f(\mathbf{x}{i})])^{2}$, is the squared difference between the expected value of the classifier $f$ and the perfect, true classifier $f^{}$.  This difference is known as the bias of a classifier.  For instance, a linear model has a strong bias since its functional form is rather simple (unless the optimial classifier is also a linear function).  The second term, $\text{Var}[f(\mathbf{x}_{i})]$, is the variance of our classifier.  Basically, this term captures the variability in outputs.  The main point is that if a classifier has low bias, meaning it is a very powerful function, then it will usually have high variance since this power allows it to generate a wide range of outputs.  And vice versa.  What I just said can be represented graphically as 

Ensembling classifiers all but always produces better performance because it reduces variance without incurring bias*.
2.  Bootstrap Aggregating ('Bagging')
In section 1, when I mentioned using multiple classifiers, you probably thought I was talking about training and combining several different kinds of classifiers.  We will do that.  But first we'll discuss something simpler: training the same classifier on multiple datasets.  
2.1 Bootstrap Resampling
We (always) want more training data, but unfortunately, it is not available.  We can use the training data we do have and resample it with replacement to generate additional 'fake' datasets.  Formally, our original dataset is $${\mathbf{y}, \mathbf{X}} \sim p(\mathcal{D}),$$ where $p(\mathcal{D})$ is the unknown population distribution.  We then treat the original data as a substitute for the population, writing $${\mathbf{\tilde y}, \mathbf{ \tilde X}}{1} , {\mathbf{\tilde y}, \mathbf{ \tilde X}}{2}, ... \sim {\mathbf{y}, \mathbf{X}}.$$  ${\mathbf{\tilde y}, \mathbf{ \tilde X}}$ are called bootstrap (re)samples.  Usually, they contain the same number of instances as the original training set.  A diagram showing sampling with replacement is below  

Now let's write some Python code to generate Bootstrap samples from a given dataset...
End of explanation
"""


def fit_and_predict_on_bootstrap_samples(model, bootstrapped_features, bootstrapped_labels, \
                                         test_features, n_bootstrap_samples):
    n_test = test_features.shape[0]
    ensemb_probs = np.zeros((n_test,))
    ensemb_preds = np.zeros((n_test,))
    for idx in xrange(n_bootstrap_samples):
        print "training model #%d" %(idx+1)
        model.fit(bootstrapped_features[idx], bootstrapped_labels[idx])
        ensemb_probs += model.predict_proba(test_features)[:,1]
        ensemb_preds += model.predict(test_features)
    ensemb_probs /= n_bootstrap_samples
    ensemb_preds /= n_bootstrap_samples
    ensemb_probs = np.around(ensemb_probs)
    ensemb_preds = np.around(ensemb_preds)
    return ensemb_probs, ensemb_preds

"""
Explanation: 2.2 Training on Bootstrap Samples
Given the Bootstrap datasets, we next train a classifier on each dataset and then combine their predictions.  These classifiers can all be instances of the model, as is the case in the code below, or different ones.  Once the models are trained, we can combine them in two way: by averaging the probabilities or the predictions.  Both methods are shown below.
End of explanation
"""


np.random.seed(182)

n_bootstrap_samples = 7
bootstrapped_features, bootstrapped_labels = bootstrap_resample(train_features, train_labels, n_bootstrap_samples)

ensembled_probs, ensembled_preds = \
fit_and_predict_on_bootstrap_samples(lr_model, bootstrapped_features, bootstrapped_labels, \
                                     test_features, n_bootstrap_samples)

print
print "Averaging probabilities: classification error on test set is %.2f%%" \
%(calc_classification_error(ensembled_probs, test_labels)*100)
print "Averaging predictions: classification error on test set is %.2f%%" \
%(calc_classification_error(ensembled_preds, test_labels)*100)
print
print "One logistic regression model error: %.2f%%"%(one_model_test_error_rate*100)
# compute the baseline error since the classes are imbalanced
print "Baseline error: %.2f%%" %((sum(test_labels)*100.)/len(test_labels))

"""
Explanation: Now run the function...
End of explanation
"""


np.random.seed(182)

# import the three new classifiers
from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB

# initialize models
d_tree_model = DecisionTreeClassifier()
knn_model = KNeighborsClassifier(n_neighbors=5)
nb_model = GaussianNB()

# fit models
d_tree_model.fit(train_features, train_labels)
knn_model.fit(train_features, train_labels)
nb_model.fit(train_features, train_labels)

# predict on test data
tree_predictions = d_tree_model.predict(test_features)
knn_predictions = knn_model.predict(test_features)
nb_predictions = nb_model.predict(test_features)

# average predictions
# add in the logistic regression predictions calcuated previously
avg_predictions = np.around((tree_predictions + knn_predictions + nb_predictions + lr_predictions)/4.)

print "Averaging predictions: classification error on test set is %.2f%%" \
%(calc_classification_error(avg_predictions, test_labels)*100)
print
print "One logistic regression model error: %.2f%%"%(one_model_test_error_rate*100)
# compute the baseline error since the classes are imbalanced
print "Baseline error: %.2f%%" %((sum(test_labels)*100.)/len(test_labels))

"""
Explanation: 2.  Voting
Now we'll consider combining a diverse set of classifers, each trained on an identitical copy of the data.  But first, we need to quickly introduce three new types of classifiers.  Unfortunately, we don't have enough time to cover each in detail and will have to use them as somewhat of a black-box.
2.1 Overview of Three Classifiers: Decision Tree, k-Nearest Neighbors, and Naive Bayes

SciKit-Learn Documentation:
- Decision Tree
- k-Nearest Neighbors
- Naive Bayes
2.2 Voting via Averaging Predictions
$$ \hat f(\mathbf{x}{i}) = \frac{1}{4} \hat y{\text{DT}} + \frac{1}{4} \hat y_{\text{kNN}} + \frac{1}{4} \hat y_{\text{NB}} + \frac{1}{4} \hat y_{\text{LogReg}}.$$
End of explanation
"""


np.random.seed(182)

### TRAINING
# calculate probabilities on the training data
tree_probs = d_tree_model.predict_proba(train_features)[:,1][np.newaxis].T
knn_probs = knn_model.predict_proba(train_features)[:,1][np.newaxis].T
nb_probs = nb_model.predict_proba(train_features)[:,1][np.newaxis].T
logReg_probs = lr_model.predict_proba(train_features)[:,1][np.newaxis].T

# combine into a new 'feature' matrix
train_probs_matrix = np.hstack([tree_probs, knn_probs, nb_probs, logReg_probs])

# train logistic regression
meta_classifier = LogisticRegression()
meta_classifier.fit(train_probs_matrix, train_labels)

# plot the weights learned for each classifier
f,ax = plt.subplots()
ticks = np.arange(4)
ax.bar(ticks, meta_classifier.coef_[0])
ax.set_xticks(ticks+.4)
ax.set_xticklabels(['Decision Tree', 'kNN', 'Naive Bayes', 'Log. Regression'])
ax.set_title('Weights Learned for Each Classifier')
plt.show()

### TESTING
# calculate probabilities on the test data
tree_probs = d_tree_model.predict_proba(test_features)[:,1][np.newaxis].T
knn_probs = knn_model.predict_proba(test_features)[:,1][np.newaxis].T
nb_probs = nb_model.predict_proba(test_features)[:,1][np.newaxis].T
logReg_probs = lr_model.predict_proba(test_features)[:,1][np.newaxis].T

# combine into a new 'feature' matrix
test_probs_matrix = np.hstack([tree_probs, knn_probs, nb_probs, logReg_probs])

stacked_predictions = meta_classifier.predict(test_probs_matrix)

print "Averaging predictions: classification error on test set is %.2f%%" %(calc_classification_error(stacked_predictions, test_labels)*100)
print
print "One logistic regression model error: %.2f%%"%(one_model_test_error_rate*100)
# compute the baseline error since the classes are imbalanced
print "Baseline error: %.2f%%" %((sum(test_labels)*100.)/len(test_labels))

train_stacked_predictions = meta_classifier.predict(train_probs_matrix)

print "Stacking train error: %.2f%%" %(calc_classification_error(train_stacked_predictions, train_labels))

"""
Explanation: 3.  Stacking Models
When we performed the averaging above, we used this formula: $$ \hat f(\mathbf{x}{i}) = \frac{1}{4} \hat y{\text{DT}} + \frac{1}{4} \hat y_{\text{kNN}} + \frac{1}{4} \hat y_{\text{NB}} + \frac{1}{4} \hat y_{\text{LogReg}}.$$  That is, we gave each classifier equal weighting.  While this approach is reasonable, probably it would be better to give an unequal weighting to each classifer, allowing the 'smartest' model to contribute most to the decision.  We can accomplish this by training a second-level logistic regression classifier on the predicted probabilities: $$ \hat f(\mathbf{x}{i}) = \sigma ( \alpha{1} f_{\text{DT}}(\mathbf{x}{i}) + \alpha{2}f_{\text{kNN}}(\mathbf{x}{i}) + \alpha{3}f_{\text{NB}}(\mathbf{x}{i}) + \alpha{4}f_{\text{LogReg}}(\mathbf{x}_{i})) \text{ where } \sigma (\cdot) \text{ is the logistic function}.$$  A depiction of the pipeline is below:
End of explanation
"""


from sklearn.datasets import fetch_lfw_pairs
lfw_train = fetch_lfw_pairs(subset='train')
lfw_test = fetch_lfw_pairs(subset='test')
lfw_train_pairs = lfw_train['pairs']
lfw_train_targets = lfw_train['target']
lfw_test_pairs = lfw_test['pairs']
lfw_test_targets = lfw_test['target']

print "The training data is of size: %d instances x %d faces x %d pixels x %d pixels" %(lfw_train_pairs.shape)
print "The test data is of size: %d instances x %d faces x %d pixels x %d pixels" %(lfw_test_pairs.shape)

"""
Explanation: Ah!  The training error is zero percent--a dead giveaway that the added power of stacking caused us to overfit.  One way of prevent this overfitting is to split the training data and fit the base learners and meta-classifier on different subsets.
<span style="color:red">STUDENT ACTIVITY (until end of session)</span>
4.  Mini Competition
For the remainder of the session, we'll have a mini predictive modeling competition on the Labeled Faces in the Wild dataset.  Your task to get as low of an error on the test set as possible by empolying all the techniques we covered today: dimensionality reduction (PCA), Bagging, stacking, etc.  Be mindful to save some data as a validation set.
The code below will load the data, display an image, convert it to feature vectors, and train a logistic regression classifier.
End of explanation
"""


face_idx=0

# subplot containing first image
ax1 = plt.subplot(1,2,1)
ax1.imshow(lfw_train_pairs[face_idx,0,:,:],cmap='Greys_r')

# subplot containing second image
ax2 = plt.subplot(1,2,2)
ax2.imshow(lfw_train_pairs[face_idx,1,:,:],cmap='Greys_r')
plt.show()

"""
Explanation: Let's visualize the images...
End of explanation
"""


train_x = np.reshape(lfw_train_pairs, (2200, 5828))
train_y = lfw_train_targets
test_x = np.reshape(lfw_test_pairs, (1000, 5828))
test_y = lfw_test_targets

# print the shapes just to check its what we expect
print train_x.shape
print train_y.shape
print test_x.shape
print test_y.shape

"""
Explanation: As a last step of pre-processing, let's flatten the data tensor...
End of explanation
"""


np.random.seed(182)

# initialize and train a logistic regression model
lr_model = LogisticRegression()
lr_model.fit(train_x, train_y)

# compute error on test data
lr_predictions = lr_model.predict(test_x)
one_model_test_error_rate = calc_classification_error(lr_predictions, test_y)

print "Classification error on test set: %.2f%%" %(one_model_test_error_rate*100)
# compute the baseline error since the classes are imbalanced
print "Baseline Error: %.2f%%" %((sum(test_y)*100.)/len(test_y))

"""
Explanation: And just to get you started, here's code to train a logistic regression classifier...
End of explanation
"""


