{"cells":[{"metadata":{},"cell_type":"markdown","source":"## frequency domain feature from : https://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":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":{"_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\")\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":"print(train.shape)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.isna().sum(axis=0).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train['time_to_eruption'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### time domain feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"# generate feature\n# collect mean / std / 5 / 10 / 20 / 40 percentile / min / max / +5000 / +10000 / +20000 self-corr\ndef 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":{},"cell_type":"markdown","source":"### frequency domain feature"},{"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":"with open('train_part_fea_time_domain.pkl', 'wb') as f:\n    pickle.dump(train_part_fea, f)\n    \nwith open('test_part_fea_time_domain.pkl', 'wb') as f:\n    pickle.dump(test_part_fea, f)","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":"with open('train_part_fea_freq_domain.pkl', 'wb') as f:\n    pickle.dump(train_part_fea_freq, f)\n    \nwith open('test_part_fea_freq_domain.pkl', 'wb') as f:\n    pickle.dump(test_part_fea_freq, f)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"***MODEL***"},{"metadata":{"trusted":true},"cell_type":"code","source":"#train_part_fea = pd.read_pickle('../input/ingv-eda-basemodel/train_part_fea.pkl')\n#test_part_fea = pd.read_pickle('../input/ingv-eda-basemodel/test_part_fea.pkl')","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":"train.to_pickle('train.pkl')\nsample_submission.to_pickle('test.pkl')","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, \n                                                    train['time_to_eruption'].values, \n                                                    test_size=0.25, random_state=42)","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#params = {'objective': 'mae', \n#          'num_iterations': ,\n#          'learning_rate': , \n#          'num_leaves': ,\n#          'seed': ,\n#          'metric': 'mae'}\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":{},"cell_type":"markdown","source":"## predict"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame({'segment_id': sample_submission['segment_id'].values, \n    'time_to_eruption': model.predict(sample_submission.iloc[:, 2:].values)})\nsubmission.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}