{"cells":[{"metadata":{},"cell_type":"markdown","source":"Thanks for the notebook from https://www.kaggle.com/ajcostarino/ingv-volcanic-eruption-prediction-lgbm-baseline @ajcostarino\n\nIn this notebook, I resampled time series from 60000 rows to 6000 rows by rolling mean (window=10) in order to save the time of processing data.\n\nHowever, out of my expectation, it also helps improve the performance of the prediction."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport plotly.graph_objs as go\nimport plotly\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nimport cufflinks as cf\ncf.set_config_file(offline=True)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Take segment_id:2037160701 as an example, you can compare the patterns before and after reampling:"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_2037160701 = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/2037160701.csv')\ndf_2037160701.iplot(title='before resampling') #before resampling\ndf_2037160701.rolling(10).mean().iloc[list(np.arange(10,60001,10))].iplot(title='after resampling') #after resampling","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Baseline model for Volcanic time to eruption\n\nThis notebook is a simpler, single model, version of what I used in the Lanl Earthqauke competition here (https://www.kaggle.com/c/LANL-Earthquake-Prediction). We collapsed every observation into a single row, calculate a couple of overall metrics on each of the sensors that are provided (1-10). And then we use Light GBM regressor with some early stopping conditions to calculate time to eruption.\n\nA couple of functions follow below: Basic statistics will calculate distribution values for each sensor. Quantiles will calculate quantiles for each sensor. Linear regression builds a linear regression model for each sensor and returns the coefficients. I have provided some other useful functions such as Fast-Fourier transforms but I do not use to make this prediction."},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\n\nfrom sklearn.linear_model import LinearRegression\nimport scipy.stats as spstats\n\n\ndef basic_statistics(t_X, x, s, sensor, postfix=''):\n    \"\"\"Computes basic statistics for the training feature set.\n    \n    Args:\n        t_X (pandas.DataFrame): The feature set being built.\n        x (pandas.Series): The signal values.\n        s (int): The integer number of the segment.\n        postfix (str): The postfix string value.\n    Return:\n        t_X (pandas.DataFrame): The feature set being built.\n    \"\"\"\n\n    t_X.loc[s, f'{sensor}_sum{postfix}']       = x.sum()\n    t_X.loc[s, f'{sensor}_mean{postfix}']      = x.mean()\n    t_X.loc[s, f'{sensor}_std{postfix}']       = x.std()\n    t_X.loc[s, f'{sensor}_var{postfix}']       = x.var() \n    t_X.loc[s, f'{sensor}_max{postfix}']       = x.max()\n    t_X.loc[s, f'{sensor}_min{postfix}']       = x.min()\n    t_X.loc[s, f'{sensor}_median{postfix}']    = x.median()\n    t_X.loc[s, f'{sensor}_skew{postfix}']      = x.skew()\n    t_X.loc[s, f'{sensor}_mad{postfix}']       = x.mad()\n    t_X.loc[s, f'{sensor}_kurtosis{postfix}']  = x.kurtosis()\n\n    return t_X\n\n\n\ndef quantiles(t_X, x, s, sensor, postfix=''):\n    \"\"\"Calculates quantile features for the training feature set.\n    Args:\n        t_X (pandas.DataFrame): The feature set being built.\n        x (pandas.Series): The signal values.\n        s (int): The integer number of the segment.\n        postfix (str): The postfix string value.\n    Return:\n        t_X (pandas.DataFrame): The feature set being built.\n    \"\"\"\n    t_X.loc[s, f'{sensor}_q999{postfix}']     = np.quantile(x ,0.999)\n    t_X.loc[s, f'{sensor}_q99{postfix}']      = np.quantile(x, 0.99)\n    t_X.loc[s, f'{sensor}_q95{postfix}']      = np.quantile(x, 0.95)\n    t_X.loc[s, f'{sensor}_q87{postfix}']      = np.quantile(x, 0.87)\n    t_X.loc[s, f'{sensor}_q13{postfix}']      = np.quantile(x, 0.13)  \n    t_X.loc[s, f'{sensor}_q05{postfix}']      = np.quantile(x, 0.05)\n    t_X.loc[s, f'{sensor}_q01{postfix}']      = np.quantile(x, 0.01)\n    t_X.loc[s, f'{sensor}_q001{postfix}']     = np.quantile(x ,0.001)\n    \n    x_abs = np.abs(x)\n    t_X.loc[s, f'{sensor}_q999_abs{postfix}'] = np.quantile(x_abs, 0.999)\n    t_X.loc[s, f'{sensor}_q99_abs{postfix}']  = np.quantile(x_abs, 0.99)\n    t_X.loc[s, f'{sensor}_q95_abs{postfix}']  = np.quantile(x_abs, 0.95)\n    t_X.loc[s, f'{sensor}_q87_abs{postfix}']  = np.quantile(x_abs, 0.87)\n    t_X.loc[s, f'{sensor}_q13_abs{postfix}']  = np.quantile(x_abs, 0.13)\n    t_X.loc[s, f'{sensor}_q05_abs{postfix}']  = np.quantile(x_abs, 0.05)\n    t_X.loc[s, f'{sensor}_q01_abs{postfix}']  = np.quantile(x_abs, 0.01)\n    t_X.loc[s, f'{sensor}_q001_abs{postfix}'] = np.quantile(x_abs, 0.001)\n    \n    t_X.loc[s, f'{sensor}_iqr']     = np.subtract(*np.percentile(x, [75, 25]))\n    t_X.loc[s, f'{sensor}_iqr_abs'] = np.subtract(*np.percentile(x_abs, [75, 25]))\n\n    return t_X\n\n\ndef __linear_regression(arr, abs_v=False):\n    \"\"\"\n    \"\"\"\n    idx = np.array(range(len(arr)))\n    if abs_v:\n        arr = np.abs(arr)\n    lr = LinearRegression()\n    fit_X = idx.reshape(-1, 1)\n    lr.fit(fit_X, arr)\n    return lr.coef_[0]\n\n\ndef __classic_sta_lta(x, length_sta, length_lta):\n    sta = np.cumsum(x ** 2)\n    # Convert to float\n    sta = np.require(sta, dtype=np.float)\n    # Copy for LTA\n    lta = sta.copy()\n    # Compute the STA and the LTA\n    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n    sta /= length_sta\n    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n    lta /= length_lta\n    # Pad zeros\n    sta[:length_lta - 1] = 0\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    return sta / lta\n\n\ndef linear_regression(t_X, x, s, sensor, postfix=''):\n    t_X.loc[s, f'{sensor}_lr_coef{postfix}'] = __linear_regression(x)\n    t_X.loc[s, f'{sensor}_lr_coef_abs{postfix}'] = __linear_regression(x, True)\n    return t_X\n\n\ndef classic_sta_lta(t_X, x, sensor, s):\n    t_X.loc[s, f'{sensor}_classic_sta_lta1_mean'] = __classic_sta_lta(x, 500, 10000).mean()\n    t_X.loc[s, f'{sensor}_classic_sta_lta2_mean'] = __classic_sta_lta(x, 5000, 100000).mean()\n    t_X.loc[s, f'{sensor}_classic_sta_lta3_mean'] = __classic_sta_lta(x, 3333, 6666).mean()\n    t_X.loc[s, f'{sensor}_classic_sta_lta4_mean'] = __classic_sta_lta(x, 10000, 25000).mean()\n    return t_X\n\n\ndef fft(t_X, x, s, sensor, postfix=''):\n    \"\"\"Generates basic statistics over the fft of the signal\"\"\"\n    z = np.fft.fft(x)\n    fft_real = np.real(z)\n    fft_imag = np.imag(z)\n\n    t_X.loc[s, f'fft_A0']             = abs(z[0])\n    \n    t_X.loc[s, f'{sensor}_fft_real_mean{postfix}']      = fft_real.mean()\n    t_X.loc[s, f'{sensor}_fft_real_std{postfix}']       = fft_real.std()\n    t_X.loc[s, f'{sensor}_fft_real_max{postfix}']       = fft_real.max()\n    t_X.loc[s, f'{sensor}_fft_real_min{postfix}']       = fft_real.min()\n    t_X.loc[s, f'{sensor}_fft_real_median{postfix}']    = np.median(fft_real)\n    t_X.loc[s, f'{sensor}_fft_real_skew{postfix}']      = spstats.skew(fft_real)\n    t_X.loc[s, f'{sensor}_fft_real_kurtosis{postfix}']  = spstats.kurtosis(fft_real)\n    \n    t_X.loc[s, f'{sensor}_fft_imag_mean{postfix}']      = fft_imag.mean()\n    t_X.loc[s, f'{sensor}_fft_imag_std{postfix}']       = fft_imag.std()\n    t_X.loc[s, f'{sensor}_fft_imag_max{postfix}']       = fft_imag.max()\n    t_X.loc[s, f'{sensor}_fft_imag_min{postfix}']       = fft_imag.min()\n    t_X.loc[s, f'{sensor}_fft_imag_median{postfix}']    = np.median(fft_imag)\n    t_X.loc[s, f'{sensor}_fft_imag_skew{postfix}']      = spstats.skew(fft_imag)\n    t_X.loc[s, f'{sensor}_fft_imag_kurtosis{postfix}']  = spstats.kurtosis(fft_imag)\n    \n    return t_X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Build out the features for the training data"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/train.csv')\ntrain_set = pd.DataFrame()\ntrain_set['segment_id'] = train.segment_id\ntrain_set = train_set.set_index('segment_id')\n\nj = 0\nfor seg in train.segment_id:\n    signals = pd.read_csv(f'/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/{seg}.csv')\n    signals = signals.rolling(10).mean().iloc[list(np.arange(10,60001,10))]\n    for i in range(1, 11):\n        sensor_id = f'sensor_{i}'\n        train_set = basic_statistics(train_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        train_set = quantiles(train_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        train_set = linear_regression(train_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        train_set = fft(train_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        \ntrain_set = pd.merge(train_set.reset_index(), train, on=['segment_id'], how='left').set_index('segment_id')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Build out the features for the testing data"},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv')\ntest_set = pd.DataFrame()\ntest_set['segment_id'] = test.segment_id\ntest_set = test_set.set_index('segment_id')\n\n\nfor seg in test.segment_id:\n    signals = pd.read_csv(f'/kaggle/input/predict-volcanic-eruptions-ingv-oe/test/{seg}.csv')\n    signals = signals.rolling(10).mean().iloc[list(np.arange(10,60001,10))]\n    \n    for i in range(1, 11):\n        sensor_id = f'sensor_{i}'\n        test_set = basic_statistics(test_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        test_set = quantiles(test_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        test_set = linear_regression(test_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')\n        test_set = fft(test_set, signals[sensor_id].fillna(0), seg, sensor_id, postfix='')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Build out the Light GBM Regression model using MAE as the eval metric"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgbm\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\n\ny = train_set['time_to_eruption']\nfeature_df = train_set.drop(['time_to_eruption'], axis = 1)\n\nscaler = StandardScaler()\nscaler.fit(feature_df)\nscaled_feature_df = pd.DataFrame(scaler.transform(feature_df), columns=feature_df.columns)\nscaled_test_df    = pd.DataFrame(scaler.transform(test_set), columns=test_set.columns)\n\nprint(scaled_feature_df.shape)\nprint(scaled_test_df.shape)\n\n\nn_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\nscaled_feature_df_columns = scaled_feature_df.columns.values\n\n\nparams = {\n    'num_leaves': 85,\n    'min_data_in_leaf': 10, \n    'objective':'regression',\n    'max_depth': -1,\n    'learning_rate': 0.001,\n    'max_bins': 2048,\n    \"boosting\": \"gbdt\",\n    \"feature_fraction\": 0.91,\n    \"bagging_freq\": 1,\n    \"bagging_fraction\": 0.91,\n    \"bagging_seed\": 42,\n    \"metric\": 'mae',\n    \"lambda_l1\": 0.1,\n    \"verbosity\": -1,\n    \"nthread\": -1,\n    \"random_state\": 42\n}\n\n\noof = np.zeros(len(scaled_feature_df))\npredictions = np.zeros(len(scaled_test_df))\nfeature_importance_df = pd.DataFrame()\n\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(scaled_feature_df, y.values)):\n    \n    strLog = \"fold {}\".format(fold_)\n    print(strLog)\n    \n    X_tr, X_val = scaled_feature_df.iloc[trn_idx], scaled_feature_df.iloc[val_idx]\n    y_tr, y_val = y.iloc[trn_idx], y.iloc[val_idx]\n\n    model = lgbm.LGBMRegressor(**params, n_estimators = 20000, n_jobs = -1)\n    model.fit(X_tr, y_tr, \n              eval_set=[(X_tr, y_tr), (X_val, y_val)], eval_metric='mae',\n              verbose=1000, early_stopping_rounds=400)\n    \n    oof[val_idx] = model.predict(X_val, num_iteration=model.best_iteration_)\n\n    fold_importance_df = pd.DataFrame()\n    fold_importance_df[\"Feature\"] = scaled_feature_df_columns\n    fold_importance_df[\"importance\"] = model.feature_importances_[:len(scaled_feature_df_columns)]\n    fold_importance_df[\"fold\"] = fold_ + 1\n    feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n    #predictions\n    predictions += model.predict(scaled_test_df, num_iteration=model.best_iteration_) / folds.n_splits","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Display some feature importances"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\n\ncols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:3014].index)\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\n\nplt.figure(figsize=(14,26*3))\nsns.barplot(x=\"importance\", y=\"Feature\", data=best_features.sort_values(by=\"importance\",ascending=False))\nplt.title('LightGBM Features (averaged over folds)')\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submit"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame()\nsubmission['segment_id'] = test_set.index\nsubmission['time_to_eruption'] = predictions\nsubmission.to_csv('submission_recent.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_df.to_pickle('feature_df.pickle')\ntest_set.to_pickle('test_set.pickle')\nscaled_feature_df.to_pickle('scaled_feature_df.pickle')\nscaled_test_df.to_pickle('scaled_test_df.pickle')","execution_count":null,"outputs":[]}],"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}