{"cells":[{"metadata":{"_uuid":"3b7efa15014e822ca40ee55d2e7e712a6ca63563"},"cell_type":"markdown","source":"<h1><center><font size=\"6\">LANL Earthquake EDA and Ensemble 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\nMy goal with this kernel is to have an end-to-end model which encorporates:\n\n* Extensible feature engineering\n* Ability to try and measure many regression models\n* Fexlibility of both a blended and stacked model production to produce testable ensemble models\n\nThis kernel is based on a fork from Gabriel Preda's excellent kernel, 'LANLEarthquake EDA and Prediction.\" I have left most of his EDA and feature engineering in place, while adding a few features of my own. The modeling portion is based partially on his models. In addition, I use ensembling approaches from Serigne's excellent model in the House Price Advanced Regression Techniques competition, https://www.kaggle.com/serigne/stacked-regressions-top-4-on-leaderboard.\n\nThis model produces 1.488 leader board entry with the averaged model, using RF, lgb, cat, and lasso.\n\n\n## Simulated earthquake experiment\n\nThis introduction is from Gabrial Prada's kernel:\n\n\"The 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 is based on a fork from Gabriel Preda's excellent kernel, 'LANLEarthquake EDA and Prediction.\" His 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\nIt also includes the elegant stacking model from Serigne's regression model (see reference above.)\n\nThis kernel extends the existing model by\n\n* normalizing the values to a mean of zero\n* using log values to remove outliers\n* removal of segments in training which contain the quake\n* additional prediction models\n* combining results from multiple models through various ensemble techniques\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":false},"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.linear_model import ElasticNet, Lasso,  BayesianRidge, LassoLarsIC\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.ensemble import RandomForestRegressor,  GradientBoostingRegressor\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nfrom sklearn.model_selection import cross_val_score, train_test_split, cross_val_predict\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import RobustScaler\nfrom sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone\n\n\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":false,"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"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":false,"trusted":true,"_uuid":"bb8916c762c26e98e31a721e4aae673aa12e6c5f"},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(os.path.join(PATH,'train.csv'), \n                       dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df_save = train_df.copy","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. We normalize the data to a 0 mean value."},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_acoustic = np.mean(train_df.acoustic_data)\ntrain_df.acoustic_data = train_df.acoustic_data - mean_acoustic\nprint(mean_acoustic)\nprint (round(np.mean(train_df.acoustic_data),2))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"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":false,"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 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 2% of the data. For this we will sample every 50 points of data.  "},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"206af56767581d8e98777da3b0670e11d4e444bd","scrolled":true},"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[::50]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[::50]\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 2% 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. Because the data is all positive and log, the peaks are not as severe.\n\nWe 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 two earthquakes of the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"845e5a241e545d73cfb1dbaea2cfa464faf8148a"},"cell_type":"code","source":"gc.collect()\n#train_ad_sample_df = train_df['acoustic_data'].values[:6291455]\n#train_ttf_sample_df = (train_df['time_to_failure'].values[:6291455])\ntrain_ad_sample_df = train_df['acoustic_data'].values[:50580000]\ntrain_ttf_sample_df = (train_df['time_to_failure'].values[:50580000])\n\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1st 2 quakes\")\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":false,"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":false,"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":{"trusted":true},"cell_type":"code","source":"gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"c65536bf1bf32037175518c3cc8fc7439e217e9d"},"cell_type":"code","source":"train_X = pd.DataFrame(index=range(segments), dtype=np.float64)\n\ntrain_y = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\n# These may be needed later\nneed_aggregated_features = False\nif need_aggregated_features:\n    total_mean = train_df['acoustic_data'].mean()\n    total_std = train_df['acoustic_data'].std()\n    total_max = train_df['acoustic_data'].max()\n    total_min = train_df['acoustic_data'].min()\n    total_sum = train_df['acoustic_data'].sum()\n    total_abs_sum = np.abs(train_df['acoustic_data']).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X.shape, train_y.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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    X.loc[seg_id, 'mean'] = xc.mean()\n    X.loc[seg_id, 'std'] = xc.std()\n    X.loc[seg_id, 'max'] = xc.max()\n    \n    #FFT transform values\n    realFFT = np.real(zc)\n    imagFFT = np.imag(zc)\n    \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    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    \n    X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n    X.loc[seg_id, 'abs_min'] = np.abs(xc).min()\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    \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    \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 [5, 10, 50, 100, 500, 1000, 5000, 10000]:\n    \n        x_roll_std = xc.rolling(windows).std().dropna().values\n        x_roll_mean = xc.rolling(windows).mean().dropna().values\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()  \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a446450bfd25296d6b98dde3e750e11d96250a4"},"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    # the y value is the last entry in the time to failure in the segment\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 and train_y."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X_save = train_X.copy\ntrain_y_save = train_y.copy\ntrain_y.head(5)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We will not train on the segments with a quake, because there are likely outliers\ntrain_y_quake = np.nonzero(np.diff(train_y.time_to_failure) > 0)[0] + 1\nprint(len(train_y_quake))\nprint (len(train_y))\n\nfor idx in train_y_quake: \n    train_y.drop([idx],inplace=True)\n    train_X.drop([idx],inplace = True)\n#np.abs(train_X.corrwith(train_y)).sort_values(ascending=False).head(12)\ntrain_X.to_csv('train_features.csv', index=False)\ntrain_y.to_csv('train_y.csv', index=False)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"ff98b8f65290cfa14a5fed507a5567a589109fc9"},"cell_type":"code","source":"train_X.shape, train_y.shape\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X.head(), train_y.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"58196ccf11c59b8bca96279f6bdc7e66eaad6fbe"},"cell_type":"markdown","source":"We scale the data."},{"metadata":{"_kg_hide-input":false,"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)\n#scaled_train_X = train_X","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"885370309aca6e12d22a811d3ccc126205de468d"},"cell_type":"markdown","source":"Let's check the obtained dataframe."},{"metadata":{"_kg_hide-input":false,"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":false,"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":false,"trusted":true,"_uuid":"1f06fee3edb09ac05890fa109808487de395d1de"},"cell_type":"code","source":"submission.shape, test_X.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e3b54cd654438c2af63060b03ddb767eacbd70c6"},"cell_type":"code","source":"for seg_id in tqdm_notebook(test_X.index):\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    # convert to mean 0 of the training dataset \n    # seg_mean = np.mean(seg.acoustic_data)\n    seg.acoustic_data = seg.acoustic_data - mean_acoustic\n    create_features(seg_id, seg, test_X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# save before scaling\ntest_X.to_csv('test_features.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5630b8d9234af698ed83af94980538e9ed8eaa59"},"cell_type":"markdown","source":"We scale the test data."},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"a18c26409ced2a71c3b7b4926d05455afd5fa26a"},"cell_type":"code","source":"scaled_test_X = pd.DataFrame(scaler.transform(test_X), columns=test_X.columns)\n#scaled_test_X = test_X\nscaled_test_X.values[1117]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"9c5701e9cc1c9076bde818af54fa67ec4ad3eb4b"},"cell_type":"code","source":"scaled_test_X.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7bf4d31bea2f1fd82fbfdc1a4a84f36834daf5c3","_kg_hide-input":false},"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. First we define a validation function to evaluate model performance.\n\nWe next define a set of models which can later be blended and stacked.\n\n\n## Define Validation\n\nFirst define a cross validation measurement routine."},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 5\ndef mae_cv (model):\n    folds = KFold(n_splits=n_fold, shuffle=True, random_state=42).get_n_splits(scaled_train_X.values)\n    mae = -cross_val_score (model, scaled_train_X.values, train_y, scoring=\"neg_mean_absolute_error\",\n                           verbose=0,\n                           cv=folds)\n    return mae\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d506ac30835fbfd5ed324bab3a3370e01739f907"},"cell_type":"markdown","source":"### LGB Model\nWe define the model parameters and the definition for the first model, an LGB regression."},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"b3fc19ff1bf1bb0897e29747421514f211ea0edc"},"cell_type":"code","source":"%%time\n\nlgb_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}\n\n\nlgb_model = lgb.LGBMRegressor(objective='regression',num_leaves=5,\n                              learning_rate=0.01, n_estimators=720,\n                              bagging_freq = 5, feature_fraction = 0.2319,\n                              feature_fraction_seed=9, bagging_seed=9,\n                              min_data_in_leaf =6, min_sum_hessian_in_leaf = 11, n_jobs = -1)\n\nscore = mae_cv(lgb_model)\nprint(\"LGBM score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nlgb_model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lgb_gamma_model = lgb.LGBMRegressor(objective='gamma',num_leaves=5,\n                              learning_rate=0.01, n_estimators=720,\n                              bagging_freq = 5, feature_fraction = 0.2319,\n                              feature_fraction_seed=9, bagging_seed=9,\n                              min_data_in_leaf =6, min_sum_hessian_in_leaf = 11, n_jobs = -1)\n\nscore = mae_cv(lgb_gamma_model)\nprint(\"LGBM - gamma score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nlgb_gamma_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### XGB Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time \nxgb_params = {'eta': 0.03,\n              'max_depth': 9,\n              'subsample': 0.85,\n              'objective': 'reg:linear',\n              'eval_metric': 'mae',\n              'silent': True,\n              'nthread': 4}\n    \nxgb_model = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, \n                             learning_rate=0.05, max_depth=3, \n                             min_child_weight=1.7817, n_estimators=2200,\n                             reg_alpha=0.4640, reg_lambda=0.8571,\n                             subsample=0.5213, silent=1,\n                             random_state =7, nthread = -1, eval_metric = 'mae',)\n\nscore = mae_cv(xgb_model)\nprint(\"XGB score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nxgb_model\n\n#    xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, \n#                          verbose_eval=500, params=xgb_params)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f5e691c0c52f7adfb13064794c7bf2600f07f39"},"cell_type":"markdown","source":"### Random Forest"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nrf_model = RandomForestRegressor(n_estimators=120, n_jobs=-1, min_samples_leaf=1, \n                           max_features = \"auto\",max_depth=15, )\nscore = mae_cv(rf_model)\nprint(\"Random Forest score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nrf_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Cat Boost\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nparams = {'loss_function':'MAE',}\ncat_model = CatBoostRegressor(iterations=1000,  eval_metric='MAE', verbose=False, **params)\n\nscore = mae_cv(cat_model)\nprint(\"Cat Boost score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\ncat_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Kernel Ridge"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nKRR_model = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)\nscore = mae_cv(KRR_model)\nprint(\"Kernel Ridge score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nprint (score)\nKRR_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Elastic Net\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n#ENet_model = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=0.9, random_state=3,max_iter=5000))\nENet_model = ElasticNet(alpha=0.0005, l1_ratio=0.9, random_state=3,max_iter=5000)\nscore = mae_cv(ENet_model)\nprint(\"Elastic Net score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nENet_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Lasso"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nlasso_model = Lasso(alpha =0.0005, random_state=1)\nscore = mae_cv(lasso_model)\nprint(\"Lasso score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nlasso_model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):\n    def __init__(self, models):\n        self.models = models\n        \n    # we define clones of the original models to fit the data in\n    def fit(self, X, y):\n        self.models_ = [clone(x) for x in self.models]\n        \n        # Train cloned base models\n        for model in self.models_:\n            model.fit(X, y)\n\n        return self\n    \n    #Now we do the predictions for cloned models and average them\n    def predict(self, X):\n        \n        predictions = np.column_stack([\n            model.predict(X) for model in self.models_\n        ])\n        return np.mean(predictions, axis=1)   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):\n    def __init__(self, base_models, meta_model, n_folds=5):\n        self.base_models = base_models\n        self.meta_model = meta_model\n        self.n_folds = n_folds\n   \n    # We again fit the data on clones of the original models\n    def fit(self, X, y):\n        print (type(X))\n        self.base_models_ = [list() for x in self.base_models]\n        self.meta_model_ = clone(self.meta_model)\n        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)\n        print (KFold)\n        # Train cloned base models then create out-of-fold predictions\n        # that are needed to train the cloned meta-model\n        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))\n        for i, model in enumerate(self.base_models):\n            for train_index, holdout_index in kfold.split(X, y):\n                instance = clone(model)\n                self.base_models_[i].append(instance)\n                instance.fit(X.iloc[train_index], y.iloc[train_index])\n                y_pred = instance.predict(X.iloc[holdout_index])\n                out_of_fold_predictions[holdout_index, i] = y_pred\n                \n        # Now train the cloned  meta-model using the out-of-fold predictions as new feature\n        self.meta_model_.fit(out_of_fold_predictions, y)\n        return self\n   \n    #Do the predictions of all base models on the test data and use the averaged predictions as \n    #meta-features for the final prediction which is done by the meta-model\n    def predict(self, X):\n        meta_features = np.column_stack([\n            np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)\n            for base_models in self.base_models_ ])\n        return self.meta_model_.predict(meta_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class StackingCVRegressorRetrained(BaseEstimator, RegressorMixin, TransformerMixin):\n    def __init__(self, regressors, meta_regressor, n_folds=5, use_features_in_secondary=False):\n        self.regressors = regressors\n        self.meta_regressor = meta_regressor\n        self.n_folds = n_folds\n        self.use_features_in_secondary = use_features_in_secondary\n\n    def fit(self, X, y):\n        self.regr_ = [clone(x) for x in self.regressors]\n        self.meta_regr_ = clone(self.meta_regressor)\n\n        kfold = KFold(n_splits=self.n_folds, shuffle=True)\n\n        out_of_fold_predictions = np.zeros((X.shape[0], len(self.regressors)))\n\n        # Create out-of-fold predictions for training meta-model\n        for i, regr in enumerate(self.regr_):\n            for train_idx, holdout_idx in kfold.split(X, y):\n                instance = clone(regr)\n                instance.fit(X[train_idx], y[train_idx])\n                out_of_fold_predictions[holdout_idx, i] = instance.predict(X[holdout_idx])\n\n        # Train meta-model\n        if self.use_features_in_secondary:\n            self.meta_regr_.fit(np.hstack((X, out_of_fold_predictions)), y)\n        else:\n            self.meta_regr_.fit(out_of_fold_predictions, y)\n        \n        # Retrain base models on all data\n        for regr in self.regr_:\n            regr.fit(X, y)\n\n        return self\n\n    def predict(self, X):\n        meta_features = np.column_stack([\n            regr.predict(X) for regr in self.regr_\n        ])\n\n        if self.use_features_in_secondary:\n            return self.meta_regr_.predict(np.hstack((X, meta_features)))\n        else:\n            return self.meta_regr_.predict(meta_features)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n#averaged_models = AveragingModels(models = (rf_model, xgb_model, KRR_model, lgb_model, ENet_model, cat_model, lasso_model))\n#averaged_models = AveragingModels(models = (rf_model, lgb_model,  cat_model, lasso_model))\naveraged_models = AveragingModels(models = (rf_model,cat_model))\n\nscore =mae_cv(averaged_models)\nprint(\" Averaged base models score: {:.4f} ({:.4f})\\n\".format(score.mean(), score.std()))\n#averaged_models.fit (scaled_train_X.values, train_y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"averaged_models.fit (scaled_train_X.values, train_y)\naveraged2_train_predict = averaged_models.predict(scaled_train_X.values)\nprint(mean_absolute_error(train_y, averaged2_train_predict))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\naveraged_prediction = np.zeros(len(scaled_test_X))\naveraged_prediction += averaged_models.predict(scaled_test_X.values)\naveraged_prediction","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nstacked_predict = StackingAveragedModels(base_models =(rf_model, xgb_model, lgb_model, cat_model,ENet_model), \n                                          meta_model =lasso_model) \nstacked_predict.fit(scaled_train_X, train_y)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stacked_train_pred = stacked_predict.predict(scaled_train_X)\n\nprint(mean_absolute_error(train_y, stacked_train_pred))\n\nstacked_prediction = np.zeros(len(scaled_test_X))\nstacked_prediction += stacked_predict.predict(scaled_test_X)**1.0\nstacked_prediction[0:4]\n","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":false},"cell_type":"code","source":"submission.time_to_failure = averaged_prediction\nsubmission.to_csv('submissionV30_averaged_cat_rf.csv',index=True)\nsubmission.time_to_failure = stacked_prediction\nsubmission.to_csv('submissionV30_stacked.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}