# -*- coding: utf-8 -*-
"""DMML Final Project.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/15spI7M-iNDjDFKOLMx5P2roJ7fFvytPS
"""

# Commented out IPython magic to ensure Python compatibility.
# Packages
import pprint
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from cv2 import equalizeHist
from skimage import filters, morphology
from sklearn.metrics import confusion_matrix, classification_report, recall_score, f1_score
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.preprocessing import normalize
from sklearn.feature_selection import chi2

# %matplotlib inline

pp = pprint.PrettyPrinter(indent=2)

"""# Importing the data

## Importing Functions
"""

def import_data(seed: object = 42) -> object:
    """ Import data from files.

        Note: The data must be in the root of the current working directory
    """

    # Read input data
    df = pd.read_csv("x_train_gr_smpl.csv").astype(int)

    # label data-frame rows based on sample data
    for x in range(10):
        index = ~pd.read_csv("y_train_smpl_%s.csv" % x, squeeze=True).astype(bool)  # reversed flags (~)
        df.loc[index, 'label'] = str(x)

    input_data_ordered = df.iloc[:, 0:2304].to_numpy()
    output_data_ordered = df.iloc[:, 2304].to_numpy()

    # Randomise instance order (forcing the same result each time)
    np.random.seed(seed)
    permutation = np.random.permutation(df.shape[0])

    # Create base input and output arrays
    input_data = input_data_ordered[permutation]
    output_data = output_data_ordered[permutation]

    return input_data, output_data, df, input_data_ordered, output_data_ordered

"""## Implementation

### Import the data
"""

input_data, output_data, df, input_data_ordered, output_data_ordered = import_data()

# Image labels with classes as dictionary
labels = {
    '0': 'Speed limit 60',
    '1': 'Speed limit 80',
    '2': 'Speed limit 80 lifter',
    '3': 'Right of way at crossing',
    '4': 'Right of way in general',
    '5': 'Give way',
    '6': 'Stop',
    '7': 'No speed limit general',
    '8': 'Turn right down',
    '9': 'Turn left down'
}

# Show first 5 rows to make sure data imported
df.head()

"""# Data analysis

## Display Functions
"""

import numpy as np

def image_as_square(image: np.ndarray):
    image_width = 48
    image_height = image_width

    return np.resize(image, (image_width, image_height)).astype(np.uint8)


def image_as_array(image: np.ndarray):
    image_width = 48
    image_height = image_width
    return np.resize(image, (image_width * image_height,)).astype(np.uint8)

def is_data_unique(data: np.array) -> bool:
    """
      Compare the shape of the data set to a version which only contains
      unique images.

      If the shapes are equal, then the data only consists of unique images.
    """

    data_shape = data.shape
    unique_data_shape = np.unique(data, axis=0).shape

    if data_shape == unique_data_shape:
        return True

    return False

def display_images_as_grid(images: np.ndarray):
    image_count = len(images)
    column_count = 8
    fig_width = 11

    # Calculate the number of rows needed
    row_count = image_count // column_count

    if image_count % column_count > 0:
        row_count += 1

    # Calulate the figure size so each image isn't tiny
    fig_height = fig_width / column_count * row_count

    # Setup the plot
    fig, axes = plt.subplots(
        nrows=row_count,
        ncols=column_count,
        figsize=(fig_width, fig_height)
    )

    ax = axes.ravel()

    #  Disable axis on all subplots, even those without images
    for i in range(0, len(ax)):
        ax[i].axis('off')

    # Resize images to shape of (48, 48)
    images = np.apply_along_axis(image_as_square, 1, images)

    # Plot the images
    for i in range(0, len(images)):
        image = images[i]
        ax[i].imshow(image, cmap='gray', vmin=0, vmax=255)
        ax[i].set_title(i, color='black')

    fig.tight_layout()

def plot_frequency_distribution(data):
    data_labels = [v for v in labels.values()]

    data = pd.DataFrame({
        'count': pd.Series([
            np.sum((data == '0')),
            np.sum((data == '1')),
            np.sum((data == '2')),
            np.sum((data == '3')),
            np.sum((data == '4')),
            np.sum((data == '5')),
            np.sum((data == '6')),
            np.sum((data == '7')),
            np.sum((data == '8')),
            np.sum((data == '9')),
        ]),
        'class': pd.Categorical(data_labels)
    })

    # Set style
    sns.set_style("white")
    sns.set_style("ticks")

    # Colour palette
    sns.color_palette("hls", 10)

    # Create new figure
    plt.subplots(figsize=(11, 3))
    # Create plot
    sns.barplot(y='class', x='count', data=data, palette='hls')
    plt.xticks(np.linspace(0, 2250, 10))
    # Remove spine
    sns.despine(bottom=True)
    #  Add grid lines
    plt.grid(b=True, which='major', axis='x')
    # Remove axis labels
    plt.ylabel("")
    plt.xlabel("")
    plt.show()

def compare_class_distribution_with_normalized(input_data, output_data):
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(11, 2))

    sns.set_style("white")
    sns.set_style("ticks")
    palette = sns.color_palette('hls', 10)

    input_data_normalized = np.zeros(input_data.shape)
    
    for label in range(10):
        x = input_data[output_data.astype(np.int) == label]
        normalized = normalize(x)
        input_data_normalized[output_data.astype(np.int) == label] = normalized
        sns.distplot(x.flatten(), ax=ax[0], hist=False, kde=True, color=palette[int(label)])
        sns.distplot(normalized.flatten(), ax=ax[1], hist=False, kde=True, color=palette[int(label)])

    sns.despine(ax=ax[0], left=True, trim=True)
    sns.despine(ax=ax[1], left=True, trim=True)

    ax[0].get_yaxis().set_ticks([])
    ax[1].get_yaxis().set_ticks([])

    ax[0].grid(b=True, which='major', axis='x')
    ax[1].grid(b=True, which='major', axis='x')

    fig.tight_layout()
    plt.show()

"""## Implementation

### Check instances within data
"""

def instance_count_per_class(data):
    unique, counts = np.unique(data, return_counts=True)
    count_dict = dict(zip(unique, counts))

    return dict((labels[key], value) for (key, value) in count_dict.items())


# Number of instances in the data set
instance_count = np.shape(input_data)[0]
print(f"Number of images: {instance_count}")

# Check if each instance (row/image) is unique
print(f"Is data unique? {is_data_unique(input_data)}")

# Print number of instances per class
pp.pprint(instance_count_per_class(output_data))

"""### Plot class frequency distribution"""

plot_frequency_distribution(output_data)

"""### Display the images as a grid"""

display_images_as_grid(input_data[0:200])

"""# Naive Bayes

Applying Naive Bayes
1. Split the data in 70/30 for train and test sets respectively
2. Apply Naive Bayes Gaussian variant
3. Predict classification using trained model
4. Display classification report and performance statistics

Notes:
- Using non-random data model predicts 955 out of 3798 signs
- With randomized data model predicts 1125 out of 3798 signs 
- With normalization model predicts 2067 out of 3798 signs
- When overtly enhancing the contrast, model predicts 3257 out of 3798 signs
- Selecting the top featurns improves the model to 3383 of 3798 signs

## Functions
"""

from sklearn.metrics import confusion_matrix

def plot_confusion_matrix(y_true, y_pred):

  title = 'Normalized confusion matrix'
  classes = [v for v in labels.values()]

  cm = confusion_matrix(y_true, y_pred)
  cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
  cm_df = pd.DataFrame(cm, index=classes, columns=classes)

  flg, ax = plt.subplots(figsize=(12, 8))
  ax = sns.heatmap(cm_df, annot=True, cmap="Blues")
  bottom, top = ax.get_ylim()
  ax.set_ylim(bottom + 0.5, top - 0.5)
  ax.set_title(f'{title}')
  ax.set_ylabel('True label')
  ax.set_xlabel('Predicted label')
  plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")

  plt.show()

# Commented out IPython magic to ensure Python compatibility.
from sklearn.metrics import confusion_matrix, classification_report, recall_score, f1_score

def print_outcomes(x_test, y_test, y_predictions):
  cm = pd.DataFrame(data=confusion_matrix(y_test, y_predictions), index=labels, columns=labels)
  cr = classification_report(y_test, y_predictions, target_names=list(labels.values()))
  rc = recall_score(y_test, y_predictions, average='weighted')
  f1 = f1_score(y_test, y_predictions, average='weighted')

  print("Outcome:\nNaive Bayes correctly predicted the classification for %s of %s signs using %s attributes\n" 
#         % (sum(y_predictions == y_test), len(y_test), x_test.shape[1]))
  print("Classification report:\n%s\n" % cr)
  print("Confusion matrix:\n%s\n" % cm)
  print("Weighted Recall: %.3f" % rc)
  print("Weighted F1: %.3f\n" % f1)

"""## Running Naive Bayes with different filters

#### Without pre-processing
"""

from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# split the data into training and test sets for each class
x_train = []
x_test = []
y_train = []
y_test = []

for label in range(10):
    x = input_data_ordered[output_data_ordered.astype(np.int)==label]
    y = output_data_ordered[output_data_ordered.astype(np.int)==label]
    x_train_tmp, x_test_tmp, y_train_tmp, y_test_tmp = train_test_split(x, y, test_size=0.3, shuffle=False)
    x_train.extend(x_train_tmp)
    x_test.extend(x_test_tmp)
    y_train.extend(y_train_tmp)
    y_test.extend(y_test_tmp)
    

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(np.array(x_train), np.array(y_train))
y_predictions = model.predict(np.array(x_test))

print_outcomes(np.array(x_test),np.array(y_test), y_predictions)
plot_confusion_matrix(np.array(y_test), y_predictions)

"""#### With Randomizing"""

from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# Run naive bayes on randomized data set
x = input_data
y = output_data

# split the data into training and test sets using stratified random sampling
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(x_train, y_train)
y_predictions = model.predict(x_test)

print_outcomes(x_test, y_test, y_predictions)
plot_confusion_matrix(y_test, y_predictions)

"""#### With Normalization"""

from sklearn.preprocessing import normalize
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# Normalize the input data
input_data_normalized = np.zeros(input_data.shape)

for label in range(10):
    x = input_data[output_data.astype(np.int)==label]
    normalized = normalize(x)
    input_data_normalized[output_data.astype(np.int)==label] = normalized

compare_class_distribution_with_normalized(input_data, output_data)

# Run naive bayes on data without image preprocessing
x = input_data_normalized
y = output_data

# split the data into training and test sets using stratified random sampling
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(x_train, y_train)
y_predictions = model.predict(x_test)

print_outcomes(x_test, y_test, y_predictions)
plot_confusion_matrix(y_test, y_predictions)

"""### Interesting find (found too late)"""

# *** late interesting find resulting in f1-score of 0.97

from sklearn.preprocessing import PowerTransformer

pt = PowerTransformer(method='box-cox', standardize=False)
input_data_transformed = np.zeros(input_data.shape)

for label in range(10):
    x = input_data[output_data.astype(np.int)==label]
    transformed = pt.fit_transform(x)
    input_data_transformed[output_data.astype(np.int)==label] = transformed
    sns.set_style('darkgrid')
    sns.distplot(transformed.flatten(), hist=False, norm_hist=True)

"""## Feature selection

Using Chi2 to determine feature importance for each of the labels of enhanced and normalized datasets
"""

from sklearn.feature_selection import chi2

def display_image_from_array(image: np.array, title):
    # resize the image to 48 x 48 shape
    image = image_as_square(image)

    # Display the image
    plt.axis('off')
    plt.tight_layout()
    plt.imshow(image, cmap='gray')
    plt.title(title)
    plt.show()

def vizualize_important_features(input_data, output_data):
    X = input_data
    Y = output_data
    feature_importances = chi2(X=X, y=Y)[0]

    N, M = X.shape

    out = {}
    for c in set(Y):
        out[c] = dict(zip(range(N), np.mean(X[Y==c, :], axis=0)*feature_importances))

    for item in out.items():
        # this transformation is only done for visualisation purposes
        scores = np.fromiter(item[1].values(), dtype=float)
        scores *= 255.0/max(scores)
        display_image_from_array(image=scores, title=f"Label: {item[0]}, Sign: {labels[item[0]]}")

"""### Visualise important features for normalized data"""

from sklearn.feature_selection import chi2

X = input_data_normalized
Y = output_data
feature_importances = chi2(X=X, y=Y)[0]

N, M = X.shape

out = {}
for c in set(Y):
    out[c] = dict(zip(range(N), np.mean(X[Y==c, :], axis=0)*feature_importances))

for item in out.items():
    # print("\nLabel: %s, Sign: %s" % (item[0], labels[item[0]]))
    # this transformation is only done for visualisation purposes
    scores = np.fromiter(item[1].values(), dtype=float)
    scores *= 255.0/max(scores)
    display_image_from_array(image=scores, title=f"Label: {item[0]}, Sign: {labels[item[0]]}")

"""### Top features for each class"""

from heapq import nlargest


def get_n_best_fetures(X, Y, N, show_: bool = False, print_: bool = False):

    idx = []
    for item in out.items():
        top = nlargest(N, item[1], key=item[1].get)
        idx.extend(top)

        if print_:
            print("\nLabel: %s, Sign: %s" % (item[0], labels[item[0]]))
            print("Selected %s features: %s" % (len(top), top))

    selected_features = list(set(idx))

    if print_:
        print("\n%s selected feature indexes:\n%s" % (len(selected_features), selected_features))

    if show_:
        all_pixels = np.repeat(0, X.shape[1])
        np.put(all_pixels, selected_features, 255)
        display_image_from_array(all_pixels, title="")

    return X[:, selected_features]



from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

def run_naive_bayes(x, y, print_: bool = False):

    # split the data into training and test sets using stratified random sampling
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

    # train data sets and use model for predictions 
    model = GaussianNB()
    model.fit(x_train, y_train)
    y_predictions = model.predict(x_test)

    f1 = f1_score(y_test, y_predictions, average='weighted')
    
    if print_:
        print("\nUsing %s best attributes, f1 score: %.3f\n" % (x.shape[1], f1))
        print_outcomes(x_test, y_test, y_predictions)
        plot_confusion_matrix(y_test, y_predictions)

    return f1

"""### 2 best attributes for each class

Selecting top 2 features for each class results in total of 9 pixels selected for the model due to overlap of best features between classes. Seemingly random features perform suprisingly good predicting 2481 of 3798 signs (65.3%) in the test data.
"""

input_data_2_best  = get_n_best_fetures(input_data_normalized, output_data, 2, True, True)
run_naive_bayes(input_data_2_best, output_data, True)

"""### 5 best attributes for each class

Increasing the number of best attributes to 5 per each class results in a total of 26 attributes selected across all classes. The performance improves to correctly predicting 2868 of 3798 signs (75.5%).
"""

input_data_5_best  = get_n_best_fetures(input_data_normalized, output_data, 5, True, True)
run_naive_bayes(input_data_5_best, output_data, True)

"""### 10 best attributes for each class

Selecting 10 attributes for each of the classes results in a total of 43 attributes. The performance is slightly improved and the model predicts correctly 3021 of 3798 signs (79.5%).
"""

input_data_10_best = get_n_best_fetures(input_data_normalized, output_data, 10, True, True)
run_naive_bayes(input_data_10_best, output_data, True)

"""### Best number of attributes for each class

Looking for the overall best number of attributes to be selected for each class the model was run on a range of 1 to 2000 best attributes. This is considered a slow method as model needs to be built and assesed for each of the 2000 iterations.

As a result, the best overall balanced accuracy score of 0.891 indicated 136 as the best number of attributes to be selected for each class resulting in a total of 376 attributes. The performance seems to be best around 100 best pixels and drops ever so slightly when increased up to the maxiumum available pixels. This shows that filtering for the most informative attributes, in this case pixels, can not only improve performance and reduce the size of dataset processed and stored, but also is expected to have a positive impact on overal predictive power of a naive bayes model.
"""

x_array = []
y_array = []
best = [0, 0]
print("Looking for best number of attributes for each class...")

for idx in range(1, 2000, 1):
    x = get_n_best_fetures(input_data_normalized, output_data, idx)
    f1 = run_naive_bayes(x, output_data)
    x_array.append(idx)
    y_array.append(f1)
    if f1 > best[1]:
        best = [idx, f1]

print("Best number of attributes: %s, F1 score: %s\n" % (best[0], best[1]))

plt.plot(x_array, y_array)
plt.show()

input_data_overall_best = get_n_best_fetures(input_data_normalized, output_data, 3, True)
run_naive_bayes(input_data_overall_best, output_data, True)

"""### KFold cross validation

Data was split into 5 equally sized bins to perform kfold cross validation on each of the combinations of 4 train and 1 test set - total of 5 runs to validate the model and best number of attributes selected in previous step. The resulted balanced accuracy scores are comparable which indicates a high validity of the selected model and assumptions made.
"""

# Commented out IPython magic to ensure Python compatibility.
from sklearn.model_selection import StratifiedKFold

x = input_data_overall_best
y = output_data

skf = StratifiedKFold(n_splits=5)

for train, test in skf.split(x, y):
    f1 = run_naive_bayes(x[train], y[train], False)
    print("Train: %s, Test: %s, F1: %.3f" 
#           % (len(x[train]), len(x[test]), f1))

"""# Complex Bayesian Networks

1. Data Transformations to csv
2. Running various classifiers in Weka
3. Comparing the results based on f1-score (visualise with seaborn)
"""

# Stratify function
from sklearn.model_selection import train_test_split

def stratify(x, y, filename, selected_features):
    # split the data into training and test sets using stratified random sampling
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)
    
    # make them dataframes
    df_x70 = pd.DataFrame(x_train)
    df_x30 = pd.DataFrame(x_test)
    df_y70 = pd.DataFrame(y_train)
    df_y30 = pd.DataFrame(y_test)

    # Join them all up
    df_x = df_x70.append(df_x30)
    df_y = df_y70.append(df_y30)
    df_y.columns = ['Class']
    df_x.columns = selected_features

    df_stratified = pd.concat([df_x, df_y], axis=1, ignore_index=True)

    # export to csv
    df_stratified.to_csv(f"Weka_{filename}.csv", index=False)

# build csv files for Weka

# FULL DATASET:
stratify(input_data, output_data, "Full_Original")
stratify(enhanced_input_data_hist, output_data, "Full_Enhanced")

# BEST_2:
original_input_data_2_best = get_n_best_fetures(input_data, output_data, 2, True)
enhanced_input_data_2_best = get_n_best_fetures(enhanced_input_data_hist, output_data, 2, True)
stratify(original_input_data_2_best, output_data, "2_Best_Original")
stratify(enhanced_input_data_2_best, output_data, "2_Best_Enhanced")

# BEST_5:
original_input_data_5_best = get_n_best_fetures(input_data, output_data, 5, True)
enhanced_input_data_5_best = get_n_best_fetures(enhanced_input_data_hist, output_data, 5, True)
stratify(original_input_data_5_best, output_data, "5_Best_Original")
stratify(enhanced_input_data_5_best, output_data, "5_Best_Enhanced")

# BEST_10:
original_input_data_10_best = get_n_best_fetures(input_data, output_data, 10, True)
enhanced_input_data_10_best = get_n_best_fetures(enhanced_input_data_hist, output_data, 10, True)
stratify(original_input_data_10_best, output_data, "10_Best_Original")
stratify(enhanced_input_data_10_best, output_data, "10_Best_Enhanced")

"""# Clustering

## Evaluation of clustering algorithms

**Evaluation of Clustering Algorithms**

The evaluation of clustering algorithms is quite different from normal cases (of supervised learning) where there is a known (confirmed) value associated with the data passed. This is not always the case with clustering. Two classes of evaluation metrics exist depending on weather the 'class' (or label) associated with data is known or not. 

(source: https://scikit-learn.org/stable/modules/clustering.html#clustering-performance-evaluation)

1.   Ground truth dependent metrics

*   **Adjusted Rand Index** ('ari'): This is a measure of the similarity between the actual labels assigned to the instances and to the predicted labels. Are all the instances in one present in the other. It is symmetric in the sense that order of passing variables does not matter to the function. 1.0 is considered the perfect score that all the instances belonging to one actual label are predicted to be in the same cluster. A score 0.0 implies there is no correlation between the actual and predicted labels, the clustering is more or less random. 

*   **Adjusted mutual Info Score** ('ami'), **normalized mutual info score** (nmi): Mutual Information is a function that measures the agreement of the two assignments, ignoring permutations. NMI is the more prominent one present in literature, AMI was more recently proposed. Mutual Information as a mathematical concept is a function of the entropy of the class. A higher value (bound by 1.0) means more similarity of distribution of instances between actual and predicted labels, a lower value indicates that the instances are spread differently in predicted labels from actual values.

*   **homogenity**: Each cluster contains only members of a single class. (1.0 indicates completely homogenus class with only members of the same class, 0.0 means predicted label contains members from different actual classes). 

*   **completeness**: All the members of the actual class belong to the same cluster (1.0 for proper completeness i.e the statement is true for all classes. 0.0 implies all classes are spread through out different clusters). 

2.   Ground truth independent metrics:

*   **Silhouette Coefficient** ('sc'): A measure of how well defined the clusters are. A high score (upper bound 1.0) implies that the clusters are well separated. A score of 0.0 implies possible overlap between clusters. A further lower score (bound -1.0) implies poor clustering. The Silhouette Coefficient is defined for each sample and is composed of two scores: **a**: The mean distance between one instance and all other instances of the same class. **b**: The mean distance between one instance and all other instances in the nearest neighbouring cluster. 

![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAQwAAABuCAYAAADSxy12AAAflklEQVR4Ae1dB3RURRe+tNAVlK4QIAkQOoEkQOhNeofQe1M6oiiCglhpBkQFREURCU26gkgNEJLQISQhEAgl1ICA4A8J7H++SWbysnmb7Ibspuydc/a8eTN32vdm7t65c2cmm8FgMBA7RoARYATMQCC7GTRMwggwAoyAQIAZBncERoARMBsBZhhmQ8WEjAAjwAyD+wAjwAiYjQAzDLOhYkJGgBFghsF9gBFgBMxGgBmG2VAxISPACDDD4D7ACDACZiPADMNsqJiQEWAEmGFwH2AEGAGzEWCGYTZUTMgIMALMMLgPMAKMgNkIMMMwGyomZAQYAWYY3AcYAUbAbASYYZgNFRMyAowAMwzuA4wAI2A2AjnNpmRCu0AA5ylpz1R6/vw55czJ3cQuPr4ZjeSeYAZI9kDy7Nkzmv+VDz18+JCuXYuiBw8eULZs2Sg6Opp2/rWdsmdnYdQe+kFKbWSGkRJCdhIfdi6cvl70bZLWNmncKEkYB9gvAvy3Yb/fPlHLK7tWovBzIbRr5w4a/dYoFefs7MzShUKDPcwwuA8oBHI7OJCTU3k6cuSoCvPyqqf87GEEmGFwH0iEAJScZ4LPijCHXLnIza12onh+sW8EmGHY9/dP0voLFyLo0aNHIrxuXU8qWCB/EhoOsF8EmGHY77fXbXnQkSMqHAyDl1QVHOwhImYY3A0UApiOHD4coN696tdXfvYwAkCAl1W5HygEwDBOnTot3l1dK1G1alUJYbDRwPV4ObJnF7YZWc0mA4ZqsbGxot2wPUH7slob1Ud+QQ8zjBcEMCslDwkNpUuRkaJJHu51KCAgkH748Sc6GxJCd6PvUunSr1PzFi1o8qQJ5ODgkCWa/vjxY9qwcROdPnWaDh7yp9y5HcjLy4vGjnmLChUuTGFhYVS4cGEqVbJklmjvizYiG9+t+qIQZp30P/y4nGZ+PEs0qFKlihQaGkbdu3Wlxo0a0uP//qO9e/fRn9t3ULu2beibRQsz9b8wpIrw8PM0/aMZ5O9/mPr26U3VqlWj8PBztGnzFnJ2ciIPD3da+PU3Ao8L4aGUK1eurPOxU9kSljBSCVxWS4aph7+/v2oWmMVCn/nUsWMHxRggroNhbPvjT2q6/nfq2aO7os9MHjCL/fv9aOq06XTlytUk7ezQvj31HziYDgcEimbV9fRQGGSmdlqjrqz0tAaqmTDP27fvCJFcVh3MonPnTokGStMmTWQ0rV69Rug3VEAm8kCiAEMAs4CuRssU0Qw3t1o0eNAA1aIyZcokwkFF2KGHGYYdfnS9JmM5VdpfdO3Smdq3b5eErEiRV1XY9es3MiXDgHJzydJlqh3jxo1NwgwggVy9ek3RODs7CWWvCrBjD09J7Pjja5seEhKiXiu7uuraX2C1RLqr165RbOwzXTpJo30+jYmh7NmyaYNe2J8jRw6LBjIYwS8rfqU9e/eKstu0foPatWmdpB6Ynp08eUqF16nD1q4SDGYYEgk7fuJfNzAwSCHg6emu/FoPVk2kK1vWkXLkME9AffDwIXl796bgswlMSebzIs85s78g7549zM7i0ePHtFQjXVSvVk037dOnT+lCRISIK1a0KNWqWVOXzh4DmWHY41c3avMzzf6R/Pnzk6urqxFF3Oup03E2Gngr6+ho9r97gfz5acq77xAkgrR09etbtjEuMjKSoq5fV1WoWbOG8ms9Z84Eq1d39zpmt1MlysIeZhhZ+OOa27ToO3eU/qJqlcq60wyI6Vor0FKlSpmbvdARNGnS2Gx6axFevBRnY4L8YYRWq5a+5HDs+HFVBUxH2IhLwcGm4QlQ2K9Pq+CrWrWK7gCJiIigffv9FEj4d85M+0yePzfQ9agoVf8aNapT3rx51bv0gDFqp2ewxYCD/kNag0pae3yaNwm1R2TsqM0YJNKVK1dOehM9d++JUxQisJd3T+rRvVui+Iz+An0rdBPSgeHBrsTY3X/wgAIC43Q1FSq4EBTAcOs3bKTyzhXtnmkwwzDuMXb4XqxYUdVqxzJllF96YmJiyc/vgHit4OJMkyaOT3N9hCzLWk8whyKvJiwLFyhQQLeoAwcO0oMHD0WcW61aStrasnkLde7UUTeNPQUyw7Cnr22irUWKFFExp0+fUX54IIr7LFggpiM5c+agCRPGUfHixRPRZJaXEiVKqKqiXcYOe2Z8fBaoYBhwgdFcvHiJ9uzdR61atshU0zDVkDT0MMNIQzAza1YFChSkwYMGiurv8/Mj2EzAYaoy/cOPxOHAMNr6fsliatumja4onxna3qCBl5hOoa7YaHb9xg3VTpi7jx07npo3b6aaki9fPoHBT8uXU/16dallyxYqzl49vPnMXr+8Ubuh0Bs9ZpzYK4K9E56eHnTuXLh4dypfjmZ/+QXVru2mRHSj5JnmNSYmhrx79xXnlhYsWJCGDB5IwcFn6e9du2nUyBE0aeIE+vyLL+in5b8QTkwvWrQorV23ntat9SX3OrzEygwj03R161cUlpzY6r3f7wA9uH+fsJJQo0YNcqtVkwoVKmT9CtioBDBHWHueOHGSQkJCBSOE7QkYBJZQYcG6Zu1aOnbsBF2+HEkdOrSnPr17ZTq9jTXgZIZhDVQzeZ6YimCOj/l7VrZBkKbuaKPeiolcPcrKGFjaVZlhWIoY0zMCdowAKz3t+ONz0xkBSxFghmEpYkzPCNgxAsww7Pjjc9MZAUsRYIZhKWJMzwjYMQLMMOz443PTGQFLEUjX7e1YD4cpcmhYGF29epVKlixJr7/2GjVs2IDXvC39kkzPCNgAgXRjGLC4m/bhDFq1ylfsTShX1lEc+x599y7hTMl5c2cz07BBB+AiGAFLEEgXhgFmMWbseGF2PH3aVBoyeJAwnImJfUZHgoKod9/+gmk0atTQkrYIWkgtekY4Fmekk8CUgY8OKQcxAlkSAZszDFjXSWbRsmVzGjZ0iBrgBw7sE5udgPQ3334npiaWDH5Y5n3ls4CCghIuFE7Lr4bLbjrxFue0hJTzymQI2Jxh3LlzR0gWwGlAv36KWeD92LFj6nj3/v37WQwlmEvlypUtSIcDVJJuczaVQW0+PdoUNBxuJwjYnGEcOpRwu9a69evJ2cWZSpYoIfYsTJwwnsLCzpGjYxlq4FU/ETMx53uAYbRt01r8zKEX+yWwj8AcYqZhBBgBsvlekr92/k3Dho9U0IM5NGrYkOrVq0te9esRthxbet+EyiyDe2Z+/InYVp3Bq8nVy2QIvFrkVVr+Y8LlTNasvs0ljKZNGhPutLgUf4JzZORlWhG5klb8upJwPNy7704WEkJaH0lvTRAtyfvkqYQLcixJx7SMgCkEmjVraioqzcNtLmFgGrB79x5xES6u29Nz8+fNoW5du1g8JdHLKyOF7fx7V0aqDtclCyHQskVzm7TG5gwDrQLTOH3mjDjO/cyZs3T02FGCpCFdq5Ytacnib1Jlh4FVGORvDYdlVT4bwRrIcp6ZBQGbMYx79+6R38FDdDkyUlw95+VVX2AkDynxXb2G3nv/AxFWqNDLdDQogHLlymURjshr8jtTaN363y1KZy7x5599QlhaZccI2CsCNmEY+Nfv0KkLySvoZn08kwb075toygGpwKtBY8Ilv+3btaUFPvMtZhj4iMFnz9Ld6LtW+Z5QzGamy3usAgJnatcI2ETpef7CBcUsgHZZxzKJmAXCnj17Tq+9VkowDBxCa6l0Ib9iFYvsMGQqfjICjIA5CNhkt+qrr7yi6oLbtvUu0fVdvZoCAoMIV/V5e/dU9OxhBBiBjIOATRjGK6+8QgMHxFluQs9w9Ogxun37Nv3zzz906tRpWrhwEU39YDrVdnOj9evWUO7cuTMOQlwTRoARUAjYRIeB0rApzGfBQlq9ei3dvHVLVQC3aTk6OopbpSa/PSnVUxGVIXuyNALYuPjVgoVUrGhRGjRwQJZqK/5MZ8+dRx3at6OMOrW2GcPAlwUg//33H926dUswjcKFClOZMqUFk+AlyyzV963SGCjPsZKG7QVfL/QhXGWYlRwU/x9+NIN+/uVXWrtmlbg4KaMt49uUYWSlj8ttsS0CYBYTJr5NmzZvIZ/5c6lr1y62rYAFpaGuhw4dphs3b9ClS5dUyrqenmIHtgrQ8cD8oFeffuKCpV1/7yAXZ2cdqvQLsskqSfo1j0vOCghAMl28ZKlgFv36ZvwjBuZ/5SPuozXGHgwjJVe4cGEaP24sjXpzNA0fPpLWr19L2kWDlNJbO94mSk9rN4Lzz9oIHDlylL6cPVfoLYYPG5oqC2BbIlSpYkUaN3Y0jRn9ZqJiXVxcEr2besGO60ED+1PExUs0Zcr7BIklozhmGBnlS3A9dBGAdPHzLytE3JAhg6hcuXK6dBkpEHexQoHfrl1bVa0ypUtT0aJF1HtKHhz14OpaibC7e+++fSmR2yyeGYbNoOaCUoMABsyWrdvEGSkjRwxPTRbplubokaOqbJzvYokCE1MT7OyGW79+Q4aRMphhqE/KnoyGAJbi//jjT1Gtxo0aZvipiBa/uA2WwSrIq4FXEutmFWnC07BBAxGzddsfVjt20kTRJoOZYZiEhiPSG4Hz5y/Q1m3bRDUaN26S3tWxqPznzw106vRplSY1S8CwiMYZMXBr1/0uzBJUhunkYYaRTsBzsSkjsGXrVoqNfUYFCxQQRzamnCLjUNy//49YGkWNYGRWonhxiyuHIycbNPAS6fb77aenT59anEdaJ2CGkdaIpiI/KPbkTy+5jMMzOY05RHhJm9ozQZBOlpPaPFAP+UuuPcnlj3aGh58XyWvWrEF58uTRy8qssBfFxKxCjIj8DweokFq1aorpiGwv6iP9isiEBzu34W7evJUhjndkOwwTH8oWwTjmBze/7du3T+yvcXZyojZtWlPt2m6ieAya7Tt20sWICAoJDRUWsqBxq+0m7m2RxxhiS/+JEycpMjJSHKKMvTulSpUkKAlfeukls5qCTowpAA42OnHiBN29e49w9Funjh2Esg7xcMlt70d99+zdR2FhYaI92H1crWpV6ty5Mzk45BJMZM3adXT8+AnKmzcPNWvaVPyD6uWJAXUmOE4H4OHhTtksPKkZ6XGiG9pz+vRpCggIjCuvoRdVrVJFMF7QpHZXdHKgIt+AgASGASkhMDCIIDFduXKV8ubNS6Vff53Gjh0tvk9yV2nUretJDg4OQrrA5kxMUyxRniZXz9TEsYSRGtTSKM3Zs2dp4sRJYus/NtwtXfYDdenWQ2zIu3v3Lr3z7nv05lujaemyZZQ3X1568vQprVzlS29Pfpd+/XWl+Jfy8ztAI0aMovenTqNdu3aL3b679+wVhkO9evelmzdvplhbSAOfff4ltWjVmpZ+v4wePXpM0dHRNHHSZJoxcxZt2LCRWrdpR9179hKDXi9D5IE7YYYMHS6YBZYE/f0P0+R336P33p9K//vf/2jKe1Pps8+/oGLFitLZsyE0ZNgIGj1mnK7UdOPmTXXlBK7QtMSBcf22ypfq1m9Aw0eMEmV5eHrQom+/o5Gj3qJtf/xJ3Xp4U+OmzWnjxk2WZG0WLZgrBrd0OMe1Z68+dPt2NNWoUV1cs4FvPXDQUAo7d06S6T7BHEqUiJvO7N+/X5fGpoEGdumCQExMjGHY8JGGfgMGGWJiYw1RUdcNpR3Li9/sOfMMHTt1MTRr3tKwdt16Q2zsM1HH2NhYQ/cevQRN8xatDBs3bjJUrFTFMHvOXENAQKDh+fPngu7UqdMqrw+mTTc8exaXXq+hT58+Nbzz7hRBP2jwUEN0dLTKJyAwUOUj6/a/J0+SZIN6/fDjT4L2t1W+qrzQ0DCDS0VXET5v/leGOu51DSdPnjSgzE8+/VzlfSY4OEmeO3f+reL9/A4kiTcVgLp8t3iJShsREaHqc/nKFUOr1m0NTs4VVbyv72pTWaU6HBhKvPB0q+1h8PVdY0Dd4PA9mjRtLmjGjB2nwk0V2L2Ht8rviQ7+ptJZI5wlDJuy54TCjp84QTv+2kktmzennDly0K1bCZLA14u+oX///ZcWLVpI3bt1pRw54j4T/m3c3WuLTM6Fn6dlP/wo9lS8PWkixYntcXJ75cquqqB9+/1MzpchFYwbP4F8V68V9MOHDyVMZ6SI7OHuTmPHvKXyqlmjuqirCog/n3XLlq1CEnlz1Ajq3ctbicwuLs5UwaWCIPdZ8DX17NmDqlevTpCAliz9XmWTL29e5ZcerWSE6ZU5Dv/sMMuGtAQHa0sYekkRXkwDxoympzExKjtrnLh9WKO/qODiTMt/+oG8vXuoZWHURxp1bdq8lfbu81P10fNICQNxDx8+1COxWRgzDJtBnbig4DPBQvvfqlVLEfHgQUJHcHF2EszCtVKlRIkwIC5ciFBhZcuWpU9mzVQDQkZcv35deunJkye6DAN5QWzf9sd2QTtwQH+qX6+eSic9VatWlV6xcUoOPhkI8f/LOXOpWtUqNGH8OBksnog7fz5OcYmATp06iHBMd6QDEypdurR8FU8wMu0RCOZOSY4dP05Ll8bdz4Ep0aSJExLlixd3ze11ONkNDDItHfQX0ClJ17r1G1StWgKGMlwbFhoaovuNJK22/Q8f/iuD0+XJSs90gB2dCvNoHISM+TzeD2uUZJMnv03GzEJWU56Lmj9/fpoy5Z0kzAJ0UVEJDKNWzTgNvUwvn9euRanB9VLBgjRs2BAZpZ6oF/QQ0kEZK6UPGQZlHvKaNvV9ocyT4fL5+L//hBeHI1WI30vRo3s3IalEXr5M/Qf0T6JIBVPCqgAcNl5B6ZeSA5NZsuR7oecBLZiXMXND+IGDh1RWTk5OujSKIBUeMGLolaSrp8OEESdXgOA/duy46APG2Mo8SpQoIb307yNmGAoMe/GgY/iuWimai06Nf+IDBxI6cp34VRJjPC5fvkJXrl4VwfhH1/7zaGkDAgPVKwa58cBBp8buz8tXrgi6OnVqE/Y6GDvQHTkaZ96MVQ33OnWMScSNdRHnw5IMehBKK034PT3dVVqsTGB6kpyTEoZ2sCRHv2bNWjHFA81LLxWkhvH2C8Zp/DVXdWLKZGqQGqcz9x2nyZ04GXdZVb58+QgSlLEDIz4ZT4M47DFJrh64SlS69J6SsIQhv4SNn3JJFMVei4oi6DTgcKapKTH5yJGEW+mrVKlM2XXWGtEZ9+9PmBNDijHujGBQu/fsUS12d6+ThAaR+BeUEk2dOnUoj46uAXmbWhbduy9Bq6833VEVMPIgz4RV1JTvmAFj00oOnh4eBAnM2KHd/ocTJCYsUae1C9dMwcCssYRq7FAPLKdLV7VK1WQlDC0C2bOlrxYhfUuXiNn5MygoYQkOR7MZSwSAB4zgyNFjCikwFmNGgMjw8HC1pAdJpWKFOKWjSkgk9CCYRkgHhameC9RIKu61aydReOqlkWGYIhw9FlffkiVLUKNGDWWUWU8pWUSZuB1PmwkGoNYMG8paPWxu3b5NV65eE0lz5cxJnp4e2mxe2I9vBBsU6YCZ3re8ePEiRWn0TNC36NHJfKKiEr5VgYIFZHC6PJlhpAvsCYXi3/Ho0eMqAFMNvc6OQXEsfgCCuFrVaiqN1oMVCOm00gXKQYeGOxOc8O8G8d2tVpyhmEyHJ2gDgxIkGq0UIvPR0hv7MYDl/bnQX+g55GMqL7kygBOocKxjcg42G5iuSWeKAYaEhEoSwjRM6kaSq4dKYIYH+WgZcas34hTaxknllAXhYKYpnZOhvVI0f758xtnZ9J0Zhk3hTloYOtnx4wkMo0aNGkmJiIQxV9i5cBGHm+HKlSubhA5MYdfuhKlG82bNxFIeyvDu1YfmzJ0nBiisDaWD+J49e8IEQIYjL6m/QJj8N4bFZE/vPsrYCnk/fvxYKPru378vkyfaXYmNV1it0TowwI9mzBSnxWvD4QfDlBIG3q9plLjGtHi/Fa8ghR8MULsCoaUPCQlRr5hioRy0E9jMnTdfxaXWg7wePHggkuMb6U15QANFsXSDBw6kgilIDVoJ4+VChWTSdHkyw0gX2BMKvXMnmoLPxnVkR8cyQoeREJvg2++XoJdo1LChrt4Ay4oY0HBVKruqgQOTcVgeuji7iEGCzVDSQRzWk2iCgo6o1RbUS4rMuIYSSlVcPAW3ync1tW3Xkfr2H0jDho8UAxBM5KBmNaJG9epJro7YufNvWv7zCqpfr66sSqJnKY11Z1RU3DQiEYHmxbWyKznEX6tpakqHKRKsS6WrVzfuuLzz5yPocEAgaZePQYOzOLHsvPWPP+nJE/M2fQEj7BuBK+voqOwuZJl4wix+9Zo4uxfQYnVKD39tmowkYbDSU/tl0sGvXdHAmQ9yYGqrYjwAGzdupNvJtJp3GCRBsYq0Py7/mTAA27VrI7ItVqyYyj5a51pJ/Pt/v+wHRYPlUNQL/46wMcAx+DAmwzsGQMTFi4IWS6go7+rVqwSDMenKlnWUXvHErsvvFi+h6tWqir0ziSLjX7DhTDrtMrEM0z5zOzgIRnvs+Am6c+eONkr4UScognEQD1zxYsWUxLR9x3ZyrVSRmjVN2D6Pdn3z7WI1sKEv2rThd8qVK/nhAozk9EKu8mgrA1xhbCfdkEEDdRm/jMcTdbl+I26ZHP1DT8Gspbe2nyUMayOcTP7oyFqFZ6NG+owAneZSZKTICaKuPInJOOvCGnFVnjZ98uRJWr78F3pvyjuqs4GZ9O3dSyQPDg5W0wsE4J+4S9futPPvXSp7OeAPHvIXqyYdO7RXzOjRo0eKbsSwoYKxyH9QGaGds0MygSXr8RMnadq0qapOklY+HRxyU614pqEVyWW89omB2ju+PeHnLySSboDdT8t/pkFDhlHzZk1FMugN8K+OAQx7mBYtmifahIYt9do2YKXIZ8ECbZEm/Tj0BhIZpAIthqjHgoVfi/Iw6D/7dBZ17BhnyGYys3hdkpQw6tdPuuKVXFprxCXPMq1RIuepEECHDQ2N06pDhG3SuJGK03pgGSkNmTAdKWzCOhFnSc6ZO18s02J5b/nyn2nFrytp3tzZ1LZtGyW9QPKYOfMj2rJtm7AZ6NipCw0dOljsUP3rr50ECWTH9m00fsIkUT8oL/EPjYE+7YP3qWXLFqJ6GHSdO3cinAj19qQJhGPlxo6bQJu3bBVXAcBmYPpHM2nWJ5+KHbCvv/66MATbuGmzMDn38DC9SoF/c2zUAmOBshLM1ZTojvBuXbsqKWL8xElil22FChUIZtpBR47QJ7M+Fub30PGA+YKJQN8AZSlO6dZKdlitdnYqT+c1VrVocHJ1kN+rePFiNHjQQGEqj412165dE8Znf27fITadQdk69b0pYkeyqfbIvPDcvuMv9eoWv01eBaSDJ8eMGTNmpEO5XGQ8Art376F7d+/RqJHDxcU8ep0IO1n9DhwUpuTo3Lj8Sc+h05ct5yg6Kf7dYOQF0/Mhgwcl+gdFWjANnMANReD1Gzfi95NkIze3mvTpJx+LQ19wFgOmDxcuXKBvFy8RS7QfTp+mpALUFdLHyy+/JJjG4iXfU/ny5WjcmNHUqVNHql69Gj1//owMz+MsW2HUhEE35Z3JNHTIYF07Etku5A2xXv5LD+jfL9GglnTyiba3fqMVFSiQn65HXRfMAybyTk7lacZHHwpmjGP+Y2NjKEf2HPTLil/pbnS0uBAJEoHWARswxUKFColvgis9oRuB4lfv+2jTIh6MDlvog4NDxIAPDAoS5u/9+vSmqe9PEQrrlPKReS5a9C2FhoUJK1kw5fSekvBFRvLLpNMTUgYGETq89l/OuDqgg0NnTsmBFiIwOiXok+ucoMM0ROaLOmjpES/zQrheHVF/TDWw2gJabadGnCyDKJvQA+jlodemy5cvU6fO3Sj67l3atmWTUuLq0WrDUB7KhTOus6wPwuGX7daml37kg5+HZ32aPHkS9fLuqdt+SW/8FN8BklF8PZIryzgt3lE/z7pehGXjiRPGJZGE9NJYO4x1GNZGOIX80YkwwFIaRKAzt8OBDubXyFc7+PWqgnJhjyDzN6ZHPPJBvKk6Ig3uyJW02nIQh7SQknLndjCZhzaN9GMK061b3A1nv2/YIINTfKIesj3GdZb1kTTJZQYaSCl3oqOFotQ4r+TSIg51gIGYxC8leuN4WMqCWeTJnVsomi0t3zi/tHhnhpEWKHIeVkEAA6R7t26UJ09u2rx5ayLrSKsUqJMprjmAotlSS1WdrCwO8vf3F2mgHMXO5IzgmGFkhK/AdTCJQIUKLkKhefvOHVq5cpVJOmtE3Lhxk378aTnJVSFrlGEqTywP74vfi9O3T28hpZiitWU4Mwxbos1lWYwApIw+fXqJYwCwIxUb9WzhoH/w8VkgzuDs2qWLRVOpF60fdBez58yjkNAw8u7ZQyhRXzTPtErPDCOtkOR8rIYADhL+YOr7YtXkt99WJbIbsVahsMP4zXc1TZ82lbRGZNYqT5uv7+o1hJ9T+XL06aezbMqstPXQ8/MqiR4qHJbhEMBqBaYHH8/6lJYtXSyMraypBITx1tatWwm2LVBe2srBLqdv/wEUExNDQQH+SUzqbVUPU+UwwzCFDIdnOATANObN9xEGZCtX/CyODMxwlXyBCqF92I9z9949WvDVPHJ0TGxS/wJZp1lSZhhpBiVnZAsEML+H9WaRokVImr/bolxblIG2HThwUOgszL1Pxhb10pbBDEOLBvsZAUYgWQRY6ZksPBzJCDACWgSYYWjRYD8jwAgkiwAzjGTh4UhGgBHQIsAMQ4sG+xkBRiBZBJhhJAsPRzICjIAWAWYYWjTYzwgwAskiwAwjWXg4khFgBLQIMMPQosF+RoARSBaB/wPY49vByXrRaAAAAABJRU5ErkJggg==)

*  **Davies-Bouldin Index** ('dbi'): An alternate and slightly more efficient way of measuring the spearation between clusters.

## K-means
"""

from sklearn import metrics
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
from sklearn.cluster import KMeans
from scipy.stats import mode
from sklearn.metrics import f1_score

# Function to create a Kmeans classfier
# which has been trained with appropriate data.

# Generic function to create a kmeans classfier. 
# This was kept separate as it initially took the input data and 
# trained it as well. Keeping this separate for now. 

def create_trained_classifier(k):
  clf = KMeans(n_clusters=k)
  print("Classifier built with {} Clusters".format(k))
  return clf

# Function to map predicted labels to actual labels 
# Labels predicted as 1, 2 may not correspond to actual '1', '2' in the output class
# Creating a function to map predicted labels to actual labels such that 
# if kmeans is specifying an instance belonging to label 1 (cluster 1) it is implying its 
# predicted actual output as 0 (for instance) and so on. 
# this would enable us to obtain accuracy score and confusion matrix related metrics as if k-means 
# (or any other clustering algorithm) was a supervised learning algorithm. 
def map_predicted_to_actual(labels_true, labels_predicted):
  mapped_labels = np.zeros_like(labels_predicted)
  unique_predicted_labels = list(set(list(labels_predicted)))
  for label in unique_predicted_labels:
    mapped_labels[labels_predicted == label] = mode(labels_true[labels_predicted == label])[0]
  return mapped_labels

# Function to generate the evaluation metrics for clustering algorithms.
def generate_metrics(labels_true, labels_predicted, data, metric='euclidean' ):
  metrics_dict = {}
  metrics_dict['ari'] = metrics.adjusted_rand_score(labels_true, labels_predicted)
  metrics_dict['ami'] = metrics.adjusted_mutual_info_score(labels_true, labels_predicted)
  metrics_dict['nmi'] = metrics.normalized_mutual_info_score(labels_true, labels_predicted)
  metrics_dict['homogenity'] = metrics.homogeneity_score(labels_true, labels_predicted)
  metrics_dict['completeness'] = metrics.completeness_score(labels_true, labels_predicted)
  # metrics independent of ground truth (actual labels assigned to data instances before hand)
  metrics_dict['sc'] = metrics.silhouette_score(data, labels_predicted, metric=metric)
  metrics_dict['dbi'] = metrics.davies_bouldin_score(data, labels_predicted)
  # supervised learning metrics generated as if this was a supervised learning algorithm
  mapped_labels = map_predicted_to_actual(labels_true, labels_predicted)
  metrics_dict['f1-score'] = f1_score(labels_true, mapped_labels, average='weighted')
  return metrics_dict

# Creating a dictionary to store the metrics obtained in different setups. 
# This will be used in the end to create a dataFrame of evaluation metrics 
# to plot and compare the perfromance of kmeans in different setups.

metrics_dict = {}

"""### K-Means without class field"""

# Setting number of clusters K as 10 for now as there are 10 classes.
k = 10 
x = input_data
y= output_data.astype(np.int32)
x = x.astype('float')/255

raw_k_means = create_trained_classifier(k)

labels = raw_k_means.fit_predict(x)
metrics_raw = generate_metrics(y, labels, x)
print(metrics_raw)
metrics_dict["Kmeans on Raw data"] = metrics_raw

"""### K-means including the class field"""

input_df = pd.DataFrame(x)
input_df["output"] = y
x_new = input_df.to_numpy()

k_means_clf = create_trained_classifier(k)
labels_new = k_means_clf.fit_predict(x_new)
dependent_metrics_raw_with_output = generate_metrics(y, labels_new, x)

print(dependent_metrics_raw_with_output)

metrics_dict["Kmeans on Raw data with output"] = dependent_metrics_raw_with_output

"""### K-means with PCA (dimensionality reduction)"""

x = input_data
y = output_data.astype(np.int32)

# reducing to 1000 dimensions (1000 being arbitarary here)
pca = PCA(n_components=1000)

x_reduced = pca.fit_transform(x)

k_means_clf = create_trained_classifier(k)

labels_with_pca = k_means_clf.fit_predict(x_reduced)

dependent_metrics_pca = generate_metrics(y, labels_with_pca, x_reduced)

print(dependent_metrics_pca)

metrics_dict["Kmeans with PCA on raw data"] = dependent_metrics_pca

"""### K-means with PCA (dimensionality reduction) (with labels)"""

x = input_data
y = output_data.astype(np.int32)
pca = PCA(n_components=1000)
x_reduced = pd.DataFrame(pca.fit_transform(x))

x_reduced['output'] = y
x_reduced = x_reduced.to_numpy()

k_means_clf = create_trained_classifier(k)

labels_with_pca_including_output = k_means_clf.fit_predict(x_reduced)

dependent_metrics_pca_including_output = generate_metrics(y, labels_with_pca_including_output, x_reduced)

print(dependent_metrics_pca_including_output)

metrics_dict["Kmeans with PCA on raw data with output"] = dependent_metrics_pca_including_output

"""## Gaussian Mixture Models (GMM)

### GMM without class field
"""

from sklearn import mixture

x = input_data
y = output_data.astype(np.int32)

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x)
gmm_metrics = generate_metrics(y, labels, x)
print(gmm_metrics)
metrics_dict["GMM on Raw Data"] = gmm_metrics

"""### GMM with class field"""

# appending output to input class

x_df = pd.DataFrame(input_data)
x_df['output'] = y
x = x_df.to_numpy()

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x)

gmm_with_output_metrics = generate_metrics(y, labels, x)

metrics_dict["GMM on Raw data with Output"] = gmm_with_output_metrics

print(gmm_with_output_metrics)

"""### GMM with PCA (dimensionality reduction)"""

# running PCA before GMM

x = input_data
y = output_data.astype(np.int32)

# reducing to 1000 dimensions (1000 being arbitarary here)
pca = PCA(n_components=1000)

x_reduced = pca.fit_transform(x)

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x)

gmm_pca_metrics = generate_metrics(y, labels, x)
metrics_dict["GMM with PCA on raw data"] = gmm_pca_metrics
print(gmm_pca_metrics)

"""## Visulising evaluation metrics"""

# Converting the nested dict of metrics to a dataframe for convenient plotting

index_of_experiments = list(metrics_dict.keys())

index_of_metrics = list(metrics_dict[index_of_experiments[0]].keys())

metrics_df = pd.DataFrame({'experiment': index_of_experiments})
# need one name and list of corresponding values to initialize the df

for name in index_of_metrics:
  list_of_values = [metrics_dict[x][name] for x in index_of_experiments]
  metrics_df[name] = list_of_values

metrics_df.head(10)

sns.barplot('f1-score', 'experiment', data=metrics_df)

sns.barplot('sc', 'experiment', data=metrics_df)

sns.barplot('dbi', 'experiment', data=metrics_df)

"""# Research Question"""

def image_equalise_hist(image: np.ndarray): 
    """ Equalise the histogram of the image globally.

        Reference:
        - https://en.wikipedia.org/wiki/Histogram_equalization
        - 
    """
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    #  Equalize the histogram of the image
    image = equalizeHist(image)

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)


def image_thresholding(image: np.ndarray): 
    """ Threshold the image using the Otsu thresholding
        algorithm. 

        Reference: 
        - https://en.wikipedia.org/wiki/Thresholding_(image_processing)
        - https://en.wikipedia.org/wiki/Otsu%27s_method
        - https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_thresholding.html
    """
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    # Find threshold using Otsu filter
    threshold: float = filters.threshold_otsu(image)
    binary = image > threshold

    binary_image = np.where(image, binary, 0) * 255

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def image_local_enhance_contrast(image: np.ndarray):
    """ Locally enhance the contrast of the image

        Reference: 
        - https://en.wikipedia.org/wiki/Histogram_equalization
        - https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
        - https://scikit-image.org/docs/dev/api/skimage.filters.rank.html#enhance-contrast
    """
    
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    image = filters.rank.enhance_contrast(image, morphology.disk(2))

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def image_local_autolevel(image: np.ndarray):
    """ Stretch histogram to normalise levels

        Reference:
        - https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
        - https://scikit-image.org/docs/dev/api/skimage.filters.rank.html#autolevel
    """
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    image = filters.rank.autolevel(image, morphology.disk(4))

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def image_enhance(image: np.ndarray):
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    image = filters.rank.enhance_contrast(image, morphology.disk(2))
    image = filters.rank.autolevel(image, morphology.disk(4))

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def enhance_all_images(data): 
    return np.apply_along_axis(image_enhance, 1, input_data)

def display_image_enhancement_comparison(image: np.array):
    def add_plot_line(im, axis):
        ax2 = axis.twinx()
        sns.distplot(im.flatten(), hist=True, bins=30, rug=False, ax=axis, kde=False)
        sns.distplot(im.flatten(), hist=False, rug=False, ax=ax2, kde=True, kde_kws={
            "alpha": 0.8, "cumulative":
                True, 'color': 'r'})
        sns.despine(ax=axis, left=True, right=True)
        sns.despine(ax=ax2, left=True, right=True)
        axis.set_yticks([])
        ax2.set_yticks([])
        axis.set_xticks([0, 50, 100, 150, 200, 250])
        ax2.set_xticks([0, 50, 100, 150, 200, 250])
        axis.set_xlim([0, 255])
        ax2.set_xlim([0, 255])

    hist_equalized_image = image_equalise_hist(image)
    autoleveled_image = image_enhance(image)
    thres_img = image_thresholding(image)
    loc_enhance = image_local_enhance_contrast(image)
    loc_autolev = image_local_autolevel(image)

    widths = [1, 1, 1, 1, 1, 1]
    heights = [4, 1]
    spec = dict(width_ratios=widths, height_ratios=heights)
    fig, axes = plt.subplots(nrows=2, ncols=6, figsize=(16, 5), gridspec_kw=spec)

    axes[0, 0].imshow(image_as_square(image), cmap='gray')
    axes[0, 0].set_axis_off()
    axes[0, 0].set_title("Original")
    add_plot_line(image, axes[1, 0])

    axes[0, 1].imshow(image_as_square(hist_equalized_image), cmap='gray')
    axes[0, 1].set_axis_off()
    axes[0, 1].set_title("Histogram equalization")
    add_plot_line(hist_equalized_image, axes[1, 1])

    axes[0, 2].imshow(image_as_square(thres_img), cmap='gray')
    axes[0, 2].set_axis_off()
    axes[0, 2].set_title("Thresholding")
    add_plot_line(thres_img, axes[1, 2])

    axes[0, 3].imshow(image_as_square(loc_enhance), cmap='gray')
    axes[0, 3].set_axis_off()
    axes[0, 3].set_title("Local contrasting")
    add_plot_line(loc_enhance, axes[1, 3])

    axes[0, 4].imshow(image_as_square(loc_autolev), cmap='gray')
    axes[0, 4].set_axis_off()
    axes[0, 4].set_title("Local equalisation")
    add_plot_line(loc_autolev, axes[1, 4])

    axes[0, 5].imshow(image_as_square(autoleveled_image), cmap='gray')
    axes[0, 5].set_axis_off()
    add_plot_line(autoleveled_image, axes[1, 5])

    fig.tight_layout()
    plt.show()

"""## Comparing different image enhancements"""

display_image_enhancement_comparison(input_data[0])

"""## Comparing models"""

enhanced_input_data = enhance_all_images(input_data)

"""### Naive Bayes with normalization and image mutations"""

# Normalize the data
enhanced_input_data_normalized = np.zeros(enhanced_input_data.shape)

for label in range(10):
    x = enhanced_input_data[output_data.astype(np.int)==label]
    normalized = normalize(x)
    enhanced_input_data_normalized[output_data.astype(np.int)==label] = normalized

compare_class_distribution_with_normalized(input_data, output_data)

# Run naive bayes on data after image preprocessing
x = enhanced_input_data_normalized
y = output_data

# split the data into training and test sets using stratified random sampling
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(x_train, y_train)
y_predictions = model.predict(x_test)

print_outcomes(x_test, y_test, y_predictions)
plot_confusion_matrix(y_test, y_predictions)

"""### Naive Bayes feature selection with image mutations"""

vizualize_important_features(enhanced_input_data_normalized, output_data)

"""### K-Means with image mutations

#### Without labels
"""

x = enhanced_input_data
y = output_data.astype(np.int32)

k_means_clf = create_trained_classifier(k)
labels_with_enhanced_input = k_means_clf.fit_predict(x)

dependent_metrics_enhanced = generate_metrics(y, labels_with_enhanced_input, x)

print(dependent_metrics_enhanced)

metrics_dict["Enhanced"] = dependent_metrics_enhanced

"""#### With labels"""

enhanced_df = pd.DataFrame(enhanced_input_data)
enhanced_df['output'] = y
x_enhanced_with_labels = enhanced_df.to_numpy()

k_means_clf = create_trained_classifier(k)
labels_enhanced_including_actual = k_means_clf.fit_predict(x_enhanced_with_labels)

dependent_metrics_enhanced_with_actual = generate_metrics(y, labels_enhanced_including_actual, x_enhanced_with_labels)


print(dependent_metrics_enhanced_with_actual)

metrics_dict["Enhanced with output"] = dependent_metrics_enhanced_with_actual

"""#### With PCA"""

# Running PCA on the enhanced data set
x = enhanced_input_data
y = output_data.astype(np.int32)

pca = PCA(n_components=1000)
x_reduced = pca.fit_transform(x)

k_means_clf = create_trained_classifier(k)

labels = k_means_clf.fit_predict(x_reduced)

enhanced_pca_metrics = generate_metrics(y, labels, x)

print(enhanced_pca_metrics)
metrics_dict["Kmeans with PCA on enhanced data"] = enhanced_pca_metrics

"""### Gaussian Mixture Models with Image Mutations"""

# enhanced input
x = enhanced_input_data

gmm = mixture.GaussianMixture(n_components=10)
labels = gmm.fit_predict(x)
print("fit complete")
gmm_enhanced_metrics = generate_metrics(y, labels, x)
metrics_dict["GMM on Enhanced data"] = gmm_enhanced_metrics

print(gmm_enhanced_metrics)

"""#### with PCA"""

# GMM with PCA on enhanced
x = enhanced_input_data
y = output_data.astype(np.int32)

pca = PCA(n_components=1000)
x_reduced = pca.fit_transform(x)

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x_reduced)

gmm_pca_enhanced = generate_metrics(y, labels, x_reduced)

print(gmm_pca_enhanced)

metrics_dict['GMM with PCA on enhanced data'] = gmm_pca_enhanced

"""### Comparing performance of image mutations on clustering algorithms"""

# Converting the nested dict of metrics to a dataframe for convenient plotting

index_of_experiments = list(metrics_dict.keys())

index_of_metrics = list(metrics_dict[index_of_experiments[0]].keys())

metrics_df = pd.DataFrame({'experiment': index_of_experiments})
# need one name and list of corresponding values to initialize the df

for name in index_of_metrics:
  list_of_values = [metrics_dict[x][name] for x in index_of_experiments]
  metrics_df[name] = list_of_values

metrics_df.head(10)

"""### Comparing performance of models before and after mutations"""

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

sns.set()
sns.set(style='whitegrid')
palette = dict(Yes="#69c597", No="#F08080")

# data sets: raw, enhanced
# models: naive full, naive class-best, bayes k2, bayes tan, kmeans, gmm

model_comparison = pd.DataFrame({
    'F-Score': pd.Series([
        0.562,
        0.86,
        0.618,
        0.880,
        0.249131,
        0.723776,
        0.247867,
        0.677946,
        0.219816,
        0.732641,
        0.847,
        0.945,
        0.842,
        0.935,
        0.200787,
        0.637211,
    ]),
    'Mutated': pd.Categorical([
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes'
    ]),
    'Model': pd.Categorical([
        'Naive Bayes',
        'Naive Bayes',
        'Naive Bayes\n(Class Best)',
        'Naive Bayes\n(Class Best)',
        'K-Means',
        'K-Means',
        'K-Means\n(with PCA)',
        'K-Means\n(with PCA)',
        'Gaussian Mixture',
        'Gaussian Mixture',
        'TAN',
        'TAN',
        'K2 (max.\n2 parents)',
        'K2 (max.\n2 parents)',
        'Gaussian Mixture\n(with PCA)',
        'Gaussian Mixture\n(with PCA)'
    ])
})

fig, ax = plt.subplots(figsize=(11, 4))

sns.barplot(x="Model", y="F-Score", data=model_comparison, hue="Mutated", palette=palette, ax=ax)
sns.despine(left=True)
plt.yticks(np.arange(0, 11) / 10)
plt.ylim([0, 1])
plt.xlabel('')

fig.tight_layout()
plt.show()
66ae80c43918b281a18f9098eba694b0bdcce092
# -*- coding: utf-8 -*-
"""DMML Final Project.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/15spI7M-iNDjDFKOLMx5P2roJ7fFvytPS
"""

# Commented out IPython magic to ensure Python compatibility.
# Packages
import pprint
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from cv2 import equalizeHist
from skimage import filters, morphology
from sklearn.metrics import confusion_matrix, classification_report, recall_score, f1_score
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.preprocessing import normalize
from sklearn.feature_selection import chi2

# %matplotlib inline

pp = pprint.PrettyPrinter(indent=2)

"""# Importing the data

## Importing Functions
"""

def import_data(seed: object = 42) -> object:
    """ Import data from files.

        Note: The data must be in the root of the current working directory
    """

    # Read input data
    df = pd.read_csv("x_train_gr_smpl.csv").astype(int)

    # label data-frame rows based on sample data
    for x in range(10):
        index = ~pd.read_csv("y_train_smpl_%s.csv" % x, squeeze=True).astype(bool)  # reversed flags (~)
        df.loc[index, 'label'] = str(x)

    input_data_ordered = df.iloc[:, 0:2304].to_numpy()
    output_data_ordered = df.iloc[:, 2304].to_numpy()

    # Randomise instance order (forcing the same result each time)
    np.random.seed(seed)
    permutation = np.random.permutation(df.shape[0])

    # Create base input and output arrays
    input_data = input_data_ordered[permutation]
    output_data = output_data_ordered[permutation]

    return input_data, output_data, df, input_data_ordered, output_data_ordered

"""## Implementation

### Import the data
"""

input_data, output_data, df, input_data_ordered, output_data_ordered = import_data()

# Image labels with classes as dictionary
labels = {
    '0': 'Speed limit 60',
    '1': 'Speed limit 80',
    '2': 'Speed limit 80 lifter',
    '3': 'Right of way at crossing',
    '4': 'Right of way in general',
    '5': 'Give way',
    '6': 'Stop',
    '7': 'No speed limit general',
    '8': 'Turn right down',
    '9': 'Turn left down'
}

# Show first 5 rows to make sure data imported
df.head()

"""# Data analysis

## Display Functions
"""

import numpy as np

def image_as_square(image: np.ndarray):
    image_width = 48
    image_height = image_width

    return np.resize(image, (image_width, image_height)).astype(np.uint8)


def image_as_array(image: np.ndarray):
    image_width = 48
    image_height = image_width
    return np.resize(image, (image_width * image_height,)).astype(np.uint8)

def is_data_unique(data: np.array) -> bool:
    """
      Compare the shape of the data set to a version which only contains
      unique images.

      If the shapes are equal, then the data only consists of unique images.
    """

    data_shape = data.shape
    unique_data_shape = np.unique(data, axis=0).shape

    if data_shape == unique_data_shape:
        return True

    return False

def display_images_as_grid(images: np.ndarray):
    image_count = len(images)
    column_count = 8
    fig_width = 11

    # Calculate the number of rows needed
    row_count = image_count // column_count

    if image_count % column_count > 0:
        row_count += 1

    # Calulate the figure size so each image isn't tiny
    fig_height = fig_width / column_count * row_count

    # Setup the plot
    fig, axes = plt.subplots(
        nrows=row_count,
        ncols=column_count,
        figsize=(fig_width, fig_height)
    )

    ax = axes.ravel()

    #  Disable axis on all subplots, even those without images
    for i in range(0, len(ax)):
        ax[i].axis('off')

    # Resize images to shape of (48, 48)
    images = np.apply_along_axis(image_as_square, 1, images)

    # Plot the images
    for i in range(0, len(images)):
        image = images[i]
        ax[i].imshow(image, cmap='gray', vmin=0, vmax=255)
        ax[i].set_title(i, color='black')

    fig.tight_layout()

def plot_frequency_distribution(data):
    data_labels = [v for v in labels.values()]

    data = pd.DataFrame({
        'count': pd.Series([
            np.sum((data == '0')),
            np.sum((data == '1')),
            np.sum((data == '2')),
            np.sum((data == '3')),
            np.sum((data == '4')),
            np.sum((data == '5')),
            np.sum((data == '6')),
            np.sum((data == '7')),
            np.sum((data == '8')),
            np.sum((data == '9')),
        ]),
        'class': pd.Categorical(data_labels)
    })

    # Set style
    sns.set_style("white")
    sns.set_style("ticks")

    # Colour palette
    sns.color_palette("hls", 10)

    # Create new figure
    plt.subplots(figsize=(11, 3))
    # Create plot
    sns.barplot(y='class', x='count', data=data, palette='hls')
    plt.xticks(np.linspace(0, 2250, 10))
    # Remove spine
    sns.despine(bottom=True)
    #  Add grid lines
    plt.grid(b=True, which='major', axis='x')
    # Remove axis labels
    plt.ylabel("")
    plt.xlabel("")
    plt.show()

def compare_class_distribution_with_normalized(input_data, output_data):
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(11, 2))

    sns.set_style("white")
    sns.set_style("ticks")
    palette = sns.color_palette('hls', 10)

    input_data_normalized = np.zeros(input_data.shape)
    
    for label in range(10):
        x = input_data[output_data.astype(np.int) == label]
        normalized = normalize(x)
        input_data_normalized[output_data.astype(np.int) == label] = normalized
        sns.distplot(x.flatten(), ax=ax[0], hist=False, kde=True, color=palette[int(label)])
        sns.distplot(normalized.flatten(), ax=ax[1], hist=False, kde=True, color=palette[int(label)])

    sns.despine(ax=ax[0], left=True, trim=True)
    sns.despine(ax=ax[1], left=True, trim=True)

    ax[0].get_yaxis().set_ticks([])
    ax[1].get_yaxis().set_ticks([])

    ax[0].grid(b=True, which='major', axis='x')
    ax[1].grid(b=True, which='major', axis='x')

    fig.tight_layout()
    plt.show()

"""## Implementation

### Check instances within data
"""

def instance_count_per_class(data):
    unique, counts = np.unique(data, return_counts=True)
    count_dict = dict(zip(unique, counts))

    return dict((labels[key], value) for (key, value) in count_dict.items())


# Number of instances in the data set
instance_count = np.shape(input_data)[0]
print(f"Number of images: {instance_count}")

# Check if each instance (row/image) is unique
print(f"Is data unique? {is_data_unique(input_data)}")

# Print number of instances per class
pp.pprint(instance_count_per_class(output_data))

"""### Plot class frequency distribution"""

plot_frequency_distribution(output_data)

"""### Display the images as a grid"""

display_images_as_grid(input_data[0:200])

"""# Naive Bayes

Applying Naive Bayes
1. Split the data in 70/30 for train and test sets respectively
2. Apply Naive Bayes Gaussian variant
3. Predict classification using trained model
4. Display classification report and performance statistics

Notes:
- Using non-random data model predicts 955 out of 3798 signs
- With randomized data model predicts 1125 out of 3798 signs 
- With normalization model predicts 2067 out of 3798 signs
- When overtly enhancing the contrast, model predicts 3257 out of 3798 signs
- Selecting the top featurns improves the model to 3383 of 3798 signs

## Functions
"""

from sklearn.metrics import confusion_matrix

def plot_confusion_matrix(y_true, y_pred):

  title = 'Normalized confusion matrix'
  classes = [v for v in labels.values()]

  cm = confusion_matrix(y_true, y_pred)
  cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
  cm_df = pd.DataFrame(cm, index=classes, columns=classes)

  flg, ax = plt.subplots(figsize=(12, 8))
  ax = sns.heatmap(cm_df, annot=True, cmap="Blues")
  bottom, top = ax.get_ylim()
  ax.set_ylim(bottom + 0.5, top - 0.5)
  ax.set_title(f'{title}')
  ax.set_ylabel('True label')
  ax.set_xlabel('Predicted label')
  plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")

  plt.show()

# Commented out IPython magic to ensure Python compatibility.
from sklearn.metrics import confusion_matrix, classification_report, recall_score, f1_score

def print_outcomes(x_test, y_test, y_predictions):
  cm = pd.DataFrame(data=confusion_matrix(y_test, y_predictions), index=labels, columns=labels)
  cr = classification_report(y_test, y_predictions, target_names=list(labels.values()))
  rc = recall_score(y_test, y_predictions, average='weighted')
  f1 = f1_score(y_test, y_predictions, average='weighted')

  print("Outcome:\nNaive Bayes correctly predicted the classification for %s of %s signs using %s attributes\n" 
#         % (sum(y_predictions == y_test), len(y_test), x_test.shape[1]))
  print("Classification report:\n%s\n" % cr)
  print("Confusion matrix:\n%s\n" % cm)
  print("Weighted Recall: %.3f" % rc)
  print("Weighted F1: %.3f\n" % f1)

"""## Running Naive Bayes with different filters

#### Without pre-processing
"""

from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# split the data into training and test sets for each class
x_train = []
x_test = []
y_train = []
y_test = []

for label in range(10):
    x = input_data_ordered[output_data_ordered.astype(np.int)==label]
    y = output_data_ordered[output_data_ordered.astype(np.int)==label]
    x_train_tmp, x_test_tmp, y_train_tmp, y_test_tmp = train_test_split(x, y, test_size=0.3, shuffle=False)
    x_train.extend(x_train_tmp)
    x_test.extend(x_test_tmp)
    y_train.extend(y_train_tmp)
    y_test.extend(y_test_tmp)
    

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(np.array(x_train), np.array(y_train))
y_predictions = model.predict(np.array(x_test))

print_outcomes(np.array(x_test),np.array(y_test), y_predictions)
plot_confusion_matrix(np.array(y_test), y_predictions)

"""#### With Randomizing"""

from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# Run naive bayes on randomized data set
x = input_data
y = output_data

# split the data into training and test sets using stratified random sampling
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(x_train, y_train)
y_predictions = model.predict(x_test)

print_outcomes(x_test, y_test, y_predictions)
plot_confusion_matrix(y_test, y_predictions)

"""#### With Normalization"""

from sklearn.preprocessing import normalize
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# Normalize the input data
input_data_normalized = np.zeros(input_data.shape)

for label in range(10):
    x = input_data[output_data.astype(np.int)==label]
    normalized = normalize(x)
    input_data_normalized[output_data.astype(np.int)==label] = normalized

compare_class_distribution_with_normalized(input_data, output_data)

# Run naive bayes on data without image preprocessing
x = input_data_normalized
y = output_data

# split the data into training and test sets using stratified random sampling
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(x_train, y_train)
y_predictions = model.predict(x_test)

print_outcomes(x_test, y_test, y_predictions)
plot_confusion_matrix(y_test, y_predictions)

"""### Interesting find (found too late)"""

# *** late interesting find resulting in f1-score of 0.97

from sklearn.preprocessing import PowerTransformer

pt = PowerTransformer(method='box-cox', standardize=False)
input_data_transformed = np.zeros(input_data.shape)

for label in range(10):
    x = input_data[output_data.astype(np.int)==label]
    transformed = pt.fit_transform(x)
    input_data_transformed[output_data.astype(np.int)==label] = transformed
    sns.set_style('darkgrid')
    sns.distplot(transformed.flatten(), hist=False, norm_hist=True)

"""## Feature selection

Using Chi2 to determine feature importance for each of the labels of enhanced and normalized datasets
"""

from sklearn.feature_selection import chi2

def display_image_from_array(image: np.array, title):
    # resize the image to 48 x 48 shape
    image = image_as_square(image)

    # Display the image
    plt.axis('off')
    plt.tight_layout()
    plt.imshow(image, cmap='gray')
    plt.title(title)
    plt.show()

def vizualize_important_features(input_data, output_data):
    X = input_data
    Y = output_data
    feature_importances = chi2(X=X, y=Y)[0]

    N, M = X.shape

    out = {}
    for c in set(Y):
        out[c] = dict(zip(range(N), np.mean(X[Y==c, :], axis=0)*feature_importances))

    for item in out.items():
        # this transformation is only done for visualisation purposes
        scores = np.fromiter(item[1].values(), dtype=float)
        scores *= 255.0/max(scores)
        display_image_from_array(image=scores, title=f"Label: {item[0]}, Sign: {labels[item[0]]}")

"""### Visualise important features for normalized data"""

from sklearn.feature_selection import chi2

X = input_data_normalized
Y = output_data
feature_importances = chi2(X=X, y=Y)[0]

N, M = X.shape

out = {}
for c in set(Y):
    out[c] = dict(zip(range(N), np.mean(X[Y==c, :], axis=0)*feature_importances))

for item in out.items():
    # print("\nLabel: %s, Sign: %s" % (item[0], labels[item[0]]))
    # this transformation is only done for visualisation purposes
    scores = np.fromiter(item[1].values(), dtype=float)
    scores *= 255.0/max(scores)
    display_image_from_array(image=scores, title=f"Label: {item[0]}, Sign: {labels[item[0]]}")

"""### Top features for each class"""

from heapq import nlargest


def get_n_best_fetures(X, Y, N, show_: bool = False, print_: bool = False):

    idx = []
    for item in out.items():
        top = nlargest(N, item[1], key=item[1].get)
        idx.extend(top)

        if print_:
            print("\nLabel: %s, Sign: %s" % (item[0], labels[item[0]]))
            print("Selected %s features: %s" % (len(top), top))

    selected_features = list(set(idx))

    if print_:
        print("\n%s selected feature indexes:\n%s" % (len(selected_features), selected_features))

    if show_:
        all_pixels = np.repeat(0, X.shape[1])
        np.put(all_pixels, selected_features, 255)
        display_image_from_array(all_pixels, title="")

    return X[:, selected_features]



from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

def run_naive_bayes(x, y, print_: bool = False):

    # split the data into training and test sets using stratified random sampling
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

    # train data sets and use model for predictions 
    model = GaussianNB()
    model.fit(x_train, y_train)
    y_predictions = model.predict(x_test)

    f1 = f1_score(y_test, y_predictions, average='weighted')
    
    if print_:
        print("\nUsing %s best attributes, f1 score: %.3f\n" % (x.shape[1], f1))
        print_outcomes(x_test, y_test, y_predictions)
        plot_confusion_matrix(y_test, y_predictions)

    return f1

"""### 2 best attributes for each class

Selecting top 2 features for each class results in total of 9 pixels selected for the model due to overlap of best features between classes. Seemingly random features perform suprisingly good predicting 2481 of 3798 signs (65.3%) in the test data.
"""

input_data_2_best  = get_n_best_fetures(input_data_normalized, output_data, 2, True, True)
run_naive_bayes(input_data_2_best, output_data, True)

"""### 5 best attributes for each class

Increasing the number of best attributes to 5 per each class results in a total of 26 attributes selected across all classes. The performance improves to correctly predicting 2868 of 3798 signs (75.5%).
"""

input_data_5_best  = get_n_best_fetures(input_data_normalized, output_data, 5, True, True)
run_naive_bayes(input_data_5_best, output_data, True)

"""### 10 best attributes for each class

Selecting 10 attributes for each of the classes results in a total of 43 attributes. The performance is slightly improved and the model predicts correctly 3021 of 3798 signs (79.5%).
"""

input_data_10_best = get_n_best_fetures(input_data_normalized, output_data, 10, True, True)
run_naive_bayes(input_data_10_best, output_data, True)

"""### Best number of attributes for each class

Looking for the overall best number of attributes to be selected for each class the model was run on a range of 1 to 2000 best attributes. This is considered a slow method as model needs to be built and assesed for each of the 2000 iterations.

As a result, the best overall balanced accuracy score of 0.891 indicated 136 as the best number of attributes to be selected for each class resulting in a total of 376 attributes. The performance seems to be best around 100 best pixels and drops ever so slightly when increased up to the maxiumum available pixels. This shows that filtering for the most informative attributes, in this case pixels, can not only improve performance and reduce the size of dataset processed and stored, but also is expected to have a positive impact on overal predictive power of a naive bayes model.
"""

x_array = []
y_array = []
best = [0, 0]
print("Looking for best number of attributes for each class...")

for idx in range(1, 2000, 1):
    x = get_n_best_fetures(input_data_normalized, output_data, idx)
    f1 = run_naive_bayes(x, output_data)
    x_array.append(idx)
    y_array.append(f1)
    if f1 > best[1]:
        best = [idx, f1]

print("Best number of attributes: %s, F1 score: %s\n" % (best[0], best[1]))

plt.plot(x_array, y_array)
plt.show()

input_data_overall_best = get_n_best_fetures(input_data_normalized, output_data, 3, True)
run_naive_bayes(input_data_overall_best, output_data, True)

"""### KFold cross validation

Data was split into 5 equally sized bins to perform kfold cross validation on each of the combinations of 4 train and 1 test set - total of 5 runs to validate the model and best number of attributes selected in previous step. The resulted balanced accuracy scores are comparable which indicates a high validity of the selected model and assumptions made.
"""

# Commented out IPython magic to ensure Python compatibility.
from sklearn.model_selection import StratifiedKFold

x = input_data_overall_best
y = output_data

skf = StratifiedKFold(n_splits=5)

for train, test in skf.split(x, y):
    f1 = run_naive_bayes(x[train], y[train], False)
    print("Train: %s, Test: %s, F1: %.3f" 
#           % (len(x[train]), len(x[test]), f1))

"""# Complex Bayesian Networks

1. Data Transformations to csv
2. Running various classifiers in Weka
3. Comparing the results based on f1-score (visualise with seaborn)
"""

# Stratify function
from sklearn.model_selection import train_test_split

def stratify(x, y, filename, selected_features):
    # split the data into training and test sets using stratified random sampling
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)
    
    # make them dataframes
    df_x70 = pd.DataFrame(x_train)
    df_x30 = pd.DataFrame(x_test)
    df_y70 = pd.DataFrame(y_train)
    df_y30 = pd.DataFrame(y_test)

    # Join them all up
    df_x = df_x70.append(df_x30)
    df_y = df_y70.append(df_y30)
    df_y.columns = ['Class']
    df_x.columns = selected_features

    df_stratified = pd.concat([df_x, df_y], axis=1, ignore_index=True)

    # export to csv
    df_stratified.to_csv(f"Weka_{filename}.csv", index=False)

# build csv files for Weka

# FULL DATASET:
stratify(input_data, output_data, "Full_Original")
stratify(enhanced_input_data_hist, output_data, "Full_Enhanced")

# BEST_2:
original_input_data_2_best = get_n_best_fetures(input_data, output_data, 2, True)
enhanced_input_data_2_best = get_n_best_fetures(enhanced_input_data_hist, output_data, 2, True)
stratify(original_input_data_2_best, output_data, "2_Best_Original")
stratify(enhanced_input_data_2_best, output_data, "2_Best_Enhanced")

# BEST_5:
original_input_data_5_best = get_n_best_fetures(input_data, output_data, 5, True)
enhanced_input_data_5_best = get_n_best_fetures(enhanced_input_data_hist, output_data, 5, True)
stratify(original_input_data_5_best, output_data, "5_Best_Original")
stratify(enhanced_input_data_5_best, output_data, "5_Best_Enhanced")

# BEST_10:
original_input_data_10_best = get_n_best_fetures(input_data, output_data, 10, True)
enhanced_input_data_10_best = get_n_best_fetures(enhanced_input_data_hist, output_data, 10, True)
stratify(original_input_data_10_best, output_data, "10_Best_Original")
stratify(enhanced_input_data_10_best, output_data, "10_Best_Enhanced")

"""# Clustering

## Evaluation of clustering algorithms

**Evaluation of Clustering Algorithms**

The evaluation of clustering algorithms is quite different from normal cases (of supervised learning) where there is a known (confirmed) value associated with the data passed. This is not always the case with clustering. Two classes of evaluation metrics exist depending on weather the 'class' (or label) associated with data is known or not. 

(source: https://scikit-learn.org/stable/modules/clustering.html#clustering-performance-evaluation)

1.   Ground truth dependent metrics

*   **Adjusted Rand Index** ('ari'): This is a measure of the similarity between the actual labels assigned to the instances and to the predicted labels. Are all the instances in one present in the other. It is symmetric in the sense that order of passing variables does not matter to the function. 1.0 is considered the perfect score that all the instances belonging to one actual label are predicted to be in the same cluster. A score 0.0 implies there is no correlation between the actual and predicted labels, the clustering is more or less random. 

*   **Adjusted mutual Info Score** ('ami'), **normalized mutual info score** (nmi): Mutual Information is a function that measures the agreement of the two assignments, ignoring permutations. NMI is the more prominent one present in literature, AMI was more recently proposed. Mutual Information as a mathematical concept is a function of the entropy of the class. A higher value (bound by 1.0) means more similarity of distribution of instances between actual and predicted labels, a lower value indicates that the instances are spread differently in predicted labels from actual values.

*   **homogenity**: Each cluster contains only members of a single class. (1.0 indicates completely homogenus class with only members of the same class, 0.0 means predicted label contains members from different actual classes). 

*   **completeness**: All the members of the actual class belong to the same cluster (1.0 for proper completeness i.e the statement is true for all classes. 0.0 implies all classes are spread through out different clusters). 

2.   Ground truth independent metrics:

*   **Silhouette Coefficient** ('sc'): A measure of how well defined the clusters are. A high score (upper bound 1.0) implies that the clusters are well separated. A score of 0.0 implies possible overlap between clusters. A further lower score (bound -1.0) implies poor clustering. The Silhouette Coefficient is defined for each sample and is composed of two scores: **a**: The mean distance between one instance and all other instances of the same class. **b**: The mean distance between one instance and all other instances in the nearest neighbouring cluster. 

![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAQwAAABuCAYAAADSxy12AAAflklEQVR4Ae1dB3RURRe+tNAVlK4QIAkQOoEkQOhNeofQe1M6oiiCglhpBkQFREURCU26gkgNEJLQISQhEAgl1ICA4A8J7H++SWbysnmb7Ibspuydc/a8eTN32vdm7t65c2cmm8FgMBA7RoARYATMQCC7GTRMwggwAoyAQIAZBncERoARMBsBZhhmQ8WEjAAjwAyD+wAjwAiYjQAzDLOhYkJGgBFghsF9gBFgBMxGgBmG2VAxISPACDDD4D7ACDACZiPADMNsqJiQEWAEmGFwH2AEGAGzEWCGYTZUTMgIMALMMLgPMAKMgNkIMMMwGyomZAQYAWYY3AcYAUbAbASYYZgNFRMyAowAMwzuA4wAI2A2AjnNpmRCu0AA5ylpz1R6/vw55czJ3cQuPr4ZjeSeYAZI9kDy7Nkzmv+VDz18+JCuXYuiBw8eULZs2Sg6Opp2/rWdsmdnYdQe+kFKbWSGkRJCdhIfdi6cvl70bZLWNmncKEkYB9gvAvy3Yb/fPlHLK7tWovBzIbRr5w4a/dYoFefs7MzShUKDPcwwuA8oBHI7OJCTU3k6cuSoCvPyqqf87GEEmGFwH0iEAJScZ4LPijCHXLnIza12onh+sW8EmGHY9/dP0voLFyLo0aNHIrxuXU8qWCB/EhoOsF8EmGHY77fXbXnQkSMqHAyDl1QVHOwhImYY3A0UApiOHD4coN696tdXfvYwAkCAl1W5HygEwDBOnTot3l1dK1G1alUJYbDRwPV4ObJnF7YZWc0mA4ZqsbGxot2wPUH7slob1Ud+QQ8zjBcEMCslDwkNpUuRkaJJHu51KCAgkH748Sc6GxJCd6PvUunSr1PzFi1o8qQJ5ODgkCWa/vjxY9qwcROdPnWaDh7yp9y5HcjLy4vGjnmLChUuTGFhYVS4cGEqVbJklmjvizYiG9+t+qIQZp30P/y4nGZ+PEs0qFKlihQaGkbdu3Wlxo0a0uP//qO9e/fRn9t3ULu2beibRQsz9b8wpIrw8PM0/aMZ5O9/mPr26U3VqlWj8PBztGnzFnJ2ciIPD3da+PU3Ao8L4aGUK1eurPOxU9kSljBSCVxWS4aph7+/v2oWmMVCn/nUsWMHxRggroNhbPvjT2q6/nfq2aO7os9MHjCL/fv9aOq06XTlytUk7ezQvj31HziYDgcEimbV9fRQGGSmdlqjrqz0tAaqmTDP27fvCJFcVh3MonPnTokGStMmTWQ0rV69Rug3VEAm8kCiAEMAs4CuRssU0Qw3t1o0eNAA1aIyZcokwkFF2KGHGYYdfnS9JmM5VdpfdO3Smdq3b5eErEiRV1XY9es3MiXDgHJzydJlqh3jxo1NwgwggVy9ek3RODs7CWWvCrBjD09J7Pjja5seEhKiXiu7uuraX2C1RLqr165RbOwzXTpJo30+jYmh7NmyaYNe2J8jRw6LBjIYwS8rfqU9e/eKstu0foPatWmdpB6Ynp08eUqF16nD1q4SDGYYEgk7fuJfNzAwSCHg6emu/FoPVk2kK1vWkXLkME9AffDwIXl796bgswlMSebzIs85s78g7549zM7i0ePHtFQjXVSvVk037dOnT+lCRISIK1a0KNWqWVOXzh4DmWHY41c3avMzzf6R/Pnzk6urqxFF3Oup03E2Gngr6+ho9r97gfz5acq77xAkgrR09etbtjEuMjKSoq5fV1WoWbOG8ms9Z84Eq1d39zpmt1MlysIeZhhZ+OOa27ToO3eU/qJqlcq60wyI6Vor0FKlSpmbvdARNGnS2Gx6axFevBRnY4L8YYRWq5a+5HDs+HFVBUxH2IhLwcGm4QlQ2K9Pq+CrWrWK7gCJiIigffv9FEj4d85M+0yePzfQ9agoVf8aNapT3rx51bv0gDFqp2ewxYCD/kNag0pae3yaNwm1R2TsqM0YJNKVK1dOehM9d++JUxQisJd3T+rRvVui+Iz+An0rdBPSgeHBrsTY3X/wgAIC43Q1FSq4EBTAcOs3bKTyzhXtnmkwwzDuMXb4XqxYUdVqxzJllF96YmJiyc/vgHit4OJMkyaOT3N9hCzLWk8whyKvJiwLFyhQQLeoAwcO0oMHD0WcW61aStrasnkLde7UUTeNPQUyw7Cnr22irUWKFFExp0+fUX54IIr7LFggpiM5c+agCRPGUfHixRPRZJaXEiVKqKqiXcYOe2Z8fBaoYBhwgdFcvHiJ9uzdR61atshU0zDVkDT0MMNIQzAza1YFChSkwYMGiurv8/Mj2EzAYaoy/cOPxOHAMNr6fsliatumja4onxna3qCBl5hOoa7YaHb9xg3VTpi7jx07npo3b6aaki9fPoHBT8uXU/16dallyxYqzl49vPnMXr+8Ubuh0Bs9ZpzYK4K9E56eHnTuXLh4dypfjmZ/+QXVru2mRHSj5JnmNSYmhrx79xXnlhYsWJCGDB5IwcFn6e9du2nUyBE0aeIE+vyLL+in5b8QTkwvWrQorV23ntat9SX3OrzEygwj03R161cUlpzY6r3f7wA9uH+fsJJQo0YNcqtVkwoVKmT9CtioBDBHWHueOHGSQkJCBSOE7QkYBJZQYcG6Zu1aOnbsBF2+HEkdOrSnPr17ZTq9jTXgZIZhDVQzeZ6YimCOj/l7VrZBkKbuaKPeiolcPcrKGFjaVZlhWIoY0zMCdowAKz3t+ONz0xkBSxFghmEpYkzPCNgxAsww7Pjjc9MZAUsRYIZhKWJMzwjYMQLMMOz443PTGQFLEUjX7e1YD4cpcmhYGF29epVKlixJr7/2GjVs2IDXvC39kkzPCNgAgXRjGLC4m/bhDFq1ylfsTShX1lEc+x599y7hTMl5c2cz07BBB+AiGAFLEEgXhgFmMWbseGF2PH3aVBoyeJAwnImJfUZHgoKod9/+gmk0atTQkrYIWkgtekY4Fmekk8CUgY8OKQcxAlkSAZszDFjXSWbRsmVzGjZ0iBrgBw7sE5udgPQ3334npiaWDH5Y5n3ls4CCghIuFE7Lr4bLbjrxFue0hJTzymQI2Jxh3LlzR0gWwGlAv36KWeD92LFj6nj3/v37WQwlmEvlypUtSIcDVJJuczaVQW0+PdoUNBxuJwjYnGEcOpRwu9a69evJ2cWZSpYoIfYsTJwwnsLCzpGjYxlq4FU/ETMx53uAYbRt01r8zKEX+yWwj8AcYqZhBBgBsvlekr92/k3Dho9U0IM5NGrYkOrVq0te9esRthxbet+EyiyDe2Z+/InYVp3Bq8nVy2QIvFrkVVr+Y8LlTNasvs0ljKZNGhPutLgUf4JzZORlWhG5klb8upJwPNy7704WEkJaH0lvTRAtyfvkqYQLcixJx7SMgCkEmjVraioqzcNtLmFgGrB79x5xES6u29Nz8+fNoW5du1g8JdHLKyOF7fx7V0aqDtclCyHQskVzm7TG5gwDrQLTOH3mjDjO/cyZs3T02FGCpCFdq5Ytacnib1Jlh4FVGORvDYdlVT4bwRrIcp6ZBQGbMYx79+6R38FDdDkyUlw95+VVX2AkDynxXb2G3nv/AxFWqNDLdDQogHLlymURjshr8jtTaN363y1KZy7x5599QlhaZccI2CsCNmEY+Nfv0KkLySvoZn08kwb075toygGpwKtBY8Ilv+3btaUFPvMtZhj4iMFnz9Ld6LtW+Z5QzGamy3usAgJnatcI2ETpef7CBcUsgHZZxzKJmAXCnj17Tq+9VkowDBxCa6l0Ib9iFYvsMGQqfjICjIA5CNhkt+qrr7yi6oLbtvUu0fVdvZoCAoMIV/V5e/dU9OxhBBiBjIOATRjGK6+8QgMHxFluQs9w9Ogxun37Nv3zzz906tRpWrhwEU39YDrVdnOj9evWUO7cuTMOQlwTRoARUAjYRIeB0rApzGfBQlq9ei3dvHVLVQC3aTk6OopbpSa/PSnVUxGVIXuyNALYuPjVgoVUrGhRGjRwQJZqK/5MZ8+dRx3at6OMOrW2GcPAlwUg//33H926dUswjcKFClOZMqUFk+AlyyzV963SGCjPsZKG7QVfL/QhXGWYlRwU/x9+NIN+/uVXWrtmlbg4KaMt49uUYWSlj8ttsS0CYBYTJr5NmzZvIZ/5c6lr1y62rYAFpaGuhw4dphs3b9ClS5dUyrqenmIHtgrQ8cD8oFeffuKCpV1/7yAXZ2cdqvQLsskqSfo1j0vOCghAMl28ZKlgFv36ZvwjBuZ/5SPuozXGHgwjJVe4cGEaP24sjXpzNA0fPpLWr19L2kWDlNJbO94mSk9rN4Lzz9oIHDlylL6cPVfoLYYPG5oqC2BbIlSpYkUaN3Y0jRn9ZqJiXVxcEr2besGO60ED+1PExUs0Zcr7BIklozhmGBnlS3A9dBGAdPHzLytE3JAhg6hcuXK6dBkpEHexQoHfrl1bVa0ypUtT0aJF1HtKHhz14OpaibC7e+++fSmR2yyeGYbNoOaCUoMABsyWrdvEGSkjRwxPTRbplubokaOqbJzvYokCE1MT7OyGW79+Q4aRMphhqE/KnoyGAJbi//jjT1Gtxo0aZvipiBa/uA2WwSrIq4FXEutmFWnC07BBAxGzddsfVjt20kTRJoOZYZiEhiPSG4Hz5y/Q1m3bRDUaN26S3tWxqPznzw106vRplSY1S8CwiMYZMXBr1/0uzBJUhunkYYaRTsBzsSkjsGXrVoqNfUYFCxQQRzamnCLjUNy//49YGkWNYGRWonhxiyuHIycbNPAS6fb77aenT59anEdaJ2CGkdaIpiI/KPbkTy+5jMMzOY05RHhJm9ozQZBOlpPaPFAP+UuuPcnlj3aGh58XyWvWrEF58uTRy8qssBfFxKxCjIj8DweokFq1aorpiGwv6iP9isiEBzu34W7evJUhjndkOwwTH8oWwTjmBze/7du3T+yvcXZyojZtWlPt2m6ieAya7Tt20sWICAoJDRUWsqBxq+0m7m2RxxhiS/+JEycpMjJSHKKMvTulSpUkKAlfeukls5qCTowpAA42OnHiBN29e49w9Funjh2Esg7xcMlt70d99+zdR2FhYaI92H1crWpV6ty5Mzk45BJMZM3adXT8+AnKmzcPNWvaVPyD6uWJAXUmOE4H4OHhTtksPKkZ6XGiG9pz+vRpCggIjCuvoRdVrVJFMF7QpHZXdHKgIt+AgASGASkhMDCIIDFduXKV8ubNS6Vff53Gjh0tvk9yV2nUretJDg4OQrrA5kxMUyxRniZXz9TEsYSRGtTSKM3Zs2dp4sRJYus/NtwtXfYDdenWQ2zIu3v3Lr3z7nv05lujaemyZZQ3X1568vQprVzlS29Pfpd+/XWl+Jfy8ztAI0aMovenTqNdu3aL3b679+wVhkO9evelmzdvplhbSAOfff4ltWjVmpZ+v4wePXpM0dHRNHHSZJoxcxZt2LCRWrdpR9179hKDXi9D5IE7YYYMHS6YBZYE/f0P0+R336P33p9K//vf/2jKe1Pps8+/oGLFitLZsyE0ZNgIGj1mnK7UdOPmTXXlBK7QtMSBcf22ypfq1m9Aw0eMEmV5eHrQom+/o5Gj3qJtf/xJ3Xp4U+OmzWnjxk2WZG0WLZgrBrd0OMe1Z68+dPt2NNWoUV1cs4FvPXDQUAo7d06S6T7BHEqUiJvO7N+/X5fGpoEGdumCQExMjGHY8JGGfgMGGWJiYw1RUdcNpR3Li9/sOfMMHTt1MTRr3tKwdt16Q2zsM1HH2NhYQ/cevQRN8xatDBs3bjJUrFTFMHvOXENAQKDh+fPngu7UqdMqrw+mTTc8exaXXq+hT58+Nbzz7hRBP2jwUEN0dLTKJyAwUOUj6/a/J0+SZIN6/fDjT4L2t1W+qrzQ0DCDS0VXET5v/leGOu51DSdPnjSgzE8+/VzlfSY4OEmeO3f+reL9/A4kiTcVgLp8t3iJShsREaHqc/nKFUOr1m0NTs4VVbyv72pTWaU6HBhKvPB0q+1h8PVdY0Dd4PA9mjRtLmjGjB2nwk0V2L2Ht8rviQ7+ptJZI5wlDJuy54TCjp84QTv+2kktmzennDly0K1bCZLA14u+oX///ZcWLVpI3bt1pRw54j4T/m3c3WuLTM6Fn6dlP/wo9lS8PWkixYntcXJ75cquqqB9+/1MzpchFYwbP4F8V68V9MOHDyVMZ6SI7OHuTmPHvKXyqlmjuqirCog/n3XLlq1CEnlz1Ajq3ctbicwuLs5UwaWCIPdZ8DX17NmDqlevTpCAliz9XmWTL29e5ZcerWSE6ZU5Dv/sMMuGtAQHa0sYekkRXkwDxoympzExKjtrnLh9WKO/qODiTMt/+oG8vXuoZWHURxp1bdq8lfbu81P10fNICQNxDx8+1COxWRgzDJtBnbig4DPBQvvfqlVLEfHgQUJHcHF2EszCtVKlRIkwIC5ciFBhZcuWpU9mzVQDQkZcv35deunJkye6DAN5QWzf9sd2QTtwQH+qX6+eSic9VatWlV6xcUoOPhkI8f/LOXOpWtUqNGH8OBksnog7fz5OcYmATp06iHBMd6QDEypdurR8FU8wMu0RCOZOSY4dP05Ll8bdz4Ep0aSJExLlixd3ze11ONkNDDItHfQX0ClJ17r1G1StWgKGMlwbFhoaovuNJK22/Q8f/iuD0+XJSs90gB2dCvNoHISM+TzeD2uUZJMnv03GzEJWU56Lmj9/fpoy5Z0kzAJ0UVEJDKNWzTgNvUwvn9euRanB9VLBgjRs2BAZpZ6oF/QQ0kEZK6UPGQZlHvKaNvV9ocyT4fL5+L//hBeHI1WI30vRo3s3IalEXr5M/Qf0T6JIBVPCqgAcNl5B6ZeSA5NZsuR7oecBLZiXMXND+IGDh1RWTk5OujSKIBUeMGLolaSrp8OEESdXgOA/duy46APG2Mo8SpQoIb307yNmGAoMe/GgY/iuWimai06Nf+IDBxI6cp34VRJjPC5fvkJXrl4VwfhH1/7zaGkDAgPVKwa58cBBp8buz8tXrgi6OnVqE/Y6GDvQHTkaZ96MVQ33OnWMScSNdRHnw5IMehBKK034PT3dVVqsTGB6kpyTEoZ2sCRHv2bNWjHFA81LLxWkhvH2C8Zp/DVXdWLKZGqQGqcz9x2nyZ04GXdZVb58+QgSlLEDIz4ZT4M47DFJrh64SlS69J6SsIQhv4SNn3JJFMVei4oi6DTgcKapKTH5yJGEW+mrVKlM2XXWGtEZ9+9PmBNDijHujGBQu/fsUS12d6+ThAaR+BeUEk2dOnUoj46uAXmbWhbduy9Bq6833VEVMPIgz4RV1JTvmAFj00oOnh4eBAnM2KHd/ocTJCYsUae1C9dMwcCssYRq7FAPLKdLV7VK1WQlDC0C2bOlrxYhfUuXiNn5MygoYQkOR7MZSwSAB4zgyNFjCikwFmNGgMjw8HC1pAdJpWKFOKWjSkgk9CCYRkgHhameC9RIKu61aydReOqlkWGYIhw9FlffkiVLUKNGDWWUWU8pWUSZuB1PmwkGoNYMG8paPWxu3b5NV65eE0lz5cxJnp4e2mxe2I9vBBsU6YCZ3re8ePEiRWn0TNC36NHJfKKiEr5VgYIFZHC6PJlhpAvsCYXi3/Ho0eMqAFMNvc6OQXEsfgCCuFrVaiqN1oMVCOm00gXKQYeGOxOc8O8G8d2tVpyhmEyHJ2gDgxIkGq0UIvPR0hv7MYDl/bnQX+g55GMqL7kygBOocKxjcg42G5iuSWeKAYaEhEoSwjRM6kaSq4dKYIYH+WgZcas34hTaxknllAXhYKYpnZOhvVI0f758xtnZ9J0Zhk3hTloYOtnx4wkMo0aNGkmJiIQxV9i5cBGHm+HKlSubhA5MYdfuhKlG82bNxFIeyvDu1YfmzJ0nBiisDaWD+J49e8IEQIYjL6m/QJj8N4bFZE/vPsrYCnk/fvxYKPru378vkyfaXYmNV1it0TowwI9mzBSnxWvD4QfDlBIG3q9plLjGtHi/Fa8ghR8MULsCoaUPCQlRr5hioRy0E9jMnTdfxaXWg7wePHggkuMb6U15QANFsXSDBw6kgilIDVoJ4+VChWTSdHkyw0gX2BMKvXMnmoLPxnVkR8cyQoeREJvg2++XoJdo1LChrt4Ay4oY0HBVKruqgQOTcVgeuji7iEGCzVDSQRzWk2iCgo6o1RbUS4rMuIYSSlVcPAW3ync1tW3Xkfr2H0jDho8UAxBM5KBmNaJG9epJro7YufNvWv7zCqpfr66sSqJnKY11Z1RU3DQiEYHmxbWyKznEX6tpakqHKRKsS6WrVzfuuLzz5yPocEAgaZePQYOzOLHsvPWPP+nJE/M2fQEj7BuBK+voqOwuZJl4wix+9Zo4uxfQYnVKD39tmowkYbDSU/tl0sGvXdHAmQ9yYGqrYjwAGzdupNvJtJp3GCRBsYq0Py7/mTAA27VrI7ItVqyYyj5a51pJ/Pt/v+wHRYPlUNQL/46wMcAx+DAmwzsGQMTFi4IWS6go7+rVqwSDMenKlnWUXvHErsvvFi+h6tWqir0ziSLjX7DhTDrtMrEM0z5zOzgIRnvs+Am6c+eONkr4UScognEQD1zxYsWUxLR9x3ZyrVSRmjVN2D6Pdn3z7WI1sKEv2rThd8qVK/nhAozk9EKu8mgrA1xhbCfdkEEDdRm/jMcTdbl+I26ZHP1DT8Gspbe2nyUMayOcTP7oyFqFZ6NG+owAneZSZKTICaKuPInJOOvCGnFVnjZ98uRJWr78F3pvyjuqs4GZ9O3dSyQPDg5W0wsE4J+4S9futPPvXSp7OeAPHvIXqyYdO7RXzOjRo0eKbsSwoYKxyH9QGaGds0MygSXr8RMnadq0qapOklY+HRxyU614pqEVyWW89omB2ju+PeHnLySSboDdT8t/pkFDhlHzZk1FMugN8K+OAQx7mBYtmifahIYt9do2YKXIZ8ECbZEm/Tj0BhIZpAIthqjHgoVfi/Iw6D/7dBZ17BhnyGYys3hdkpQw6tdPuuKVXFprxCXPMq1RIuepEECHDQ2N06pDhG3SuJGK03pgGSkNmTAdKWzCOhFnSc6ZO18s02J5b/nyn2nFrytp3tzZ1LZtGyW9QPKYOfMj2rJtm7AZ6NipCw0dOljsUP3rr50ECWTH9m00fsIkUT8oL/EPjYE+7YP3qWXLFqJ6GHSdO3cinAj19qQJhGPlxo6bQJu3bBVXAcBmYPpHM2nWJ5+KHbCvv/66MATbuGmzMDn38DC9SoF/c2zUAmOBshLM1ZTojvBuXbsqKWL8xElil22FChUIZtpBR47QJ7M+Fub30PGA+YKJQN8AZSlO6dZKdlitdnYqT+c1VrVocHJ1kN+rePFiNHjQQGEqj412165dE8Znf27fITadQdk69b0pYkeyqfbIvPDcvuMv9eoWv01eBaSDJ8eMGTNmpEO5XGQ8Art376F7d+/RqJHDxcU8ep0IO1n9DhwUpuTo3Lj8Sc+h05ct5yg6Kf7dYOQF0/Mhgwcl+gdFWjANnMANReD1Gzfi95NkIze3mvTpJx+LQ19wFgOmDxcuXKBvFy8RS7QfTp+mpALUFdLHyy+/JJjG4iXfU/ny5WjcmNHUqVNHql69Gj1//owMz+MsW2HUhEE35Z3JNHTIYF07Etku5A2xXv5LD+jfL9GglnTyiba3fqMVFSiQn65HXRfMAybyTk7lacZHHwpmjGP+Y2NjKEf2HPTLil/pbnS0uBAJEoHWARswxUKFColvgis9oRuB4lfv+2jTIh6MDlvog4NDxIAPDAoS5u/9+vSmqe9PEQrrlPKReS5a9C2FhoUJK1kw5fSekvBFRvLLpNMTUgYGETq89l/OuDqgg0NnTsmBFiIwOiXok+ucoMM0ROaLOmjpES/zQrheHVF/TDWw2gJabadGnCyDKJvQA+jlodemy5cvU6fO3Sj67l3atmWTUuLq0WrDUB7KhTOus6wPwuGX7daml37kg5+HZ32aPHkS9fLuqdt+SW/8FN8BklF8PZIryzgt3lE/z7pehGXjiRPGJZGE9NJYO4x1GNZGOIX80YkwwFIaRKAzt8OBDubXyFc7+PWqgnJhjyDzN6ZHPPJBvKk6Ig3uyJW02nIQh7SQknLndjCZhzaN9GMK061b3A1nv2/YIINTfKIesj3GdZb1kTTJZQYaSCl3oqOFotQ4r+TSIg51gIGYxC8leuN4WMqCWeTJnVsomi0t3zi/tHhnhpEWKHIeVkEAA6R7t26UJ09u2rx5ayLrSKsUqJMprjmAotlSS1WdrCwO8vf3F2mgHMXO5IzgmGFkhK/AdTCJQIUKLkKhefvOHVq5cpVJOmtE3Lhxk378aTnJVSFrlGEqTywP74vfi9O3T28hpZiitWU4Mwxbos1lWYwApIw+fXqJYwCwIxUb9WzhoH/w8VkgzuDs2qWLRVOpF60fdBez58yjkNAw8u7ZQyhRXzTPtErPDCOtkOR8rIYADhL+YOr7YtXkt99WJbIbsVahsMP4zXc1TZ82lbRGZNYqT5uv7+o1hJ9T+XL06aezbMqstPXQ8/MqiR4qHJbhEMBqBaYHH8/6lJYtXSyMraypBITx1tatWwm2LVBe2srBLqdv/wEUExNDQQH+SUzqbVUPU+UwwzCFDIdnOATANObN9xEGZCtX/CyODMxwlXyBCqF92I9z9949WvDVPHJ0TGxS/wJZp1lSZhhpBiVnZAsEML+H9WaRokVImr/bolxblIG2HThwUOgszL1Pxhb10pbBDEOLBvsZAUYgWQRY6ZksPBzJCDACWgSYYWjRYD8jwAgkiwAzjGTh4UhGgBHQIsAMQ4sG+xkBRiBZBJhhJAsPRzICjIAWAWYYWjTYzwgwAskiwAwjWXg4khFgBLQIMMPQosF+RoARSBaB/wPY49vByXrRaAAAAABJRU5ErkJggg==)

*  **Davies-Bouldin Index** ('dbi'): An alternate and slightly more efficient way of measuring the spearation between clusters.

## K-means
"""

from sklearn import metrics
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
from sklearn.cluster import KMeans
from scipy.stats import mode
from sklearn.metrics import f1_score

# Function to create a Kmeans classfier
# which has been trained with appropriate data.

# Generic function to create a kmeans classfier. 
# This was kept separate as it initially took the input data and 
# trained it as well. Keeping this separate for now. 

def create_trained_classifier(k):
  clf = KMeans(n_clusters=k)
  print("Classifier built with {} Clusters".format(k))
  return clf

# Function to map predicted labels to actual labels 
# Labels predicted as 1, 2 may not correspond to actual '1', '2' in the output class
# Creating a function to map predicted labels to actual labels such that 
# if kmeans is specifying an instance belonging to label 1 (cluster 1) it is implying its 
# predicted actual output as 0 (for instance) and so on. 
# this would enable us to obtain accuracy score and confusion matrix related metrics as if k-means 
# (or any other clustering algorithm) was a supervised learning algorithm. 
def map_predicted_to_actual(labels_true, labels_predicted):
  mapped_labels = np.zeros_like(labels_predicted)
  unique_predicted_labels = list(set(list(labels_predicted)))
  for label in unique_predicted_labels:
    mapped_labels[labels_predicted == label] = mode(labels_true[labels_predicted == label])[0]
  return mapped_labels

# Function to generate the evaluation metrics for clustering algorithms.
def generate_metrics(labels_true, labels_predicted, data, metric='euclidean' ):
  metrics_dict = {}
  metrics_dict['ari'] = metrics.adjusted_rand_score(labels_true, labels_predicted)
  metrics_dict['ami'] = metrics.adjusted_mutual_info_score(labels_true, labels_predicted)
  metrics_dict['nmi'] = metrics.normalized_mutual_info_score(labels_true, labels_predicted)
  metrics_dict['homogenity'] = metrics.homogeneity_score(labels_true, labels_predicted)
  metrics_dict['completeness'] = metrics.completeness_score(labels_true, labels_predicted)
  # metrics independent of ground truth (actual labels assigned to data instances before hand)
  metrics_dict['sc'] = metrics.silhouette_score(data, labels_predicted, metric=metric)
  metrics_dict['dbi'] = metrics.davies_bouldin_score(data, labels_predicted)
  # supervised learning metrics generated as if this was a supervised learning algorithm
  mapped_labels = map_predicted_to_actual(labels_true, labels_predicted)
  metrics_dict['f1-score'] = f1_score(labels_true, mapped_labels, average='weighted')
  return metrics_dict

# Creating a dictionary to store the metrics obtained in different setups. 
# This will be used in the end to create a dataFrame of evaluation metrics 
# to plot and compare the perfromance of kmeans in different setups.

metrics_dict = {}

"""### K-Means without class field"""

# Setting number of clusters K as 10 for now as there are 10 classes.
k = 10 
x = input_data
y= output_data.astype(np.int32)
x = x.astype('float')/255

raw_k_means = create_trained_classifier(k)

labels = raw_k_means.fit_predict(x)
metrics_raw = generate_metrics(y, labels, x)
print(metrics_raw)
metrics_dict["Kmeans on Raw data"] = metrics_raw

"""### K-means including the class field"""

input_df = pd.DataFrame(x)
input_df["output"] = y
x_new = input_df.to_numpy()

k_means_clf = create_trained_classifier(k)
labels_new = k_means_clf.fit_predict(x_new)
dependent_metrics_raw_with_output = generate_metrics(y, labels_new, x)

print(dependent_metrics_raw_with_output)

metrics_dict["Kmeans on Raw data with output"] = dependent_metrics_raw_with_output

"""### K-means with PCA (dimensionality reduction)"""

x = input_data
y = output_data.astype(np.int32)

# reducing to 1000 dimensions (1000 being arbitarary here)
pca = PCA(n_components=1000)

x_reduced = pca.fit_transform(x)

k_means_clf = create_trained_classifier(k)

labels_with_pca = k_means_clf.fit_predict(x_reduced)

dependent_metrics_pca = generate_metrics(y, labels_with_pca, x_reduced)

print(dependent_metrics_pca)

metrics_dict["Kmeans with PCA on raw data"] = dependent_metrics_pca

"""### K-means with PCA (dimensionality reduction) (with labels)"""

x = input_data
y = output_data.astype(np.int32)
pca = PCA(n_components=1000)
x_reduced = pd.DataFrame(pca.fit_transform(x))

x_reduced['output'] = y
x_reduced = x_reduced.to_numpy()

k_means_clf = create_trained_classifier(k)

labels_with_pca_including_output = k_means_clf.fit_predict(x_reduced)

dependent_metrics_pca_including_output = generate_metrics(y, labels_with_pca_including_output, x_reduced)

print(dependent_metrics_pca_including_output)

metrics_dict["Kmeans with PCA on raw data with output"] = dependent_metrics_pca_including_output

"""## Gaussian Mixture Models (GMM)

### GMM without class field
"""

from sklearn import mixture

x = input_data
y = output_data.astype(np.int32)

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x)
gmm_metrics = generate_metrics(y, labels, x)
print(gmm_metrics)
metrics_dict["GMM on Raw Data"] = gmm_metrics

"""### GMM with class field"""

# appending output to input class

x_df = pd.DataFrame(input_data)
x_df['output'] = y
x = x_df.to_numpy()

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x)

gmm_with_output_metrics = generate_metrics(y, labels, x)

metrics_dict["GMM on Raw data with Output"] = gmm_with_output_metrics

print(gmm_with_output_metrics)

"""### GMM with PCA (dimensionality reduction)"""

# running PCA before GMM

x = input_data
y = output_data.astype(np.int32)

# reducing to 1000 dimensions (1000 being arbitarary here)
pca = PCA(n_components=1000)

x_reduced = pca.fit_transform(x)

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x)

gmm_pca_metrics = generate_metrics(y, labels, x)
metrics_dict["GMM with PCA on raw data"] = gmm_pca_metrics
print(gmm_pca_metrics)

"""## Visulising evaluation metrics"""

# Converting the nested dict of metrics to a dataframe for convenient plotting

index_of_experiments = list(metrics_dict.keys())

index_of_metrics = list(metrics_dict[index_of_experiments[0]].keys())

metrics_df = pd.DataFrame({'experiment': index_of_experiments})
# need one name and list of corresponding values to initialize the df

for name in index_of_metrics:
  list_of_values = [metrics_dict[x][name] for x in index_of_experiments]
  metrics_df[name] = list_of_values

metrics_df.head(10)

sns.barplot('f1-score', 'experiment', data=metrics_df)

sns.barplot('sc', 'experiment', data=metrics_df)

sns.barplot('dbi', 'experiment', data=metrics_df)

"""# Research Question"""

def image_equalise_hist(image: np.ndarray): 
    """ Equalise the histogram of the image globally.

        Reference:
        - https://en.wikipedia.org/wiki/Histogram_equalization
        - 
    """
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    #  Equalize the histogram of the image
    image = equalizeHist(image)

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)


def image_thresholding(image: np.ndarray): 
    """ Threshold the image using the Otsu thresholding
        algorithm. 

        Reference: 
        - https://en.wikipedia.org/wiki/Thresholding_(image_processing)
        - https://en.wikipedia.org/wiki/Otsu%27s_method
        - https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_thresholding.html
    """
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    # Find threshold using Otsu filter
    threshold: float = filters.threshold_otsu(image)
    binary = image > threshold

    binary_image = np.where(image, binary, 0) * 255

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def image_local_enhance_contrast(image: np.ndarray):
    """ Locally enhance the contrast of the image

        Reference: 
        - https://en.wikipedia.org/wiki/Histogram_equalization
        - https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
        - https://scikit-image.org/docs/dev/api/skimage.filters.rank.html#enhance-contrast
    """
    
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    image = filters.rank.enhance_contrast(image, morphology.disk(2))

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def image_local_autolevel(image: np.ndarray):
    """ Stretch histogram to normalise levels

        Reference:
        - https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
        - https://scikit-image.org/docs/dev/api/skimage.filters.rank.html#autolevel
    """
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    image = filters.rank.autolevel(image, morphology.disk(4))

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def image_enhance(image: np.ndarray):
    #  Resize image to a shape of (48, 48)
    image = image_as_square(image)

    image = filters.rank.enhance_contrast(image, morphology.disk(2))
    image = filters.rank.autolevel(image, morphology.disk(4))

    #  Resize the iamge back to a shape of (2304, )
    return image_as_array(image)

def enhance_all_images(data): 
    return np.apply_along_axis(image_enhance, 1, input_data)

def display_image_enhancement_comparison(image: np.array):
    def add_plot_line(im, axis):
        ax2 = axis.twinx()
        sns.distplot(im.flatten(), hist=True, bins=30, rug=False, ax=axis, kde=False)
        sns.distplot(im.flatten(), hist=False, rug=False, ax=ax2, kde=True, kde_kws={
            "alpha": 0.8, "cumulative":
                True, 'color': 'r'})
        sns.despine(ax=axis, left=True, right=True)
        sns.despine(ax=ax2, left=True, right=True)
        axis.set_yticks([])
        ax2.set_yticks([])
        axis.set_xticks([0, 50, 100, 150, 200, 250])
        ax2.set_xticks([0, 50, 100, 150, 200, 250])
        axis.set_xlim([0, 255])
        ax2.set_xlim([0, 255])

    hist_equalized_image = image_equalise_hist(image)
    autoleveled_image = image_enhance(image)
    thres_img = image_thresholding(image)
    loc_enhance = image_local_enhance_contrast(image)
    loc_autolev = image_local_autolevel(image)

    widths = [1, 1, 1, 1, 1, 1]
    heights = [4, 1]
    spec = dict(width_ratios=widths, height_ratios=heights)
    fig, axes = plt.subplots(nrows=2, ncols=6, figsize=(16, 5), gridspec_kw=spec)

    axes[0, 0].imshow(image_as_square(image), cmap='gray')
    axes[0, 0].set_axis_off()
    axes[0, 0].set_title("Original")
    add_plot_line(image, axes[1, 0])

    axes[0, 1].imshow(image_as_square(hist_equalized_image), cmap='gray')
    axes[0, 1].set_axis_off()
    axes[0, 1].set_title("Histogram equalization")
    add_plot_line(hist_equalized_image, axes[1, 1])

    axes[0, 2].imshow(image_as_square(thres_img), cmap='gray')
    axes[0, 2].set_axis_off()
    axes[0, 2].set_title("Thresholding")
    add_plot_line(thres_img, axes[1, 2])

    axes[0, 3].imshow(image_as_square(loc_enhance), cmap='gray')
    axes[0, 3].set_axis_off()
    axes[0, 3].set_title("Local contrasting")
    add_plot_line(loc_enhance, axes[1, 3])

    axes[0, 4].imshow(image_as_square(loc_autolev), cmap='gray')
    axes[0, 4].set_axis_off()
    axes[0, 4].set_title("Local equalisation")
    add_plot_line(loc_autolev, axes[1, 4])

    axes[0, 5].imshow(image_as_square(autoleveled_image), cmap='gray')
    axes[0, 5].set_axis_off()
    add_plot_line(autoleveled_image, axes[1, 5])

    fig.tight_layout()
    plt.show()

"""## Comparing different image enhancements"""

display_image_enhancement_comparison(input_data[0])

"""## Comparing models"""

enhanced_input_data = enhance_all_images(input_data)

"""### Naive Bayes with normalization and image mutations"""

# Normalize the data
enhanced_input_data_normalized = np.zeros(enhanced_input_data.shape)

for label in range(10):
    x = enhanced_input_data[output_data.astype(np.int)==label]
    normalized = normalize(x)
    enhanced_input_data_normalized[output_data.astype(np.int)==label] = normalized

compare_class_distribution_with_normalized(input_data, output_data)

# Run naive bayes on data after image preprocessing
x = enhanced_input_data_normalized
y = output_data

# split the data into training and test sets using stratified random sampling
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=52, stratify=y)

# train data sets and use model for predictions 
model = GaussianNB()
model.fit(x_train, y_train)
y_predictions = model.predict(x_test)

print_outcomes(x_test, y_test, y_predictions)
plot_confusion_matrix(y_test, y_predictions)

"""### Naive Bayes feature selection with image mutations"""

vizualize_important_features(enhanced_input_data_normalized, output_data)

"""### K-Means with image mutations

#### Without labels
"""

x = enhanced_input_data
y = output_data.astype(np.int32)

k_means_clf = create_trained_classifier(k)
labels_with_enhanced_input = k_means_clf.fit_predict(x)

dependent_metrics_enhanced = generate_metrics(y, labels_with_enhanced_input, x)

print(dependent_metrics_enhanced)

metrics_dict["Enhanced"] = dependent_metrics_enhanced

"""#### With labels"""

enhanced_df = pd.DataFrame(enhanced_input_data)
enhanced_df['output'] = y
x_enhanced_with_labels = enhanced_df.to_numpy()

k_means_clf = create_trained_classifier(k)
labels_enhanced_including_actual = k_means_clf.fit_predict(x_enhanced_with_labels)

dependent_metrics_enhanced_with_actual = generate_metrics(y, labels_enhanced_including_actual, x_enhanced_with_labels)


print(dependent_metrics_enhanced_with_actual)

metrics_dict["Enhanced with output"] = dependent_metrics_enhanced_with_actual

"""#### With PCA"""

# Running PCA on the enhanced data set
x = enhanced_input_data
y = output_data.astype(np.int32)

pca = PCA(n_components=1000)
x_reduced = pca.fit_transform(x)

k_means_clf = create_trained_classifier(k)

labels = k_means_clf.fit_predict(x_reduced)

enhanced_pca_metrics = generate_metrics(y, labels, x)

print(enhanced_pca_metrics)
metrics_dict["Kmeans with PCA on enhanced data"] = enhanced_pca_metrics

"""### Gaussian Mixture Models with Image Mutations"""

# enhanced input
x = enhanced_input_data

gmm = mixture.GaussianMixture(n_components=10)
labels = gmm.fit_predict(x)
print("fit complete")
gmm_enhanced_metrics = generate_metrics(y, labels, x)
metrics_dict["GMM on Enhanced data"] = gmm_enhanced_metrics

print(gmm_enhanced_metrics)

"""#### with PCA"""

# GMM with PCA on enhanced
x = enhanced_input_data
y = output_data.astype(np.int32)

pca = PCA(n_components=1000)
x_reduced = pca.fit_transform(x)

gmm = mixture.GaussianMixture(n_components=10)

labels = gmm.fit_predict(x_reduced)

gmm_pca_enhanced = generate_metrics(y, labels, x_reduced)

print(gmm_pca_enhanced)

metrics_dict['GMM with PCA on enhanced data'] = gmm_pca_enhanced

"""### Comparing performance of image mutations on clustering algorithms"""

# Converting the nested dict of metrics to a dataframe for convenient plotting

index_of_experiments = list(metrics_dict.keys())

index_of_metrics = list(metrics_dict[index_of_experiments[0]].keys())

metrics_df = pd.DataFrame({'experiment': index_of_experiments})
# need one name and list of corresponding values to initialize the df

for name in index_of_metrics:
  list_of_values = [metrics_dict[x][name] for x in index_of_experiments]
  metrics_df[name] = list_of_values

metrics_df.head(10)

"""### Comparing performance of models before and after mutations"""

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

sns.set()
sns.set(style='whitegrid')
palette = dict(Yes="#69c597", No="#F08080")

# data sets: raw, enhanced
# models: naive full, naive class-best, bayes k2, bayes tan, kmeans, gmm

model_comparison = pd.DataFrame({
    'F-Score': pd.Series([
        0.562,
        0.86,
        0.618,
        0.880,
        0.249131,
        0.723776,
        0.247867,
        0.677946,
        0.219816,
        0.732641,
        0.847,
        0.945,
        0.842,
        0.935,
        0.200787,
        0.637211,
    ]),
    'Mutated': pd.Categorical([
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes',
        'No', 'Yes'
    ]),
    'Model': pd.Categorical([
        'Naive Bayes',
        'Naive Bayes',
        'Naive Bayes\n(Class Best)',
        'Naive Bayes\n(Class Best)',
        'K-Means',
        'K-Means',
        'K-Means\n(with PCA)',
        'K-Means\n(with PCA)',
        'Gaussian Mixture',
        'Gaussian Mixture',
        'TAN',
        'TAN',
        'K2 (max.\n2 parents)',
        'K2 (max.\n2 parents)',
        'Gaussian Mixture\n(with PCA)',
        'Gaussian Mixture\n(with PCA)'
    ])
})

fig, ax = plt.subplots(figsize=(11, 4))

sns.barplot(x="Model", y="F-Score", data=model_comparison, hue="Mutated", palette=palette, ax=ax)
sns.despine(left=True)
plt.yticks(np.arange(0, 11) / 10)
plt.ylim([0, 1])
plt.xlabel('')

fig.tight_layout()
plt.show()