{"metadata":{"accelerator":"GPU","colab":{"gpuType":"A100","machine_shape":"hm","provenance":[]},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":87230,"databundleVersionId":9911932,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# Mount Google Drive\n# This is just for Google Colab\nfrom google.colab import drive\ndrive.mount('/content/drive')\n\n# Set the directory at MyDrive\n# This is just when running on Google colab\nimport os\nos.chdir('/content/drive/MyDrive/')","metadata":{"id":"6hC4261swlwJ","outputId":"0f9ea675-a0f7-484b-cb24-439ddaff26a6"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nfrom matplotlib import pyplot as plt\nimport cv2\nfrom sklearn.model_selection import train_test_split\nimport copy\nimport tensorflow as tf\nfrom tensorflow.keras.layers import Input, Conv2D, MaxPooling2D, UpSampling2D, concatenate, Cropping2D, Dropout\nfrom tensorflow.keras.models import Model, load_model\nfrom tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom scipy import interpolate\nfrom scipy.interpolate import NearestNDInterpolator\n\n# Load the training data (BScans)\n# np.shape(train) = (number of samples, 230, 230)\n# The Bscans are already normalised\ntrain = np.load('Bscan_Kaggle.npy')\n\n# Load the labels (permittivity ground truth)\n# np.shape(labels) = (number of samples, 224, 224)\nlabels = np.load('Labels_Kaggle.npy')\n\n\nfor i, c in enumerate(train):\n    # Go through my training data i.e. BScans\n    # Normalise them from -1 to 1\n    train[i] = c/np.max(np.abs(c))\n    # Normalise the labels from -0.5 to 0.5\n    labels[i] = labels[i]/10 - 0.5\n\n    # plt.imshow(labels[i])\n    # plt.colorbar()\n    # plt.show()\n","metadata":{"id":"tmwajFejwsZj"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras.layers import ZeroPadding2D, BatchNormalization\nfrom tensorflow.keras.initializers import he_normal\n\n\n\ndef unet_model(input_size=(230, 230, 1)):\n    # This is a vanilla typical U-net that takes an image as input with 230x230 dimensions and gives another image as an output with 224x224 dimensions\n    # The hyperparameters are not tuned.\n\n\n    inputs = Input(input_size)\n\n    # Size of the kernels\n    ra=4\n    # Encoding path\n    c1 = Conv2D(16, (ra, ra), activation='relu', padding='same', kernel_initializer=he_normal())(inputs)\n    # c1 = BatchNormalization()(c1)\n    c1 = Dropout(0.1)(c1)\n    c1 = Conv2D(16, (ra, ra), activation='relu', padding='same', kernel_initializer=he_normal())(c1)\n    # c1 = BatchNormalization()(c1)\n    p1 = MaxPooling2D((2, 2))(c1)\n\n    c2 = Conv2D(32, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(p1)\n    # c2 = BatchNormalization()(c2)\n    c2 = Dropout(0.1)(c2)\n    c2 = Conv2D(32, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c2)\n    # c2 = BatchNormalization()(c2)\n    p2 = MaxPooling2D((2, 2))(c2)\n\n    c3 = Conv2D(64, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(p2)\n    # c3 = BatchNormalization()(c3)\n    c3 = Dropout(0.2)(c3)\n    c3 = Conv2D(64, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c3)\n    # c3 = BatchNormalization()(c3)\n    p3 = MaxPooling2D((2, 2))(c3)\n\n    c4 = Conv2D(128, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(p3)\n    # c4 = BatchNormalization()(c4)\n    c4 = Dropout(0.2)(c4)\n    c4 = Conv2D(128, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c4)\n    # c4 = BatchNormalization()(c4)\n    p4 = MaxPooling2D(pool_size=(2, 2))(c4)\n\n    c5 = Conv2D(256, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(p4)\n    # c5 = BatchNormalization()(c5)\n\n    c5 = Dropout(0.3)(c5)\n    c5 = Conv2D(256, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c5)\n    # c5 = BatchNormalization()(c5)\n\n    # Decoding path\n    u6 = UpSampling2D((2, 2))(c5)\n    # c4 = Cropping2D(cropping=((1, 0), (1, 0)))(c4)\n    u6 = concatenate([u6, c4])\n    c6 = Conv2D(128, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(u6)\n    # c6 = BatchNormalization()(c6)\n\n    c6 = Dropout(0.2)(c6)\n    c6 = Conv2D(128, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c6)\n    # c6 = BatchNormalization()(c6)\n\n    u7 = UpSampling2D((2, 2))(c6)\n\n    c3 = Cropping2D(cropping=((1, 0), (1, 0)))(c3)\n    u7 = concatenate([u7, c3])\n    c7 = Conv2D(64, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(u7)\n    # c7 = BatchNormalization()(c7)\n\n    c7 = Dropout(0.2)(c7)\n    c7 = Conv2D(64, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c7)\n    # c7 = BatchNormalization()(c7)\n\n    u8 = UpSampling2D((2, 2))(c7)\n    c2 = Cropping2D(cropping=((2, 1), (2, 1)))(c2)\n    u8 = concatenate([u8, c2])\n    c8 = Conv2D(32, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(u8)\n    # c8 = BatchNormalization()(c8)\n\n    c8 = Dropout(0.1)(c8)\n    c8 = Conv2D(32, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c8)\n    # c8 = BatchNormalization()(c8)\n\n    u9 = UpSampling2D((2, 2))(c8)\n    c1 = Cropping2D(cropping=((3, 3), (3, 3)))(c1)\n    u9 = concatenate([u9, c1], axis=3)\n    c9 = Conv2D(32, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(u9)\n    # c9 = BatchNormalization()(c9)\n\n    c9 = Dropout(0.1)(c9)\n    c9 = Conv2D(16, (ra, ra), activation='relu', padding='same',kernel_initializer=he_normal())(c9)\n    # c9 = BatchNormalization()(c9)\n\n    outputs = Conv2D(1, (1, 1), activation='linear')(c9)\n    # outputs = UpSampling2D((2, 2))(c10)\n\n\n    model = Model(inputs=[inputs], outputs=[outputs])\n    return model","metadata":{"id":"x_BpwQ5qyj-q"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Divide the data in training and testing. We use 10% for testing purposes.\nx_train, x_test, y_train, y_test = train_test_split(train[:,:,:], labels[:,:,:], test_size = 0.1)\n\n# We increase the dimensions by 1 to be consisted with the tensorflow format\nx_train = x_train[...,np.newaxis]\ny_train = y_train[...,np.newaxis]\nx_test = x_test[...,np.newaxis]\ny_test = y_test[...,np.newaxis]\n\n# The size of the Bscans i.e. the inputs of the U-net\ninput_size = (230, 230, 1)\n# Initialise the unet\nmodel = unet_model(input_size)\n\n# You can play with the loss, instead of using a typucal mse or mae you can add e.g. a total variational constrain to force the model to be blocky\ndef custom_loss(y_true, y_pred):\n    # Mean Absolute Error\n    mae = tf.math.reduce_mean(tf.math.abs(y_true - y_pred))\n\n    # Mean Absolute Differences\n    # reg=tf.math.reduce_mean(tf.math.abs(y_pred[1:]-y_pred[:-1]))\n\n    # return mae + reg*0.001\n\n    return mae\n\n\n\n# Generate 10 models to apply ensemble U-nets\nfor i in range(0,10):\n    model = unet_model(input_size)\n    # Compile the model\n    model.compile(optimizer='adam', loss=custom_loss, metrics=['mae'])\n\n    # You can play with this parameter. Usually learning rate of adam should decrease as the iterations increase\n    model.optimizer.learning_rate = 5e-4\n\n    # The model with the best performance in the validation set will be saved in best_model_2_1.hp5\n    p1 = EarlyStopping(monitor=\"val_loss\", mode= \"min\", verbose=1, patience=30)\n    p2 = ModelCheckpoint('best_model_2_{}.hp5'.format(i), monitor=\"val_loss\", mode=\"min\", save_best_only=True )\n\n    # Train the model with 5 batch size for 15 epochs\n    history = model.fit(x_train, y_train, epochs=55, batch_size=5, validation_split=0.1, callbacks=[p1, p2])\n\n    # Decrease the learning rate and train for another 10 epochs\n    model.optimizer.learning_rate = 1e-4\n    history = model.fit(x_train, y_train, epochs=30, batch_size=5, validation_split=0.1, callbacks=[p1, p2])\n\n    # Decrease the learning rate and traing for another 10 epochs\n    model.optimizer.learning_rate = 5e-5\n    history = model.fit(x_train, y_train, epochs=30, batch_size=5, validation_split=0.1, callbacks=[p1, p2])\n\n\n\n","metadata":{"id":"NagrQBoOyzu_","outputId":"63650af7-8b5e-45b4-8285-05c02052bd34"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from mpl_toolkits.axes_grid1 import make_axes_locatable\n\n\nsaved_model = load_model('best_model_2_0.hp5', custom_objects= {'custom_loss': custom_loss})\n\n\n# Predictions using x_test. x_test are unknown data that were not included in the training process\npredictions = saved_model.predict(x_test)\n\n\n# Load the best model\nfor i in range(1,10):\n    saved_model = load_model('best_model_2_{}.hp5'.format(i), custom_objects= {'custom_loss': custom_loss})\n\n\n    # Predictions using x_test. x_test are unknown data that were not included in the training process\n    predictions = predictions + saved_model.predict(x_test)\n\npredictions = predictions/10\n\n# plot the results\nfor i, c in enumerate(predictions):\n    sh = np.shape(c)\n\n    # Actual size of the model i.e. 50x22 meters\n    xc = np.linspace(0,50,sh[0])\n    yc = np.linspace(0,22,sh[1])\n\n    plt.subplot(1,2,1)\n    ax = plt.gca()\n    c = (c+0.5)*10\n    im = ax.imshow(c, extent=[0, 50, -22, 0],vmin = c.min(), vmax = c.max())\n    plt.title('ML Inversion')\n    divider = make_axes_locatable(ax)\n    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n    plt.colorbar(im, cax=cax)\n    plt.xlabel('X-axis (m)')\n    plt.ylabel('Depth (m)')\n\n    plt.subplot(1,2,2)\n    ax = plt.gca()\n    im = ax.imshow((y_test[i,:,:]+0.5)*10, extent=[0, 50, -22, 0], vmin = c.min(), vmax = c.max())\n    plt.title('Ground_Truth')\n\n    divider = make_axes_locatable(ax)\n    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n    plt.colorbar(im, cax=cax)\n    plt.xlabel('X-axis (m)')\n    plt.ylabel('Depth (m)')\n    plt.tight_layout()\n\n    plt.show()\n","metadata":{"id":"EiuAiQFH1DUV","outputId":"3c600452-01e4-4300-a9ca-9aedcd0f46d2"},"execution_count":null,"outputs":[]}]}