{"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":"code","source":"!pip install faiss-gpu\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\n%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport os\nimport sys\nimport shutil\nimport math\nimport random\nimport heapq \nimport time\nimport copy\nimport itertools  \nfrom typing import Dict, List\nfrom PIL import Image\nfrom io import StringIO,BytesIO \nfrom scipy.spatial.distance import pdist\nfrom scipy.signal import butter, lfilter\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import label_binarize\nfrom sklearn.metrics import confusion_matrix,roc_curve,accuracy_score,auc,roc_auc_score \nfrom sklearn.decomposition import PCA\nfrom functools import reduce\nfrom scipy.io import loadmat\nimport cv2\nimport faiss \nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.autograd import Variable\nimport torchvision.ops as ops\n# torch.cuda.set_device(3)\nprint (torch.cuda.current_device())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\n# import numpy as np # linear algebra\n# import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# # Input data files are available in the read-only \"../input/\" directory\n# # For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tstart = time.time()\nroot_dir = '/kaggle/input/aptos2019-blindness-detection/train_images' #the path of images\nData = np.array(pd.read_csv(\"/kaggle/input/aptos2019-blindness-detection/train.csv\" , sep=','))\n#trainset \nnp.random.shuffle(Data)\ntrData=Data[0:int(0.8*len(Data))]\nteData=Data[int(0.8*len(Data)):]\n# trData=Data[0:40]\n# teData=Data[40:60]\ntrN, trI, trM, trY = [],[],[],[]\nfor iname, itype in trData:\n    try:\n        trN.append(iname)\n        trY.append(itype) #0 refer to Benign, and 1 refers to malignant\n        image_path = os.path.join(root_dir, iname + \".png\")\n        img = cv2.resize(cv2.imread(image_path).astype(np.float32), (256, 256))#(256,256,3)\n        trI.append(img)\n#         mask_path = os.path.join(root_dir,'mask', iname)\n#         mask = cv2.resize(cv2.imread(mask_path,0).astype(np.float32), (256, 256))/255#(256,256)\n#         trM.append(np.where(mask == 0.0, 0, 1)) #0, 1\n    except:\n        print(iname+\":\"+str(image_path))\n#     sys.stdout.write('\\r{} / {} '.format(len(trN),trData.shape[0]))\n#     sys.stdout.flush()\nprint('The length of trainset is %d'%len(trN))\n# testset\nteN, teI, teM, teY = [],[],[],[]\nfor iname, itype in teData:\n    try:\n        teN.append(iname)\n        teY.append(itype) #0 refer to Benign, and 1 refers to malignant\n        image_path = os.path.join(root_dir, iname + \".png\")\n        img = cv2.resize(cv2.imread(image_path).astype(np.float32), (256, 256))#(256,256,3)\n        teI.append(img)\n#         mask_path = os.path.join(root_dir,'mask', iname)\n#         mask = cv2.resize(cv2.imread(mask_path,0).astype(np.float32), (256, 256))/255#(256,256)\n#         teM.append(np.where(mask == 0.0, 0, 1))\n    except:\n        print(iname+\":\"+str(image_path))\n#     sys.stdout.write('\\r{} / {} '.format(len(teN),teData.shape[0]))\n#     sys.stdout.flush()\nprint('The length of testset is %d'%len(teN))\n\nelapsed = time.time() - tstart\nprint('Completed data handle in %d seconds' % int(elapsed))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"__all__ = ['ResNet', 'resnet18', 'resnet34', 'resnet50', 'resnet101',\n           'resnet152']\n\n\nmodel_urls = {\n    'resnet18': 'https://download.pytorch.org/models/resnet18-5c106cde.pth',\n    'resnet34': 'https://download.pytorch.org/models/resnet34-333f7ec4.pth',\n    'resnet50': 'https://download.pytorch.org/models/resnet50-19c8e357.pth',\n    'resnet101': 'https://download.pytorch.org/models/resnet101-5d3b4d8f.pth',\n    'resnet152': 'https://download.pytorch.org/models/resnet152-b121ed2d.pth',\n}\n#\n\ndef conv3x3(in_planes, out_planes, stride=1):\n    \"3x3 convolution with padding\"\n    return nn.Conv2d(in_planes, out_planes, kernel_size=3, stride=stride,\n                     padding=1, bias=False)\n\n\nclass BasicBlock(nn.Module):\n    expansion = 1\n\n    def __init__(self, inplanes, planes, stride=1, downsample=None):\n        super(BasicBlock, self).__init__()\n        self.conv1 = conv3x3(inplanes, planes, stride)\n        self.bn1 = nn.BatchNorm2d(planes)\n        self.relu = nn.ReLU(inplace=True)\n        self.conv2 = conv3x3(planes, planes)\n        self.bn2 = nn.BatchNorm2d(planes)\n        self.downsample = downsample\n        self.stride = stride\n\n    def forward(self, x):\n        residual = x\n\n        out = self.conv1(x)\n        out = self.bn1(out)\n        out = self.relu(out)\n\n        out = self.conv2(out)\n        out = self.bn2(out)\n\n        if self.downsample is not None:\n            residual = self.downsample(x)\n\n        out += residual\n        out = self.relu(out)\n\n        return out\n\n\nclass Bottleneck(nn.Module):\n    expansion = 4\n\n    def __init__(self, inplanes, planes, stride=1, downsample=None):\n        super(Bottleneck, self).__init__()\n        self.conv1 = nn.Conv2d(inplanes, planes, kernel_size=1, bias=False)\n        self.bn1 = nn.BatchNorm2d(planes)\n        self.conv2 = nn.Conv2d(planes, planes, kernel_size=3, stride=stride,\n                               padding=1, bias=False)\n        self.bn2 = nn.BatchNorm2d(planes)\n        self.conv3 = nn.Conv2d(planes, planes * 4, kernel_size=1, bias=False)\n        self.bn3 = nn.BatchNorm2d(planes * 4)\n        self.relu = nn.ReLU(inplace=True)\n        self.downsample = downsample\n        self.stride = stride\n\n    def forward(self, x):\n        residual = x\n\n        out = self.conv1(x)\n        out = self.bn1(out)\n        out = self.relu(out)\n\n        out = self.conv2(out)\n        out = self.bn2(out)\n        out = self.relu(out)\n\n        out = self.conv3(out)\n        out = self.bn3(out)\n\n        if self.downsample is not None:\n            residual = self.downsample(x)\n\n        out += residual\n        out = self.relu(out)\n\n        return out\n\n\nclass ResNet(nn.Module):\n\n    def __init__(self, block, layers, num_classes=1000):\n        self.inplanes = 64\n        super(ResNet, self).__init__()\n        self.conv1 = nn.Conv2d(3, 64, kernel_size=7, stride=2, padding=3,\n                               bias=False)\n        self.bn1 = nn.BatchNorm2d(64)\n        self.relu = nn.ReLU(inplace=True)\n        self.maxpool = nn.MaxPool2d(kernel_size=3, stride=2, padding=1)\n        self.layer1 = self._make_layer(block, 64, layers[0])\n        self.layer2 = self._make_layer(block, 128, layers[1], stride=2)\n        self.layer3 = self._make_layer(block, 256, layers[2], stride=2)\n        self.layer4 = self._make_layer(block, 512, layers[3], stride=2)\n        self.avgpool = nn.AdaptiveAvgPool2d(1)\n        self.fc = nn.Linear(512 * block.expansion, num_classes)\n\n        for m in self.modules():\n            if isinstance(m, nn.Conv2d):\n                n = m.kernel_size[0] * m.kernel_size[1] * m.out_channels\n                m.weight.data.normal_(0, math.sqrt(2. / n))\n            elif isinstance(m, nn.BatchNorm2d):\n                m.weight.data.fill_(1)\n                m.bias.data.zero_()\n\n    def _make_layer(self, block, planes, blocks, stride=1):\n        downsample = None\n        if stride != 1 or self.inplanes != planes * block.expansion:\n            downsample = nn.Sequential(\n                nn.Conv2d(self.inplanes, planes * block.expansion,\n                          kernel_size=1, stride=stride, bias=False),\n                nn.BatchNorm2d(planes * block.expansion),\n            )\n\n        layers = []\n        layers.append(block(self.inplanes, planes, stride, downsample))\n        self.inplanes = planes * block.expansion\n        for i in range(1, blocks):\n            layers.append(block(self.inplanes, planes))\n\n        return nn.Sequential(*layers)\n\n    def forward(self, x):\n        x = self.conv1(x)\n        x = self.bn1(x)\n        x = self.relu(x)\n        x = self.maxpool(x)\n\n        x = self.layer1(x)\n        x = self.layer2(x)\n        x = self.layer3(x)\n        feat = self.layer4(x)\n\n        x = self.avgpool(feat)\n        x = x.view(x.size(0), -1)\n        x = self.fc(x)\n\n        return feat, x\n\n\ndef resnet18(pretrained=False, **kwargs):\n    \"\"\"Constructs a ResNet-18 model.\n    Args:\n        pretrained (bool): If True, returns a model pre-trained on ImageNet\n    \"\"\"\n    model = ResNet(BasicBlock, [2, 2, 2, 2], **kwargs)\n    if pretrained:\n        model.load_state_dict(model_zoo.load_url(model_urls['resnet18']))\n    return model\n\n\ndef resnet34(pretrained=False, **kwargs):\n    \"\"\"Constructs a ResNet-34 model.\n    Args:\n        pretrained (bool): If True, returns a model pre-trained on ImageNet\n    \"\"\"\n    model = ResNet(BasicBlock, [3, 4, 6, 3], **kwargs)\n    if pretrained:\n        model.load_state_dict(model_zoo.load_url(model_urls['resnet34']))\n    return model\n\n\ndef resnet50(pretrained=False, **kwargs):\n    \"\"\"Constructs a ResNet-50 model.\n    Args:\n        pretrained (bool): If True, returns a model pre-trained on ImageNet\n    \"\"\"\n    model = ResNet(Bottleneck, [3, 4, 6, 3], **kwargs)\n    if pretrained:\n        model.load_state_dict(model_zoo.load_url(model_urls['resnet50']))\n    return model\n\n\ndef resnet101(pretrained=False, **kwargs):\n    \"\"\"Constructs a ResNet-101 model.\n    Args:\n        pretrained (bool): If True, returns a model pre-trained on ImageNet\n    \"\"\"\n    model = ResNet(Bottleneck, [3, 4, 23, 3], **kwargs)\n    if pretrained:\n        model.load_state_dict(model_zoo.load_url(model_urls['resnet101']))\n    return model\n\n\ndef resnet152(pretrained=False, **kwargs):\n    \"\"\"Constructs a ResNet-152 model.\n    Args:\n        pretrained (bool): If True, returns a model pre-trained on ImageNet\n    \"\"\"\n    model = ResNet(Bottleneck, [3, 8, 36, 3], **kwargs)\n    if pretrained:\n        model.load_state_dict(model_zoo.load_url(model_urls['resnet152']))\n    return model","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#define model\nmodel = resnet34(num_classes=5).cuda()#initialize model\noptimizer = torch.optim.Adam(model.parameters(), lr=0.001) #define optimizer\nce_loss  = nn.CrossEntropyLoss().cuda() #define cross-entropy loss\n#train model\nbest_net, best_loss = None, float('inf')\nbatchSize = 128\nfor epoch in range(30):#iteration\n    losses = []\n    shuffled_idx = np.random.permutation(np.arange(len(trI)))\n    trI_batch = np.array(trI)[shuffled_idx]\n    trY_batch = np.array(trY)[shuffled_idx]\n    num_batches = (len(trI)+batchSize-1) // batchSize\n    for i in range(num_batches):\n        optimizer.zero_grad()#grad vanish\n        min_idx = i * batchSize\n        max_idx = np.min([len(trI), (i+1)*batchSize])\n        X_batch = torch.from_numpy(trI_batch[min_idx:max_idx]).type(torch.FloatTensor).cuda()\n        Y_batch = torch.from_numpy(trY_batch[min_idx:max_idx]).type(torch.LongTensor).cuda()\n        #forword\n\n        _,Out_batch = model(X_batch.permute(0, 3, 1, 2))#permute the dims of matrix\n        loss = ce_loss(Out_batch,Y_batch)#loss\n        #backward\n        loss.backward()\n        optimizer.step()#update parameters\n        #show loss\n        sys.stdout.write('\\r {} / {} : loss = {}'.format(i+1, num_batches, float('%0.6f'%loss.item())))\n        sys.stdout.flush()     \n        losses.append(loss.item())\n    print(\"Eopch: %5d mean_loss = %.6f\" % (epoch + 1, np.mean(losses)))\n    if np.mean(losses) < best_loss:\n        best_loss = np.mean(losses)\n        best_net = copy.deepcopy(model)\n#     if(epoch % 5 == 0 and epoch > 0)\n#         lr\nprint(\"best_loss = %.6f\" % (best_loss))\n#release gpu memory\nmodel = model.cpu()\nce_loss = ce_loss.cpu()\ntorch.cuda.empty_cache()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#test model\nteY_pred = []\nteY_prob = []\nnum_batches = (len(teI)+batchSize-1) // batchSize\nfor i in range(num_batches):\n    min_idx = i * batchSize\n    max_idx = np.min([len(teI), (i+1)*batchSize])\n    x_batch = torch.from_numpy(np.array(teI[min_idx:max_idx])).type(torch.FloatTensor).cuda()\n\n    _,out_batch = best_net(x_batch.permute(0, 3, 1, 2))#forword\n    \n    out_batch = F.log_softmax(out_batch,dim=1)\n#     prob = out_batch.max(1,keepdim=True)[0]\n    teY_prob.extend(out_batch.cpu().data.numpy().tolist())\n    pred = out_batch.max(1,keepdim=True)[1]\n    teY_pred.extend(pred.cpu().data.numpy().flatten().tolist())\n    sys.stdout.write('\\r {} / {} '.format(i, num_batches))\n    sys.stdout.flush()\n    \n# TNR= TN / (FP+TN) ->low misdiagnosis rate->Specificity\n# TPR= TP / (TP+FN) -> low missed diagnosis rate->Sensitivity\n# ROC curves: y axis:Sensitivity, x axis:1-Specificity\n# confusion matrix\nlabels = [0,1,2,3,4]\ncm = confusion_matrix(teY, teY_pred, labels=labels) \n\n    \nprint ('Sensitivity(TPR) of Benign: %.6f'%float(cm[0][0]/np.sum(cm[0]))) \nprint ('Sensitivity(TPR) of Malignant: %.6f'%float(cm[1][1]/np.sum(cm[1])))\n# auc and roc\nteY_one_hot = label_binarize(np.array(teY), np.arange(len(labels)))\nauc_score = roc_auc_score(teY_one_hot, np.array(teY_prob), average='micro')#macro\nprint ('AUC (Area Under Curve) of Micro: %.6f'% auc_score)\n#plot roc curve\nfpr_tce, tpr_tce, thresholds = roc_curve(teY_one_hot.ravel(),np.array(teY_prob).ravel()) \n# plt.plot(fpr_ce, tpr_ce, c = 'r', ls = '--', label = u'ATH(our) AUC=%.4f' % auc_score)\nplt.plot(fpr_tce, tpr_tce, c = 'b', ls = '--', label = u'R-MAC AUC=%.4f' % auc_score) \nplt.plot((0, 1), (0, 1), c = '#808080', lw = 1, ls = '--', alpha = 0.7)\nplt.xlim((-0.01, 1.02))\nplt.ylim((-0.01, 1.02))\nplt.xticks(np.arange(0, 1.1, 0.2))\nplt.yticks(np.arange(0, 1.1, 0.2))\nplt.xlabel('1-Specificity')\nplt.ylabel('Sensitivity')\nplt.grid(b=True, ls=':')\nplt.legend(loc='lower right')\nplt.title('TNSCUI')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#release gpu memory and save model in CPU\nmodel = model.cpu()\nce_loss = ce_loss.cpu()\nbest_net = best_net.cpu()\ntorch.cuda.empty_cache() \ntorch.save(best_net.state_dict(), '/kaggle/working/CAM_resnet_34.pkl')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n#Code: https://github.com/imatge-upc/retrieval-2017-cam\n#Paper: BMVC2017《Class-Weighted Convolutional Features for Image Retrieval》\n'''\n# Extract region of interest from CAMs\ndef extract_ROI(heatmap, threshold):\n    th = threshold * np.max(heatmap)\n    heatmap = heatmap > th\n    # Find the largest connected component\n\n    contours, hierarchy = cv2.findContours(heatmap.astype('uint8'), mode=cv2.RETR_EXTERNAL, method=cv2.CHAIN_APPROX_SIMPLE)\n\n    areas = [cv2.contourArea(ctr) for ctr in contours]\n\n    max_contour = contours[areas.index(max(areas))]\n\n    x, y, w, h = cv2.boundingRect(max_contour)\n    if w == 0:\n        w = heatmap.shape[1]\n    if h == 0:\n        h = heatmap.shape[0]\n    return x, y, w, h\n\ndef compute_crow_channel_weight(X):\n    \"\"\"\n    Given a tensor of features, compute channel weights as the\n    log of inverse channel sparsity.\n    :param ndarray X:\n        3d tensor of activations with dimensions (channels, height, width)\n    :returns ndarray:\n        a channel weight vector\n    \"\"\"\n    K, w, h = X.shape\n    area = float(w * h)\n    nonzeros = np.zeros(K, dtype=np.float32)\n    for i, x in enumerate(X):\n        nonzeros[i] = np.count_nonzero(x) / area\n\n    nzsum = nonzeros.sum()\n    for i, d in enumerate(nonzeros):\n        nonzeros[i] = np.log(nzsum / d) if d > 0. else 0.\n\n    return nonzeros\n\ndef compute_pca(descriptors, pca_dim=512, whiten=True):\n    print (descriptors.shape)\n    t1 = time.time()\n    print( 'Computing PCA with dimension reduction to: ', pca_dim)\n    sys.stdout.flush()\n    pca = PCA(n_components=pca_dim, whiten=whiten)\n    pca.fit(descriptors)\n    print (pca.components_.shape)\n    print ('PCA finished!')\n    print ('Time elapsed computing PCA: ', time.time() - t1)\n    return pca\n\n\ndef sum_pooling(features):\n    num_samples = features.shape[0]\n    num_features = features.shape[1]\n    sys.stdout.flush()\n    descriptors = np.zeros((num_samples, num_features), dtype=np.float32)\n    for i in range(0, num_samples):\n        #print 'Image: ', i\n        #sys.stdout.flush()\n        for f in range(0, num_features):\n            descriptors[i, f] = features[i, f].sum()\n    descriptors /= np.linalg.norm(descriptors, axis=1)[:, None]\n    return descriptors\n\ndef weighted_cam_pooling(features, cams, max_pool=False, channel_weights=True):\n    '''\n    :param features: Feature Maps\n    :param cams: Class Activation Maps\n    :param max_pool: Perform also Max pooling\n    :param channel_weights: Channel Weighting as in Crow\n    :return: A descriptor for each CAM.\n    '''\n    t = time.time()\n    num_samples = features.shape[0]\n    num_features = features.shape[1]\n    num_classes = cams.shape[1]\n\n    wp_regions = np.zeros((num_features, num_classes), dtype=np.float32)\n    wsp_descriptors_reg = np.zeros((num_samples * num_classes, num_features), dtype=np.float32)\n    wmp_descriptors_reg = np.zeros((num_samples * num_classes, num_features), dtype=np.float32)\n\n    if max_pool:\n        mp_regions = np.zeros((num_features, num_classes), dtype=np.float32)\n\n    for i in range(0, num_samples):\n        #CROW\n        if channel_weights:\n            C = np.array(compute_crow_channel_weight(features[i]))\n\n        for f in range(0, num_features):\n            for k in range(0, num_classes):\n                # For each region compute avg weighted sum of activations and l2 normalize\n                if max_pool:\n                        mp_regions[f, k] = np.amax(np.multiply(features[i, f], cams[i, k]))\n                wp_regions[f, k] = np.multiply(features[i, f], cams[i, k]).sum()\n\n        if channel_weights:\n            wp_regions = wp_regions * C[:, None]\n        wp_regions /= np.linalg.norm(wp_regions, axis=0)\n\n        if max_pool:\n            if channel_weights:\n                mp_regions = mp_regions * C[:, None]\n            mp_regions /= np.linalg.norm(mp_regions, axis=0)\n\n        wsp_descriptors_reg[num_classes*i:num_classes*(i+1)] = np.transpose(wp_regions)\n\n        if max_pool:\n            wmp_descriptors_reg[num_classes*i:num_classes*(i+1)] = np.transpose(mp_regions)\n\n    #print 'Time elapsed computing image representations for the batch: ', time.time() - t\n\n    if max_pool:\n        return wsp_descriptors_reg, wmp_descriptors_reg\n    else:\n        return wsp_descriptors_reg\n    \n# General Descriptor Aggregation : PCA + Aggregation\ndef descriptor_aggregation(descriptors_cams, num_images, num_classes, pca=None):\n\n    num_classes_ori = int(descriptors_cams.shape[0] / num_images)\n    descriptors = np.zeros((num_images, descriptors_cams.shape[1]), dtype=np.float32)\n\n    if pca is not None:\n        # Sometimes we may have errors during re-ranking due to bounding box generation on places where CAM=0\n        try:\n            descriptors_pca = pca.transform(descriptors_cams)\n        except:\n            print ('---------------------------->Exception')\n            desc_err = np.zeros((descriptors_cams.shape[0], descriptors_cams.shape[1]), dtype=np.float32)\n            for j in range(0, descriptors_cams.shape[0]):\n                try:\n                    desc_err[j] = pca.transform(descriptors_cams[j])\n                except:\n                    print ('---------------------------->Exception')\n                    print (j)\n                    desc_err[j] = desc_err[j-1]\n            descriptors_pca = desc_err\n\n        descriptors = np.zeros((num_images, descriptors_pca.shape[1]), dtype=np.float32)\n        #print descriptors_pca.shape\n\n    index = 0\n    for i in range(0, num_images):\n        index = num_classes_ori + index\n        if i == 0:\n            index = 0\n        if pca is not None:\n            for k in range(index, index+num_classes):\n                descriptors_pca[k] /= np.linalg.norm(descriptors_pca[k])\n                descriptors[i] += descriptors_pca[k]\n\n            descriptors[i] /= np.linalg.norm(descriptors[i])\n        else:\n            for k in range(index, index+num_classes):\n                descriptors[i] += descriptors_cams[k]\n            descriptors[i] /= np.linalg.norm(descriptors[i])\n\n    return descriptors\n\n#load model and transfer to GPU\ndevice = torch.device(\"cuda\")\nbest_net = resnet34(num_classes=5)\nbest_net.load_state_dict(torch.load( '/kaggle/working/CAM_resnet_34.pkl'))\nbest_net.to(device)\n\nweights_fc = best_net.fc.state_dict()['weight'] #tensor([2, 2048])\nweights_fc = np.transpose(weights_fc.cpu().numpy(), (1, 0)) #numpy(2048, 2)\n#extract feature-cam for trainset\nbatchSize=25\nnum_samples = len(trI)\nn_classes = 5\nclass_list = np.zeros((num_samples, n_classes), dtype=np.int32)\nscores_vec = list()\n# Region of interest for re-ranking (bounding box coordinates --> (num samples, num_thresholds, x,y,dx,dy)\nbbox_coord = np.zeros((num_samples, 5, 4), dtype=np.int16)\nnum_batches = (len(trI)+batchSize-1) // batchSize\nfor i in range(num_batches):\n    min_idx = i * batchSize\n    max_idx = np.min([len(trI), (i+1)*batchSize])\n    x_batch = torch.from_numpy(np.array(trI[min_idx:max_idx])).type(torch.FloatTensor).cuda()\n    if x_batch.shape[0] == 0:\n        continue\n    conv_outputs_cam, scores = best_net(x_batch.permute(0, 3, 1, 2))#forword\n    conv_outputs_cam = conv_outputs_cam.data.cpu().numpy()\n    scores = scores.data.cpu().numpy()\n    if i == 0:\n        features_conv = np.zeros((num_samples, conv_outputs_cam.shape[1], conv_outputs_cam.shape[2], conv_outputs_cam.shape[3]))\n        cams = np.zeros((num_samples, n_classes, conv_outputs_cam.shape[2], conv_outputs_cam.shape[3]))\n        features_conv[min_idx:max_idx, :, :, :] = conv_outputs_cam\n    else:\n        features_conv[min_idx:max_idx, :, :, :] = conv_outputs_cam\n    \n    for ii in range(0, max_idx-min_idx):\n        indexed_scores = scores[ii].argsort()[::-1]\n        scores_vec.append(scores[ii])\n        for k in range(0, n_classes):\n            w_class = weights_fc[:, indexed_scores[k]]\n            cam = np.zeros(dtype=np.float32, shape=conv_outputs_cam.shape[2:4])\n            for ind, w in enumerate(w_class):\n                cam += w * conv_outputs_cam[ii, ind, :, :]\n            cam /= np.max(cam)\n            cam[np.where(cam < 0)] = 0\n            cams[min_idx+ii, k, :, :] = cam\n            class_list[min_idx+ii, k] = indexed_scores[k]\n            \n        heatmap = cams[min_idx+ii, 0]\n        bbox_coord[min_idx+ii, 0, :] = extract_ROI(heatmap=heatmap, threshold=0.01)# Full Image\n        bbox_coord[min_idx+ii, 1, :] = extract_ROI(heatmap=heatmap, threshold=0.1)\n        bbox_coord[min_idx+ii, 2, :] = extract_ROI(heatmap=heatmap, threshold=0.2)\n        bbox_coord[min_idx+ii, 3, :] = extract_ROI(heatmap=heatmap, threshold=0.3)\n        bbox_coord[min_idx+ii, 4, :] = extract_ROI(heatmap=heatmap, threshold=0.4)\n#output: features_conv, cams, class_list, scores_vec, bbox_coord\n#extract features    \nd_wp = weighted_cam_pooling(features_conv, cams)\n#d_sp = sum_pooling(features_conv)\n# Compute Query Descriptor\ntrF = descriptor_aggregation(d_wp, num_samples, n_classes)\n# trF_pca = compute_pca(trF)\n\n#extract feature-cam for testset\nbatchSize=25\nnum_samples = len(teI)\nn_classes = 5\nclass_list = np.zeros((num_samples, n_classes), dtype=np.int32)\nscores_vec = list()\n# Region of interest for re-ranking (bounding box coordinates --> (num samples, num_thresholds, x,y,dx,dy)\nbbox_coord = np.zeros((num_samples, 5, 4), dtype=np.int16)\nnum_batches = (len(teI)+batchSize-1) // batchSize\nfor i in range(num_batches):\n    min_idx = i * batchSize\n    max_idx = np.min([len(teI), (i+1)*batchSize])\n    x_batch = torch.from_numpy(np.array(teI[min_idx:max_idx])).type(torch.FloatTensor).cuda()\n    if x_batch.shape[0] == 0:\n        continue\n    conv_outputs_cam, scores = best_net(x_batch.permute(0, 3, 1, 2))#forword\n    conv_outputs_cam = conv_outputs_cam.data.cpu().numpy()\n    scores = scores.data.cpu().numpy()\n    if i == 0:\n        features_conv = np.zeros((num_samples, conv_outputs_cam.shape[1], conv_outputs_cam.shape[2], conv_outputs_cam.shape[3]))\n        cams = np.zeros((num_samples, n_classes, conv_outputs_cam.shape[2], conv_outputs_cam.shape[3]))\n        features_conv[min_idx:max_idx, :, :, :] = conv_outputs_cam\n    else:\n        features_conv[min_idx:max_idx, :, :, :] = conv_outputs_cam\n    \n    for ii in range(0, max_idx-min_idx):\n        indexed_scores = scores[ii].argsort()[::-1]\n        scores_vec.append(scores[ii])\n        for k in range(0, n_classes):\n            w_class = weights_fc[:, indexed_scores[k]]\n            cam = np.zeros(dtype=np.float32, shape=conv_outputs_cam.shape[2:4])\n            for ind, w in enumerate(w_class):\n                cam += w * conv_outputs_cam[ii, ind, :, :]\n            cam /= np.max(cam)\n            cam[np.where(cam < 0)] = 0\n            cams[min_idx+ii, k, :, :] = cam\n            class_list[min_idx+ii, k] = indexed_scores[k]\n            \n        heatmap = cams[min_idx+ii, 0]\n        bbox_coord[min_idx+ii, 0, :] = extract_ROI(heatmap=heatmap, threshold=0.01)# Full Image\n        bbox_coord[min_idx+ii, 1, :] = extract_ROI(heatmap=heatmap, threshold=0.1)\n        bbox_coord[min_idx+ii, 2, :] = extract_ROI(heatmap=heatmap, threshold=0.2)\n        bbox_coord[min_idx+ii, 3, :] = extract_ROI(heatmap=heatmap, threshold=0.3)\n        bbox_coord[min_idx+ii, 4, :] = extract_ROI(heatmap=heatmap, threshold=0.4)\n#output: features_conv, cams, class_list, scores_vec, bbox_coord\n#extract features    \nd_wp = weighted_cam_pooling(features_conv, cams)\n#d_sp = sum_pooling(features_conv)\n# Compute Query Descriptor\nteF = descriptor_aggregation(d_wp, num_samples, n_classes)\n#teF_pca = compute_pca(teF)\n\n#compute the size of lesion\ndef Func_IOU_size(pred,target):\n    ious = []\n    # ignore IOU for background class\n    pred_inds = pred != 0\n    pred_sum = pred_inds.sum()\n    target_inds = target != 0\n    target_sum = target_inds.sum()\n    ious.append(round(float(min(pred_sum,target_sum)/max(pred_sum,target_sum)),4))\n    return np.mean(ious)\n# buliding index of trainset\ntstart = time.time()\ncpu_index = faiss.IndexFlatL2(np.shape(trF)[1]) #\ngpu_index = faiss.index_cpu_to_all_gpus(cpu_index) #make all gpu usable\nnp.savetxt(\"/kaggle/working/index.txt\", trF,fmt='%f',delimiter=',')\ngpu_index.add(np.ascontiguousarray(trF, dtype=np.float32)) #add data(must be float32) to index\nelapsed = time.time() - tstart    \nprint('Completed buliding index in %d seconds' % int(elapsed))\nfor topk in [5, 10,25,50]:\n    mAP = [] #mean average precision\n    mIoU = []\n    scores, neighbors = gpu_index.search(np.array(teF).astype('float32'), k=topk)\n    for i, teVal in enumerate(teF):\n        stype = teY[i]\n        #perfromance\n        pos_len = 0\n        rank_len = len(neighbors[i])\n        #for j in ranklist:\n        for j in neighbors[i].tolist():\n            dtype = trY[j]\n            rank_len=rank_len+1\n            if stype==dtype:  #hit\n                pos_len += 1\n        mAP.append(pos_len/rank_len)\n#             mIoU.append(Func_IOU_size(teM[i],trM[j]))\n    print(\"mAP of top {} ={:.4f}\".format(topk,np.mean(mAP)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}