{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Submission notebook","metadata":{}},{"cell_type":"markdown","source":"This notebook takes a pretrained Keras model (UNet with an EfficientNet-B0 backend) and builds the submission file. The process to create the submission is:\n\n1. Get the prediction from the model, output is sigmoid, so we need to find a good threshold.\n2. Find the best \"threshold-kernel size\" relation to improve the score on the three original fragments. \"Kernel size\" is used in step 4.\n3. Apply the threshold to the prediction.\n4. Apply classical image processing morphological \"open-close\" transformations to improve the quality of the predicted labels.\n5. Resize to original size.\n6. Generate submission file.\n\n4 different models are provided, which can be selected using the VAL_FOLD config:\n\n* \"1\": model was trained with the entire volume 1 as validation data\n* \"2\": model was trained with the entire volume 2 as validation data\n* \"3\": model was trained with the entire volume 3 as validation data\n* \"CUSTOM\": Validation data is build from 3 random patches taken from the middle of each of the 3 original volumes. The rest of the volumes were used as training data\n\nI am making this public for three reasons:\n\n1. Using classical open-close transformations to improve the predicted label perhaps can help people squeeze some decimals in the public leaderboard score.\n2. Perhaps someone wants to play around with model ensembling.\n3. The public leaderboard for these 4 models is extremely low compared with the score obtained on the validation set. I am still trying to understand if this is a bug on the submission code, or an overfitting problem in the model.","metadata":{}},{"cell_type":"code","source":"import os\nimport seaborn as sns\nimport pandas as pd\nfrom sklearn.metrics import fbeta_score","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:46:47.185978Z","iopub.execute_input":"2023-04-23T18:46:47.186503Z","iopub.status.idle":"2023-04-23T18:46:47.194249Z","shell.execute_reply.started":"2023-04-23T18:46:47.186460Z","shell.execute_reply":"2023-04-23T18:46:47.192327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PROD = False\n\nfrom tensorflow.keras.layers import Conv2D, BatchNormalization, Activation, MaxPool2D, Conv2DTranspose, Concatenate, Input\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.applications import EfficientNetB0\nimport tensorflow as tf\nfrom tensorflow import keras\nimport keras.backend as K\nfrom keras import layers\nfrom keras.utils import get_file\n\nfrom tensorflow.keras.metrics import Metric\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nimport PIL.Image as Image\n\nimport random\nimport gc\nimport cv2\nimport time\nfrom tqdm import tqdm\nimport re\nimport math\nfrom collections import namedtuple\nfrom io import StringIO\nimport os\n\ntf.keras.utils.set_random_seed(1234)\n\nDATA_DIR = \"/kaggle/input/vesuvius-challenge-ink-detection\"\nPATCH_SIZE = 256  # e.g. 128x128\nPATCH_HALFSIZE = PATCH_SIZE // 2\nDOWNSAMPLING = 0.5 # Setting this to e.g. 0.5 means images will be loaded as 2x smaller. 1 does nothing.\nZ_DIM = 8   # Number of slices in the z direction. Max value is 65 - Z_START\nZ_START = 28  # Offset of slices in the z direction\nBATCH_SIZE = 8\nTHRESHOLD = 0.5\nLEARNING_RATE = 0.0001\nWARMUP_EPOCHS = 1\nDECAY_EPOCHS = 125\n\nCONFIG = {\n    \"PATCH_SIZE\": PATCH_SIZE,\n    \"PATCH_HALFSIZE\": PATCH_HALFSIZE,\n    \"DOWNSAMPLING\": DOWNSAMPLING,\n    \"Z_DIM\": Z_DIM,\n    \"Z_START\": Z_START,\n    \"BATCH_SIZE\": BATCH_SIZE,\n    \"learning_rate\": LEARNING_RATE,\n    \"epochs\": 50,\n    \"steps_per_epoch\": 100,\n    \"THRESHOLD\": THRESHOLD,\n    \"LOSS\": \"MODIFIED_DICE\",\n    \"VAL_FOLD\": \"CUSTOM\",\n    \"SIGMOID_OUTPUT\": True,\n    'model_name': 'efficientnetb0-unet',\n    'DATA_AUG': True,\n    'WARMUP_EPOCHS': WARMUP_EPOCHS,\n    'DECAY_EPOCHS': DECAY_EPOCHS\n}\n\nprint(\"TF Version: \", tf.__version__)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:46:47.196617Z","iopub.execute_input":"2023-04-23T18:46:47.197053Z","iopub.status.idle":"2023-04-23T18:46:47.252585Z","shell.execute_reply.started":"2023-04-23T18:46:47.196993Z","shell.execute_reply":"2023-04-23T18:46:47.251213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data import","metadata":{}},{"cell_type":"code","source":"def resize(img,downsampling=DOWNSAMPLING):\n    if downsampling != 1.:\n        size = int(img.shape[1] * downsampling), int(img.shape[0] * downsampling)\n        img = cv2.resize(img, size)\n    return img\n\ndef resize_to_original(img,original_mask):\n    size = original_mask.shape[1], original_mask.shape[0]\n    img = cv2.resize(img, size)\n    return img\n\ndef load_mask(split, index):\n    img = cv2.imread(f\"{DATA_DIR}/{split}/{index}/mask.png\", 0)\n    img = resize(img)\n    return img.astype(\"bool\")\n\n\ndef load_labels(split, index):\n    img = cv2.imread(f\"{DATA_DIR}/{split}/{index}/inklabels.png\", 0)\n    img = resize(img)\n    return np.expand_dims(img, axis=-1)\n\n\ndef load_volume(split, index):\n    # A more memory-efficient volune loader\n    fnames = [f\"{DATA_DIR}/{split}/{index}/surface_volume/{i:02}.tif\"\n             for i in range(Z_START, Z_START + Z_DIM)]\n\n    batch_size = 8\n    fname_batches = [fnames[i :i + batch_size] for i in range(0, len(fnames), batch_size)]\n    volumes = []\n    for fname_batch in fname_batches:\n        z_slices = []\n        for fname in tqdm(fname_batch):\n            img = cv2.imread(fname, 0)\n            img = resize(img)\n            z_slices.append(img)\n        volumes.append(np.stack(z_slices, axis=-1))\n        del z_slices\n    return np.concatenate(volumes, axis=-1)\n\n\ndef load_sample(split, index):\n    print(f\"Loading '{split}/{index}'...\")\n    gc.collect()\n    if split == \"train\":\n        return load_volume(split, index), load_mask(split, index), load_labels(split, index)\n    return load_volume(split, index), load_mask(split, index), None","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:46:47.253924Z","iopub.execute_input":"2023-04-23T18:46:47.255167Z","iopub.status.idle":"2023-04-23T18:46:47.272032Z","shell.execute_reply.started":"2023-04-23T18:46:47.255114Z","shell.execute_reply":"2023-04-23T18:46:47.270194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"volume_1, mask_1, labels_1 = load_sample(split=\"train\", index=1)\nvolume_2, mask_2, labels_2 = load_sample(split=\"train\", index=2)\nvolume_3, mask_3, labels_3 = load_sample(split=\"train\", index=3)\ngc.collect()\nprint(\"Loading complete.\")","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:46:47.274952Z","iopub.execute_input":"2023-04-23T18:46:47.275393Z","iopub.status.idle":"2023-04-23T18:47:10.044981Z","shell.execute_reply.started":"2023-04-23T18:46:47.275363Z","shell.execute_reply":"2023-04-23T18:47:10.044086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prod_data  = {\n    \"train_volumes\": [volume_1, volume_2, volume_3],\n    \"train_labels\": [labels_1, labels_2, labels_3],\n    \"train_masks\": [mask_1, mask_2, mask_3],\n}","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:10.046137Z","iopub.execute_input":"2023-04-23T18:47:10.046942Z","iopub.status.idle":"2023-04-23T18:47:10.052497Z","shell.execute_reply.started":"2023-04-23T18:47:10.046906Z","shell.execute_reply":"2023-04-23T18:47:10.050840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Utilities to create dataset","metadata":{}},{"cell_type":"code","source":"def sample_random_location(shape):\n    x = random.randint(PATCH_HALFSIZE, shape[0] - PATCH_HALFSIZE - 1)\n    y = random.randint(PATCH_HALFSIZE, shape[1] - PATCH_HALFSIZE - 1)\n    return (x, y)\n\n\ndef list_all_locations(mask, stride=PATCH_SIZE):\n    locations = []\n    for x in range(PATCH_HALFSIZE, mask.shape[0] - PATCH_HALFSIZE, stride):\n        for y in range(PATCH_HALFSIZE, mask.shape[1] - PATCH_HALFSIZE, stride):\n            if mask[x, y]:  \n                locations.append((x, y))\n    return locations\n\n\ndef extract_patch(location, volume):\n    x = location[0]\n    y = location[1]\n    patch = volume[x - PATCH_HALFSIZE :x + PATCH_HALFSIZE,\n                   y - PATCH_HALFSIZE :y + PATCH_HALFSIZE, :]\n    return patch.astype(\"float32\") / 255.\n\n\ndef extract_labels(location, labels):\n    x = location[0]\n    y = location[1]\n    \n    label = labels[x - PATCH_HALFSIZE :x + PATCH_HALFSIZE,\n                   y - PATCH_HALFSIZE :y + PATCH_HALFSIZE, :]\n    return label.astype(\"float32\") / 255.\n\n\ndef make_random_data_generator(volume, mask, labels):\n    def data_generator():\n        while True:\n            loc = sample_random_location(mask.shape)\n            x = loc[0]\n            y = loc[1]\n            if mask[x, y]:    \n                patch = extract_patch(loc, volume)\n                label = extract_labels(loc, labels)\n                yield patch, label\n    return data_generator\n\n\ndef make_iterated_data_generator(volume, mask, labels=None, return_locations=False):\n    locations = list_all_locations(mask)\n    def data_generator():\n        for loc in locations:\n            patch = extract_patch(loc, volume)\n            if labels is None:\n                if return_locations:\n                    yield patch, loc\n                else:\n                    yield patch\n            else:\n                label = extract_labels(loc, labels)\n                if return_locations:\n                    yield patch, label, loc\n                else:\n                    yield patch, label\n    return data_generator\n\ndef make_random_data_generator_from_list(volume_list, mask_list, labels_list):\n    def data_generator():\n        while True:\n            dataset_idx = random.randint(0,len(volume_list)-1)\n            loc = sample_random_location(mask_list[dataset_idx].shape)\n            x = loc[0]\n            y = loc[1]\n            if mask_list[dataset_idx][x, y]:    \n                patch = extract_patch(loc, volume_list[dataset_idx])\n                label = extract_labels(loc, labels_list[dataset_idx])\n                yield patch, label\n    return data_generator\n    \ndef make_iterated_data_generator_from_list(volume_list, mask_list, labels_list=None, return_locations=False):\n    locations = []\n    list_length = len(mask_list)\n    for i in range(list_length):\n        locations.append(list_all_locations(mask_list[i]))\n    def data_generator():\n        for i in range(list_length):\n            for loc in locations[i]:\n                patch = extract_patch(loc, volume_list[i])\n                if labels_list is None:\n                    if return_locations:\n                        yield patch, loc\n                    else:\n                        yield patch\n                else:\n                    label = extract_labels(loc, labels_list[i])\n                    if return_locations:\n                        yield patch, label, loc\n                    else:\n                        yield patch, label\n    return data_generator\n\ndef make_tf_dataset(gen_fn, labeled=True, return_locations=False):\n    if labeled:\n        if return_locations:\n            output_signature = (\n                tf.TensorSpec(shape=(PATCH_SIZE, PATCH_SIZE, Z_DIM), dtype=tf.float32),\n                tf.TensorSpec(shape=(PATCH_SIZE, PATCH_SIZE, 1), dtype=tf.float32),\n                tf.TensorSpec(shape=(2,), dtype=tf.float32),\n            )\n        else:\n            output_signature = (\n                tf.TensorSpec(shape=(PATCH_SIZE, PATCH_SIZE, Z_DIM), dtype=tf.float32),\n                tf.TensorSpec(shape=(PATCH_SIZE, PATCH_SIZE, 1), dtype=tf.float32),\n            )\n    else:\n        if return_locations:\n            output_signature = (\n                tf.TensorSpec(shape=(PATCH_SIZE, PATCH_SIZE, Z_DIM), dtype=tf.float32),\n                tf.TensorSpec(shape=(2,), dtype=tf.float32),\n            )\n        else:\n            output_signature = tf.TensorSpec(shape=(PATCH_SIZE, PATCH_SIZE, Z_DIM), dtype=tf.float32)\n    ds = tf.data.Dataset.from_generator(\n        gen_fn,\n        output_signature=output_signature,\n    )\n    return ds.prefetch(tf.data.AUTOTUNE).batch(BATCH_SIZE)\n\ndef make_datasets_for_fold(fold, train_augment_fn=None, return_locations=False):\n    train_volumes = fold[\"train_volumes\"]\n    train_masks = fold[\"train_masks\"]\n    train_labels = fold[\"train_labels\"]\n    \n    include_validation = \"validation_volume\" in fold\n    if include_validation:\n        validation_volume = fold[\"validation_volume\"]\n        validation_mask = fold[\"validation_mask\"]\n        validation_labels = fold[\"validation_labels\"]\n\n    train_ds = make_tf_dataset(\n        make_random_data_generator_from_list(train_volumes, train_masks, train_labels),\n        labeled=True\n    )\n    \n    if train_augment_fn:\n        train_ds = train_ds.map(train_augment_fn, num_parallel_calls=tf.data.AUTOTUNE)\n    train_ds = train_ds.prefetch(tf.data.AUTOTUNE)\n\n    if not include_validation:\n        return train_ds\n\n    val_ds = make_tf_dataset(\n            make_iterated_data_generator_from_list(validation_volume, validation_mask, validation_labels, return_locations=return_locations),\n            labeled=True,\n            return_locations=return_locations\n        )\n    \n    return train_ds, val_ds","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:10.054367Z","iopub.execute_input":"2023-04-23T18:47:10.055336Z","iopub.status.idle":"2023-04-23T18:47:10.086778Z","shell.execute_reply.started":"2023-04-23T18:47:10.055291Z","shell.execute_reply":"2023-04-23T18:47:10.085487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model definition","metadata":{}},{"cell_type":"code","source":"class Swish(layers.Layer):\n    def __init__(self, name=None, **kwargs):\n        super().__init__(name=name, **kwargs)\n\n    def call(self, inputs, **kwargs):\n        return tf.nn.swish(inputs)\n\n    def get_config(self):\n        config = super().get_config()\n        config['name'] = self.name\n        return config\n    \nclass DropConnect(layers.Layer):\n\n    def __init__(self, drop_connect_rate, **kwargs):\n        super().__init__(**kwargs)\n        self.drop_connect_rate = drop_connect_rate\n\n    def call(self, inputs, **kwargs):\n        def drop_connect():\n            keep_prob = 1.0 - self.drop_connect_rate\n\n            # Compute drop_connect tensor\n            batch_size = tf.shape(inputs)[0]\n            random_tensor = keep_prob\n            random_tensor += tf.random.uniform([batch_size, 1, 1, 1], dtype=inputs.dtype)\n            binary_tensor = tf.floor(random_tensor)\n            output = tf.math.divide(inputs, keep_prob) * binary_tensor\n            return output\n\n        return K.in_train_phase(drop_connect(), inputs, training=None)\n\n    def get_config(self):\n        config = super().get_config()\n        config['drop_connect_rate'] = self.drop_connect_rate\n        return config","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:10.088535Z","iopub.execute_input":"2023-04-23T18:47:10.088883Z","iopub.status.idle":"2023-04-23T18:47:10.105421Z","shell.execute_reply.started":"2023-04-23T18:47:10.088851Z","shell.execute_reply":"2023-04-23T18:47:10.104004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def binary_fbeta(ytrue , ypred, beta=1, epsilon=1e-7):\n    # epsilon is set so as to avoid division by zero error\n    beta_squared = beta**2 # squaring beta\n\n    # casting ytrue and ypred as float dtype\n    ytrue = tf.cast(ytrue, tf.float32)\n    ypred = tf.cast(ypred, tf.float32)\n\n    tp = tf.reduce_sum(ytrue*ypred) # calculating true positives\n    predicted_positive = tf.reduce_sum(ypred) # calculating predicted positives\n    actual_positive = tf.reduce_sum(ytrue) # calculating actual positives\n    \n    precision = tp/(predicted_positive+epsilon) # calculating precision\n    recall = tp/(actual_positive+epsilon) # calculating recall\n    \n    # calculating fbeta\n    fb = (1+beta_squared)*precision*recall / (beta_squared*precision + recall + epsilon)\n\n    return fb\n\ndef modified_dice_loss(y_true, y_pred, beta=0.5, epsilon=1e-7):\n    \n    beta_squared = beta**2 # squaring beta\n    \n    y_true_f = K.flatten(y_true)\n    if CONFIG[\"SIGMOID_OUTPUT\"]:\n        y_pred_f = K.flatten(y_pred)\n    else:\n        y_pred_f = tf.keras.activations.sigmoid(K.flatten(y_pred))\n    \n    tp = tf.reduce_sum(y_true_f*y_pred_f) # calculating true positives\n    predicted_positive = tf.reduce_sum(y_pred_f) # calculating predicted positives\n    actual_positive = tf.reduce_sum(y_true_f) # calculating actual positives\n    \n    precision = tp/(predicted_positive+epsilon) # calculating precision\n    recall = tp/(actual_positive+epsilon) # calculating recall\n    \n    fb = (1+beta_squared)*precision*recall / (beta_squared*precision + recall + epsilon)\n    \n    return 1-fb","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:10.106977Z","iopub.execute_input":"2023-04-23T18:47:10.107500Z","iopub.status.idle":"2023-04-23T18:47:10.125175Z","shell.execute_reply.started":"2023-04-23T18:47:10.107455Z","shell.execute_reply":"2023-04-23T18:47:10.124176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class StatefullBinaryFBeta(Metric):\n    def __init__(self, name='state_full_binary_fbeta', beta=0.5, threshold=THRESHOLD, epsilon=1e-7, **kwargs):\n        # initializing an object of the super class\n        super(StatefullBinaryFBeta, self).__init__(name=name, **kwargs)\n\n        # initializing state variables\n        self.tp = self.add_weight(name='tp', initializer='zeros') # initializing true positives \n        self.actual_positive = self.add_weight(name='fp', initializer='zeros') # initializing actual positives\n        self.predicted_positive = self.add_weight(name='fn', initializer='zeros') # initializing predicted positives\n\n        # initializing other atrributes that wouldn't be changed for every object of this class\n        self.beta_squared = beta**2 \n        self.threshold = threshold\n        self.epsilon = epsilon\n\n    def update_state(self, ytrue, ypred, sample_weight=None):\n        # casting ytrue and ypred as float dtype\n        ytrue = tf.cast(ytrue, tf.float32)\n        if CONFIG[\"SIGMOID_OUTPUT\"]:\n            ypred = tf.cast(ypred, tf.float32)\n        else:\n            ypred = tf.keras.activations.sigmoid(tf.cast(ypred, tf.float32))\n\n        # setting values of ypred greater than the set threshold to 1 while those lesser to 0\n        ypred = tf.cast(tf.greater_equal(ypred, tf.constant(self.threshold)), tf.float32)\n\n        self.tp.assign_add(tf.reduce_sum(ytrue*ypred)) # updating true positives atrribute\n        self.predicted_positive.assign_add(tf.reduce_sum(ypred)) # updating predicted positive atrribute\n        self.actual_positive.assign_add(tf.reduce_sum(ytrue)) # updating actual positive atrribute\n\n    def result(self):\n        self.precision = self.tp/(self.predicted_positive+self.epsilon) # calculates precision\n        self.recall = self.tp/(self.actual_positive+self.epsilon) # calculates recall\n\n        # calculating fbeta\n        self.fb = (1+self.beta_squared)*self.precision*self.recall / (self.beta_squared*self.precision + self.recall + self.epsilon)\n\n        return self.fb\n\n    def reset_state(self):\n        self.tp.assign(0) # resets true positives to zero\n        self.predicted_positive.assign(0) # resets predicted positives to zero\n        self.actual_positive.assign(0) # resets actual positives to zero","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:10.126455Z","iopub.execute_input":"2023-04-23T18:47:10.127304Z","iopub.status.idle":"2023-04-23T18:47:10.142874Z","shell.execute_reply.started":"2023-04-23T18:47:10.127245Z","shell.execute_reply":"2023-04-23T18:47:10.141444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"custom_objects = {\n    \"Swish\": Swish, \n    \"DropConnect\": DropConnect, \n    \"modified_dice_loss\": modified_dice_loss,\n    \"StatefullBinaryFBeta\": StatefullBinaryFBeta\n}\nwith keras.utils.custom_object_scope(custom_objects):\n    model_name = CONFIG[\"VAL_FOLD\"].lower()\n    model = keras.models.load_model(f\"/kaggle/input/model-baseline/model_{model_name}.keras\")","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:10.146382Z","iopub.execute_input":"2023-04-23T18:47:10.147575Z","iopub.status.idle":"2023-04-23T18:47:13.403667Z","shell.execute_reply.started":"2023-04-23T18:47:10.147509Z","shell.execute_reply":"2023-04-23T18:47:13.402300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Sanity check","metadata":{}},{"cell_type":"code","source":"thresholds = np.arange(0.1,1,0.05)\nkernel_sizes = range(1,50,2)\n\ndef grid_search_th_kernel(predicted_labels,true_labels):\n    results = []\n    for th in tqdm(thresholds):\n        for k in kernel_sizes:\n            if CONFIG[\"SIGMOID_OUTPUT\"]:\n                predicted_labels_th = np.where(predicted_labels > th, 1, 0).astype(np.uint8)\n            else:\n                predicted_labels_th = np.where(tf.keras.activations.sigmoid(predicted_labels) > th, 1, 0).astype(np.uint8)\n            kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (k,k))\n            predicted_labels_closed = cv2.morphologyEx(predicted_labels_th, cv2.MORPH_CLOSE, kernel)\n            beta_score = float(binary_fbeta(ytrue=true_labels,ypred=np.expand_dims(predicted_labels_closed,axis=-1),beta=0.5))\n            results.append({\n                'threshold': th,\n                'kernel_size': k,\n                'beta_score': beta_score\n            })\n    df = pd.DataFrame(results)\n    return df\n\ndef best_th_kernel_comb(df):\n    \n    best_threshold = 0\n    best_kernel_size = 0\n    \n    x = []\n    y = []\n    z = []\n    \n    for th in tqdm(thresholds):\n        for k in kernel_sizes:\n            x.append(th)\n            y.append(k)\n            mean_beta_score = df[(df['threshold'] == th) & (df['kernel_size'] == k)]['beta_score'].mean()\n            z.append(mean_beta_score)\n                \n    data = pd.DataFrame(data={'x':x, 'y':y, 'z':z})\n    data = data.pivot(index='x', columns='y', values='z')\n    sns.heatmap(data)\n    plt.show()\n    \n    max_idx = np.argmax(z)\n    best_threshold = x[max_idx]\n    best_kernel_size = y[max_idx]\n    \n    return best_threshold, best_kernel_size","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:13.405969Z","iopub.execute_input":"2023-04-23T18:47:13.406405Z","iopub.status.idle":"2023-04-23T18:47:13.420571Z","shell.execute_reply.started":"2023-04-23T18:47:13.406361Z","shell.execute_reply":"2023-04-23T18:47:13.418856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_best = None\n\nfor i in range(len(prod_data[\"train_volumes\"])):\n    print(\"===============================\")\n    print(f\"Training volume {i}:\")\n    train_ds = make_tf_dataset(\n            make_iterated_data_generator(\n                prod_data[\"train_volumes\"][i], \n                prod_data[\"train_masks\"][i],\n                prod_data[\"train_labels\"][i],\n                return_locations=True),\n            labeled=True,\n            return_locations=True\n    )\n        \n    predicted_labels = np.zeros(prod_data[\"train_volumes\"][i].shape[:2] + (1,), dtype=\"float32\")\n    predictions_map_counts = np.zeros(prod_data[\"train_volumes\"][i].shape[:2] + (1,), dtype=\"int8\")\n    true_labels = np.zeros(prod_data[\"train_labels\"][i].shape[:2] + (1,), dtype=\"float32\")\n    for patch_batch, y_true_batch, loc_batch in tqdm(train_ds):\n        predictions = model.predict_on_batch(patch_batch)\n        for (x, y), pred, y_true in zip(loc_batch, predictions, y_true_batch):\n            x = int(x.numpy())\n            y = int(y.numpy())\n            predicted_labels[x-PATCH_HALFSIZE:x+PATCH_HALFSIZE,y-PATCH_HALFSIZE:y+PATCH_HALFSIZE,:] += pred\n            predictions_map_counts[x - PATCH_HALFSIZE : x + PATCH_HALFSIZE, y - PATCH_HALFSIZE : y + PATCH_HALFSIZE, :] += 1  \n            true_labels[x-PATCH_HALFSIZE:x+PATCH_HALFSIZE,y-PATCH_HALFSIZE:y+PATCH_HALFSIZE,:] += y_true\n\n    print(f\"\\tMin/Max predictions_map_counts : {np.min(predictions_map_counts)}/{np.max(predictions_map_counts)}\")\n\n    predicted_labels /= (predictions_map_counts + 1e-7)\n    true_labels /= (predictions_map_counts + 1e-7)\n    true_labels = np.where(true_labels > THRESHOLD, 1, 0)\n\n    print(f\"\\tMin/Max predicted_labels : {np.min(predicted_labels)}/{np.max(predicted_labels)}\")\n    fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4)\n    ax1.set_title(f\"True \")\n    ax1.imshow(true_labels, cmap='gray')\n    ax2.set_title(f\"Pred\")\n    ax2.imshow(predicted_labels, cmap='gray')\n    \n    df = grid_search_th_kernel(predicted_labels,true_labels)\n    if df_best is None:\n        df_best = df\n    else:\n        df_best = pd.concat([df_best,df])\n    \n    best_row = df.iloc[df['beta_score'].idxmax()]\n    best_threshold = best_row['threshold']\n    best_kernel_size = int(best_row['kernel_size'])\n    \n    ax3.set_title(f\"th = {best_threshold:.2f}\")\n    if CONFIG[\"SIGMOID_OUTPUT\"]:\n        predicted_labels_th = np.where(predicted_labels > best_threshold, 1, 0).astype(np.uint8)\n    else:\n        predicted_labels_th = np.where(tf.keras.activations.sigmoid(predicted_labels) > best_threshold, 1, 0).astype(np.uint8)\n    print(f\"\\tMin/Max predicted_labels_th : {np.min(predicted_labels_th)}/{np.max(predicted_labels_th)}\")\n    ax3.imshow(predicted_labels_th, cmap='gray')\n\n    beta_score = binary_fbeta(ytrue=true_labels,ypred=predicted_labels_th,beta=0.5)\n    \n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (best_kernel_size,best_kernel_size))\n    predicted_labels_closed = cv2.morphologyEx(predicted_labels_th, cv2.MORPH_CLOSE, kernel)\n    beta_score_closed = binary_fbeta(ytrue=true_labels,ypred=np.expand_dims(predicted_labels_closed,axis=-1),beta=0.5)\n\n    ax4.set_title(f\"k = {best_kernel_size}\")\n    ax4.imshow(predicted_labels_closed, cmap='gray')\n\n    orig_labels = cv2.imread(DATA_DIR + f\"/train/{i+1}/inklabels.png\", 0)\n    orig_labels = np.asarray(orig_labels) / 255.0\n    predicted_labels_closed_orig_size = resize_to_original(predicted_labels_closed,orig_labels)\n    beta_score_orig_size = binary_fbeta(ytrue=orig_labels,ypred=predicted_labels_closed_orig_size,beta=0.5)\n    beta_score_sklearn = fbeta_score(y_true=K.flatten(orig_labels),y_pred=K.flatten(predicted_labels_closed_orig_size),beta=0.5)\n\n    print(f\"Beta score (DOWNSAMPLING = {DOWNSAMPLING}): {beta_score}\")\n    print(f\"Beta score (CLOSED): {beta_score_closed}\")\n    print(f\"Beta score (ORIGINAL SIZE): {beta_score_orig_size}\")\n    print(f\"Beta score (sklearn): {beta_score_sklearn}\")\n\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-23T18:47:13.422488Z","iopub.execute_input":"2023-04-23T18:47:13.422948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Selecting best threshold and kernel size for morphological closing...\")\nbest_threshold, best_kernel_size = best_th_kernel_comb(df_best)\nprint(f\"Best threshold: {best_threshold}\")\nprint(f\"Best kernel size: {best_kernel_size}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predicting test","metadata":{}},{"cell_type":"code","source":"def post_process(predicted_labels,threshold,index):\n    test_mask = Image.open(DATA_DIR + f\"/test/{index}/mask.png\")\n    test_mask = np.asarray(test_mask)\n    \n    assert np.max(test_mask) == 1\n    \n    if CONFIG[\"SIGMOID_OUTPUT\"]:\n        pred_th = np.where(predicted_labels > threshold, 1, 0).astype(np.uint8)\n    else:\n        pred_th = np.where(tf.keras.activations.sigmoid(predicted_labels) > threshold, 1, 0).astype(np.uint8)\n    \n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (best_kernel_size,best_kernel_size))\n    pred_th_closed = cv2.morphologyEx(pred_th, cv2.MORPH_CLOSE, kernel)\n    \n    pred_th_resized = resize_to_original(pred_th_closed,test_mask)\n    \n    pred_th_resized_masked = pred_th_resized*test_mask\n    \n    return pred_th_resized_masked\n    \ndef compute_test_output(model,index):\n    test_volume, test_mask, _ = load_sample(split=\"test\", index=index)\n    test_ds = make_tf_dataset(\n        make_iterated_data_generator(test_volume, test_mask,return_locations=True),\n        labeled=False,\n        return_locations=True\n    )\n    \n    predicted_labels = np.zeros(test_volume.shape[:2], dtype=\"float32\")\n    predictions_map_counts = np.zeros(test_volume.shape[:2], dtype=\"int8\")\n    \n    for patch_batch, loc_batch in tqdm(test_ds):\n        predictions = model.predict_on_batch(patch_batch)\n        for (x, y), pred in zip(loc_batch, predictions):\n            x = int(x.numpy())\n            y = int(y.numpy())\n            predicted_labels[x-PATCH_HALFSIZE:x+PATCH_HALFSIZE,y-PATCH_HALFSIZE:y+PATCH_HALFSIZE] += pred[:,:,0]\n            predictions_map_counts[x - PATCH_HALFSIZE : x + PATCH_HALFSIZE, y - PATCH_HALFSIZE : y + PATCH_HALFSIZE] += 1  \n    \n    predicted_labels /= (predictions_map_counts + 1e-7)\n    del test_volume\n    del test_ds\n    del predictions_map_counts\n    \n    gc.collect()\n    return predicted_labels\n\ndef rle(img):\n    '''\n    img: numpy array, 1 - mask, 0 - background\n    Returns run length as string formated\n    '''\n    pixels = img.flatten()\n    \n    pixels = np.concatenate([[0], pixels, [0]])\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 1\n    runs[1::2] -= runs[::2]\n    return ' '.join(str(x) for x in runs)\n\ndef fast_rle(img):\n    flat_img = img.flatten().astype(np.uint8)\n\n    starts = np.array((flat_img[:-1] == 0) & (flat_img[1:] == 1))\n    ends = np.array((flat_img[:-1] == 1) & (flat_img[1:] == 0))\n    starts_ix = np.where(starts)[0] + 2\n    ends_ix = np.where(ends)[0] + 2\n    lengths = ends_ix - starts_ix\n    predicted_arr = np.stack([starts_ix, lengths]).T.flatten()\n    f = StringIO()\n    np.savetxt(f, predicted_arr.reshape(1, -1), delimiter=\" \", fmt=\"%d\")\n    predicted = f.getvalue().strip()\n    return predicted","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"indexes = sorted(os.listdir(f\"{DATA_DIR}/test\"))\nindexes","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results = []\n\nfor index in indexes:\n    predicted = compute_test_output(model,index=index)\n    predicted_postprocess = post_process(predicted,threshold=best_threshold,index=index)\n    print(f\"Min predicted: {np.min(predicted)}\")\n    print(f\"Max predicted: {np.max(predicted)}\")\n    print(f\"Min predicted: {np.min(predicted_postprocess)}\")\n    print(f\"Max predicted: {np.max(predicted_postprocess)}\")\n    fig, (ax1, ax2) = plt.subplots(1, 2)\n    ax1.set_title(f\"Test: {index}\")\n    ax1.imshow(predicted, cmap='gray')\n    ax2.set_title(f\"Test: {index} with postprocess\")\n    ax2.imshow(predicted_postprocess, cmap='gray')\n    plt.show()\n    predicted_rle = rle(predicted_postprocess)\n    predicted_rle_fast = fast_rle(predicted_postprocess)\n    assert predicted_rle == predicted_rle_fast\n    results.append((index,predicted_rle))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = pd.DataFrame(results, columns=['Id', 'Predicted'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_sub = pd.read_csv(DATA_DIR + '/sample_submission.csv')\nsample_sub = pd.merge(sample_sub[['Id']], sub, on='Id', how='left')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_sub.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_sub.to_csv(\"submission.csv\", index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}