import pandas as pd
import numpy as np
from collections import OrderedDict, defaultdict

survey_data = pd.read_csv("https://raw.githubusercontent.com/samtalasila/e-bike-survey-response-results/master/E-Bike_Survey_Responses.csv")
print(survey_data.shape)
survey_data.head()

"""
Explanation: E-Bike survey analysis
This notebook contains my analysis for how to charactrize the survey data gathered by the City of Toronto regarding E-Bike usage, and to discern if there were any patterns in the responses which could predict a particular answer to the question "Does your household have access to any of the following types of private motorized vehicles?" 
The conclusions from this data analysis are presented in the Conclusions cell, at the bottom of the notebook.  My responses to the challenge question are presented therein.  I began by importing the tools needed to load the data and perform EDA:
End of explanation
"""


response_q = 'Does your household have access to any of the following private motorized vehicles?'
positive_label = 'No - I do not have access to a private motorized vehicle'

labels = survey_data[response_q]
predictors = survey_data.drop(columns=[response_q, 'Timestamp'])
missing_values_per_predictor = predictors.isna().sum()

print("Original data was shape: ", survey_data.shape)
print("Predictors have shape: ", predictors.shape)

print("Found ", labels.isna().sum(), " rows with response missing")
print("Found ", missing_values_per_predictor.sum(), " number of missing values total among all predictor columns")

"""
Explanation: First impression of the data and assumptions
Based on my first examination of the data, the responses in each column can be modeled as one of three types:
- ordinal; ordered string responses
- binary; yes / no, male / female responses
- categorical; un-ordered string responses
So the way forward to begin attacking the problem is to re-encode each column appropriately.  The appropriate encoder for each response type (in scikit-learn anyhow) is:
- CategoricalEncoder(encoding='ordinal') for ordinal data
- LabelEncoder(n=2) for binary data
- CategoricalEncoder(encoding='onehot') for categorical data
Except I don't have a version of scikit-learn that has ordinal encoding (it's in pre-release alpha version 0.20, I have 0.19.1)  So I'll assume that all non-binary features are categorical features, ignoring the ordinal nature of many of the responses.  
However I choose to represent the data, I should make sure to isolate the responses my model will try to predict.  Also, we want to drop any rows where the response is missing:
End of explanation
"""


predictors = predictors.fillna(value='NA')
print("Found ", len(set(labels)), "different responses")

positives = 0
negatives = 0
for label in labels:
    if label.strip().startswith(positive_label):
        positives += 1
    else:
        negatives += 1
print("Got ", positives, " matches to specific answer")
print("Got ", negatives, " non-mathces to specific answer")
print("Expected total of ", labels.shape[0], " responses, got ", positives + negatives)

"""
Explanation: Let's do a bit of simple cleaning.  We will fill missing values with 'NA', and make sure that when we break down the label distribution, we get the expected number of positive and negative cases that together form the total:
End of explanation
"""


question_response_complexity_dict = OrderedDict()
for k in predictors.keys():
    question_response_complexity_dict[k] = ""
    
print("Number of questions: ", len(predictors.keys()))  
print("Cardinality of responses to each survey question:\n")    
for k in question_response_complexity_dict.keys():
    data = predictors[k]
    cardinality = len(set(data))
    question_response_complexity_dict[k] = cardinality
    print(k, "\t", cardinality)


"""
Explanation: Initial assumptions about the data are wrong
We should pause here to reflect on a two points revealed by our EDA so far.


There are 123 different answers to the question we are trying to predict as response.  We cannot hope to predict 123 separate classes with only 2238 data points, so we will treat this as a binary classification problem.


This makes me want to check the variability in the answers to each other question.  I inferred from a scan of the first five records that the responses to each question were codified; that is, I could associate a categorical or ordinal code to each text response.  But counting the number of different responses to each question shows that this is not the case:
End of explanation
"""


%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-white')    
    
### look at hist of responses to each predictor
fig, ax = plt.subplots(5, 4, figsize=(14, 14))
fig.subplots_adjust(hspace=0.4, wspace=0.4)
r = 0
c = 0

for i,k in enumerate(predictors.keys()):
    data = [d for _, d in predictors[k].value_counts().items()]
    names = [n for n, _ in predictors[k].value_counts().items()]
    
    r = i % 5
    c = (i - r) % 4
    title_length = min(len(k.split()),8)
    title_line_one = ' '.join(k.split()[0 : title_length // 2])
    title_line_two = ' '.join(k.split()[title_length // 2 : title_length])
    title = title_line_one + "\n" + title_line_two
    ax[r,c].bar([i for i in range(len(data))], data, color='gray')
    ax[r,c].set_title(title)

"""
Explanation: Revising my assumptions of the data
It looks like only a small minority of questions have low cardinality: age, income, average commute distance, length of daily commute, aware of speed limits.  What about the distribution of responses to each question?  Let's look at a quick histogram of how frequently each response to the corresponding question appears in the data.  We can do this for each question using pd.Series.value_counts() and arranging the histograms into a grid:
End of explanation
"""


# Go through the predictors, find those answers in each question to be grouped into 'other'
# N.B: using a poorly chosen heuristic of 4 is not great, but based on my reading of
# the data manually, it seems an okay heuristic. 

# One of the ways to improve this would be to automatically determine an appropriate threshold via an information theoretic measure
# for each predictor, but this is out of scope for a first attempt.

replace_as_other = {k: [] for k in predictors.keys()}
for k in predictors.keys():
    replace_as_other[k] = [n for n, d in predictors[k].value_counts().items() if d < 4]

# perform the in-place replacement
for k in predictors.keys():
    predictors.loc[predictors[k].isin(replace_as_other[k]), k] = 'unusual answer'
    

"""
Explanation: Obervations about the data, deciding on an approach to predict the response
The histograms reveal a few things about the data:
1) Most of the questions have a few common responses (likely provided in the survey), followed by a long tail of principly unique responses that seem provided free-form. If I were to drop or group the free-form answers into one factor level (e.g 'other'), I could model most questions as categorical factors.
2) A second way would be to model each answer in the responses as a bag of words, and turn each answer in to a sparse bag of vectorized features.
Without wanting to get too deep into reverse engineering the survey, I'll try to model the collection of responses first as a set of categorical variables, and then later as a stemmed bag of words.
Modeling predictors as categorical variables
End of explanation
"""


%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-white')    
    
### look at hist of responses to each predictor
fig, ax = plt.subplots(5, 4, figsize=(14, 14))
fig.subplots_adjust(hspace=0.4, wspace=0.4)
r = 0
c = 0

for i,k in enumerate(predictors.keys()):
    data = [d for _, d in predictors[k].value_counts().items()]
    names = [n for n, _ in predictors[k].value_counts().items()]
    
    r = i % 5
    c = (i - r) % 4
    title_length = min(len(k.split()),8)
    title_line_one = ' '.join(k.split()[0 : title_length // 2])
    title_line_two = ' '.join(k.split()[title_length // 2 : title_length])
    title = title_line_one + "\n" + title_line_two
    ax[r,c].bar([i for i in range(len(data))], data, color='gray')
    ax[r,c].set_title(title)

"""
Explanation: Let's re-examine the distribution of answers to each question now:
End of explanation
"""


from sklearn.feature_extraction import DictVectorizer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import Imputer

# Encode labels in the response
label_encoder = {positive_label: 1}
encoded_labels = [label_encoder.get(r,0) for r in labels]
y = np.array(encoded_labels)

# Encode the features in the predictors, imputing missing values
encoder_pipeline = Pipeline([
    ("vectorizer", DictVectorizer(sparse=False)), 
    ("imputer", Imputer(missing_values=np.nan, 
                    strategy="mean",
                    axis=0))
                    ])

features = predictors.to_dict('records')
X = encoder_pipeline.fit_transform(features)

# Make sure we get sensible features
print("shape of X: ", X.shape)
print("type of X: ", type(X))
print("Any NaN or infinity? ", np.isnan(X).sum())

from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.linear_model import LogisticRegressionCV
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier

# CV on logistic regression
lr_clf = LogisticRegressionCV()
lr_clf.fit(X,y)

# CV on n_estimators, max_depth for both GBC, RFC
param_grid_dict = {'n_estimators': [5,10,20,50,100,200], 'max_depth': [2,4,6,8,10], 'min_samples_leaf': [1,3,5,10,15]}

rfc = RandomForestClassifier(n_estimators=10, max_depth=None, min_samples_split=2, random_state=0)
gbc = GradientBoostingClassifier(n_estimators=10, max_depth=None, min_samples_split=2, random_state=0)

random_forest_clf = GridSearchCV(estimator=rfc, param_grid=param_grid_dict, n_jobs=2)
random_forest_clf.fit(X, y) 

gbc_clf = GridSearchCV(estimator=gbc, param_grid=param_grid_dict, n_jobs=2)
gbc_clf.fit(X,y)

"""
Explanation: Great.  Now, we encode the answers to each question as categotical 'one-hot' answers, and train models to accurately predict our response. Also we encode the labels to be predicted.
End of explanation
"""


print("Model: Logistic Regression")
print("Model best params: ", lr_clf.C_)
print("Model best score: ", lr_clf.score(X,y))
print("\n")

for name, model in zip(["Random Forest", "Gradient Boosted Classifier"],[random_forest_clf, gbc_clf]):
    print("Model: ", name)
    print("Model best params: ", model.best_params_)
    print("Model best score: ", model.best_score_)
    print("\n")
 

"""
Explanation: Examining the model performances and the best features
Let's take a look at a crude breakdown of the best scoring models, and what hyperparameters define them:
End of explanation
"""


# Unpack the Logistic Regression classifier
best = lr_clf.coef_
feature_importance = 100.0 * (best / best.max())

# Retrieve the feature names
feature_names = encoder_pipeline.named_steps['vectorizer'].get_feature_names()
indices = np.argsort(np.abs(feature_importance))

top_n = 10
print("\nTop ", top_n, " feature ranking:\n")
for f in range(top_n):
    print(f + 1, feature_names[indices[0,f]], feature_importance[0,indices[0,f]])



fig, ax = plt.subplots(1, 2, figsize=(12, 2))
fig.subplots_adjust(hspace=0.4, wspace=0.4)
from sklearn.ensemble import partial_dependence

for i, (name, model) in enumerate(zip(["Random Forest", "Gradient Boosted Classifier"],[random_forest_clf, gbc_clf])):

    best = model.best_estimator_
    importances = best.feature_importances_
    #std = np.std([tree.feature_importances_ for tree in best.estimators_],
    #             axis=0)
    indices = np.argsort(importances)[::-1]

    # Print the feature ranking
    print("\nTop ", top_n, name, " feature ranking:")
    print("(indexed feature, question and response, feature weight) \n")
    for f in range(top_n):
        print("*",indices[f],"*", feature_names[indices[f]], importances[indices[f]])

    # Plot the feature ranking
    title = name + " feature importances"
    ax[i].bar(range(top_n), importances[indices[0:top_n]],
           color="r", align="center")
    ax[i].set_xticks(range(top_n))
    ax[i].set_xticklabels(indices[0:top_n])
    ax[i].set_title(title)



"""
Explanation: General comments
Since we aren't training this model to generalize to new data but rather to understand our survey data set, I cheated and did not follow best data hygene, which is to make sure training, test and validation data are strictly exclusive.  Despite this we can still learn a lot from examining the distribution of model scores.  The reported best score is accuracy over a held out fold of training data, which for binary classification is the proportion of correctly classified elements.  
Each classifier I trained performed about equally, with perhaps a minor difference between the random forests and the logistic regression classifier, but on a data set this small with the cross-validation scheme used, we cannot say with certainty.  Let's proceed to unpack the most important features by classifier.
End of explanation
"""


# shorten the names of the features, removing non-informative words
import nltk
nltk.download('stopwords')
from nltk.corpus import stopwords
stops = set(stopwords.words('english'))

def strip_to_three_max(q_or_a):
    ''' strip a multi word question or answer down to a max
    length of its first three words '''
    words = [word for word in q_or_a.split() if word not in stops]
    new_len = min(3, len(words))
    return " ".join([w for w in words[0:new_len]])

def repro(name):
    ''' shorten the feature name '''
    
    parts = name.strip().split("=")
    if len(parts) != 2:
        return name
    
    q, a = parts
    return strip_to_three_max(q) + " = " + strip_to_three_max(a)
    
    
reprocessed_features = [repro(f) for f in feature_names]

from sklearn.ensemble import partial_dependence as pdep

top_n = 5
model_name = "Gradient Boosted Regression Tree"
model = gbc_clf
    
best = model.best_estimator_
importances = best.feature_importances_
indices = np.argsort(importances)[::-1] 

top_n_features = [indices[f] for f in range(top_n)]

# Plot the feature ranking
title = name + " feature importances"
fig, axs = pdep.plot_partial_dependence(best, X, top_n_features, feature_names=reprocessed_features, n_cols=top_n)
fig.set_size_inches(16, 6)
fig.suptitle('Q: ' + response_q + '\n' + "A: " + positive_label)
plt.subplots_adjust(wspace=0.75)  


"""
Explanation: Feature importances
The plots above reveal the most important 10 features for the random forest model as well as the gradient boosted regression tree classifier model.  The overlap between both sets of features is large, which is encouraging.  However, it only identifies those features which are most informative for constructing the trees used to separate positive from negative responders, and not the importance of the directionality of the response.  
To find out which features drive a responsdent to answer affirmatively (i.e "No - I do not have access to a private motorized vehicle" to the question "'Does your household have access to any of the following private motorized vehicles?'"), we need to resort to a different analysis.  Fortunately, the gradient boosted regression tree classifier model allows us to efficiently calculate just this relationship, termed partial dependence.  We do so below for the top five features:
End of explanation
"""


