{"cells":[{"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\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\nimport os\nfor 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 20GB 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":{"trusted":true},"cell_type":"code","source":"# more necessary imports\nfrom tqdm import tqdm, trange\nimport time\nfrom pathlib import Path\nimport dask.dataframe as dd\n\nimport matplotlib.pyplot as plt\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler\n\n\n#import lightgbm as lgb\nfrom optuna.integration import lightgbm as lgb\n\nimport glob","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# number of train and test files\ntrains = []\ntests = []\nfor file in os.listdir(\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/\"):\n    trains.append(file)\n    \nfor file in os.listdir(\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/test/\"):\n    tests.append(file)\n    \nprint(f\"Number of train files: {len(trains)}\")\nprint(f\"Number of test files: {len(tests)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Feature Extraction**\nCode from [Volcano](https://www.kaggle.com/ymdhryk/volcano/notebook#Feature-Extraction) notebook by [ymdhryk](https://www.kaggle.com/ymdhryk)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define feature extraction function\n\nfs = [\"_mean\",\"_std\",\"_max\",\"_min\",\"_mad\",\"_skew\",\"_kurt\",\"_nunique\",\n      \"_quantile_05\",\"_quantile_10\",\"_quantile_30\",\"_quantile_70\",\"_quantile_90\",\"_quantile_95\",\n      \"_fft_power_mean\",\"_fft_power_std\",\"_fft_power_min\",\"_fft_power_max\",\n      \"_fft_power_sum_low\",\"_fft_power_sum_middle\",\"_fft_power_sum_high\",\n      \"_fft_power_mad\",\"_fft_power_skew\",\"_fft_power_kurt\",\"_fft_power_nunique\",\n      \"_fft_power_quantile_05\",\"_fft_power_quantile_10\",\"_fft_power_quantile_30\",\"_fft_power_quantile_70\",\"_fft_power_quantile_90\",\"_fft_power_quantile_95\",\n      \"_cross_0_count\",\n      \"_roll_mean_min\",\"_roll_mean_max\",\"_roll_dist_min\",\"_roll_dist_max\",\"_roll_dist_diff_min\",\"_roll_dist_diff_max\",\n      #\"_first_005\",\"_last_005\",\"_first_010\",\"_last_010\",\"_first_030\",\"_last_030\",\"_first_070\",\"_last_070\",\"_first_090\",\"_last_090\",\"_first_095\",\"_last_095\",\n      #\"_abs_0250_min\",\"_abs_0250_max\",\"_abs_0500_min\",\"_abs_0500_max\",\"_abs_0750_min\",\"_abs_0750_max\",\"_abs_1000_min\",\"_abs_1000_max\",\"_abs_1250_min\",\"_abs_1250_max\",\"_abs_1500_min\",\"_abs_1500_max\",\n     ]\n\ndef extract(segment_id,dir_=\"train\"):\n    \"\"\"\n    Extract statistical features for each sensor signal\n    \n    - Mean\n    - Standard Deviation\n    - Maximum\n    - Minimum\n    - Mean Absolute Deviation\n    - Skewness\n    - Kurtosis\n    - Median\n    - Mode\n    - (Unbiased) Standard Error of the Mean\n    - Number of Unique Values\n    \"\"\"\n    #display(segment_id)\n    f = pd.read_csv(f\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/{dir_}/{segment_id}.csv\")\n    \n    # Fill NaN\n    f.interpolate(axis=0,inplace=True)\n    #display(f)\n    \n    # Quantile\n    q005 = f.quantile(0.05,axis=0)\n    q010 = f.quantile(0.1 ,axis=0)\n    q030 = f.quantile(0.3 ,axis=0)\n    q070 = f.quantile(0.7 ,axis=0)\n    q090 = f.quantile(0.9 ,axis=0)\n    q095 = f.quantile(0.95,axis=0)\n    \n    # Rolling\n    roll = f.rolling(500)\n    roll_mean = roll.mean()\n    roll_max = roll.max()\n    roll_min = roll.min()\n    roll_dist = roll_max - roll_min\n    roll_dist_diff = roll_dist.diff()\n\n    # FFT power\n    # Remove under flowed 0 frequency and mirrored higher half.\n    fft = pd.DataFrame(np.fft.fft(f.fillna(0)),columns=f.columns).abs().iloc[1:30001,:]\n    fft[f.iloc[1:30001,:].isnull()] = np.nan\n    #display(fft)\n    \n    # Timing information inside 10 minute.\n    f005 = f.where(f < q005)\n    f010 = f.where(f < q010)\n    f030 = f.where(f < q030)\n\n    f070 = f.where(f > q070)\n    f090 = f.where(f > q090)\n    f095 = f.where(f > q095)\n    \n    f_abs = f.abs()\n    f_abs_0250 = f_abs.where(f_abs >  250)\n    f_abs_0500 = f_abs.where(f_abs >  500)\n    f_abs_0750 = f_abs.where(f_abs >  750)\n    f_abs_1000 = f_abs.where(f_abs > 1000)\n    f_abs_1250 = f_abs.where(f_abs > 1250)\n    f_abs_1500 = f_abs.where(f_abs > 1500)\n\n\n    return pd.concat((f.mean(axis=0).add_suffix(\"_mean\"),\n                      f.std(axis=0).add_suffix(\"_std\"),\n                      f.max(axis=0).add_suffix(\"_max\"),\n                      f.min(axis=0).add_suffix(\"_min\"),\n                      f.mad(axis=0).add_suffix(\"_mad\"),\n                      f.skew(axis=0).add_suffix(\"_skew\"),\n                      f.kurt(axis=0).add_suffix(\"_kurt\"),\n                      f.nunique(axis=0).add_suffix(\"_nunique\"),\n                      q005.add_suffix(\"_quantile_05\"),\n                      q010.add_suffix(\"_quantile_10\"),\n                      q030.add_suffix(\"_quantile_30\"),\n                      q070.add_suffix(\"_quantile_70\"),\n                      q090.add_suffix(\"_quantile_90\"),\n                      q095.add_suffix(\"_quantile_95\"),\n                      fft.mean(axis=0).add_suffix(\"_fft_power_mean\"),\n                      fft.std(axis=0).add_suffix(\"_fft_power_std\"),\n                      fft.min(axis=0).add_suffix(\"_fft_power_min\"),\n                      fft.max(axis=0).add_suffix(\"_fft_power_max\"),\n                      fft.iloc[:10000,:].sum(axis=0).add_suffix(\"_fft_power_sum_low\"),\n                      fft.iloc[10000:20000,:].sum(axis=0).add_suffix(\"_fft_power_sum_middle\"),\n                      fft.iloc[20000:,:].sum(axis=0).add_suffix(\"_fft_power_sum_high\"),\n                      fft.mad(axis=0).add_suffix(\"_fft_power_mad\"),\n                      fft.skew(axis=0).add_suffix(\"_fft_power_skew\"),\n                      fft.kurt(axis=0).add_suffix(\"_fft_power_kurt\"),\n                      fft.nunique(axis=0).add_suffix(\"_fft_power_nunique\"),\n                      fft.quantile(0.05,axis=0).add_suffix(\"_fft_power_quantile_05\"),\n                      fft.quantile(0.1,axis=0).add_suffix(\"_fft_power_quantile_10\"),\n                      fft.quantile(0.3,axis=0).add_suffix(\"_fft_power_quantile_30\"),\n                      fft.quantile(0.7,axis=0).add_suffix(\"_fft_power_quantile_70\"),\n                      fft.quantile(0.9,axis=0).add_suffix(\"_fft_power_quantile_90\"),\n                      fft.quantile(0.95,axis=0).add_suffix(\"_fft_power_quantile_95\"),\n                      ((f * f.shift()) < 0).sum(axis=0).add_suffix(\"_cross_0_count\"),\n                      roll_mean.min(axis=0).add_suffix(\"_roll_mean_min\"),\n                      roll_mean.max(axis=0).add_suffix(\"_roll_mean_max\"),\n                      roll_dist.min(axis=0).add_suffix(\"_roll_dist_min\"),\n                      roll_dist.max(axis=0).add_suffix(\"_roll_dist_max\"),\n                      roll_dist_diff.min(axis=0).add_suffix(\"_roll_dist_diff_min\"),\n                      roll_dist_diff.max(axis=0).add_suffix(\"_roll_dist_diff_max\"),\n                      f005.idxmin().add_suffix(\"_first_005\"),\n                      f005.idxmax().add_suffix(\"_last_005\"),\n                      f010.idxmin().add_suffix(\"_first_010\"),\n                      f010.idxmax().add_suffix(\"_last_010\"),\n                      f030.idxmin().add_suffix(\"_first_030\"),\n                      f030.idxmax().add_suffix(\"_last_030\"),\n                      f070.idxmin().add_suffix(\"_first_070\"),\n                      f070.idxmax().add_suffix(\"_last_070\"),\n                      f090.idxmin().add_suffix(\"_first_090\"),\n                      f090.idxmax().add_suffix(\"_last_090\"),\n                      f095.idxmin().add_suffix(\"_first_095\"),\n                      f095.idxmax().add_suffix(\"_last_095\"),\n                      f_abs_0250.idxmin().add_suffix(\"_abs_0250_min\"),\n                      f_abs_0250.idxmax().add_suffix(\"_abs_0250_max\"),\n                      f_abs_0500.idxmin().add_suffix(\"_abs_0500_min\"),\n                      f_abs_0500.idxmax().add_suffix(\"_abs_0500_max\"),\n                      f_abs_0750.idxmin().add_suffix(\"_abs_0750_min\"),\n                      f_abs_0750.idxmax().add_suffix(\"_abs_0750_max\"),\n                      f_abs_1000.idxmin().add_suffix(\"_abs_1000_min\"),\n                      f_abs_1000.idxmax().add_suffix(\"_abs_1000_max\"),\n                      f_abs_1250.idxmin().add_suffix(\"_abs_1250_min\"),\n                      f_abs_1250.idxmax().add_suffix(\"_abs_1250_max\"),\n                      f_abs_1500.idxmin().add_suffix(\"_abs_1500_min\"),\n                      f_abs_1500.idxmax().add_suffix(\"_abs_1500_max\"),\n                     ),\n                     axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Test with small data\n\n%time small_features = sorted_train.iloc[[0,1,2],:][\"segment_id\"].apply(extract)\n\ndisplay(small_features)\n\nframe = small_features.iloc[:0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Extract features for train data\n\ntrain = pd.read_csv(\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train.csv\")\nsorted_train = train.sort_values(\"time_to_eruption\",ascending=False)\n\n%time features = dd.from_pandas(sorted_train[\"segment_id\"],npartitions=4).apply(extract,meta=frame).compute(scheduler=\"processes\")\n\ndata = pd.concat((sorted_train,features),axis=1)\n\n# Save features to reuse\ndata.to_csv(\"train_data.csv\")\n\ndata","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Extract features for test data\n\ntest_ = pd.DataFrame([os.path.basename(f)[:-4] for f in glob.glob('/kaggle/input/predict-volcanic-eruptions-ingv-oe/test/*')], columns=[\"segment_id\"])\n\n%time _test = dd.from_pandas(test_[\"segment_id\"],npartitions=4).apply(extract,dir_=\"test\",meta=frame).compute(scheduler=\"processes\")\ntest = pd.concat((test_,_test),axis=1)\n\n# Save features to reuse\ntest.to_csv(\"test_data.csv\")\n\ntest","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **LightGBM Regressor**\nCode from [Introduction to Volcanology, Seismograms and LGBM](https://www.kaggle.com/jesperdramsch/introduction-to-volcanology-seismograms-and-lgbm#Train-a-LightGBM-Regressor) notebook by [jesperdramsch](https://www.kaggle.com/jesperdramsch)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgbm\nfrom sklearn.model_selection import KFold\nimport gc\n\nn_fold = 7\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=101)\nfeatures = list(data.drop([\"segment_id\", \"time_to_eruption\"], axis=1).columns)\ntarget_name = [\"time_to_eruption\"]\n\nparams = {\n    \"n_estimators\": 2000,\n    \"boosting_type\": \"gbdt\",\n    \"metric\": \"mae\",\n    \"num_leaves\": 66,\n    \"learning_rate\": 0.005,\n    \"feature_fraction\": 0.9,\n    \"bagging_fraction\": 0.8,\n    \"agging_freq\": 3,\n    \"max_bins\": 2048,\n    \"verbose\": 0,\n    \"random_state\": 101,\n    \"nthread\": -1,\n    \"device\": \"gpu\",\n}\n\nsub_preds = np.zeros(test.shape[0])\nfeature_importance = pd.DataFrame(index=list(range(n_fold)), columns=features)\n\nfor n_fold, (trn_idx, val_idx) in enumerate(folds.split(data)):\n    print(f\"Fold {n_fold}:\")\n    trn_x, trn_y = data[features].iloc[trn_idx], data[target_name].iloc[trn_idx]\n    val_x, val_y = data[features].iloc[val_idx], data[target_name].iloc[val_idx]\n    \n    model = lgbm.LGBMRegressor(**params)\n    \n    model.fit(trn_x, trn_y, \n            eval_set= [(trn_x, trn_y), (val_x, val_y)], \n            eval_metric=\"mae\", verbose=200, early_stopping_rounds=50\n           )\n\n    feature_importance.iloc[n_fold, :] = model.feature_importances_\n    \n    sub_preds += model.predict(test[features], num_iteration=model.best_iteration_) / folds.n_splits","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submission file"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame()\nsubmission['segment_id'] = test[\"segment_id\"]\nsubmission['time_to_eruption'] = sub_preds\nsubmission.to_csv('submission.csv', header=True, index=False)\n\nsubmission","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}