{"cells":[{"metadata":{},"cell_type":"markdown","source":"  \n **stft code copied from** [https://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft](https://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft) "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nfrom scipy.stats import skew, kurtosis\nfrom scipy.signal import stft\nimport os\nfrom tqdm.notebook import tqdm\nimport scipy\nimport re\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error as mse\nimport matplotlib.pyplot as plt\nfrom xgboost import plot_importance\nimport xgboost as xgb\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for dirname, _, filenames_train in os.walk('/kaggle/input/predict-volcanic-eruptions-ingv-oe/train'): \n    continue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for dirname, _, filenames_test in os.walk('/kaggle/input/predict-volcanic-eruptions-ingv-oe/test'): \n    continue","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**COPIED**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# STFT(Short Time Fourier Transform) Specifications\nfs = 100                # sampling frequency \nN = 60001     # 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\n\nDIR = '../input/predict-volcanic-eruptions-ingv-oe'","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')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_set = pd.merge(train, feature_df, on = 'segment_id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_columns(df,result):\n    df = df.fillna(0)\n    for column in df:\n        result.at[index,'sum_'+column] = sum(df[column])\n        result.at[index,'med_'+column] = np.median(df[column])\n        result.at[index,'permiss_'+column] =  df[column].isnull().sum() / df[column].size\n        result.at[index, 'skew_'+ column] = skew(df[column])\n        result.at[index, 'kurtosis'+ column] = kurtosis(df[column])\n        result.at[index,'max_'+column] = df[column].max()\n        result.at[index,'min_'+column] = df[column].min()\n        result.at[index, 'std_'+ column] = np.std(df[column])\n        result.at[index, 'var_'+ column] = np.var(df[column])\n        result.at[index,'quan_0.05'+ column] = np.quantile(df[column],0.05)\n        result.at[index,'quan_0.1'+ column] = np.quantile(df[column],0.1)\n        result.at[index,'quan_0.15'+ column] = np.quantile(df[column],0.15)\n        result.at[index,'quan_0.2'+ column] = np.quantile(df[column],0.2)\n        result.at[index,'quan_0.025'+ column] = np.quantile(df[column],0.25)\n        result.at[index,'quan_0.3'+ column] = np.quantile(df[column],0.3)\n        result.at[index,'quan_0.35'+ column] = np.quantile(df[column],0.35)\n        result.at[index,'quan_0.4'+ column] = np.quantile(df[column],0.4)\n        result.at[index,'quan_0.45'+ column] = np.quantile(df[column],0.45)\n        result.at[index,'quan_0.5'+ column] = np.quantile(df[column],0.5)\n        result.at[index,'quan_0.55'+ column] = np.quantile(df[column],0.55)\n        result.at[index,'quan_0.6'+ column] = np.quantile(df[column],0.6)\n        result.at[index,'quan_0.65'+ column] = np.quantile(df[column],0.65)\n        result.at[index,'quan_0.7'+ column] = np.quantile(df[column],0.7)\n        result.at[index,'quan_0.75'+ column] = np.quantile(df[column],0.75)\n        result.at[index,'quan_0.8'+ column] = np.quantile(df[column],0.8)\n        result.at[index,'quan_0.85'+ column] = np.quantile(df[column],0.85)\n        result.at[index,'quan_0.9'+ column] = np.quantile(df[column],0.9)\n        result.at[index,'quan_0.95'+ column] = np.quantile(df[column],0.95)\n        \n    return result\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for row in  train_set.itertuples():\n    index = row[0]\n    segmentid = row[1]\n    if str(segmentid)+\".csv\" in filenames_train:\n    \n        df_segement = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/'+str(segmentid)+\".csv\")\n       \n        result = pd.concat([df_segement, df_segement.abs().add_suffix(\"_abs\")], axis=1, join=\"inner\")\n        df_segemnt = result\n        train_set = create_columns(df_segement,train_set)\n    if index % 100 == 0:\n        print(index)\n           \n            \n        \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.DataFrame([(re.findall(\"[0-9]+\",a)) for a in filenames_test] , columns=['segment_id'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_df_test = make_features('test')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_df_test.drop('segment_id', axis = 1, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_set = test.join(feature_df_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_set","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for row in test_set.itertuples():\n    index = row[0]\n    segementId = row[1]\n    \n    if str(segementId)+\".csv\" in filenames_test:\n        df_segement = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/test/'+str(segementId)+\".csv\")\n       \n        result = pd.concat([df_segement, df_segement.abs().add_suffix(\"_abs\")], axis=1, join=\"inner\")\n        df_segemnt = result\n        test_set = create_columns(df_segement,test_set)\n    if index % 100 == 0:\n        print(index)\n           ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_set","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_set","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_set.to_csv('test_work.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_set.to_csv('train_work.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train = train_set['time_to_eruption']\ntrain_df = train_set.drop(['time_to_eruption','segment_id'], axis = 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def report(results, n_top=3):\n    for i in range(1, n_top + 1):\n        candidates = np.flatnonzero(results['rank_test_score'] == i)\n        for candidate in candidates:\n            print(\"Model with rank: {0}\".format(i))\n            print(\"Mean validation score: {0:.3f} (std: {1:.3f})\"\n                  .format(results['mean_test_score'][candidate],\n                          results['std_test_score'][candidate]))\n            print(\"Parameters: {0}\".format(results['params'][candidate]))\n            print(\"\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nfrom sklearn.model_selection import RandomizedSearchCV\nimport random\nparam = {'eta': [0.05,0.1,0.2,0.3],\n        'max_depth': [4,5,6,7,8,9,10],\n         'subsample ': [0.5,0.75,1],\n         'gamma': [0.05,0.075,0.09,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]\n         \n        }\nsearch = RandomizedSearchCV(xgb.XGBRegressor(),param,n_iter = 30,cv = 5, scoring = 'neg_mean_squared_error')\nsearch.fit(train_df,y_train)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"report(search.cv_results_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'subsample': 0.75, 'max_depth': 9, 'gamma': 0.1, 'eta': 0.1}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xgb_reg = xgb.XGBRegressor(**params)\nxgb_reg.fit(train_df, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sorted_idx = xgb_reg.feature_importances_.argsort()\nimp = xgb_reg.feature_importances_\ncol = train_df.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.barh(col[sorted_idx[:50]], imp[sorted_idx[:50]])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"todrop = col[sorted_idx[:50]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.drop(todrop, axis = 1, inplace = True)\ntest_set.drop(todrop, axis = 1, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xgb_reg = xgb.XGBRegressor(**params)\nxgb_reg.fit(train_df, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = test_set.drop('segment_id', axis = 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame(test_set['segment_id'], columns=['segment_id'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission['time_to_eruption'] = xgb_reg.predict(test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.set_index('segment_id', inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nsubmission.to_csv('out.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}