{"cells":[{"metadata":{},"cell_type":"markdown","source":"# INGV - Volcano prediction\n\n![](https://media.metrolatam.com/2018/01/23/mordorsauron-1200x600.jpg)\n### **Motivation**\n#### 10,000 Orcs now stand between Frodo and Mount Doom. There's still hope for Frodo. He needs time, and safe passage across the Plains of Gorgoroth. He needs to know when Mount Doom will erupt.\n\n#### In this notebook I tried to tackle the problem of volcano eruption prediction via converting sensor readings to images using techniques like Gramian Angular Field and CNNs. A Gramian Angular Field is an image obtained from a time series, representing some temporal correlation between each time point ([source](https://pyts.readthedocs.io/en/stable/auto_examples/image/plot_gaf.html)):\n\n\n![](https://pyts.readthedocs.io/en/stable/_images/sphx_glr_plot_gaf_001.png)\n\n\n#### A similar approach was used [here](https://ieeexplore.ieee.org/document/8700425). I wasn't able to find the right button to make it work (bad approach? bad preprocessing? not enough data? not enough expertise with CNNs?) and I don't want to dedicate more time to this either. Maybe somebody finds it interesting and is able to achieve good results. In that case please let me know (and maybe upvote :D)! Also any feedback is welcome."},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install pyts","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# IMPORTINGS\nimport os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nimport gc\ngc.enable()\n\n# Sklearn\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_absolute_error\n\nfrom pyts.image import GramianAngularField\n\n# Keras\nfrom keras import applications\nfrom keras.models import Sequential\nfrom keras.layers import Conv2D, BatchNormalization, MaxPooling2D, Dropout, Flatten, Dense\nfrom keras.optimizers import Adam","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# GLOBAL VARIABLES\nPATH = '/kaggle/input/predict-volcanic-eruptions-ingv-oe'\nN_SPLITS = 5\nIMG_SIZE = 28\nCHANNELS = 10\nLR = 1e-5\nEPOCHS = 30\nBS = 64\nDROP = 0.5\nSEED = 42\n\n\n# HELPING FUNCTIONS\ndef tseries_to_img(df_list, type_d, method='difference'):\n    \"\"\"\n    Function that transforms segments time series into images using GramianAngularField\n    \n    Segments are preprocessed as follows:\n    - Missing data in channels are filled with the mean\n    - If a channel is completely missing --> fill with zeroes\n    \n    Once the data of all sensors are converted to an image, all images are stacked horizontally\n    \n    \"\"\"\n    \n    df_list.index = df_list['segment_id']\n    df_img = np.zeros((df_list.shape[0], IMG_SIZE, IMG_SIZE*CHANNELS, 1))\n\n    \n    for r,seg in enumerate(tqdm(df_list['segment_id'].values.tolist())):\n        seg_df = pd.read_csv(os.path.join(PATH, type_d, str(seg)+'.csv'))\n        \n        for r2,sens in enumerate(range(CHANNELS)):\n            seg_sens = seg_df['sensor_'+str(sens+1)].values.reshape(1,-1)\n            \n            if np.isnan(seg_sens).sum() < seg_sens.shape[1]:\n                seg_sens[np.isnan(seg_sens)] = np.mean(seg_sens[~np.isnan(seg_sens)])\n            else:\n                seg_sens[np.isnan(seg_sens)] = 0\n                \n            gadf = GramianAngularField(image_size=IMG_SIZE, method=method)\n            seg_sens_gadf = gadf.fit_transform(seg_sens)\n            \n            df_img[r, :, 0+IMG_SIZE*r2:IMG_SIZE*(r2+1), 0] = seg_sens_gadf[0,:,:]\n        \n    return df_img\n    \n    \ndef create_cnn():\n    \"\"\"\n    Function that creates a CNN model using Keras\n    \"\"\"\n    model = Sequential()\n    model.add(Conv2D(32, (3,3), activation='relu', kernel_initializer='he_uniform',\n                     input_shape=(IMG_SIZE, IMG_SIZE*CHANNELS, 1)))\n    model.add(BatchNormalization())\n    model.add(MaxPooling2D())\n    model.add(Dropout(DROP))\n    model.add(Flatten())\n    model.add(Dense(16, activation='relu', kernel_initializer='he_uniform'))\n    model.add(Dense(1, activation='linear'))\n    model.compile(optimizer=Adam(lr=LR), loss='mae', metrics=['mae'])\n    print(model.summary())\n    return model\n    \n\ndef crossvalidate_model(kf, X_train_f, y_train_f, X_test_f, plt_hist=True):\n    \"\"\"\n    Function to crossvalidate model\n    \"\"\"\n    \n    iofs_preds = np.zeros((X_train_f.shape[0],1))\n    oofs_preds = np.zeros((X_train_f.shape[0],1))\n    test_preds = np.zeros((X_test_f.shape[0],1))\n    \n    for k, (trn_idx, val_idx) in enumerate(kf.split(X_train_f, y_train_f)):\n\n        X_trn_cv, y_trn_cv = X_train_f[trn_idx,:], y_train_f[trn_idx,:]\n        X_val_cv, y_val_cv = X_train_f[val_idx,:], y_train_f[val_idx,:]\n        \n        sc = MinMaxScaler(feature_range=(-1,1))\n        sc.fit(y_trn_cv)\n        y_trn_cv_sc = sc.transform(y_trn_cv)\n        y_val_cv_sc = sc.transform(y_val_cv)\n        \n        model = create_cnn()\n        history = model.fit(X_trn_cv, y_trn_cv_sc, \n                            epochs=EPOCHS, batch_size=BS, \n                            validation_data=(X_val_cv, y_val_cv_sc),\n                            verbose=1)\n        \n        if plt_hist:\n            fig, ax = plt.subplots(1,2, figsize=(10,5))\n            ax[0].plot(history.history['loss'])\n            ax[0].plot(history.history['val_loss'])\n            ax[1].plot(history.history['mae'])\n            ax[1].plot(history.history['val_mae'])\n            plt.show()\n        \n        y_pred_trn_cv = sc.inverse_transform(model.predict(X_trn_cv))\n        y_pred_val_cv = sc.inverse_transform(model.predict(X_val_cv))\n        y_pred_test = sc.inverse_transform(model.predict(X_test_f))\n        \n        print(\"Fold {} train MAE: {}\".format(k+1, mean_absolute_error(y_trn_cv, y_pred_trn_cv)))\n        print(\"Fold {} val MAE: {}\".format(k+1, mean_absolute_error(y_val_cv, y_pred_val_cv)))\n        \n        iofs_preds[trn_idx] = y_pred_trn_cv\n        oofs_preds[val_idx] = y_pred_val_cv\n        test_preds += y_pred_test / kf.get_n_splits()\n        \n    print(\"Overall train MAE: {}\".format(mean_absolute_error(y_train_f, iofs_preds)))\n    print(\"Overall val MAE: {}\".format(mean_absolute_error(y_train_f, oofs_preds)))\n    \n    return iofs_preds, oofs_preds, test_preds","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load data\ntrain_df = pd.read_csv(os.path.join(PATH,'train.csv'))\nsub = pd.read_csv(os.path.join(PATH,'sample_submission.csv'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Prepare data\ndf_list_train = train_df.copy()\ny_train = df_list_train['time_to_eruption'].values.reshape(-1,1)\nX_train = tseries_to_img(df_list_train, type_d='train', method='difference')\n\ndf_list_test = sub.copy()\nX_test = tseries_to_img(df_list_test, type_d='test', method='difference')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot some training and test images\nfig, ax = plt.subplots(4,1, figsize=(20,10))\nax[0].imshow(X_train[0,:,:,0])\nax[0].set_title('Segment '+str(df_list_train['segment_id'].iloc[0]))\nax[1].imshow(X_train[1,:,:,0])\nax[1].set_title('Segment '+str(df_list_train['segment_id'].iloc[1]))\nax[2].imshow(X_test[0,:,:,0])\nax[2].set_title('Segment '+str(df_list_test['segment_id'].iloc[0]))\nax[3].imshow(X_test[1,:,:,0])\nax[3].set_title('Segment '+str(df_list_test['segment_id'].iloc[1]))\nplt.show()\n\ndel df_list_train, df_list_test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Crossvalidate model\nkf = KFold(n_splits=N_SPLITS, shuffle=True, random_state=SEED)\niofs_preds, oofs_preds, test_preds = crossvalidate_model(kf, X_train, y_train, X_test, plt_hist=False)\n\n# Save predictions for test set\nsub['time_to_eruption'] = test_preds\nsub.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# nulls_df = pd.DataFrame(index=train_df['segment_id'], columns=['channels_nulls'])\n# nulls_seg = []\n# nulls_nulls = []\n# for s in tqdm(train_df['segment_id']):\n#     s_df = pd.read_csv(os.path.join(PATH, 'train', str(s)+'.csv'))\n#     if s_df.isnull().sum().any():\n#         tmp = s_df.isnull().sum()\n#         nulls_df.loc[s,'channels_nulls'] = str(tmp[tmp!=0].to_dict())\n#         nulls_seg.append(tmp[tmp!=0].index.tolist())\n#         nulls_nulls.append(tmp[tmp!=0].tolist())\n        \n# nulls_df = nulls_df.loc[nulls_df['channels_nulls'].notnull(),:]\n# nulls_seg = [c1 for c2 in n_seg for c1 in c2]\n# nulls_nulls = [c1 for c2 in n_nulls for c1 in c2]","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}