{"cells":[{"metadata":{"trusted":true,"_uuid":"f02146926f9478e8a005abf19c743f1be3de39f6"},"cell_type":"code","source":"from fastai.tabular import *","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"30c1f20e07a15cd91f3056208f3c8d3eb238ebf3"},"cell_type":"code","source":"path = Path(\"../input\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"890d3950a6132214736c4074922521737aca80b5"},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(path/'train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9796412b9d4fc2361ffb229e126071c24756fc2e"},"cell_type":"markdown","source":"# Preprocessing data.\nI used https://www.kaggle.com/gpreda/lanl-earthquake-eda-and-prediction"},{"metadata":{"trusted":true,"_uuid":"89261a53cabd2c2069f627ff6c4d47746c259424"},"cell_type":"code","source":"rows = 150000\nsegments = int(np.floor(train_df.shape[0] / rows))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5cdb7f28e08186acd27e57849c69393620954fa4"},"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'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c019e349542faf34c2993789d205cec326b1130c"},"cell_type":"code","source":"import warnings\nfrom scipy import stats\nfrom scipy.signal import hann\nimport matplotlib.pyplot as plt\nfrom scipy.signal import hilbert\nfrom scipy.signal import convolve\nfrom sklearn.svm import NuSVR, SVR\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 tqdm import tqdm\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0365dffe580dce46c07c8f361283c59076b02739"},"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":{"trusted":true,"_uuid":"8d8c8db2a98dac84c0046d33bb675cfa328308e9"},"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    \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, '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    \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    \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, '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":"fc5a25685c5af8a12d4478dbe75e6a9bad0c19d1"},"cell_type":"code","source":"# iterate over all segments\nfor seg_id in tqdm(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":{"trusted":true,"_uuid":"c8b2fa785da89673e3a76b6eab4c0032789e7b43"},"cell_type":"code","source":"submission = pd.read_csv(path/'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":{"trusted":true,"_uuid":"a391483ba9565d6e789c4052f088ea12bcf27093"},"cell_type":"code","source":"for seg_id in tqdm(test_X.index):\n    seg = pd.read_csv(path/f'test/{seg_id}.csv')\n    create_features(seg_id, seg, test_X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4ec1f654a82d5403d759393fa25d01f4c50625a0"},"cell_type":"code","source":"test_X.reset_index(inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b4d47c592cd5cc0c90b96a1710e92f51d07bdfd9"},"cell_type":"code","source":"df = pd.concat([train_X, train_y], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"82b1d711b7b2bc5b9d1e59bd394572c94747ca9c"},"cell_type":"code","source":"df.to_feather(\"df\")\ntest_X.to_hdf('df_test.h5', key=\"df_test\")\ndf = pd.read_feather(\"df\")\ntest_df = pd.read_hdf('df_test.h5', key=\"df_test\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"871ecc3b4af29928cceb2daa125e9c4f5d8e0034"},"cell_type":"code","source":"df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b214aadd4763eccdca19fc9bb03ad9a22ed458c9"},"cell_type":"code","source":"dep_var = 'time_to_failure'\ncont_names = None\ncat_names = None\n\nprocs = [FillMissing, Categorify, Normalize]\n# procs = [FillMissing, Categorify]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6655de447bc4e83ccde7fa866d67af80349ae06e"},"cell_type":"code","source":"data = TabularDataBunch.from_df(path, df, dep_var, valid_idx=range(0,800), procs=procs, test_df=test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e5d9cabc42831249a3f8ff69f0bdee1fdb015ade"},"cell_type":"code","source":"data.show_batch(rows=10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7256e16869575e415a4a8dd36b8157d901cf2066"},"cell_type":"code","source":"learn = tabular_learner(data, layers=[200,100], model_dir=Path(\"/\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e70251a46e1c3ac27cd21bcdccb2a380f091ddf6"},"cell_type":"code","source":"learn.lr_find(); learn.recorder.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a5f6ab60d06f00f23ff10662270367bcf833af68"},"cell_type":"code","source":"lr = 1e-02","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"681b27183682e2f427507afa00b4d203e21ab508"},"cell_type":"code","source":"learn.fit_one_cycle(5, slice(lr))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e7c0908a4c8d13b4b3bd84b0a1dbfa904c53919b"},"cell_type":"code","source":"learn.save('stage-1')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"13346f6333929353a21e2ad32e70bdfd77b946a8"},"cell_type":"code","source":"learn.load('stage-1');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"44b76938db2044c5b04fc07f3781c04dfd258b95"},"cell_type":"code","source":"learn.unfreeze()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f23b8b61ca5fe9b6404cbefab85ff50f6e1765a2"},"cell_type":"code","source":"learn.lr_find(); learn.recorder.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d564eb83b9720b21f0b760742b4b250605c849e9"},"cell_type":"code","source":"learn.fit_one_cycle(10, slice(2e-05))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"905c0162f2783bcfa41c65bd2c650fc411581d44"},"cell_type":"code","source":"preds, _ = learn.get_preds(ds_type=DatasetType.Test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c86df9491b846d161e261be4e8eea8dd3094e666"},"cell_type":"code","source":"subm_df = pd.DataFrame(test_df[\"seg_id\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"568869e012a92d9999924efd1dcc18ae0ca8e0b0"},"cell_type":"code","source":"subm_df[\"time_to_failure\"] = preds.numpy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fbca742892d1e1c85697cc5736f8cf4ae28b8b88"},"cell_type":"code","source":"subm_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"83c396a4704275cb04a5070955d32f60be6c7bd1"},"cell_type":"code","source":"subm_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"06dfb1099bcf2c576b62c0e0bcbb6fc513757aba"},"cell_type":"code","source":"!kaggle competitions submit LANL-Earthquake-Prediction -f {'submission.csv'} -m \"More features\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9b428235ad4901c142966a4d1ede10ea9fe0ff9e"},"cell_type":"code","source":"","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}