# modules
from keras.layers import Input, Dense, Dropout
from keras.models import Model
from keras.datasets import mnist
from keras.models import Sequential, load_model
from keras.optimizers import RMSprop
from keras.callbacks import TensorBoard
from __future__ import print_function
from keras.utils import plot_model
from IPython.display import SVG
from keras.utils.vis_utils import model_to_dot
from sklearn import preprocessing
from keras import layers
from keras import initializers
from matplotlib import axes
from matplotlib import rc

import keras
import matplotlib.pyplot as plt
import numpy as np
import math
import pydot
import graphviz
import pandas as pd
import IPython

%matplotlib inline
font = {'family' : 'monospace',
        'weight' : 'bold',
        'size'   : 20}

rc('font', **font)

"""
Explanation: Title of Database: Wall-Following navigation task with mobile robot SCITOS-G5
The data were collected as the SCITOS G5 navigates through the room following the wall in a clockwise
direction, for 4 rounds. To navigate, the robot uses 24 ultrasound sensors arranged circularly around its "waist". 
The numbering of the ultrasound sensors starts at the front of the robot and increases in clockwise direction.
End of explanation
"""


# import
data_raw = pd.read_csv('data/sensor_readings_24.csv', sep=",", header=None)
data = data_raw.copy()

"""
Explanation: Import and basic data inspection
End of explanation
"""


data.head()

"""
Explanation: The dataframe consists of only positive values and the classes are encoded as strings in the variable with index 24
End of explanation
"""


df_tab = data_raw
df_tab[24] = df_tab[24].astype('category')
tab = pd.crosstab(index=df_tab[24], columns="frequency")
tab.index.name = 'Class/Direction'
tab/tab.sum()

"""
Explanation: Whats the distribution of the classes?
End of explanation
"""


mapping = {key: value for (key, value) in zip(data[24].unique(), range(len(data[24].unique())))}
print(mapping)
data.replace({24:mapping}, inplace=True)

data[24].unique()

"""
Explanation: The Move_Forward and the Sharp-Right-Turn Class combine nearly 80% of all observated classes. So it might happen, that the accuracy may still be high with around 75% although most of the features are eliminated.
Preprocessing
0. Mapping integer values to the classes.
End of explanation
"""


data_train = data.sample(frac=0.9, random_state=42)
data_val = data.drop(data_train.index)

df_x_train = data_train.iloc[:,:-1]
df_y_train = data_train.iloc[:,-1]

df_x_val = data_val.iloc[:,:-1]
df_y_val = data_val.iloc[:,-1]

"""
Explanation: 1. Take a random sample of 90% of the rows from the dataframe. To ensure reproducability the random_state variable is set. The other 10% are placed aside for validation after training. The last column is the class column and is stored in the y variables respectively.
End of explanation
"""


x_train = df_x_train.values
x_train = (x_train - x_train.min()) / (x_train.max() - x_train.min())
y_train = df_y_train.values

x_val = df_x_val.values
x_val = (x_val - x_val.min()) / (x_val.max() - x_val.min())
y_val = df_y_val.values
y_eval = y_val

"""
Explanation: 2. Normalization between 0 and 1
End of explanation
"""


y_train = keras.utils.to_categorical(y_train, 4)
y_val = keras.utils.to_categorical(y_val, 4)

"""
Explanation: 3. Make useful categorical variables out of the single column data by one-hot encoding it.
End of explanation
"""


epochsize = 150
batchsize = 24
shuffle = False
dropout = 0.1
num_classes = 4
input_dim = x_train.shape[1]
hidden1_dim = 30
hidden2_dim = 30
class_names = mapping.keys()

"""
Explanation: 4. Set Global Parameters
End of explanation
"""


input_data = Input(shape=(input_dim,), dtype='float32', name='main_input')
hidden_layer1 = Dense(hidden1_dim, activation='relu', input_shape=(input_dim,), kernel_initializer='normal')(input_data)
dropout1 = Dropout(dropout)(hidden_layer1)
hidden_layer2 = Dense(hidden2_dim, activation='relu', input_shape=(input_dim,), kernel_initializer='normal')(dropout1)
dropout2 = Dropout(dropout)(hidden_layer2)
output_layer = Dense(num_classes, activation='softmax', kernel_initializer='normal')(dropout2)

model = Model(inputs=input_data, outputs=output_layer)

model.compile(loss='categorical_crossentropy',
              optimizer=RMSprop(),
              metrics=['accuracy'])

plot_model(model, to_file='images/robo1_nn.png', show_shapes=True, show_layer_names=True)

IPython.display.Image("images/robo1_nn.png")

model.fit(x_train, y_train, 
          batch_size=batchsize,
          epochs=epochsize,
          verbose=0,
          shuffle=shuffle,
          validation_split=0.05)
nn_score = model.evaluate(x_val, y_val)[1]
print(nn_score)

fig = plt.figure(figsize=(20,10))
plt.plot(model.history.history['val_acc'])
plt.plot(model.history.history['acc'])
plt.axhline(y=nn_score, c="red")
plt.text(0, nn_score, "test: " + str(round(nn_score, 4)), fontdict=font)
plt.title('model accuracy for neural net with 2 hidden layers')
plt.ylabel('accuracy')
plt.xlabel('epochs')
plt.legend(['valid', 'train'], loc='lower right')
plt.show()

import itertools

from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix

def plot_confusion_matrix(cm, classes,
                          normalize=False,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    print(cm)

    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45)
    plt.yticks(tick_marks, classes)

    fmt = '.2f' if normalize else 'd'
    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, format(cm[i, j], fmt),
                 horizontalalignment="center",
                 color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')

# Compute confusion matrix
cnf_matrix = confusion_matrix(y_eval, model.predict(x_val).argmax(axis=-1))
np.set_printoptions(precision=2)

# Plot non-normalized confusion matrix
plt.figure(figsize=(20,10))
plot_confusion_matrix(cnf_matrix, classes=class_names,
                      title='Confusion matrix, without normalization')

# Plot normalized confusion matrix
plt.figure(figsize=(20,10))
plot_confusion_matrix(cnf_matrix, classes=class_names, normalize=True,
                      title='Normalized confusion matrix')

"""
Explanation: Train Neural Net
Due to a tight schedule we will not perform any cross validation. So it might happen that our accuracy estimators lack a little bit in potential of generalization. We shall live with that. Another setup of experiments would be, that we loop over some different dataframes samples up in the preprocessing steps and repeat all the steps below to finally average the results.
The dimension of the hidden layers are set arbitrarily but some runs have shown that 30 is a good number. The input_dim Variable is set to 24 because initially there are 24 features. The aim is to build the best possible neural net.
Optimizer
RMSprop is a mini batch gradient descent algorithm which divides the gradient by a running average of the learning rate. More information: http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf
The weights are initialized by a normal distribution with mean 0 and standard deviation of 0.05.
End of explanation
"""


IPython.display.Image("images/2018-01-25 18_44_01-PubMed Central, Table 2_ Sensors (Basel). 2017 Mar; 17(3)_ 549. Published online.png")

"""
Explanation: Comparison
The following data is from a paper published in March 2017. You can find that here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5375835/
End of explanation
"""


encoder_dim = 8
hidden1_dim = 30
hidden2_dim = 30

main_input = Input(shape=(input_dim,), dtype='float32', name='main_input')

encoding_layer = Dense(encoder_dim, activation='relu', kernel_initializer='normal')
encoding_layer_output = encoding_layer(main_input)
decoding_layer_output = Dense(input_dim
                              ,activation='sigmoid'
                              ,name='decoder_output'
                              ,kernel_initializer='normal')(encoding_layer_output)

x = Dense(hidden1_dim, activation='relu', kernel_initializer='normal')(encoding_layer_output)
x = Dropout(dropout)(x)
x = Dense(hidden2_dim, activation='relu', kernel_initializer='normal')(x)
x = Dropout(dropout)(x)

classifier_output = Dense(num_classes
                          ,activation='softmax'
                          ,name='main_output'
                          ,kernel_initializer='normal')(x)

auto_classifier = Model(inputs=main_input, outputs=[classifier_output, decoding_layer_output])

auto_classifier.compile(optimizer=RMSprop(),
                        loss={'main_output': 'categorical_crossentropy', 'decoder_output': 'mean_squared_error'},
                        loss_weights={'main_output': 1., 'decoder_output': 1.},
                        metrics=['accuracy'])

plot_model(auto_classifier, to_file='images/robo4_auto_class_LR.png', show_shapes=True, show_layer_names=True)

IPython.display.Image("images/robo4_auto_class_LR.png")

auto_classifier.fit({'main_input': x_train},
                    {'main_output': y_train, 'decoder_output': x_train},
                    epochs=epochsize, 
                    batch_size=batchsize,
                    shuffle=shuffle,
                    validation_split=0.05,
                    verbose=0)

score = auto_classifier.evaluate(x=x_val, y=[y_val, x_val], verbose=1)[3]
print(score)

fig = plt.figure(figsize=(20,10))
plt.plot(auto_classifier.history.history['val_main_output_acc'])
plt.plot(auto_classifier.history.history['main_output_acc'])
plt.axhline(y=score, c="red")
plt.text(0, score, "test: " + str(round(score, 4)), fontdict=font)
plt.title('model accuracy for ' + str(round(input_dim/encoder_dim, 2)) +  ' x compression with single layer autoencoder')
plt.ylabel('accuracy')
plt.xlabel('epochs')
plt.legend(['valid', 'train'], loc='lower right')
plt.show()

"""
Explanation: One can easily see that our results are better.
Dimensionality Reduction with Single Hidden Layer Autoencoder
End of explanation
"""


encoder_dim1 = 16
encoder_dim2 = 8
decoder_dim1 = 16

main_input = Input(shape=(input_dim,), dtype='float32', name='main_input')

encoding_layer1 = Dense(encoder_dim1, activation='relu', kernel_initializer='normal')(main_input)
encoding_layer2 = Dense(encoder_dim2, activation='relu', kernel_initializer='normal')(encoding_layer1)

decoding_layer1 = Dense(decoder_dim1
                        ,activation='relu'
                        ,kernel_initializer='normal')(encoding_layer2)
decoding_layer2 = Dense(input_dim
                        ,activation='sigmoid'
                        ,name='decoder_output'
                        ,kernel_initializer='normal')(decoding_layer1)

x = Dense(hidden1_dim, activation='relu', kernel_initializer='normal')(encoding_layer2)
x = Dropout(dropout)(x)
x = Dense(hidden2_dim, activation='relu', kernel_initializer='normal')(x)
x = Dropout(dropout)(x)

classifier_output = Dense(num_classes
                          ,activation='softmax'
                          ,name='main_output'
                          ,kernel_initializer='normal')(x)

stacked_auto_classifier = Model(inputs=main_input, outputs=[classifier_output, decoding_layer2])

stacked_auto_classifier.compile(optimizer=RMSprop(),
                        loss={'main_output': 'categorical_crossentropy', 'decoder_output': 'mean_squared_error'},
                        loss_weights={'main_output': 1., 'decoder_output': 1.},
                        metrics=['accuracy'])

plot_model(stacked_auto_classifier, to_file='images/stacked__auto_class.png', show_shapes=True, show_layer_names=True)

IPython.display.Image("images/stacked__auto_class.png")

stacked_auto_classifier.fit({'main_input': x_train},
                            {'main_output': y_train, 'decoder_output': x_train},
                            epochs=epochsize, 
                            batch_size=batchsize,
                            shuffle=shuffle,
                            validation_split=0.05,
                            verbose=0)

stacked_score = stacked_auto_classifier.evaluate(x=x_val, y=[y_val, x_val], verbose=1)[3]
print(stacked_score)

fig = plt.figure(figsize=(20,10))
plt.plot(stacked_auto_classifier.history.history['val_main_output_acc'])
plt.plot(stacked_auto_classifier.history.history['main_output_acc'])
plt.axhline(y=stacked_score, c="red")
plt.text(0, stacked_score, "test: " + str(round(stacked_score, 4)), fontdict=font)
plt.title('model accuracy for ' + str(round(input_dim/encoder_dim, 2)) +  ' x compression with stacked autoencoder')
plt.ylabel('accuracy')
plt.xlabel('epochs')
plt.legend(['valid', 'train'], loc='lower right')
plt.show()

"""
Explanation: Dimensionality reduction with stacked (multi hidden layer) autoencoder
End of explanation
"""


# the initial coding dimension s.t. there is no feature selection at the beginning
encoding_dim = input_dim
result3 = {'encoding_dim': []
           ,'auto_classifier_acc': []}

while encoding_dim > 0:
    main_input = Input(shape=(input_dim,), dtype='float32', name='main_input')

    encoding_layer = Dense(encoding_dim, activation='relu', name='encoder', kernel_initializer='normal')
    encoding_layer_output = encoding_layer(main_input)
    decoding_layer_output = Dense(input_dim, activation='sigmoid'
                                  ,name='decoder_output'
                                  ,kernel_initializer='normal')(encoding_layer_output)

    x = Dense(hidden1_dim, activation='relu', kernel_initializer='normal')(encoding_layer_output)
    x = Dropout(dropout)(x)
    x = Dense(hidden2_dim, activation='relu', kernel_initializer='normal')(x)
    x = Dropout(dropout)(x)

    classifier_output = Dense(num_classes, activation='softmax', name='main_output', kernel_initializer='normal')(x)

    auto_classifier = Model(inputs=main_input, outputs=[classifier_output, decoding_layer_output])

    auto_classifier.compile(optimizer=RMSprop(),
                            loss={'main_output': 'categorical_crossentropy', 'decoder_output': 'mean_squared_error'},
                            loss_weights={'main_output': 1., 'decoder_output': 1.},
                            metrics=['accuracy'])

    auto_classifier.fit({'main_input': x_train},
                        {'main_output': y_train, 'decoder_output': x_train},
                        epochs=epochsize, 
                        batch_size=batchsize,
                        shuffle=shuffle,
                        validation_split=0.05,
                        verbose=0)

    accuracy = auto_classifier.evaluate(x=x_val, y=[y_val, x_val], verbose=1)[3]
    result3['encoding_dim'].append(encoding_dim)
    result3['auto_classifier_acc'].append(accuracy)

    encoding_dim -=1
    print(result3)

result_df = pd.DataFrame(result3)
result_df['neural_net_acc'] = nn_score
result_df

fig = plt.figure(figsize=(20,10))
plt.bar(result_df['encoding_dim'], result_df['auto_classifier_acc'])
plt.axhline(y=result_df['neural_net_acc'][0], c="red")
plt.text(0, result_df['neural_net_acc'][0], "best neural net: " + str(round(result_df['neural_net_acc'][0], 4))
         ,fontdict=font)
plt.title('model accuracy for different encoding dimensions')
plt.ylabel('accuracy')
plt.xlabel('dimension')
plt.ylim(0.6, 1)

"""
Explanation: Whats the best dimensionality reduction with single autoencoder?
End of explanation
"""


# Compute confusion matrix
cnf_matrix = confusion_matrix(y_eval, auto_classifier.predict(x_val)[0].argmax(axis=-1))
np.set_printoptions(precision=2)

# Plot non-normalized confusion matrix
plt.figure(figsize=(20,10))
plot_confusion_matrix(cnf_matrix, classes=class_names,
                      title='Confusion matrix, without normalization')

# Plot normalized confusion matrix
plt.figure(figsize=(20,10))
plot_confusion_matrix(cnf_matrix, classes=class_names, normalize=True,
                      title='Normalized confusion matrix')

result_df.to_csv('results/robo_results.csv')

"""
Explanation: Prediction for a classifier with a dimension of 16 (not for report)
End of explanation
"""


encoding_weights = encoding_layer.get_weights()
sum_of_weights = {index: item.sum() for (index, item) in enumerate(encoding_weights[0])}

weights = sum_of_weights
features = []
for i in range(encoder_dim1):
    max_key = max(weights, key=lambda key: weights[key])
    features.append(max_key)
    del weights[max_key]
    
print(features)

"""
Explanation: Experimental area (not for presentation)
The indices of the extracted features are:
End of explanation
"""


x_train_selected = np.array([x[features] for x in x_train])
x_val_selected = np.array([x[features] for x in x_val])

input_dim = x_train_selected.shape[1]
hidden1_dim = 26
hidden2_dim = 26

result3 = []

for i in range(1,4):
    model_new = Sequential()
    model_new.add(Dense(hidden1_dim, activation='relu', input_shape=(input_dim,)))
    model_new.add(Dense(hidden2_dim, activation='relu'))
    model_new.add(Dense(num_classes, activation='softmax'))

    model_new.compile(loss='categorical_crossentropy',
                  optimizer=RMSprop(),
                  metrics=['accuracy'])

    model_new.fit(x_train_selected, y_train, 
              batch_size=batchsize,
              epochs=epochsize,
              verbose=0,
              shuffle=shuffle,
              validation_split=0.1)
    score = model_new.evaluate(x_val_selected, y_val)[1]
    result3.append(score)
    print(result3)
print(np.mean(result3))

"""
Explanation: Check the neural net performance with new selected features
End of explanation
"""


# the initial coding dimension s.t. there is no feature selection at the beginning
encoding_dim = 24
# dimension of the neural net layer1
hidden1_dim = 30
# dimension of the second neural net layer
hidden2_dim = 30
epoch_size = 150
batch_size = 24
shuffle = False
result2 = {'encoding_dim/features': []
           ,'compression_level': []
           ,'auto_classifier_acc': []
           ,'selected_classifier_acc': []
           ,'features': []}

"""
Explanation: Finding good features
"Is there a good number of features for the Robo- Dataset?" For that we create a loop over some numbers (i.e. dimension of hidden layer from the single autoencoder) and get the respective result.
The architecture within those loops is quite special. For that we recommend to have a look at the neural net visualization below.
The parameters for the ensemble training
End of explanation
"""


while encoding_dim > 0:
    main_input = Input(shape=(input_dim,), dtype='float32', name='main_input')

    encoding_layer = Dense(encoding_dim, activation='relu', name='encoder', kernel_initializer='normal')
    encoding_layer_output = encoding_layer(main_input)
    decoding_layer_output = Dense(input_dim, activation='sigmoid'
                                  ,name='decoder_output'
                                  ,kernel_initializer='normal')(encoding_layer_output)

    x = Dense(hidden1_dim, activation='relu', kernel_initializer='normal')(encoding_layer_output)
    x = Dense(hidden2_dim, activation='relu', kernel_initializer='normal')(x)

    classifier_output = Dense(num_classes, activation='softmax', name='main_output', kernel_initializer='normal')(x)

    auto_classifier = Model(inputs=main_input, outputs=[classifier_output, decoding_layer_output])

    auto_classifier.compile(optimizer=RMSprop(),
                            loss={'main_output': 'categorical_crossentropy', 'decoder_output': 'mean_squared_error'},
                            loss_weights={'main_output': 1., 'decoder_output': 1.},
                            metrics=['accuracy'])

    auto_classifier.fit({'main_input': x_train},
                        {'main_output': y_train, 'decoder_output': x_train},
                        epochs=epoch_size, 
                        batch_size=batch_size,
                        shuffle=shuffle,
                        validation_split=0.1,
                        verbose=0)

    accuracy = auto_classifier.evaluate(x=x_val, y=[y_val, x_val], verbose=1)[3]
    result2['encoding_dim/features'].append(encoding_dim)
    result2['compression_level'].append(1 - encoding_dim/24)
    result2['auto_classifier_acc'].append(accuracy)

    encoding_weights = encoding_layer.get_weights()
    sum_of_weights = {index: item.sum() for (index, item) in enumerate(encoding_weights[0])}

    weights = sum_of_weights
    features = []
    for i in range(encoding_dim):
        max_key = max(weights, key=lambda key: weights[key])
        features.append(max_key)
        del weights[max_key]
    
    result2['features'].append(features)
    x_train_selected = np.array([x[features] for x in x_train])
    x_val_selected = np.array([x[features] for x in x_val])

    input_dim_new = x_train_selected.shape[1]
    
    accuracy_list = []
    for i in range(1):
        model_new = Sequential()
        model_new.add(Dense(hidden1_dim, activation='relu', input_shape=(input_dim_new,), kernel_initializer='normal'))
        model_new.add(Dense(hidden2_dim, activation='relu', kernel_initializer='normal'))
        model_new.add(Dense(num_classes, activation='softmax', kernel_initializer='normal'))

        model_new.compile(loss='categorical_crossentropy',
                      optimizer=RMSprop(),
                      metrics=['accuracy'])

        model_new.fit(x_train_selected, y_train, 
                  batch_size=batch_size,
                  epochs=epoch_size,
                  verbose=0,
                  shuffle=shuffle,
                  validation_split=0.1)
        score = model_new.evaluate(x_val_selected, y_val)[1]
        accuracy_list.append(score)

    result2['selected_classifier_acc'].append(np.mean(accuracy_list))
    encoding_dim -=1
    print(result2)

"""
Explanation: The idea is:
- reduce the number of features encoding_dim iteratively by training an auto_classifier network
- after training, the weights of the first layer are analyzed.
- therefore sum up the weights for each input node and take the encoding_dim number as new features
- test the new features with the pre- defined neural net and see which selection is best.
End of explanation
"""


auto_classifier.summary()

"""
Explanation: Summary
End of explanation
"""


plot_model(auto_classifier, to_file='images/robo2_auto_class_LR.png', show_shapes=True, show_layer_names=True, rankdir='LR')

IPython.display.Image("images/robo2_auto_class_LR.png")

"""
Explanation: Architecture
End of explanation
"""


result_df = pd.DataFrame(result2)
result_df['neural_net_acc'] = 0.949938949939
result_df

result_df.to_csv('results/robo_results.csv')

result_df.plot(x='encoding_dim/features', y=['selected_classifier_acc', 'neural_net_acc'], kind='bar', figsize=(20,10))

"""
Explanation: Results
End of explanation
"""


