{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport time\n\nfrom tqdm import tqdm\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Analysis"},{"metadata":{},"cell_type":"markdown","source":"Data Exploration\n\nWe are going to load the entire data in memory even though the dataset is fairly large. This is being done so that we can run some analysis on the original dataset."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain = pd.read_csv(\"../input/train.csv\", dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train.shape)\nprint(train.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,1, figsize=(20,10))\nax[0].plot(train['acoustic_data'].values[::100], color='g')\nax[0].set_title(\"Acoustic data for 1% sample data\")\nax[0].set_xlabel(\"Index\")\nax[0].set_ylabel(\"Acoustic Data Signal\");\nax[1].plot(train['time_to_failure'].values[::100], color='b')\nax[1].set_title(\"Time to Failure for 1% sample data\")\nax[1].set_xlabel(\"Index\")\nax[1].set_ylabel(\"Time to Failure in ms\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can further look into some more statistical features of the reduced dataset (2% sample of the total data) to understand the nature of the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"train.iloc[::50].describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Exploratory Visualization"},{"metadata":{},"cell_type":"markdown","source":"Looking at the above Time To Failure plot, we see that there are 16 earthquakes that occur (Time To Failue = 0). We can explore those specific points of earthquakes further in the next few plots, by investigating each of the quake points further by taking an approximate sample of 30,000,000 consecutive points (in steps of 50, so 2% of the 30,000,000 points) during which the quake occurs."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plotAroundPoints(start, end, ith):\n    fig, ax1 = plt.subplots(figsize=(8, 4))\n    plt.title(\"Trends of acoustic_data and time_to_failure around the {} earthquake\".format(ith))\n    plt.plot(train['acoustic_data'].values[start:end:50], color='b')\n    ax1.set_ylabel('acoustic_data', color='b')\n    plt.legend(['acoustic_data'])\n    ax2 = ax1.twinx()\n    plt.plot(train['time_to_failure'].values[start:end:50], color='g')\n    ax2.set_ylabel('time_to_failure', color='g')\n    plt.legend(['time_to_failure'], loc=(0.75, 0.1))\n\nplotAroundPoints(0, 30000000, \"first\")\nplotAroundPoints(30000000, 60000000, \"second\")\nplotAroundPoints(90000000, 120000000, \"third\")\nplotAroundPoints(125000000, 155000000, \"fourth\")\nplotAroundPoints(170000000, 200000000, \"fifth\")\nplotAroundPoints(200000000, 230000000, \"sixth\")\nplotAroundPoints(225000000, 255000000, \"seventh\")\nplotAroundPoints(285000000, 315000000, \"eigth\")\nplotAroundPoints(325000000, 355000000, \"ninth\")\nplotAroundPoints(360000000, 390000000, \"tenth\")\nplotAroundPoints(405000000, 455000000, \"eleventh\")\nplotAroundPoints(440000000, 470000000, \"twelvth\")\nplotAroundPoints(480000000, 510000000, \"thirteenth\")\nplotAroundPoints(510000000, 540000000, \"fourteenth\")\nplotAroundPoints(560000000, 590000000, \"fifteenth\")\nplotAroundPoints(605000000, 635000000, \"sixteenth\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observing the above plots, we can thus conclude that there always is a big spike in the acoustic signal value just before the earthquake. This could be useful information and can be used as a random visual verification metric (as a secondary form of verification, in addition to our metric scores) once our model has been implemented."},{"metadata":{},"cell_type":"markdown","source":"We also need to look at the provided test data folder, which contain the `acoustic_data` column against which we need to predict the `time_to_failure` values."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_dir = \"../input/test\"\ntest_files = os.listdir(test_dir)\nprint(test_files[0:5])\nprint(\"Number of test files: {}\".format(len(test_files)))\ntest_file_0 = pd.read_csv('../input/test/' + test_files[0])\nprint(\"Dimensions of the first test file: {}\".format(test_file_0.shape))\ntest_file_0.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we also need to check the submission file which needs to be populated with our predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(\"../input/sample_submission.csv\", index_col='seg_id', dtype={\"time_to_failure\": np.float32})\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(submission)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This verifies that there are `2624` rows in the `sample_submission.csv` file, which corresponds to the number of the test files in the test folder. We need to provide a value for the `time_to_failure` column in the `sample_submission.csv` file, each row of the file coinciding with a test file containing `150,000` values for the `acoustic_data` column."},{"metadata":{},"cell_type":"markdown","source":"Algorithms and Techniques"},{"metadata":{},"cell_type":"markdown","source":"Looking at the nature of the data, the following algorithms could be tried to determinepatterns in the test data (data in the test folder) that match the data in thetraining set (some or all may be implemented):\n* Polynomial Regression  \n  Since the nature of the data is non-linear, Polynomial Regression model could be used to check how the model performs against the given data. L2 regularisation can be set with a small lambda value to not overly punish the model.\n* Neural Network\n  I am also going to try to use a Neural Network model as the data is non-linear."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(train['acoustic_data'].values[::100], train['time_to_failure'].values[::100], s=10)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Benchmark"},{"metadata":{},"cell_type":"markdown","source":"Let us try to use a Linear Regression model to fit the training data and compute a Mean Absolute Error (MAE) score, so that we can get a benchmark against a simple model. Since the data set is massive, and is contiguous in nature, we will take a sample of the entire data i.e. 2% and try to derive a linear regression model on it."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import make_scorer\nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import median_absolute_error\n\nX_train, X_test, y_train, y_test = train_test_split(train['acoustic_data'].values[::25].reshape(-1, 1), train['time_to_failure'].values[::25], test_size=0.2)\n\nquake_linear_model = LinearRegression()\nquake_linear_model.fit(X_train, y_train)\n\ny_train_pred = quake_linear_model.predict(X_train)\ny_test_pred = quake_linear_model.predict(X_test)\n\nr2_train_score = r2_score(y_train, y_train_pred)\nmae_train_score = median_absolute_error(y_train, y_train_pred)\nr2_test_score = r2_score(y_test, y_test_pred)\nmae_test_score = median_absolute_error(y_test, y_test_pred)\n\nprint(\"R2 score for training data: {} and for the test data: {}\".format(r2_train_score, r2_test_score))\nprint(\"Mean Absolute Error score for training data: {} and for the test data: {}\".format(mae_train_score, mae_test_score))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**III Methodology**"},{"metadata":{},"cell_type":"markdown","source":"**Data Preprocessing**"},{"metadata":{},"cell_type":"markdown","source":"Firstly, we need to add some more features to this data set, which contains only 1 column for the input signal.\n\nThis can be done by generating statistical information over chunks of the original dataset of length 150,000. This is the length of the test files that we need to make the predictions for. (Ref: Basic Feature Benchmark by Inversion: https://www.kaggle.com/inversion/basic-feature-benchmark)\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"chunk_size = 150000\n\nchunks = int(np.floor(train.shape[0]/chunk_size))\n\nX_data = pd.DataFrame(index=range(chunks), dtype=np.float32, columns=['min','max','std', 'avg', 'sum', 'median', 'mean_diff', \n                                                                       'q05', 'q25', 'q75', 'q95'])\ny_data = pd.DataFrame(index=range(chunks), dtype=np.float32, columns=['ttf'])\n\ndef create_features(data_chunk, X_df, chunk_no, col_name='acoustic_data'):\n    x = data_chunk[col_name]\n    X_df.loc[chunk_no, 'min'] = x.min()\n    X_df.loc[chunk_no, 'max'] = x.max()\n    X_df.loc[chunk_no, 'std'] = x.std()\n    X_df.loc[chunk_no, 'avg'] = x.mean()\n    X_df.loc[chunk_no, 'sum'] = x.sum()\n    X_df.loc[chunk_no, 'median'] = x.median()\n    X_df.loc[chunk_no, 'mean_diff'] = np.mean(np.diff(x))\n    X_df.loc[chunk_no, 'q05'] = np.quantile(x, 0.05)\n    X_df.loc[chunk_no, 'q25'] = np.quantile(x, 0.25)\n    X_df.loc[chunk_no, 'q75'] = np.quantile(x, 0.75)\n    X_df.loc[chunk_no, 'q95'] = np.quantile(x, 0.95)\n    return X_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for chunk_no in tqdm(range(chunks)):\n    data_chunk = train.iloc[chunk_no*chunk_size:chunk_no*chunk_size+chunk_size]\n    X_data = create_features(data_chunk, X_data, chunk_no)\n    y = data_chunk['time_to_failure'].values[-1]\n    y_data.loc[chunk_no, 'ttf'] = y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(X_data.shape)\nprint(y_data.shape)\nprint(X_data.shape[1])\nX_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n**Implementation**"},{"metadata":{},"cell_type":"markdown","source":"Now that we have created a new dataset with additional features, the next step is to divide this new dataset into training and test sets."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(X_data.values, y_data.values, test_size=0.2)\n# X_test\n# X_data.values\nX_train.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we shall try and build a RNN (Recurrent Neural Network) model with Keras, as the training data is sequential in nature. "},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.models import Sequential\nfrom keras.layers import Dense, Dropout, CuDNNGRU, CuDNNLSTM, Flatten\nfrom keras.optimizers import adam\nfrom keras.callbacks import ModelCheckpoint\n\nmodel = Sequential()\n# model.add(CuDNNLSTM(64, kernel_initializer=\"RandomUniform\", input_shape= (X_train.shape[1], 1)))\nmodel.add(CuDNNGRU(64, kernel_initializer=\"RandomUniform\", input_shape= (X_train.shape[1], 1)))\nmodel.add(Dropout(0.2))\nmodel.add(Dense(64, activation=\"relu\"))\nmodel.add(Dropout(0.2))\nmodel.add(Dense(64, activation=\"relu\"))\nmodel.add(Dropout(0.2))\nmodel.add(Dense(64, activation=\"relu\"))\nmodel.add(Dropout(0.2))\n# model.add(Flatten())\nmodel.add(Dense(1))\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.callbacks import ModelCheckpoint\n\n# Reshaping for fit\n# X_train_array = np.reshape(X_train, (X_train.shape[0], 1, X_train.shape[1]))\nX_train_array = np.reshape(X_train, (X_train.shape[0], X_train.shape[1], 1))\ny_train_array = np.reshape(y_train, (y_train.shape[0], y_train.shape[1], 1))\n\n# model.compile(loss=\"mean_squared_error\", optimizer=\"rmsprop\", metrics=[\"mse\"])\nmodel.compile(loss=\"mean_absolute_error\", optimizer=\"adam\", metrics=[\"mae\", \"mse\"])\n\ncheckpointer = ModelCheckpoint(\"model.weights.hdf5\", save_best_only=True, verbose=1)\n\nbuild = model.fit(X_train_array, y_train, epochs=200, batch_size=30, validation_split = 0.20, callbacks=[checkpointer], verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(build.history.keys())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# summarize history for loss\nplt.plot(build.history['loss'])\nplt.plot(build.history['val_loss'])\nplt.title('model loss')\nplt.ylabel('loss')\nplt.xlabel('epoch')\nplt.legend(['train', 'validation'], loc='upper right')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_array = np.reshape(X_train, (X_train.shape[0], X_train.shape[1], 1))\npredictions = model.predict(X_train_array)\n\nr2_pred_score = r2_score(y_train, predictions)\nmae_pred_score = median_absolute_error(y_train, predictions)\n\nprint(\"R2 score for training data: {}\".format(r2_pred_score))\nprint(\"Mean Absolute Error score for training data: {}\".format(mae_pred_score))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test_array = np.reshape(X_test, (X_test.shape[0], X_test.shape[1], 1))\npredictions_test = model.predict(X_test_array)\n\nr2_pred_test_score = r2_score(y_test, predictions_test)\nmae_pred_test_score = median_absolute_error(y_test, predictions_test)\n\nprint(\"R2 score for test data: {}\".format(r2_pred_test_score))\nprint(\"Mean Absolute Error score for test data: {}\".format(mae_pred_test_score))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Refinement**\n\nAs mentioned above, RNN was selected as the model (RNN layers) given the nature of the data. \nThe various layers available in RNN (Keras) are: \n* RNN\n* GRU\n* LSTM\n* CuDNNGRU\n* CuDNNLSTM\netc.\n\nThe first model was implemented using the following parameter values:\n* epochs = 30\n* Number of units in the Dense layers = 32\n* Dropout probability of 0.1\n* CuDNNLSTM as in the input layer for the RNN implementation in Keras\n\nThese parameters were then continually refined to obtain better results: \n* The results obtained using these parameters were good with mean validation loss (mean squared error) around 7.85. The number of epochs were then incrementally increased, which kept reducing the validation loss and hence the number of epochs were then finally set to 200 as increasing the epochs further did not result in further improvement over the validation loss, causing overfitting.\n* The number of units in the Dense layers of the network were increased to 64, which also improved the model results with a lower validation loss.\n* Increase in the Dropout value from 0.1 to 0.2 also yielded better results.\n* And finally, switching the *CuDNNLSTM* input layer with *CuDNNGRU* implementation further improved the result.\n"},{"metadata":{},"cell_type":"markdown","source":"**Results**\n\n**Model Evaluation and Validation**\n\nOnce we built our final model, we plotted a learning curve for the validation loss over the training and validation sets (refer to Figure).\n\nThis plot shows that the training loss reduces continually with the validation loss closely following the training loss. This trend in the loss plot confirms a good fit of our model.\n\n**Justification**\n\nR2 score for training data: 2.9657491629198063e-07 and for the test data: 4.92912790694966e-07\nMean Absolute Error score for training data: 2.804914134725992 and for the test data: 2.8038857678803177\n\nLets tabulate our final results in terms of the metrics i.e. the R2 score and the Mean Absolute Error against the benchmark of our sample data.\n\n____________________________________________________________________________\n|    |Training Set Score Benchmark/Final  | Test Set Score Benchmark/Final |\n|____|____________________________________|________________________________|\n| R2 | 2.9657e-07/0.4471                  | 4.9291e-07/0.4147              |\n|____|____________________________________|________________________________|\n| MAE| 2.8049/1.6096                      | 2.8038/1.6573                  |\n____________________________________________________________________________\n\nLooking at the comparison of the benchmark vs the implemented model metrics, there is a massive improvement for the R2 score of our final model over the benchmark model, as well as a big reduction in the mean absolute error as well. This comparison gives us confidence in the implemented model to predict results for the given test set (different from the test set used above).\n"},{"metadata":{},"cell_type":"markdown","source":"**V. Conclusion**\n\n**Free Form Visualization**"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_sub = pd.DataFrame(columns=X_data.columns, dtype=np.float32)\n\nfor i, seg_id in enumerate(tqdm(submission.index)):\n    seg = pd.read_csv('../input/test/' + str(seg_id) + '.csv')\n    X_seg = create_features(seg, X_sub, i)\n    # print(X_seg)\n    # X_seg_array = np.reshape(X_seg.values, (X_seg.shape[0], X_seg.shape[1], 1))\n    # pred_seg = model.predict(X_seg_array)\n    # print(pred_seg)\n    # submission.time_to_failure[i] = pred_seg","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(X_sub.shape)\nX_sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_seg_array = np.reshape(X_seg.values, (X_seg.shape[0], X_seg.shape[1], 1))\npred_final = model.predict(X_seg_array)\nsubmission['time_to_failure'] = pred_final\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After running the predictions on the test files using the model implemented, we generate an output submission.csv by outputting the time_to_failure values against each test file. \n\nFurthermore, on doing a simple search for all the time_to_failure values that are <1.0, we get a list of the following test files and their predicted time_to_failure values:\n\n|seg_26a2a0|0.88940096|\n|__________|__________|\n|seg_724df9|0.61797214|\n|seg_7a9f2b|0.92327976|\n|seg_7fa6ec|0.96871   |\n|seg_aa98cc|0.8351287 |\n|seg_c80857|0.6413747 |\n|seg_e3d751|0.8323933 |\n\nWe next plot the acoustic_data values for each of these files to visualise the trends in the input signal values."},{"metadata":{"trusted":true},"cell_type":"code","source":"possible_eq = submission.loc[submission[\"time_to_failure\"] < 1.0]\nprint(possible_eq)\nprint(type(possible_eq))\nprint(possible_eq.columns)\n# segments = [\"seg_26a2a0\", \"seg_724df9\", \"seg_7a9f2b\", \"seg_7fa6ec\", \"seg_aa98cc\", \"seg_b35174\", \"seg_c80857\", \"seg_e3d751\"]\n\nfor seg_id in possible_eq.index:\n    print(seg_id)\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    fig, ax1 = plt.subplots(figsize=(8, 4))\n    plt.title(\"Trends of acoustic_data for test file {}\".format(seg_id))\n    plt.plot(seg['acoustic_data'].values, color='b')\n    ax1.set_ylabel('acoustic_data', color='b')\n    plt.legend(['acoustic_data'])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The above plots of the acoustic_data values for all the test files (each containing 150,000 samples) that were predicted with a time_to_failure of < 1.0 show a massive peak in the acoustic_data values, in the range of 3000-6000. This observation is in line with the trend also observed in the input data sample where each of the 16 earthquakes (time_to_failure=0) were preceded by a massive spike just before the occurrence of the earthquake (refer to Figure). \n"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}