Datasets
Models
Downloads: 1
Diff of
/baseline.ipynb
[000000]
..
[bc293e]
Switch to side-by-side view
--- a +++ b/baseline.ipynb @@ -0,0 +1,1680 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# GI Tract Segmentation Competition" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load And Prepare" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:47.438209Z", + "start_time": "2022-06-28T03:46:44.722879Z" + } + }, + "outputs": [], + "source": [ + "from fastai.vision.all import *\n", + "import matplotlib.patches as mpatches\n", + "import albumentations as A\n", + "import cv2\n", + "import pynvml\n", + "from scipy.spatial.distance import directed_hausdorff\n", + "import timm\n", + "import segmentation_models_pytorch as smp\n", + "\n", + "\n", + "# !cp kaggle.json /home/kgeorgio/.kaggle\n", + "# !kaggle competitions download -c uw-madison-gi-tract-image-segmentation\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load CSV" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:47.672193Z", + "start_time": "2022-06-28T03:46:47.439755Z" + }, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "train_df = pd.read_csv('dataset/train.csv', low_memory=False)\n", + "train_df = train_df.pivot(index='id', columns='class', values='segmentation').reset_index()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:47.739865Z", + "start_time": "2022-06-28T03:46:47.730425Z" + }, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "<div>\n", + "<style scoped>\n", + " .dataframe tbody tr th:only-of-type {\n", + " vertical-align: middle;\n", + " }\n", + "\n", + " .dataframe tbody tr th {\n", + " vertical-align: top;\n", + " }\n", + "\n", + " .dataframe thead th {\n", + " text-align: right;\n", + " }\n", + "</style>\n", + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: right;\">\n", + " <th>class</th>\n", + " <th>id</th>\n", + " <th>large_bowel</th>\n", + " <th>small_bowel</th>\n", + " <th>stomach</th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <th>0</th>\n", + " <td>case101_day20_slice_0001</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " </tr>\n", + " <tr>\n", + " <th>1</th>\n", + " <td>case101_day20_slice_0002</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " </tr>\n", + " <tr>\n", + " <th>2</th>\n", + " <td>case101_day20_slice_0003</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " </tr>\n", + " <tr>\n", + " <th>3</th>\n", + " <td>case101_day20_slice_0004</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " </tr>\n", + " <tr>\n", + " <th>4</th>\n", + " <td>case101_day20_slice_0005</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " <td>NaN</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>\n", + "</div>" + ], + "text/plain": [ + "class id large_bowel small_bowel stomach\n", + "0 case101_day20_slice_0001 NaN NaN NaN\n", + "1 case101_day20_slice_0002 NaN NaN NaN\n", + "2 case101_day20_slice_0003 NaN NaN NaN\n", + "3 case101_day20_slice_0004 NaN NaN NaN\n", + "4 case101_day20_slice_0005 NaN NaN NaN" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "train_df.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:47.950258Z", + "start_time": "2022-06-28T03:46:47.742286Z" + }, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dataset/train/case101/case101_day26/scans/slice_0121_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0081_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0065_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0068_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0059_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0064_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0002_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0088_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0098_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0074_266_266_1.50_1.50.png\n", + "dataset/train/case101/case101_day26/scans/slice_0021_266_266_1.50_1.50.png\n", + "...\n", + "And 38486 more lines.\n" + ] + } + ], + "source": [ + "path = Path('dataset/train')\n", + "\n", + "fnames = get_image_files(path)\n", + "\n", + "for ind, file_name in enumerate(fnames):\n", + " print(file_name)\n", + " if ind>9:\n", + " print('...')\n", + " break\n", + "print(f\"And {len(fnames)-10} more lines.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Functions" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def get_slice_id(fname):\n", + " return fname.parts[3] + '_' + fname.parts[5][:10]\n", + "\n", + "\n", + "def rle_decode(mask_rle, shape, color=1):\n", + " \"\"\" TBD\n", + "\n", + " Args:\n", + " mask_rle (str): run-length as string formated (start length)\n", + " shape (tuple of ints): (height,width) of array to return \n", + "\n", + " Returns: \n", + " Mask (np.array)\n", + " - 1 indicating mask\n", + " - 0 indicating background\n", + "\n", + " \"\"\"\n", + " # Split the string by space, then convert it into a integer array\n", + " s = np.array(mask_rle.split(), dtype=int)\n", + "\n", + " # Every even value is the start, every odd value is the \"run\" length\n", + " starts = s[0::2] - 1\n", + " lengths = s[1::2]\n", + " ends = starts + lengths\n", + "\n", + " # The image image is actually flattened since RLE is a 1D \"run\"\n", + " if len(shape) == 3:\n", + " h, w, d = shape\n", + " img = np.zeros((h * w, d), dtype=np.float32)\n", + " else:\n", + " h, w = shape\n", + " img = np.zeros((h * w,), dtype=np.float32)\n", + "\n", + " # The color here is actually just any integer you want!\n", + " for lo, hi in zip(starts, ends):\n", + " img[lo: hi] = color\n", + "\n", + " # Don't forget to change the image back to the original shape\n", + " return img.reshape(shape)\n", + "\n", + "\n", + "def label_func(fname):\n", + " # First we need to get the slice row\n", + " slice_id = get_slice_id(fname)\n", + " slice_row = train_df.query('id == @slice_id')\n", + "\n", + " # Then we need to extract the slice width and height which are provided in the fname last part\n", + " # Typically the height is the first part of a dimension, but for some reason the slices have\n", + " # widths provided first\n", + " w, h = map(lambda x: int(x), fname.parts[-1].split('_')[2:4])\n", + "\n", + " # Create mask array (It needs to have 3 channels but fastai will only keep the first one anyways)\n", + " mask = np.zeros((h, w, 3), dtype=np.uint8)\n", + "\n", + " # If the segmentation mask is str\n", + " # Each mask should have it's own code (color) where fastai will use them for identification\n", + " if isinstance(slice_row['large_bowel'].item(), str):\n", + " mask[:, :, 0] = rle_decode(slice_row['large_bowel'].item(), shape=(h, w), color=255)\n", + "\n", + " if isinstance(slice_row['small_bowel'].item(), str):\n", + " mask[:, :, 1] = rle_decode(slice_row['small_bowel'].item(), shape=(h, w), color=255)\n", + "\n", + " if isinstance(slice_row['stomach'].item(), str):\n", + " mask[:, :, 2] = rle_decode(slice_row['stomach'].item(), shape=(h, w), color=255)\n", + "\n", + " return mask\n", + "\n", + "\n", + "# This was the code available in fastai\n", + "@ToTensor\n", + "def encodes(self, o: PILMask): return o._tensor_cls(image2tensor(o)[0])\n", + "\n", + "\n", + "# And this is how we customize it to suit our needs\n", + "@ToTensor\n", + "def encodes(self, o: PILMask): return o._tensor_cls(image2tensor(o))\n", + "\n", + "\n", + "@typedispatch\n", + "def show_batch(x: TensorImage, y: TensorMask, samples, ctxs=None, max_n=6, nrows=None, ncols=2,\n", + " figsize=None, **kwargs):\n", + " if figsize is None: figsize = (ncols * 3, max_n // ncols * 3)\n", + " if ctxs is None: ctxs = get_grid(max_n, nrows=nrows, ncols=ncols, figsize=figsize)\n", + " for i, ctx in enumerate(ctxs):\n", + " x_i = x[i] / x[i].max()\n", + " show_image(x_i, ctx=ctx, cmap='gray', **kwargs)\n", + " show_image(y[i], ctx=ctx, cmap='Spectral_r', alpha=0.35, **kwargs)\n", + " red_patch = mpatches.Patch(color='red', label='large_bowel')\n", + " green_patch = mpatches.Patch(color='green', label='small_bowel')\n", + " blue_patch = mpatches.Patch(color='blue', label='stomach')\n", + " ctx.legend(handles=[red_patch, green_patch, blue_patch], fontsize=figsize[0] / 2)\n", + "\n", + "\n", + "def pad_img(img, up_size=None):\n", + " if up_size is None:\n", + " return img\n", + " shape0 = np.array(img.shape[:2])\n", + " resize = np.array(up_size)\n", + " if np.any(shape0 != resize):\n", + " diff = resize - shape0\n", + " pad0 = diff[0]\n", + " pad1 = diff[1]\n", + " pady = [pad0 // 2, pad0 // 2 + pad0 % 2]\n", + " padx = [pad1 // 2, pad1 // 2 + pad1 % 2]\n", + " img = np.pad(img, [pady, padx])\n", + " img = img.reshape((*resize))\n", + " return img\n", + "\n", + "\n", + "def unpad_img(img, up_size, org_size):\n", + " shape0 = np.array(org_size)\n", + " resize = np.array(up_size)\n", + " if np.any(shape0 != resize):\n", + " diff = resize - shape0\n", + " pad0 = diff[0]\n", + " pad1 = diff[1]\n", + " pady = [pad0 // 2, pad0 // 2 + pad0 % 2]\n", + " padx = [pad1 // 2, pad1 // 2 + pad1 % 2]\n", + " img = img[pady[0]:-pady[1], padx[0]:-padx[1], :]\n", + " img = img.reshape((*shape0, 3))\n", + " return img\n", + "\n", + "\n", + "def load_image(fname, up_size=None):\n", + " img = np.array(Image.open(fname))\n", + " img = np.interp(img, [np.min(img), np.max(img)], [0, 255])\n", + " return pad_img(img, up_size)\n", + "\n", + "\n", + "def get_25D_image(row, up_size=None):\n", + " if up_size:\n", + " imgs = np.zeros((*up_size, len(row['fnames'])))\n", + " else:\n", + " imgs = np.zeros((row['slice_h'], row['slice_w'], len(row['fnames'])))\n", + "\n", + " for i, fname in enumerate(row['fnames']):\n", + " img = load_image(fname, up_size)\n", + " imgs[..., i] += img\n", + " return imgs.astype(np.uint8)\n", + "\n", + "\n", + "def get_mask(row, up_size=None):\n", + " if up_size:\n", + " mask = np.zeros((*up_size, 3))\n", + " else:\n", + " mask = np.zeros((row['slice_h'], row['slice_w'], 3))\n", + "\n", + " if isinstance(row['large_bowel'], str):\n", + " mask[..., 0] += pad_img(\n", + " rle_decode(row['large_bowel'], shape=(row['slice_h'], row['slice_w']), color=255), up_size)\n", + " if isinstance(row['small_bowel'], str):\n", + " mask[..., 1] += pad_img(\n", + " rle_decode(row['small_bowel'], shape=(row['slice_h'], row['slice_w']), color=255), up_size)\n", + " if isinstance(row['stomach'], str):\n", + " mask[..., 2] += pad_img(\n", + " rle_decode(row['stomach'], shape=(row['slice_h'], row['slice_w']), color=255), up_size)\n", + "\n", + " return mask.astype(np.uint8)\n", + "\n", + "\n", + "def get_train_aug(img_size, crop=0.9, p=0.4):\n", + " crop_size = round(img_size[0] * crop)\n", + " return A.Compose([\n", + " A.RandomCrop(height=crop_size, width=crop_size, always_apply=True),\n", + " A.HorizontalFlip(p=p),\n", + " A.OneOf([\n", + " A.GridDistortion(num_steps=5, distort_limit=0.05, p=1.0),\n", + " A.ElasticTransform(alpha=1, sigma=50, alpha_affine=50, p=1.0)\n", + " ], p=p),\n", + " A.CoarseDropout(\n", + " max_holes=8, min_holes=8,\n", + " max_height=crop_size // 10, max_width=crop_size // 10,\n", + " min_height=4, min_width=4, mask_fill_value=0, p=0.2 * p),\n", + " A.ShiftScaleRotate(\n", + " shift_limit=0.0625, scale_limit=0.2, rotate_limit=25,\n", + " interpolation=cv2.INTER_AREA, p=p),\n", + " A.HorizontalFlip(p=0.5 * p),\n", + " A.OneOf([\n", + " A.MotionBlur(p=0.2 * p),\n", + " A.MedianBlur(blur_limit=3, p=0.1 * p),\n", + " A.Blur(blur_limit=3, p=0.1 * p),\n", + " ], p=0.2 * p),\n", + " A.GaussNoise(var_limit=0.001, p=0.2 * p),\n", + " A.OneOf([\n", + " A.OpticalDistortion(p=0.3 * p),\n", + " A.GridDistortion(p=0.1 * p),\n", + " A.PiecewiseAffine(p=0.3 * p),\n", + " ], p=0.2 * p),\n", + " A.OneOf([\n", + " A.Sharpen(p=0.2 * p),\n", + " A.Emboss(p=0.2 * p),\n", + " A.RandomBrightnessContrast(p=0.2 * p),\n", + " ]),\n", + " ])\n", + "\n", + "\n", + "def get_test_aug(img_size, crop=0.9):\n", + " crop_size = round(crop * img_size[0])\n", + " return A.Compose([\n", + " A.CenterCrop(height=crop_size, width=crop_size),\n", + " ])\n", + "\n", + "\n", + "class AlbumentationsTransform(ItemTransform, RandTransform):\n", + " split_idx, order = None, 2\n", + "\n", + " def __init__(self, train_aug, valid_aug):\n", + " store_attr()\n", + "\n", + " def before_call(self, b, split_idx):\n", + " self.idx = split_idx\n", + "\n", + " def encodes(self, x):\n", + " if len(x) > 1:\n", + " img, mask = x\n", + " if self.idx == 0:\n", + " aug = self.train_aug(image=np.array(img), mask=np.array(mask))\n", + " else:\n", + " aug = self.valid_aug(image=np.array(img), mask=np.array(mask))\n", + " return PILImage.create(aug[\"image\"]), PILMask.create(aug[\"mask\"])\n", + " else:\n", + " img = x[0]\n", + " aug = self.valid_aug(image=np.array(img))\n", + " return PILImage.create(aug[\"image\"])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Preprocess" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Split name metada into different columns" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:48.134095Z", + "start_time": "2022-06-28T03:46:47.951449Z" + }, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "train_df['partial_fname'] = train_df.id\n", + "fname_df = pd.DataFrame({'partial_fname': [f'{fname.parts[-3]}_slice_{fname.parts[-1][6:10]}' for fname in fnames],\n", + " 'fname': fnames})\n", + "\n", + "train_df = train_df.merge(fname_df, on='partial_fname').drop('partial_fname', axis=1)\n", + "\n", + "train_df['case_id'] = train_df.id.apply(lambda x: x.split('_')[0])\n", + "train_df['day_num'] = train_df.id.apply(lambda x: x.split('_')[1])\n", + "\n", + "train_df['slice_w'] = train_df[\"fname\"].apply(lambda x: int(str(x)[:-4].rsplit(\"_\",4)[1]))\n", + "train_df['slice_h'] = train_df[\"fname\"].apply(lambda x: int(str(x)[:-4].rsplit(\"_\",4)[2]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Add Fnames" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:48.186750Z", + "start_time": "2022-06-28T03:46:48.135151Z" + }, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "channels = 3\n", + "stride = 2\n", + "for j, i in enumerate(range(-1*(channels-channels//2-1), channels//2+1)):\n", + " method = 'ffill'\n", + " if i <= 0: method = 'bfill'\n", + " train_df[f'fname_{j:02}'] = train_df.groupby(['case_id', 'day_num'])['fname'].shift(stride*-i).fillna(method=method)\n", + " \n", + "train_df['fnames'] = train_df[[f'fname_{j:02d}' for j in range(channels)]].values.tolist()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "source": [ + "# Training" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "source": [ + "## Metrics" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:49.086943Z", + "start_time": "2022-06-28T03:46:49.086930Z" + } + }, + "outputs": [], + "source": [ + "def dice_coeff_adj(inp, targ):\n", + " inp = np.where(sigmoid(inp).cpu().detach().numpy() > 0.5, 1, 0)\n", + " targ = targ.cpu().detach().numpy()\n", + " eps = 1e-5\n", + " dice_scores = []\n", + " for i in range(targ.shape[0]):\n", + " dice_i = []\n", + " for j in range(targ.shape[1]):\n", + " if inp[i, j].sum() == targ[i, j].sum() == 0:\n", + " continue\n", + " I = (targ[i, j] * inp[i, j]).sum()\n", + " U = targ[i, j].sum() + inp[i, j].sum()\n", + " dice_i.append((2.*I)/(U+eps))\n", + " if dice_i:\n", + " dice_scores.append(np.mean(dice_i))\n", + " \n", + " if dice_scores:\n", + " return np.mean(dice_scores)\n", + " else:\n", + " return 0\n", + " \n", + " \n", + "def hd_dist_per_slice(inp, targ, seed): \n", + " inp = np.argwhere(inp) / np.array(inp.shape)\n", + " targ = np.argwhere(targ) / np.array(targ.shape)\n", + " haussdorf_dist = 1 - directed_hausdorff(inp, targ, seed)[0]\n", + " return haussdorf_dist if haussdorf_dist > 0 else 0\n", + "\n", + "def hd_dist_adj(inp, targ, seed=42):\n", + " inp = np.where(sigmoid(inp).cpu().detach().numpy() > 0.5, 1, 0)\n", + " targ = targ.cpu().detach().numpy()\n", + " hd_scores = []\n", + " for i in range(targ.shape[0]):\n", + " hd_i = []\n", + " for j in range(targ.shape[1]):\n", + " if inp[i, j].sum() == targ[i, j].sum() == 0:\n", + " continue\n", + " hd_i.append(hd_dist_per_slice(inp[i, j], targ[i, j], seed))\n", + " if hd_i:\n", + " hd_scores.append(np.mean(hd_i))\n", + " if hd_scores:\n", + " return np.mean(hd_scores)\n", + " else:\n", + " return 0\n", + "\n", + "def custom_metric_adj(inp, targ, seed=42):\n", + " hd_score_per_batch = hd_dist_adj(inp, targ, seed)\n", + " dice_score_per_batch = dice_coeff_adj(inp, targ)\n", + " \n", + " return 0.4*dice_score_per_batch + 0.6*hd_score_per_batch" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Loss Functions" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:49.088279Z", + "start_time": "2022-06-28T03:46:49.088271Z" + } + }, + "outputs": [], + "source": [ + "class DiceBCEModule(Module):\n", + " def __init__(self, eps:float=1e-5, from_logits=True):\n", + " store_attr()\n", + " \n", + " def forward(self, inp:Tensor, targ:Tensor) -> Tensor:\n", + " inp = inp.view(-1)\n", + " targ = targ.view(-1)\n", + " \n", + " if self.from_logits: \n", + " bce_loss = nn.BCEWithLogitsLoss()(inp, targ)\n", + " inp = torch.sigmoid(inp)\n", + " \n", + " \n", + " intersection = (inp * targ).sum() \n", + " dice = (2.*intersection + self.eps)/(inp.sum() + targ.sum() + self.eps) \n", + " \n", + " return 0.5*(1 - dice) + 0.5*bce_loss\n", + "\n", + "\n", + "class DiceBCELoss(BaseLoss):\n", + " def __init__(self, *args, eps:float=1e-5, from_logits=True, thresh=0.5, **kwargs):\n", + " super().__init__(DiceBCEModule, *args, eps=eps, from_logits=from_logits, flatten=False, is_2d=True, floatify=True, **kwargs)\n", + " self.thresh = thresh\n", + " \n", + " def decodes(self, x:Tensor) -> Tensor:\n", + " \"Converts model output to target format\"\n", + " return (x>self.thresh).long()\n", + "\n", + " def activation(self, x:Tensor) -> Tensor:\n", + " \"`nn.BCEWithLogitsLoss`'s fused activation function applied to model output\"\n", + " return torch.sigmoid(x)\n", + "\n", + "# Source: https://www.kaggle.com/code/thedrcat/focal-multilabel-loss-in-pytorch-explained/notebook\n", + "def focal_binary_cross_entropy(logits, targets, gamma=2, n=3):\n", + " p = torch.sigmoid(logits)\n", + " p = torch.where(targets >= 0.5, p, 1-p)\n", + " logp = - torch.log(torch.clamp(p, 1e-4, 1-1e-4))\n", + " loss = logp*((1-p)**gamma)\n", + " loss = n*loss.mean()\n", + " return loss\n", + "\n", + "class DiceFocalModule(Module):\n", + " def __init__(self, eps:float=1e-5, from_logits=True, ws=[0.5, 0.5], gamma=2, n=3):\n", + " store_attr()\n", + " \n", + " def forward(self, inp:Tensor, targ:Tensor) -> Tensor:\n", + " inp = inp.view(-1)\n", + " targ = targ.view(-1)\n", + " \n", + " if self.from_logits: \n", + " focal_loss = focal_binary_cross_entropy(inp, targ, self.gamma, self.n)\n", + " inp = torch.sigmoid(inp)\n", + " \n", + " \n", + " intersection = (inp * targ).sum() \n", + " dice = (2.*intersection + self.eps)/(inp.sum() + targ.sum() + self.eps) \n", + " \n", + " return self.ws[0]*(1 - dice) + self.ws[1]*focal_loss\n", + " \n", + "class DiceFocalLoss(BaseLoss):\n", + " def __init__(self, *args, eps:float=1e-5, from_logits=True, ws=[0.5, 0.5], gamma=2, n=3, thresh=0.5, **kwargs):\n", + " super().__init__(DiceFocalModule, *args, eps=eps, from_logits=from_logits, ws=ws, gamma=gamma, n=n, flatten=False, is_2d=True, floatify=True, **kwargs)\n", + " self.thresh = thresh\n", + " \n", + " def decodes(self, x:Tensor) -> Tensor:\n", + " \"Converts model output to target format\"\n", + " return (x>self.thresh).long()\n", + "\n", + " def activation(self, x:Tensor) -> Tensor:\n", + " \"`nn.BCEWithLogitsLoss`'s fused activation function applied to model output\"\n", + " return torch.sigmoid(x)\n", + "\n", + "class FocalTverskyLossModule(Module):\n", + " def __init__(self, eps:float=1e-5, from_logits=True, alpha=0.3, beta=0.7, gamma=3/4):\n", + " store_attr()\n", + " \n", + " def forward(self, inp:Tensor, targ:Tensor) -> Tensor:\n", + " inp = inp.view(-1)\n", + " targ = targ.view(-1)\n", + " \n", + " if self.from_logits: \n", + " inp = torch.sigmoid(inp)\n", + " \n", + " inp_0, inp_1 = inp, 1 - inp\n", + " targ_0, targ_1 = targ, 1 - targ\n", + " \n", + " num = (inp_0 * targ_0).sum() \n", + " denom = num + (self.alpha * (inp_0 * targ_1).sum()) + (self.beta * (inp_1 * targ_0).sum()) + self.eps\n", + " loss = 1 - (num / denom)\n", + " return loss**self.gamma \n", + " \n", + "class FocalTverskyLoss(BaseLoss):\n", + " def __init__(self, *args, eps:float=1e-5, from_logits=True, alpha=0.3, beta=0.7, gamma=3/4, thresh=0.5, **kwargs):\n", + " super().__init__(FocalTverskyLossModule, *args, eps=eps, from_logits=from_logits, alpha=alpha, beta=beta, gamma=gamma, flatten=False, is_2d=True, floatify=True, **kwargs)\n", + " self.thresh = thresh\n", + " \n", + " def decodes(self, x:Tensor) -> Tensor:\n", + " \"Converts model output to target format\"\n", + " return (x>self.thresh).long()\n", + "\n", + " def activation(self, x:Tensor) -> Tensor:\n", + " \"`nn.BCEWithLogitsLoss`'s fused activation function applied to model output\"\n", + " return torch.sigmoid(x)\n", + "\n", + "def focal_binary_cross_entropy(logits, targets, gamma=2, n=3):\n", + " p = torch.sigmoid(logits)\n", + " p = torch.where(targets >= 0.5, p, 1-p)\n", + " logp = - torch.log(torch.clamp(p, 1e-4, 1-1e-4))\n", + " loss = logp*((1-p)**gamma)\n", + " loss = n*loss.mean()\n", + " return loss\n", + "\n", + "class ComboModule(Module):\n", + " def __init__(self, eps:float=1e-5, from_logits=True, ws=[2, 3, 1], gamma=2, n=3):\n", + " store_attr()\n", + " \n", + " def forward(self, inp:Tensor, targ:Tensor) -> Tensor:\n", + " inp = inp.view(-1)\n", + " targ = targ.view(-1)\n", + " \n", + " if self.from_logits: \n", + " focal_loss = focal_binary_cross_entropy(inp, targ, self.gamma, self.n)\n", + " bce_loss = nn.BCEWithLogitsLoss()(inp, targ)\n", + " inp = torch.sigmoid(inp)\n", + " \n", + " intersection = (inp * targ).sum() \n", + " dice = (2.*intersection + self.eps)/(inp.sum() + targ.sum() + self.eps) \n", + " \n", + " return self.ws[0]*(1 - dice) + self.ws[1]*focal_loss + self.ws[2]*bce_loss\n", + " \n", + "class ComboLoss(BaseLoss):\n", + " def __init__(self, *args, eps:float=1e-5, from_logits=True, ws=[2, 3, 1], gamma=2, n=3, thresh=0.5, **kwargs):\n", + " super().__init__(ComboModule, *args, eps=eps, from_logits=from_logits, ws=ws, gamma=gamma, n=n, flatten=False, is_2d=True, floatify=True, **kwargs)\n", + " self.thresh = thresh\n", + " \n", + " def decodes(self, x:Tensor) -> Tensor:\n", + " \"Converts model output to target format\"\n", + " return (x>self.thresh).long()\n", + "\n", + " def activation(self, x:Tensor) -> Tensor:\n", + " \"`nn.BCEWithLogitsLoss`'s fused activation function applied to model output\"\n", + " return torch.sigmoid(x)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set up CUDA" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "device_name: GeForce RTX 2080 Ti\n", + "device_capability: (7, 5)\n", + "device_properties: _CudaDeviceProperties(name='GeForce RTX 2080 Ti', major=7, minor=5, total_memory=11019MB, multi_processor_count=68)\n", + "current_device: 0\n", + "\n", + "Available Device IDs: (0, 1, 2, 3)\n", + "Device 0: 10620 MB free\n", + "Device 1: 7899 MB free\n", + "Device 2: 16274 MB free\n", + "Device 3: 9241 MB free\n" + ] + } + ], + "source": [ + "# Select Cuda GPU device\n", + "def get_memory_free_MiB(gpu_index):\n", + " pynvml.nvmlInit()\n", + " handle = pynvml.nvmlDeviceGetHandleByIndex(int(gpu_index))\n", + " mem_info = pynvml.nvmlDeviceGetMemoryInfo(handle)\n", + " return mem_info.free // 1024 ** 2\n", + "\n", + "device_count = torch.cuda.device_count()\n", + "current_device = torch.cuda.current_device()\n", + "device_name = torch.cuda.get_device_name(current_device)\n", + "device_capability = torch.cuda.get_device_capability(current_device)\n", + "device_properties = torch.cuda.get_device_properties(current_device)\n", + "is_available = torch.cuda.is_available()\n", + "device_cuda = torch.device(\"cuda\")\n", + "devices_tup = tuple(range(device_count))\n", + "\n", + "print('device_name: {device_name}'.format(device_name=device_name))\n", + "print('device_capability: {device_capability}'.format(device_capability=device_capability))\n", + "print('device_properties: {device_properties}'.format(device_properties=device_properties))\n", + "print('current_device: {current_device}'.format(current_device=current_device))\n", + "print(\"\\nAvailable Device IDs: \", devices_tup)\n", + "\n", + "for device_id in devices_tup:\n", + " print(f\"Device {device_id}: {get_memory_free_MiB(device_id)} MB free\")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Successfully selected device 2\n" + ] + } + ], + "source": [ + "device_id = 2\n", + "torch.backends.cudnn.benchmark = True\n", + "\n", + "if device_id is not None:\n", + " torch.cuda.set_device(device_id)\n", + " current_device = torch.cuda.current_device()\n", + " if current_device == device_id:\n", + " print(f\"\\nSuccessfully selected device {device_id}\")\n", + " else:\n", + " print(f\"Error: Couldn't change device from {current_device} to {device_id}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Baseline train using SMP" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:49.089093Z", + "start_time": "2022-06-28T03:46:49.089085Z" + } + }, + "outputs": [], + "source": [ + "def build_model(encoder_name, in_c=3, classes=3, weights=\"imagenet\"):\n", + " model = smp.Unet(\n", + " encoder_name=encoder_name, \n", + " encoder_weights=weights, \n", + " in_channels=in_c, \n", + " classes=classes, \n", + " activation=None\n", + " )\n", + " return model\n", + "\n", + "# Split any of model parameters from smp into encoder and decoder, \n", + "# so that we can freeze and unfreeze encoder layers.\n", + "def smp_splitter(model):\n", + " model_layers = list(model.children())\n", + " encoder_params = params(model_layers[0])\n", + " decoder_params = params(model_layers[1]) + params(model_layers[2])\n", + " return L(encoder_params, decoder_params)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/kgeorgio/miniconda3/envs/gi-tract/lib/python3.7/site-packages/torch/_tensor.py:575: UserWarning: floor_divide is deprecated, and will be removed in a future version of pytorch. It currently rounds toward 0 (like the 'trunc' function NOT 'floor'). This results in incorrect rounding for negative values.\n", + "To keep the current behavior, use torch.div(a, b, rounding_mode='trunc'), or for actual floor division, use torch.div(a, b, rounding_mode='floor'). (Triggered internally at /opt/conda/conda-bld/pytorch_1631630742027/work/aten/src/ATen/native/BinaryOps.cpp:467.)\n", + " return torch.floor_divide(self, other)\n" + ] + } + ], + "source": [ + "dev_df = train_df.sample(frac=0.2)\n", + "up_size = (320, 384)\n", + "tfms = [[partial(get_25D_image, up_size=up_size), PILImage.create],\n", + " [partial(get_mask, up_size=up_size), PILMask.create]] # the pipeline\n", + "\n", + "splits = RandomSplitter()(dev_df)\n", + "# https://docs.fast.ai/tutorial.albumentations.html\n", + "albu_aug = AlbumentationsTransform(get_train_aug(up_size), \n", + " get_test_aug(up_size))\n", + " \n", + "# https://docs.fast.ai/data.core.html#Datasets\n", + "dsets = Datasets(train_df, tfms, splits=splits)\n", + "# https://docs.fast.ai/data.load.html#DataLoader\n", + "dls = dsets.dataloaders(bs=16, \n", + " after_item=[albu_aug, ToTensor], \n", + " after_batch=[IntToFloatTensor(div_mask=255), \n", + " Normalize.from_stats(*imagenet_stats)],\n", + " device=device_id)\n", + "# https://smp.readthedocs.io/en/latest/encoders.html\n", + "model = build_model('efficientnet-b0', in_c=3, classes=3, weights=\"imagenet\")\n", + "model = model.cuda(device_id)\n", + "metrics = [\n", + " dice_coeff_adj, \n", + " hd_dist_adj, \n", + " custom_metric_adj\n", + " ]\n", + "loss_func = ComboLoss()\n", + "splitter = smp_splitter" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "ExecuteTime": { + "end_time": "2022-06-28T03:46:49.091707Z", + "start_time": "2022-06-28T03:46:49.091699Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "<style>\n", + " /* Turns off some styling */\n", + " progress {\n", + " /* gets rid of default border in Firefox and Opera. */\n", + " border: none;\n", + " /* Needs to be in here for Safari polyfill so background images work as expected. */\n", + " background-size: auto;\n", + " }\n", + " .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n", + " background: #F44336;\n", + " }\n", + "</style>\n" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: left;\">\n", + " <th>epoch</th>\n", + " <th>train_loss</th>\n", + " <th>valid_loss</th>\n", + " <th>dice_coeff_adj</th>\n", + " <th>hd_dist_adj</th>\n", + " <th>custom_metric_adj</th>\n", + " <th>time</th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <td>0</td>\n", + " <td>1.383521</td>\n", + " <td>1.162622</td>\n", + " <td>0.383626</td>\n", + " <td>0.797801</td>\n", + " <td>0.632131</td>\n", + " <td>01:41</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "<style>\n", + " /* Turns off some styling */\n", + " progress {\n", + " /* gets rid of default border in Firefox and Opera. */\n", + " border: none;\n", + " /* Needs to be in here for Safari polyfill so background images work as expected. */\n", + " background-size: auto;\n", + " }\n", + " .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n", + " background: #F44336;\n", + " }\n", + "</style>\n" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: left;\">\n", + " <th>epoch</th>\n", + " <th>train_loss</th>\n", + " <th>valid_loss</th>\n", + " <th>dice_coeff_adj</th>\n", + " <th>hd_dist_adj</th>\n", + " <th>custom_metric_adj</th>\n", + " <th>time</th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <td>0</td>\n", + " <td>0.814299</td>\n", + " <td>0.605600</td>\n", + " <td>0.621061</td>\n", + " <td>0.814855</td>\n", + " <td>0.737338</td>\n", + " <td>01:59</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# https://docs.fast.ai/learner.html#Learner\n", + "learn = Learner(dls, \n", + " model,\n", + " metrics=metrics, \n", + " loss_func=loss_func, \n", + " splitter=splitter).to_fp16()\n", + "\n", + "#https://docs.fast.ai/callback.schedule.html#Learner.fine_tune\n", + "learn.freeze()\n", + "learn.fine_tune(1, 1e-2)\n", + "learn.export('test_model.pkl')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "source": [ + "## Use DynamicUnet" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "def timm_model_sizes(encoder, img_size):\n", + " sizes = []\n", + " for layer in encoder.feature_info:\n", + " sizes.append(torch.Size([1, layer['num_chs'], img_size[0]//layer['reduction'], img_size[1]//layer['reduction']]))\n", + " return sizes\n", + "\n", + "\n", + "def get_timm_output_layers(encoder):\n", + " outputs = []\n", + " for layer in encoder.feature_info:\n", + " # Converts 'blocks.0.0' to ['blocks', '0', '0']\n", + " attrs = layer['module'].split('.')\n", + " output_layer = getattr(encoder, attrs[0])[int(attrs[1])][int(attrs[2])]\n", + " outputs.append(output_layer)\n", + " return outputs\n", + "\n", + "\n", + "class DynamicTimmUnet(SequentialEx):\n", + " \"Create a U-Net from a given architecture in timm.\"\n", + " def __init__(self, encoder, n_out, img_size, blur=False, blur_final=True, self_attention=False,\n", + " y_range=None, last_cross=True, bottle=False, act_cls=defaults.activation,\n", + " init=nn.init.kaiming_normal_, norm_type=None, **kwargs):\n", + " imsize = img_size\n", + " sizes = timm_model_sizes(encoder, img_size)\n", + " sz_chg_idxs = list(reversed(range(len(sizes))))\n", + " outputs = list(reversed(get_timm_output_layers(encoder)))\n", + " self.sfs = hook_outputs(outputs, detach=False)\n", + " \n", + " # cut encoder\n", + " encoder = nn.Sequential(*list(encoder.children()))[:-5]\n", + " \n", + " x = dummy_eval(encoder, imsize).detach()\n", + "\n", + " ni = sizes[-1][1]\n", + " middle_conv = nn.Sequential(ConvLayer(ni, ni*2, act_cls=act_cls, norm_type=norm_type, **kwargs),\n", + " ConvLayer(ni*2, ni, act_cls=act_cls, norm_type=norm_type, **kwargs)).eval()\n", + " x = middle_conv(x)\n", + " layers = [encoder, BatchNorm(ni), nn.ReLU(), middle_conv]\n", + "\n", + " for i,idx in enumerate(sz_chg_idxs):\n", + " not_final = i!=len(sz_chg_idxs)-1\n", + " up_in_c, x_in_c = int(x.shape[1]), int(sizes[idx][1])\n", + " do_blur = blur and (not_final or blur_final)\n", + " sa = self_attention and (i==len(sz_chg_idxs)-3)\n", + " unet_block = UnetBlock(up_in_c, x_in_c, self.sfs[i], final_div=not_final, blur=do_blur, self_attention=sa,\n", + " act_cls=act_cls, init=init, norm_type=norm_type, **kwargs).eval()\n", + " layers.append(unet_block)\n", + " x = unet_block(x)\n", + "\n", + " ni = x.shape[1]\n", + " if imsize != sizes[0][-2:]: layers.append(PixelShuffle_ICNR(ni, act_cls=act_cls, norm_type=norm_type))\n", + " layers.append(ResizeToOrig())\n", + " if last_cross:\n", + " layers.append(MergeLayer(dense=True))\n", + " ni += in_channels(encoder)\n", + " layers.append(ResBlock(1, ni, ni//2 if bottle else ni, act_cls=act_cls, norm_type=norm_type, **kwargs))\n", + " layers += [ConvLayer(ni, n_out, ks=1, act_cls=None, norm_type=norm_type, **kwargs)]\n", + " apply_init(nn.Sequential(layers[3], layers[-2]), init)\n", + " #apply_init(nn.Sequential(layers[2]), init)\n", + " if y_range is not None: layers.append(SigmoidRange(*y_range))\n", + " layers.append(ToTensorBase())\n", + " super().__init__(*layers)\n", + "\n", + " def __del__(self):\n", + " if hasattr(self, \"sfs\"): self.sfs.remove()\n", + " \n", + " \n", + "def dynamic_unet_splitter(model):\n", + " return L(model[0], model[1:]).map(params)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "<style>\n", + " /* Turns off some styling */\n", + " progress {\n", + " /* gets rid of default border in Firefox and Opera. */\n", + " border: none;\n", + " /* Needs to be in here for Safari polyfill so background images work as expected. */\n", + " background-size: auto;\n", + " }\n", + " .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n", + " background: #F44336;\n", + " }\n", + "</style>\n" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "SuggestedLRs(valley=0.00013182566908653826)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtpklEQVR4nO3deXxcdbn48c8zk8m+NUu3pCVt09J9syBYkNqC7AVFQC+uPxSRe0FQ4cq9ioD7dYGLioKieEXBWlDKorIWFBDbQlvpQps0aZtumSTNMpM9eX5/zKSEkKRJO2fOLM/79ZpXZ845c87z7bTzzHc536+oKsYYY5KXx+0AjDHGuMsSgTHGJDlLBMYYk+QsERhjTJKzRGCMMUnOEoExxiS5FLcDGK2ioiItKytzOwxjjIkrGzZsqFPV4sH2xV0iKCsrY/369W6HYYwxcUVEdg+1z5qGjDEmyVkiMMaYJGeJwBhjklzc9REMpquri5qaGtrb290OxTXp6emUlpbi8/ncDsUYE2cSIhHU1NSQk5NDWVkZIuJ2OFGnqtTX11NTU8OUKVPcDscYE2cSommovb2dwsLCpEwCACJCYWFhUteIjDHHLiESAZC0SaBPspffmESmqjy2aT/7G9scOX/CJIJ4kp2dDUB1dTVz5851ORpjTKw70NTOtQ++znPbax05f3Imgs2r4I65cGt+6M/Nq9yOyBhjhrSzNgDA9LHZjpw/+RLB5lXw2HXQtBfQ0J+PXXdcyeDLX/4yP/nJT468vvXWW/nGN77BihUrWLx4MfPmzePRRx8d9hw9PT3ceOONnHTSScyfP5977rkHgI9//OP86U9/OnLcFVdccdRzGWMSy85DLQBMH5fjyPmTLxE8ezt0DWhn62oLbT9Gl19+OatWvZVIVq1axSc+8Qn++Mc/8tprr/H888/zxS9+keGWBb3vvvvIy8tj3bp1rFu3jp///OdUVVVx5ZVXcv/99wPQ1NTEyy+/zPnnn3/MsRpj4k9FbYDCrFQKslIdOX9CDB8dlaaa0W0fgUWLFlFbW8v+/fvx+/2MGTOG8ePHc8MNN/Diiy/i8XjYt28fhw4dYvz48YOe46mnnmLz5s2sXr06FE5TEzt37uT9738/11xzDX6/n4cffphLLrmElJTk+9iMSWY7awOUO9QsBMmYCPJKw81Cg2w/DpdeeimrV6/m4MGDXH755fz2t7/F7/ezYcMGfD4fZWVlww7vVFV+9KMfcfbZZ79j38c//nEeeOABHnroIX71q18dV5zGmPiiquw41MJFCyc6do3kaxpacQv4Mt6+zZcR2n4cLr/8ch566CFWr17NpZdeSlNTE2PHjsXn8/H888+ze/eQE/8BcPbZZ/PTn/6Urq4uAHbs2EEwGATgk5/8JHfeeScAs2fPPq44jTHxpbalg5b2bqaPdaZ/AJKxRjD/stCfz94eag7KKw0lgb7tx2jOnDm0tLRQUlLChAkTuOKKK7jwwguZN28eS5YsYebMmcO+/9Of/jTV1dUsXrwYVaW4uPhIJ/G4ceOYNWsWF1988XHFaIyJPzsPOTtiCECG68CMRUuWLNGB6xFs27aNWbNmuRSR81pbW5k3bx6vvfYaeXl5Qx6X6H8PxiSjX71UxW2PbWXdf59JcU7aMZ9HRDao6pLB9iVf01CceeaZZ5g1axbXXnvtsEnAGJOYdtYGyM/0UZTtzIghSMamoThz5plnHrV/wRiTuCoOBZg+NtvRaWSsRmCMMTFKVdlR20K5gx3FkECJIN76OiIt2ctvTCKqD3bS2NrlaEcxJEgiSE9Pp76+Pmm/DPvWI0hPT3c7FGNMBB0ZMTTO2UTgeB+BiHiB9cA+Vb1giGMuAVYDJ6nq+sGOGU5paSk1NTX4/f7jCzaO9a1QZoxJHBW14TmGHG4aikZn8eeBbUDuYDtFJCd8zKvHegGfz2crcxljEs7O2gA5aSmMyz32YaMj4WjTkIiUAucDvxjmsK8D3wVseS1jjOln56EA5eOcHTEEztcI7gRuAgat14jIYmCSqj4hIjc6HMuwenqV1/ccpqouSEOwk4ZgJ/XBTnp7ldIxGZSOyTzy57i8NNJSvG97v6rSEOxkV10Qj8CJ43PJTrPRucaYY7ezNsDymcWOX8exbyoRuQCoVdUNIrJskP0e4IfAJ0dwrquAqwAmT54csRi7enp5pbKev2w5yFNbDlIX6DyyLy3FQ2FWKiLCnza20TugH7ooO40JeemMy02jLtDJLn+A5vbutx1TVpjJrAm5lI/NJjsthcxULxmpKaSmePC3dFBzuJV9h9vY19hGZqqXKUVZTC3OZkpRFvkZPg63dlIX6KQ+0ElrVzflxdnMLcmjfGw2Pm/kKnO9vUpjWxcNwQ46unspyEplTGYq6b53JrvOnl56ewc/T9+PlhSPkBLB+IxJRoeDndQFOhzvHwBnawRLgZUich6QDuSKyAOq+tHw/hxgLrA2XO0ZD6wRkZUDO4xV9V7gXghNMXEswVT6A6yramBPQyu7G1rZU99KVV2QQEc3male3jdzLOfOHc/8knwKslPJSvUeqY519fRysKmdvYdbqTncxoHGdg42t7G/sZ2aw20UZKVy4YKJTC3OZmpxFt09yrYDzUcef9lykMEGNGWmeikdk0FJfgbBzh6e2+5n1frBp8P2eYWuntBJ0lI8zByfw5isVHxeD6leDylewesRBMEjoS9lQVA0/HcIvQqtnd0EOrppbu8m0N5FY2sXh1s735HoALJSveRm+Ojs7qWtq4f2rp5BjxtMZqqXvAwfuek+8jJ8ZKenkJ2W8taf/R/pKWSkevF5PPi8oSTi8woeEUTAI6HnmalecsLvt0RjEl2FPzRiqNzhEUPgYCJQ1ZuBmwHCNYIv9UsCqGoTUNT3WkTWho8Z9aihkXh22yG+9eR2UjxC6ZgMJhdmsWhyPqdPL+b06UXv+PXbn8/rYVJBJpMKMkd8vbNmjzvyvLdXaevqobWzh9bObjq6eynOTiM/0/eOtr/m9i6q/EGa2roozE6lKDuNMZmpeD1CVV2QLfubeGNfE1sPNNMQ7KSzu5fuXqWrp5funr4vfaVXQVEEOfJLve/LNDs9hdz0FCbmpZOfmUphViqF2aFFL9JSPDQEQ8mhIdhJc1sXqSkeMnxeMlK9pPu8eD3vbK/U8PVUobtHaW7vormti+b2LprauqhtaaeqroeW9m5a2rvo6B6iWjFCmalexmSG4g7Fn0ZxThrjctIYn5fOuNx0xuelMzYnfdB4jYl1O/pWJXP4HgJwYYoJEbkdWK+qa6J53UsWl3Lu3AlMyEuP+q9Jj0fISkshKy0FGL73Pzfdx4JJ+YPuKx+bTfnYbC5aWBL5IKOsu6eXYEcPLR1dBDq6ae3sobtH6e7ppTOc1Ho1VJ9RVXp6IdjZfSSRtLR309jaRV2gA3+gg20HWqgLdNA9oMqS6vWEE38mkwsyGZuTRm64ppKTnkJBViqTCzIpCDcDGhMrdh4KkJXqpSQ/4+gHH6eoJAJVXQusDT8fdOJ/VV3mZAyF2WkUOnkBMyopXg95mR7yMn0RO2dvr9LQ2snBpnYONbdzoCnUdLenIciehlY27D5My4B+nD5ZqV4mFYSSxZSirLcexVkUZ6dZkjBRVxFelSwa//ZsWItJGB6PUJSdRlF2GnNLBp+ptbO7l5b2Lprbu2lu68Lf0sHew63srm9lb0Mru+qCrH3TT2fPW01XxTlpLJqUz6LJY1g0OZ8FpflkpA7dlGhMJOysbeG0cudHDIElApNkUlM8odph9tBNdD29yv7GNnbVBdnlD/CvmiZe39vIU1sPAaHO+tOnF3HmrHEsnzWWsTk2tYeJrKa2Lg41dzg+tUQfSwTGDOD1yJHBAWfMeOsXWUOwk417D/Pijjqe3nqIZ7bVArBocj4rF0zk/PkTLCmYiKiodX5Vsv4sERgzQgVZqSyfOY7lM8fxtQtns/1gC89sPcSTbxzktse28vXHt7K0vIgLF0zk1KmFlI7JsL4Fc0z65hgqt0RgTOwSEWZNyGXWhFyuXTGdnYdaWLNpP2s27eem1ZsByM/0Ma8kj7kleZxWXsQpUwttKKsZkYraAKkpHkrHjHzI+vFIiDWLjYkVqsrWA81s3NvIG/ua2FzTxJsHW+juVcbmpHHhgolctHAi80ryrLZghvT/7l/H/sY2/nL9eyN2zuHWLLYagTERJCLMmZjHnIlvjVpq7+rhue21/On1ffzmld3c9/cqZo7P4Svnz+a0tufg2duhqQbySmHFLTD/MhdLYGJBRW2A+aXRW6PcEoExDkv3eTlv3gTOmzeBptYunnzjAD9dW8mqX/2Ak9LuI007Qgc27YXHrgs9t2SQtNq7eth7uJUPLIrejaM2YYsxUZSX6eMjJ0/mqRveyzdyHnkrCfTpagvVEEzSqqoLohq9jmKwGoExrkj3eUnvODT4zqbBJx40yaFv6Oi04uglAqsRGOOWvMGXFu3Njf+5pMyxq/QHEIGpxVlRu6YlAmPcsuIW8L19QrFWTeUO/QiNrZ1DvMkkuoraAJPGZA47I3KkWSIwxi3zL4ML74K8SYBA3iR2vvtb3NPwLj587z/wt3Qc9RQm8VT6g0yLYm0ArI/AGHfNv+xtI4QWAPdN93PV/23gs79Zz0NXnUpqiv1eSxY9vcouf4DTyqM7V7L9CzMmxpw+vZjvXTqf1/Y08q0nt7kdjomi/Y1tdHT3RrWjGCwRGBOTLpg/kU8tLeP+l6t5bNN+t8MxUdI3YiiaQ0fBEoExMevmc2fxrhPG8J8Pbz4yCZlJbG4MHQVLBMbErNQUDz/5t8Vk+Lxc/cBrBDsGX13NJI5Kf4DCrFTGZKVG9bqWCIyJYePz0rnrI4vY5Q9w08ObibdJIs3oVNQGmBblZiGwRGBMzFtaXsSNZ8/kic0HuHttpdvhGAdV+gNRbxYCSwTGxIWrz5jKygUT+d5f3+SpLQfdDsc4oD7QweHWrqh3FIMlAmPigojwPx+az/zSPG74/UbePGidx4mm0h8EiPrNZGCJwJi4ke7zcu/HlpCZlsKn/28dh4M2DUUicWvoKFgiMCaujM9L596PvYtDzR1c89vX6O7pdTskEyEVtQEyfF4m5mUc/eAIs0RgTJxZNHkM3/rAPF7ZVc89L+5yOxwTIZX+AFOLs/C4sK61JQJj4tAli0s4f94E7nxmB9sONLsdjomAitqAK81CYInAmLgkInz94rnkZaTyhVWb6Oy2JqJ41tbZw77GNleGjkIUEoGIeEXkdRF5fJB9XxCRrSKyWUSeFZETnI7HmERRkJXKtz84j20Hmrnr2Z1uh2OOQ6XfvY5iiE6N4PPAUFMovg4sUdX5wGrgf6IQjzEJ46zZ4/jQu0q5e20Fr+857HY45hj1JYKErBGISClwPvCLwfar6vOq2hp++Q9g8LX7jDFDuuXC2YzPTeeLf9hEe1eP2+GYY1BZG8AjUFaU6cr1na4R3AncBIykAfNK4M+ORmNMAspN9/E/H1rALn+Q7/x5u9vhmGNQ4Q8wuSCTtJToLU/Zn2OJQEQuAGpVdcMIjv0osAT43hD7rxKR9SKy3u/3RzhSY+LfadOL+OR7QusXvLDD/o/Em8raoGv9A+BsjWApsFJEqoGHgOUi8sDAg0TkTOC/gZWqOugirap6r6ouUdUlxcXFDoZsTPz68rkzmTEumy/9YRMNdtdx3Oju6aWqLujKrKN9HEsEqnqzqpaqahnwYeA5Vf1o/2NEZBFwD6EkUOtULMYkg3SflzsvX0RTaxf/9ci/bMrqOLH3cBudPdFfnrK/qN9HICK3i8jK8MvvAdnAH0Rko4isiXY8xiSS2RNz+dLZM/jLloP8YUON2+GYEah0cY6hPinRuIiqrgXWhp/f0m/7mdG4vjHJ5NOnTeX57X5uW7OFd08p4ITC6M9maUauwuWho2B3FhuTcDwe4QeXLcDrEa5+4DVa2rvcDskMo6I2QHFOGnkZPtdisERgTAKamJ/Bj/9tMTsPtXD1AxtsCooYVukPUO5ibQAsERiTsN47o5jvXjKflyrquWn1Jnp7rfM41qiqq5PN9YlKH4Exxh2XvKuUg83tfO+vbzIuL52bz53ldkimH39LBy3t3a6sStafJQJjEtw1y6ZxsKmde17YxfjcdD61dIrbIZmwiiOTzeW4Goc1DRmT4ESEW1fO4azZ4/jGE9vYXNPodkgmLBaGjoIlAmOSgtcjfP9DCyjMSuWm1Zut8zhGVNQGyE5LYVxumqtxWCIwJknkZfr45gfmsf1gC3evrXA7HANU+oNMK85CJPrLU/ZnicCYJHLW7HFctHAiP36uwpa4jAEVtQFXbyTrY4nAmCTztQvnkJ/p46bVm+nusSYit7S0d3Gwud3Vyeb6WCIwJskUZKVy+0Vz+de+Ju792y63w0lalf4g4H5HMVgiMCYpnTdvAufOHc+dT++kui7odjhJqW/EkDUNGWNcc9vKOYjAT9dWuh1KUqrwB0jxCCcUurM8ZX+WCIxJUmNz07n8pEk88noN+xvb3A4n6VTUBigrysLndf9r2P0IjDGu+ewZ01CFe1+0voJoq/QHXJ9aoo8lAmOSWEl+Bh9YVMKD/9yDv2XQlWKNAzq7e9ld3xoTHcVgicCYpHf1sml09vTyy5eq3A4laeyuD9LTq5YIjDGxYVpxNufNm8BvXtlNU6stYhMNlTGwKll/lgiMMfz7snICHd38+pVqt0NJChUxNHQULBEYYwgter9i5lh++VIVwY5ut8NJeBW1ASbkpZOVFhsrAVgiMMYA8O/Ly2ls7eLBf+5xO5SEV+kPxkz/AFgiMMaELZ48hlOmFvCLv1XZNNUOUtXw0FFLBMaYGHT1GdM42NzOoxv3uR1KwjrQ1E5rZ09MTDbXxxKBMeaIM2YUM3N8Dve+uMsWu3fIWyOGYuNmMrBEYIzpR0S4+oxp7KwN8Nz2WrfDSUixsjxlf5YIjDFvc/78CZTkZ/CzF2wyOidU+oPkpKdQnO3u8pT9WSIwxryNz+vh06dPYf3uw6yvbnA7nITT11Hs9vKU/TmeCETEKyKvi8jjg+xLE5Hfi0iFiLwqImVOx2OMObrLT5pEfqaPn71gk9FFWqyNGILo1Ag+D2wbYt+VwGFVLQfuAL4bhXiMMUeRmZrCx08t45lth6iobXE7nITR3N7FoeYOpo2NnY5icDgRiEgpcD7wiyEOuQj4dfj5amCFxFJ9yZgk9olTTyDd57FaQQTt6lueMslqBHcCNwFD3Z1SAuwFUNVuoAkodDgmY8wIFGan8ZGTJ/PIazXsPGS1gkg4sjxlDI0YAgcTgYhcANSq6oYInOsqEVkvIuv9fn8EojPGjMS1y6eTlZrCd/683e1QEkJleHnKyQXuL0/Z34gSgYhkiYgn/HyGiKwUEd9R3rYUWCki1cBDwHIReWDAMfuASeHzpgB5QP3AE6nqvaq6RFWXFBcXjyRkY0wEFGSlcs37ynl2ey0vV9a5HU7cq/QHOKEwMyaWp+xvpNG8CKSLSAnwFPAx4P7h3qCqN6tqqaqWAR8GnlPVjw44bA3wifDzD4WPsdsZjYkhn1paxsS8dL715Da72/g4VfqDMTdiCEaeCERVW4EPAner6qXAnGO5oIjcLiIrwy/vAwpFpAL4AvDlYzmnMcY56T4vN55zIm/sa+bRTTYH0bHq6umlui62Zh3tM+JEICKnAlcAT4S3eUd6EVVdq6oXhJ/foqprws/bVfVSVS1X1ZNV1YYnGBODLlpQwtySXL7/1x20d/W4HU5c2tPQSnevxnWN4HrgZuCPqrpFRKYCzzsWlTEmpng8wn+dO4t9jW3c/3K12+HEpVgdMQQjTASq+oKqrlTV74Y7jetU9TqHYzPGxJD3lBexfOZYfvJcBY2tnW6HE3cqw/cQTI2hWUf7jHTU0O9EJFdEsoA3gK0icqOzoRljYs1N55xIS0c3v3yp2u1Q4k6lP8DYnDRy04824DL6Rto0NFtVm4GLgT8DUwiNHDLGJJGZ43M5e844fvVSFc3tXW6HE1cq/YGY7CiGkScCX/i+gYuBNaraBdg4MmOS0LXLp9PS3s2vrVYwYqpKZW3sTTbXZ6SJ4B6gGsgCXhSRE4Bmp4IyxsSuuSV5rJg5lvteqiLQ0e12OHHBH+igub07plYl62+kncV3qWqJqp6nIbuB9zkcmzEmRl27YjqNrV385pXdbocSFyprQx3FsThiCEbeWZwnIj/sm+9HRH5AqHZgjElCCyfl894Zxfzib7to7bRawdG8tU5xHCcC4JdAC3BZ+NEM/MqpoIwxse+65eXUBzv53at73A4l5lX6A2SmepmQl+52KIMaaSKYpqpfU9Vd4cdtwFQnAzPGxLYlZQWcOrWQe17cZXcbH0XfHEOxutzKSBNBm4ic1vdCRJYCbc6EZIyJF9etmI6/pYOH/mm1guGERgzFbmv6SBPB1cBPRKQ6PK30j4HPOhaVMSYunDK1gJOnFHD32kqrFQyhrbOHfY1tMds/ACMfNbRJVRcA84H5qroIWO5oZMaYmCci3HDmDGpbOnjQagWDOtJRHKMjhmCUK5SpanP4DmMITRttjElyp04r5N1WKxhSVV3szjHU53iWyYnNXg9jTNRdf+YM/C0dNoJoEH2JoKwwMROBTTFhjAFCtYJTphbw0xesVjBQdV2QiXnppPtGvIRL1A2bCESkRUSaB3m0ABOjFKMxJg701Qp+a7WCt9lVF2RKDDcLwVESgarmqGruII8cVU2JVpDGmNh3ytRCTp1ayE/XVtLWabWCPtX1wZhuFoLjaxoyxpi3uf7M6dQFOnjgHzYHEcDhYCeNrV1MKbJEYIxJEu+eWsjp04u467md1AU63A7HdVX1oY5iSwTGmKTytQvn0N7Vw7ef3O52KK6r8lsiMMYkofKx2Xz69Kk8/FoN66ob3A7HVdX1QbweYVJBptuhDMsSgTEm4q5dXs7EvHS++qc36O7pdTsc1+yqCzJpTAY+b2x/1cZ2dMaYuJSZmsItF85m+8EWfp3Ei9dU1wUpi/FmIbBEYIxxyNlzxnPGjGLueHoHh5rb3Q4n6lSVqrpgzPcPgCUCY4xDRITbVs6hs7uXbz6xze1wos7f0kFrZ48lAmNMcisryuLqM6ayZtN+/rGr3u1wompXXXyMGAIHE4GIpIvIP0Vkk4hsEZHbBjlmsog8LyKvi8hmETnPqXiMMe743LJySvIzuHXNlqTqOK6Og8nm+jhZI+gAlofXMVgInCMipww45ivAqvD6Bh8G7nYwHmOMCzJSvXz1gllsP9iSVPMQVdUFSU3xMDE/w+1QjsqxRKAhgfBLX/gxcMZSBXLDz/OA/U7FY4xxz9lzxnNaeRE/eOpN6pPkjuOquiAnFGTi9cT+jP2O9hGIiFdENgK1wNOq+uqAQ24FPioiNcCTwLVOxmOMcYeIcOvK2bR29vC9v77pdjhRES8jhsDhRKCqPaq6ECgFThaRuQMO+Qhwv6qWAucBvxGRd8QkIleJyHoRWe/3+50M2RjjkPKxOXxqaRm/X7+XzTWNbofjqJ5eZXd9qyWC/lS1EXgeOGfAriuBVeFjXgHSgaJB3n+vqi5R1SXFxcUOR2uMccp1K6ZTlJ3GLY9uobc3cde22t/YRmdPryUCESkWkfzw8wzgLGDgLFR7gBXhY2YRSgT2k9+YBJWT7uPL58xk495GnvjXAbfDccyR5SmTPREAE4DnRWQzsI5QH8HjInK7iKwMH/NF4DMisgl4EPikqibuzwRjDB9YVML0sdn86LmdCVsrqA5PPz01ThKBY6uMqepmYNEg22/p93wrsNSpGIwxscfjEa5dMZ3rHnydv2w5yHnzJrgdUsTt8gfJSvVSnJPmdigjYncWG2Oi7vx5E5hanMVdzyZmraC6PjTZnEjsDx0FSwTGGBd4PcK1y8vZfrCFp7cdcjuciKuKk1lH+1giMMa44sL5EykrzOSuZ3eSSF2Dnd291Bxui5v+AbBEYIxxSYrXw7+/r5wt+5t5dlut2+FEzN7DrfT0alzMMdTHEoExxjUXLyphckEmdz2XOLWCvsnmphRbIjDGmKPyeT38+/umsbmmibVvJsYtRH33EEyxGoExxozMBxaVUpKfwR3P7EiIWkFVXZC8DB9jslLdDmXELBEYY1yVmuLhuhXlbK5pSoi+gr6ho/HEEoExxnUfXFzKCYWZ/PDpHXF/X0F1XStTCjPdDmNULBEYY1zn83q4bvl0th5o5qmtB90O55i1d/Wwv6nNagTGGHMsLlo4kanFWdzxdPzebbynoRXV+FinuD9LBMaYmJDi9fD5FdN581BL3M5MWhVH6xT3Z4nAGBMzLpw/kRnjsrnzmR30xGGtoDrOpp/uY4nAGBMzPB7hhjNnUOkPsmbTPrfDGbXq+iAFWankZfjcDmVULBEYY2LK2XPGM2tCLnc8vZPO7l63wxmVqrogZXE2YggsERhjYozHI9x09onsaWjloXV73A5nVKrrWuOuWQgsERhjYtCyE4t595QC7np2J8GObrfDGZG2zh4ONrfH1dQSfSwRGGNijojwn+fOpC7QyS/+VuV2OCPStzyl1QiMMSZCFk8ewzlzxnPvi5XUBTrcDueojsw6aonAGGMi50tnn0hbVw8/fq7C7VCOqspqBMYYE3nlY7O5/KRJ/PbV3eypb3U7nGFV1wUpyk4jOy3F7VBGzRKBMSamfX7FDLwe4QdPv+l2KMOqrmtlSlH8DR0FSwTGmBg3Pi+dTy2dwqMb97O5ptHtcIZUVR+Mu6kl+lgiMMbEvM8tm0ZRdiq3PbY1JhevCXR042/piMv+AbBEYIyJA7npPm48+0Q27D7Mmk373Q7nHeJ5xBBYIjDGxIkPvWsScybm8p0/b6e1M7ZuMjtyD4E1DRljjHO8HuFrF87hQFM797ywy+1w3uatWUets/htRCRdRP4pIptEZIuI3DbEcZeJyNbwMb9zKh5jTPw7eUoBF8yfwM9eqGRfY5vb4RxRVdfKuNw0MlPjb+goOFsj6ACWq+oCYCFwjoic0v8AEZkO3AwsVdU5wPUOxmOMSQA3nzcLgG8/uc3lSN5SHccjhsDBRKAhgfBLX/gxsLv/M8BPVPVw+D21TsVjjEkMJfkZfPaMaTy++QCv7qp3Oxwg1DQUrx3F4HAfgYh4RWQjUAs8raqvDjhkBjBDRF4SkX+IyDlOxmOMSQxXnzGVkvwMvvroG3T1uLtmQXN7F/XBTk6wGsHgVLVHVRcCpcDJIjJ3wCEpwHRgGfAR4Ocikj/wPCJylYisF5H1fr/fyZCNMXEgMzWFW1fOYcehAPf93d3ZSd8aOhqfHcUQpVFDqtoIPA8M/MVfA6xR1S5VrQJ2EEoMA99/r6ouUdUlxcXFjsdrjIl9Z80ex1mzx3HnMzvY2+DePERVcbpOcX9Ojhoq7vt1LyIZwFnA9gGH/YlQbQARKSLUVBRb48KMMTHr1pVzEITbHtviWgzVdaEkdEKBJYLBTACeF5HNwDpCfQSPi8jtIrIyfMxfgXoR2UqoxnCjqsZG748xJuaV5Gdww1nTeWZbLU9tOehKDNX1QSbkpZOR6nXl+pHg2KBXVd0MLBpk+y39nivwhfDDGGNG7VNLp/Dwhn3cumYLS8uLyIryNNChBevjtzYAdmexMSbO+bwevvmBuexvaueHT++I+vX3NLTG7R3FfSwRGGPi3pKyAq5492R++VIV/6xqiNp1m9u7aAh2MjmO+wfAEoExJkH813mzmDQmky/9YRPBjuhMSte3atoJhVYjMMYY12WlpfD9Sxew93Ar34zS9BO7w4lgcoElAmOMiQknTyngM6dP5Xev7mHtm87PWLO7IXQPgdUIjDEmhnzhrBnMGJfNfz68mabWLkevtae+lcKsVHLSfY5ex2mWCIwxCSXd5+WHly2kPtDJVx99w9Fr7a5vZXKc1wbAEoExJgHNLcnj8yums2bTfh7duM+x6+xpaOWEOO8fAEsExpgE9bll01hywhi+8sc3HJmLqKO7h/1NbUyO85vJwBKBMSZBpXg93HH5QgBu+P1GuiM8XXXN4TZUsRqBMcbEskkFmXz94rms332Yn66tjOi5E+UeArBEYIxJcBcvKuHihRO589mdvLbncMTOu7s+NHTUOouNMSYO3H7xXMbnpnP9QxtpaY/MkNLdDa1kpnopzk6LyPncZInAGJPwctN9/O+HF7KvsY0vP/IvQhMfH5899a1MLshERCIQobssERhjksKSsgK+9P4TeWLzAf7vld3Hfb7dDa1xP7VEH0sExpik8dn3TmXFzLF844mtvH4c/QW9vRqefjr+h46CJQJjTBLxeIQfXLaAsTnp/MfvXudwsPOYznOopZ3O7l6rERhjTDzKz0zl7isW42/p4AurNtLbO/r+gt0JNHQULBEYY5LQgkn5fOWCWTz/pp//ePA1quuCo3r/kXsI4nxBmj7RXdzTGGNixMdOOYGGYCc/e6GSp7Yc4tIlk7huRTkT8jKO+t7dDUFSPMLE/PQoROo8qxEYY5KSiHD9mTN48cb3ccW7J7N6w17O+N5a7l5bcdT37q5vpWRMBinexPgKTYxSGGPMMRqbm85tF83luS8u430nFvM/f3mTF3f4h33P7vrEGToKlgiMMQYIzUv0vx9eRPnY8KI2bUPfgby7PpgwHcVgicAYY44ILWqzgNqWDm5/bOugxzS2dtLc3p0wHcVgicAYY95mfmk+1yybxsOv1fD01kPv2H9kwXqrERhjTOK6dvl0Zk3I5eZH/vWOm852NyTWPQRgicAYY94hNcXDDy5dQFNbJ1959I23TVK3p2/6aessPjoRSReRf4rIJhHZIiK3DXPsJSKiIrLEqXiMMWY0Zk/M5fozZ/DE5gPc/Mi/6AqvcLa7vpWxOWlkpibObVhOlqQDWK6qARHxAX8XkT+r6j/6HyQiOcDngVcdjMUYY0btmmXTaO/q4UfPVbCvsY27r1jM7obWhGoWAgdrBBoSCL/0hR+DTerxdeC7QLtTsRhjzLEQEb74/hP57iXzeLmynkt/9goVtQEmJ9CIIXC4j0BEvCKyEagFnlbVVwfsXwxMUtUnnIzDGGOOx+UnTeb+T53EvsNtNAQ7rUYwGqrao6oLgVLgZBGZ27dPRDzAD4EvHu08InKViKwXkfV+//B3/BljjBNOn17MHz53KqeVF7Fi1li3w4koicSSbSO6kMgtQKuqfj/8Og+oBPqaj8YDDcBKVV0/1HmWLFmi69cPudsYY8wgRGSDqg46IMfJUUPFIpIffp4BnAVs79uvqk2qWqSqZapaBvyDoyQBY4wxkedk09AE4HkR2QysI9RH8LiI3C4iKx28rjHGmFFwbPioqm4GFg2y/ZYhjl/mVCzGGGOGZncWG2NMkrNEYIwxSc4SgTHGJDlLBMYYk+QsERhjTJKL2g1lkSIifmB3+GUe0NRvd//Xgz0vAuqO4/IDr3csxw22b7hyDHzd97z/tmiU63jLNHCb05/VUDGM5phIfVb9n7v9b3Ck20fzWUF8/BtMhs9q4Ov+3xf5qlo8aASqGrcP4N6hXg/2HFgfyesdy3GD7RuuHMOUpf82x8t1vGUayecTyc8qWuWKt3+DI90+ms8qWuWyz+rYy3W0WOK9aeixYV4P9TyS1zuW4wbbN1w5Br5+bIhjjsdIznW8ZRq4zenPaqTnisZnNdJYRuJ4/w2OdLt9VsfPjc9q4OsRfV/EXdPQ8RCR9TrEXBvxLBHLlYhlAitXPEnEMg0l3msEo3Wv2wE4JBHLlYhlAitXPEnEMg0qqWoExhhj3inZagTGGGMGsERgjDFJzhKBMcYkOUsEYSJyuoj8TER+ISIvux1PpIiIR0S+KSI/EpFPuB1PJIjIMhH5W/jzWuZ2PJEkIlnhZVkvcDuWSBCRWeHPabWIfM7teCJFRC4WkZ+LyO9F5P1ux3O8EiIRiMgvRaRWRN4YsP0cEXlTRCpE5MvDnUNV/6aqVwOPA792Mt6RikS5gIsIrRndBdQ4FetIRahMSmiJ03RioEwQsXIB/CewypkoRydC/6+2hf9fXQYsdTLekYpQuf6kqp8BrgYudzLeaEiIUUMi8l5CXwz/p6pzw9u8wA5CS2TWEFol7SOAF/j2gFP8P1WtDb9vFXClqrZEKfwhRaJc4cdhVb1HRFar6oeiFf9gIlSmOlXtFZFxwA9V9YpoxT+UCJVrAVBIKMHVqerj0Yl+cJH6fxVekfBzwG9U9XfRin8oEf6++AHwW1V9LUrhO8KxFcqiSVVfFJGyAZtPBipUdReAiDwEXKSq3wYGrXaLyGSgKRaSAESmXCJSA3SGX/Y4GO6IROqzCjsMpDkS6ChF6LNaBmQBs4E2EXlSVXudjHs4kfqsVHUNsEZEngBcTwQR+qwE+A7w53hPApAgiWAIJcDefq9rgHcf5T1XAr9yLKLIGG25HgF+JCKnAy86GdhxGFWZROSDwNlAPvBjRyM7PqMql6r+N4CIfJJwrcfR6I7NaD+rZcAHCSXsJ50M7DiN9v/VtcCZQJ6IlKvqz5wMzmmJnAhGTVW/5nYMkaaqrYQSXMJQ1UcIJbiEpKr3ux1DpKjqWmCty2FEnKreBdzldhyRkhCdxUPYB0zq97o0vC3eJWK5ErFMkJjlSsQyQeKWa0QSORGsA6aLyBQRSQU+DKxxOaZISMRyJWKZIDHLlYhlgsQt14gkRCIQkQeBV4ATRaRGRK5U1W7gP4C/AtuAVaq6xc04RysRy5WIZYLELFcilgkSt1zHIyGGjxpjjDl2CVEjMMYYc+wsERhjTJKzRGCMMUnOEoExxiQ5SwTGGJPkLBEYY0ySs0RgEoKIBKJ8vYisWSGhtRWaRGSjiGwXke+P4D0Xi8jsSFzfGLBEYMygRGTYebhU9T0RvNzfVHUhsAi4QESONm//xYRmKDUmIiwRmIQlItNE5C8iskFCK5rNDG+/UEReFZHXReSZ8LoGiMitIvIbEXkJ+E349S9FZK2I7BKR6/qdOxD+c1l4/+rwL/rfhqcoRkTOC2/bICJ3iciw6wuoahuwkdBMmIjIZ0RknYhsEpGHRSRTRN4DrAS+F65FTBuqnMaMlCUCk8juBa5V1XcBXwLuDm//O3CKqi4CHgJu6vee2cCZqvqR8OuZhKa8Phn4moj4BrnOIuD68HunAktFJB24Bzg3fP3iowUrImOA6bw1XfgjqnqSqi4gNO3Blar6MqE5cG5U1YWqWjlMOY0ZEZuG2iQkEckG3gP8IfwDHd5axKYU+L2ITABSgap+b10T/mXe5wlV7QA6RKQWGMc7l8f8p6rWhK+7ESgjtALWLlXtO/eDwFVDhHu6iGwilATuVNWD4e1zReQbhNZdyCY0D85oymnMiFgiMInKAzSG294H+hGhJS7XhBdOubXfvuCAYzv6Pe9h8P8zIzlmOH9T1QtEZArwDxFZpaobgfuBi1V1U3ixmmWDvHe4chozItY0ZBKSqjYDVSJyKYSWFhSRBeHdebw11/wnHArhTWBqvyURj7rAebj28B1CC9gD5AAHws1R/ddlbgnvO1o5jRkRSwQmUWSGpxTue3yB0JfnleFmly3AReFjbyXUlLIBqHMimHDz0jXAX8LXaQGaRvDWnwHvDSeQrwKvAi8B2/sd8xBwY7izexpDl9OYEbFpqI1xiIhkq2ogPIroJ8BOVb3D7biMGchqBMY45zPhzuMthJqj7nE3HGMGZzUCY4xJclYjMMaYJGeJwBhjkpwlAmOMSXKWCIwxJslZIjDGmCRnicAYY5Lc/wcltH5GZboWDgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "dls = dsets.dataloaders(bs=16, after_item=[albu_aug, ToTensor],\n", + " after_batch=[IntToFloatTensor(div_mask=255), \n", + " Normalize.from_stats(*imagenet_stats)])\n", + "img_size = [round(0.9*320) for _ in range(2)]\n", + "\n", + "encoder = timm.create_model('efficientnet_b0', pretrained=True)\n", + "\n", + "# Let's use self attentions and Mish activation function \n", + "model = DynamicTimmUnet(encoder, 3, img_size, self_attention=True, act_cls=Mish)\n", + "\n", + "# We'll also use ranger optimizer with is RAdam with Lookahead\n", + "learn = Learner(dls, model, metrics=metrics, loss_func=loss_func, splitter=dynamic_unet_splitter, opt_func=ranger).to_fp16()\n", + "learn.freeze()\n", + "learn.lr_find()" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "<style>\n", + " /* Turns off some styling */\n", + " progress {\n", + " /* gets rid of default border in Firefox and Opera. */\n", + " border: none;\n", + " /* Needs to be in here for Safari polyfill so background images work as expected. */\n", + " background-size: auto;\n", + " }\n", + " .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n", + " background: #F44336;\n", + " }\n", + "</style>\n" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: left;\">\n", + " <th>epoch</th>\n", + " <th>train_loss</th>\n", + " <th>valid_loss</th>\n", + " <th>dice_coeff_adj</th>\n", + " <th>hd_dist_adj</th>\n", + " <th>custom_metric_adj</th>\n", + " <th>time</th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <td>0</td>\n", + " <td>0.991883</td>\n", + " <td>0.699887</td>\n", + " <td>0.562718</td>\n", + " <td>0.762979</td>\n", + " <td>0.682874</td>\n", + " <td>02:08</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Let's also use flat cosine annealing lr shceduler\n", + "lr = 1e-3\n", + "learn.fit_flat_cos(1, lr)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "<style>\n", + " /* Turns off some styling */\n", + " progress {\n", + " /* gets rid of default border in Firefox and Opera. */\n", + " border: none;\n", + " /* Needs to be in here for Safari polyfill so background images work as expected. */\n", + " background-size: auto;\n", + " }\n", + " .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n", + " background: #F44336;\n", + " }\n", + "</style>\n" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: left;\">\n", + " <th>epoch</th>\n", + " <th>train_loss</th>\n", + " <th>valid_loss</th>\n", + " <th>dice_coeff_adj</th>\n", + " <th>hd_dist_adj</th>\n", + " <th>custom_metric_adj</th>\n", + " <th>time</th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <td>0</td>\n", + " <td>0.805488</td>\n", + " <td>0.588406</td>\n", + " <td>0.645035</td>\n", + " <td>0.856163</td>\n", + " <td>0.771712</td>\n", + " <td>02:24</td>\n", + " </tr>\n", + " <tr>\n", + " <td>1</td>\n", + " <td>0.660503</td>\n", + " <td>0.477541</td>\n", + " <td>0.698127</td>\n", + " <td>0.851142</td>\n", + " <td>0.789936</td>\n", + " <td>02:23</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>" + ], + "text/plain": [ + "<IPython.core.display.HTML object>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Let's unfreeze the encoder layers and train with discriminative learning rates.\n", + "learn.unfreeze()\n", + "learn.fit_flat_cos(2, slice(lr/400, lr/4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Inference" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "def create_df(df, fnames):\n", + " df = df.copy()\n", + " df = df.pivot(index='id', columns='class', values='segmentation').reset_index()\n", + " \n", + " df['partial_fname'] = df.id\n", + " fname_df = pd.DataFrame({'partial_fname': [f'{fname.parts[-3]}_slice_{fname.parts[-1][6:10]}' for fname in fnames],\n", + " 'fname': fnames})\n", + "\n", + " df = df.merge(fname_df, on='partial_fname').drop('partial_fname', axis=1)\n", + "\n", + " df['case_id'] = df.id.apply(lambda x: x.split('_')[0])\n", + " df['day_num'] = df.id.apply(lambda x: x.split('_')[1])\n", + "\n", + " df['slice_w'] = df[\"fname\"].apply(lambda x: int(str(x)[:-4].rsplit(\"_\",4)[1]))\n", + " df['slice_h'] = df[\"fname\"].apply(lambda x: int(str(x)[:-4].rsplit(\"_\",4)[2]))\n", + " \n", + " channels = 3\n", + " stride = 2\n", + " for j, i in enumerate(range(-1*(channels-channels//2-1), channels//2+1)):\n", + " method = 'ffill'\n", + " if i <= 0: method = 'bfill'\n", + " df[f'fname_{j:02}'] = df.groupby(['case_id', 'day_num'])['fname'].shift(stride*-i).fillna(method=method)\n", + "\n", + " df['fnames'] = df[[f'fname_{j:02d}' for j in range(channels)]].values.tolist()\n", + " \n", + " return df\n", + "\n", + "def mask2rle(mask):\n", + " \"\"\"\n", + " img: numpy array, 1 - mask, 0 - background\n", + " Returns run length as string formated\n", + " \"\"\"\n", + " mask = np.array(mask)\n", + " pixels = mask.flatten()\n", + " pad = np.array([0])\n", + " pixels = np.concatenate([pad, pixels, pad])\n", + " runs = np.where(pixels[1:] != pixels[:-1])[0] + 1\n", + " runs[1::2] -= runs[::2]\n", + "\n", + " return \" \".join(str(x) for x in runs)\n", + "\n", + "def resize_img_to_org_size(img, org_size):\n", + " shape0 = np.array(img.shape[:2])\n", + " diff = org_size - shape0\n", + " if np.any(diff < 0):\n", + " img = pad_img_nc(img, (320, 384))\n", + " resized = unpad_img_nc(img, org_size)\n", + " else:\n", + " resized = pad_img_nc(img, org_size)\n", + " return resized\n", + "\n", + "def get_rle_masks(preds, df):\n", + " rle_masks = []\n", + " for pred, width, height in zip(preds, df['slice_w'], df['slice_h']):\n", + " upsized_mask = resize_img_to_org_size(pred, (height, width))\n", + " for i in range(3):\n", + " rle_mask = mask2rle(upsized_mask[:, :, i])\n", + " rle_masks.append(rle_mask)\n", + " return rle_masks\n", + "\n", + "def unpad_img_nc(img, org_size):\n", + " shape0 = np.array(org_size)\n", + " resize = np.array(img.shape[:2])\n", + " if np.any(shape0!=resize):\n", + " diff = resize - shape0\n", + " pad0 = diff[0]\n", + " pad1 = diff[1]\n", + " pady = [pad0//2, pad0//2 + pad0%2]\n", + " padx = [pad1//2, pad1//2 + pad1%2]\n", + " \n", + " if pady[0] != 0:\n", + " img = img[pady[0]:-pady[1], :, :]\n", + " \n", + " if padx[0] != 0:\n", + " img = img[:, padx[0]:-padx[1], :]\n", + " \n", + " img = img.reshape((*shape0, img.shape[-1]))\n", + " return img\n", + "\n", + "def pad_img_nc(img, up_size=None):\n", + " if up_size is None:\n", + " return img\n", + " shape0 = np.array(img.shape[:2])\n", + " resize = np.array(up_size)\n", + " if np.any(shape0!=resize):\n", + " diff = resize - shape0\n", + " pad0 = diff[0]\n", + " pad1 = diff[1]\n", + " pady = [pad0//2, pad0//2 + pad0%2]\n", + " padx = [pad1//2, pad1//2 + pad1%2]\n", + " padz = [0, 0]\n", + " img = np.pad(img, [pady, padx, padz])\n", + " img = img.reshape((*resize, img.shape[-1]))\n", + " return img\n", + "\n", + "def get_rle_masks(preds, df):\n", + " rle_masks = []\n", + " for pred, width, height in zip(preds, df['slice_w'], df['slice_h']):\n", + " upsized_mask = resize_img_to_org_size(pred, (height, width))\n", + " for i in range(3):\n", + " rle_mask = mask2rle(upsized_mask[:, :, i])\n", + " rle_masks.append(rle_mask)\n", + " return rle_masks" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "data_path = 'dataset/'\n", + "\n", + "train_path = Path(data_path+'train')\n", + "test_path = Path(data_path+'test')\n", + "\n", + "train_fnames = get_image_files(train_path)\n", + "test_fnames = get_image_files(test_path)\n", + "\n", + "sample_submission = pd.read_csv(data_path+'sample_submission.csv')\n", + "\n", + "if sample_submission.shape[0] > 0: \n", + " test = sample_submission.copy()\n", + "else:\n", + " test_fnames = train_fnames\n", + " test_path = train_path\n", + " train = pd.read_csv('dataset/train.csv', low_memory=False)\n", + " test = train.copy()\n", + " test = test.sample(frac=1.0, random_state=42)\n", + "\n", + "test_df = create_df(test, test_fnames)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "learn = load_learner('test_model.pkl')\n", + "# Sample the test set demonstartion purposes\n", + "test_df = test_df.sample(frac=0.1)\n", + "bs = learn.dls.bs\n", + "test_dl = learn.dls.test_dl(test_df, bs=bs, shuffle=False).to('cuda')" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████████████████████████████████████████████████████████████████████████████| 241/241 [00:50<00:00, 4.81it/s]\n" + ] + } + ], + "source": [ + "from tqdm import tqdm\n", + "import gc\n", + "\n", + "learn.model = learn.model.cuda()\n", + "learn.model.eval()\n", + "masks = []\n", + "\n", + "with torch.no_grad():\n", + " for i, b in enumerate(tqdm(test_dl)):\n", + " b.to('cuda')\n", + " b_preds = (sigmoid(learn.model(b)) > 0.5).permute(0, 2, 3, 1).cpu().detach().numpy().astype(np.uint8)\n", + "\n", + " masks.extend(get_rle_masks(b_preds, test_df.iloc[i*bs:i*bs+bs]))\n", + "\n", + " # test_preds[i*bs:i*bs+bs] = b_preds\n", + " del b_preds\n", + " torch.cuda.empty_cache()\n", + " gc.collect()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Submission" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "def get_case_id(fname):\n", + " return fname.parts[3] + '_' + fname.parts[5][:10]" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "from itertools import chain\n", + "\n", + "submission = pd.DataFrame({\n", + " 'id': chain.from_iterable([[get_case_id(fname)]*3 for fname in test_df['fname']]),\n", + " 'class': chain.from_iterable([['large_bowel', 'small_bowel', 'stomach'] for _ in test_df['fname']]),\n", + " 'predicted': masks,\n", + " })\n", + " \n", + "\n", + "# Merge with sample submission to preserve order to slices during scoring and avoid 0 scores\n", + "if sample_submission.shape[0] > 0:\n", + " del sample_submission['segmentation']\n", + " submission = sample_submission.merge(submission, on=['id', 'class'])\n", + "\n", + "submission.to_csv('submission.csv', index=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.11" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}