{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Earthquake time prediction using XGB\n\nThis kernel manipulates with a dataset that was generated using the kernel https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples\n\n* [Additional methods](#addmethods)\n* [Datasets loading and preparation](#data)\n* [Feature selection](#feature_selection)\n* [Cross-validation strategy](#CV)\n* [Validation and Learning Curves](#curves)\n* [Parameters tuning](#tuning)\n    * [Exhaustive tuning with RandomizedSearchCV](#rand)\n    * [Fine tuning with GridSearcCV](#grid)\n* [Make a prediction using the chosen CV strategy and tuned algorithm](#preds)    \n\n22.06.2019\n    "},{"metadata":{"trusted":true},"cell_type":"code","source":"import multiprocessing\nn_jobs = multiprocessing.cpu_count()-1\n\nimport numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA as RandomizedPCA\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import r2_score\n\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.model_selection import RandomizedSearchCV\nfrom sklearn.model_selection import cross_val_score, cross_validate\nfrom sklearn.model_selection import learning_curve, validation_curve\nfrom sklearn.model_selection import KFold, ShuffleSplit, StratifiedShuffleSplit\n\nfrom time import time, ctime\n\nimport xgboost","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Additional methods** <a class=\"anchor\" id=\"addmethods\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"def changeparam(value, param):\n    if param == 'max_depth' or param == 'min_child_weight':\n        if value >= 1 and value <= 10:\n            return list(filter(lambda x: x != 0, [ value + i for i in [-1, 0, 1]]))\n        else: \n            return value\n    elif param == 'n_estimators':\n        if value >= 1 and value < 10:\n            return list(filter(lambda x: x != 0, [ value + i for i in [-1, 0, 1]]))\n        elif value >= 10 and value < 100:\n            return [ value + i for i in [-5, 0, 5]]\n        elif value >= 100 and value < 1000:\n            return [ value + i for i in [-20, -10, 0, 20, 10]]\n        elif value >= 1000 and value < 10000:\n            return [ value + i for i in [-200, -100, 0, 100, 200]]\n        elif value >= 10000:\n            return [ value + i for i in [-2000, -1000, 0, 1000, 2000]]\n\ndef plotCurves(model, X, Y, param_range, param_name, scoring):\n    best_values = []\n \n    for i, st_i in enumerate(scoring):\n            fig, ax = plt.subplots( nrows = 1, ncols = 2, figsize=(18, 5))\n\n            train_scores, test_scores = validation_curve(\n                model, X, Y.values.ravel(), param_name=param_name, param_range=param_range,\n                cv=cv, scoring=scoring[i], n_jobs = n_jobs, error_score=0)\n\n            train_scores_mean = np.mean(train_scores, axis=1)\n            train_scores_std = np.std(train_scores, axis=1)\n\n            test_scores_mean = np.mean(test_scores, axis=1)\n            test_scores_std = np.std(test_scores, axis=1)\n\n            ax[0].plot(param_range, train_scores_mean, label=\"Training score\", color=\"darkorange\")\n            ax[0].fill_between(param_range, train_scores_mean - train_scores_std,\n                         train_scores_mean + train_scores_std, alpha=0.2,\n                         color=\"darkorange\")\n\n            ax[0].plot(param_range, test_scores_mean, label=\"Cross-validation score\", color=\"navy\")\n            ax[0].fill_between(param_range, test_scores_mean - test_scores_std,\n                         test_scores_mean + test_scores_std, alpha=0.2,\n                         color=\"navy\")\n            ax[0].legend(loc=\"best\")\n\n            ax[0].set_title('Validation Curve for '+scoring[i])\n            ax[0].set_ylabel(scoring[i])\n            ax[0].set_xlabel(param_name)\n\n            ax[0].set_ylim(test_scores_mean[-1] - abs(6*test_scores_std[-1]), \n                           test_scores_mean[-1] + abs(6*test_scores_std[-1]))\n            ax[0].grid(True)\n\n            best_param = list(param_range)[np.argmax(test_scores_mean)]\n            \n            if param_name == 'n_estimators':\n                model.n_estimators = best_param\n            elif param_name == 'max_depth':\n                model.max_depth = best_param\n            elif param_name == 'min_child_weight':\n                model.min_child_weight = best_param\n\n            best_values.append(best_param)\n\n            sizes = np.linspace(.1, 1.0, 10)\n            train_sizes, train_scores, test_scores = learning_curve(model, \n                                                                    X, Y.values.ravel(), \n                                                                    cv = cv, \n                                                                    scoring=scoring[i],\n                                                                    n_jobs = n_jobs, \n                                                                    train_sizes = sizes,\n                                                                    error_score=0)\n            train_scores_mean = np.mean(train_scores, axis=1)\n            train_scores_std = np.std(train_scores, axis=1)\n            test_scores_mean = np.mean(test_scores, axis=1)\n            test_scores_std = np.std(test_scores, axis=1)\n\n            ax[1].fill_between(sizes, train_scores_mean - train_scores_std,\n                             train_scores_mean + train_scores_std, alpha=0.1,\n                             color=\"r\")\n            ax[1].fill_between(sizes, test_scores_mean - test_scores_std,\n                             test_scores_mean + test_scores_std, alpha=0.1, color=\"g\")\n            ax[1].plot(sizes, train_scores_mean, 'o-', color=\"r\",\n                     label=\"Training score\")\n            ax[1].plot(sizes, test_scores_mean, 'o-', color=\"g\",\n                     label=\"Cross-validation score\")\n\n            ax[1].set_title('Learning curve for '+scoring[i] +' for ' + param_name + ' = '+ str(best_param) )\n            ax[1].set_ylabel(scoring[i])\n            ax[1].set_xlabel('Number training observations')\n            ax[1].grid(True)\n            ax[1].legend(loc=\"best\")\n\n            plt.pause(0.01)\n    fig.tight_layout()\n    return best_values\n\ndef plotfig (ypred, yactual, strtitle):\n    plt.scatter(ypred, yactual.values.ravel())\n    plt.title(strtitle)\n    plt.plot([(0, 0), (20, 20)], [(0, 0), (20, 20)])\n    plt.xlim(0, 20)\n    plt.ylim(0, 20)\n    plt.xlabel('Predicted', fontsize=12)\n    plt.ylabel('Actual', fontsize=12)\n    plt.show()\n\n# Utility function to report best scores\ndef report(results, n_top=3):\n    for i in range(1, n_top + 1):\n        candidates = np.flatnonzero(results['rank_test_score'] == i)\n        for candidate in candidates:\n            print(\"Model with rank: {0}\".format(i))\n            print(\"Mean validation score: {0:.3f} (std: {1:.3f})\".format(\n                  results['mean_test_score'][candidate],\n                  results['std_test_score'][candidate]))\n            print(\"Parameters: {0}\".format(results['params'][candidate]))\n            print(\"\")\n            \ndef mode_custom(List): \n    dict = {} \n    count, itm = 0, '' \n    for item in reversed(List): \n        dict[item] = dict.get(item, 0) + 1\n        if dict[item] >= count : \n            count, itm = dict[item], item \n    return(itm) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Datasets loading and preparation** <a class=\"anchor\" id=\"data\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"XY = pd.read_csv('../input/LANL_train.csv')\nX_TEST = pd.read_csv('../input/LANL_test.csv')\ncol = [c for c in XY.columns if c not in ['time_to_failure']]\nX = XY[col]\nY = XY['time_to_failure']\nprint(X.shape, X_TEST.shape, Y.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(X.isnull().values.any())\nprint(X_TEST.isnull().values.any())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check the number of features with one unique value \nallunique = X.nunique().reset_index()\nlst = [allunique.loc[i,'index'] for i in range(len(allunique)) if allunique.loc[i,0] == 1]\nprint(len(lst), lst)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Drop colums with one unique value (in this case this is null)\nX.drop(lst, axis = 1, inplace = True)\nX_TEST.drop(lst, axis = 1, inplace = True)\nprint(X.shape, X_TEST.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = StandardScaler()\n\nscaler.fit(X)\nX = scaler.transform(X)\n\nscaler.fit(X_TEST)\nX_TEST = scaler.transform(X_TEST)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Feature selection** <a class=\"anchor\" id=\"feature_selection\"></a>"},{"metadata":{},"cell_type":"markdown","source":"**Bi-objective feature selection**\n\nThe features were selected using a metaheuristic optimization approach which idea had been presented in https://iopscience.iop.org/article/10.1088/1742-6596/1210/1/012086/meta. We were considering two objectives: maximization of R-squared and features number minimization. In order to check the idea an get a solution faster, we performed the experiment for Linear Regression model. So, we got 79 features that allow getting almost the same private score, but with less time consumption. \n\nFor demonstration purpose, in the figure below you can see the results of the first phase of the optimization procedure that allows concluding the best R-squared is about 0.52 with approximately 80 features.  \nhttps://github.com/LyubAlex/kaggle/blob/master/LANL%20Earthquake%20Prediction/phase1.png\n<img src=\"phase1.png\" width=\"500\" class=\"image left\" />\n\n \n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# features = [5,6,7,9,10,11,12,16,18,19,21,22,24,26,28,29,32,35,37,38,39,41,42,44,45,\n#             46,47,48,49,50,51,52,54,55,56,57,58,62,64,66,68,69,71,72,73,74,76,78,79,\n#             82,86,89,90,91,96,97,98,99,100,102,103,104,105,107,108,109,110,113,120,\n#             121,123,126,127,128,129,130,131,133,134]\n# X = X.iloc[:,features]\n# X_TEST = X_TEST.iloc[:,features]\n# print(X.shape, X_TEST.shape, Y.shape)\n\n# X = scaler.fit(X).transform(X)\n# X_TEST = scaler.fit(X_TEST).transform(X_TEST)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**PCA**\n\nThis approach of feature selection worsen the score. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# pca = RandomizedPCA(copy = True, \n#                     iterated_power = 3,\n#                     n_components = 55, \n#                     svd_solver='randomized', \n#                     random_state = 0, \n#                     whiten=False).fit(X)\n\n# plt.plot(np.cumsum(pca.explained_variance_ratio_))\n# plt.grid()\n\n# X = pca.transform(X)\n# print('Dimension of train dataset', X.shape)\n\n# X_TEST = pca.transform(X_TEST)\n# print('Dimension of test dataset', X_TEST.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Cross-validation strategy** <a class=\"anchor\" id=\"CV\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 5\n\ncv = ShuffleSplit(n_splits=n_fold, test_size=0.4, random_state = 0)\n# cv = KFold(n_splits=n_fold, shuffle=True, random_state = 0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**XGBoost initialization**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# scoring = ['r2', 'neg_mean_absolute_error', 'neg_mean_squared_error', \n#            'neg_mean_squared_log_error','neg_median_absolute_error',\n#            'explained_variance']\n\nscoring = ['r2', 'neg_mean_absolute_error', 'neg_mean_squared_error', \n           'neg_median_absolute_error', 'explained_variance']\n\nn_estimators_range = range(10, 5000, 500) \nmax_depth_range = np.arange(1, 11, 2)\nmin_child_weight_range = np.arange(1, 11, 2)\n\nmodel = xgboost.XGBRegressor(learning_rate = 0.01, \n                            criterion = 'mae',\n                            subsample = 1,\n                            colsample_bytree = 1,\n                            max_depth = 1,\n                            min_child_weight = 1\n                            )\n\nmodel","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Validation and Learning Curves <a class=\"anchor\" id=\"curves\"></a>"},{"metadata":{},"cell_type":"markdown","source":"1. On **n_estimators** parameter"},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"best_n_estimators = plotCurves(model, X, Y, n_estimators_range, 'n_estimators', scoring)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.n_estimators = mode_custom(best_n_estimators)\nmodel","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"2. On **max_depth** parameter"},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"best_max_depth = plotCurves(model, X, Y, max_depth_range, 'max_depth', scoring)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.max_depth = mode_custom(best_max_depth)\nmodel","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"3. On **min_child_weight** parameter"},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"best_min_child_weight = plotCurves(model, X, Y, min_child_weight_range, 'min_child_weight', scoring)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.min_child_weight = mode_custom(best_min_child_weight)\nch_algorithm = model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Parameters tuning <a class=\"anchor\" id=\"tuning\"></a>"},{"metadata":{},"cell_type":"markdown","source":"**Exhaustive tuning with RandomizedSearchCV** <a class=\"anchor\" id=\"rand\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"# random_grid = {'n_estimators': n_estimators_range,\n#                'learning_rate': [0.01],\n#                'max_depth': max_depth_range,\n#                'min_child_weight': min_child_weight_range,\n#                'subsample' : [1],\n#                'colsample_bytree' : [1],\n#                'criterion' : ['mae'],\n#               }\n\n# model = xgboost.XGBRegressor()\n\n# n_iter_search = 20\n# random_search = RandomizedSearchCV(estimator = model, \n#                                     param_distributions = random_grid, \n#                                     n_iter = n_iter_search, \n#                                     cv = cv, \n#                                     scoring = 'neg_mean_absolute_error',\n#                                     verbose = 10, \n#                                     error_score = 0,\n#                                     random_state = 42,\n#                                     n_jobs = n_jobs)\n\n# start = time()\n# random_search.fit(X, Y)\n# print(\"RandomizedSearchCV took %.2f seconds for %d candidates\"\n#       \" parameter settings.\" % ((time() - start), n_iter_search))\n# report(random_search.cv_results_)\n\n# ch_algorithm = random_search.best_estimator_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Fine tuning with GridSearcCV** <a class=\"anchor\" id=\"grid\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nsearch_grid = {'n_estimators': changeparam(ch_algorithm.n_estimators, 'n_estimators'),\n               'learning_rate': [0.01],\n               'max_depth': changeparam(ch_algorithm.max_depth, 'max_depth'),\n               'min_child_weight': changeparam(ch_algorithm.min_child_weight, 'min_child_weight'),\n               'subsample' : [0.5, 1],\n               'colsample_bytree' : [0.5, 1],\n               'criterion' : ['mae'],\n              }\n\nn_iter_search = sum([len(value) for key, value in search_grid.items()])\nmodel = xgboost.XGBRegressor()\n\ngrid = GridSearchCV(model, \n                    search_grid, \n                    cv=cv, \n                    scoring = 'neg_mean_absolute_error', \n                    verbose = 10,\n                    error_score = 0,\n                    n_jobs = n_jobs)\n\nstart = time()\ngrid.fit(X, Y)\nprint(\"RandomizedSearchCV took %.2f seconds for %d candidates\"\n      \" parameter settings.\" % ((time() - start), n_iter_search))\n\nreport(grid.cv_results_)\n\nch_algorithm = grid.best_estimator_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ch_algorithm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ch_algorithm.n_estimators =ch_algorithm.n_estimators*2\n# ch_algorithm.learning_rate = ch_algorithm.learning_rate / 2\nch_algorithm.n_jobs = n_jobs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ch_algorithm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Make a prediction using the chosen CV strategy and tuned algorithm <a class=\"anchor\" id=\"preds\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nsubmission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\noof = np.zeros(len(X))\nprediction = np.zeros(len(submission))\nmae, r2 = [], []\n\nfor fold_n, (train_index, valid_index) in enumerate(cv.split(X)):\n    print('\\nFold', fold_n, 'started at', ctime())\n\n    X_train = X[train_index]\n    X_valid = X[valid_index]\n    Y_train = Y.values.ravel()[train_index]\n    Y_valid = Y.values.ravel()[valid_index]\n       \n    best_model = ch_algorithm.fit(X_train, Y_train)\n    y_pred = best_model.predict(X_valid)   \n  \n    oof[valid_index] = y_pred\n\n    mae.append(mean_absolute_error(Y_valid, y_pred))\n    r2.append(r2_score(Y_valid, y_pred))\n\n    print('MAE: ', mean_absolute_error(Y_valid, y_pred))\n    print('R2: ', r2_score(Y_valid, y_pred))\n\n    prediction += best_model.predict(X_TEST)\n        \nprediction /= n_fold\n\nprint('='*45)\nprint('CV mean MAE: {0:.4f}, std: {1:.4f}.'.format(np.mean(mae), np.std(mae)))\nprint('CV mean R2:  {0:.4f}, std: {1:.4f}.'.format(np.mean(r2), np.std(r2)))\n\nplotfig(best_model.predict(X), Y, 'Predicted vs. Actual responses for XGB')\n\n# non_zeros_ind = np.argwhere(oof != 0)[:, 0]\n# oof     = oof[non_zeros_ind]\n# Y_fixed = pd.DataFrame(Y.values[non_zeros_ind])\n# plotfig(oof, Y_fixed, 'Predicted vs. Actual responses for XGB')\n# print('\\nMAE: ', mean_absolute_error(Y_fixed, oof))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# scores = cross_validate(best_model, X, Y, cv=cv, scoring = {'neg_mean_absolute_error', 'r2'})\n\n# print('CV mean MAE: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores['test_neg_mean_absolute_error']), \n#                                                    np.std(scores['test_neg_mean_absolute_error'])))\n# print('CV mean R2:  {0:.4f}, std: {1:.4f}.'.format(np.mean(scores['test_r2']), \n#                                                    np.std(scores['test_r2'])))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission['time_to_failure'] = prediction \nprint(submission.head())\nsubmission.to_csv('submission.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.1"}},"nbformat":4,"nbformat_minor":1}