{"cells":[{"metadata":{},"cell_type":"markdown","source":"\nThis code is mainly taken from:\n\nhttps://github.com/danielegrattarola/spektral/blob/master/examples/other/graph_signal_classification_mnist.py\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install --upgrade pip","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install spektral","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Refer to an example, https://www.kaggle.com/delai50/ingv-imaging-ts-for-mt-doom-eruption-prediction\n!pip install pyts   \nfrom pyts.image import GramianAngularField","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n       os.path.join(dirname, filename)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport gc\ngc.enable()\nfrom tqdm.notebook import tqdm\n\n\nimport tensorflow as tf\nfrom tensorflow.keras import Model\nfrom tensorflow.keras.layers import Dense, Flatten\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.regularizers import l2\nfrom tensorflow.keras.losses import MeanAbsoluteError\n#from tensorflow.keras.metrics import MeanAbsoluteError\n\nfrom spektral.data import PackedBatchLoader\nfrom spektral.data import Dataset, Graph, DisjointLoader \nfrom spektral.layers import GCNConv\nfrom spektral.layers.ops import sp_matrix_to_sp_tensor\n\nimport scipy.sparse as sp\nfrom sklearn.neighbors import kneighbors_graph\nfrom sklearn.metrics import mean_absolute_error, r2_score\n\nimport scipy as scp  \nfrom scipy.signal import butter,filtfilt,freqz\nimport scipy.fftpack","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Detect hardware, return appropriate distribution strategy\ndef get_strategy():\n    gpu = \"\"\n    try:\n        tpu = tf.distribute.cluster_resolver.TPUClusterResolver()\n        print('Running on TPU ', tpu.master())     \n    except ValueError:\n        tpu = None\n        os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"0\"\n        gpu = tf.config.list_physical_devices(\"GPU\")\n        if len(gpu) == 1:\n            print('Running on GPU ', gpu)\n    if tpu:\n        tf.config.experimental_connect_to_cluster(tpu)\n        tf.tpu.experimental.initialize_tpu_system(tpu)\n        strategy = tf.distribute.experimental.TPUStrategy(tpu)\n    elif len(gpu) == 1:\n        strategy = tf.distribute.OneDeviceStrategy(device=\"/gpu:0\")\n        tf.config.optimizer.set_experimental_options({\"auto_mixed_precision\":True})\n    else:\n        strategy = tf.distribute.get_strategy()\n    print(\"REPLICAS: \", strategy.num_replicas_in_sync)\n    return strategy \n\nstrategy = get_strategy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# GLOBAL VARIABLES\nPATH = '/kaggle/input/predict-volcanic-eruptions-ingv-oe'\n\n#IMG_SIZE=10 is corr among sensors/10min, 20=corr among sensors/5min, else=GramianAngularField\nIMG_SIZE = 10          #10, 20, 28\nkk= int(0.3*IMG_SIZE)  #number of neighbours of each node.\n\n# Parameters\nCHANNELS = 1 \nbatch_size = 128  \nepochs = 1500        # Number of training epochs\npatience = epochs    # Patience for early stopping\nl2_reg = 5e-4        # Regularization rate for l2\n\n#BPF\nlow_cutoff = 1       #Hz, to remove dc\nhigh_cutoff = 35     #Hz\norder_bpf = 6        #filter order\nfs = 100             #sample rate, 100 Hz\nnyq = 0.5 * fs       #Nyquist Frequency","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DEBUG = True   \n\nif DEBUG:\n    batch_size = 32  \n    epochs = 1500","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def butter_bandpass_filter(data, cutoff_low, cutoff_high, fs, order):\n    normal_cutoff_low = cutoff_low / nyq\n    normal_cutoff_high = cutoff_high / nyq    \n    # Get the filter coefficients \n    b, a = butter(order, [normal_cutoff_low,normal_cutoff_high], btype='band', analog=False)\n    y = filtfilt(b, a, data)\n    return y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# HELPING FUNCTIONS\nLIMIT_PRECISION = 2**15-1     #16 bits precision on sensors\nMINUTES_10 = 60000            #60k=10min\nMINUTES_5  = 30000\nSPLIT_MINUTES = int(MINUTES_10/MINUTES_5)\n\ndef tseries_to_imgAvg(df_list, type_d, method='difference'):   \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    cols = ['s'+str(i)+str(j) for i in range(10) for j in range(SPLIT_MINUTES)]\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        seg_df = seg_df.iloc[:MINUTES_10,:]\n        df_rand = pd.DataFrame(np.random.randn(seg_df.shape[0],seg_df.shape[1]), \n                               columns=seg_df.columns, index=seg_df.index)/LIMIT_PRECISION\n        #seg_df = seg_df.fillna(0)                      \n        seg_df[pd.isna(seg_df)] = df_rand[pd.isna(seg_df)]               \n        \n        #BPF filtering\n        seg_columns = seg_df.columns\n        for i in seg_columns:\n            seg_df[i] = butter_bandpass_filter(seg_df[i].values, low_cutoff, high_cutoff, fs, order_bpf) \n         \n        if IMG_SIZE == 10:     #corr sensors/10min\n            A = seg_df.corr()\n            A = A.fillna(0)    #NaN if sensor =0\n            A[np.abs(A)<0.01]=0\n            np.fill_diagonal(A.values, 1)    #missing sensor\n            df_img[r, :, 0:IMG_SIZE, 0] = A  #(10,10) \n        elif IMG_SIZE == 20:    #corr sensors/5min\n            np_seg = seg_df.values\n            seg_sens = np_seg.reshape([MINUTES_5,-1],order='F')  #(30k, 20)\n            A = pd.DataFrame(seg_sens, columns=cols)\n            A = A.corr()\n            A = A.fillna(0)\n            A[np.abs(A)<0.01] = 0\n            np.fill_diagonal(A.values, 1)    #missing sensor\n            df_img[r, :, 0:IMG_SIZE, 0] = A  #(20,20)\n        else:       #GramianAngularField\n            #norm\n            train_input_mean = seg_df.mean()\n            train_input_sigma = seg_df.std()\n            seg_df_avg = (seg_df - train_input_mean) / train_input_sigma\n            seg_df_avg = seg_df_avg.fillna(0)\n            seg_df_avg['avg'] = seg_df_avg.mean(axis=1)           \n            seg_sens = seg_df_avg['avg'].values.reshape(1,-1)  #(1, 60000)              \n\n            gadf = GramianAngularField(image_size=IMG_SIZE, method=method)\n            seg_sens_gadf = gadf.fit_transform(seg_sens)\n            df_img[r, :, 0:IMG_SIZE, 0] = seg_sens_gadf[0,:,:]  \n            \n    return df_img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#BandPass Filter Frequency Response\n\nb2, a2 = butter(order_bpf, [low_cutoff/nyq, high_cutoff/nyq], btype='band', analog=False)\nw2,h2 = freqz(b2,a2, fs=fs)\nprint(len(w2))  #len(w2) = 50Hz\n\nplt.figure(figsize=(18,4))\nplt.subplot(131)\nplt.plot(w2, 20 * np.log10(abs(h2)+2**-31), 'b')  #add 2**-31 to avoid log(0)\nplt.ylabel('Amplitude [dB]', color='b')\nplt.xlabel('Frequency [Hz]')\nplt.grid()\nplt.title('BPF, passband @ 0.01Hz~35Hz')\n\nbins=20   # 1.9Hz = 50* 20/len(w2)\nplt.subplot(132)\nplt.plot(w2[:bins], 20 * np.log10(abs(h2[:bins])+2**-31), 'b')\nplt.ylabel('Amplitude [dB]', color='b')\nplt.xlabel('Frequency [Hz]')\nplt.grid()\nplt.title('Zoom-In BPF')\n\nbins=256  # 25Hz = 50* 256/len(w2)\nplt.subplot(133)\nplt.plot(w2[-bins:], 20 * np.log10(abs(h2[-bins:])+2**-31), 'b')\nplt.ylabel('Amplitude [dB]', color='b')\nplt.xlabel('Frequency [Hz]')\nplt.grid()\nplt.title('Zoom-In BPF')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train/601524801.csv')  #shortest erupttime=6250\ntrain_df.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#add random noise to NaN\ndt = train_df.copy()\ndf_rand = pd.DataFrame(np.random.randn(dt.shape[0],dt.shape[1]), columns=dt.columns, index=dt.index)/(2**15-1)\ndt[pd.isna(dt)] = df_rand[pd.isna(dt)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Original Sensors\nfig, ax = plt.subplots(nrows=10, ncols=2, figsize=(24, 30))\nfig.subplots_adjust(hspace = .5)\ncolors = plt.rcParams[\"axes.prop_cycle\"]()\n  \n# Sensor#1 ~ #10\nfor j in range(1,11):\n        \n    fft = np.fft.fft(dt[f'sensor_{j}'].values)\n    psd = np.abs(fft) ** 2\n    fftfreq = scp.fftpack.fftfreq(len(psd),1/fs)\n    i = fftfreq > 0\n\n    c = next(colors)[\"color\"]\n    ax[j-1,0].plot(dt[f'sensor_{j}'].values,color=c)\n    ax[j-1,0].set_title('Sensor_'+str(j))\n    ax[j-1,0].set_xlabel('time-domain samples')\n    ax[j-1,0].set_ylabel('Amp')\n\n    ax[j-1,1].plot(fftfreq[i], 10 * np.log10(psd[i]),color=c)\n    ax[j-1,1].set_title('Sensor_'+str(j))\n    ax[j-1,1].set_xlabel('Frequency (Hz)')\n    ax[j-1,1].set_ylabel('PSD (dB)')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sensor 4:Very strong low frequency signal close to DC. This created a non-stationary bias. Strong DC at sensors 1,4,6,7,9,10\n\nSensor 1,5,7,9,10: Impulsive signals\n\nSensor 2,3,8: Small random noises\n\nSensor 9: Peak PSD 125dB, noise floor PSD 50dB, dynamic range of sensor 125-50 = 75dB about 13bits precision"},{"metadata":{"trusted":true},"cell_type":"code","source":"#BPF on all sensors to remove DC\nseg_columns = dt.columns\nfor i in seg_columns:\n    dt[i] = butter_bandpass_filter(dt[i].values, low_cutoff, high_cutoff, fs, order_bpf) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#norm\ndt_mean = dt.mean()\ndt_sigma = dt.std()\ndt_avg = (dt - dt_mean) / dt_sigma\n#avg all sensors\ndt_avg['avg'] = dt_avg.mean(axis=1)  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=11, ncols=2, figsize=(24, 30))\nfig.subplots_adjust(hspace = .5)\ncolors = plt.rcParams[\"axes.prop_cycle\"]()\n  \n# Sensor#1 ~ #10\nfor j in range(1,11):\n        \n    fft = np.fft.fft(dt_avg[f'sensor_{j}'].values)\n    psd = np.abs(fft) ** 2\n    fftfreq = scp.fftpack.fftfreq(len(psd),1/fs)\n    i = fftfreq > 0\n\n    c = next(colors)[\"color\"]\n    ax[j-1,0].plot(dt_avg[f'sensor_{j}'].values,color=c)\n    ax[j-1,0].set_title('Sensor_'+str(j))\n    ax[j-1,0].set_xlabel('time-domain samples')\n    ax[j-1,0].set_ylabel('Amp')\n\n    ax[j-1,1].plot(fftfreq[i], 10 * np.log10(psd[i]),color=c)\n    ax[j-1,1].set_title('Sensor_'+str(j))\n    ax[j-1,1].set_xlabel('Frequency (Hz)')\n    ax[j-1,1].set_ylabel('PSD (dB)')\n\n#Average Sensor signals\nj +=1\nfft = np.fft.fft(dt_avg.avg)\npsd = np.abs(fft) ** 2\nfftfreq = scp.fftpack.fftfreq(len(psd),1/fs)\ni = fftfreq > 0\n\nc = next(colors)[\"color\"]\nax[j-1,0].plot(dt_avg.avg,color=c)\nax[j-1,0].set_title('Sensor_Avg')\nax[j-1,0].set_xlabel('time-domain samples')\nax[j-1,0].set_ylabel('Amp')\n\nax[j-1,1].plot(fftfreq[i], 10 * np.log10(psd[i]),color=c)\nax[j-1,1].set_title('Sensor_Avg')\nax[j-1,1].set_xlabel('Frequency (Hz)')\nax[j-1,1].set_ylabel('PSD (dB)')    \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Freq below 1Hz are removed.  Very obvious in sensors 4 time-domain signal.\n\nSensors 2,3,8 noises are amplified due to normalization.  Not sure this is a good idea.  For correlations usage, correlations with random noise should resulted in zero, this should be fine.  Another way is to average the sensors first, then normalize."},{"metadata":{"trusted":true},"cell_type":"code","source":"#Lets look at sensor 5 on every minute.\n\nfig, ax = plt.subplots(nrows=10, ncols=2, figsize=(24, 30))\nfig.subplots_adjust(hspace = .5)\ncolors = plt.rcParams[\"axes.prop_cycle\"]()\n  \n# Minute#1 ~ #10 on sensor 5\nfor j in range(10):\n        \n    fft = np.fft.fft(dt_avg.sensor_5[6000*(j):6000*(j+1)])\n    psd = np.abs(fft) ** 2\n    fftfreq = scp.fftpack.fftfreq(len(psd),1/fs)\n    i = fftfreq > 0\n\n    c = next(colors)[\"color\"]\n    ax[j,0].plot(dt_avg.sensor_5[6000*(j):6000*(j+1)],color=c)\n    ax[j,0].set_title('Minute'+str(j+1))\n    ax[j,0].set_xlabel('time-domain samples')\n    ax[j,0].set_ylabel('Amp')\n\n    ax[j,1].plot(fftfreq[i], 10 * np.log10(psd[i]),color=c)\n    ax[j,1].set_title('Minute'+str(j+1))\n    ax[j,1].set_xlabel('Frequency (Hz)')\n    ax[j,1].set_ylabel('PSD (dB)')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_df,dt,dt_avg\ngc.collect()","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":"train_df.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_list_train = train_df.copy()\ndf_list_test = sub.copy()\n\nif DEBUG:\n    df_list_train = df_list_train[df_list_train[\"time_to_eruption\"] < 600_000]  #debug\n    df_list_test = df_list_test[:32]  #debug\n    \ndf_list_train.shape, df_list_test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Prepare data\nX_train = tseries_to_imgAvg(df_list_train, type_d='train', method='difference')\ny_train = df_list_train['time_to_eruption'].values.reshape(-1,1)\nX_test = tseries_to_imgAvg(df_list_test, type_d='test', method='difference')\ny_test = df_list_test['time_to_eruption'].values.reshape(-1,1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.shape, y_train.shape, X_test.shape, y_test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#norm train dataset\ntrain_mean = X_train.mean()\ntrain_sigma = X_train.std()\ntrain_norm = (X_train - train_mean) / train_sigma\n\n#scale\nx_train_max = train_norm.max()\nx_train_min = train_norm.min()\nX_train_std = (train_norm - x_train_min) / (x_train_max - x_train_min)\n\ny_train_max = y_train.max()\ny_train_min = y_train.min()\ny_train_std = (y_train - y_train_min) / (y_train_max - y_train_min)\n#y_train_scaled = y_train_std * (y_train_max - y_train_min) + y_train_min\n\n#Apply train mean, std to the test dataset\n#norm test dataset\ntest_norm = (X_test - train_mean) / train_sigma\n#scale\nx_test_max = test_norm.max()\nx_test_min = test_norm.min()\nX_test_std = (test_norm - x_test_min) / (x_test_max - x_test_min)\ny_test_std = y_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mean, train_sigma, x_train_max, x_train_min, y_train_max,y_train_min, x_test_max,x_test_min","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_std.max(), X_train_std.min(), y_train_std.max(), y_train_std.min(), X_test_std.max(), X_test_std.min()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_list_train.shape, y_train.shape, X_train_std.shape, y_train.shape, X_test_std.shape, y_test.shape ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot some training and test images\nfig, ax = plt.subplots(1,4, figsize=(20,14))\n\n#train\nax[0].imshow(X_train_std[0,:,:,0])\nax[0].set_title('Segment '+str(df_list_train['segment_id'].iloc[0]))\nax[1].imshow(X_train_std[1,:,:,0])\nax[1].set_title('Segment '+str(df_list_train['segment_id'].iloc[1]))\n#test\nax[2].imshow(X_test_std[0,:,:,0])\nax[2].set_title('Segment '+str(df_list_test['segment_id'].iloc[0]))\nax[3].imshow(X_test_std[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":"features = np.squeeze(X_train_std)\nlabels = np.squeeze(y_train_std)\nX_test = np.squeeze(X_test_std)\ny_test = np.squeeze(y_test)\n        \nprint(features.shape, labels.shape, X_test.shape, y_test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class INGV(Dataset):   \n    def __init__(self, features=features, labels=labels, p_flip=0., k=8, **kwargs):\n        self.a = None\n        self.features = features\n        self.labels = labels\n        self.k = k\n        self.p_flip = p_flip\n        super().__init__(**kwargs)\n    \n    def read(self):\n        self.a = _mnist_grid_graph(self.k)\n        self.a = _flip_random_edges(self.a, self.p_flip)\n\n        x = self.features.reshape(-1, IMG_SIZE**2, 1)\n        y = self.labels\n        \n        return [Graph(x=x_, y=y_) for x_, y_ in zip(x, y)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def _grid_coordinates(side):\n    M = side ** 2\n    x = np.linspace(0, 1, side, dtype=np.float32)\n    y = np.linspace(0, 1, side, dtype=np.float32)\n    xx, yy = np.meshgrid(x, y)\n    z = np.empty((M, 2), np.float32)\n    z[:, 0] = xx.reshape(M)\n    z[:, 1] = yy.reshape(M)\n    return z\n\n\ndef _get_adj_from_data(X, k, **kwargs):\n    A = kneighbors_graph(X, k, **kwargs).toarray()\n    A = sp.csr_matrix(np.maximum(A, A.T))\n    return A\n\n\ndef _mnist_grid_graph(k):\n    X = _grid_coordinates(IMG_SIZE)\n    A = _get_adj_from_data(X, k, \n                           mode='connectivity', \n                           metric='euclidean', \n                           include_self=False)\n    return A\n\n\ndef _flip_random_edges(A, percent):\n    if not A.shape[0] == A.shape[1]:\n        raise ValueError('A must be a square matrix.')\n    \n    dtype = A.dtype\n    A = sp.lil_matrix(A).astype(np.bool)\n    \n    n_elem = A.shape[0] ** 2\n    n_elem_to_flip = round(percent * n_elem)\n    unique_idx = np.random.choice(n_elem, replace=False, size=n_elem_to_flip)\n    row_idx = unique_idx // A.shape[0]\n    col_idx = unique_idx % A.shape[0]\n    idxs = np.stack((row_idx, col_idx)).T\n    \n    for i in idxs:\n        i = tuple(i)\n        A[i] = np.logical_not(A[i])\n    A = A.tocsr().astype(dtype)\n    A.eliminate_zeros()\n    return A","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Build model\nclass Net(Model):\n    def __init__(self, **kwargs):\n        super().__init__(**kwargs)\n        self.conv1 = GCNConv(32, activation='elu', kernel_regularizer=l2(l2_reg))\n        self.conv2 = GCNConv(32, activation='elu', kernel_regularizer=l2(l2_reg))\n        self.flatten = Flatten()\n        self.fc1 = Dense(512, activation='relu')  \n        #self.fc2 = Dense(1, activation='linear')\n        self.fc2 = Dense(1, activation='sigmoid') \n        \n    def call(self, inputs):\n        x, a = inputs\n        x = self.conv1([x, a])\n        x = self.conv2([x, a])\n        output = self.flatten(x)\n        output = self.fc1(output)\n        output = self.fc2(output)\n\n        return output","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create model\n\nwith strategy.scope():\n    model = Net()\n    optimizer = Adam(lr=0.0001) \n    loss_fn = MeanAbsoluteError()\n    acc_fn = tf.keras.metrics.MeanAbsoluteError()    \n    \n    #model.compile(optimizer=optimizer,\n    #              loss='mae', \n    #              metrics=['mae'])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Training function\n@tf.function\ndef train_on_batch(inputs, target):\n    with tf.GradientTape() as tape:\n        predictions = model(inputs, training=True)\n        loss = loss_fn(target, predictions) + sum(model.losses)\n        acc = acc_fn(target, predictions)\n\n    gradients = tape.gradient(loss, model.trainable_variables)\n    optimizer.apply_gradients(zip(gradients, model.trainable_variables))\n    return loss, acc\n\n\n# Evaluation function\ndef evaluate(loader):\n    step = 0\n    results = []\n    for batch in loader:\n        step += 1\n        x, target = batch\n        predictions = model([x, adj], training=False)\n        loss = loss_fn(target, predictions)\n        acc = acc_fn(target, predictions)\n        results.append((loss, acc, len(target)))  # Keep track of batch size\n        if step == loader.steps_per_epoch:\n            results = np.array(results)\n            return np.average(results[:, :-1], 0, weights=results[:, -1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# The adjacency matrix is stored as an attribute of the dataset.\n# Create filter for GCN and convert to sparse tensor.\n\ndata = INGV(features, labels, k=kk, p_flip=0 )\n\nadj = data.a\nadj = GCNConv.preprocess(adj)\nadj = sp_matrix_to_sp_tensor(adj)\n\n# Train/valid split\np_split = int(0.20 * len(data))  #20 percent split\nnp.random.shuffle(data)\ndata_tr, data_va = data[:-p_split], data[-p_split:]\n\nlen(data_tr), len(data_va)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setup training# Setup training\nbest_val_loss = 99999\ncurrent_patience = patience\nstep = 0\n\n# We can use PackedBatchLoader because we only need to create batches of node\n# features with the same dimensions.\nloader_tr = PackedBatchLoader(data_tr, batch_size=batch_size, epochs=epochs)\nloader_va = PackedBatchLoader(data_va, batch_size=batch_size)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Training loop\nresults_tr = []\n\nloss_history = []\nepoch_cnt = 0\nepoch_mod =1\n\nif (100 < epochs < 1000):\n    epoch_mod=10\nelif (epochs > 1000):\n    epoch_mod=20\n\nfor batch in loader_tr:\n    step += 1\n\n    # Training step\n    x, y = batch                              #x(batch_size, IMG**2), y(batch_size,)\n    #l, a = model.train_on_batch([x, adj], y) #TPUs do not support DT_VARIANT\n    l, a = train_on_batch([x, adj], y)        #l=loss, a=acc\n    results_tr.append((l, a, len(y)))\n\n    if step == loader_tr.steps_per_epoch:\n        results_va = evaluate(loader_va)      #Validate\n        if results_va[0] < best_val_loss:\n            best_val_loss = results_va[0]\n            current_patience = patience\n            #results_te = evaluate(loader_te)\n        else:\n            current_patience -= 1\n            if current_patience == 0:\n                print('Early stopping')\n                break\n\n        # Print results\n        results_tr = np.array(results_tr)\n        results_tr = np.average(results_tr[:, :-1], 0, weights=results_tr[:, -1])  \n        MAE = results_va[1] * (y_train_max - y_train_min) + y_train_min\n        epoch_cnt +=1\n        \n        if (epoch_cnt%epoch_mod ==0):               \n            print('Epoch {:3}/{:<3} | '\n                  'Train loss: {:.5f}, mae: {:.5f} | '\n                  'Val loss: {:.5f}, mae: {:.5f} | '\n                  'MAE: {:10.1f}'\n                  .format(epochs, epoch_cnt, *results_tr, *results_va, MAE))\n        \n        loss_history.append((results_tr[0], results_tr[1], \n                             results_va[0], results_va[1]))\n        \n        # Reset epoch\n        results_tr = []\n        step = 0\n        \nloss_history = np.array(loss_history)         ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plots\nplt.figure(figsize=(10, 5))\n\nplt.subplot(121)\nplt.plot(loss_history[:, 0], label='Train loss')\nplt.plot(loss_history[:, 2], label='Val loss')\nplt.legend()\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\n\nplt.subplot(122)\nplt.plot(loss_history[:, 1], label='Train acc')\nplt.plot(loss_history[:, 3], label='Val acc')\nplt.legend()\nplt.ylabel('Acc')\nplt.xlabel('Epoch')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Validate model')\ny_pred=[]\ny_true=[]\n\nloader_va = PackedBatchLoader(data_va, batch_size=batch_size, epochs=1)\n#batches = [b for b in loader_va]\n#x,y = batches[-1]\n\nfor batch in loader_va:\n    x, y = batch   \n    p_va = model([x, adj], training=False)  #predict label per batch  \n    y_pred.append(p_va)   \n    y_va = np.vstack(y)   #True label per batch\n    y_true.append(y_va)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x.shape, y.shape, np.shape(y_true),np.shape(y_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Last batch of Validate images\n\nval_images = np.squeeze(x)     #(batch,IMG_SIZE**2) <- (batch,IMG_SIZE**2,1)\nval_images = np.reshape(val_images,[val_images.shape[0], IMG_SIZE,-1])\n\ny_va = np.squeeze(y_va)\np_va = np.squeeze(p_va)\n\nnum_plt = int(np.sqrt(len(val_images)))\nnum_plt = np.amin([5,num_plt])\nnum_fig = num_plt**2\n\nplt.figure(figsize=(13,14))\nfor i in range(num_fig):\n    plt.subplot(num_plt,num_plt,i+1)\n    plt.xticks([])\n    plt.yticks([])\n    plt.grid(False)\n    plt.imshow(val_images[i], cmap=plt.cm.binary)\n    plt.xlabel(p_va[i]) #Predict label\n    plt.title(y_va[i])  #True label\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Testing model')\n\ndata_testset = INGV(X_test_std, y_test, k=kk, p_flip=0 )\n\nadj_tst = data_testset.a\nadj_tst = GCNConv.preprocess(adj_tst)\nadj_tst = sp_matrix_to_sp_tensor(adj_tst)\n\ndata_te = data_testset\nlen(data_te)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#loader_te = PackedBatchLoader(data_te, batch_size=data_te.n_graphs, epochs=1)\nloader_te = PackedBatchLoader(data_te, batch_size=batch_size, epochs=1)\n\ny_pred=[]\nfor batch in loader_te:\n    x_te, y_te = batch    \n    p_te = model([x_te, adj_tst], training=False)    \n    y_pred.append(p_te)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot last batch of Test images\n\ntest_images = np.squeeze(x_te)\ntest_images = np.reshape(test_images,[test_images.shape[0], IMG_SIZE,-1])\np_te = np.squeeze(p_te)\n\n#Plot Test images\nnum_plt = int(np.sqrt(len(test_images)))\nnum_plt = np.amin([5,num_plt])\nnum_fig = num_plt**2\n\nplt.figure(figsize=(10,10))\nfor i in range(num_fig):\n    plt.subplot(num_plt,num_plt,i+1)\n    plt.xticks([])\n    plt.yticks([])\n    plt.grid(False)\n    plt.imshow(test_images[i], cmap=plt.cm.binary)\n    plt.xlabel(p_te[i])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Submission:**"},{"metadata":{"trusted":true},"cell_type":"code","source":"flat_list=[item for sublist in y_pred for item in sublist]       \nflat_list= np.squeeze(flat_list)\nnp.shape(flat_list)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred_scaled = flat_list * (y_train_max - y_train_min) + y_train_min\nlen(y_pred_scaled)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.shape, y_pred_scaled.shape ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if not DEBUG:\n    sub['time_to_eruption'] = y_pred_scaled\n    sub.to_csv('submission.csv', header=True, index=False)\n    !cat submission.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}