{"cells":[{"metadata":{},"cell_type":"markdown","source":"# INGV Volcanic : Basic_solution (STFT)\n**Volcanic earthquakes** have various characteristics depending on the volcano, and there are various classifications, but the typical classifications can be classified as follows.\n \n**Volcanic earthquake**\n * type A (10Hz or higher): Earthquake caused by destruction of rocks surrounding magma chambers and conduits\n * type BH (5-8Hz) : Earthquake caused by magma intruding into the conduit and destroying the conduit and rocks around the conduit\n * type BL (1.5-2.5Hz) : An earthquake around the conduit due to gas etc. ejecting from the crater prior to the explosive eruption and reducing the pressure inside the conduit.\n\n**Volcanic tremor**\n * type C (0.5-1.2Hz) : Vibration due to increase in gas pressure in the cavity along with BH\n * type D (2-4Hz) : Vibration due to gas ejection along with BL\n \n> I am a complete amateur about volcanoes. The jargon may be wrong, but please forgive me."},{"metadata":{},"cell_type":"markdown","source":"**Version 3**\n * At the beginning, briefly add the mechanism of earthquake and tremor\n * Corrected the unit of `time_to_eruption` from millisecond to centisecond (Thanks Alex V B)\n * Changed LightGBM parameters (Thanks [Dave E](https://www.kaggle.com/davidedwards1/volcano-stft-data-optimisation))"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport datetime\nfrom tqdm.notebook import tqdm\n\nimport numpy as np\nimport pandas as pd\nimport scipy\nimport scipy.signal\n\nimport matplotlib\nimport matplotlib.pyplot as plt\npd.options.display.max_columns = None    # disp all columns\n%matplotlib inline\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error as mse\n\n# from lightgbm import LGBMRegressor\nimport lightgbm as lgb","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"### Kaggle or Local-PC ###\nKAGGLE = True       # <==== SET ============\n\nif KAGGLE:\n    DIR = '../input/predict-volcanic-eruptions-ingv-oe'\nelse:              # local PC\n    DIR = './predict-volcanic-eruptions-ingv-oe/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(os.path.join(DIR, 'train.csv'))\ntest = pd.read_csv(os.path.join(DIR, 'sample_submission.csv'))\n\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Convert 'time_to_eruption'to hours:minutes:seconds (Just for reference)\ntrain['h:m:s'] = (train['time_to_eruption']\n                  .apply(lambda x:datetime.timedelta(seconds = x/100)))\ntrain","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Observe sample data"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# plot utility function\ndef plot(ax, x, y, xlabel=None, ylabel=None, legend=None):\n    ax.plot(x, y, label = legend)\n    if xlabel != None:\n        ax.set_xlabel(xlabel)\n    if ylabel != None:\n        ax.set_ylabel(ylabel)\n    if legend != None:\n        ax.legend()\n    ax.grid(True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Serect sample segment\nsample_df = (train.sort_values('time_to_eruption')\n             .reset_index()\n             .rename(columns={'index': 'train_id'}))\nsample_df = sample_df[sample_df.index % (len(train) // 5) == 5].reset_index(drop = True)\nsample_ids = sample_df['segment_id'].values\nsample_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Time Domain"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor = 4      #### 1 ～ 10\n\nfig, ax = plt.subplots(len(sample_ids), 1, figsize = (12, len(sample_ids)*2))\nfor i, segment_id in enumerate(sample_ids):\n    segment_df = pd.read_csv(os.path.join(DIR, f'train/{segment_id}.csv')).fillna(0)\n    ax[i].plot(range(len(segment_df)), segment_df[f'sensor_{sensor}'])\n    ax[i].set_title(f'segment_id : {segment_id},  sensor : {sensor}')\n\nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Time-Frequency Domain (STFT)\nSTFT : Short Time Fourier Transform"},{"metadata":{"trusted":true},"cell_type":"code","source":"fs = 100                # sampling frequency \nN = len(segment_df)     # data size\nn = 256                 # FFT segment size\n\nfig, ax = plt.subplots(len(sample_ids), 1, figsize = (12, len(sample_ids)*2))\nfor i, segment_id in enumerate(sample_ids):\n    segment_df = pd.read_csv(os.path.join(DIR, f'train/{segment_id}.csv')).fillna(0)\n    \n    x = segment_df[f'sensor_{sensor}'][:N]\n    f, t, Z = scipy.signal.stft(x, fs = fs, window = 'hann', nperseg = n)\n    Z = np.abs(Z)\n\n    ax[i].pcolormesh(t, f, Z, vmin = 0, vmax = Z.mean()*10)\n    ax[i].set_ylim(0, 20)\n    ax[i].set_ylabel('Frequency [Hz]'); plt.xlabel('Time [s]')\n    ax[i].set_title(f'segment_id : {segment_id},  sensor : {sensor}')\nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"# STFT(Short Time Fourier Transform) Specifications\nfs = 100                # sampling frequency \nN = len(segment_df)     # data size\nn = 256                 # FFT segment size\nmax_f = 20              # ～20Hz\n\ndelta_f = fs / n        # 0.39Hz\ndelta_t = n / fs / 2    # 1.28s","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_features(tgt):\n    tgt_df = train if tgt == 'train' else test\n    feature_set = []\n    for segment_id in tqdm(tgt_df['segment_id']):\n        segment_df = pd.read_csv(os.path.join(DIR,f'{tgt}/{segment_id}.csv'))\n        segment = [segment_id]\n        for sensor in segment_df.columns:\n            x = segment_df[sensor][:N]\n            if x.isna().sum() > 1000:     ##########\n                segment += ([np.NaN] * 10)\n                continue\n            f, t, Z = scipy.signal.stft(x.fillna(0), fs = fs, window = 'hann', nperseg = n)\n            f = f[:round(max_f/delta_f)+1]\n            Z = np.abs(Z[:round(max_f/delta_f)+1]).T    # ～max_f, row:time,col:freq\n\n            th = Z.mean() * 1     ##########\n            Z_pow = Z.copy()\n            Z_pow[Z < th] = 0\n            Z_num = Z_pow.copy()\n            Z_num[Z >= th] = 1\n\n            Z_pow_sum = Z_pow.sum(axis = 0)\n            Z_num_sum = Z_num.sum(axis = 0)\n\n            A_pow = Z_pow_sum[round(10/delta_f):].sum()\n            A_num = Z_num_sum[round(10/delta_f):].sum()\n            BH_pow = Z_pow_sum[round(5/delta_f):round(8/delta_f)].sum()\n            BH_num = Z_num_sum[round(5/delta_f):round(8/delta_f)].sum()\n            BL_pow = Z_pow_sum[round(1.5/delta_f):round(2.5/delta_f)].sum()\n            BL_num = Z_num_sum[round(1.5/delta_f):round(2.5/delta_f)].sum()\n            C_pow = Z_pow_sum[round(0.6/delta_f):round(1.2/delta_f)].sum()\n            C_num = Z_num_sum[round(0.6/delta_f):round(1.2/delta_f)].sum()\n            D_pow = Z_pow_sum[round(2/delta_f):round(4/delta_f)].sum()\n            D_num = Z_num_sum[round(2/delta_f):round(4/delta_f)].sum()\n            segment += [A_pow, A_num, BH_pow, BH_num, BL_pow, BL_num, C_pow, C_num, D_pow, D_num]\n\n        feature_set.append(segment)\n\n    cols = ['segment_id']\n    for i in range(10):\n        for j in ['A_pow', 'A_num','BH_pow', 'BH_num','BL_pow', 'BL_num','C_pow', 'C_num','D_pow', 'D_num']:\n            cols += [f's{i+1}_{j}']\n    feature_df = pd.DataFrame(feature_set, columns = cols)\n    feature_df['segment_id'] = feature_df['segment_id'].astype('int')\n    return feature_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_df = make_features('train')\ntrain_set = pd.merge(train, feature_df, on = 'segment_id')\ntrain_set","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2, 5, figsize = (12, 6))\nx = train_set['time_to_eruption']\nfor i,type in enumerate(['A_pow','A_num','BH_pow','BH_num','BL_pow','BL_num','C_pow','C_num','D_pow','D_num']):\n    y = np.zeros(len(x))\n    for j in range(10):\n        y += train_set[f's{j+1}_{type}']\n    y /= 10\n    x1 = np.polyfit(x, y.fillna(y.mean()), 2)\n    y1 = np.poly1d(x1)(x)\n    ax[i%2, i//2].plot(x, y,'.')\n    ax[i%2, i//2].plot(x, y1,'.')\n    ax[i%2, i//2].set_ylim(0,)\n    ax[i%2, i//2].set_title(type)\nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Modeling and Predicting"},{"metadata":{"trusted":true},"cell_type":"code","source":"df = train_set.drop(['segment_id', 'time_to_eruption','h:m:s'], axis=1)\ny = train_set['time_to_eruption']\n\nX_train, X_val, y_train, y_val = train_test_split(df, y,\n                                                  random_state = 42,\n                                                  test_size = 0.2,\n                                                  shuffle = True)\n\nfeatures = X_train.columns.tolist()\ncat_features = {}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# lgb = LGBMRegressor(random_state = 42,\n#                     max_depth = 7,\n#                     n_estimators = 250,       ######### \n#                     learning_rate = 0.05)\n# lgb.fit(X_train, y_train)\n# preds = lgb.predict(X_val)\n\n# print('RMSE: ', np.sqrt(mse(y_val, preds)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def do_lgb(X_train, y_train, X_val, y_val):\n    params = {'objective': 'rmse',\n              'metric': 'rmse',\n              'max_depth':14,\n              'min_data_in_leaf':5,         # = min_child_samples\n              'num_leaves': 2**7 - 1,\n              'learning_rate': 0.05,\n              'feature_fraction': 0.7,      # = colsample_bytree\n              'bagging_fraction': 0.5,      # = subsample\n              'bagging_freq': 5,\n              'lambda_l1':80,               # = reg_alpha\n              'num_iterations': 10000,      # = n_estimators\n              'seed': 42,\n              'verbose': 1\n             }\n\n    lgb_train = lgb.Dataset(X_train, y_train)\n    lgb_eval = lgb.Dataset(X_val, y_val, reference=lgb_train)\n\n    evals_result = {}\n    model = lgb.train(\n        params,\n        lgb_train,\n        valid_sets = (lgb_train, lgb_eval), \n        feature_name = features,\n        categorical_feature = cat_features,\n        verbose_eval = 100,\n        evals_result = evals_result,\n        early_stopping_rounds = 200)\n\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lgb_model = do_lgb(X_train, y_train, X_val, y_val)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict test data"},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_df = make_features('test')\ntest_set = pd.merge(test, feature_df, on = 'segment_id')\ntest_set","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Predict test data\npreds = lgb_model.predict(test_set.drop(['segment_id', 'time_to_eruption'], axis=1))\ntest['time_to_eruption'] = preds\ntest[['segment_id','time_to_eruption']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test[['segment_id','time_to_eruption']].to_csv('submission.csv', index=False)","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}