{"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":"This is a Tensorflow implementation of the [Vesuvius Challenge: Ink Detection tutorial](https://www.kaggle.com/code/jpposma/vesuvius-challenge-ink-detection-tutorial). ","metadata":{}},{"cell_type":"code","source":"import gc\nimport os\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nimport PIL.Image as Image\n\nimport tensorflow as tf","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-03-23T13:22:13.527989Z","iopub.execute_input":"2023-03-23T13:22:13.528468Z","iopub.status.idle":"2023-03-23T13:22:25.440383Z","shell.execute_reply.started":"2023-03-23T13:22:13.528428Z","shell.execute_reply":"2023-03-23T13:22:25.438802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.random.set_seed(10)\nnp.random.seed(10)\n\npath = '/kaggle/input/vesuvius-challenge-ink-detection/train/1'\ndepth = (27, 37)\nradius = 30\nrect = ((1100, 3500), 700, 950)\nmax_steps = 30_000\nlr = 0.03\nbatch_size = 32","metadata":{"execution":{"iopub.status.busy":"2023-03-23T13:22:25.442771Z","iopub.execute_input":"2023-03-23T13:22:25.44355Z","iopub.status.idle":"2023-03-23T13:22:25.451578Z","shell.execute_reply.started":"2023-03-23T13:22:25.443507Z","shell.execute_reply":"2023-03-23T13:22:25.450139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create dataset","metadata":{}},{"cell_type":"code","source":"def select_pixels(path: str, radius: int, rect: tuple = None, include: bool = False) -> (np.ndarray):\n    \"\"\"\n    Select pixels.\n    \n    :param path: Location of the data.\n    :param radius: Size of the area around the pixel.\n    :param depth: Range of slices used.\n    :param rect: Rectangular area defined by the tuple ((x0, y0), width, height). \n    :param include: Whether to return the pixels inside (True) or outside (False) the rectangular area.\n    \"\"\"\n    \n    mask = np.array(Image.open(os.path.join(path, 'mask.png')), dtype='int')\n    height, width = mask.shape\n    \n    p = np.argwhere(mask == 1)\n    p = p[(p[:, 0] > radius) & (p[:, 1] > radius) & (p[:, 0] + radius + 1 < height - radius) & (p[:, 1] + radius + 1 < width)]\n    \n    if rect is None:\n        return p\n    \n    (xmin, ymin), width, height = rect\n    xmax, ymax = (xmin + width), (ymin + height)\n    \n    if include:\n        return p[(p[:, 1] > xmin) & (p[:, 1] < xmax) & (p[:, 0] > ymin) & (p[:, 0] < ymax)]\n    \n    return p[(p[:, 1] < xmin) | (p[:, 1] > xmax) | (p[:, 0] < ymin) | (p[:, 0] > ymax)]\n\n\ndef load_data(path: str, depth: tuple):\n    \"\"\"\n    Load data from disc.\n    \n    :param path: Location of the data.\n    :param depth: Range of slices used.\n    \"\"\"\n    \n    targets = np.array(Image.open(os.path.join(path, 'inklabels.png')), dtype='int')\n    height, width = targets.shape\n    \n    images = []\n    \n    for d in range(*depth):\n        image = np.array(Image.open(f'{path}/surface_volume/{d :02}.tif'), dtype='float32').reshape(height, width, 1)\n        images.append(image / 65535.0)\n        \n    images = np.stack(images)\n    \n    return images, targets\n\n\nclass Dataset(tf.keras.utils.Sequence):\n    \n    def __init__(\n        self, \n        images: np.ndarray, \n        targets: np.ndarray, \n        pixels: np.ndarray, \n        radius: int, \n        batch_size: int, \n        max_steps: int = None, \n        shuffle: bool = False\n    ):\n        \"\"\"\n        :param images: Array of images.\n        :param targets: Array of targets.\n        :param pixels: Array of pixels.\n        :param radius: Size of the area around the pixel.\n        :param batch_size: Batch size.\n        :param max_steps: Total number of batches.\n        :param shuffle: Wheter to shuffle the data.\n        \"\"\"\n        \n        self.images = images\n        self.targets = targets\n        self.pixels = pixels\n        self.radius = radius\n        self.batch_size = batch_size\n        \n        if shuffle:\n            np.random.shuffle(self.pixels)\n            \n        if max_steps is not None:\n            self.pixels = self.pixels[:self.batch_size * max_steps]\n            \n        \n    def _shape(self):\n        return (self.images.shape[0], int(self.radius * 2 + 1), int(self.radius * 2 + 1), 1)\n\n        \n    def __getitem__(self, index):\n        \n        radius = self.radius\n        pixels = self.pixels[(index * self.batch_size):((index + 1) * self.batch_size)]\n        \n        images, targets = [], []\n        \n        for p in pixels:\n            i, j = p[0], p[1]\n            \n            images.append(self.images[:, (i - radius):(i + radius + 1), (j - radius):(j + radius + 1), :])\n            targets.append(self.targets[i, j])\n        \n        return np.array(images), np.array(targets)\n    \n    def __len__(self):\n        return np.ceil(self.pixels.shape[0] / self.batch_size).astype('int')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-23T13:22:25.453769Z","iopub.execute_input":"2023-03-23T13:22:25.454194Z","iopub.status.idle":"2023-03-23T13:22:25.478094Z","shell.execute_reply.started":"2023-03-23T13:22:25.454116Z","shell.execute_reply":"2023-03-23T13:22:25.476618Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"images, targets = load_data(path, depth)\n\ntrain_pixels = select_pixels(path, radius, rect, include=False)\ntrain = Dataset(images, targets, train_pixels, radius, batch_size=batch_size, shuffle=True, max_steps=max_steps)\n\nvalid_pixels = select_pixels(path, radius, rect, include=True)\nvalid = Dataset(images, targets, valid_pixels, radius, batch_size=batch_size, shuffle=False)\n\n\nfig, ax = plt.subplots()\nax.imshow(targets, cmap='gray')\nax.add_patch(patches.Rectangle(*rect, edgecolor='r', facecolor='none'))\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-23T13:22:25.480773Z","iopub.execute_input":"2023-03-23T13:22:25.481233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train model","metadata":{}},{"cell_type":"code","source":"class OneCycleScheduler(tf.keras.callbacks.Callback):\n    \"\"\"\n    https://www.avanwyk.com/tensorflow-2-super-convergence-with-the-1cycle-policy/\n    \"\"\"\n    \n    def __init__(self, lr_max, steps, mom_min=0.85, mom_max=0.95, phase_1_pct=0.3, div_factor=25.):\n        super(OneCycleScheduler, self).__init__()\n        \n        lr_min = lr_max / div_factor\n        final_lr = lr_max / (div_factor * 1e4)\n        phase_1_steps = steps * phase_1_pct\n        phase_2_steps = steps - phase_1_steps\n        \n        self.phase_1_steps = phase_1_steps\n        self.phase_2_steps = phase_2_steps\n        self.phase = 0\n        self.step = 0\n        \n        self.phases = [\n            [CosineAnnealer(lr_min, lr_max, phase_1_steps), CosineAnnealer(mom_max, mom_min, phase_1_steps)], \n            [CosineAnnealer(lr_max, final_lr, phase_2_steps), CosineAnnealer(mom_min, mom_max, phase_2_steps)]\n        ]\n        \n        self.lrs = []\n        self.moms = []\n\n    def on_train_begin(self, logs=None):\n        self.phase = 0\n        self.step = 0\n\n        self.set_lr(self.lr_schedule().start)\n        self.set_momentum(self.mom_schedule().start)\n        \n    def on_train_batch_begin(self, batch, logs=None):\n        self.lrs.append(self.get_lr())\n        self.moms.append(self.get_momentum())\n\n    def on_train_batch_end(self, batch, logs=None):\n        self.step += 1\n        if self.step >= self.phase_1_steps:\n            self.phase = 1\n            \n        self.set_lr(self.lr_schedule().step())\n        self.set_momentum(self.mom_schedule().step())\n        \n    def get_lr(self):\n        try:\n            return tf.keras.backend.get_value(self.model.optimizer.lr)\n        except AttributeError:\n            return None\n        \n    def get_momentum(self):\n        try:\n            return tf.keras.backend.get_value(self.model.optimizer.momentum)\n        except AttributeError:\n            return None\n        \n    def set_lr(self, lr):\n        try:\n            tf.keras.backend.set_value(self.model.optimizer.lr, lr)\n        except AttributeError:\n            pass # ignore\n        \n    def set_momentum(self, mom):\n        try:\n            tf.keras.backend.set_value(self.model.optimizer.momentum, mom)\n        except AttributeError:\n            pass # ignore\n\n    def lr_schedule(self):\n        return self.phases[self.phase][0]\n    \n    def mom_schedule(self):\n        return self.phases[self.phase][1]\n    \n    \nclass CosineAnnealer:\n    \n    def __init__(self, start, end, steps):\n        self.start = start\n        self.end = end\n        self.steps = steps\n        self.n = 0\n        \n    def step(self):\n        self.n += 1\n        cos = np.cos(np.pi * (self.n / self.steps)) + 1\n        return self.end + (self.start - self.end) / 2. * cos\n    \n    \nclass LogBatchMetricsCallback(tf.keras.callbacks.Callback):\n    \n    def __init__(self):\n        super().__init__()\n        self._metrics = {'loss': [], 'lr': []}\n        \n    \n    def on_train_batch_end(self, batch, logs=None): \n        self._metrics['loss'].append(logs['loss'])\n        self._metrics['lr'].append(logs['lr'])\n        \n        \nclass TrackLRCallback(tf.keras.callbacks.Callback):\n    \n    def __init__(self, optimizer):\n        super().__init__()\n        self.optimizer = optimizer\n    \n    def on_batch_end(self, batch, logs):\n        logs.update({'lr' : self.optimizer.lr.numpy()})\n        \n    def on_epoch_end(self, epoch, logs):\n        logs.update({'lr' : self.optimizer.lr.numpy()})","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.keras.backend.clear_session()\n\nmodel = tf.keras.models.Sequential([\n    tf.keras.layers.InputLayer(input_shape=train._shape()),\n    tf.keras.layers.Conv3D(filters=16, kernel_size=3, activation=None, padding='same'),\n    tf.keras.layers.MaxPool3D(),\n    tf.keras.layers.Conv3D(filters=32, kernel_size=3, activation=None, padding='same'),\n    tf.keras.layers.MaxPool3D(),\n    tf.keras.layers.Conv3D(filters=64, kernel_size=3, activation=None, padding='same'),\n    tf.keras.layers.MaxPool3D(),\n    tf.keras.layers.Flatten(),\n    tf.keras.layers.Dense(128, 'relu'),\n    tf.keras.layers.Dense(1, 'sigmoid'),\n])\n\nmodel.summary()\n\nloss = tf.keras.losses.BinaryCrossentropy()\noptimizer = tf.keras.optimizers.SGD(lr)\n\ntrack_lr = TrackLRCallback(optimizer)\nbatch_metrics = LogBatchMetricsCallback()\nlr_schedule = OneCycleScheduler(lr, max_steps)\n\nmodel.compile(loss=loss, optimizer=optimizer)\n_ = model.fit(train, callbacks=[track_lr, batch_metrics, lr_schedule])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = pd.DataFrame(batch_metrics._metrics)\n\nfig, ax = plt.subplots(ncols=2, figsize=(12, 4))\nhistory.loss.plot(ax=ax[0]).set_title('loss')\nhistory.lr.plot(ax=ax[1]).set_title('lr')\nfig.tight_layout()\nfig.show()\n\n\ndel train\n_ = gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predict","metadata":{}},{"cell_type":"code","source":"y_pred = model.predict(valid)\n\ndel valid, model\n_ = gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, width, height = rect\n\nfig, ax = plt.subplots(figsize=(4, 4))\nax.imshow(y_pred.reshape((height - 1), (width - 1)), cmap='gray')\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}