{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"scrolled":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 dask.dataframe as dd\n\nimport matplotlib.pyplot as plt\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler\n\n\n#import lightgbm as lgb\nfrom optuna.integration import lightgbm as lgb\n\nimport glob\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":{},"cell_type":"markdown","source":"# Read Data"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"scrolled":false},"cell_type":"code","source":"# meta data\ntrain = pd.read_csv(\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train.csv\")\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Sort data by time to eruption\nsorted_train = train.sort_values(\"time_to_eruption\",ascending=False)\nsorted_train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# meta data for test\ntest_ = pd.DataFrame([os.path.basename(f)[:-4] for f in glob.glob('/kaggle/input/predict-volcanic-eruptions-ingv-oe/test/*')], columns=[\"segment_id\"])\ntest_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# EDA"},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# See time to eruption distribution -> Roughly uniform distribution\n\nsorted_train[[\"time_to_eruption\"]].plot(kind=\"hist\",bins=100,figsize=(10,7))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Max and Min data\ndisplay(sorted_train.iloc[[0,-1],:])\n\nmax_id = 1923243961\nmin_id =  601524801","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Plot the first (probably normal) and the last (probably abnormal) raw signal data -> Looks different\n\npd.read_csv(f\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/{max_id}.csv\").plot(figsize=(15,10),title=\"max time_to_eruption (probably normal)\",subplots=True,ylim=(-10000,10000))\nplt.show()\n\n\npd.read_csv(f\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/{min_id}.csv\").plot(figsize=(15,10),title=\"min time_to_eruption (probably abnormal)\",subplots=True,ylim=(-10000,10000))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Check NaN data -> Some Sensors are completely empty.\n\ndisplay(pd.read_csv(f\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/{max_id}.csv\").isnull().sum())\n\ndisplay(pd.read_csv(f\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/{min_id}.csv\").isnull().sum())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature Extraction"},{"metadata":{"trusted":true,"scrolled":false},"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,"scrolled":false},"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,"scrolled":false},"cell_type":"code","source":"#Extract features for train data\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 resuse\ndata.to_csv(\"train_data.csv\")\n\ndata","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Extract features for test data\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":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# Check Features\n\nfor _fs in fs:\n    data.plot(x=\"time_to_eruption\",\n              y=[f\"sensor_{i}\" + _fs for i in range(1,11)],\n              marker=\".\",linestyle=\"\",figsize=(15,15),subplots=True)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# Check Features\n\n\nfor _fs in fs:\n    for i in range(1,11):\n        _c = f\"sensor_{i}{_fs}\"\n        plot_data = data[_c].sort_values()\n        _nd = plot_data.notna().sum()\n        plt.plot(plot_data,np.arange(1,plot_data.shape[0]+1)/_nd,\n                 color=\"tab:blue\",label=\"train data\",marker=\".\",linestyle=\":\",alpha=0.5)\n\n        plot_test = test[_c].sort_values()\n        _nt = plot_test.notna().sum()\n        plt.plot(plot_test,np.arange(1,plot_test.shape[0]+1)/_nt,\n                 color=\"tab:red\" ,label=\"test data\" ,marker=\".\",linestyle=\":\",alpha=0.5)\n\n        plt.title(_c)\n        plt.legend()\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train"},{"metadata":{"trusted":true},"cell_type":"code","source":"scale_X = StandardScaler()\nscale = StandardScaler()\n\n\ntrain_X, test_X, train_y, test_y = train_test_split(data.drop(columns=[\"segment_id\",\"time_to_eruption\"]),\n                                                    data[[\"time_to_eruption\"]],\n                                                    test_size=0.05)\n\n\ncols = [f\"sensor_{i}{_fs}\" for i in range(1,11) for _fs in fs]\n\n\n# Scale X\nscaled_train_X = pd.DataFrame(scale_X.fit_transform(train_X[cols]),index=train_X.index,columns=cols)\nscaled_test_X  = pd.DataFrame(scale_X.transform(test_X[cols])     ,index=test_X.index ,columns=cols)\nscaled_test    = pd.DataFrame(scale_X.transform(test[cols])       ,index=test.index   ,columns=cols)\n\n# Scale y\nscaled_train_y = scale.fit_transform(train_y)[:,0]\nscaled_test_y  = scale.transform(test_y)[:,0]\n\n\n# LightGBM parameters\nparams = {'task' : 'train',\n          'boosting_type' : 'gbdt',\n          'objective' : 'regression',\n          'metric' : 'mae',\n          'verbose' : 0}\n\n    \n# Create GBM\n# create dataset for lightgbm\nlgb_train = lgb.Dataset(scaled_train_X, scaled_train_y)\nlgb_eval  = lgb.Dataset(scaled_test_X , scaled_test_y , reference=lgb_train)\n\n# Train LightGBM\ngbm = lgb.train(params,\n                lgb_train,\n                num_boost_round=1000,\n                valid_sets=[lgb_train, lgb_eval],\n                valid_names=[\"train\",\"test\"],\n                early_stopping_rounds=30)\nprint(f\"Best Iteration: {gbm.best_iteration}\")\nprint(f\"Best Score: {gbm.best_score}\")\n\n# Plot Training Results\nplt.figure(figsize=(12,7))\nplt.plot(scaled_train_y,gbm.predict(scaled_train_X),marker=\".\",linestyle=\"\",color=\"tab:blue\",label=\"train data\")\nplt.plot(scaled_test_y ,gbm.predict(scaled_test_X) ,marker=\".\",linestyle=\"\",color=\"tab:red\" ,label=\"test data\")\nplt.plot(np.arange(-2.0,2.0,0.1),np.arange(-2.0,2.0,0.1),color=\"tab:green\")\nplt.legend()\nplt.xlabel(\"(Scaled) True time_to_eruption\")\nplt.ylabel(\"(Scaled) Pred time_to_eruption\")\nplt.show()\n\n# Plot Feature Importance\nlgb.plot_importance(gbm,figsize=(10,70))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create submission file"},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"submit = pd.DataFrame(scale.inverse_transform(gbm.predict(scaled_test)),\n                      index=test[\"segment_id\"],\n                      columns=[\"time_to_eruption\"],\n                      dtype=\"int\")\nsubmit.clip(lower=0,inplace=True)\n\ndisplay(submit)\n\nsubmit.to_csv(\"submission.csv\")\n\n!cat submission.csv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}