{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"## General information\n\nCorrectly predicting earthquakes is very important for preventing deaths and damage to infrastructure. In this competition we try to predict time left to the next laboratory earthquake based on seismic signal data.\nTraining data represents one huge signal, but in test data we have many separate chunks, for each of which we need to predict time to failure.\n\nThis is my second kernel for this competition, here is the [link](https://www.kaggle.com/artgor/seismic-data-eda-and-baseline) to the first one.\n\nIn this kernel I'll try to create more useful features and generate more data for training.\n\n![](https://i.cbc.ca/1.4972912.1547133821!/fileImage/httpImage/image.jpg_gen/derivatives/16x9_780/new-brunswick-earthquake.jpg)"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm_notebook\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.metrics import mean_absolute_error\npd.options.display.precision = 15\n\nimport lightgbm as lgb\nimport xgboost as xgb\nimport time\nimport datetime\nfrom catboost import CatBoostRegressor\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.linear_model import LinearRegression\nimport gc\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nfrom scipy.signal import hilbert\nfrom scipy.signal import hann\nfrom scipy.signal import convolve\nfrom scipy import stats\nfrom sklearn.kernel_ridge import KernelRidge","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bc314853a64c3ee3b7b14dd2ad99c06e77e16561"},"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":{"_kg_hide-input":true,"trusted":true,"_uuid":"029a1f0d374a39d7fac9a7c81d417973bce1c884"},"cell_type":"code","source":"train_acoustic_data_small = train['acoustic_data'].values[::50]\ntrain_time_to_failure_small = train['time_to_failure'].values[::50]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"70a699101fd731a3d0963a7c71b8229510bfe67f"},"cell_type":"markdown","source":"* We can see that usually acoustic data shows huge fluctuations just before the failure and the nature of data is cyclical;\n* Another important point: visually failures can be predicted as cases when huge fluctuations in signal are followes by small signal values. This could be useful for predicting \"time_to_failure\" changes from 0 to high values;\n* I thought that comparing max values of signal in a segment to some threshold value (1000 or 2000) could be useful, but it didn't work;"},{"metadata":{"_uuid":"8fbeddf32663d5a5dd3399c3ed19271bbc937c92"},"cell_type":"markdown","source":"## Feature generation\n\nI create several groups of features:\n* Usual aggregations: mean, std, min and max;\n* Average difference between the consequitive values in absolute and percent values;\n* Absolute min and max vallues;\n* Aforementioned aggregations for first and last 10000 and 50000 values - I think these data should be useful;\n* Max value to min value and their differencem also count of values bigger that 500 (arbitrary threshold);\n* Quantile features from this kernel: https://www.kaggle.com/andrekos/basic-feature-benchmark-with-quantiles\n* Trend features from this kernel: https://www.kaggle.com/jsaguiar/baseline-with-abs-and-trend-features\n* Rolling features from this kernel: https://www.kaggle.com/wimwim/rolling-quantiles"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"e75947cef3bff46114dd09ce164625871d54b30e","scrolled":true},"cell_type":"code","source":"# Create a training file with simple derived features\nrows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\ndef 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    \n    sta = np.cumsum(x ** 2)\n\n    # Convert to float\n    sta = np.require(sta, dtype=np.float)\n\n    # Copy for LTA\n    lta = sta.copy()\n\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\n    # Pad zeros\n    sta[:length_lta - 1] = 0\n\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\n    return sta / lta\n\nX_tr = pd.DataFrame(index=range(segments), dtype=np.float64)\n\ny_tr = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\n\ntotal_mean = train['acoustic_data'].mean()\ntotal_std = train['acoustic_data'].std()\ntotal_max = train['acoustic_data'].max()\ntotal_min = train['acoustic_data'].min()\ntotal_sum = train['acoustic_data'].sum()\ntotal_abs_sum = np.abs(train['acoustic_data']).sum()\n\ndef calc_change_rate(x):\n    change = (np.diff(x) / x[:-1]).values\n    change = change[np.nonzero(change)[0]]\n    change = change[~np.isnan(change)]\n    change = change[change != -np.inf]\n    change = change[change != np.inf]\n    return np.mean(change)\n\nfor segment in tqdm_notebook(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n    x = pd.Series(seg['acoustic_data'].values)\n    y = seg['time_to_failure'].values[-1]\n    \n    y_tr.loc[segment, 'time_to_failure'] = y\n    X_tr.loc[segment, 'mean'] = x.mean()\n    X_tr.loc[segment, 'std'] = x.std()\n    X_tr.loc[segment, 'max'] = x.max()\n    X_tr.loc[segment, 'min'] = x.min()\n    \n    X_tr.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x))\n    X_tr.loc[segment, 'mean_change_rate'] = calc_change_rate(x)\n    X_tr.loc[segment, 'abs_max'] = np.abs(x).max()\n    X_tr.loc[segment, 'abs_min'] = np.abs(x).min()\n    \n    X_tr.loc[segment, 'std_first_50000'] = x[:50000].std()\n    X_tr.loc[segment, 'std_last_50000'] = x[-50000:].std()\n    X_tr.loc[segment, 'std_first_10000'] = x[:10000].std()\n    X_tr.loc[segment, 'std_last_10000'] = x[-10000:].std()\n    \n    X_tr.loc[segment, 'avg_first_50000'] = x[:50000].mean()\n    X_tr.loc[segment, 'avg_last_50000'] = x[-50000:].mean()\n    X_tr.loc[segment, 'avg_first_10000'] = x[:10000].mean()\n    X_tr.loc[segment, 'avg_last_10000'] = x[-10000:].mean()\n    \n    X_tr.loc[segment, 'min_first_50000'] = x[:50000].min()\n    X_tr.loc[segment, 'min_last_50000'] = x[-50000:].min()\n    X_tr.loc[segment, 'min_first_10000'] = x[:10000].min()\n    X_tr.loc[segment, 'min_last_10000'] = x[-10000:].min()\n    \n    X_tr.loc[segment, 'max_first_50000'] = x[:50000].max()\n    X_tr.loc[segment, 'max_last_50000'] = x[-50000:].max()\n    X_tr.loc[segment, 'max_first_10000'] = x[:10000].max()\n    X_tr.loc[segment, 'max_last_10000'] = x[-10000:].max()\n    \n    X_tr.loc[segment, 'max_to_min'] = x.max() / np.abs(x.min())\n    X_tr.loc[segment, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n    X_tr.loc[segment, 'count_big'] = len(x[np.abs(x) > 500])\n    X_tr.loc[segment, 'sum'] = x.sum()\n    \n    X_tr.loc[segment, 'mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n    X_tr.loc[segment, 'mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])\n    X_tr.loc[segment, 'mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n    X_tr.loc[segment, 'mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])\n    \n    X_tr.loc[segment, 'q95'] = np.quantile(x, 0.95)\n    X_tr.loc[segment, 'q99'] = np.quantile(x, 0.99)\n    X_tr.loc[segment, 'q05'] = np.quantile(x, 0.05)\n    X_tr.loc[segment, 'q01'] = np.quantile(x, 0.01)\n    \n    X_tr.loc[segment, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n    X_tr.loc[segment, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n    X_tr.loc[segment, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n    X_tr.loc[segment, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    X_tr.loc[segment, 'trend'] = add_trend_feature(x)\n    X_tr.loc[segment, 'abs_trend'] = add_trend_feature(x, abs_values=True)\n    X_tr.loc[segment, 'abs_mean'] = np.abs(x).mean()\n    X_tr.loc[segment, 'abs_std'] = np.abs(x).std()\n    \n    X_tr.loc[segment, 'mad'] = x.mad()\n    X_tr.loc[segment, 'kurt'] = x.kurtosis()\n    X_tr.loc[segment, 'skew'] = x.skew()\n    X_tr.loc[segment, 'med'] = x.median()\n    \n    X_tr.loc[segment, 'Hilbert_mean'] = np.abs(hilbert(x)).mean()\n    X_tr.loc[segment, 'Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n    X_tr.loc[segment, 'classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n    X_tr.loc[segment, 'classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n    X_tr.loc[segment, 'classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n    X_tr.loc[segment, 'classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n    X_tr.loc[segment, 'classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean()\n    X_tr.loc[segment, 'classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n    X_tr.loc[segment, 'classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean()\n    X_tr.loc[segment, 'classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n    X_tr.loc[segment, 'Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    X_tr.loc[segment, 'exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n    X_tr.loc[segment, 'exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n    X_tr.loc[segment, 'exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n    no_of_std = 3\n    X_tr.loc[segment, 'MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n    X_tr.loc[segment,'MA_700MA_BB_high_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] + no_of_std * X_tr.loc[segment, 'MA_700MA_std_mean']).mean()\n    X_tr.loc[segment,'MA_700MA_BB_low_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] - no_of_std * X_tr.loc[segment, 'MA_700MA_std_mean']).mean()\n    X_tr.loc[segment, 'MA_400MA_std_mean'] = x.rolling(window=400).std().mean()\n    X_tr.loc[segment,'MA_400MA_BB_high_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] + no_of_std * X_tr.loc[segment, 'MA_400MA_std_mean']).mean()\n    X_tr.loc[segment,'MA_400MA_BB_low_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] - no_of_std * X_tr.loc[segment, 'MA_400MA_std_mean']).mean()\n    X_tr.loc[segment, 'MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n    X_tr.drop('Moving_average_700_mean', axis=1, inplace=True)\n    \n    X_tr.loc[segment, 'iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n    X_tr.loc[segment, 'q999'] = np.quantile(x,0.999)\n    X_tr.loc[segment, 'q001'] = np.quantile(x,0.001)\n    X_tr.loc[segment, 'ave10'] = stats.trim_mean(x, 0.1)\n\n    for windows in [10, 100, 1000]:\n        x_roll_std = x.rolling(windows).std().dropna().values\n        x_roll_mean = x.rolling(windows).mean().dropna().values\n        \n        X_tr.loc[segment, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X_tr.loc[segment, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X_tr.loc[segment, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X_tr.loc[segment, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X_tr.loc[segment, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X_tr.loc[segment, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X_tr.loc[segment, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X_tr.loc[segment, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X_tr.loc[segment, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X_tr.loc[segment, 'av_change_rate_roll_std_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n        X_tr.loc[segment, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X_tr.loc[segment, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X_tr.loc[segment, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X_tr.loc[segment, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X_tr.loc[segment, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X_tr.loc[segment, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X_tr.loc[segment, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X_tr.loc[segment, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X_tr.loc[segment, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X_tr.loc[segment, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X_tr.loc[segment, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n        X_tr.loc[segment, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6ad94de98384e95f8a769ce9670d8b55f636eec3"},"cell_type":"code","source":"print(f'{X_tr.shape[0]} samples in new train data and {X_tr.shape[1]} columns.')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25734169d50198f51f5fe6fff44caac7967d171d"},"cell_type":"markdown","source":"### On sampling\n\nI tried to randomly sample 150000 rows 1k-10k times and add these samples to training data, but it severely decreased my score."},{"metadata":{"trusted":true,"_uuid":"b3e87358c38c3dc2a675be51a241d8bfa7171bbd"},"cell_type":"code","source":"np.abs(X_tr.corrwith(y_tr['time_to_failure'])).sort_values(ascending=False).head(12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a2bd1c5cd7178d083b99b19a226e89486aa02d8b"},"cell_type":"code","source":"plt.figure(figsize=(44, 24))\ncols = list(np.abs(X_tr.corrwith(y_tr['time_to_failure'])).sort_values(ascending=False).head(24).index)\nfor i, col in enumerate(cols):\n    plt.subplot(6, 4, i + 1)\n    plt.plot(X_tr[col], color='blue')\n    plt.title(col)\n    ax1.set_ylabel(col, color='b')\n\n    ax2 = ax1.twinx()\n    plt.plot(y_tr, color='g')\n    ax2.set_ylabel('time_to_failure', color='g')\n    plt.legend([col, 'time_to_failure'], loc=(0.875, 0.9))\n    plt.grid(False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"02a2b2d6f64393ea928e60e4e660e47bec1afc16"},"cell_type":"code","source":"means_dict = {}\nfor col in X_tr.columns:\n    if X_tr[col].isnull().any():\n        print(col)\n        mean_value = X_tr.loc[X_tr[col] != -np.inf, col].mean()\n        X_tr.loc[X_tr[col] == -np.inf, col] = mean_value\n        X_tr[col] = X_tr[col].fillna(mean_value)\n        means_dict[col] = mean_value","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a154abe0bd5d3e16b216d9375f09166aabba7d36"},"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_tr)\nX_train_scaled = pd.DataFrame(scaler.transform(X_tr), columns=X_tr.columns)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bb0b51806d92a3e4c788e439ac17d5a5685fdfa8"},"cell_type":"markdown","source":"### Reading test dat"},{"metadata":{"trusted":true,"_uuid":"a1b46d0a370c3817a1cc98c022f9c62ebd135331","_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(columns=X_tr.columns, dtype=np.float64, index=submission.index)\nplt.figure(figsize=(22, 16))\n\nfor i, seg_id in enumerate(tqdm_notebook(X_test.index)):\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    \n    x = pd.Series(seg['acoustic_data'].values)\n    X_test.loc[seg_id, 'mean'] = x.mean()\n    X_test.loc[seg_id, 'std'] = x.std()\n    X_test.loc[seg_id, 'max'] = x.max()\n    X_test.loc[seg_id, 'min'] = x.min()\n        \n    X_test.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(x))\n    X_test.loc[seg_id, 'mean_change_rate'] = calc_change_rate(x)\n    X_test.loc[seg_id, 'abs_max'] = np.abs(x).max()\n    X_test.loc[seg_id, 'abs_min'] = np.abs(x).min()\n    \n    X_test.loc[seg_id, 'std_first_50000'] = x[:50000].std()\n    X_test.loc[seg_id, 'std_last_50000'] = x[-50000:].std()\n    X_test.loc[seg_id, 'std_first_10000'] = x[:10000].std()\n    X_test.loc[seg_id, 'std_last_10000'] = x[-10000:].std()\n    \n    X_test.loc[seg_id, 'avg_first_50000'] = x[:50000].mean()\n    X_test.loc[seg_id, 'avg_last_50000'] = x[-50000:].mean()\n    X_test.loc[seg_id, 'avg_first_10000'] = x[:10000].mean()\n    X_test.loc[seg_id, 'avg_last_10000'] = x[-10000:].mean()\n    \n    X_test.loc[seg_id, 'min_first_50000'] = x[:50000].min()\n    X_test.loc[seg_id, 'min_last_50000'] = x[-50000:].min()\n    X_test.loc[seg_id, 'min_first_10000'] = x[:10000].min()\n    X_test.loc[seg_id, 'min_last_10000'] = x[-10000:].min()\n    \n    X_test.loc[seg_id, 'max_first_50000'] = x[:50000].max()\n    X_test.loc[seg_id, 'max_last_50000'] = x[-50000:].max()\n    X_test.loc[seg_id, 'max_first_10000'] = x[:10000].max()\n    X_test.loc[seg_id, 'max_last_10000'] = x[-10000:].max()\n    \n    X_test.loc[seg_id, 'max_to_min'] = x.max() / np.abs(x.min())\n    X_test.loc[seg_id, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n    X_test.loc[seg_id, 'count_big'] = len(x[np.abs(x) > 500])\n    X_test.loc[seg_id, 'sum'] = x.sum()\n    \n    X_test.loc[seg_id, 'mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n    X_test.loc[seg_id, 'mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])\n    X_test.loc[seg_id, 'mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n    X_test.loc[seg_id, 'mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])\n    \n    X_test.loc[seg_id, 'q95'] = np.quantile(x,0.95)\n    X_test.loc[seg_id, 'q99'] = np.quantile(x,0.99)\n    X_test.loc[seg_id, 'q05'] = np.quantile(x,0.05)\n    X_test.loc[seg_id, 'q01'] = np.quantile(x,0.01)\n    \n    X_test.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n    X_test.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n    X_test.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n    X_test.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    X_test.loc[seg_id, 'trend'] = add_trend_feature(x)\n    X_test.loc[seg_id, 'abs_trend'] = add_trend_feature(x, abs_values=True)\n    X_test.loc[seg_id, 'abs_mean'] = np.abs(x).mean()\n    X_test.loc[seg_id, 'abs_std'] = np.abs(x).std()\n    \n    X_test.loc[seg_id, 'mad'] = x.mad()\n    X_test.loc[seg_id, 'kurt'] = x.kurtosis()\n    X_test.loc[seg_id, 'skew'] = x.skew()\n    X_test.loc[seg_id, 'med'] = x.median()\n    \n    X_test.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(x)).mean()\n    X_test.loc[seg_id, 'Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n    X_test.loc[seg_id, 'classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n    X_test.loc[seg_id, 'classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n    X_test.loc[seg_id, 'classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n    X_test.loc[seg_id, 'classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n    X_test.loc[seg_id, 'classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean()\n    X_test.loc[seg_id, 'classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n    X_test.loc[seg_id, 'classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean()\n    X_test.loc[seg_id, 'classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n    X_test.loc[seg_id, 'Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    X_test.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n    X_test.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n    X_test.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n    no_of_std = 3\n    X_test.loc[seg_id, 'MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n    X_test.loc[seg_id,'MA_700MA_BB_high_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X_test.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X_test.loc[seg_id,'MA_700MA_BB_low_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X_test.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X_test.loc[seg_id, 'MA_400MA_std_mean'] = x.rolling(window=400).std().mean()\n    X_test.loc[seg_id,'MA_400MA_BB_high_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X_test.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X_test.loc[seg_id,'MA_400MA_BB_low_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X_test.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X_test.loc[seg_id, 'MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n    X_test.drop('Moving_average_700_mean', axis=1, inplace=True)\n    \n    X_test.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n    X_test.loc[seg_id, 'q999'] = np.quantile(x,0.999)\n    X_test.loc[seg_id, 'q001'] = np.quantile(x,0.001)\n    X_test.loc[seg_id, 'ave10'] = stats.trim_mean(x, 0.1)\n    \n    for windows in [10, 100, 1000]:\n        x_roll_std = x.rolling(windows).std().dropna().values\n        x_roll_mean = x.rolling(windows).mean().dropna().values\n        \n        X_test.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X_test.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X_test.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X_test.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X_test.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X_test.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X_test.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X_test.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X_test.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X_test.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_test.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X_test.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X_test.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X_test.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X_test.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X_test.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X_test.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X_test.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X_test.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X_test.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X_test.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_test.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()\n    \n    if i < 12:\n        plt.subplot(6, 4, i + 1)\n        plt.plot(seg['acoustic_data'])\n        plt.title(seg_id)\n\nfor col in X_test.columns:\n    if X_test[col].isnull().any():\n        X_test.loc[X_test[col] == -np.inf, col] = means_dict[col]\n        X_test[col] = X_test[col].fillna(means_dict[col])\n        \nX_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"759245da1d8ecdcda5f78cf6c866ef9ccb4f3aff"},"cell_type":"markdown","source":"## Building models"},{"metadata":{"trusted":true,"_uuid":"b1e871cec60e445e82260cb6c3d7b553e290417a"},"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b0f5152fc9a9a851838d0d1a9b134f3bcde03d0d","_kg_hide-input":true},"cell_type":"code","source":"def train_model(X=X_train_scaled, X_test=X_test_scaled, y=y_tr, params=None, folds=folds, model_type='lgb', plot_feature_importance=False, model=None):\n\n    oof = np.zeros(len(X))\n    prediction = np.zeros(len(X_test))\n    scores = []\n    feature_importance = pd.DataFrame()\n    for fold_n, (train_index, valid_index) in enumerate(folds.split(X)):\n        print('Fold', fold_n, 'started at', time.ctime())\n        X_train, X_valid = X.iloc[train_index], X.iloc[valid_index]\n        y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]\n        \n        if model_type == 'lgb':\n            model = lgb.LGBMRegressor(**params, n_estimators = 50000, n_jobs = -1)\n            model.fit(X_train, y_train, \n                    eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric='mae',\n                    verbose=10000, early_stopping_rounds=200)\n            \n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test, num_iteration=model.best_iteration_)\n            \n        if model_type == 'xgb':\n            train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)\n            valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)\n\n            watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]\n            model = xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, verbose_eval=500, params=params)\n            y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n            y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n        \n        if model_type == 'sklearn':\n            model = model\n            model.fit(X_train, y_train)\n            \n            y_pred_valid = model.predict(X_valid).reshape(-1,)\n            score = mean_absolute_error(y_valid, y_pred_valid)\n            print(f'Fold {fold_n}. MAE: {score:.4f}.')\n            print('')\n            \n            y_pred = model.predict(X_test).reshape(-1,)\n        \n        if model_type == 'cat':\n            model = CatBoostRegressor(iterations=20000,  eval_metric='MAE', **params)\n            model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)\n\n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test)\n        \n        oof[valid_index] = y_pred_valid.reshape(-1,)\n        scores.append(mean_absolute_error(y_valid, y_pred_valid))\n\n        prediction += y_pred    \n        \n        if model_type == 'lgb':\n            # feature importance\n            fold_importance = pd.DataFrame()\n            fold_importance[\"feature\"] = X.columns\n            fold_importance[\"importance\"] = model.feature_importances_\n            fold_importance[\"fold\"] = fold_n + 1\n            feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n    prediction /= n_fold\n    \n    print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n    \n    if model_type == 'lgb':\n        feature_importance[\"importance\"] /= n_fold\n        if plot_feature_importance:\n            cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n                by=\"importance\", ascending=False)[:50].index\n\n            best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n            plt.figure(figsize=(16, 12));\n            sns.barplot(x=\"importance\", y=\"feature\", data=best_features.sort_values(by=\"importance\", ascending=False));\n            plt.title('LGB Features (avg over folds)');\n        \n            return oof, prediction, feature_importance\n        return oof, prediction\n    \n    else:\n        return oof, prediction","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f368d0dc3a116451650e878e81da909ef4e942ae"},"cell_type":"code","source":"params = {'num_leaves': 128,\n          'min_data_in_leaf': 79,\n          'objective': 'huber',\n          'max_depth': -1,\n          'learning_rate': 0.01,\n          \"boosting\": \"gbdt\",\n          \"bagging_freq\": 5,\n          \"bagging_fraction\": 0.8126672064208567,\n          \"bagging_seed\": 11,\n          \"metric\": 'mae',\n          \"verbosity\": -1,\n          'reg_alpha': 0.1302650970728192,\n          'reg_lambda': 0.3603427518866501\n         }\noof_lgb, prediction_lgb, feature_importance = train_model(params=params, model_type='lgb', plot_feature_importance=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1dbc733abb7ef3315410d987689bbc0c181d974a"},"cell_type":"code","source":"top_cols = list(feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n                by=\"importance\", ascending=False)[:50].index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"389f82c089d500c80644862d7e55acdaf4bae508"},"cell_type":"code","source":"# Taking less columns seriously decreases score.\n# X_train_scaled = X_train_scaled[top_cols]\n# X_test_scaled = X_test_scaled[top_cols]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bd51243ddfa1b86b6f713f64c5005800380afa04"},"cell_type":"code","source":"oof_lgb, prediction_lgb, feature_importance = train_model(X=X_train_scaled, X_test=X_test_scaled, params=params, model_type='lgb', plot_feature_importance=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6f52e5cea8fea8f58c10e888616997fbbfe5c749"},"cell_type":"code","source":"xgb_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}\noof_xgb, prediction_xgb = train_model(X=X_train_scaled, X_test=X_test_scaled, params=xgb_params, model_type='xgb')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"48a35d7fe14b3ab9c4778571928fce9b7c2cd042","scrolled":true},"cell_type":"code","source":"model = NuSVR(gamma='scale', nu=0.9, C=10.0, tol=0.01)\noof_svr, prediction_svr = train_model(X=X_train_scaled, X_test=X_test_scaled, params=None, model_type='sklearn', model=model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a3b4ec5e4a8a8b824a6f1a56106b1e8e7ab6ebc0"},"cell_type":"code","source":"model = NuSVR(gamma='scale', nu=0.7, tol=0.01, C=1.0)\noof_svr1, prediction_svr1 = train_model(X=X_train_scaled, X_test=X_test_scaled, params=None, model_type='sklearn', model=model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f66e18dbce5690fde671799b809e40560991de90"},"cell_type":"code","source":"params = {'loss_function':'MAE'}\noof_cat, prediction_cat = train_model(X=X_train_scaled, X_test=X_test_scaled, params=params, model_type='cat')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5fecb618f872c6991afe86a8af17544ec8b76665"},"cell_type":"code","source":"model = KernelRidge(kernel='rbf', alpha=0.1, gamma=0.01)\noof_r, prediction_r = train_model(X=X_train_scaled, X_test=X_test_scaled, params=None, model_type='sklearn', model=model)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6adbcb019ff54fc86e276fc69392d24e8a6b46ec"},"cell_type":"markdown","source":"### Stacking\nAnd now let's try stacking :) We can use the same function for it."},{"metadata":{"trusted":true,"_uuid":"6a7caafd57c26da85fd62a73a0a98275c06e15fa"},"cell_type":"code","source":"train_stack = np.vstack([oof_lgb, oof_xgb, oof_svr, oof_svr1, oof_r, oof_cat]).transpose()\ntrain_stack = pd.DataFrame(train_stack, columns = ['lgb', 'xgb', 'svr', 'svr1', 'r', 'cat'])\ntest_stack = np.vstack([prediction_lgb, prediction_xgb, prediction_svr, prediction_svr1, prediction_r, prediction_cat]).transpose()\ntest_stack = pd.DataFrame(test_stack)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"ec445fbb4ba84dbebbb73a62cb7d97294de37776"},"cell_type":"code","source":"oof_lgb_stack, prediction_lgb_stack, feature_importance = train_model(X=train_stack, X_test=test_stack, params=params, model_type='lgb', plot_feature_importance=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7b997000c2a2aba8b51dad51b66d633afbf9adb9"},"cell_type":"markdown","source":"Now let's see how do our models perform"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"29557b1a9c6619ddf71b11cd47ad5c6b720d51a0"},"cell_type":"code","source":"plt.figure(figsize=(18, 8))\nplt.subplot(2, 3, 1)\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot(oof_lgb, color='b', label='lgb')\nplt.legend(loc=(1, 0.5));\nplt.title('lgb');\nplt.subplot(2, 3, 2)\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot(oof_xgb, color='teal', label='xgb')\nplt.legend(loc=(1, 0.5));\nplt.title('xgb');\nplt.subplot(2, 3, 3)\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot(oof_svr, color='red', label='svr')\nplt.legend(loc=(1, 0.5));\nplt.title('svr');\nplt.subplot(2, 3, 4)\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot(oof_cat, color='b', label='cat')\nplt.legend(loc=(1, 0.5));\nplt.title('cat');\nplt.subplot(2, 3, 5)\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot(oof_lgb_stack, color='gold', label='stack')\nplt.legend(loc=(1, 0.5));\nplt.title('blend');\nplt.legend(loc=(1, 0.5));\nplt.suptitle('Predictions vs actual');\nplt.subplot(2, 3, 6)\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot((oof_lgb + oof_xgb + oof_svr + oof_svr1 + oof_r + oof_cat) / 6, color='gold', label='blend')\nplt.legend(loc=(1, 0.5));\nplt.title('blend');\nplt.legend(loc=(1, 0.5));\nplt.suptitle('Predictions vs actual');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2916c7f8c188b5455bfd55cbbfba1c72ca1ef86a"},"cell_type":"markdown","source":"It turned out that stacking is much worse than blending on LB."},{"metadata":{"trusted":true,"_uuid":"c5d60a48d8c9071236815b0931111e86d20730fa"},"cell_type":"code","source":"submission['time_to_failure'] = (prediction_lgb + prediction_xgb + prediction_svr + prediction_svr1 + prediction_cat + prediction_r) / 6\n# submission['time_to_failure'] = prediction_lgb_stack\nprint(submission.head())\nsubmission.to_csv('submission.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Saving features to csv, so that they can be easily used."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_tr.to_csv('train_features.csv', index=False)\nX_test.to_csv('test_features.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# From https://www.kaggle.com/rtatman/download-a-csv-file-from-a-kernel\nfrom IPython.display import HTML\nimport base64\n\n# function that takes in a dataframe and creates a text link to  \n# download it (will only work for files < 2MB or so)\ndef create_download_link(df, title = \"Download CSV file\", filename = \"submission.csv\"):  \n    csv = df.to_csv()\n    b64 = base64.b64encode(csv.encode())\n    payload = b64.decode()\n    html = '<a download=\"{filename}\" href=\"data:text/csv;base64,{payload}\" target=\"_blank\">{title}</a>'\n    html = html.format(payload=payload,title=title,filename=filename)\n    return HTML(html)\n\n# create a link to download the dataframe\ncreate_download_link(Earthquakes FE.More features and samples 8d9773)\n\n\n","execution_count":null,"outputs":[]}],"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}