{"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\nfrom bayes_opt import BayesianOptimization","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')","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","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)","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')","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":"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)","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":"test_df = test_set.drop('segment_id', axis = 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dtrain = xgb.DMatrix(train_df, label=y_train)\ndtest = xgb.DMatrix(test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def xgb_evaluate(max_depth, gamma, colsample_bytree, min_child_weight, subsample):\n    params = {'eval_metric': 'mae',\n              'max_depth': int(max_depth),\n              'subsample': 0.8,\n              'eta': 0.1,\n              'gamma': gamma,\n              'colsample_bytree': colsample_bytree,\n              'min_child_weight': min_child_weight,\n              'subsample': subsample\n             }\n    # Used around 1000 boosting rounds in the full model\n    cv_result = xgb.cv(params, dtrain, num_boost_round=100, nfold=3)    \n    \n    # Bayesian optimization only knows how to maximize, not minimize, so return the negative RMSE\n    return -1.0 * cv_result['test-mae-mean'].iloc[-1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xgb_bo = BayesianOptimization(xgb_evaluate, {'max_depth': (3, 20), \n                                             'gamma': (0, 1),\n                                             'colsample_bytree': (0.3, 1),\n                                             'min_child_weight': (0.5, 2.5),\n                                             'subsample': (0.3, 1)\n                                            })\n# Use the expected improvement acquisition function to handle negative numbers\n# Optimally needs quite a few more initiation points and number of iterations\nxgb_bo.maximize(init_points=20, n_iter=100, acq='ei')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'colsample_bytree': 0.7074388259225073,\n         'gamma': 0.596567557549089,\n         'max_depth': 11,\n         'min_child_weight': 0.7825572308779551,\n         'subsample': 0.977879394816535}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = xgb.train(params, dtrain)\npred = model.predict(dtest)","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'] = pred","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":"submission.to_csv('out.csv')","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}