{"cells":[{"metadata":{"_uuid":"1bc99ff52d68285950e6b52535b03351449f5eef"},"cell_type":"markdown","source":"---\n\n<h1><span style=\"color:red;\"><strong>LANL Earthquake Prediction</strong></span></h1>\n\n---\n\n<h3><span style=\"color:Blue;\"><strong><i>Can you predict upcoming laboratory earthquakes?</i></strong></span></h3>\n\n---\n\n> ## **_Objective_**\n> In this competition, you will address when the earthquake will take place. Specifically, you’ll predict the time remaining before laboratory earthquakes occur from real-time seismic data.\n> ## **_Solution thought by me_**\n> _In this kernel, I tried to apply lightgbm with initial parameter and kfold validation and also try to apply xgboost and other ensemble model and stacking is also applied._\n\n---\n> ## **_Outline_**\n* [**1.Load library**](#1.Load-library)\n* [**2.Read Data**](#2.Read-Data)\n* [**3.Feature Engineering**](#3.Feature-Engineering)\n* [**4.Data transformation**](#4.Data-transformation)\n* [**5.Test Data**](#5.Test-Data)\n* [**6.Model Training**](#6.Model-Training)\n    * [**1. Lightgbm**](#1.-Lightgbm)\n    * [**2. XGboost**](#2.-XGboost)\n* [**7.Stacking**](#7.Stacking)\n> * [**8.Final Prediction**](#8.Final-Prediction)\n---"},{"metadata":{"_uuid":"8d519229ccbe835bd270f584bfddd280271efa41"},"cell_type":"markdown","source":"## **1.Load library**"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm_notebook as tqdm\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import NuSVR\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.metrics import mean_absolute_error, make_scorer\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import StratifiedKFold, KFold\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.linear_model import BayesianRidge\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ac832eb60aa3be659710024160d992360896913d"},"cell_type":"markdown","source":"## **2.Read Data**"},{"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.int16, 'time_to_failure': np.float64})\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b0107b91790b46fb42d0093a1aa8a22f61ccc905"},"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]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4e37e5f6d74968b3d96ce7ade89fac4cf0df587b"},"cell_type":"markdown","source":"## **3.Feature Engineering**"},{"metadata":{"trusted":true,"_uuid":"e04936d09afde6091a302337f4947f8116c568a5"},"cell_type":"code","source":"rows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\nX_train = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       columns=['ave', 'std', 'max', 'min','q95','q99', 'q05','q01',\n                                'abs_max', 'abs_mean', 'abs_std', 'trend', 'abs_trend'])\ny_train = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       columns=['time_to_failure'])\n\nfor segment in tqdm(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_train.loc[segment, 'time_to_failure'] = y\n    \n    X_train.loc[segment, 'ave'] = x.mean()\n    X_train.loc[segment, 'std'] = x.std()\n    X_train.loc[segment, 'max'] = x.max()\n    X_train.loc[segment, 'min'] = x.min()\n    X_train.loc[segment, 'q95'] = np.quantile(x,0.95)\n    X_train.loc[segment, 'q99'] = np.quantile(x,0.99)\n    X_train.loc[segment, 'q05'] = np.quantile(x,0.05)\n    X_train.loc[segment, 'q01'] = np.quantile(x,0.01)\n    \n    X_train.loc[segment, 'abs_max'] = np.abs(x).max()\n    X_train.loc[segment, 'abs_mean'] = np.abs(x).mean()\n    X_train.loc[segment, 'abs_std'] = np.abs(x).std()\n    X_train.loc[segment, 'trend'] = add_trend_feature(x)\n    X_train.loc[segment, 'abs_trend'] = add_trend_feature(x, abs_values=True)\n    \nX_train.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9e7b32d237b4fd402e72ee691d3f36e8b06ed4cd"},"cell_type":"markdown","source":"## **4.Data transformation**"},{"metadata":{"trusted":true,"_uuid":"f3bb36700ff5b557f559029a42cff0e9f9cba557"},"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_train)\nX_train_scaled = scaler.transform(X_train)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"26bf2e5c080e51900681a2c8bfdc4a78325ac66a"},"cell_type":"markdown","source":"## **5.Test Data**"},{"metadata":{"trusted":true,"_uuid":"523e73c6f68c40e5694dbb796cb0da5ccf117ab4"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(columns=X_train.columns, dtype=np.float64, index=submission.index)\nfor seg_id in tqdm(X_test.index):\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    \n    x = seg['acoustic_data'].values\n    \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    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_max'] = np.abs(x).max()\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    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\nX_test_scaled = scaler.transform(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2e7b3abef5702bb48e852d6c5231c5cfddff7eb4"},"cell_type":"code","source":"X_train_scaled = pd.DataFrame(X_train_scaled,columns=X_train.columns)\nX_train_scaled.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"09325f7b837e9b34568fa3e68c9204f9456cd2e5"},"cell_type":"code","source":"X_test_scaled = pd.DataFrame(X_test_scaled,columns=X_test.columns)\nX_test_scaled.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a4e4d8d1f06a9d7ad8f3b2c2a9ae1d9a7f032677"},"cell_type":"markdown","source":"## **6.Model Training**\n---\n## **1. Lightgbm**"},{"metadata":{"trusted":true,"_uuid":"0b091a4082ec1061069c19a245d08c9378851f1d"},"cell_type":"code","source":"param = {'num_leaves': 31,\n         'min_data_in_leaf': 32, \n         'objective':'regression',\n         'max_depth': -1,\n         'learning_rate': 0.001,\n         \"min_child_samples\": 20,\n         \"boosting\": \"gbdt\",\n         \"feature_fraction\": 0.9,\n         \"bagging_freq\": 1,\n         \"bagging_fraction\": 0.9 ,\n         \"bagging_seed\": 11,\n         \"metric\": 'rmse',\n         \"lambda_l1\": 0.1,\n         \"nthread\": 4,\n         \"verbosity\": -1}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a380c449b0a45db979085c3fa0fbb4f69ab4769e"},"cell_type":"code","source":"features = X_train_scaled.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"5ed1d1c1fedafe49976eee6f18dba91cf0ffc969"},"cell_type":"code","source":"import time\nimport lightgbm as lgb\n\nfolds = KFold(n_splits=5, shuffle=True, random_state=15)\noof = np.zeros(len(X_train_scaled))\npredictions = np.zeros(len(X_test_scaled))\nstart = time.time()\nfeature_importance_df = pd.DataFrame()\n\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(X_train_scaled.values, y_train.values)):\n    print(\"fold n°{}\".format(fold_))\n    trn_data = lgb.Dataset(X_train_scaled.iloc[trn_idx][features], label=y_train.iloc[trn_idx])\n    val_data = lgb.Dataset(X_train_scaled.iloc[val_idx][features], label=y_train.iloc[val_idx])\n\n    num_round = 10000\n    clf = lgb.train(param, trn_data, num_round, valid_sets = [trn_data, val_data], verbose_eval=100, early_stopping_rounds = 200)\n    oof[val_idx] = clf.predict(X_train_scaled.iloc[val_idx][features], num_iteration=clf.best_iteration)\n    \n    fold_importance_df = pd.DataFrame()\n    fold_importance_df[\"feature\"] = features\n    fold_importance_df[\"importance\"] = clf.feature_importance()\n    fold_importance_df[\"fold\"] = fold_ + 1\n    feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n    \n    predictions += clf.predict(X_test_scaled[features], num_iteration=clf.best_iteration) / folds.n_splits\n\nprint(\"CV score: {:<8.5f}\".format(mean_squared_error(oof, y_train)**0.5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f8f90b78f9ddc400cb336e52c72a6f4bc15e8c99"},"cell_type":"code","source":"cols = (feature_importance_df[[\"feature\", \"importance\"]]\n        .groupby(\"feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:1000].index)\n\nbest_features = feature_importance_df.loc[feature_importance_df.feature.isin(cols)]\n\nplt.figure(figsize=(14,16))\nsns.barplot(x=\"importance\",\n            y=\"feature\",\n            data=best_features.sort_values(by=\"importance\",\n                                           ascending=False))\nplt.title('LightGBM Features (avg over folds)')\nplt.tight_layout()\nplt.savefig('lgbm_importances.png')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"84916e0bcf15a7bbc8ce622bf3def97ac62672ab"},"cell_type":"markdown","source":"## **2. XGboost**"},{"metadata":{"trusted":true,"_uuid":"85f8070cb4a9f7599bfe68679923933e8469c65c"},"cell_type":"code","source":"%%time\nimport xgboost as xgb\n\nxgb_params = {'eta': 0.001, 'max_depth': 5, 'subsample': 0.8, 'colsample_bytree': 0.8, 'alpha':0.1,\n          'objective': 'reg:linear', 'eval_metric': 'mae', 'silent': True, 'random_state':folds}\n\n\nfolds = KFold(n_splits=5, random_state=4520)\noof_xgb = np.zeros(len(X_train_scaled))\npredictions_xgb = np.zeros(len(X_test_scaled))\n\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(X_train_scaled.values, y_train.values)):\n    print(\"fold n°{}\".format(fold_ + 1))\n    trn_data = xgb.DMatrix(data=X_train_scaled.iloc[trn_idx][features], label=y_train.iloc[trn_idx])\n    val_data = xgb.DMatrix(data=X_train_scaled.iloc[val_idx][features], label=y_train.iloc[val_idx])\n    watchlist = [(trn_data, 'train'), (val_data, 'valid')]\n    print(\"-\" * 10 + \"Xgboost \" + str(fold_) + \"-\" * 10)\n    num_round = 11000\n    xgb_model = xgb.train(xgb_params, trn_data, num_round, watchlist, early_stopping_rounds=50, verbose_eval=1000)\n    oof_xgb[val_idx] = xgb_model.predict(xgb.DMatrix(X_train_scaled.iloc[val_idx][features]), ntree_limit=xgb_model.best_ntree_limit+50)\n\n    predictions_xgb += xgb_model.predict(xgb.DMatrix(X_test_scaled[features]), ntree_limit=xgb_model.best_ntree_limit+50) / folds.n_splits\n    \nnp.save('oof_xgb', oof_xgb)\nnp.save('predictions_xgb', predictions_xgb)\nprint(\"CV score: {:<8.5f}\".format(mean_squared_error(oof_xgb, y_train)**0.5))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"fa69f7acab27e2b6190e47fb2ab9880ed5bf8c85"},"cell_type":"code","source":"# %%time\n# from catboost import CatBoostRegressor\n# folds = KFold(n_splits=5, random_state=4520)\n# oof_cat = np.zeros(len(X_train_scaled))\n# predictions_cat = np.zeros(len(X_test_scaled))\n\n# for fold_, (trn_idx, val_idx) in enumerate(folds.split(X_train_scaled.values, y_train.values)):\n#     print(\"fold n°{}\".format(fold_ + 1))\n#     trn_data, trn_y = X_train_scaled.iloc[trn_idx][features], y_train.iloc[trn_idx]\n#     val_data, val_y = X_train_scaled.iloc[val_idx][features], y_train.iloc[val_idx]\n#     print(\"-\" * 10 + \"Catboost \" + str(fold_) + \"-\" * 10)\n#     cb_model = CatBoostRegressor(iterations=8000, learning_rate=0.01, depth=8, l2_leaf_reg=20, bootstrap_type='Bernoulli',  eval_metric='RMSE', metric_period=50, od_type='Iter', od_wait=45, random_seed=17, allow_writing_files=False)\n#     cb_model.fit(trn_data, trn_y, eval_set=(val_data, val_y), use_best_model=True, verbose=True,)\n    \n#     oof_cat[val_idx] = cb_model.predict(val_data)\n#     predictions_cat += cb_model.predict(X_test_scaled[features]) / folds.n_splits\n    \n# np.save('oof_cat', oof_cat)\n# np.save('predictions_cat', predictions_cat)\n# np.sqrt(mean_squared_error(y_train.values, oof_cat))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2765dbed40f1f34bf9507b8567473bc793cad263"},"cell_type":"markdown","source":"## **7.Stacking**"},{"metadata":{"trusted":true,"_uuid":"4d590837dfedc76f5b2448b19009147f7f5c5793"},"cell_type":"code","source":"train_stack = np.vstack([oof, oof_xgb]).transpose()\ntest_stack = np.vstack([predictions,predictions_xgb]).transpose()\n\nfolds = KFold(n_splits=5, shuffle=True, random_state=15)\noof_stack = np.zeros(train_stack.shape[0])\npredictions_stack = np.zeros(test_stack.shape[0])\n\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(train_stack, y_train)):\n    print(\"fold n°{}\".format(fold_))\n    trn_data, trn_y = train_stack[trn_idx], y_train.iloc[trn_idx].values\n    val_data, val_y = train_stack[val_idx], y_train.iloc[val_idx].values\n\n    print(\"-\" * 10 + \"Ridge Regression\" + str(fold_) + \"-\" * 10)\n#     cb_model = CatBoostRegressor(iterations=3000, learning_rate=0.1, depth=8, l2_leaf_reg=20, bootstrap_type='Bernoulli',  eval_metric='RMSE', metric_period=50, od_type='Iter', od_wait=45, random_seed=17, allow_writing_files=False)\n#     cb_model.fit(trn_data, trn_y, eval_set=(val_data, val_y), cat_features=[], use_best_model=True, verbose=True)\n    clf = BayesianRidge()\n    clf.fit(trn_data, trn_y)\n    \n    oof_stack[val_idx] = clf.predict(val_data)\n    predictions_stack += clf.predict(test_stack) / 5\n\n\nprint(\"CV score: {:<8.5f}\".format(mean_squared_error(oof, y_train)**0.5))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3a5b1849731681cb075e90537daebdcb775e1ed4"},"cell_type":"markdown","source":"## **8.Final Prediction**"},{"metadata":{"trusted":true,"_uuid":"dd46c9e2c3ec68f45900cc49f95aee82459137a1"},"cell_type":"code","source":"sample_submission = pd.read_csv('../input/sample_submission.csv')\nsample_submission['time_to_failure'] = predictions_stack\nsample_submission.to_csv('Bayesian_Ridge_Stacking.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b2f170d0b38206487afc2b8337dd2f1038924483"},"cell_type":"code","source":"sample_submission.shape","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}