{"cells":[{"metadata":{},"cell_type":"markdown","source":"
\n

DISCLAIMER

\n

The approach that will be shown here is simply a fun experiment. It can provide idea or encourage ...

\n
\n\n\n---\n\n### Version Update 2\n\n[WIP]: Working to add **Model interpretability with Integrated Gradients**. FYI, the model `LamdaDenseNet` seems can able to take the leverage of TPU; the MXU seems remain to high. That's really interesting. \n\nReferences:\n\n- Integrated Gradients original [paper](https://arxiv.org/abs/1703.01365)\n- Original [implementation](https://github.com/ankurtaly/Integrated-Gradients)\n- Reusing code from [here](https://keras.io/examples/vision/integrated_gradients/) by [Aakash Kumar Nain](https://www.kaggle.com/aakashnain)\n\n\n### Version Update 1\n\nPreviously we've built a custom **mixed depth-wise group convolutional** layer in `tf.keras` and added it at the top of the pre trained model. In this version we will extend it and try something new. We will modify the pre-trained model **internally** and add some custom layer. Mainly, we will be tweaking with the following stuff. \n\n```\n- Dense Net 121\n- Mixed Depth-wise Group Convolution\n- Lamda-Layer: Long-Range Interaction without Attention.\n```\n\nThe end architecture will be look something like as follows, namely **LamdaDenseNet-121**: \n\n![LamdaDenseNet](https://user-images.githubusercontent.com/17668390/103923527-68532100-513f-11eb-97f5-d16806865f08.png)\n\n\nThe features flow will be something like this:\n\n- First we init the convolution part of `DenseNet-121` with a custom `MixDepthGroupConvolution2D` layer.\n- Second we iterate over the consecutive Dense Blocks and add a custom `LamdaLayer` in between the Dense Block and Transition Layer until Dense Block 5. "},{"metadata":{},"cell_type":"markdown","source":"## -)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"!pip install efficientnet -q\n!pip install einops -q\n\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom glob import glob\nimport albumentations as A \nfrom pylab import rcParams\nimport matplotlib.pyplot as plt\nimport plotly.graph_objects as go\nimport os, gc, cv2, random, warnings, math, sys, json, pprint\n\n# sklearn\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import roc_auc_score\n\n# tf \nimport tensorflow as tf\nimport efficientnet.tfkeras as efn\nfrom tensorflow.keras import backend as K\n\nos.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' \nwarnings.simplefilter('ignore')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SEED = 101\nTRAIN_DF = '../input/ranzcr-clip-catheter-line-classification/train.csv'\nTRAIN_IMG_PATH = '../input/ranzcr-clip-catheter-line-classification/train/'\nTEST_IMG_PATH = '../input/ranzcr-clip-catheter-line-classification/test/'\nCLASS_MAP = '../input/ranzcr-clip-catheter-line-classification/train_annotations.csv'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv(TRAIN_DF)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## RANZCR-CLiP Dataloader"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from kaggle_datasets import KaggleDatasets\n\ndef auto_select_accelerator():\n \"\"\"\n Reference: \n * https://www.kaggle.com/mgornergoogle/getting-started-with-100-flowers-on-tpu\n * https://www.kaggle.com/xhlulu/ranzcr-efficientnet-tpu-training\n \"\"\"\n try:\n tpu = tf.distribute.cluster_resolver.TPUClusterResolver()\n tf.config.experimental_connect_to_cluster(tpu)\n tf.tpu.experimental.initialize_tpu_system(tpu)\n strategy = tf.distribute.experimental.TPUStrategy(tpu)\n print(\"Running on TPU:\", tpu.master())\n except ValueError:\n strategy = tf.distribute.get_strategy()\n print(f\"Running on {strategy.num_replicas_in_sync} replicas\")\n return strategy\n\n\ndef build_decoder(with_labels=True, target_size=(256, 256), ext='jpg'):\n def decode(path):\n file_bytes = tf.io.read_file(path)\n if ext == 'png':\n img = tf.image.decode_png(file_bytes, channels=3)\n elif ext in ['jpg', 'jpeg']:\n img = tf.image.decode_jpeg(file_bytes, channels=3)\n else:\n raise ValueError(\"Image extension not supported\")\n\n img = tf.cast(img, tf.float32) / 255.0\n img = tf.image.resize(img, target_size)\n\n return img\n \n def decode_with_labels(path, label):\n return decode(path), label\n \n return decode_with_labels if with_labels else decode\n\n\ndef build_augmenter(with_labels=True):\n def augment(img):\n img = tf.image.random_flip_left_right(img)\n img = tf.image.random_flip_up_down(img)\n return img\n \n def augment_with_labels(img, label):\n return augment(img), label\n \n return augment_with_labels if with_labels else augment\n\n\ndef build_dataset(paths, labels=None, bsize=32, cache=True,\n decode_fn=None, augment_fn=None,\n augment=True, repeat=True, shuffle=1024, \n cache_dir=\"\"):\n if cache_dir != \"\" and cache is True:\n os.makedirs(cache_dir, exist_ok=True)\n \n if decode_fn is None:\n decode_fn = build_decoder(labels is not None)\n \n if augment_fn is None:\n augment_fn = build_augmenter(labels is not None)\n \n AUTO = tf.data.experimental.AUTOTUNE\n slices = paths if labels is None else (paths, labels)\n \n dset = tf.data.Dataset.from_tensor_slices(slices)\n dset = dset.map(decode_fn, num_parallel_calls=AUTO)\n dset = dset.cache(cache_dir) if cache else dset\n dset = dset.map(augment_fn, num_parallel_calls=AUTO) if augment else dset\n dset = dset.repeat() if repeat else dset\n dset = dset.shuffle(shuffle) if shuffle else dset\n dset = dset.batch(bsize).prefetch(AUTO)\n \n return dset\n\nstrategy = auto_select_accelerator()\nBATCH_SIZE = strategy.num_replicas_in_sync * 16\nGCS_DS_PATH = KaggleDatasets().get_gcs_path('ranzcr-clip-catheter-line-classification')\n\nprint(BATCH_SIZE)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\npaths = GCS_DS_PATH + \"/train/\" + df['StudyInstanceUID'] + '.jpg'\nsub_df = pd.read_csv('../input/ranzcr-clip-catheter-line-classification/sample_submission.csv')\ntest_paths = GCS_DS_PATH + \"/test/\" + sub_df['StudyInstanceUID'] + '.jpg'\n\n# Get the multi-labels\nlabel_cols = sub_df.columns[1:]\nlabels = df[label_cols].values\n\n# Train test split\n(\n train_paths, valid_paths, \n train_labels, valid_labels\n) = train_test_split(paths, labels, test_size=0.15, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Build the tensorflow datasets\nIMSIZES = (224, 240, 260, 300, 380, 456, 528, 600, 850)\nim_size = IMSIZES[3]\n\ndecoder = build_decoder(with_labels=True, target_size=(im_size, im_size))\ntest_decoder = build_decoder(with_labels=False, target_size=(im_size, im_size))\n\ntrain_dataset = build_dataset(\n train_paths, train_labels, bsize=BATCH_SIZE, decode_fn=decoder\n)\n\nvalid_dataset = build_dataset(\n valid_paths, valid_labels, bsize=BATCH_SIZE, decode_fn=decoder,\n repeat=False, shuffle=False, augment=False\n)\n\ntest_dataset = build_dataset(\n test_paths, cache=False, bsize=BATCH_SIZE, decode_fn=test_decoder,\n repeat=False, shuffle=False, augment=False\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Mixed Depthwise Convolution \n\nor **MDConv** is a core building block of [MixNet](https://arxiv.org/abs/1907.09595) architecture. This is the layer that we will add before the dense blocks. Generally, Mixed Convolution is a group of convolutions with varying filter sizes on each layer. Typically, we use constant filter dimention on per layer i.e `Conv2D`. But the group convolution is a concept where we use multiple size kernel on per layer. \n\n

\n \n

\n\nHowever, unlike `pytorch`, currently in `tf.keras`, still there is no such `group` parameter inside any convolutional layers that conveniently do this job. Instead we need to do something like as follows: \n\n

\n \n

\n\n\nThe original implementation I've found a bit messy to use for custom data set. However, here I will reuse their code bases and write a simple `tf.keras` layer for mix depth-wise convolution for easy use. It's actually pretty simple and with reasonable effort the whole `MixNet` is achievable. \n\n- [Official Code](https://github.com/tensorflow/tpu/tree/master/models/official/mnasnet/mixnet).\n\n---"},{"metadata":{},"cell_type":"markdown","source":"## 1. Mix-Depth Group Convolution\n\n\nThe code is adopted from official [code](https://github.com/tensorflow/tpu/tree/master/models/official/mnasnet/mixnet) and modified accordingly. To understand the operation, we must first understand what depthwise convolution is. The following diagram stated it at best. \n\n

\n \n

\n\n\nMixed Depthwise is simply an extend of it. Unlike a constant kernel size, above which is `5 x 5`, the `MixConv` uses multiple kernel sizes to the splitted channels, like below: `3 x 3`, `5 x 5` etc.\n\n

\n \n

"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"class MixDepthGroupConvolution2D(tf.keras.layers.Layer):\n def __init__(self, kernels=[3, 5],\n conv_kwargs=None,\n **kwargs):\n super(MixDepthGroupConvolution2D, self).__init__(**kwargs)\n\n if conv_kwargs is None:\n conv_kwargs = {\n 'strides': (1, 1),\n 'padding': 'same',\n 'dilation_rate': (1, 1),\n 'use_bias': False,\n }\n self.channel_axis = -1 \n self.kernels = kernels\n self.groups = len(self.kernels)\n self.strides = conv_kwargs.get('strides', (1, 1))\n self.padding = conv_kwargs.get('padding', 'same')\n self.dilation_rate = conv_kwargs.get('dilation_rate', (1, 1))\n self.use_bias = conv_kwargs.get('use_bias', False)\n self.conv_kwargs = conv_kwargs or {}\n\n self.layers = [tf.keras.layers.DepthwiseConv2D(kernels[i],\n strides=self.strides,\n padding=self.padding,\n activation=tf.nn.relu, \n dilation_rate=self.dilation_rate,\n use_bias=self.use_bias,\n kernel_initializer='he_normal')\n for i in range(self.groups)]\n\n def call(self, inputs, **kwargs):\n if len(self.layers) == 1:\n return self.layers[0](inputs)\n filters = K.int_shape(inputs)[self.channel_axis]\n splits = self.split_channels(filters, self.groups)\n x_splits = tf.split(inputs, splits, self.channel_axis)\n x_outputs = [c(x) for x, c in zip(x_splits, self.layers)]\n return tf.keras.layers.concatenate(x_outputs, \n axis=self.channel_axis)\n\n def split_channels(self, total_filters, num_groups):\n split = [total_filters // num_groups for _ in range(num_groups)]\n split[0] += total_filters - sum(split)\n return split\n\n def get_config(self):\n config = {\n 'kernels': self.kernels,\n 'groups': self.groups,\n 'strides': self.strides,\n 'padding': self.padding,\n 'dilation_rate': self.dilation_rate,\n 'conv_kwargs': self.conv_kwargs,\n }\n base_config = super(MixDepthGroupConvolution2D, self).get_config()\n return dict(list(base_config.items()) + list(config.items()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. Lamda Layer\n\nLamdaNetworks, basically an approach for the alternatives of Attention Mechanism. The main intuition is to capture the long-range interaction between an input and structured contextual information. \n\n\n

\n \n

\n\nHowever, currently the paper is under in double blind review and that's why there's no official implementation of it. If I'm not wrong they've showed a performance boost after adding the Lamda layer with `ResNet` family and `MobileNet`. Check the experimental details in page 17 and also 19 for hybrid lamda network. In the paper, it's claimed that the lamda layer can capture such long-range interaction by transforming available contexts into linear functions. Here, we'll be using an un-official implementaion of `LamdaLayer` from [here](https://github.com/lucidrains/lambda-networks). You can also find `PyTorch` code too. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import tensorflow as tf \nfrom tensorflow.keras import backend as K\n\nfrom einops.layers.tensorflow import Rearrange\nfrom tensorflow.keras import initializers\nfrom tensorflow import einsum, nn, meshgrid\nfrom tensorflow.keras.layers import Layer, InputSpec\nfrom tensorflow.keras import initializers as initializations\n\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.layers import Input, concatenate, ZeroPadding2D\nfrom tensorflow.keras.layers import Dense, Dropout, Activation\nfrom tensorflow.keras.layers import Conv2D, Conv3D, Softmax, Lambda, Add, Layer\nfrom tensorflow.keras.layers import AveragePooling2D, GlobalAveragePooling2D, MaxPooling2D\nfrom tensorflow.keras.layers import BatchNormalization","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# helpers functions\ndef exists(val):\n return val is not None\n\ndef default(val, d):\n return val if exists(val) else d\n\ndef calc_rel_pos(n):\n pos = tf.stack(meshgrid(tf.range(n), tf.range(n), indexing = 'ij'))\n pos = Rearrange('n i j -> (i j) n')(pos) # [n*n, 2] pos[n] = (i, j)\n rel_pos = pos[None, :] - pos[:, None] # [n*n, n*n, 2] rel_pos[n, m] = (rel_i, rel_j)\n rel_pos += n - 1 # shift value range from [-n+1, n-1] to [0, 2n-2]\n return rel_pos\n\n# lambda layer\nclass LambdaLayer(Layer):\n def __init__(\n self,\n *,\n dim_k,\n n = None,\n r = None,\n heads = 4,\n dim_out = None,\n dim_u = 1):\n super(LambdaLayer, self).__init__()\n '''\n Ref: https://github.com/lucidrains/lambda-networks/blob/main/lambda_networks/tfkeras.py\n '''\n\n self.out_dim = dim_out\n self.u = dim_u # intra-depth dimension\n self.heads = heads\n\n assert (dim_out % heads) == 0, 'values dimension must be divisible by number of heads for multi-head query'\n self.dim_v = dim_out // heads\n self.dim_k = dim_k\n self.heads = heads\n\n self.to_q = Conv2D(self.dim_k * heads, 1, use_bias=False)\n self.to_k = Conv2D(self.dim_k * dim_u, 1, use_bias=False)\n self.to_v = Conv2D(self.dim_v * dim_u, 1, use_bias=False)\n\n self.norm_q = BatchNormalization()\n self.norm_v = BatchNormalization()\n\n self.local_contexts = exists(r)\n if exists(r):\n assert (r % 2) == 1, 'Receptive kernel size should be odd'\n self.pos_conv = Conv3D(dim_k, (1, r, r), padding='same')\n else:\n assert exists(n), 'You must specify the window length (n = h = w)'\n rel_length = 2 * n - 1\n self.rel_pos_emb = self.add_weight(name='pos_emb',\n shape=(rel_length, rel_length, dim_k, dim_u),\n initializer=initializers.random_normal,\n trainable=True)\n self.rel_pos = calc_rel_pos(n)\n\n def call(self, x, **kwargs):\n b, hh, ww, c, u, h = *x.get_shape().as_list(), self.u, self.heads\n\n q = self.to_q(x)\n k = self.to_k(x)\n v = self.to_v(x)\n\n q = self.norm_q(q)\n v = self.norm_v(v)\n\n q = Rearrange('b hh ww (h k) -> b h k (hh ww)', h=h)(q)\n k = Rearrange('b hh ww (u k) -> b u k (hh ww)', u=u)(k)\n v = Rearrange('b hh ww (u v) -> b u v (hh ww)', u=u)(v)\n\n k = nn.softmax(k)\n\n Lc = einsum('b u k m, b u v m -> b k v', k, v)\n Yc = einsum('b h k n, b k v -> b n h v', q, Lc)\n\n if self.local_contexts:\n v = Rearrange('b u v (hh ww) -> b v hh ww u', hh=hh, ww=ww)(v)\n Lp = self.pos_conv(v)\n Lp = Rearrange('b v h w k -> b v k (h w)')(Lp)\n Yp = einsum('b h k n, b v k n -> b n h v', q, Lp)\n else:\n rel_pos_emb = tf.gather_nd(self.rel_pos_emb, self.rel_pos)\n Lp = einsum('n m k u, b u v m -> b n k v', rel_pos_emb, v)\n Yp = einsum('b h k n, b n k v -> b n h v', q, Lp)\n\n Y = Yc + Yp\n out = Rearrange('b (hh ww) h v -> b hh ww (h v)', hh = hh, ww = ww)(Y)\n return out\n\n def compute_output_shape(self, input_shape):\n return (*input_shape[:2], self.out_dim)\n\n def get_config(self):\n config = {'output_dim': (*self.input_shape[:2], self.out_dim)}\n base_config = super(LambdaLayer, self).get_config()\n return dict(list(base_config.items()) + list(config.items()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Dense Net 121\n\nFor such pre-trained model, direct source code from `tf.keras.app` would be nice but for concised implementaion I will be using from [this](https://github.com/flyyufelix/DenseNet-Keras) un-official implementaion. But the main problem is the author last updated his repo four years ago. So, lots of codes are old and won't fit with today's API. We need to modify many things. \n\nHere I'm porting the required code with proper upgradation according to the recent API. I'm choosing `DenseNet 121`, feel free to take others. However we remember our end network\n\n![LamdaDenseNet](https://user-images.githubusercontent.com/17668390/103923527-68532100-513f-11eb-97f5-d16806865f08.png)\n\nLet's build piece by piece. "},{"metadata":{},"cell_type":"markdown","source":"### Dense Block"},{"metadata":{"trusted":true},"cell_type":"code","source":"def dense_block(x, stage, nb_layers, nb_filter, \n growth_rate, dropout_rate=None, weight_decay=1e-4, grow_nb_filters=True):\n ''' Build a dense_block where the output of each conv_block is fed to subsequent ones\n # Arguments\n x: input tensor\n stage: index for dense block\n nb_layers: the number of layers of conv_block to append to the model.\n nb_filter: number of filters\n growth_rate: growth rate\n dropout_rate: dropout rate\n weight_decay: weight decay factor\n grow_nb_filters: flag to decide to allow number of filters to grow\n '''\n\n eps = 1.1e-5\n concat_feat = x\n\n for i in range(nb_layers):\n branch = i+1\n x = conv_block(concat_feat, stage, branch, \n growth_rate, dropout_rate, weight_decay)\n concat_feat = concatenate([concat_feat, x], \n axis=concat_axis, name='concat_'+str(stage)+'_'+str(branch))\n\n if grow_nb_filters:\n nb_filter += growth_rate\n\n return concat_feat, nb_filter","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Conv Block"},{"metadata":{"trusted":true},"cell_type":"code","source":"def conv_block(x, stage, branch, nb_filter, dropout_rate=None, weight_decay=1e-4):\n '''Apply BatchNorm, Relu, bottleneck 1x1 Conv2D, 3x3 Conv2D, and option dropout\n # Arguments\n x: input tensor \n stage: index for dense block\n branch: layer index within each dense block\n nb_filter: number of filters\n dropout_rate: dropout rate\n weight_decay: weight decay factor\n '''\n eps = 1.1e-5\n conv_name_base = 'conv' + str(stage) + '_' + str(branch)\n relu_name_base = 'relu' + str(stage) + '_' + str(branch)\n\n # 1x1 Convolution (Bottleneck layer)\n inter_channel = nb_filter * 4 \n x = BatchNormalization(epsilon=eps, axis=concat_axis, name=conv_name_base+'_x1_bn')(x)\n x = Scale(axis=concat_axis, name=conv_name_base+'_x1_scale')(x)\n x = Activation('relu', name=relu_name_base+'_x1')(x)\n x = Conv2D(inter_channel, 1, 1, name=conv_name_base+'_x1', use_bias=False)(x)\n\n if dropout_rate:\n x = Dropout(dropout_rate)(x)\n\n # 3x3 Convolution\n x = BatchNormalization(epsilon=eps, axis=concat_axis, name=conv_name_base+'_x2_bn')(x)\n x = Scale(axis=concat_axis, name=conv_name_base+'_x2_scale')(x)\n x = Activation('relu', name=relu_name_base+'_x2')(x)\n x = Conv2D(nb_filter, 1, 1, name=conv_name_base+'_x2', use_bias=False)(x) # 3, 3\n\n if dropout_rate:\n x = Dropout(dropout_rate)(x)\n\n return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Transition Block"},{"metadata":{"trusted":true},"cell_type":"code","source":"def transition_block(x, stage, nb_filter, compression=1.0, dropout_rate=None, weight_decay=1E-4):\n ''' Apply BatchNorm, 1x1 Convolution, averagePooling, optional compression, dropout \n # Arguments\n x: input tensor\n stage: index for dense block\n nb_filter: number of filters\n compression: calculated as 1 - reduction. Reduces the number of feature maps in the transition block.\n dropout_rate: dropout rate\n weight_decay: weight decay factor\n '''\n eps = 1.1e-5\n conv_name_base = 'conv' + str(stage) + '_blk'\n relu_name_base = 'relu' + str(stage) + '_blk'\n pool_name_base = 'pool' + str(stage) \n\n x = BatchNormalization(epsilon=eps, axis=concat_axis, name=conv_name_base+'_bn')(x)\n x = Scale(axis=concat_axis, name=conv_name_base+'_scale')(x)\n x = Activation('relu', name=relu_name_base)(x)\n x = Conv2D(int(nb_filter * compression), 1, 1, name=conv_name_base, use_bias=False)(x)\n\n if dropout_rate:\n x = Dropout(dropout_rate)(x)\n\n x = AveragePooling2D((2, 2), strides=(2, 2), name=pool_name_base)(x)\n\n return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**A Custom Layer for DenseNet used for BatchNormalization**"},{"metadata":{"trusted":true},"cell_type":"code","source":"class Scale(Layer):\n '''Custom Layer for DenseNet used for BatchNormalization.\n \n Learns a set of weights and biases used for scaling the input data.\n the output consists simply in an element-wise multiplication of the input\n and a sum of a set of constants:\n out = in * gamma + beta,\n where 'gamma' and 'beta' are the weights and biases larned.\n '''\n def __init__(self, weights=None, axis=-1, momentum = 0.9, beta_init='zero', gamma_init='one', **kwargs):\n self.momentum = momentum\n self.axis = axis\n self.beta_init = initializations.get(beta_init)\n self.gamma_init = initializations.get(gamma_init)\n self.initial_weights = weights\n super(Scale, self).__init__(**kwargs)\n\n def build(self, input_shape):\n self.input_spec = [InputSpec(shape=input_shape)]\n shape = (int(input_shape[self.axis]),)\n\n # Tensorflow >= 1.0.0 compatibility\n self.gamma = K.variable(self.gamma_init(shape), name='{}_gamma'.format(self.name))\n self.beta = K.variable(self.beta_init(shape), name='{}_beta'.format(self.name))\n #self.gamma = self.gamma_init(shape, name='{}_gamma'.format(self.name))\n #self.beta = self.beta_init(shape, name='{}_beta'.format(self.name))\n self._trainable_weights = [self.gamma, self.beta]\n\n if self.initial_weights is not None:\n self.set_weights(self.initial_weights)\n del self.initial_weights\n\n def call(self, x, mask=None):\n input_shape = self.input_spec[0].shape\n broadcast_shape = [1] * len(input_shape)\n broadcast_shape[self.axis] = input_shape[self.axis]\n out = K.reshape(self.gamma, broadcast_shape) * x + K.reshape(self.beta, broadcast_shape)\n return out\n\n def get_config(self):\n config = {\"momentum\": self.momentum, \"axis\": self.axis}\n base_config = super(Scale, self).get_config()\n return dict(list(base_config.items()) + list(config.items()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## MixDepthConvLamdaDenseNet\n\nLet's combine the above chunks and build the whole model. In the following function we will add `Dense-Block + Lamda-Layer + Transition-Layer` iteratively; just like shown in the below diagram.\n\n![LamdaDenseNet](https://user-images.githubusercontent.com/17668390/103923527-68532100-513f-11eb-97f5-d16806865f08.png)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def MixDepthConvLamdaDenseNet(image_size, \n nb_dense_block=3, # should 4 but set 3 for fast prototype \n growth_rate=32, \n nb_filter=64, reduction=0.0, \n dropout_rate=0.0, weight_decay=1e-4, \n classes=1000, weights_path=None):\n '''Instantiate the DenseNet 121 architecture,\n # Arguments\n nb_dense_block: number of dense blocks to add to end\n growth_rate: number of filters to add per dense block\n nb_filter: initial number of filters\n reduction: reduction factor of transition blocks.\n dropout_rate: dropout rate\n weight_decay: weight decay factor\n classes: optional number of classes to classify images\n weights_path: path to pre-trained weights\n # Returns\n A Keras model instance.\n '''\n eps = 1.1e-5\n\n # compute compression factor\n compression = 1.0 - reduction\n\n # Handle Dimension Ordering for different backends\n global concat_axis\n if K.image_data_format() == 'channels_last':\n concat_axis = -1\n img_input = Input(shape=(image_size, image_size, 3), name='data')\n else:\n concat_axis = 1\n img_input = Input(shape=(3, image_size, image_size), name='data')\n\n # From architecture for ImageNet (Table 1 in the paper)\n nb_filter = 64\n nb_layers = [6,12,24,16] # For DenseNet-121\n\n # Initial convolution\n x = MixDepthGroupConvolution2D(kernels=[3,5,7])(img_input)\n x = ZeroPadding2D((3, 3), name='conv1_zeropadding')(x)\n x = Conv2D(nb_filter, 7, 2, name='conv1', use_bias=False)(x)\n x = BatchNormalization(epsilon=eps, axis=concat_axis, name='conv1_bn')(x)\n x = Scale(axis=concat_axis, name='conv1_scale')(x)\n x = Activation('relu', name='relu1')(x)\n x = ZeroPadding2D((1, 1), name='pool1_zeropadding')(x)\n x = MaxPooling2D((3, 3), strides=(2, 2), name='pool1')(x)\n\n # Add dense blocks\n for block_idx in range(nb_dense_block - 1):\n stage = block_idx+2\n x, nb_filter = dense_block(x, stage, nb_layers[block_idx], nb_filter, \n growth_rate, dropout_rate=dropout_rate, weight_decay=weight_decay)\n \n # add lamda layer\n x = LambdaLayer(\n dim_out = nb_filter, # channels out\n r = 23, # the receptive field for relative positional encoding (23 x 23)\n dim_k = 16, # key dimension\n heads = 8, # number of heads, for multi-query; values dimension must be divisible by number of heads for multi-head query\n dim_u = 1 # 'intra-depth' dimension\n )(x)\n\n # Add transition_block\n x = transition_block(x, stage, nb_filter, compression=compression, \n dropout_rate=dropout_rate, weight_decay=weight_decay)\n nb_filter = int(nb_filter * compression)\n\n final_stage = stage + 1\n x, nb_filter = dense_block(x, final_stage, nb_layers[-1], nb_filter, \n growth_rate, dropout_rate=dropout_rate, weight_decay=weight_decay)\n\n x = BatchNormalization(epsilon=eps, axis=concat_axis, name='conv'+str(final_stage)+'_blk_bn')(x)\n x = Scale(axis=concat_axis, name='conv'+str(final_stage)+'_blk_scale')(x)\n x = Activation('relu', name='relu'+str(final_stage)+'_blk')(x)\n x = GlobalAveragePooling2D(name='pool'+str(final_stage))(x)\n \n x = Dense(classes, name='fc6')(x)\n x = Activation('sigmoid', name='prob')(x)\n\n model = Model(img_input, x, name='densenet')\n\n return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with strategy.scope():\n model = MixDepthConvLamdaDenseNet(\n image_size=im_size,\n reduction=0.5, \n classes=labels.shape[1])\n \n model.compile(\n optimizer=tf.keras.optimizers.Adam(learning_rate=0.1),\n loss='binary_crossentropy',\n metrics=[tf.keras.metrics.AUC(multi_label=True)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import time\n\nclass TimeHistory(tf.keras.callbacks.Callback):\n def on_train_begin(self, logs={}):\n self.times = []\n\n def on_epoch_begin(self, batch, logs={}):\n self.epoch_time_start = time.time()\n\n def on_epoch_end(self, batch, logs={}):\n self.times.append(time.time() - self.epoch_time_start)\n \n\ndef set_callback():\n checkpoint = tf.keras.callbacks.ModelCheckpoint(\n 'model.h5', save_best_only=True, \n save_weights_only=True,\n monitor='val_auc', mode='max')\n\n lr_reducer = tf.keras.callbacks.ReduceLROnPlateau(\n monitor=\"val_auc\", patience=3, min_lr=1e-6, mode='max')\n\n csv_loger = tf.keras.callbacks.CSVLogger('his.csv')\n \n time_his = TimeHistory()\n \n return [checkpoint, lr_reducer, csv_loger, time_his]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# print out the model params\ntrainable_count = np.sum([K.count_params(w) \\\n for w in model.trainable_weights]) \nnon_trainable_count = np.sum([K.count_params(w) \\\n for w in model.non_trainable_weights])\nprint('Total params: {:,}'.format(trainable_count + non_trainable_count))\nprint('Trainable params: {:,}'.format(trainable_count))\nprint('Non-trainable params: {:,}'.format(non_trainable_count))\n\nsteps_per_epoch = train_paths.shape[0] // BATCH_SIZE\n\nmodel.fit(\n train_dataset, \n epochs=10,\n verbose=1,\n callbacks=set_callback(),\n steps_per_epoch=steps_per_epoch,\n validation_data=valid_dataset)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = pd.read_csv('his.csv') \n\n# find the lowest validation loss score\nprint(history.loc[history['val_auc'].idxmin()])\nhistory.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(19,6))\n\nplt.subplot(121)\nplt.plot(history.epoch, history.loss, label=\"loss\")\nplt.plot(history.epoch, history.val_loss, label=\"val_loss\")\nplt.legend()\n\n\nplt.subplot(122)\nplt.plot(history.epoch, history.auc, label=\"auc\")\nplt.plot(history.epoch, history.val_auc, label=\"val_auc\")\nplt.legend()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model interpretability with Integrated Gradients"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_img_array(img_path, size=(300, 300)):\n # `img` is a PIL image of size 300\n img = tf.keras.preprocessing.image.load_img(img_path, target_size=size)\n # `array` is a float32 Numpy array of shape (300, 300, 3)\n array = tf.keras.preprocessing.image.img_to_array(img)\n # We add a dimension to transform our array into a \"batch\"\n # of size (1, 300, 300, 3)\n array = np.expand_dims(array, axis=0)\n return array","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.applications import densenet\n\ndef get_gradients(img_input, top_pred_idx):\n \"\"\"Computes the gradients of outputs w.r.t input image.\n\n Args:\n img_input: 4D image tensor\n top_pred_idx: Predicted label for the input image\n\n Returns:\n Gradients of the predictions w.r.t img_input\n \"\"\"\n images = tf.cast(img_input, tf.float32)\n\n with tf.GradientTape() as tape:\n tape.watch(images)\n preds = model(images)\n top_class = preds[:, top_pred_idx]\n\n grads = tape.gradient(top_class, images)\n return grads\n\n\ndef get_integrated_gradients(img_input, top_pred_idx, baseline=None, num_steps=50):\n \"\"\"Computes Integrated Gradients for a predicted label.\n\n Args:\n img_input (ndarray): Original image\n top_pred_idx: Predicted label for the input image\n baseline (ndarray): The baseline image to start with for interpolation\n num_steps: Number of interpolation steps between the baseline\n and the input used in the computation of integrated gradients. These\n steps along determine the integral approximation error. By default,\n num_steps is set to 50.\n\n Returns:\n Integrated gradients w.r.t input image\n \"\"\"\n # If baseline is not provided, start with a black image\n # having same size as the input image.\n if baseline is None:\n baseline = np.zeros(img_size).astype(np.float32)\n else:\n baseline = baseline.astype(np.float32)\n\n # 1. Do interpolation.\n img_input = img_input.astype(np.float32)\n interpolated_image = [\n baseline + (step / num_steps) * (img_input - baseline)\n for step in range(num_steps + 1)\n ]\n interpolated_image = np.array(interpolated_image).astype(np.float32)\n\n # 2. Preprocess the interpolated images\n interpolated_image = densenet.preprocess_input(interpolated_image)\n\n # 3. Get the gradients\n grads = []\n for i, img in enumerate(interpolated_image):\n img = tf.expand_dims(img, axis=0)\n grad = get_gradients(img, top_pred_idx=top_pred_idx)\n grads.append(grad[0])\n grads = tf.convert_to_tensor(grads, dtype=tf.float32)\n\n # 4. Approximate the integral using the trapezoidal rule\n grads = (grads[:-1] + grads[1:]) / 2.0\n avg_grads = tf.reduce_mean(grads, axis=0)\n\n # 5. Calculate integrated gradients and return\n integrated_grads = (img_input - baseline) * avg_grads\n return integrated_grads","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def random_baseline_integrated_gradients(\n img_input, top_pred_idx, num_steps=50, num_runs=2\n):\n \"\"\"Generates a number of random baseline images.\n\n Args:\n img_input (ndarray): 3D image\n top_pred_idx: Predicted label for the input image\n num_steps: Number of interpolation steps between the baseline\n and the input used in the computation of integrated gradients. These\n steps along determine the integral approximation error. By default,\n num_steps is set to 50.\n num_runs: number of baseline images to generate\n\n Returns:\n Averaged integrated gradients for `num_runs` baseline images\n \"\"\"\n # 1. List to keep track of Integrated Gradients (IG) for all the images\n integrated_grads = []\n\n # 2. Get the integrated gradients for all the baselines\n for run in range(num_runs):\n baseline = np.random.random(img_size) * 255\n igrads = get_integrated_gradients(\n img_input=img_input,\n top_pred_idx=top_pred_idx,\n baseline=baseline,\n num_steps=num_steps,\n )\n integrated_grads.append(igrads)\n\n # 3. Return the average integrated gradients for the image\n integrated_grads = tf.convert_to_tensor(integrated_grads)\n return tf.reduce_mean(integrated_grads, axis=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Gradient Visualization**"},{"metadata":{"trusted":true},"cell_type":"code","source":"class GradVisualizer:\n \"\"\"Plot gradients of the outputs w.r.t an input image.\"\"\"\n\n def __init__(self, positive_channel=None, negative_channel=None):\n if positive_channel is None:\n self.positive_channel = [0, 255, 0]\n else:\n self.positive_channel = positive_channel\n\n if negative_channel is None:\n self.negative_channel = [255, 0, 0]\n else:\n self.negative_channel = negative_channel\n\n def apply_polarity(self, attributions, polarity):\n if polarity == \"positive\":\n return np.clip(attributions, 0, 1)\n else:\n return np.clip(attributions, -1, 0)\n\n def apply_linear_transformation(\n self,\n attributions,\n clip_above_percentile=99.9,\n clip_below_percentile=70.0,\n lower_end=0.2,\n ):\n # 1. Get the thresholds\n m = self.get_thresholded_attributions(\n attributions, percentage=100 - clip_above_percentile\n )\n e = self.get_thresholded_attributions(\n attributions, percentage=100 - clip_below_percentile\n )\n\n # 2. Transform the attributions by a linear function f(x) = a*x + b such that\n # f(m) = 1.0 and f(e) = lower_end\n transformed_attributions = (1 - lower_end) * (np.abs(attributions) - e) / (\n m - e\n ) + lower_end\n\n # 3. Make sure that the sign of transformed attributions is the same as original attributions\n transformed_attributions *= np.sign(attributions)\n\n # 4. Only keep values that are bigger than the lower_end\n transformed_attributions *= transformed_attributions >= lower_end\n\n # 5. Clip values and return\n transformed_attributions = np.clip(transformed_attributions, 0.0, 1.0)\n return transformed_attributions\n\n def get_thresholded_attributions(self, attributions, percentage):\n if percentage == 100.0:\n return np.min(attributions)\n\n # 1. Flatten the attributions\n flatten_attr = attributions.flatten()\n\n # 2. Get the sum of the attributions\n total = np.sum(flatten_attr)\n\n # 3. Sort the attributions from largest to smallest.\n sorted_attributions = np.sort(np.abs(flatten_attr))[::-1]\n\n # 4. Calculate the percentage of the total sum that each attribution\n # and the values about it contribute.\n cum_sum = 100.0 * np.cumsum(sorted_attributions) / total\n\n # 5. Threshold the attributions by the percentage\n indices_to_consider = np.where(cum_sum >= percentage)[0][0]\n\n # 6. Select the desired attributions and return\n attributions = sorted_attributions[indices_to_consider]\n return attributions\n\n def binarize(self, attributions, threshold=0.001):\n return attributions > threshold\n\n def morphological_cleanup_fn(self, attributions, structure=np.ones((4, 4))):\n closed = ndimage.grey_closing(attributions, structure=structure)\n opened = ndimage.grey_opening(closed, structure=structure)\n return opened\n\n def draw_outlines(\n self, attributions, percentage=90, connected_component_structure=np.ones((3, 3))\n ):\n # 1. Binarize the attributions.\n attributions = self.binarize(attributions)\n\n # 2. Fill the gaps\n attributions = ndimage.binary_fill_holes(attributions)\n\n # 3. Compute connected components\n connected_components, num_comp = ndimage.measurements.label(\n attributions, structure=connected_component_structure\n )\n\n # 4. Sum up the attributions for each component\n total = np.sum(attributions[connected_components > 0])\n component_sums = []\n for comp in range(1, num_comp + 1):\n mask = connected_components == comp\n component_sum = np.sum(attributions[mask])\n component_sums.append((component_sum, mask))\n\n # 5. Compute the percentage of top components to keep\n sorted_sums_and_masks = sorted(component_sums, key=lambda x: x[0], reverse=True)\n sorted_sums = list(zip(*sorted_sums_and_masks))[0]\n cumulative_sorted_sums = np.cumsum(sorted_sums)\n cutoff_threshold = percentage * total / 100\n cutoff_idx = np.where(cumulative_sorted_sums >= cutoff_threshold)[0][0]\n if cutoff_idx > 2:\n cutoff_idx = 2\n\n # 6. Set the values for the kept components\n border_mask = np.zeros_like(attributions)\n for i in range(cutoff_idx + 1):\n border_mask[sorted_sums_and_masks[i][1]] = 1\n\n # 7. Make the mask hollow and show only the border\n eroded_mask = ndimage.binary_erosion(border_mask, iterations=1)\n border_mask[eroded_mask] = 0\n\n # 8. Return the outlined mask\n return border_mask\n\n def process_grads(\n self,\n image,\n attributions,\n polarity=\"positive\",\n clip_above_percentile=99.9,\n clip_below_percentile=0,\n morphological_cleanup=False,\n structure=np.ones((3, 3)),\n outlines=False,\n outlines_component_percentage=90,\n overlay=True,\n ):\n if polarity not in [\"positive\", \"negative\"]:\n raise ValueError(\n f\"\"\" Allowed polarity values: 'positive' or 'negative'\n but provided {polarity}\"\"\"\n )\n if clip_above_percentile < 0 or clip_above_percentile > 100:\n raise ValueError(\"clip_above_percentile must be in [0, 100]\")\n\n if clip_below_percentile < 0 or clip_below_percentile > 100:\n raise ValueError(\"clip_below_percentile must be in [0, 100]\")\n\n # 1. Apply polarity\n if polarity == \"positive\":\n attributions = self.apply_polarity(attributions, polarity=polarity)\n channel = self.positive_channel\n else:\n attributions = self.apply_polarity(attributions, polarity=polarity)\n attributions = np.abs(attributions)\n channel = self.negative_channel\n\n # 2. Take average over the channels\n attributions = np.average(attributions, axis=2)\n\n # 3. Apply linear transformation to the attributions\n attributions = self.apply_linear_transformation(\n attributions,\n clip_above_percentile=clip_above_percentile,\n clip_below_percentile=clip_below_percentile,\n lower_end=0.0,\n )\n\n # 4. Cleanup\n if morphological_cleanup:\n attributions = self.morphological_cleanup_fn(\n attributions, structure=structure\n )\n # 5. Draw the outlines\n if outlines:\n attributions = self.draw_outlines(\n attributions, percentage=outlines_component_percentage\n )\n\n # 6. Expand the channel axis and convert to RGB\n attributions = np.expand_dims(attributions, 2) * channel\n\n # 7.Superimpose on the original image\n if overlay:\n attributions = np.clip((attributions * 0.8 + image), 0, 255)\n return attributions\n\n def visualize(\n self,\n image,\n gradients,\n integrated_gradients,\n polarity=\"positive\",\n clip_above_percentile=99.9,\n clip_below_percentile=0,\n morphological_cleanup=False,\n structure=np.ones((3, 3)),\n outlines=False,\n outlines_component_percentage=90,\n overlay=True,\n figsize=(15, 8),\n ):\n # 1. Make two copies of the original image\n img1 = np.copy(image)\n img2 = np.copy(image)\n\n # 2. Process the normal gradients\n grads_attr = self.process_grads(\n image=img1,\n attributions=gradients,\n polarity=polarity,\n clip_above_percentile=clip_above_percentile,\n clip_below_percentile=clip_below_percentile,\n morphological_cleanup=morphological_cleanup,\n structure=structure,\n outlines=outlines,\n outlines_component_percentage=outlines_component_percentage,\n overlay=overlay,\n )\n\n # 3. Process the integrated gradients\n igrads_attr = self.process_grads(\n image=img2,\n attributions=integrated_gradients,\n polarity=polarity,\n clip_above_percentile=clip_above_percentile,\n clip_below_percentile=clip_below_percentile,\n morphological_cleanup=morphological_cleanup,\n structure=structure,\n outlines=outlines,\n outlines_component_percentage=outlines_component_percentage,\n overlay=overlay,\n )\n\n _, ax = plt.subplots(1, 3, figsize=figsize)\n ax[0].imshow(image)\n ax[1].imshow(grads_attr.astype(np.uint8))\n ax[2].imshow(igrads_attr.astype(np.uint8))\n\n ax[0].set_title(\"Input\")\n ax[1].set_title(\"Normal gradients\")\n ax[2].set_title(\"Integrated gradients\")\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_path = [\n '../input/ranzcr-clip-catheter-line-classification/test/1.2.826.0.1.3680043.8.498.10003659706701445041816900371598078663.jpg',\n '../input/ranzcr-clip-catheter-line-classification/test/1.2.826.0.1.3680043.8.498.10003890246067211044742686138544513464.jpg'\n]\n\n# 1. Convert the image to numpy array\nimg = get_img_array(img_path[0])\n\n# 2. Keep a copy of the original image\norig_img = np.copy(img[0]).astype(np.uint8)\n\n# 3. Preprocess the image\nimg_processed = tf.cast(densenet.preprocess_input(img), dtype=tf.float32)\n\n# 4. Get model predictions\npreds = model.predict(img_processed)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_size = (300, 300, 3)\n\ntop_pred_idx = tf.argmax(preds[0])\n\n# 5. Get the gradients of the last layer for the predicted label\ngrads = get_gradients(img_processed, top_pred_idx=top_pred_idx)\n\n# 6. Get the integrated gradients\nigrads = random_baseline_integrated_gradients(\n np.copy(orig_img), top_pred_idx=top_pred_idx, num_steps=50, num_runs=2\n)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}