{"cells":[{"metadata":{},"cell_type":"markdown","source":"In this notebook, I am using the global statistical values (mean, max, min etc), spectral density of the sensor signals and shorttime fourier transforms as the features. \nI am using the principal components (first 21 components that explain upto 99.9% of the variance)of the spectral density features to reduce overfitting. \nFinally I am training an XGB model with five fold validation(Ordering the data in increasing order of time to eruption so as to ensure a similar distribution in each of the folds)"},{"metadata":{},"cell_type":"markdown","source":"I referred to the following notebooks to obtain ideas and some code.Thanks to the respective authors.\n1. https://www.kaggle.com/soheild91/ingv-nn-xgb-baseline (Getting started)\n2. https://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft (STFT features)\n3. https://www.kaggle.com/kylesnyder/ingv-spectral-density-w-randomforest (Spectral density features)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport scipy\nfrom scipy import signal\nimport random\nimport pickle\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.decomposition import PCA\n\nimport xgboost","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Spectral density features"},{"metadata":{"trusted":true},"cell_type":"code","source":"# #spectral density along with statistical features\n# input_df = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\n# # data = np.empty((input_df.shape[0],1410))#129 features per sensor 10 sensors =1290 + 120 stat ftrs(12*10)\n# data = np.empty((input_df.shape[0],1290))\n# time = np.empty((input_df.shape[0],1))\n\n# for i,segment in enumerate(input_df['segment_id']):\n#     temp = pd.read_csv(f'../input/predict-volcanic-eruptions-ingv-oe/train/{segment}.csv').fillna(0)\n#     temp_arr = np.empty((0,))\n#     for col in temp.columns:\n#         freq,psd =signal.welch(temp[col],100)\n#         temp_arr= np.concatenate((temp_arr,psd))  \n# #     temp_arr = np.concatenate((temp_arr,temp.abs().mean().to_numpy(),\n# #                                 temp.std().to_numpy(),\n# #                                 temp.mean().to_numpy(),\n# #                                 temp.var().to_numpy(),\n# #                                 temp.min().to_numpy(),\n# #                                 temp.max().to_numpy(),\n# #                                 temp.median().to_numpy(),\n# #                                 temp.quantile([0.1,0.25,0.5,0.75,0.9]).to_numpy().reshape(1,-1)[0]))\n#     temp_arr = temp_arr.reshape((1,-1))\n#     data[i,:] = temp_arr\n#     time[i,0] = input_df.loc[i,'time_to_eruption']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# sample_submission_df=pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv',nrows=5)\n# #data_test=np.empty((sample_submission_df.shape[0],1410))\n# data_test=np.empty((sample_submission_df.shape[0],1290))\n# for i,segment in enumerate(sample_submission_df['segment_id']):\n#     temp=pd.read_csv(f'../input/predict-volcanic-eruptions-ingv-oe/test/{segment}.csv').fillna(0)\n#     temp_arr = np.empty((0,))\n#     for col in temp.columns:\n#         freq,psd =signal.welch(temp[col],100)\n#         temp_arr= np.concatenate((temp_arr,psd))\n# #     temp_arr = np.concatenate((temp_arr,temp.abs().mean().to_numpy(),\n# #                                 temp.std().to_numpy(),\n# #                                 temp.mean().to_numpy(),\n# #                                 temp.var().to_numpy(),\n# #                                 temp.min().to_numpy(),\n# #                                 temp.max().to_numpy(),\n# #                                 temp.median().to_numpy(),\n# #                                 temp.quantile([0.1,0.25,0.5,0.75,0.9]).to_numpy().reshape(1,-1)[0]))\n#     temp_arr = temp_arr.reshape((1,-1))\n#     data_test[i,:] = temp_arr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train_df_sd = pd.DataFrame(data)\n# test_df_sd = pd.DataFrame(data_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# with open('train_df_sd.pkl', 'wb') as f:\n#     pickle.dump(train_df_sd, f)\n    \n# with open('test_df_sd.pkl', 'wb') as f:\n#     pickle.dump(test_df_sd, f)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# STFT Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"# # STFT(Short Time Fourier Transform) Specifications\n# fs = 100                # sampling frequency \n# #N = len(segment_df)     # data size\n# N = 60000\n# n = 256                 # FFT segment size\n# max_f = 20              # ～20Hz\n\n# delta_f = fs / n        # 0.39Hz\n# delta_t = n / fs / 2    # 1.28s","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# #path ../input/predict-volcanic-eruptions-ingv-oe/test\n# def make_features(input_df,path):\n#     feature_set = []\n#     for i,segment_id in enumerate(input_df['segment_id']):\n#         temp=pd.read_csv(f'{path}/{segment_id}.csv')#.fillna(0)\n#         segment = [segment_id]\n#         for sensor in temp.columns:\n#             x = temp[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#     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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# input_df = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\n# path_train = \"../input/predict-volcanic-eruptions-ingv-oe/train\"\n# train_df_stft = make_features(input_df,path_train)\n# train_df_stft = pd.merge(input_df, train_df_stft, on = 'segment_id')\n# #test\n# sample_submission_df=pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv')\n# path_test = \"../input/predict-volcanic-eruptions-ingv-oe/test\"\n# test_df_stft = make_features(sample_submission_df,path_test)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# with open('train_df_stft.pkl', 'wb') as f:\n#     pickle.dump(train_df_stft, f)\n    \n# with open('test_df_stft.pkl', 'wb') as f:\n#     pickle.dump(test_df_stft, f)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df_stft = pd.read_pickle('../input/volcano-features/train_df_stft.pkl')\ntime = train_df_stft['time_to_eruption']\ntrain_df_stft = train_df_stft.drop(['segment_id', 'time_to_eruption'], axis=1)\n\n# train_df_sd = pd.read_pickle('../input/volcano-features/train_df_sd_only.pkl')\ntrain_df_sd = pd.read_pickle('../input/volcano-features/train_df_sd.pkl')\n\ntest_df_stft = pd.read_pickle('../input/volcano-features/test_df_stft.pkl')\ntest_df_stft = test_df_stft.drop(['segment_id'],axis = 1)\n\n# test_df_sd = pd.read_pickle('../input/volcano-features/test_df_sd_only.pkl')\ntest_df_sd = pd.read_pickle('../input/volcano-features/test_df_sd.pkl')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PCA"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.decomposition import PCA\npca = PCA(n_components=0.999,svd_solver=\"full\") #all components that explain upto 99.9% of the variance\ntrain_pca = pca.fit_transform(train_df_sd)\ntest_pca = pca.transform(test_df_sd)\ntrain_pca = pd.DataFrame(train_pca)\ntest_pca = pd.DataFrame(test_pca)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.concat([train_df_stft,train_pca] ,axis=1)\ntest_df = pd.concat([test_df_stft,test_pca],axis = 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# XGB with stratified k fold "},{"metadata":{"trusted":true},"cell_type":"code","source":"input_df = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\ninput_df = input_df.sort_values(\"time_to_eruption\")\nfold_list = [1,2,3,4,5]\nfolds = []\nfor i in range(int((input_df.shape[0]-1)/5)):\n    random.shuffle(fold_list)\n    folds.extend(fold_list)\nfolds = folds + [1] #adding a remaining solitary record to fold 1\ninput_df['fold'] = folds\n#input_df.head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = np.zeros(len(test_df))\nfor fold in range(1,6):\n    train_index_list = input_df[input_df['fold'] != fold].index\n    test_index_list = input_df[input_df['fold'] == fold].index\n\n    X_train = train_df.iloc[train_index_list]\n    y_train = time[train_index_list]\n    X_val = train_df.iloc[test_index_list]\n    y_val = time[test_index_list]\n\n    model = xgboost.XGBRegressor(n_estimators=100000,tree_method='gpu_hist',max_depth=8,learning_rate=0.05,alpha=0.1,SUBSAMPLE=0.6)#,colsample_bytree=0.5)\n    eval_set = [(X_val, y_val)]\n    model.fit(X_train, y_train,early_stopping_rounds=5,eval_metric='mae', eval_set=eval_set, verbose=False)\n    print(model.evals_result()['validation_0']['mae'][-5:])\n    predictions += model.predict(test_df)\npredictions = predictions/5\n\nsample_submission_df1=pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv')\nsample_submission_df1['time_to_eruption']=predictions\nsample_submission_df1.to_csv('xgb_5fldst_ft_sdpca_stft.csv',index=False)","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}