{"cells":[{"metadata":{"_uuid":"8e49039a682501b5ffaa1f4459accab7912a0d37"},"cell_type":"markdown","source":"![](https://pp.userapi.com/c848636/v848636381/10387a/82TkN23uVpQ.jpg)"},{"metadata":{"_uuid":"e40474c5aeb9bea8604b63ea1af9c1109febd466"},"cell_type":"markdown","source":"# Intro\nThis kernel is dedicated to exploration of LANL Earthquake Prediction Challenge. \nThis kernel is different in that I suggest trying out different methods and functions that are used to process signals for feature extraction:\n* [Hilbert transform](http:/en.wikipedia.org/wiki/Analytic_signal) \n* Smooth a pulse using a [Hann](http://en.wikipedia.org/wiki/Hann_function) window \n* Use trigger [classic STA/LTA](http://docs.obspy.org/tutorial/code_snippets/trigger_tutorial.html#available-methods)\n\nThank [ Andrew Lukyanenko](http://www.kaggle.com/artgor) for his [kernel](http://www.kaggle.com/artgor/seismic-data-eda-and-baseline/output)  - my decision is based on his.\n\nThank [Vishy](http://www.kaggle.com/viswanathravindran) for his [discuss](http://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/77267) and links.\n\nYou can view and discuss these features in my [kernel](https://www.kaggle.com/nikitagribov/analysis-function-for-signal-data) or [discuss](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/77267#455024) [Vishy](http://www.kaggle.com/viswanathravindran)."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"import numpy as np\nimport pandas as pd \nfrom tqdm import tqdm_notebook\n\nimport plotly.offline as py\npy.init_notebook_mode(connected=True)\nimport plotly.graph_objs as go\nfrom plotly import tools\nimport plotly.figure_factory as ff\n\nfrom scipy.signal import hilbert, hann, convolve\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm_notebook\nfrom sklearn.preprocessing import StandardScaler\npd.options.display.precision = 15\n\nimport lightgbm as lgb\nimport time\nimport datetime\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\nimport gc\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"%%time\ntrain = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.float32, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9e739445936bee2fce7ec27b44e6d59bf64fb279"},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ef5da6b59666d4d0737c423d37f2dd8b7562075c"},"cell_type":"markdown","source":"# New Features example"},{"metadata":{"trusted":true,"_uuid":"c2fe9f2b49d2e02775628e85686d7913323a8708"},"cell_type":"code","source":"def classic_sta_lta_py(a, nsta, nlta):\n    \"\"\"\n    Computes the standard STA/LTA from a given input array a. The length of\n    the STA is given by nsta in samples, respectively is the length of the\n    LTA given by nlta in samples. Written in Python.\n    .. note::\n        There exists a faster version of this trigger wrapped in C\n        called :func:`~obspy.signal.trigger.classic_sta_lta` in this module!\n    :type a: NumPy :class:`~numpy.ndarray`\n    :param a: Seismic Trace\n    :type nsta: int\n    :param nsta: Length of short time average window in samples\n    :type nlta: int\n    :param nlta: Length of long time average window in samples\n    :rtype: NumPy :class:`~numpy.ndarray`\n    :return: Characteristic function of classic STA/LTA\n    \"\"\"\n    # The cumulative sum can be exploited to calculate a moving average (the\n    # cumsum function is quite efficient)\n    sta = np.cumsum(a ** 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[nsta:] = sta[nsta:] - sta[:-nsta]\n    sta /= nsta\n    lta[nlta:] = lta[nlta:] - lta[:-nlta]\n    lta /= nlta\n\n    # Pad zeros\n    sta[:nlta - 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f23d017cb1fcd826c736c4422ad578331411ab6a"},"cell_type":"code","source":"per = 0.000005","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e0958ef8cb3e9abaaa23fc24e0e257057e8102f8"},"cell_type":"code","source":"#Calculate Hilbert transform\nsignal = train.acoustic_data[:int(len(train)*per)]\nanalytic_signal = hilbert(signal)\namplitude_envelope = np.abs(analytic_signal)\n\n#Calculate Hann func\nwin = hann(50)\nfiltered = convolve(signal, win, mode='same') / sum(win)\n\n#Calculate STA/LTA\nsta_lta = classic_sta_lta_py(signal, 50, 1000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9622001199494d19772ad23a467a50bf12ffafa2","_kg_hide-input":false},"cell_type":"code","source":"trace0 = go.Scatter(\n    y = signal,\n    name = 'signal'\n)\n\ntrace1 = go.Scatter(\n    y = amplitude_envelope,\n    name = 'amplitude_envelope'\n)\n\n\ntrace3 = go.Scatter(\n    y = filtered,\n    name= 'filtered'\n) \n\ntrace4 = go.Scatter(\n    y = sta_lta,\n    name= 'sta_lta'\n) \n\n\ntrace_time = go.Scatter(\n    y = train.time_to_failure[:int(len(train)*per)],\n)\n\ndata = [trace0, trace1, trace3,trace4]\n\nlayout = go.Layout(\n    title = \"Part acoustic_data\"\n)\n\nfig = go.Figure(data=data,layout=layout)\npy.iplot(fig, filename = \"Part acoustic_data\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e61b68b039700c06b936a1a3ca9638fb484d5c6c"},"cell_type":"code","source":"col_feat = ['ave','med', 'std', 'max', 'min',\n                               'av_change_abs', 'av_change_rate', 'abs_max', 'abs_min',\n                               'std_first_50000', 'std_last_50000', 'std_first_10000', 'std_last_10000',\n                               'avg_first_50000', 'avg_last_50000', 'avg_first_10000', 'avg_last_10000',\n                               'min_first_50000', 'min_last_50000', 'min_first_10000', 'min_last_10000',\n                               'max_first_50000', 'max_last_50000', 'max_first_10000', 'max_last_10000',\n                               ]\nname_lines = ['','amp_env','filt','sta_lta']\ncolumns = []\nfor cf in col_feat:\n    for name in name_lines:\n        columns.append(cf+'_'+name)\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5ff9e421dba8db91df738e9102263019c6d87f0e"},"cell_type":"code","source":"len(columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bf279eae37e8b41fda3de7c67e31ed1b7bd46ccd"},"cell_type":"code","source":"rows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\ny_tr = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       columns=['time_to_failure'])\n\nX_tr = pd.DataFrame()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4240fa4a96a403dbc623dc81306955f258d6f00f"},"cell_type":"code","source":"\ndef get_feature(x,segment, name_lines):\n    X_tr1 = pd.DataFrame(dtype=np.float64)\n    X_tr1.loc[segment, 'ave'+'_'+name_lines] = x.mean()\n    X_tr1.loc[segment, 'ave'+'_'+name_lines] = np.median(x)\n    X_tr1.loc[segment, 'std'+'_'+name_lines] = x.std()\n    X_tr1.loc[segment, 'max'+'_'+name_lines] = x.max()\n    X_tr1.loc[segment, 'min'+'_'+name_lines] = x.min()\n    \n    X_tr1.loc[segment, 'av_change_abs'+'_'+name_lines] = np.mean(np.diff(x))\n    X_tr1.loc[segment, 'av_change_rate'+'_'+name_lines] = np.mean(np.nonzero((np.diff(x) / x[:-1]))[0])\n    X_tr1.loc[segment, 'abs_max'+'_'+name_lines] = np.abs(x).max()\n    X_tr1.loc[segment, 'abs_min'+'_'+name_lines] = np.abs(x).min()\n    \n    X_tr1.loc[segment, 'std_first_50000'+'_'+name_lines] = x[:50000].std()\n    X_tr1.loc[segment, 'std_last_50000'+'_'+name_lines] = x[50000:].std()\n    X_tr1.loc[segment, 'std_first_10000'+'_'+name_lines] = x[:10000].std()\n    X_tr1.loc[segment, 'std_last_10000'+'_'+name_lines] = x[10000:].std()\n    \n    X_tr1.loc[segment, 'avg_first_50000'+'_'+name_lines] = x[:50000].mean()\n    X_tr1.loc[segment, 'avg_last_50000'+'_'+name_lines] = x[50000:].mean()\n    X_tr1.loc[segment, 'avg_first_10000'+'_'+name_lines] = x[:10000].mean()\n    X_tr1.loc[segment, 'avg_last_10000'+'_'+name_lines] = x[10000:].mean()\n    \n    X_tr1.loc[segment, 'min_first_50000'+'_'+name_lines] = x[:50000].min()\n    X_tr1.loc[segment, 'min_last_50000'+'_'+name_lines] = x[50000:].min()\n    X_tr1.loc[segment, 'min_first_10000'+'_'+name_lines] = x[:10000].min()\n    X_tr1.loc[segment, 'min_last_10000'+'_'+name_lines] = x[10000:].min()\n    \n    X_tr1.loc[segment, 'max_first_50000'+'_'+name_lines] = x[:50000].max()\n    X_tr1.loc[segment, 'max_last_50000'+'_'+name_lines] = x[50000:].max()\n    X_tr1.loc[segment, 'max_first_10000'+'_'+name_lines] = x[:10000].max()\n    X_tr1.loc[segment, 'max_last_10000'+'_'+name_lines] = x[10000:].max()\n    return X_tr1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d04f6dcfc7caceafb556dc1612cbd0d71e262cb2"},"cell_type":"code","source":"np.random.seed(42)\nfor segment in tqdm_notebook(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n    x = seg['acoustic_data'].values\n    y = seg['time_to_failure'].values[-1]\n    y_tr.loc[segment, 'time_to_failure'] = y\n    \n    x_ac_signal = get_feature(x,segment, 'data')\n    \n    analytic_signal = hilbert(x)\n    amplitude_envelope = np.abs(analytic_signal)\n    x_amp_env = get_feature(amplitude_envelope,segment, 'amp_env')\n    \n    win = hann(100)\n    filtered = convolve(x, win, mode='same') / sum(win)\n    x_filt_hann = get_feature(filtered,segment, 'filt')\n    \n    sta_lta = classic_sta_lta_py(x, 1000, 5000)\n    x_sta_lta = get_feature(sta_lta,segment, 'sta_lta')\n    \n    df_loc = pd.concat([x_ac_signal, x_amp_env, x_filt_hann, x_sta_lta], axis = 1)\n    X_tr = X_tr.append(df_loc)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1e0e633b8b6c0dcb357888514515d23aaa506370"},"cell_type":"code","source":"segments = 10000\n\ny_tr_more = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       columns=['time_to_failure'])\nX_tr_more = pd.DataFrame()\nnp.random.seed(42)\nfor segment in tqdm_notebook(range(segments)):\n    ind = np.random.randint(0, train.shape[0]-150001)\n    seg = train.iloc[ind:ind+rows]\n    x = seg['acoustic_data'].values\n    y = seg['time_to_failure'].values[-1]\n    \n    y_tr_more.loc[segment, 'time_to_failure'] = y\n    x_ac_signal = get_feature(x,segment, 'data')\n    \n    analytic_signal = hilbert(x)\n    amplitude_envelope = np.abs(analytic_signal)\n    x_amp_env = get_feature(amplitude_envelope,segment, 'amp_env')\n    \n    win = hann(100)\n    filtered = convolve(x, win, mode='same') / sum(win)\n    x_filt_hann = get_feature(filtered,segment, 'filt')\n    \n    sta_lta = classic_sta_lta_py(x, 1000, 5000)\n    x_sta_lta = get_feature(sta_lta,segment, 'sta_lta')\n    \n    df_loc = pd.concat([x_ac_signal, x_amp_env, x_filt_hann, x_sta_lta], axis = 1)\n    X_tr_more = X_tr_more.append(df_loc)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8d9303fc98f6c3a54ea347d97c2d64025e9a52ea"},"cell_type":"code","source":"X_tr = X_tr.append(X_tr_more)\ny_tr = y_tr.append(y_tr_more)\nprint(f'{X_tr.shape[0]} samples in new train data now.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3fddbdaad2cc3ef8de920f270bf4bb202a7052e3"},"cell_type":"code","source":"X_tr.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7a0fff93c666a89d863ff0203fb8ecbe62277142"},"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":{"trusted":true,"_uuid":"0a5bf0079f24a690c8a8b351ca87eb2f90fb23e9"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0703f9d39a2468002148b4be15c1c1bdf5fec199"},"cell_type":"code","source":"for i, seg_id in enumerate(tqdm_notebook(submission.index)):\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    \n    x = seg['acoustic_data'].values\n    \n    x_ac_signal = get_feature(x,seg_id, 'data')\n    \n    analytic_signal = hilbert(x)\n    amplitude_envelope = np.abs(analytic_signal)\n    x_amp_env = get_feature(amplitude_envelope,seg_id, 'amp_env')\n    \n    win = hann(100)\n    filtered = convolve(x, win, mode='same') / sum(win)\n    x_filt_hann = get_feature(filtered,seg_id, 'filt')\n    \n    sta_lta = classic_sta_lta_py(x, 1000, 5000)\n    x_sta_lta = get_feature(sta_lta,seg_id, 'sta_lta')\n    \n    df_loc = pd.concat([x_ac_signal, x_amp_env, x_filt_hann, x_sta_lta], axis = 1)\n    X_test = X_test.append(df_loc)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"129b13cc3e3dff173edd3ff3cd2036cca3690012"},"cell_type":"code","source":"X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7375e626331de1fa5e3ca960c501af53095d85e6"},"cell_type":"code","source":"len(X_test_scaled)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1c207a5ae53ecdb62fde971b1167ac8e392f9b34"},"cell_type":"code","source":"len(X_train_scaled)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cd3e78a90a0067a3cb83f2f03002eb54a0f33827"},"cell_type":"code","source":"X_train_scaled = X_train_scaled.fillna(0)\nX_test_scaled = X_test_scaled.fillna(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"df647dbb9a4d174b10c61f485c5a14815a9ee1ef"},"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":"b53a895f1080dfe53b8e4452e4d85e8a2f0900da"},"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 = 20000, nthread = 4, 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        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":"cdb99f45594345290fb66745e60d48062772f9a4"},"cell_type":"code","source":"params = {'num_leaves': 54,\n         'min_data_in_leaf': 79,\n         'objective': 'huber',\n         'max_depth': -1,\n         'learning_rate': 0.01,\n         \"boosting\": \"gbdt\",\n         # \"feature_fraction\": 0.8354507676881442,\n         \"bagging_freq\": 3,\n         \"bagging_fraction\": 0.8126672064208567,\n         \"bagging_seed\": 11,\n         \"metric\": 'mae',\n         \"verbosity\": -1,\n         'reg_alpha': 1.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":"ffedf1ee55702637ac1cbf818c6152a05c4e0aac"},"cell_type":"code","source":"len(prediction_lgb)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9f7ddf4237383ad25675d3dd7edac77a39fc60f4"},"cell_type":"code","source":"submission['time_to_failure'] = prediction_lgb\nprint(submission.head())\nsubmission.to_csv('submission.csv')","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}