{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Siamese Network with KNN outliers detector","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow import keras\nfrom keras.layers import Input, Conv2D, MaxPooling2D, Flatten, Dense, Lambda, BatchNormalization, Activation, AveragePooling2D, Dropout\nfrom keras.models import Sequential, Model\nfrom tensorflow.keras.optimizers import Adam, SGD\nimport keras.backend as K\nimport cv2\nimport numpy as np\nimport random\nfrom scipy.io import wavfile\nfrom scipy import signal\nimport matplotlib.pyplot as plt\nimport os\nimport numpy as np\nimport pandas as pd\nimport librosa\nimport librosa.display\nfrom tensorflow.keras import layers\nfrom sklearn.model_selection import train_test_split\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2022-05-24T16:06:22.056784Z","iopub.execute_input":"2022-05-24T16:06:22.057197Z","iopub.status.idle":"2022-05-24T16:06:28.958131Z","shell.execute_reply.started":"2022-05-24T16:06:22.057089Z","shell.execute_reply":"2022-05-24T16:06:28.956945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install librosa\n!pip install pyod","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dev_data_path =  \"../input/eurecom-aml-2022-challenge-2/dev_data/dev_data/slider/train/\"\neval_data_path =  \"../input/eurecom-aml-2022-challenge-2/eval_data/eval_data/slider/train/\"\ndev_files = os.listdir(dev_data_path)\neval_files = os.listdir(eval_data_path)\n# add root path\ndev_files = [dev_data_path + file for file in dev_files]\neval_files = [eval_data_path + file for file in eval_files]\nprint(f\"There are {len(dev_files)} dev train audio\")\nprint(f\"There are {len(eval_files)} eval train audio\")\n\n# merge in a unique train\ntrain_files = [*dev_files, *eval_files]\n\n# get test set (from training set) - This will be used for assessing performances with AUC score\ntest_data_path =  \"../input/eurecom-aml-2022-challenge-2/dev_data/dev_data/slider/test/\"\ntest_files = os.listdir(test_data_path)\ntest_files = [test_data_path + file for file in test_files]\nprint(f\"There are {len(test_files)} dev test audio\")\n\n# get final test set (the one of the leaderboard)\nfinal_test_data_path =  \"../input/eurecom-aml-2022-challenge-2/eval_data/eval_data/slider/test/\"\nfinal_test_files = os.listdir(final_test_data_path)\nfinal_test_files = [final_test_data_path + file for file in final_test_files]\nprint(f\"There are {len(final_test_files)} final test audio\")","metadata":{"execution":{"iopub.status.busy":"2022-05-24T16:32:32.682673Z","iopub.execute_input":"2022-05-24T16:32:32.683232Z","iopub.status.idle":"2022-05-24T16:32:32.706353Z","shell.execute_reply.started":"2022-05-24T16:32:32.683101Z","shell.execute_reply":"2022-05-24T16:32:32.705619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train Data preprocessing","metadata":{}},{"cell_type":"code","source":"import warnings\n# Used to suppress warning caused by deprecated API of NumPy version 1.19 or higher\nwarnings.filterwarnings(\"ignore\", category=np.VisibleDeprecationWarning) \n\nwav_files = np.array([librosa.load(file, sr=None) for file in train_files]) \nsample_rates = np.array(list(map(lambda tup: tup[1], wav_files)))\ntrain_audios = np.array(list(map(lambda tup: tup[0], wav_files)))\n\nwav_files = np.array([librosa.load(file, sr=None) for file in test_files])\nsample_rates = np.array(list(map(lambda tup: tup[1], wav_files)))\ntest_audios = np.array(list(map(lambda tup: tup[0], wav_files)))\n\nwav_files = np.array([librosa.load(file, sr=None) for file in final_test_files])\nsample_rates = np.array(list(map(lambda tup: tup[1], wav_files)))\nfinal_test_audios = np.array(list(map(lambda tup: tup[0], wav_files)))","metadata":{"execution":{"iopub.status.busy":"2022-05-24T16:08:35.497468Z","iopub.execute_input":"2022-05-24T16:08:35.498102Z","iopub.status.idle":"2022-05-24T16:09:48.685766Z","shell.execute_reply.started":"2022-05-24T16:08:35.498066Z","shell.execute_reply":"2022-05-24T16:09:48.684955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12,7))\nlibrosa.display.waveshow(train_audios[9], sr=16000)\nplt.title(\"Slider machine sound wave\", fontsize=25, pad=20)\nplt.xticks(fontsize=15)\nplt.yticks(fontsize=15)\nplt.xlabel(\"Time\", fontsize=20, labelpad=15)\nplt.ylabel(\"Intensity\", fontsize=20, labelpad=15)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-24T16:32:41.397809Z","iopub.execute_input":"2022-05-24T16:32:41.399304Z","iopub.status.idle":"2022-05-24T16:32:41.403207Z","shell.execute_reply.started":"2022-05-24T16:32:41.39926Z","shell.execute_reply":"2022-05-24T16:32:41.402165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get spectrograms\nsample_rate = sample_rates[0]\n\nX_train_val = librosa.power_to_db(librosa.feature.melspectrogram(y=train_audios,\n                      sr=sample_rate,\n                      n_fft=1024,\n                      hop_length=417,\n                      n_mels=384,\n))\n\nX_test = librosa.power_to_db(librosa.feature.melspectrogram(y=test_audios,\n                      sr=sample_rate,\n                      n_fft=1024,\n                      hop_length=417,\n                      n_mels=384,\n))\n\nX_test_final = librosa.power_to_db(librosa.feature.melspectrogram(y=final_test_audios,\n                      sr=sample_rate,\n                      n_fft=1024,\n                      hop_length=417,\n                      n_mels=384,\n))","metadata":{"execution":{"iopub.status.busy":"2022-05-24T16:10:22.690657Z","iopub.execute_input":"2022-05-24T16:10:22.690917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# display a spectrogram\nplt.figure(figsize=(15,5))\nplt.imshow(X_train_val[0])\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# resize image for Siamese Neural Network - 128x128 input for a fast training\nfrom skimage.transform import resize\n\nX_train_val = resize(X_train_val, (X_train_val.shape[0],128,128))\nX_test = resize(X_test, (X_test.shape[0],128,128))\nX_test_final = resize(X_test_final, (X_test_final.shape[0],128,128))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_machine_id_from_file_name(file_name):\n    # Siamese Network solves a supervised problem, hence we use machine IDs as labels in the first part of our work\n    return int(file_name.split(\"/\")[-1].split(\"_\")[2])\n\ndef get_labels(files):\n    labels = []\n    for file in files:\n    labels.append(int(get_machine_id_from_file_name(file)))\n    return np.array(labels)\n\ny_train_val = get_labels(train_files)\ny_test = get_labels(test_files)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize training slider machine IDs distribution\n\nplt.figure(figsize=(12,7))\nsns.countplot(x=y_train_val)\nplt.title(\"Training machine IDs distribution\", fontsize=25, pad=20)\nplt.xticks(fontsize=15)\nplt.yticks(fontsize=15)\nplt.xlabel(\"Machine ID\", fontsize=20, labelpad=15)\nplt.ylabel(\"Count\", fontsize=20, labelpad=15)\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use 20% of train_val as validation set\nX_train, X_val, y_train, y_val = train_test_split(X_train_val, y_train_val, train_size=0.8, random_state=42, stratify=y_train_val)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Normalize\ndef normalize_sample(sample):\n    return (sample - sample.mean()) / sample.std()\n\ndef norm_set(X):\n    for i, element in enumerate(X):\n        X[i] = normalize_sample(X[i])\n    return X\n\nX_train = norm_set(X_train)\nX_val = norm_set(X_val)\nX_test = norm_set(X_test)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generate pairs for input to Siamese Network. \n# If ID is either 4 or 5 generate double pairs since those are the half of other machine IDs (we balance the sets)\ndef make_pairs(x, y, test=False):\n    # num classes, pay attention to class 0 (Slider machine IDs range from 0 to 5)\n    num_classes = max(y) + 1\n    # list of ndarrays, each containing the indices of a specific class in the dataset\n    digit_indices = [np.where(y == i)[0] for i in range(num_classes)]\n    pairs = []\n    labels = []\n\n    for idx1 in range(len(x)):\n        # add 5 matching examples - if classes 4 or 5 add 10 matching examples\n        if idx1 == 4 or idx1 == 5:\n            num_pairs = 10\n        else:\n            num_pairs = 5\n\n        x1 = x[idx1]\n        label1 = y[idx1]\n        for j in range(num_pairs):\n            idx2 = random.choice(digit_indices[label1])\n            x2 = x[idx2]\n          \n            pairs += [[x1, x2]]  # add positive pair\n            labels += [1]\n\n        # add 5 non-matching example (or 10 if classes are either 4 or 5)\n        for j in range(num_pairs):\n            label2 = random.randint(0, num_classes - 1)\n            if test:  # in case it is test we have to skip classes 1,3,5 since dev_test contains only 0,2,4 machine IDs\n                while label2 == label1 or label2 in [1,3,5]:\n                    label2 = random.randint(0, num_classes - 1)\n            else:\n                while label2 == label1:\n                    label2 = random.randint(0, num_classes - 1)\n\n            idx2 = random.choice(digit_indices[label2])\n            x2 = x[idx2]\n\n            pairs += [[x1, x2]]  # add negative pair\n            labels += [0]\n\n    return np.array(pairs), np.array(labels).astype(\"float32\")\n\n\n# make train pairs\npairs_train, labels_train = make_pairs(X_train, y_train)\n# make validation pairs\npairs_val, labels_val = make_pairs(X_val, y_val)\n# make test pairs\npairs_test, labels_test = make_pairs(X_test, y_test, test=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split train, validation and test pairs\nx_train_1 = pairs_train[:, 0]\nx_train_2 = pairs_train[:, 1]\n\nx_val_1 = pairs_val[:, 0]\nx_val_2 = pairs_val[:, 1]\n\nx_test_1 = pairs_test[:, 0]\nx_test_2 = pairs_test[:, 1]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize pairs if wanted\ndef visualize(pairs, labels, to_show=6, num_col=3, predictions=None, test=False):\n    \n    num_row = to_show // num_col if to_show // num_col != 0 else 1\n    to_show = num_row * num_col\n\n    # Plot the images\n    fig, axes = plt.subplots(num_row, num_col, figsize=(5, 5))\n    for i in range(to_show):\n        if num_row == 1:\n            ax = axes[i % num_col]\n        else:\n            ax = axes[i // num_col, i % num_col]\n\n        ax.imshow(tf.concat([pairs[i][0], pairs[i][1]], axis=1), cmap=\"gray\")\n        ax.set_axis_off()\n        if test:\n            ax.set_title(\"True: {} | Pred: {:.5f}\".format(labels[i], predictions[i][0]))\n        else:\n            ax.set_title(\"Label: {}\".format(labels[i]))\n    if test:\n        plt.tight_layout(rect=(0, 0, 1.9, 1.9), w_pad=0.0)\n    else:\n        plt.tight_layout(rect=(0, 0, 1.5, 1.5))\n    plt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize(pairs_train[:-1], labels_train[:-1], to_show=8, num_col=4)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Siamese Network with euclidean distance\ndef euclidean_distance(vects):\n    x, y = vects\n    sum_square = tf.math.reduce_sum(tf.math.square(x - y), axis=1, keepdims=True)\n    return tf.math.sqrt(tf.math.maximum(sum_square, tf.keras.backend.epsilon()))\n\n# pretty fast network, few layers\ninput = layers.Input((128, 128, 1))\nx = tf.keras.layers.BatchNormalization()(input)\nx = layers.Conv2D(4, (5, 5), activation=\"tanh\")(x)\nx = layers.AveragePooling2D(pool_size=(2, 2))(x)\nx = layers.Dropout(0.5)(x)\nx = layers.Conv2D(16, (5, 5), activation=\"tanh\")(x)\nx = layers.AveragePooling2D(pool_size=(2, 2))(x)\nx = layers.Dropout(0.5)(x)\nx = layers.Flatten()(x)\n\nx = tf.keras.layers.BatchNormalization()(x)\nx = layers.Dense(64, activation=\"tanh\")(x)  # 64 features that will be used as input to KNN outlier detector\nembedding_network = keras.Model(input, x)\n\n\ninput_1 = layers.Input((128, 128, 1))\ninput_2 = layers.Input((128, 128, 1))\n\ntower_1 = embedding_network(input_1)\ntower_2 = embedding_network(input_2)\n\nmerge_layer = layers.Lambda(euclidean_distance)([tower_1, tower_2])\nnormal_layer = tf.keras.layers.BatchNormalization()(merge_layer)\noutput_layer = layers.Dense(1, activation=\"sigmoid\")(normal_layer)\nsiamese = keras.Model(inputs=[input_1, input_2], outputs=output_layer)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def loss(margin=1):\n    # contrastive loss with margin 1 by default\n    def contrastive_loss(y_true, y_pred):\n        square_pred = tf.math.square(y_pred)\n        margin_square = tf.math.square(tf.math.maximum(margin - (y_pred), 0))\n        return tf.math.reduce_mean(\n            (1 - y_true) * square_pred + (y_true) * margin_square\n        )\n\n    return contrastive_loss","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"opt = SGD(lr=0.05, momentum=0.7)\nsiamese.compile(loss=loss(margin=1), optimizer=opt, metrics=[\"accuracy\"])\nsiamese.summary()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = siamese.fit(\n    [x_train_1, x_train_2],\n    labels_train,\n    validation_data=([x_val_1, x_val_2], labels_val),\n    batch_size=128,\n    epochs=35,\n)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# function to plot history of loss and accuracy - training and validation\ndef plot_history(losses, val_losses, accuracies, val_accuracies):\n    plt.plot(losses)\n    plt.plot(val_losses)\n    plt.legend([\"train_loss\", \"val_loss\"])\n    plt.xlabel(\"Epochs\")\n    plt.ylabel(\"Loss\")\n    plt.show()\n\n    plt.plot(accuracies)\n    plt.plot(val_accuracies)\n    plt.legend([\"train_accuracy\", \"val_accuracy\"])\n    plt.xlabel(\"Epochs\")\n    plt.ylabel(\"Accuracy\")\n    plt.show()\n\nplot_history(\n    history.history[\"loss\"],\n    history.history[\"val_loss\"],\n    history.history[\"accuracy\"],\n    history.history[\"val_accuracy\"],\n)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract encoder which is used as feature extractor for the KNN outlier detector\ninput = layers.Input((128, 128, 1))\nx = tf.keras.layers.BatchNormalization()(input)\nx = layers.Conv2D(4, (5, 5), activation=\"tanh\")(x)\nx = layers.AveragePooling2D(pool_size=(2, 2))(x)\nx = layers.Dropout(0.5)(x)\nx = layers.Conv2D(16, (5, 5), activation=\"tanh\")(x)\nx = layers.AveragePooling2D(pool_size=(2, 2))(x)\nx = layers.Dropout(0.5)(x)\nx = layers.Flatten()(x)\nx = tf.keras.layers.BatchNormalization()(x)\nx = layers.Dense(64, activation=\"tanh\")(x)\nembedding_network = keras.Model(input, x)\n\n# set weights of the trained Siamese network encoder with weights from training\nembedding_network.layers[1].set_weights(siamese.layers[2].get_weights()[:4])  # BN params\nembedding_network.layers[2].set_weights(siamese.layers[2].get_weights()[4:6])  # conv2d weights + biases\nembedding_network.layers[5].set_weights(siamese.layers[2].get_weights()[6:8])  # conv2d weights + biases\nembedding_network.layers[9].set_weights(siamese.layers[2].get_weights()[8:12])  # BN params\nembedding_network.layers[10].set_weights(siamese.layers[2].get_weights()[12:])  # Dense matrix + biases","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now extract features to feed to KNN outliers detector\nX_train = np.concatenate((X_train, X_val), axis=0)\nfeat_vect_train = embedding_network.predict(X_train)\nfeat_vect_test = embedding_network.predict(X_test)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# If you want to apply PCA\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA \n\nsc = StandardScaler()\nfeat_vect_train_sc = sc.fit_transform(feat_vect_train)\nfeat_vect_test_sc = sc.transform(feat_vect_test)\npca = PCA(0.7)\nfeat_vect_train_pca = pca.fit_transform(feat_vect_train_sc)\nfeat_vect_test_pca = pca.transform(feat_vect_test_sc)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# KNN outliers detector\nfrom pyod.models.knn import KNN\n\nknn = KNN(contamination=0.01)\nknn.fit(feat_vect_train_pca)\ny_train_scores = knn.decision_scores_","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(y_train_scores, bins='auto')\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# compute AUC score on development test set\nfrom sklearn.metrics import roc_auc_score\n\n# Extract test labels\ny_true = []\nfor file in test_files:\n    if \"anomaly\" in file:\n        y_true.append(1)\n    else:\n        y_true.append(0)\ny_true = np.array(y_true)\n\nroc_auc_score(y_true, knn.decision_function(feat_vect_test_pca))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Final test submission\nfinal_test_files = [file.split(\"/\")[-1] for file in final_test_files]\nX_test_final = norm_set(X_test_final)\nfeat_vect_test_final = embedding_network.predict(X_test_final)\nanomaly_scores = knn.decision_function(pca.transform(sc.transform(feat_vect_test_final)))\nsubmission_df = pd.DataFrame(data={'file_name': final_test_files,\n                                   'anomaly_score': anomaly_scores.squeeze()})\nsubmission_df_sorted = submission_df.sort_values(by=\"file_name\", ignore_index=True)\nsubmission_df_sorted.to_csv(\"submission.csv\", index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df_sorted","metadata":{},"execution_count":null,"outputs":[]}]}