{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### Important notice\n\nDue to some weird bug I can't answer the comments on this kernel - my answers aren't posted. And I can't upvote comments on this kernel. So, I'm sorry, but I can't answer :(\nBut it seems I can answer in other kernels, so I can reply, if you post a comment to my second kernel.\n\n## 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.\n\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 kernel is dedicated to exploration of LANL Earthquake Prediction Challenge.\n\nIn my second [kernel](https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples) I pay more attention to feature engineering and data generation.\n\n![](https://images.spot.im/v1/production/vbwfwtqv6krdqcmxu2k9)","metadata":{"_uuid":"041eccd08aafc3e46f18b1b7a1b80a2b28a0003f"}},{"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\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\nfrom sklearn.linear_model import Ridge, RidgeCV\nimport gc\nfrom catboost import CatBoostRegressor\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ntrain = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"_uuid":"d1d33084c48168284d75d0da861009789cb081d4","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.shape","metadata":{"_uuid":"7ad2361064a108547c6c9e7277dd920583cbe165","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have 629 million rows! Huge data, so let's plot a sample of it.","metadata":{"_uuid":"a2c436d0c0abbefac651799f63ad52554a5bcbd4"}},{"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","metadata":{"_uuid":"95b8d85ba4348b7ff4ce389d3a44cebf6e549e6a","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. First 2% of data\")\nplt.plot(train['acoustic_data'].values[:12582910], color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train['time_to_failure'].values[:12582910], color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)","metadata":{"_kg_hide-input":true,"_uuid":"4149effc6e32ee4a5c4bab817af080c3696a0288","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"On the first plot you can see 2% of all data (all data skipping each 50 rows). We can see that usually acoustic data shows huge fluctuations just before the failure and the nature of data is cyclical.\n\nOn the second plot we see first 2% of the data. It seems that at first the signal has huge fluctuations for a short time, then it lowers and after some time the earthquake occurs. I think it will be quite diffucult to distinguish target values properly","metadata":{"_uuid":"584a3f891dbfee0b82ab8c4dfdab5fb35c88d07b"}},{"cell_type":"markdown","source":"### Feature engineering\n\nLet's create some new features.\n\nWhy 150000? Test segments are 150000 each. At first I create features similar to baseline kernel, but with more aggregations.","metadata":{"_uuid":"c4972c4a8bc047aa05d909f0b447340f3c5353b3"}},{"cell_type":"code","source":"# Create a training file with simple derived features\nrows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\nX_tr = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       columns=['ave', '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'])\ny_tr = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       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_max = np.abs(train['acoustic_data']).sum()\n\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    \n    y_tr.loc[segment, 'time_to_failure'] = y\n    X_tr.loc[segment, 'ave'] = 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    \n    X_tr.loc[segment, 'av_change_abs'] = np.mean(np.diff(x))\n    X_tr.loc[segment, 'av_change_rate'] = np.mean(np.nonzero((np.diff(x) / x[:-1]))[0])\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()","metadata":{"_uuid":"2b369f9462919c0c6e29f8a5b6183cb3c57f057a","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'{X_tr.shape[0]} samples in new train data.')","metadata":{"_uuid":"c03e361f5390471f0281367b9520933c7c999223","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see all new features","metadata":{"_uuid":"d4539c37e6282084d760678b1d70a2b76e7fe3d1"}},{"cell_type":"code","source":"plt.figure(figsize=(26, 24))\nfor i, col in enumerate(X_tr.columns):\n    plt.subplot(7, 4, i + 1)\n    plt.plot(X_tr[col], color='blue')\n    plt.title(col)\n    ax1.set_ylabel(col, color='b')\n    # plt.legend([col])\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)","metadata":{"_uuid":"b0f0a3d919074b0630a655dc560e7fc2287e8017","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_tr)\nX_train_scaled = pd.DataFrame(scaler.transform(X_tr), columns=X_tr.columns)","metadata":{"_uuid":"33d7dc2509af6eb317111b9b6b0021c0af3c75e0","trusted":true},"execution_count":null,"outputs":[]},{"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 = seg['acoustic_data'].values\n    X_test.loc[seg_id, 'ave'] = 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, 'av_change_abs'] = np.mean(np.diff(x))\n    X_test.loc[seg_id, 'av_change_rate'] = np.mean(np.nonzero((np.diff(x) / x[:-1]))[0])\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    if i < 24:\n        plt.subplot(7, 4, i + 1)\n        plt.plot(seg['acoustic_data'])\n        plt.title(seg_id)\n    \nX_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns)","metadata":{"_uuid":"8a8fefba8af693620aa3598ac3de8f2ad521bbf4","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","metadata":{"_uuid":"9172a3903b178e151d8aefbc68e01426a39e9e2e","trusted":true},"execution_count":null,"outputs":[]},{"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        if model_type == 'xgb':\n            train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X_tr.columns)\n            valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X_tr.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_tr.columns), ntree_limit=model.best_ntree_limit)\n            y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X_tr.columns), ntree_limit=model.best_ntree_limit)\n            \n        if model_type == 'rcv':\n            model = RidgeCV(alphas=(0.01, 0.1, 1.0, 10.0, 100.0), scoring='neg_mean_absolute_error', cv=3)\n            model.fit(X_train, y_train)\n            print(model.alpha_)\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 == '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","metadata":{"_uuid":"815e13e6c0a78fc3fa461d4c823c97b185838a7e","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"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)","metadata":{"_uuid":"c532c79cb6f3417c3a9eec579979aaea91a0aa10","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb_params = {'eta': 0.03, 'max_depth': 10, 'subsample': 0.85, #'colsample_bytree': 0.8, \n          'objective': 'reg:linear', 'eval_metric': 'mae', 'silent': True, 'nthread': 4}\noof_xgb, prediction_xgb = train_model(params=xgb_params, model_type='xgb')","metadata":{"_uuid":"2d01a97d68c4031b3b2d11974a17b2e7053efc88","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = NuSVR(gamma='scale', nu=0.75, C=10.0)\noof_svr, prediction_svr = train_model(params=None, model_type='sklearn', model=model)","metadata":{"scrolled":true,"_uuid":"076994056cdd612e6628f980a31415547649019a","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Interestingly, if we compare the results of training on this data with results of training on smaller dataset (as in official baseline), then we can see that MAE for XGB model became less than MAE of LGB model.\n\nNow let's see how do our models perform","metadata":{"_uuid":"4e3847381aebc9d77e1cfa12b5d16d51e09d46a0"}},{"cell_type":"code","source":"plt.figure(figsize=(16, 8))\nplt.plot(y_tr, color='g', label='y_train')\nplt.plot(oof_lgb, color='b', label='lgb')\nplt.plot(oof_xgb, color='teal', label='xgb')\nplt.plot(oof_svr, color='red', label='svr')\nplt.plot((oof_lgb + oof_xgb + oof_svr) / 3, color='gold', label='blend')\nplt.legend();\nplt.title('Predictions vs actual');","metadata":{"_uuid":"93248680af51c40093dcdc24cc2f8f1d4d70de44","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that models can't predict high values well, but additional data is predicted much better.","metadata":{"_uuid":"b917ae3354e6b1576f1cea218a94b1533bfde5dd"}},{"cell_type":"code","source":"prediction_lgb[:10], prediction_xgb[:10], prediction_svr[:10]","metadata":{"_uuid":"aed2c67cf5b7596b785bb85b0943ee9fef6e8187","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission['time_to_failure'] = (prediction_lgb + prediction_xgb + prediction_svr) / 3\nprint(submission.head())\nsubmission.to_csv('submission.csv')","metadata":{"_uuid":"acf417f3540ca69267b1088e199e4ae358b36439","trusted":true},"execution_count":null,"outputs":[]}]}