{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom scipy.signal import spectrogram\nimport pywt\nimport pywt.data\n\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\n\nimport io\nimport cv2\nfrom PIL import Image, ImageChops\n\n_input_dir = '/kaggle/input/predict-volcanic-eruptions-ingv-oe/'\n_train_dir = f'{_input_dir}train/'\n_test_dir = f'{_input_dir}test/'\n\ntrain = pd.read_csv(f'{_input_dir}train.csv')\ntest = pd.read_csv(f'{_input_dir}sample_submission.csv')\n\n# Resample starting at 25th earthquake getting every 15th earthquake in order of increasing time to eruption\ntrain_resample = train.sort_values(['time_to_eruption'])\\\n                      .reset_index(drop=True) \\\n                      .iloc[10::15, :]\ntrain_resample.head(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Regression using spectrograms and Convolutional NNs\n\n#### Methodology\nIn my other notebook: https://www.kaggle.com/ajcostarino/ingv-volcanic-eruption-prediction-lgbm-baseline. I presented a way of building a baseline prediction by aggregating sensor values from each segment and then using a classic Gradient Boosted Tree to predict time to failure. Here I tried a different approach: First we denoise the signals using a wavelet transformation. We then convert the transformed signal to a spectrogram. Each sensor has an individual spectrogram, the 10 sensors all together form a set of spectrograms for each earthquake. This can be described by a 3D tensor of size (`N` x `C` x `D` x `W` x `H`). Where `W, H` are the width and height of the spectrograms, `D` is the depth or the number of spectrograms in our case `10`. `C` is the number of channels `3` for RGB. `N` is the size of the batch. We can then use 3D convolutions to extract features from the set of spectrograms and build our regression model using those features.\n\n#### Wavelet Denoise\nHere we use some code to remove high frequency noise from our sensors."},{"metadata":{"trusted":true},"cell_type":"code","source":"def maddest(d, axis=None):\n    \"\"\"\n    Mean Absolute Deviation\n    \"\"\"\n    \n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef denoise_signal(x, wavelet='db4', level=1):\n    \"\"\"\n    1. Adapted from waveletSmooth function found here:\n    http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/\n    2. Threshold equation and using hard mode in threshold as mentioned\n    in section '3.2 denoising based on optimized singular values' from paper by Tomas Vantuch:\n    http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    \"\"\"\n    \n    # Decompose to get the wavelet coefficients\n    coeff = pywt.wavedec(x, wavelet, mode=\"per\")\n    \n    # Calculate sigma for threshold as defined in http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    # As noted by @harshit92 MAD referred to in the paper is Mean Absolute Deviation not Median Absolute Deviation\n    sigma = (1/0.6745) * maddest(coeff[-level])\n\n    # Calculate the univeral threshold\n    uthresh = sigma * np.sqrt(2*np.log(len(x)))\n    coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])\n    \n    # Reconstruct the signal using the thresholded coefficients\n    return pywt.waverec(coeff, wavelet, mode='per')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.experimental import enable_iterative_imputer\nfrom sklearn.impute import IterativeImputer\n\ndef denoise_transform(signals):\n    for i in range(1, 11):\n        signal_col = f'sensor_{i}'\n        if (signals[signal_col].isnull().sum() > 0) & (signals[signal_col].isnull().sum() < 20000):\n            # Impute the missing values if there are gaps in our data\n            imp = IterativeImputer(max_iter=10, random_state=42)\n            signals[signal_col] = imp.fit_transform(signals[signal_col].to_numpy().reshape(-1, 1))[:,0]\n        \n        signals[signal_col] = denoise_signal(signals[signal_col], level=2)[:-1]\n    \n    return signals\n\n\ndef transform_signals(signals):\n    for i in range(1, 11):\n        signal_col = f'sensor_{i}'\n        if (signals[signal_col].isnull().sum() > 0) & (signals[signal_col].isnull().sum() < 20000):\n            imp = IterativeImputer(max_iter=10, random_state=42)\n            signals[signal_col] = imp.fit_transform(signals[signal_col].to_numpy().reshape(-1, 1))[:,0]\n        \n        \n        for l in range(1,3):\n            signals[f'sensor_{i}_l{l}'] = denoise_signal(signals[signal_col], level=l)[:-1]\n            signals[f'sensor_{i}_l{l}_sum'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).sum()\n            signals[f'sensor_{i}_l{l}_mean'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).mean()\n            signals[f'sensor_{i}_l{l}_std'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).std()\n            signals[f'sensor_{i}_l{l}_var'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).var()\n            signals[f'sensor_{i}_l{l}_max'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).max()\n            signals[f'sensor_{i}_l{l}_min'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).min()\n            signals[f'sensor_{i}_l{l}_median'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).median()\n            signals[f'sensor_{i}_l{l}_skew'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).skew()\n            signals[f'sensor_{i}_l{l}_kurt'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).kurt()\n            signals[f'sensor_{i}_l{l}_q95'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).quantile(.95, interpolation='midpoint')\n            signals[f'sensor_{i}_l{l}_q87'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).quantile(.87, interpolation='midpoint')\n            signals[f'sensor_{i}_l{l}_q13'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).quantile(.13, interpolation='midpoint')\n            signals[f'sensor_{i}_l{l}_q05'] = signals[f'sensor_{i}_l{l}'].rolling(50, min_periods=1).quantile(.05, interpolation='midpoint')\n    \n    return signals","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualizing sensor signals and their respective spetrograms\nHere we visualize what a sensor signal looks like. The sensor denoised. The sensor spectrogram, and the denoised signal spectrogram. While we can't see the difference too much between the regular singal graphs, as we can see here the denoised spectrogram has less noise in the image. "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(4, 4, figsize=(40, 20))\n\ndef plot_row(temp, sensor, ax1, ax2, ax3, ax4):\n    ax1.plot(temp.index, temp[sensor], color='pink')\n    ax1.set_title(f'{sensor}')\n                  \n    ax2.plot(temp.index, denoise_transform(temp.copy())[sensor], color='pink')\n    ax2.set_title(f'{sensor} Denoised')\n    \n    # Spectrogram\n    f, t, Sxx = spectrogram(temp[sensor])\n    ax3.pcolormesh(t, f, np.log10(Sxx), shading='nearest')\n    ax3.set_title(f'{sensor} Spectrogram')\n    ax3.set_ylabel('f [Hz]')\n    ax3.set_xlabel('t [sec]')\n    ax3.set_yscale('symlog')\n    \n    # Denoise Spectrogram\n    f, t, Sxx = spectrogram(denoise_transform(temp.copy())[sensor])\n    ax4.pcolormesh(t, f, np.log10(Sxx), shading='nearest')\n    ax4.set_title(f'{sensor} Denoised Spectrogram')\n    ax4.set_ylabel('f [Hz]')\n    ax4.set_xlabel('t [sec]')\n    ax4.set_yscale('symlog')\n    \n\ntemp = pd.read_csv(f'{_train_dir}1593620672.csv')\nplot_row(temp, 'sensor_1', axes[0, 0], axes[0, 1], axes[0, 2], axes[0, 3])\nplot_row(temp, 'sensor_2', axes[1, 0], axes[1, 1], axes[1, 2], axes[1, 3])\nplot_row(temp, 'sensor_3', axes[2, 0], axes[2, 1], axes[2, 2], axes[2, 3])\nplot_row(temp, 'sensor_4', axes[3, 0], axes[3, 1], axes[3, 2], axes[3, 3])\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Spectrogram far from eruption\nSo an individual segment has 10 images associated with it, or 10 layers of information. Depicted below"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(3, 3, figsize=(40, 15))\n\ndef plot_denoise_spectrogram(temp, sensor, ax):\n    # Denoise Spectrogram\n    f, t, Sxx = spectrogram(denoise_transform(temp.copy())[sensor])\n    ax.pcolormesh(t, f, np.log10(Sxx), shading='nearest')\n    ax.set_title(f'{sensor} Denoised Spectrogram')\n    ax.set_ylabel('f [Hz]')\n    ax.set_xlabel('t [centi-seconds]')\n    ax.set_yscale('symlog')\n    \n    return ax\n\n    \ntemp = pd.read_csv(f'{_train_dir}1992733806.csv')\nplot_denoise_spectrogram(temp, 'sensor_1', axes[0, 0])\nplot_denoise_spectrogram(temp, 'sensor_2', axes[0, 1])\nplot_denoise_spectrogram(temp, 'sensor_3', axes[0, 2])\nplot_denoise_spectrogram(temp, 'sensor_4', axes[1, 0])\nplot_denoise_spectrogram(temp, 'sensor_5', axes[1, 1])\nplot_denoise_spectrogram(temp, 'sensor_6', axes[1, 2])\nplot_denoise_spectrogram(temp, 'sensor_7', axes[2, 0])\nplot_denoise_spectrogram(temp, 'sensor_8', axes[2, 1])\nplot_denoise_spectrogram(temp, 'sensor_9', axes[2, 2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot Spectrogram close to eruption\nSo an individual segment has 10 images associated with it, or 10 layers of information. Depicted below"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(3, 3, figsize=(40, 15))\n    \ntemp = pd.read_csv(f'{_train_dir}1826701813.csv')\nplot_denoise_spectrogram(temp, 'sensor_1', axes[0, 0])\nplot_denoise_spectrogram(temp, 'sensor_2', axes[0, 1])\nplot_denoise_spectrogram(temp, 'sensor_3', axes[0, 2])\nplot_denoise_spectrogram(temp, 'sensor_4', axes[1, 0])\nplot_denoise_spectrogram(temp, 'sensor_5', axes[1, 1])\nplot_denoise_spectrogram(temp, 'sensor_6', axes[1, 2])\nplot_denoise_spectrogram(temp, 'sensor_7', axes[2, 0])\nplot_denoise_spectrogram(temp, 'sensor_8', axes[2, 1])\nplot_denoise_spectrogram(temp, 'sensor_9', axes[2, 2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def spectrogram_data(temp, sensor):\n    # Denoise Spectrogram\n    fig = plt.figure(figsize=(40, 35))\n    f, t, Sxx = spectrogram(denoise_transform(temp.copy())[sensor])\n    plt.pcolormesh(t, f, np.log10(Sxx), shading='nearest')\n    plt.title(f'{sensor} Denoised Spectrogram')\n    plt.ylabel('f [Hz]')\n    plt.xlabel('t [centi-seconds]')\n    plt.yscale('symlog')\n    plt.close()\n    return fig\n\n\ndef trim(im):\n    bg = Image.new(im.mode, im.size, im.getpixel((0,0)))\n    diff = ImageChops.difference(im, bg)\n    bbox = diff.getbbox()\n    if bbox:\n        return im.crop(bbox)\n\n    \ndef get_img_from_fig(fig, dpi=0):\n    buf = io.BytesIO()\n    fig.savefig(buf, format=\"png\", dpi=None)\n    buf.seek(0)\n    img_arr = np.frombuffer(buf.getvalue(), dtype=np.uint8)\n    buf.close()\n    img = cv2.imdecode(img_arr, 1)\n    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n    \n\n    return img[305:2200][:,360:2500]\n\n\ndef get_sensor_spectrograms(signals):\n    spectrograms = []\n    for i in range(1, 11):\n        sensor_col = f'sensor_{i}'\n        fig = spectrogram_data(signals, sensor_col)\n        spectrogram = get_img_from_fig(fig)\n        spectrograms.append(spectrogram)\n        \n    return np.array(spectrograms)\n    \ngrams = get_sensor_spectrograms(temp)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## We can convert out spectrograms to Image arrays\nUsing the code here we can convert each spectrogram to an image array."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(1, 3, figsize=(40, 5))\n\naxes[0].imshow(grams[0], aspect='auto')\naxes[1].imshow(grams[1], aspect='auto')\naxes[2].imshow(grams[2], aspect='auto')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"60000/10/10/10","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class ConvNN(nn.Module):\n    \n    def __init__(self):\n        super().__init__()\n        self.conv1 = nn.Conv1d(10, 15, kernel_size=4, padding=1)\n        self.act1 =nn.Tanh()\n        self.pool1 = nn.MaxPool1d(10, ceil_mode=True)\n        self.conv2 = nn.Conv1d(15, 5, kernel_size=3, padding=1)\n        self.act2 = nn.Tanh()\n        self.pool2 = nn.MaxPool1d(10, ceil_mode=True)\n        self.conv3 = nn.Conv1d(5, 1, kernel_size=3, padding=1)\n        self.act3 = nn.Tanh()\n        self.pool3 = nn.MaxPool1d(10, ceil_mode=True)\n        \n        self.flatten1 = nn.Flatten()\n        self.flatten2 = nn.Flatten()\n        self.fc1 = nn.Linear(60, 32)\n        self.act4 = nn.Tanh()\n        self.fc2 = nn.Linear(32, 1)\n    \n    def forward(self, x):\n        out = self.pool1(self.act1(self.conv1(x)))\n        out = self.pool2(self.act2(self.conv2(out)))\n        out = self.pool3(self.act3(self.conv3(out)))\n        out = self.flatten2(self.flatten1(out))\n        #print(out.shape)\n        out = self.act4(self.fc1(out))\n        out = self.fc2(out)\n        return out\n\nmodel = ConvNN()\nloss_fn = nn.L1Loss(reduction='mean')\n\noptimizer = torch.optim.Adam(model.parameters(), lr=0.1)\nprint(model)\nprint(len(list(model.parameters())))\nfor i in range(len(list(model.parameters()))):\n    print(list(model.parameters())[i].size())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch.nn.functional as F\n\nclass TrainDataset(Dataset):\n    \n    def __init__(self, segment_ids):\n        self.segment_ids = segment_ids\n        \n    def __getitem__(self, index):\n        segment_id = self.segment_ids[index]\n        signals = pd.read_csv(f'{_train_dir}{segment_id}.csv')\n        grams = denoise_transform(signals).to_numpy()\n        del signals\n        for i in range(0, 10):\n            grams[i] = grams[i] / 32767.0\n        grams = torch.tensor(grams)\n        grams = grams.permute(1, 0)\n        ttf  = torch.tensor(train.set_index('segment_id').loc[segment_id]['time_to_eruption'])\n        return (grams, ttf)\n        \n    def __len__(self):\n        return len(self.segment_ids)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def training_loop(n_epochs, optimizer, model, loss_fn, train_loader):\n    device = torch.device('cuda')\n    model.to(device)\n    \n    for epoch in range(0, n_epochs):\n        loss_train = 0.0\n        \n        i = 1\n        for grams, ttf in train_loader:\n            grams = grams.to(device).float()\n            ttf = ttf.to(device).float()\n            outputs = model(grams)\n            loss = loss_fn(outputs, ttf)\n            optimizer.zero_grad()\n            loss.backward()\n            optimizer.step()\n            loss_train += loss.item()\n            del grams\n            del ttf\n            del outputs\n            \n            if i % 100 == 0:\n                print(f'Epoch {epoch}, Training Loss {loss_train / i}')\n            i += 1\n            \n        print(f'Epoch {epoch}, Training Loss {loss_train/i}')\n\n\ntrd = TrainDataset(list(train.segment_id[::-1]))\ntrain_loader = DataLoader(\n    trd,\n    batch_size = 1,\n    shuffle=True\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import warnings\nwith warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    training_loop(15, optimizer, model, loss_fn, train_loader)","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}