{"cells":[{"metadata":{"_uuid":"9d05929a6055563ff49aa0be47c96fb60c326698"},"cell_type":"markdown","source":"For me, this VSB power line competition was a good chance to learn how to use LSTM or RNN in general (I expect GRU should not be much different to apply with Keras..). I need a place to write things down so I remember another day and not just today. I wrote myself a [blog post](https://swenotes.wordpress.com/2019/02/22/learning-to-lstm/) to remind myself. This kernel is an attempt to put some working code somewhere.\n\nIf I got any part wrong about here, or missing something, do let me know :).\n\nI started with the [kernel](https://www.kaggle.com/braquino/5-fold-lstm-attention-fully-commented-0-694) by Bruno Marek. Then played with the data and classifiers myself, built a separate [preprocessing kernel](https://www.kaggle.com/donkeys/preprocessing-with-python-multiprocessing) as well. This kernel uses data produced by that preprocessing kernel.\n\nThere are other public kernels in the competition with better scores but I wanted to keep this simple to help myself more clearly understand the core concepts.\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport pyarrow.parquet as pq # Used to read the data\nimport os \nimport numpy as np\nfrom keras.layers import *\nfrom keras.models import Model\nfrom sklearn.model_selection import train_test_split \nfrom keras import backend as K \nfrom keras import optimizers\nimport tensorflow as tf\nfrom sklearn.model_selection import GridSearchCV, StratifiedKFold\nfrom keras.callbacks import *\n%matplotlib inline\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# select how many folds will be created\nN_SPLITS = 5\n# it is just a constant with the measurements data size\nsample_size = 800000","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ad473b0c30029b79cf45700fdc44ecf994186c0c"},"cell_type":"markdown","source":"Matthews correlation coefficient is simply the measure given in this Kaggle competition as a way to measure the score. It [seems](https://en.wikipedia.org/wiki/Matthews_correlation_coefficient) that a value of 1 would mean perfect prediction, and 0 equal to random values. So maybe 0.6-0.7 that many kernels get is not all that bad? "},{"metadata":{"trusted":true,"_uuid":"d892404a7d36bb161574aac74d42afcefcc7d44d"},"cell_type":"code","source":"def matthews_correlation_coeff(y_true, y_pred):\n    '''Calculates the Matthews correlation coefficient measure for quality\n    of binary classification problems.\n    '''\n    y_pred = tf.convert_to_tensor(y_pred, np.float32)\n    y_true = tf.convert_to_tensor(y_true, np.float32)\n\n    y_pred_pos = K.round(K.clip(y_pred, 0, 1))\n    y_pred_neg = 1 - y_pred_pos\n\n    y_pos = K.round(K.clip(y_true, 0, 1))\n    y_neg = 1 - y_pos\n\n    tp = K.sum(y_pos * y_pred_pos)\n    tn = K.sum(y_neg * y_pred_neg)\n\n    fp = K.sum(y_neg * y_pred_pos)\n    fn = K.sum(y_pos * y_pred_neg)\n\n    numerator = (tp * tn - fp * fn)\n    denominator = K.sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn))\n\n    return numerator / (denominator + K.epsilon())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6ab983e9f5733a42ec4d5e4a8b1532b1453cc244"},"cell_type":"code","source":"# load the training set metadata, defines which signals are in which order in the data\ntrain_meta = pd.read_csv('../input/vsb-power-line-fault-detection/metadata_train.csv')\n# set index, it makes the data access much faster\ntrain_meta = train_meta.set_index(['id_measurement', 'phase'])\ntrain_meta.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"04c76b951dad765be6d88a96f15aefec0a03b057"},"cell_type":"code","source":"# load the test set metadata, defines which signals are in which order in the data\ntest_meta = pd.read_csv('../input/vsb-power-line-fault-detection/metadata_train.csv')\n# set index, it makes the data access much faster\ntest_meta = test_meta.set_index(['id_measurement', 'phase'])\ntest_meta.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e938b19c039d1b0a56d6a0f1df74031e0c5c64fc"},"cell_type":"code","source":"!ls ../input/preprocessing-with-python-multiprocessing","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"350be95fee64e5501ba8dd349ee10438b7b8e773"},"cell_type":"markdown","source":"The data files produced by the preprocessing kernel is shown above. I had to compress them using gzip to fit them into the kernel 5GB output size limit. Hence the decompression and the filename suffix here."},{"metadata":{"trusted":true,"_uuid":"fb2d20e2d696000313594830cee891007c943ce8"},"cell_type":"code","source":"df_test_pre = pd.read_csv(\"../input/preprocessing-with-python-multiprocessing/my_test_combined_scaled.csv.gz\", compression=\"gzip\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8280547b2c5c3adaf28d53fa874934862687b51f"},"cell_type":"code","source":"df_test_pre.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c68738513aebfcd9f363695b7228e960b54999b0"},"cell_type":"markdown","source":"The test dataframe loaded above has 22 features calculated for each of the 3 signals per measurement id. So 66 columns. It becomes 67 when loaded, because dumping the values to disk with pandas.to_csv seems to have generated one extra column (maybe the index?).  Number of actual rows should match the number of measurement id's in the corresponding dataset:"},{"metadata":{"trusted":true,"_uuid":"34a062e44abec351a1b7351d31173df71587ab72"},"cell_type":"code","source":"1084640/160","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"90632d6de91dfd52be1fccca109ad2630ecc819e"},"cell_type":"markdown","source":"The training data-set should look about the same but with only the 2904 measurements, so fewer rows:"},{"metadata":{"trusted":true,"_uuid":"651386a3a97046f6a440eb32f2fd0629eb88d9b4"},"cell_type":"code","source":"df_train_pre = pd.read_csv(\"../input/preprocessing-with-python-multiprocessing/my_train_combined_scaled.csv.gz\", compression=\"gzip\")\ndf_train_pre.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5be9275fa98790adbcc918bd0da8a4703f58acce"},"cell_type":"markdown","source":"To drop the excess column:"},{"metadata":{"trusted":true,"_uuid":"6b3132ab94f4363bd5c67516f3115116f5e59f44"},"cell_type":"code","source":"df_train_pre.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"aeda5c10faf750695b6a2cc1bb8fc7bdd8cafbcf"},"cell_type":"code","source":"df_train_pre.drop(\"Unnamed: 0\", axis=1, inplace=True)\ndf_test_pre.drop(\"Unnamed: 0\", axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"268d81f222ed136ec2ee9821446fd17af3ba850e"},"cell_type":"markdown","source":"The preprocessed data has 160 timesteps. Number of rows should match the number of measurements times the number of timesteps:"},{"metadata":{"trusted":true,"_uuid":"51da21a0c7c4c2e5b38ed59e6d7a8622ea1f77fb"},"cell_type":"code","source":"#number of \"observations\" in test dataset\ndf_test_pre.shape[0]/160","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"70a5c35c03a05a031f650efbe792bcfbd5da2634"},"cell_type":"code","source":"df_train_pre.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"22e8d03b21ff871308d1f9a851695fa2fd043900"},"cell_type":"code","source":"#number of \"observations\" in training dataset\n464640/160","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"23447fc71d67c2a94249974b023b0ad1e1e2fa42"},"cell_type":"code","source":"train_meta.index.get_level_values('id_measurement').unique()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bd76fad6f620c23d2057324bff512420fceddc93"},"cell_type":"markdown","source":"LSTM is about timesteps, and in this case the 800k measurements per signal were summarized to 160 timesteps per signal in the pre-processing. So things like average values of 5000 measurementes (800k/5000=160). These are now in rows 0-159 for the first measurements id, where the columns 0-21 are for the first signal, columns 22-43 for second signal, and 44-65 for the third signal.\n\nThis continues for the following measurements with the 3 signals per measurement id in the columns. So the signals for the second measurement id are in rows 160-319.\n\nA look at first signal for the first measurement id:"},{"metadata":{"trusted":true,"_uuid":"2b9f2cc68d598fa14e35fd037ba7dae565b492d3"},"cell_type":"code","source":"pd.set_option('display.max_rows', 5)\ndf_train_pre.iloc[0:160,:22]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a6347010934564e651745e9cb1ddb8a69b04d440"},"cell_type":"markdown","source":"Second signal:"},{"metadata":{"trusted":true,"_uuid":"2139528960831684c24e9dadbb33a5ca701b0077"},"cell_type":"code","source":"df_train_pre.iloc[0:160,22:44]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9a76cdf4c2f1ddf74c358dec0d016d26916efdcf"},"cell_type":"markdown","source":"Third signal:"},{"metadata":{"trusted":true,"_uuid":"72bb1fbb510737a78fc5ec71347e36e5512ebe85"},"cell_type":"code","source":"df_train_pre.iloc[0:160,44:66]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a7266b17a55655d08d2e765bf431399aba9280a1"},"cell_type":"markdown","source":"And the 3 signals for the second measurement id:"},{"metadata":{"trusted":true,"_uuid":"bbbe99fa67cda6178995a315a98baff6b343955b"},"cell_type":"code","source":"df_train_pre.iloc[160:320]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7e40047ccd43b9587ba282d97c17cad5ee4b01e1"},"cell_type":"code","source":"pd.reset_option('display.max_rows')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f61162d75cc9848929bc50108613f546e28afb45"},"cell_type":"code","source":"from sklearn.metrics import matthews_corrcoef\n\n# The output of this kernel must be binary (0 or 1), but the output of the NN Model is float (0 to 1).\n# So, find the best threshold to convert float to binary is crucial to the result\n# this piece of code is a function that evaluates all the possible thresholds from 0 to 1 by 0.01\ndef threshold_search(y_true, y_proba):\n    best_threshold = 0\n    best_score = 0\n    scores = []\n    for threshold in [i * 0.01 for i in range(100)]:\n        yp_np = np.array(y_proba)\n        yp_bool = yp_np >= threshold\n        score = matthews_corrcoef(y_true, yp_bool)\n        #score = K.eval(matthews_correlation(y_true.astype(np.float64), (y_proba > threshold).astype(np.float64)))\n        scores.append(score)\n        if score > best_score:\n            print(\"found better score:\"+str(score)+\", th=\"+str(threshold))\n            best_threshold = threshold\n            best_score = score\n    search_result = {'threshold': best_threshold, 'matthews_correlation': best_score}\n    scores_df = pd.DataFrame({\"score\": scores})\n    print(\"scores plot:\")\n    scores_df.plot()\n    plt.show()\n    return search_result","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7423e955bb77b7c693e693e829a752039edbcd2a"},"cell_type":"markdown","source":"Create the actual LSTM model. For more explanations, see my [blog post](https://swenotes.wordpress.com/2019/02/22/learning-to-lstm/)."},{"metadata":{"trusted":true,"_uuid":"baf198807424261f69cd4e1dacfcb08f2b076a2e"},"cell_type":"code","source":"def create_model(input_data):\n    input_shape = input_data.shape\n    inp = Input(shape=(input_shape[1], input_shape[2],), name=\"input_signal\")\n    x = Bidirectional(CuDNNLSTM(128, return_sequences=True, name=\"lstm1\"), name=\"bi1\")(inp)\n    x = Bidirectional(CuDNNLSTM(64, return_sequences=False, name=\"lstm2\"), name=\"bi2\")(x)\n    #other kernels have used also a custom Attention layer but I leave it out for simplicity here\n#    x = Attention(input_shape[1])(x)\n    x = Dense(128, activation=\"relu\", name=\"dense1\")(x)\n    x = Dense(64, activation=\"relu\", name=\"dense2\")(x)\n    x = Dense(1, activation='sigmoid', name=\"output\")(x)\n    model = Model(inputs=inp, outputs=x)\n    model.compile(loss='binary_crossentropy', optimizer='adam', metrics=[matthews_correlation_coeff])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"76b14da058fb7aef695e84385b38c5b3120337d5"},"cell_type":"markdown","source":"Since the dataset has been combined to have all 3 phase signals per measurement id on a single row (22\\*3=66 features/columns), as combined features, I need to combine the prediction targets for all 3 signals also into one. "},{"metadata":{"trusted":true,"_uuid":"239f7df87a4fc532f5bf85ed45108498e648e07a"},"cell_type":"code","source":"#if any of the 3 signals for a measurement id is labeled as faulty, this labels the whole set of 3 as faulty\ny = (train_meta.groupby(\"id_measurement\").sum()/3 > 0)[\"target\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a5914d97cd75f8a7dfd08a1a47ad6dca5880c550"},"cell_type":"code","source":"#to see the number of targets matches the number of rows in training dataset\ny.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"845a53f037e08b3cad8d61babc8c9f9bf9d15a00"},"cell_type":"markdown","source":"LSTM requires 3-dimensional input, so this reshapes the dataframe values from 2D dataframe to a 3D numpy matrix in the required format. 2904 observations, 160 timesteps, each with 66 features. \n\nFeatures for each timestep are on a single row in the dataframe (66 on a row) as is. The dataframe here being \"df_train_pre\". The dataframe has 160 rows per measurement id as shown above (df_train_pre.iloc[0:160] for measurement id 1 and df_train_pre.iloc[160:320] for measurement id 2, and so on). Each of these measurement id sets should be its own \"observation\" in the numpy matrix used as input for the LSTM. \n\nThe following reshape creates the required input format, setting the overall input shape as (2904, 160, 66). This is 2904 observations, 160 timesteps for each of those 2904 observations, and 66 features for each of those 160 timesteps. This is the 3D format format LSTM expects as input. A timestep has been formed by splitting the sequence of signal values over time to 160 separate values on after the other, and collecting the 66 features for that timeslot."},{"metadata":{"trusted":true,"_uuid":"47572a667c8206b1c86c634829dcfd545c44b704"},"cell_type":"code","source":"#if using all signal values separately, the number of rows would be 8712, or 2904*3.\n#X = df_train_pre.values.reshape(8712, 160, 22)\n#but with the current data format I show above, it is 2904 rows, or \"observations\"\nX = df_train_pre.values.reshape(2904, 160, 66)\nX.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d718714f990b294e56b2627118bcd648ce2ba0f1"},"cell_type":"markdown","source":"Now to do the same for the test-dataset, but remembering it has more rows, so the first dimension is higher. Maybe because the people at Kaggle want to make life difficult for the competitors and so the test set is much bigger :)."},{"metadata":{"trusted":true,"_uuid":"aaa44bd8ec393e5570f5967c1ecbfca776fb07d5"},"cell_type":"code","source":"X_test = df_test_pre.values.reshape(6779, 160, 66)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e8ba0fa7e0e7e8bbb734984ceeb1dcf975fdaef7"},"cell_type":"markdown","source":"Finally, train and run a simple LSTM classifier for all this:"},{"metadata":{"trusted":true,"_uuid":"77513d34d31037a14211c0828cbe04e4f97d74c2"},"cell_type":"code","source":"eval_preds = np.zeros(X.shape[0])\nlabel_predictions = []\n\nsplits = list(StratifiedKFold(n_splits=N_SPLITS, shuffle=True, random_state=123).split(X, y))\nfor idx, (train_idx, val_idx) in enumerate(splits):\n    K.clear_session()\n    print(\"Beginning fold {}\".format(idx+1))\n    train_X, train_y, val_X, val_y = X[train_idx], y[train_idx], X[val_idx], y[val_idx]\n\n    model = create_model(X)\n    #checkpoint to save model with best validation score. keras seems to add val_xxxxx as name for metric to use here\n    ckpt = ModelCheckpoint('weights.h5', save_best_only=True, save_weights_only=True, monitor='val_matthews_correlation_coeff', verbose=1, mode='max')\n    earlystopper = EarlyStopping(patience=25, verbose=1) \n    model.fit(train_X, train_y, batch_size=128, epochs=50, validation_data=[val_X, val_y], callbacks=[ckpt, earlystopper])\n    # loads the best weights saved by the checkpoint\n    model.load_weights('weights.h5')\n\n    print(\"finding threshold\")\n    predictions = model.predict(val_X, batch_size=512)\n    best_threshold = threshold_search(val_y, predictions)['threshold']\n    \n    print(\"predicting test set\")\n    pred = model.predict(X_test, batch_size=300, verbose=1)\n    pred_bool = pred > best_threshold\n    labels = pred_bool.astype(\"int32\")\n    label_predictions.append(labels)\n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cc42a67794763aeb66832eadb82352e800d2a8d1"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d7dca4cfe67e8eccdf536e3fbb26f9055a39d6e7"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4df7db6f6454cbc7a6740cc859e26753953ea04f"},"cell_type":"markdown","source":"Convert the above predictions into suitable submission format for the competition:"},{"metadata":{"trusted":true,"_uuid":"a086b1c9aa293e62ce75930e40e68e8068327958"},"cell_type":"code","source":"label_predictions = [pred.flatten() for pred in label_predictions]\n\nimport scipy\n\n# Ensemble with voting\nlabels = np.array(label_predictions)\n#convert list of predictions into set of columns\nlabels = np.transpose(labels, (1, 0))\n#take most common value (0 or 1) or each row\nlabels = scipy.stats.mode(labels, axis=-1)[0]\nlabels = np.squeeze(labels)\n\nsubmission = pd.read_csv('../input/vsb-power-line-fault-detection/sample_submission.csv')\nlabels3 = np.repeat(labels, 3)\nsubmission['target'] = labels3\nsubmission.to_csv('_voted_submission.csv', index=False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8a38931fc7f4d22d1266874590bd43cfd2a95d57"},"cell_type":"markdown","source":"Just a quick look at how many positive (faulty power line) predictions did we get:"},{"metadata":{"trusted":true,"_uuid":"4a4b5b127f85585f6d8db1b0dfd91503bb086ef0"},"cell_type":"code","source":"sum(labels3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"76f2dda07ed573fcfe3e59769e99b0376fb76595"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}