{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# CLAM\n\nNOTE: Some of the descriptions or images are cited from: https://github.com/mahmoodlab/CLAM\n\n![img](https://github.com/mahmoodlab/CLAM/raw/master/docs/CLAM2.jpg)\n\n## TL;DR:\n\n+ CLAM is a high-throughput and interpretable method for data efficient whole slide image (WSI) classification using slide-level labels without any ROI extraction or patch-level annotations, and is capable of handling multi-class subtyping problems. Tested on three different WSI datasets, trained models adapt to independent test cohorts of WSI resections and biopsies as well as smartphone microscopy images (photomicrographs).\n+ paper: https://arxiv.org/abs/2004.09666\n\n## How to apply CLAM on the STRIP AI dataset ?\n\n+ I prepared four notebooks for pre-process, train and inference:\n\n### pre-process\n\n+ (1) image generation: https://www.kaggle.com/code/fx6300/clam-strip-ai-image-generation\n+ (2) feature extraction: https://www.kaggle.com/code/fx6300/clam-strip-ai-feature-extraction\n\n### train\n\n+ <b>&gt; THIS NOTEBOOK &lt;</b> (3) train: https://www.kaggle.com/code/fx6300/clam-strip-ai-train\n\n### inference\n\n+ (4) inference: https://www.kaggle.com/code/fx6300/clam-strip-ai-inference\n\n## How to visualize the attention generated by CLAM ?\n\n+ I prepared an example:\n  + https://www.kaggle.com/fx6300/clam-strip-ai-attention-heatmap\n\n## NOTE\n\n+ The source code from CLAM (https://github.com/mahmoodlab/CLAM) is licensed under GPLv3 and available for non-commercial academic purposes.","metadata":{}},{"cell_type":"code","source":"!conda install ../input/how-to-use-pyvips-offline/*.tar.bz2\n!git clone https://github.com/oval-group/smooth-topk.git\n!cd smooth-topk && python setup.py install","metadata":{"papermill":{"duration":42.042241,"end_time":"2022-09-04T13:26:23.047873","exception":false,"start_time":"2022-09-04T13:25:41.005632","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-09-27T04:06:50.271875Z","iopub.execute_input":"2022-09-27T04:06:50.272322Z","iopub.status.idle":"2022-09-27T04:07:00.891584Z","shell.execute_reply.started":"2022-09-27T04:06:50.272274Z","shell.execute_reply":"2022-09-27T04:07:00.890395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nsys.path.append(\"/opt/conda/lib/python3.7/site-packages/topk-1.0-py3.7.egg\")\nimport os\nimport gc\nimport cv2\nimport time\nimport random\nimport string\nimport joblib\nimport tifffile\nimport numpy as np \nimport pandas as pd \nimport torch\nfrom torch import nn\nfrom torchvision import models\nimport matplotlib.pyplot as plt\nfrom torch.utils.data import Dataset, DataLoader\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.model_selection import StratifiedKFold\nfrom tqdm.notebook import tqdm\nfrom tqdm import trange\nfrom torch.optim import lr_scheduler\nimport warnings\nimport tempfile\nfrom PIL import Image\nimport pyvips\nimport torch.nn.functional as F\nimport h5py\nfrom topk.svm import SmoothTop1SVM\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2022-09-27T04:07:00.894074Z","iopub.execute_input":"2022-09-27T04:07:00.894772Z","iopub.status.idle":"2022-09-27T04:07:00.903719Z","shell.execute_reply.started":"2022-09-27T04:07:00.894728Z","shell.execute_reply":"2022-09-27T04:07:00.902867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"debug = False\ntrain_df = pd.read_csv(\"../input/mayo-clinic-strip-ai/train.csv\")\ntest_df = pd.read_csv(\"../input/mayo-clinic-strip-ai/test.csv\")\ndirs = [\"../input/mayo-clinic-strip-ai/train/\", \"../input/mayo-clinic-strip-ai/test/\"]\nTILE_SIZE = 256\nIMG_HEIGHT = TILE_SIZE * 6\nIMG_WIDTH = TILE_SIZE * 6\nTARGET_FOLD = [0,1,2,3,4]#[0,1,2,3,4]","metadata":{"execution":{"iopub.status.busy":"2022-09-27T04:07:00.905367Z","iopub.execute_input":"2022-09-27T04:07:00.905795Z","iopub.status.idle":"2022-09-27T04:07:00.930975Z","shell.execute_reply.started":"2022-09-27T04:07:00.905759Z","shell.execute_reply":"2022-09-27T04:07:00.930091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ImgDataset(Dataset):\n    def __init__(self, df, data_dir, **kwargs):\n        self.df = df\n        self.data_dir = data_dir\n    def __len__(self):\n        return len(self.df)\n        \n    def __getitem__(self, index):\n        image_id = self.df.iloc[index].image_id\n        label = {\"CE\":0,\"LAA\":1}[self.df.iloc[index].label]\n        full_path = f\"{self.data_dir}/{self.df.iloc[index].image_id}.h5\"\n        with h5py.File(full_path,'r') as hdf5_file:\n            features = torch.stack([torch.tensor(hdf5_file[str(i)]) for i in range(64)]).view(64, 1024)\n            coords = torch.tensor([i for i in range(64)]).view(64)\n        return features, label, coords","metadata":{"execution":{"iopub.status.busy":"2022-09-27T04:07:00.933364Z","iopub.execute_input":"2022-09-27T04:07:00.93379Z","iopub.status.idle":"2022-09-27T04:07:00.941402Z","shell.execute_reply.started":"2022-09-27T04:07:00.933755Z","shell.execute_reply":"2022-09-27T04:07:00.940381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nAttention Network without Gating (2 fc layers)\nargs:\n    L: input feature dimension\n    D: hidden layer dimension\n    dropout: whether to use dropout (p = 0.25)\n    n_classes: number of classes \n\"\"\"\n\ndef initialize_weights(module):\n        for m in module.modules():\n                if isinstance(m, nn.Linear):\n                        nn.init.xavier_normal_(m.weight)\n                        m.bias.data.zero_()\n                \n                elif isinstance(m, nn.BatchNorm1d):\n                        nn.init.constant_(m.weight, 1)\n                        nn.init.constant_(m.bias, 0)\n\nclass Attn_Net(nn.Module):\n\n    def __init__(self, L = 1024, D = 256, dropout = False, n_classes = 1):\n        super(Attn_Net, self).__init__()\n        self.module = [\n            nn.Linear(L, D),\n            nn.Tanh()]\n\n        if dropout:\n            self.module.append(nn.Dropout(0.25))\n\n        self.module.append(nn.Linear(D, n_classes))\n        \n        self.module = nn.Sequential(*self.module)\n    \n    def forward(self, x):\n        return self.module(x), x # N x n_classes\n\n\"\"\"\nAttention Network with Sigmoid Gating (3 fc layers)\nargs:\n    L: input feature dimension\n    D: hidden layer dimension\n    dropout: whether to use dropout (p = 0.25)\n    n_classes: number of classes \n\"\"\"\nclass Attn_Net_Gated(nn.Module):\n    def __init__(self, L = 1024, D = 256, dropout = False, n_classes = 1):\n        super(Attn_Net_Gated, self).__init__()\n        self.attention_a = [\n            nn.Linear(L, D),\n            nn.Tanh()]\n        \n        self.attention_b = [nn.Linear(L, D),\n                            nn.Sigmoid()]\n        if dropout:\n            self.attention_a.append(nn.Dropout(0.25))\n            self.attention_b.append(nn.Dropout(0.25))\n\n        self.attention_a = nn.Sequential(*self.attention_a)\n        self.attention_b = nn.Sequential(*self.attention_b)\n        \n        self.attention_c = nn.Linear(D, n_classes)\n\n    def forward(self, x):\n        a = self.attention_a(x)\n        b = self.attention_b(x)\n        A = a.mul(b)\n        A = self.attention_c(A)  # N x n_classes\n        return A, x\n\n\"\"\"\nargs:\n    gate: whether to use gated attention network\n    size_arg: config for network size\n    dropout: whether to use dropout\n    k_sample: number of positive/neg patches to sample for instance-level training\n    dropout: whether to use dropout (p = 0.25)\n    n_classes: number of classes \n    instance_loss_fn: loss function to supervise instance-level training\n    subtyping: whether it's a subtyping problem\n\"\"\"\nclass CLAM_SB(nn.Module):\n    def __init__(self, gate = True, size_arg = \"small\", dropout = True, k_sample=8, n_classes=2,\n        instance_loss_fn=nn.CrossEntropyLoss(), subtyping=False):\n        super(CLAM_SB, self).__init__()\n        self.size_dict = {\"small\": [1024, 512, 256], \"big\": [1024, 512, 384]}\n        size = self.size_dict[size_arg]\n        fc = [nn.Linear(size[0], size[1]), nn.ReLU()]\n        if dropout:\n            fc.append(nn.Dropout(0.25))\n        if gate:\n            attention_net = Attn_Net_Gated(L = size[1], D = size[2], dropout = dropout, n_classes = 1)\n        else:\n            attention_net = Attn_Net(L = size[1], D = size[2], dropout = dropout, n_classes = 1)\n        fc.append(attention_net)\n        self.attention_net = nn.Sequential(*fc)\n        self.classifiers = nn.Linear(size[1], n_classes)\n        instance_classifiers = [nn.Linear(size[1], 2) for i in range(n_classes)]\n        self.instance_classifiers = nn.ModuleList(instance_classifiers)\n        self.k_sample = k_sample\n        self.instance_loss_fn = instance_loss_fn\n        self.n_classes = n_classes\n        self.subtyping = subtyping\n\n        initialize_weights(self)\n\n    def relocate(self):\n        device=torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n        self.attention_net = self.attention_net.to(device)\n        self.classifiers = self.classifiers.to(device)\n        self.instance_classifiers = self.instance_classifiers.to(device)\n    \n    @staticmethod\n    def create_positive_targets(length, device):\n        return torch.full((length, ), 1, device=device).long()\n    @staticmethod\n    def create_negative_targets(length, device):\n        return torch.full((length, ), 0, device=device).long()\n    \n    #instance-level evaluation for in-the-class attention branch\n    def inst_eval(self, A, h, classifier): \n        device=h.device\n        if len(A.shape) == 1:\n            A = A.view(1, -1)\n        top_p_ids = torch.topk(A, self.k_sample)[1][-1]\n        top_p = torch.index_select(h, dim=0, index=top_p_ids)\n        top_n_ids = torch.topk(-A, self.k_sample, dim=1)[1][-1]\n        top_n = torch.index_select(h, dim=0, index=top_n_ids)\n        p_targets = self.create_positive_targets(self.k_sample, device)\n        n_targets = self.create_negative_targets(self.k_sample, device)\n\n        all_targets = torch.cat([p_targets, n_targets], dim=0)\n        all_instances = torch.cat([top_p, top_n], dim=0)\n        logits = classifier(all_instances)\n        all_preds = torch.topk(logits, 1, dim = 1)[1].squeeze(1)\n        instance_loss = self.instance_loss_fn(logits, all_targets)\n        return instance_loss, all_preds, all_targets\n    \n    #instance-level evaluation for out-of-the-class attention branch\n    def inst_eval_out(self, A, h, classifier):\n        device=h.device\n        if len(A.shape) == 1:\n            A = A.view(1, -1)\n        top_p_ids = torch.topk(A, self.k_sample)[1][-1]\n        top_p = torch.index_select(h, dim=0, index=top_p_ids)\n        p_targets = self.create_negative_targets(self.k_sample, device)\n        logits = classifier(top_p)\n        p_preds = torch.topk(logits, 1, dim = 1)[1].squeeze(1)\n        instance_loss = self.instance_loss_fn(logits, p_targets)\n        return instance_loss, p_preds, p_targets\n\n    def forward(self, h, label=None, instance_eval=False, return_features=False, attention_only=False):\n        device = h.device\n        A, h = self.attention_net(h)  # NxK        \n        A = torch.transpose(A, 1, 0)  # KxN\n        if attention_only:\n            return A\n        A_raw = A\n        A = F.softmax(A, dim=1)  # softmax over N\n        total_inst_loss = torch.tensor(0.0).cuda()\n        if instance_eval:\n            all_preds = []\n            all_targets = []\n            inst_labels = F.one_hot(label, num_classes=self.n_classes).squeeze() #binarize label\n            for i in range(len(self.instance_classifiers)):\n                inst_label = inst_labels[i].item()\n                classifier = self.instance_classifiers[i]\n                if inst_label == 1: #in-the-class:\n                    instance_loss, preds, targets = self.inst_eval(A, h, classifier)\n                    all_preds.extend(preds.cpu().numpy())\n                    all_targets.extend(targets.cpu().numpy())\n                else: #out-of-the-class\n                    if self.subtyping:\n                        instance_loss, preds, targets = self.inst_eval_out(A, h, classifier)\n                        all_preds.extend(preds.cpu().numpy())\n                        all_targets.extend(targets.cpu().numpy())\n                    else:\n                        continue\n                total_inst_loss += instance_loss\n\n            if self.subtyping:\n                total_inst_loss /= len(self.instance_classifiers)\n        M = torch.mm(A, h) \n        logits = self.classifiers(M)\n        Y_hat = torch.topk(logits, 1, dim = 1)[1]\n        Y_prob = F.softmax(logits, dim = 1)\n        # print(total_inst_loss)\n        return logits, Y_prob, Y_hat, total_inst_loss, A_raw","metadata":{"execution":{"iopub.status.busy":"2022-09-27T04:07:00.943185Z","iopub.execute_input":"2022-09-27T04:07:00.943869Z","iopub.status.idle":"2022-09-27T04:07:00.978334Z","shell.execute_reply.started":"2022-09-27T04:07:00.943823Z","shell.execute_reply":"2022-09-27T04:07:00.977423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_model(model, train_loader, val_loader, criterions, optimizer, num_epochs, fold):\n    best_loss = 10000.0\n    best_acc = 0\n    for epoch in range(num_epochs):\n        val_loss = 0\n        model.cuda()\n        for phase in ['train', 'val']:\n            if phase == 'train':\n                model.train()\n            else:\n                model.eval()\n                \n            epoch_loss = 0.0\n            epoch_mcll = 0.0\n            epoch_acc = 0\n            \n            y_hat = []\n            y = []\n            for item in tqdm(train_loader if phase == \"train\" else val_loader, leave=False):\n                images = item[0][0].cuda()\n                classes = item[1].cuda().long()\n                optimizer.zero_grad()\n                \n                with torch.set_grad_enabled(phase == 'train'):\n                    logits, Y_prob, Y_hat, total_inst_loss, _ = model(images, label=classes[0], instance_eval=True)\n                    loss1 = criterions[0](Y_prob, F.one_hot(classes, num_classes=2).cuda().float())\n                    #loss2 = criterions[1](Y_prob, classes)\n                    loss3 = total_inst_loss\n                    loss = loss1 + loss3\n\n                    if phase == 'train':\n                        loss.backward()\n                        optimizer.step()\n\n                    y_hat += Y_hat\n                    y += classes.data\n                    epoch_loss += loss.item() * len(Y_hat)\n                    # print(Y_hat)\n                    # print(classes.data)\n                    epoch_acc += torch.sum(Y_hat[0] == classes.data) / len(Y_hat[0])\n            #epoch_mcll = weighted_multi_class_log_loss(s(torch.stack(y_hat).cpu()[:,:2]).detach(), torch.stack(y))\n            data_size = len(train_loader if phase == \"train\" else val_loader)\n            epoch_loss = epoch_loss / data_size\n            epoch_acc = epoch_acc.double() / data_size\n            #epoch_mcll = epoch_mcll / data_size\n            #epoch_mcll_val = epoch_mcll if phase == \"val\" else 100.0\n            if phase == 'val':\n                val_loss = epoch_loss\n\n            print(f'Epoch {epoch + 1}/{num_epochs} | {phase:^5} | Loss: {epoch_loss:.4f} | Acc: {epoch_acc:.4f}')\n        \n        if epoch_acc > best_acc:\n            traced = torch.jit.trace(model.cpu(), torch.rand(64, 1024))\n            traced.save('model_fold{}.pth'.format(fold))\n            best_acc = epoch_acc\n            del traced\n        elif abs(epoch_acc - best_acc) < 1e6 and best_loss > val_loss:\n            traced = torch.jit.trace(model.cpu(), torch.rand(64, 1024))\n            traced.save('model_fold{}.pth'.format(fold))\n            best_loss = val_loss\n            del traced\n        gc.collect()\n        torch.cuda.empty_cache()","metadata":{"papermill":{"duration":0.022161,"end_time":"2022-09-04T13:26:26.932243","exception":false,"start_time":"2022-09-04T13:26:26.910082","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-09-27T04:07:00.979935Z","iopub.execute_input":"2022-09-27T04:07:00.980304Z","iopub.status.idle":"2022-09-27T04:07:00.995656Z","shell.execute_reply.started":"2022-09-27T04:07:00.980267Z","shell.execute_reply":"2022-09-27T04:07:00.994635Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def weighted_multi_class_log_loss(y_hat, y, w):\n    # print(y)\n    # print(torch.clamp(torch.sum(y[:,0]), min=1e-15))\n    # print(torch.log(torch.clamp(y_hat, 1e-15, 1.0 - 1e-15)))\n    # print(y * torch.log(torch.clamp(y_hat, 1e-15, 1.0 - 1e-15)))\n    # print(0.5 * y * torch.log(torch.clamp(y_hat, 1e-15, 1.0 - 1e-15)))\n    a = torch.log(torch.clamp(y_hat, 1e-15, 1.0 - 1e-15)).cuda()\n    b = torch.tensor([torch.clamp(torch.sum(y[:,0]), min=1e-15), torch.clamp(torch.sum(y[:,1]), min=1e-15)]).cuda()\n    return torch.sum(-torch.sum(w * y * a * 1/b))\nweighted_multi_class_log_loss(torch.tensor([[0.5,0.5], [0.8,0.2], [0.0,1.0], [0.0,1.0]]).cuda(), F.one_hot(torch.tensor([0, 1, 0, 0]), 2).cuda(), torch.tensor([0.5,0.5]).cuda())","metadata":{"execution":{"iopub.status.busy":"2022-09-27T04:07:00.997209Z","iopub.execute_input":"2022-09-27T04:07:00.997597Z","iopub.status.idle":"2022-09-27T04:07:01.014624Z","shell.execute_reply.started":"2022-09-27T04:07:00.99754Z","shell.execute_reply":"2022-09-27T04:07:01.013598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class WeightedMultiClassLogLoss(torch.nn.Module):\n    def __init__(self, weights):\n        super(WeightedMultiClassLogLoss, self).__init__()\n        self.weights = weights\n\n    def forward(self, inputs, targets):\n        return weighted_multi_class_log_loss(inputs, targets, self.weights)","metadata":{"execution":{"iopub.status.busy":"2022-09-27T04:07:01.015849Z","iopub.execute_input":"2022-09-27T04:07:01.016275Z","iopub.status.idle":"2022-09-27T04:07:01.023616Z","shell.execute_reply.started":"2022-09-27T04:07:01.01624Z","shell.execute_reply":"2022-09-27T04:07:01.022629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"skf = StratifiedKFold(n_splits = 5, shuffle = True, random_state = 42)\n\ntrain09_df, test_df = train_test_split(train_df, test_size=0.1, random_state=42, stratify = train_df.label)\n\nfor fold, (train_idx, val_idx) in enumerate(skf.split(train09_df.index, train09_df.label)):\n    if not fold in TARGET_FOLD:\n        continue\n    loss_fn = SmoothTop1SVM(n_classes = 2).cuda()\n    model = CLAM_SB(instance_loss_fn = loss_fn, subtyping = True, k_sample = 4)\n    #model = timm.create_model(\"dm_nfnet_f0\")\n    train = train09_df.iloc[train_idx, :]\n    val = train09_df.iloc[val_idx, :]\n    ce_size = len(train[train['label'] == \"CE\"])\n    laa_size = len(train[train['label'] == \"LAA\"])\n    pre_weights = [(ce_size + laa_size) / ce_size, (ce_size + laa_size) / laa_size] \n    weights = torch.tensor(pre_weights / np.sum(pre_weights)).cuda()\n    batch_size = 1\n    train_loader = DataLoader(\n        ImgDataset(train, \"../input/my-features-v4/my_features-v4\"), \n        batch_size=batch_size, \n        shuffle=True, \n        num_workers=1\n    )\n    val_loader = DataLoader(\n        ImgDataset(val, \"../input/my-features-v4/my_features-v4\"), \n        batch_size=batch_size, \n        shuffle=False, \n        num_workers=1\n    )\n    criterions = [WeightedMultiClassLogLoss(weights), nn.CrossEntropyLoss()]\n    optimizer = torch.optim.AdamW(model.parameters(), lr=1e-5)\n    train_model(model, train_loader, val_loader, criterions, optimizer, 50, fold)\n    del model, train, val, train_loader, val_loader, criterions, optimizer\n    gc.collect()\n    torch.cuda.empty_cache()","metadata":{"papermill":{"duration":1584.98764,"end_time":"2022-09-04T13:52:52.65447","exception":false,"start_time":"2022-09-04T13:26:27.66683","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-09-27T04:07:01.024994Z","iopub.execute_input":"2022-09-27T04:07:01.025499Z","iopub.status.idle":"2022-09-27T04:13:21.559541Z","shell.execute_reply.started":"2022-09-27T04:07:01.025463Z","shell.execute_reply":"2022-09-27T04:13:21.557926Z"},"trusted":true},"execution_count":null,"outputs":[]}]}