{"cells":[{"metadata":{},"cell_type":"markdown","source":"<center><img src=\"https://images.ctfassets.net/81iqaqpfd8fy/3Wp4SEgzagcICaSqcIMOQM/5721655abf93a19521dad8a35d747f2d/Erupting_Volcano.jpg?h=620&w=1024\"></center>\n<h1><center>INGV - Volcanic Eruption Prediction</center></h1>\n<h1><center>Faster Simple Baseline Using Datatable</center></h1>"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport glob\nfrom tqdm import tqdm\n\nimport time\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.linear_model import Ridge, RidgeCV\n\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom catboost import CatBoostRegressor\n\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport matplotlib.pyplot as plt\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check Data"},{"metadata":{},"cell_type":"markdown","source":"## Check Files"},{"metadata":{"trusted":true},"cell_type":"code","source":"segment_csvs = glob.glob(\"../input/predict-volcanic-eruptions-ingv-oe/train/*\")\nlen_segment_csvs = len(segment_csvs)\nlen_segment_csvs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_csvs = glob.glob(\"../input/predict-volcanic-eruptions-ingv-oe/test/*\")\nlen_test_csvs = test_csvs\nlen(len_test_csvs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(sample_submission)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segment_csvs[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_1 = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train/2037160701.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can assume the `sampling rate` at `100Hz` because it has 60000 rows."},{"metadata":{},"cell_type":"markdown","source":"## Show Sensor values"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ndef sensor_show(df):\n    f, axes = plt.subplots(10, 1)\n    f.set_size_inches((16, 8)) \n    f.tight_layout() \n    plt.subplots_adjust(bottom=-0.4)\n    \n    # Sensor#1 ~ #10\n    for i in range(1,11):\n        axes[i-1].plot(df[f'sensor_{i}'].values)\n        axes[i-1].set_title('Sensor_'+str(i))\n        axes[i-1].set_xlabel('time')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_show(train_1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Baseline\n\n- https://www.kaggle.com/mahmoudvaziri/svm-regression-mean\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# datatable installation with internet\n#!pip install datatable==0.11.0 > /dev/null\n\n# installation without internet\n!pip install ../input/python-datatable/datatable-0.11.0-cp37-cp37m-manylinux2010_x86_64.whl\n\nimport datatable as dt","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Using datatable\n- More faster than original baseline\n"},{"metadata":{},"cell_type":"markdown","source":"<center><img src=\"https://miro.medium.com/max/446/0*w7dsjAY9CKNY7owL.png?h=620&w=1024\"></center>\n\n<p></p>\n\n- https://www.kaggle.com/rohanrao/riiid-with-blazing-fast-rid\n\nThis notebook shows how you can use [Python datatable](https://datatable.readthedocs.io/en/latest/index.html) to read the complete training data and convert it to a pandas dataframe in under a minute.\n"},{"metadata":{},"cell_type":"markdown","source":"In Original kernels, It takes about `7:23 sec`.\n\nBut using datatable, it takes only `3:17 sec`."},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv(segment_csvs[0])\ndf_mean = pd.DataFrame(df.mean()).T\ndf_mean['id'] = segment_csvs[0].split('/')[-1].split('.')[0]\nsegment_csvs.remove(segment_csvs[0])\n\nfor csv in tqdm(segment_csvs):\n    seg_name = csv.split('/')[-1].split('.')[0]\n    df = dt.fread(csv).to_jay('train.jay')\n    df = dt.fread('train.jay')\n    df_ = pd.DataFrame(df.mean().to_pandas()) # df.mean() for datatable\n    df_['id'] = csv.split('/')[-1].split('.')[0]\n    df_mean = pd.concat([df_mean,df_])\n    del df\n    \ndf_mean.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\ndf_train.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mean['id'] = df_mean['id'].astype('int64')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mean = df_mean.join(df_train.set_index('segment_id'), on='id')\ndf_mean.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = df_mean.drop(['id','time_to_eruption'],axis=1)\ny_train = df_mean['time_to_eruption']\nX_train = X_train.fillna(X_train.mean())\ndel df_mean","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_train)\nX_train_scaled = pd.DataFrame(scaler.transform(X_train), columns=X_train.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_scaled","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv(test_csvs[0])\ndf_mean_test = pd.DataFrame(df.mean()).T\ndf_mean_test['id'] = test_csvs[0].split('/')[-1].split('.')[0]\ntest_csvs.remove(test_csvs[0])\n\nfor csv in tqdm(test_csvs):\n    df = dt.fread(csv).to_jay('test.jay')\n    df = dt.fread(\"test.jay\")\n    df_ = pd.DataFrame(df.mean().to_pandas()) # df.mean() for datatable\n    df_['id'] = csv.split('/')[-1].split('.')[0]\n    df_mean_test = pd.concat([df_mean_test,df_])\n    del df\ndf_mean_test.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = df_mean_test.fillna(df_mean_test.mean())\nX_test = X_test.drop(['id'],axis=1)\n#del df_mean_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test_scaled","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_set = pd.DataFrame()\ntrain_set['segment_id'] = df_train.segment_id\ntrain_set = train_set.set_index('segment_id')\ntrain_set = pd.merge(train_set.reset_index(), df_train, on=['segment_id'], how='left').set_index('segment_id')\n\ny_train = train_set['time_to_eruption']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- https://www.kaggle.com/artgor/seismic-data-eda-and-baseline"},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_model(X=X_train_scaled, X_test=X_test_scaled, y=y_train, 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_train.columns)\n            valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X_train.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_train.columns), ntree_limit=model.best_ntree_limit)\n            y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X_train.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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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},"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')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","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":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 8))\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');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction_lgb[:10], prediction_xgb[:10], prediction_svr[:10]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame()\nsubmission['segment_id'] = sample_submission.segment_id\nsubmission['time_to_eruption'] = (prediction_lgb + prediction_xgb + prediction_svr) / 3\nprint(submission.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I think that blow kernel also can be more faster using datatable.\n- https://www.kaggle.com/ajcostarino/ingv-volcanic-eruption-prediction-lgbm-baseline"}],"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":4,"nbformat_minor":4}