{"metadata":{"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"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":11000,"databundleVersionId":875412,"sourceType":"competition"}],"dockerImageVersionId":25160,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### Overview\n\nThis kernel implements a CNN with a squeeze-and-excitation mechanism. The rationale is using sigmoid gates to reduce the effects of noise  from the original inputs. The gate parameters are learned to differentiate useful filters and less useful ones, and hopefully less useful ones are more suppressed with lower sigmoid coefficients.\n\nThe original paper for the squeeze-and-excitation mechanism:\n\nhttps://arxiv.org/abs/1709.01507\n\nSome codes were borrowed from these kernels (Great thanks!):\n\nhttps://www.kaggle.com/gpreda/lanl-earthquake-eda-and-prediction\n\nhttps://www.kaggle.com/friedchips/how-to-reduce-the-training-data-to-400mb\n\nhttps://www.kaggle.com/tandonarpit6/lanl-earthquake-prediction-fast-ai\n\nhttps://www.kaggle.com/rtatman/download-a-csv-file-from-a-kernel\n\n","metadata":{}},{"cell_type":"code","source":"## Read in training data (from https://www.kaggle.com/gpreda/lanl-earthquake-eda-and-prediction)\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm_notebook\n\n%time train = pd.read_csv(\"../input/train.csv\", dtype={\"acoustic_data\": np.int16, \"time_to_failure\": np.float32})","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample Data Generation\n\nTo evaluate the squeeze-and-excitation mechanism with a subset of data, in total 2,550 samples with 15,0000 consecutive points were generated from 17 pre-failure periods. The model was trained with 5-fold cross validation and test data was predicted with the average outputs from the five models of each cv.","metadata":{}},{"cell_type":"code","source":"# Split the whole training segments into five cross-folds.\n# There are 17 pre-failure periods in the training data, so samples for five cross folds were randomly selected \n# from pre-failure periods with index 0 - 3, 4 - 6, 7 - 9, 10 - 12, and 13 - 16.\n\nsample_num = 150\ninput_len = 150000\nnp.random.seed(7898)\n\n# From https://www.kaggle.com/friedchips/how-to-reduce-the-training-data-to-400mb\nttf = train[\"time_to_failure\"].values\nindex_start = np.nonzero(np.diff(ttf) > 0)[0] + 1\nindex_start = np.insert(index_start, 0, 0) # insert 1st period start manually\nchunk_length = np.diff(np.append(index_start, train.shape[0]))\n\nX_series = []\ny_series = []\n# Split 17 pre-failure periods into five cross-folds, where the last four folds starts at index of 4, 7, 10 and 13th pre-failure period\ncv_assign = [4, 7, 10, 13]\nX, y = None, None\nfor i in range(len(index_start)):\n    if i in cv_assign:\n        X_series.append(X)\n        y_series.append(y)\n        X, y = None, None\n    index_set = np.random.randint(low=index_start[i], high=index_start[i] + chunk_length[i] - input_len, size=sample_num)\n    ac_data = np.zeros((sample_num,input_len, 1), dtype=np.int16)\n    ac_label = np.zeros((sample_num,), dtype=np.float32)\n    for j in range(sample_num):\n        ac_data[j, :, 0] = train[\"acoustic_data\"].values[index_set[j]:index_set[j] + input_len]\n        ac_label[j] = train[\"time_to_failure\"].values[index_set[j] + input_len]\n    if X is None: X, y = ac_data, ac_label\n    else:\n        X, y = np.concatenate((X, ac_data), axis=0), np.concatenate((y, ac_label), axis=0)\nX_series.append(X)\ny_series.append(y)\ndel train # Just to save space\n\n# From https://www.kaggle.com/tandonarpit6/lanl-earthquake-prediction-fast-ai\nsubmission = pd.read_csv(\"../input/sample_submission.csv\", index_col=\"seg_id\")\nX_test = np.zeros((2624, 150000, 1), dtype=np.int32)\ni = 0\nfor seg_id in tqdm_notebook(submission.index):\n    seg = pd.read_csv(\"../input/test/{}.csv\".format(seg_id))\n    X_test[i, :, :] = seg.values\n    i += 1","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model Architecture\n\nThe input into the model is a segment with length 15,0000. The model first reduces all signals by 100 fold, and followed by two convolutional layers. Squeeze-and-excitation mechanism and a max-pooling layer are filled in between two conv layers. The final output is calculated with an intermediate dense layer from the second conv layer.","metadata":{}},{"cell_type":"code","source":"#Generate CNN model with a squeeze-and-excitation mechanism\nfrom keras.models import Sequential, Model\nfrom keras.layers import Input, Conv1D, MaxPooling1D, Dense, Flatten, Multiply, Dropout, Lambda, GlobalAveragePooling1D\nfrom keras.optimizers import Adam\n\ndef squeeze_block(x, filter_num):\n    squeeze = GlobalAveragePooling1D()(x)\n    squeeze = Dense(int(filter_num / 2), activation=\"relu\", kernel_initializer=\"he_normal\")(squeeze)\n    squeeze = Dense(filter_num, activation=\"sigmoid\", kernel_initializer=\"he_normal\", name=\"squeeze_coef\")(squeeze)\n    c = Multiply()([x, squeeze])\n    return c\n\ndef model_gen(InputDim):\n    print(\"Building model ...\")\n    inputs = Input((InputDim, 1))\n    c = Lambda(lambda x: x / 100.0)(inputs)\n    c = Conv1D(filters=16, kernel_size=(10,), strides=1, padding=\"valid\", \n               kernel_initializer=\"he_normal\", activation=\"relu\", name=\"conv_before_squeeze\")(c)\n    c = squeeze_block(c, 16)\n    c = MaxPooling1D(pool_size=10, strides=10, padding=\"valid\")(c)\n    c = Conv1D(filters=16, kernel_size=(10,), strides=5, padding=\"valid\", \n               kernel_initializer=\"he_normal\", activation='relu', name=\"conv_after_squeeze\")(c)\n    f = Flatten()(c)\n    d = Dense(64, activation=\"relu\", kernel_initializer=\"he_normal\")(f)\n    d = Dropout(rate=0.5)(d)\n    outputs = Dense(1, activation=\"elu\", kernel_initializer=\"he_normal\")(d)\n    model = Model(inputs=[inputs], outputs=[outputs])\n    optim = Adam(lr=0.005)\n    model.compile(loss=\"mean_absolute_error\", optimizer=optim)\n    return model\nmodel = model_gen(150000)\nmodel.summary()","metadata":{"_uuid":"84b9a9e92b164eceb89a232a7d689d39a59a47e2","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model Training\n\nA decaying learning rate scheduler was applied.\n\nFor each cross-fold, the model was trained for 30 epochs, and the model from the last epoch was used to predict the test data.","metadata":{}},{"cell_type":"code","source":"# Model training with 5-cv\n\nimport numpy as np\nfrom keras.callbacks import LearningRateScheduler\n\nbatch_size = 50\nepochs_number = 30\n\ndef step_decay(epoch):\n    initial_lrate = 0.005\n    drop = 0.8\n    epochs_drop = 5.0\n    lrate = initial_lrate * np.power(drop, np.floor((1 + epoch) / epochs_drop))\n    return lrate\nlrate = LearningRateScheduler(step_decay)\n\nmodel_loss = np.zeros((5,), dtype=np.float32)\ny_pred_series = []\ny_test = np.zeros((2624, 1), dtype=np.float32)\n\nfor i in range(len(X_series)):\n    model = model_gen(input_len)\n    X_fit, y_fit = None, None\n    for j in range(len(X_series)):\n        if j != i:\n            if X_fit is None: X_fit, y_fit = X_series[j], y_series[j]\n            else:\n                X_fit, y_fit = np.concatenate((X_fit, X_series[j]), axis=0), np.concatenate((y_fit, y_series[j]), axis=0)\n    X_val, y_val = X_series[i], y_series[i]\n    model.fit(X_fit, y_fit, batch_size=batch_size, epochs=epochs_number, \n              validation_data=(X_val, y_val), shuffle=True, callbacks=[lrate])\n    model_loss[i] = model.evaluate(X_val, y_val)\n    y_pred_series.append(model.predict(X_val))\n    y_test += model.predict(X_test)\nmodel.save_weights(filepath=\"model_param.hdf5\")\n    \nmodel_loss_mean = np.mean(model_loss)\nmodel_loss_std = np.std(model_loss)\nprint(\"CV Score for the model is {:.3f} +/- {:.3f}\".format(model_loss_mean, model_loss_std))\nsubmission.time_to_failure = y_test / 5\nsubmission.to_csv(\"submission.csv\", index=True)","metadata":{"_uuid":"ecbab3c71227eb649d17ec93dac57768f282c73d","collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model Evaluation\n\nModel from each cross-fold was used to predict the validation fold to have a comparison with ground truth. The model tends to overestimate small failure time and underestimate large failure time.","metadata":{}},{"cell_type":"code","source":"print(\"CV Score for the model is {:.3f} +/- {:.3f}\".format(model_loss_mean, model_loss_std))\nfig, ax = plt.subplots(5, 1, figsize=(15, 10))\nfor i in range(len(y_pred_series)):\n    index = np.argsort(y_series[i])\n    ax[i].plot(np.arange(y_series[i].shape[0]), y_series[i][index], \"b-\")\n    ax[i].plot(np.arange(y_pred_series[i].shape[0]), y_pred_series[i][index], \"r.\")\n    ax[i].set_ylabel(\"Time in CV {}\".format(i + 1))\n    if i == 0:\n        ax[i].set_title(\"Truth (Blue) versus Prediction (Red)\")\n    ax[-1].set_xlabel(\"Samples ranked with increasing failure time\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To evaluate the efficiency of the squeeze-and-excitation mechanism, one sample with the largest failure time from the last cross-fold was picked to generate the outputs of convolutional layers. Note that each conv layer has 16 filters, so the outputs of first conv layer were shown in left four columns (start from second row), and ones for second conv layer were in right four columns.\n\nIt seems all squeeze coefficients for 16 filters in first conv layer were around 0.50. While it could mean the squeeze-and-excitation didn't do the job as expected, interestingly, the CNN model without it performs consistently WORSE than with it (didn't show in current kernel).","metadata":{}},{"cell_type":"code","source":"from keras import backend as K\n\nI = y_val.argsort()[-1] # Pick some sample with the largest failure time\nx = X_val[I, :, :]\ny = y_val[I]\nx = np.expand_dims(x, axis=0)\ny_pred = model.predict(x)[0][0]\n\nget_pre_squeeze_conv_layer_output = K.function([model.layers[0].input],\n                                  [model.get_layer(\"conv_before_squeeze\").output])\nget_post_squeeze_conv_layer_output = K.function([model.layers[0].input],\n                                      [model.get_layer(\"conv_before_squeeze\").output])\nget_squeeze_coefficients = K.function([model.layers[0].input],\n                                      [model.get_layer(\"squeeze_coef\").output])\n\nconv_output = np.array(get_pre_squeeze_conv_layer_output([x])[0])\nconv2_output = np.array(get_post_squeeze_conv_layer_output([x])[0])\nsqueeze_coefs = np.array(get_squeeze_coefficients([x])[0])\n\nr, c = 5, 8\nfig, ax = plt.subplots(r, c, figsize=(30, 10))\nfor i in range(r):\n    for j in range(c):\n        ax[i][j].set_xticks([])\n        if i == 0:\n            if j == 0:\n                ax[i][j].plot(np.arange(x.shape[1]), x[0, :, 0], \"b-\")\n                ax[i][j].set_title(\"True/Pred = {:.2f}/{:.2f}\".format(y, y_pred))\n            else:\n                ax[i][j].set_visible(False)\n        else:\n            if j <= 3:\n                k = c // 2 * (i - 1) + j\n                ax[i][j].plot(np.arange(conv_output.shape[1]), conv_output[0, :, k], \"b-\")\n                ax[i][j].set_title(\"Squeeze Coef: {:.2f}\".format(squeeze_coefs[0, k]))\n            else:\n                k = c // 2 * (i - 1) + j - 4\n                ax[i][j].plot(np.arange(conv2_output.shape[1]), conv2_output[0, :, k], \"b-\")\nax[2][0].set_ylabel(\"First Conv Outputs Before Squeezing\")\nax[2][4].set_ylabel(\"Second Conv Outputs After Squeezing\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One can download the submission.csv from the following codes. The PL score with predictions from this kernel was 1.922. I expect some better scores when one trains the model with the whole set of training data.","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/rtatman/download-a-csv-file-from-a-kernel\nfrom IPython.display import HTML\nimport base64\n\n# function that takes in a dataframe and creates a text link to  \n# download it (will only work for files < 2MB or so)\ndef create_download_link(df, title = \"Download CSV file\", filename = \"submission.csv\"):  \n    csv = df.to_csv()\n    b64 = base64.b64encode(csv.encode())\n    payload = b64.decode()\n    html = '<a download=\"{filename}\" href=\"data:text/csv;base64,{payload}\" target=\"_blank\">{title}</a>'\n    html = html.format(payload=payload,title=title,filename=filename)\n    return HTML(html)\n\n# create a link to download the dataframe\ncreate_download_link(submission)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}