{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 1"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import train_test_split\nimport seaborn as sns\nimport os\nfrom tqdm import tqdm\nimport gc\nimport pickle\nimport scipy\nimport scipy.signal","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train.csv\")\nsample_submission = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def generate_feature_timedomain():\n    \n    def helper(path):\n        data = []\n        for file in tqdm(os.listdir(path)):\n            tmp = []\n            file_path = os.path.join(path, file)\n            d = pd.read_csv(file_path)\n            tmp.append(eval(file[:-4]))\n            # mean\n            tmp += d.mean(axis=0).values.astype('float32').tolist()\n            # std\n            tmp += d.std(axis=0).values.astype('float32').tolist()\n            # min\n            tmp += d.min(axis=0).values.astype('float32').tolist()\n            # max\n            tmp += d.max(axis=0).values.astype('float32').tolist()\n            # 5 percentile\n            tmp += d.quantile(0.05, axis=0).values.astype('float32').tolist()\n            # 10 percentile\n            tmp += d.quantile(0.1, axis=0).values.astype('float32').tolist()\n            # 20 percentile\n            tmp += d.quantile(0.2, axis=0).values.astype('float32').tolist()\n            # 40 percentile\n            tmp += d.quantile(0.4, axis=0).values.astype('float32').tolist()\n            # 60 percentile\n            tmp += d.quantile(0.6, axis=0).values.astype('float32').tolist()\n            # 80 percentile\n            tmp += d.quantile(0.8, axis=0).values.astype('float32').tolist()\n            # shift\n            for col in d:\n                d[col+'_5000'] = d[col].shift(5000)\n                d[col+'_10000'] = d[col].shift(10000)\n                d[col+'_20000'] = d[col].shift(20000)\n                d[col+'_30000'] = d[col].shift(30000)\n                \n            # +5000 / +10000 / +20000 / +30000 self-corr\n            for col in d.columns[:10]:\n                col1 = col+'_5000'\n                col2 = col+'_10000'\n                col3 = col+'_20000'\n                col4 = col+'_30000'\n                tmp1 = d.loc[:, [col, col1]].dropna()\n                tmp2 = d.loc[:, [col, col2]].dropna()\n                tmp3 = d.loc[:, [col, col3]].dropna()\n                tmp4 = d.loc[:, [col, col4]].dropna()\n                tmp += [tmp1[col].corr(tmp1[col1]), \n                        tmp2[col].corr(tmp2[col2]), \n                        tmp3[col].corr(tmp3[col3]),\n                        tmp4[col].corr(tmp4[col4])]\n                \n            data.append(tmp)\n        return data\n                   \n    print('train_part: ')\n    train_part_fea = helper('../input/predict-volcanic-eruptions-ingv-oe/train')\n    print('test_part: ')\n    test_part_fea = helper('../input/predict-volcanic-eruptions-ingv-oe/test')\n    \n    return train_part_fea, test_part_fea","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def generate_feature_freq_domain():\n    # STFT\n    fs = 100\n    n = 256\n    N = 60001\n    max_f = 20\n    delta_f = fs / n\n    delta_t = n / fs / 2\n    \n    def helper(fs, n, N, max_f, delta_f, path):\n        data = []\n        for file in tqdm(os.listdir(path)):\n            tmp = []\n            file_path = os.path.join(path, file)\n            d = pd.read_csv(file_path)\n            tmp.append(eval(file[:-4]))\n            \n            for i in range(d.shape[1]):\n                if d.iloc[:, i].isna().sum() > 1000:\n                    tmp += [np.nan] * 7 * 65\n                    tmp += [np.nan] * 10\n                else:\n                    # STFT\n                    f, t, Z = scipy.signal.stft(d.iloc[:, i].fillna(0).values, fs = fs, window = 'hann', nperseg = n)\n                    f = f[:round(max_f/delta_f)+1]\n                    \n                    Z_half = np.abs(Z[:round(Z.shape[0]//2)+1]).T\n                    tmp += Z_half.min(axis=0).astype('float32').tolist()\n                    tmp += Z_half.max(axis=0).astype('float32').tolist()\n                    tmp += Z_half.std(axis=0).astype('float32').tolist()\n                    tmp += Z_half.mean(axis=0).astype('float32').tolist()\n                    tmp += np.quantile(Z_half, 0.25, axis=0).astype('float32').tolist()\n                    tmp += np.quantile(Z_half, 0.5, axis=0).astype('float32').tolist()\n                    tmp += np.quantile(Z_half, 0.75, axis=0).astype('float32').tolist()\n                    \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                    tmp += [A_pow, A_num, BH_pow, BH_num, BL_pow, BL_num, C_pow, C_num, D_pow, D_num]\n            data.append(tmp)\n        return data\n    \n    print('train_part: ')\n    train_part_fea = helper(fs, n, N, max_f, delta_f, path='../input/predict-volcanic-eruptions-ingv-oe/train')\n    print('test_part: ')\n    test_part_fea = helper(fs, n, N, max_f, delta_f, path='../input/predict-volcanic-eruptions-ingv-oe/test')\n    \n    return train_part_fea, test_part_fea","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_part_fea, test_part_fea = generate_feature_timedomain()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_part_fea_freq, test_part_fea_freq = generate_feature_freq_domain()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base_colname = ['sensor_'+str(i) for i in range(1, 11)]\nfea_colname = ['segment_id'] + [j + '_mean' for j in base_colname] + [j + '_std' for j in base_colname] + \\\n                [j + '_min' for j in base_colname] + [j + '_max' for j in base_colname] + \\\n                    [j + '_5_quant' for j in base_colname] + [j + '_10_quant' for j in base_colname] + \\\n                        [j + '_20_quant' for j in base_colname] + [j + '_40_quant' for j in base_colname] + \\\n                        [j + '_60_quant' for j in base_colname] + [j + '_80_quant' for j in base_colname] + \\\n                    [j + i for j in base_colname for i in ['_5000_self_corr', '_10000_self_corr', \n                                                           '_20000_self_corr', '_30000_self_corr']]\n\ntrain = pd.merge(train, pd.DataFrame(train_part_fea, columns=fea_colname), on='segment_id', how='left')\nsample_submission = pd.merge(sample_submission, pd.DataFrame(test_part_fea, columns=fea_colname), on='segment_id', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fea_freq_colname = ['segment_id']\nfor i in base_colname:\n    for j in range(65):\n        for s in ['min','max', 'std', 'mean', '25_quant', '50_quant', '75_quant']:\n            fea_freq_colname.append(i+'_freq'+str(j)+'_'+s)\n    fea_freq_colname.extend([i + ss for ss in ['_A_pow', '_A_num', '_BH_pow', '_BH_num', '_BL_pow', \n                                               '_BL_num', '_C_pow', '_C_num', '_D_pow', '_D_num']])\n\ntrain = pd.merge(train, pd.DataFrame(train_part_fea_freq, columns=fea_freq_colname), on='segment_id', how='left')\nsample_submission = pd.merge(sample_submission, pd.DataFrame(test_part_fea_freq, \n                                                             columns=fea_freq_colname), on='segment_id', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_val, y_train, y_val = train_test_split(train.drop(['segment_id', 'time_to_eruption'], axis=1).values, train['time_to_eruption'].values, test_size=0.25, random_state=28)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\n\ntrain_data = lgb.Dataset(X_train, label=y_train)\nval_data = lgb.Dataset(X_val, y_val, reference=train_data,)\n\n\nparams = { 'num_leaves': 85,\n          'n_estimators': 1000,\n    'min_data_in_leaf': 10, \n    'objective':'mae',\n    'max_depth': -1,\n    'learning_rate': 0.01,\n    'max_bins': 2048,\n    \"boosting\": \"gbdt\",\n    \"feature_fraction\": 0.91,\n    \"bagging_freq\": 1,\n    \"bagging_fraction\": 0.91,\n    \"bagging_seed\": 42,\n    \"metric\": 'mae',\n    \"lambda_l1\": 0.1,\n    \"verbosity\": -1,\n    \"nthread\": -1,\n    \"random_state\": 42}\n\nmodel = lgb.train(params=params, train_set=train_data, valid_sets=[train_data, val_data], valid_names=['train', 'val'], \n                  early_stopping_rounds=50)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred_1 = model.predict(sample_submission.iloc[:, 2:].values)\ny_pred_1 = [x if x>=0 else 0 for x in y_pred_1]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport gc\nimport numpy as np\nimport pandas as pd\n\nfrom time import time\nfrom time import ctime\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nfrom tqdm import tqdm_notebook\nfrom tqdm import tqdm\n\nimport joblib\nfrom joblib import Parallel, delayed\nimport multiprocessing\nnum_cores = multiprocessing.cpu_count()-1\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import r2_score\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.model_selection import KFold\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time\n\nfrom sklearn.model_selection import KFold\n\nX = train.drop(['segment_id', 'time_to_eruption'], axis=1)\nY = train['time_to_eruption']\nX_test = sample_submission.iloc[:, 2:]\n\nn_fold = 5\ncv = KFold(n_splits=n_fold, shuffle=True, random_state=14)\n\noof = np.zeros(len(X))\ncat_prediction = np.zeros(len(X_test))\nmae, r2 = [], []\n\nPARAMS = {\n    \n             'random_seed': 42,\n             'eval_metric': 'MAE'\n\n        }\n\nfor fold_n, (train_index, valid_index) in enumerate(cv.split(X)):\n\n    X_train = X.iloc[train_index,:]\n    X_valid = X.iloc[valid_index,:]\n    \n    Y_train = Y.iloc[train_index]\n    Y_valid = Y.iloc[valid_index]\n          \n    best_model = CatBoostRegressor(**PARAMS, thread_count = -1)  \n    \n    train_dataset = Pool(data=X_train,\n                     label=Y_train,\n                     )\n    \n    eval_dataset = Pool(data=X_valid,\n                    label=Y_valid,\n                    )\n    \n    best_model.fit(train_dataset,\n              use_best_model=True,\n              verbose = False,\n              plot = True,\n              eval_set=eval_dataset,\n              early_stopping_rounds=100)\n\n   \n    y_pred = best_model.predict(Pool(data=X_valid))\n\n    mae.append(mean_absolute_error(Y_valid, y_pred))\n    r2.append(r2_score(Y_valid, y_pred))\n\n    print('MAE: ', mean_absolute_error(Y_valid, y_pred))\n    print('R2: ', r2_score(Y_valid, y_pred))\n\n    cat_prediction += best_model.predict(Pool(data=X_test))\n        \ncat_prediction /= n_fold","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred_2 = cat_prediction\ny_pred_2 = [x if x>=0 else 0 for x in y_pred_2]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3"},{"metadata":{"trusted":true},"cell_type":"code","source":"import 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":{"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":{"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'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['h:m:s'] = (train['time_to_eruption']\n                  .apply(lambda x:datetime.timedelta(seconds = x/100)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_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":{"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":{"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":{"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')","execution_count":null,"outputs":[]},{"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":"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":{"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":"preds = lgb_model.predict(test_set.drop(['segment_id', 'time_to_eruption'], axis=1))\ny_pred_3 = preds\ny_pred_3 = [x if x>=0 else 0 for x in y_pred_3]\ny_pred_3","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# FINAL"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame()\nsubmission['segment_id'] = test['segment_id']\nsubmission['time_to_eruption'] = [(y_pred_1[i] + y_pred_2[i] + y_pred_3[i])/3 for i in range(len(y_pred_1))]\nsubmission.to_csv('submission_recent.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission","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}