{"cells":[{"metadata":{"_uuid":"3b7efa15014e822ca40ee55d2e7e712a6ca63563"},"cell_type":"markdown","source":"# <a id='0'>Content</a>\n\n- <a href='#1'>Introduction</a>  \n- <a href='#2'>Prepare the data analysis</a>  \n- <a href='#3'>Data exploration</a>   \n- <a href='#4'>Feature engineering</a>\n- <a href='#5'>Model</a>\n- <a href='#6'>Submission</a>  \n- <a href='#7'>References</a>"},{"metadata":{"_uuid":"0bf2c2c00334bdc9f8f3d78d10d958bc521b4c7a"},"cell_type":"markdown","source":"# <a id='1'>Introduction</a>  \n\nThe data are from an experiment conducted on rock in a double direct shear geometry subjected to bi-axial loading, a classic laboratory earthquake model.\n\nTwo fault gouge layers are sheared simultaneously while subjected to a constant normal load and a prescribed shear velocity. The laboratory faults fail in repetitive cycles of stick and slip that is meant to mimic the cycle of loading and failure on tectonic faults. While the experiment is considerably simpler than a fault in Earth, it shares many physical characteristics. \n\nLos Alamos' initial work showed that the prediction of laboratory earthquakes from continuous seismic data is possible in the case of quasi-periodic laboratory seismic cycles.   \n\nIn this competition, the team has provided a much more challenging dataset with considerably more aperiodic earthquake failures.  \n\nObjective of the competition is to predict the failures for each training set."},{"metadata":{"_uuid":"3b198f4c2141e0e94e0853cea5798c0e83343716"},"cell_type":"markdown","source":"# <a id='1'>Prepare the data analysis</a>\n\n## Load packages\n\nHere we define the packages for data manipulation, feature engineering and model training."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import gc\nimport os\nimport tqdm\nimport logging\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport lightgbm as lgb\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import KFold,StratifiedKFold","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f7dc473a43cbd799a184ca71a27391bb6faa9984"},"cell_type":"markdown","source":"## Load the data\n\nLet's see first what files we have in input directory."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"45e5498b89764854036baa88eda6ead21d41f382"},"cell_type":"code","source":"IS_LOCAL = False\nif(IS_LOCAL):\n    PATH=\"../input/LANL/\"\nelse:\n    PATH=\"../input/\"\nos.listdir(PATH)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68ced52082fc77c768db9ff14fed7cbe43e097d1"},"cell_type":"markdown","source":"Let's load the train file."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"bb8916c762c26e98e31a721e4aae673aa12e6c5f"},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(PATH+'train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"Let's check the data imported."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"40c15befb5782fe1d671ae4ad0ca592aff316073"},"cell_type":"code","source":"print(\"Train: rows:{} cols:{}\".format(train_df.shape[0], train_df.shape[1]))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"0ee349f0d7f482a0cc4827a3b7515a915e4dd37f"},"cell_type":"code","source":"pd.options.display.precision = 15\ntrain_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"877ce7e2c62dc0c7b9a86c081754a7f3b19bb364"},"cell_type":"markdown","source":"# <a id='3'>Data exploration</a>  \n\nThe dimmension of the data is quite large, in excess of 600 millions rows of data.  \nThe two columns in the train dataset have the following meaning:   \n*  accoustic_data: is the accoustic signal measured in the laboratory experiment;  \n* time to failure: this gives the time until a failure will occurs.\n\nLet's plot 1% of the data. For this we will sample every 100 points of data.  "},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"206af56767581d8e98777da3b0670e11d4e444bd"},"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[::100]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[::100]\n\ndef plot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% sampled data\"):\n    fig, ax1 = plt.subplots(figsize=(12, 8))\n    plt.title(title)\n    plt.plot(train_ad_sample_df, color='r')\n    ax1.set_ylabel('acoustic data', color='r')\n    plt.legend(['acoustic data'], loc=(0.01, 0.95))\n    ax2 = ax1.twinx()\n    plt.plot(train_ttf_sample_df, color='b')\n    ax2.set_ylabel('time to failure', color='b')\n    plt.legend(['time to failure'], loc=(0.01, 0.9))\n    plt.grid(True)\n\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df)\ndel train_ad_sample_df\ndel train_ttf_sample_df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"36c363a5f2e56e54c8ffddd802aa86d91c61d97f"},"cell_type":"markdown","source":"The plot shows only 1% of the full data. \nThe acoustic data shows complex oscilations with variable amplitude. Just before each failure there is an increase in the amplitude of the acoustic data. We see that large amplitudes are also obtained at different moments in time (for example about the mid-time between two succesive failures).  \n\nLet's plot as well the first 1% of the data."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"845e5a241e545d73cfb1dbaea2cfa464faf8148a"},"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[:6291455]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[:6291455]\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% of data\")\ndel train_ad_sample_df\ndel train_ttf_sample_df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"78948b35358020156422a76c8d60dc132f17417b"},"cell_type":"markdown","source":"On this zoomed-in-time plot we can see that actually the large oscilation before the failure is not quite in the last moment. There are also trains of intense oscilations preceeding the large one and also some oscilations with smaller peaks after the large one. Then, after some minor oscilations, the failure occurs."},{"metadata":{"_uuid":"8dee473309c39d74020b37f2e4863094191c121f"},"cell_type":"markdown","source":"# <a id='4'>Features engineering</a>  \n\nThe test segments are 150,000 each.   \nWe split the train data in segments of the same dimmension with the test sets.\n\nWe will create additional aggregation features, calculated on the segments. The features we create are as following:\n* average value on the segment;\n* standard deviation on the segment;\n* variance;\n* max;\n* min;  \n* absolute max;  \n* continue component  of FFT transform and mean, standard deviation, max, min for absolute, real and imaginary part of FFT transform.  \n\n"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"abcf886f07660f12d170a2182675e15efffefde0"},"cell_type":"code","source":"rows = 150000\nsegments = int(np.floor(train_df.shape[0] / rows))\nprint(\"Number of segments: \", segments)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"263329e44ea67703470aec40ded07b9f20453fba"},"cell_type":"markdown","source":"We will create additional segments by shifting the current segment with a certain stride. Then, the total number of segments will be:\n\n$$Total_{segments} = N_{segments} * N_{shifts}$$\n\nAnd the size of a shift will be:\n\n$$ Size_{shift} = \\frac{Size_{segment}}{N_{shifts}}$$"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"fc16812019c69d35b3ba223161bfe6e8c6269e79"},"cell_type":"code","source":"shifts = 10\ntotal_segments = segments * shifts\nshift = int(rows/shifts)\nprint(\"Total number of segments: \", total_segments)\nprint(\"Shift size:\", shift)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"602534f8c1698a402a433127d0c850bd1d09b91a"},"cell_type":"code","source":"train_columns = ['mean', 'std', 'var', 'max', 'min', 'abs_max', 'A0', \n                       'Amean', 'Astd', 'Amax', 'Amin',                        \n                       'Rmean', 'Rstd', 'Rmax', 'Rmin',\n                       'Imean', 'Istd', 'Imax', 'Imin']\nX_train = pd.DataFrame(index=range(total_segments), dtype=np.float64,\n                       columns=train_columns)\ny_train = pd.DataFrame(index=range(total_segments), dtype=np.float64,\n                       columns=['time_to_failure'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0bf0608e8c5e65f43b5c885e0a5a8c87a437ae3e"},"cell_type":"code","source":"def create_features(seg_id, seg, X):\n    xc = seg['acoustic_data'].values\n    zc = np.fft.fft(xc)\n        \n    X.loc[seg_id, 'mean'] = xc.mean()\n    X.loc[seg_id, 'std'] = xc.std()\n    X.loc[seg_id, 'var'] = xc.var()\n    X.loc[seg_id, 'max'] = xc.max()\n    X.loc[seg_id, 'min'] = xc.min()\n    X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n    X.loc[seg_id, 'A0'] = abs(zc[0])\n    X.loc[seg_id, 'Amean'] = np.abs(zc).mean()\n    X.loc[seg_id, 'Astd'] = np.abs(zc).std()\n    X.loc[seg_id, 'Amax'] = np.abs(zc).max()\n    X.loc[seg_id, 'Amin'] = np.abs(zc).min()\n    X.loc[seg_id, 'Rmean'] = np.real(zc).mean()\n    X.loc[seg_id, 'Rstd'] = np.real(zc).std()\n    X.loc[seg_id, 'Rmax'] = np.real(zc).max()\n    X.loc[seg_id, 'Rmin'] = np.real(zc).min()\n    X.loc[seg_id, 'Imean'] = np.imag(zc).mean()\n    X.loc[seg_id, 'Istd'] = np.imag(zc).std()\n    X.loc[seg_id, 'Imax'] = np.imag(zc).max()\n    X.loc[seg_id, 'Imin'] = np.imag(zc).min()    ","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"200e364c28f8548519125e0cfdbd63c458e6411d"},"cell_type":"code","source":"for seg_id in tqdm(range(total_segments)):\n    seg = train_df.iloc[seg_id*shift:seg_id*shift+rows]\n    create_features(seg_id, seg, X_train)\n    y_train.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c2e064c25b00c08acd8f8973f565a0b29d09704f"},"cell_type":"markdown","source":"Let's check the result. We plot the head of X_train."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"dc6dd182bf1819c4dcfe823416efca51ff15d5bf"},"cell_type":"code","source":"X_train.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"58196ccf11c59b8bca96279f6bdc7e66eaad6fbe"},"cell_type":"markdown","source":"We scale the data."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"5e9f4c996ac9842fd3243e4f9d3a6e77b8b65cf7"},"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_train)\nX_train_scaled = pd.DataFrame(scaler.transform(X_train), columns=X_train.columns)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"885370309aca6e12d22a811d3ccc126205de468d"},"cell_type":"markdown","source":"Let's check the obtained dataframe."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"1212dd6508332a32bfb4e337f3076ee9680584d4"},"cell_type":"code","source":"X_train_scaled.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"500232cf8ff46fd1d0a9a43088c815d9624a222a"},"cell_type":"markdown","source":"# <a id='5'>Model</a>  \n\nLet's prepare the model.\n\nWe read the submission file and prepare the test file."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"38e67f8e457d7c7ce4900b454d3d4a1ab5c42d04"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(columns=X_train.columns, dtype=np.float64, index=submission.index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c33bd8377aa45cf61d6b3b340a10e81659cab4fe"},"cell_type":"markdown","source":"We apply the same processing done for the training data to the test data."},{"metadata":{"trusted":true,"_uuid":"298c9af7d50367649057a6fc00b64ab391022a37","_kg_hide-input":true},"cell_type":"code","source":"for seg_id in tqdm(X_test.index):\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    create_features(seg_id, seg, X_test)\n      \nX_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7bf4d31bea2f1fd82fbfdc1a4a84f36834daf5c3"},"cell_type":"code","source":"X_test_scaled.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dc3697a05dd856f23fb6084e82fed56915ffc32f"},"cell_type":"markdown","source":"We prepare the folds for cross-validation. \n\nWe will use 5 folds.\n"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"3ec548bc494b168f1afce0c56aefac2414b38860"},"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=2019)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"67d5d1c02c062f2e1b7afecedcff4489f414e206"},"cell_type":"markdown","source":"Model parameters."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"a9eecea2535f548b03c1b9c1af19900621e9cdcc"},"cell_type":"code","source":"params = {'num_leaves': 51,\n         'min_data_in_leaf': 10, \n         'objective':'regression',\n         'max_depth': -1,\n         'learning_rate': 0.05,\n         \"min_child_samples\": 40,\n         \"boosting\": \"gbdt\",\n         \"feature_fraction\": 0.7,\n         \"bagging_freq\": 1,\n         \"bagging_fraction\": 0.9 ,\n         \"bagging_seed\": 11,\n         \"metric\": 'rmse',\n         \"lambda_l1\": 0.1,\n         \"verbosity\": -1,\n         \"nthread\": 4,\n         \"random_state\": 2019}","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"d4eb2efbe227c2ccbfccee3470f5b62327f96103"},"cell_type":"markdown","source":"We prepare the model with cross-validation."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"76d709ed70b11a5b442d1ed17effce7c60bb9091"},"cell_type":"code","source":"oof = np.zeros(len(X_train_scaled))\npredictions = np.zeros(len(X_test_scaled))\nfeature_importance_df = pd.DataFrame()\n#run model\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(X_train_scaled,y_train.values)):\n    strLog = \"fold {}\".format(fold_)\n    print(strLog)\n    \n    X_tr, X_val = X_train_scaled.iloc[trn_idx], X_train_scaled.iloc[val_idx]\n    y_tr, y_val = y_train.iloc[trn_idx], y_train.iloc[val_idx]\n\n    model = lgb.LGBMRegressor(**params, n_estimators = 5000, n_jobs = -1)\n    model.fit(X_tr, y_tr, \n                    eval_set=[(X_tr, y_tr), (X_val, y_val)], eval_metric='mae',\n                    verbose=200, early_stopping_rounds=200)\n    oof[val_idx] = model.predict(X_val, num_iteration=model.best_iteration_)\n    #feature importance\n    fold_importance_df = pd.DataFrame()\n    fold_importance_df[\"Feature\"] = train_columns\n    fold_importance_df[\"importance\"] = model.feature_importances_[:len(train_columns)]\n    fold_importance_df[\"fold\"] = fold_ + 1\n    feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n    #predictions\n    predictions += model.predict(X_test_scaled, num_iteration=model.best_iteration_) / folds.n_splits","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"88e9e2c2b595d76cb57e41ea4cc2dd3f6702822a"},"cell_type":"markdown","source":"We plot the feature importance."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"d1808d0a30e74f610aeba21382ade235bfb3d9c5"},"cell_type":"code","source":"cols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:50].index)\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\nplt.figure(figsize=(6,5))\nsns.barplot(x=\"importance\", y=\"Feature\", data=best_features.sort_values(by=\"importance\",ascending=False))\nplt.title('LightGBM Features (averaged over folds)')\nplt.tight_layout()\nplt.savefig('lgbm_importances.png')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0fce1845bdc118dea09c917d2322c840426e34e1"},"cell_type":"markdown","source":"# <a id='6'>Submission</a>  \n\nWe set the predicted time to failure in the submission file."},{"metadata":{"trusted":true,"_uuid":"641d4feeb57e87e24c3f4bc29863c35169f9302f"},"cell_type":"code","source":"submission['time_to_failure'] = predictions\nsubmission.to_csv('submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"31de440320aa87a962deb1d166abba4ac82e4c04"},"cell_type":"markdown","source":"# <a id='7'>References</a>  \n\n[1] Fast Fourier Transform, https://en.wikipedia.org/wiki/Fast_Fourier_transform   \n[2] Shifting aperture, in Neural network for inverse mapping in eddy current testing, https://www.researchgate.net/publication/3839126_Neural_network_for_inverse_mapping_in_eddy_current_testing "}],"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}