{"cells":[{"metadata":{"_uuid":"3b7efa15014e822ca40ee55d2e7e712a6ca63563"},"cell_type":"markdown","source":"<h1><center><font size=\"6\">LANL Earthquake EDA and Prediction</font></center></h1>\n\n<h2><center><font size=\"4\">Dataset used: LANL Earthquake Prediction</font></center></h2>\n\n<img src=\"https://storage.googleapis.com/kaggle-media/competitions/LANL/nik-shuliahin-585307-unsplash.jpg\" width=\"600\"></img>\n\n<br>\n\n# <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\n## Simulated earthquake experiment\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\n## Competition \n\nIn this competition, the team has provided a much more challenging dataset with considerably more aperiodic earthquake failures.  \nObjective of the competition is to predict the failures for each test set.  \n\n## Kernel\n\nThis solution uses  Andrew's Data Munging plus some inspiration from Scirpus's [Kernel](https://www.kaggle.com/scirpus/andrews-script-plus-a-genetic-program-model/).\n\n"},{"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 time\nimport logging\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport xgboost as xgb\nimport lightgbm as lgb\nfrom scipy import stats\nfrom scipy.signal import hann\nfrom tqdm import tqdm_notebook\nimport matplotlib.pyplot as plt\nfrom scipy.signal import hilbert\nfrom scipy.signal import convolve\nfrom sklearn.svm import NuSVR, SVR\nfrom catboost import CatBoostRegressor\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nwarnings.filterwarnings(\"ignore\")","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":"95ef0a76216010678fb760844307b24719eb4bec"},"cell_type":"markdown","source":"We have two files in the **input** directory and another directory, with the **test** data.  \n\nLet's see how many files are in **test** folder."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"91bda65accbd9172a37799acd2e9198a3477060e"},"cell_type":"code","source":"print(\"There are {} files in test folder\".format(len(os.listdir(os.path.join(PATH, 'test' )))))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68ced52082fc77c768db9ff14fed7cbe43e097d1"},"cell_type":"markdown","source":"\n\nLet's load the train file."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"bb8916c762c26e98e31a721e4aae673aa12e6c5f"},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(os.path.join(PATH,'train.csv'), dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","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","scrolled":true},"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. \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":{"_uuid":"d107cd036c9a931d1e85795f11e5853811a44246"},"cell_type":"markdown","source":"Let's define some computation helper functions."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"8b46e3d6267988a865b69efc1f92ba90a20ca7b7"},"cell_type":"code","source":"def add_trend_feature(arr, abs_values=False):\n    idx = np.array(range(len(arr)))\n    if abs_values:\n        arr = np.abs(arr)\n    lr = LinearRegression()\n    lr.fit(idx.reshape(-1, 1), arr)\n    return lr.coef_[0]\n\ndef classic_sta_lta(x, length_sta, length_lta):\n    sta = np.cumsum(x ** 2)\n    # Convert to float\n    sta = np.require(sta, dtype=np.float)\n    # Copy for LTA\n    lta = sta.copy()\n    # Compute the STA and the LTA\n    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n    sta /= length_sta\n    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n    lta /= length_lta\n    # Pad zeros\n    sta[:length_lta - 1] = 0\n    # Avoid division by zero by setting zero values to tiny float\n    dtiny = np.finfo(0.0).tiny\n    idx = lta < dtiny\n    lta[idx] = dtiny\n    return sta / lta","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0cec83d75c3498355fb232c33ef43dcd2017f24e"},"cell_type":"markdown","source":"## Process train file\n\nNow let's calculate the aggregated functions for train set."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"c65536bf1bf32037175518c3cc8fc7439e217e9d"},"cell_type":"code","source":"train_X = pd.DataFrame(index=range(segments), dtype=np.float64)\ntrain_y = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\ntotal_mean = train_df['acoustic_data'].mean()\ntotal_std = train_df['acoustic_data'].std()\ntotal_max = train_df['acoustic_data'].max()\ntotal_min = train_df['acoustic_data'].min()\ntotal_sum = train_df['acoustic_data'].sum()\ntotal_abs_sum = np.abs(train_df['acoustic_data']).sum()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"0745eb336b1313d3c46c957d7b18c19c98076743"},"cell_type":"code","source":"def create_features(seg_id, seg, X):\n    xc = pd.Series(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, 'max'] = xc.max()\n    X.loc[seg_id, 'min'] = xc.min()\n    \n    #FFT transform values\n    realFFT = np.real(zc)\n    imagFFT = np.imag(zc)\n    X.loc[seg_id, 'Rmean'] = realFFT.mean()\n    X.loc[seg_id, 'Rstd'] = realFFT.std()\n    X.loc[seg_id, 'Rmax'] = realFFT.max()\n    X.loc[seg_id, 'Rmin'] = realFFT.min()\n    X.loc[seg_id, 'Imean'] = imagFFT.mean()\n    X.loc[seg_id, 'Istd'] = imagFFT.std()\n    X.loc[seg_id, 'Imax'] = imagFFT.max()\n    X.loc[seg_id, 'Imin'] = imagFFT.min()\n    X.loc[seg_id, 'Rmean_last_5000'] = realFFT[-5000:].mean()\n    X.loc[seg_id, 'Rstd__last_5000'] = realFFT[-5000:].std()\n    X.loc[seg_id, 'Rmax_last_5000'] = realFFT[-5000:].max()\n    X.loc[seg_id, 'Rmin_last_5000'] = realFFT[-5000:].min()\n    X.loc[seg_id, 'Rmean_last_15000'] = realFFT[-15000:].mean()\n    X.loc[seg_id, 'Rstd_last_15000'] = realFFT[-15000:].std()\n    X.loc[seg_id, 'Rmax_last_15000'] = realFFT[-15000:].max()\n    X.loc[seg_id, 'Rmin_last_15000'] = realFFT[-15000:].min()\n    \n    X.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(xc))\n    X.loc[seg_id, 'mean_change_rate'] = np.mean(np.nonzero((np.diff(xc) / xc[:-1]))[0])\n    X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n    X.loc[seg_id, 'abs_min'] = np.abs(xc).min()\n    \n    X.loc[seg_id, 'std_first_50000'] = xc[:50000].std()\n    X.loc[seg_id, 'std_last_50000'] = xc[-50000:].std()\n    X.loc[seg_id, 'std_first_10000'] = xc[:10000].std()\n    X.loc[seg_id, 'std_last_10000'] = xc[-10000:].std()\n    \n    X.loc[seg_id, 'avg_first_50000'] = xc[:50000].mean()\n    X.loc[seg_id, 'avg_last_50000'] = xc[-50000:].mean()\n    X.loc[seg_id, 'avg_first_10000'] = xc[:10000].mean()\n    X.loc[seg_id, 'avg_last_10000'] = xc[-10000:].mean()\n    \n    X.loc[seg_id, 'min_first_50000'] = xc[:50000].min()\n    X.loc[seg_id, 'min_last_50000'] = xc[-50000:].min()\n    X.loc[seg_id, 'min_first_10000'] = xc[:10000].min()\n    X.loc[seg_id, 'min_last_10000'] = xc[-10000:].min()\n    \n    X.loc[seg_id, 'max_first_50000'] = xc[:50000].max()\n    X.loc[seg_id, 'max_last_50000'] = xc[-50000:].max()\n    X.loc[seg_id, 'max_first_10000'] = xc[:10000].max()\n    X.loc[seg_id, 'max_last_10000'] = xc[-10000:].max()\n    \n    X.loc[seg_id, 'max_to_min'] = xc.max() / np.abs(xc.min())\n    X.loc[seg_id, 'max_to_min_diff'] = xc.max() - np.abs(xc.min())\n    X.loc[seg_id, 'count_big'] = len(xc[np.abs(xc) > 500])\n    X.loc[seg_id, 'sum'] = xc.sum()\n    \n    X.loc[seg_id, 'mean_change_rate_first_50000'] = np.mean(np.nonzero((np.diff(xc[:50000]) / xc[:50000][:-1]))[0])\n    X.loc[seg_id, 'mean_change_rate_last_50000'] = np.mean(np.nonzero((np.diff(xc[-50000:]) / xc[-50000:][:-1]))[0])\n    X.loc[seg_id, 'mean_change_rate_first_10000'] = np.mean(np.nonzero((np.diff(xc[:10000]) / xc[:10000][:-1]))[0])\n    X.loc[seg_id, 'mean_change_rate_last_10000'] = np.mean(np.nonzero((np.diff(xc[-10000:]) / xc[-10000:][:-1]))[0])\n    \n    X.loc[seg_id, 'q95'] = np.quantile(xc, 0.95)\n    X.loc[seg_id, 'q99'] = np.quantile(xc, 0.99)\n    X.loc[seg_id, 'q05'] = np.quantile(xc, 0.05)\n    X.loc[seg_id, 'q01'] = np.quantile(xc, 0.01)\n    \n    X.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(xc), 0.95)\n    X.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(xc), 0.99)\n    X.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(xc), 0.05)\n    X.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(xc), 0.01)\n    \n    X.loc[seg_id, 'trend'] = add_trend_feature(xc)\n    X.loc[seg_id, 'abs_trend'] = add_trend_feature(xc, abs_values=True)\n    X.loc[seg_id, 'abs_mean'] = np.abs(xc).mean()\n    X.loc[seg_id, 'abs_std'] = np.abs(xc).std()\n    \n    X.loc[seg_id, 'mad'] = xc.mad()\n    X.loc[seg_id, 'kurt'] = xc.kurtosis()\n    X.loc[seg_id, 'skew'] = xc.skew()\n    X.loc[seg_id, 'med'] = xc.median()\n    \n    X.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(xc)).mean()\n    X.loc[seg_id, 'Hann_window_mean'] = (convolve(xc, hann(150), mode='same') / sum(hann(150))).mean()\n    X.loc[seg_id, 'classic_sta_lta1_mean'] = classic_sta_lta(xc, 500, 10000).mean()\n    X.loc[seg_id, 'classic_sta_lta2_mean'] = classic_sta_lta(xc, 5000, 100000).mean()\n    X.loc[seg_id, 'classic_sta_lta3_mean'] = classic_sta_lta(xc, 3333, 6666).mean()\n    X.loc[seg_id, 'classic_sta_lta4_mean'] = classic_sta_lta(xc, 10000, 25000).mean()\n    X.loc[seg_id, 'Moving_average_700_mean'] = xc.rolling(window=700).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_1500_mean'] = xc.rolling(window=1500).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_3000_mean'] = xc.rolling(window=3000).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_6000_mean'] = xc.rolling(window=6000).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    X.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(xc, span=300).mean()).mean(skipna=True)\n    X.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(xc, span=3000).mean().mean(skipna=True)\n    X.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(xc, span=6000).mean().mean(skipna=True)\n    no_of_std = 2\n    X.loc[seg_id, 'MA_700MA_std_mean'] = xc.rolling(window=700).std().mean()\n    X.loc[seg_id,'MA_700MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X.loc[seg_id,'MA_700MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X.loc[seg_id, 'MA_400MA_std_mean'] = xc.rolling(window=400).std().mean()\n    X.loc[seg_id,'MA_400MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X.loc[seg_id,'MA_400MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X.loc[seg_id, 'MA_1000MA_std_mean'] = xc.rolling(window=1000).std().mean()\n    \n    X.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(xc, [75, 25]))\n    X.loc[seg_id, 'q999'] = np.quantile(xc,0.999)\n    X.loc[seg_id, 'q001'] = np.quantile(xc,0.001)\n    X.loc[seg_id, 'ave10'] = stats.trim_mean(xc, 0.1)\n    \n    for windows in [10, 100, 1000]:\n        x_roll_std = xc.rolling(windows).std().dropna().values\n        x_roll_mean = xc.rolling(windows).mean().dropna().values\n        \n        X.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X.loc[seg_id, 'av_change_rate_roll_std_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n        X.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X.loc[seg_id, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n        X.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a446450bfd25296d6b98dde3e750e11d96250a4","_kg_hide-input":true},"cell_type":"code","source":"# iterate over all segments\nfor seg_id in tqdm_notebook(range(segments)):\n    seg = train_df.iloc[seg_id*rows:seg_id*rows+rows]\n    create_features(seg_id, seg, train_X)\n    train_y.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 shape and the head of train_X."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"ff98b8f65290cfa14a5fed507a5567a589109fc9"},"cell_type":"code","source":"train_X.shape","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"dc6dd182bf1819c4dcfe823416efca51ff15d5bf"},"cell_type":"code","source":"train_X.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(train_X)\nscaled_train_X = pd.DataFrame(scaler.transform(train_X), columns=train_X.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":"scaled_train_X.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"500232cf8ff46fd1d0a9a43088c815d9624a222a"},"cell_type":"markdown","source":"## Process test data\n\nWe apply the same processing done for the training data to the test data.\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')\ntest_X = pd.DataFrame(columns=train_X.columns, dtype=np.float64, index=submission.index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b9b80332f6e9af2ee82f8e0feca4e9d1d18925ee"},"cell_type":"markdown","source":"Let's check the shape of the submission and test_X datasets."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"1f06fee3edb09ac05890fa109808487de395d1de"},"cell_type":"code","source":"submission.shape, test_X.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e3b54cd654438c2af63060b03ddb767eacbd70c6","_kg_hide-input":true},"cell_type":"code","source":"for seg_id in tqdm_notebook(test_X.index):\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    create_features(seg_id, seg, test_X)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5630b8d9234af698ed83af94980538e9ed8eaa59"},"cell_type":"markdown","source":"We scale also the test data."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"a18c26409ced2a71c3b7b4926d05455afd5fa26a"},"cell_type":"code","source":"scaled_test_X = pd.DataFrame(scaler.transform(test_X), columns=test_X.columns)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"9c5701e9cc1c9076bde818af54fa67ec4ad3eb4b"},"cell_type":"code","source":"scaled_test_X.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7bf4d31bea2f1fd82fbfdc1a4a84f36834daf5c3","_kg_hide-input":true},"cell_type":"code","source":"scaled_test_X.tail(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2bb7ba0d1ec723be641adf5f925c24237c9198c0"},"cell_type":"markdown","source":"\n# <a id='5'>Model</a>  \n\nLet's prepare the model.\n\n## Run model\n\nWe define the folds for cross-validation."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"99d43c789dde8cac37270dc16430cf98870788d3"},"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\ntrain_columns = scaled_train_X.columns.values","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d506ac30835fbfd5ed324bab3a3370e01739f907"},"cell_type":"markdown","source":"We define the model parameters."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"b3fc19ff1bf1bb0897e29747421514f211ea0edc"},"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.001,\n         \"boosting\": \"gbdt\",\n         \"feature_fraction\": 0.91,\n         \"bagging_freq\": 1,\n         \"bagging_fraction\": 0.91,\n         \"bagging_seed\": 42,\n         \"metric\": 'mae',\n         \"lambda_l1\": 0.1,\n         \"verbosity\": -1,\n         \"nthread\": -1,\n         \"random_state\": 42}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f5e691c0c52f7adfb13064794c7bf2600f07f39"},"cell_type":"markdown","source":"We run the model. During training for each fold, we validate using the validation set and also we predict using the current model for test set. The final result will be the average over the all folds for the predictions done at each fold training."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"d0b278f0a14dbd3dee610e2e27a0cc6d55844678"},"cell_type":"code","source":"oof = np.zeros(len(scaled_train_X))\npredictions = np.zeros(len(scaled_test_X))\nfeature_importance_df = pd.DataFrame()\n#run model\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(scaled_train_X,train_y.values)):\n    strLog = \"fold {}\".format(fold_)\n    print(strLog)\n    \n    X_tr, X_val = scaled_train_X.iloc[trn_idx], scaled_train_X.iloc[val_idx]\n    y_tr, y_val = train_y.iloc[trn_idx], train_y.iloc[val_idx]\n\n    model = lgb.LGBMRegressor(**params, n_estimators = 20000, n_jobs = -1)\n    model.fit(X_tr, \n              y_tr, \n              eval_set=[(X_tr, y_tr), (X_val, y_val)], \n              eval_metric='mae',\n              verbose=1000, \n              early_stopping_rounds=500)\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(scaled_test_X, num_iteration=model.best_iteration_) / folds.n_splits","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cc8966810c47ae092a300dbb3a5ec582c9740b8f"},"cell_type":"markdown","source":"## Features importance\n\nLet's print features importance."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"0f64306227bc30ffa306b5821a63f04ab6a874c0"},"cell_type":"code","source":"cols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:200].index)\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\n\nplt.figure(figsize=(14,26))\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","_kg_hide-input":true},"cell_type":"code","source":"submission.time_to_failure = predictions\nsubmission.to_csv('submission.csv',index=True)","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   \n[3] Andrews Script plus a Genetic Program Model, https://www.kaggle.com/scirpus/andrews-script-plus-a-genetic-program-model/"}],"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}