{"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":"import os\nimport pandas as pd\nimport numpy as np\nfrom sklearn.model_selection import train_test_split\nfrom pydicom import dcmread\nfrom PIL import Image\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\n\n\nimport torch\nimport torch.nn as nn\nimport torchvision\nimport torchvision.transforms as transforms\nfrom torch.utils import data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparing labels","metadata":{}},{"cell_type":"code","source":"label_data = pd.read_csv('../input/rsna-pneumonia-detection-challenge/stage_2_train_labels.csv')\ncolumns = ['patientId', 'Target']\n\nlabel_data = label_data.filter(columns)\nlabel_data.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dividing labels for train and validation set","metadata":{}},{"cell_type":"code","source":"train_labels, val_labels = train_test_split(label_data.values, test_size=0.1)\nprint(train_labels.shape)\nprint(val_labels.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'patientId: {train_labels[0][0]}, Target: {train_labels[0][1]}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparing train and validation image paths","metadata":{}},{"cell_type":"code","source":"train_f = '../input/rsna-pneumonia-detection-challenge/stage_2_train_images'\ntest_f = '../input/rsna-pneumonia-detection-challenge/stage_2_test_images'\n\ntrain_paths = [os.path.join(train_f, image[0]) for image in train_labels]\nval_paths = [os.path.join(train_f, image[0]) for image in val_labels]\n\nprint(len(train_paths))\nprint(len(val_paths))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show some samples from data","metadata":{}},{"cell_type":"code","source":"def imshow(num_to_show=9):\n    \n    plt.figure(figsize=(10,10))\n    \n    for i in range(num_to_show):\n        plt.subplot(3, 3, i+1)\n        plt.grid(False)\n        plt.xticks([])\n        plt.yticks([])\n        \n        img_dcm = dcmread(f'{train_paths[i+20]}.dcm')\n        img_np = img_dcm.pixel_array\n        plt.imshow(img_np, cmap=plt.cm.binary)\n        plt.xlabel(train_labels[i+20][1])\n\nimshow()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Composing transformations","metadata":{}},{"cell_type":"code","source":"transform = transforms.Compose([\n    transforms.RandomHorizontalFlip(),\n    transforms.Resize(224),\n    transforms.ToTensor()])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Write a custom dataset ","metadata":{}},{"cell_type":"code","source":"class Dataset(data.Dataset):\n    \n    def __init__(self, paths, labels, transform=None):\n        self.paths = paths\n        self.labels = labels\n        self.transform = transform\n    \n    def __getitem__(self, index):\n        image = dcmread(f'{self.paths[index]}.dcm')\n        image = image.pixel_array\n        image = image / 255.0\n\n        image = (255*image).clip(0, 255).astype(np.uint8)\n        image = Image.fromarray(image).convert('RGB')\n\n        label = self.labels[index][1]\n        \n        if self.transform is not None:\n            image = self.transform(image)\n            \n        return image, label\n    \n    def __len__(self):\n        \n        return len(self.paths)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check the custom dataset","metadata":{}},{"cell_type":"code","source":"train_dataset = Dataset(train_paths, train_labels, transform=transform)\nimage = iter(train_dataset)\nimg, label = next(image)\nprint(f'Tensor:{img}, Label:{label}')\nimg = np.transpose(img, (1, 2, 0))\nplt.imshow(img)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train image shape","metadata":{}},{"cell_type":"code","source":"img.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare training and validation dataloader","metadata":{}},{"cell_type":"code","source":"train_dataset = Dataset(train_paths, train_labels, transform=transform)\nval_dataset = Dataset(val_paths, val_labels, transform=transform)\ntrain_loader = data.DataLoader(dataset=train_dataset, batch_size=128, shuffle=True)\nval_loader = data.DataLoader(dataset=val_dataset, batch_size=128, shuffle=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check dataloader","metadata":{}},{"cell_type":"code","source":"batch = iter(train_loader)\nimages, labels = next(batch)\n\nimage_grid = torchvision.utils.make_grid(images[:4])\nimage_np = image_grid.numpy()\nimg = np.transpose(image_np, (1, 2, 0))\nplt.imshow(img)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Specify device object","metadata":{}},{"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Orthogonal CNN Helper (http://pwang.pw/ocnn.html)","metadata":{}},{"cell_type":"code","source":"\"\"\" helper function\noriginal author baiyu\nmodified by Peter Wang (@samaoline)\n\"\"\"\n\nimport sys\n\nimport numpy as np\n\n\ndef conv_orth_dist(kernel, stride = 1):\n    [o_c, i_c, w, h] = kernel.shape\n    assert (w == h),\"Do not support rectangular kernel\"\n    #half = np.floor(w/2)\n    assert stride<w,\"Please use matrix orthgonality instead\"\n    new_s = stride*(w-1) + w#np.int(2*(half+np.floor(half/stride))+1)\n    temp = torch.eye(new_s*new_s*i_c).reshape((new_s*new_s*i_c, i_c, new_s,new_s)).cuda()\n    out = (F.conv2d(temp, kernel, stride=stride)).reshape((new_s*new_s*i_c, -1))\n    Vmat = out[np.floor(new_s**2/2).astype(int)::new_s**2, :]\n    temp= np.zeros((i_c, i_c*new_s**2))\n    for i in range(temp.shape[0]):temp[i,np.floor(new_s**2/2).astype(int)+new_s**2*i]=1\n    return torch.norm( Vmat@torch.t(out) - torch.from_numpy(temp).float().cuda() )\n    \ndef deconv_orth_dist(kernel, stride = 2, padding = 1):\n    [o_c, i_c, w, h] = kernel.shape\n    output = torch.conv2d(kernel, kernel, stride=stride, padding=padding)\n    target = torch.zeros((o_c, o_c, output.shape[-2], output.shape[-1])).cuda()\n    ct = int(np.floor(output.shape[-1]/2))\n    target[:,:,ct,ct] = torch.eye(o_c).cuda()\n    return torch.norm( output - target )\n    \ndef orth_dist(mat, stride=None):\n    mat = mat.reshape( (mat.shape[0], -1) )\n    if mat.shape[0] < mat.shape[1]:\n        mat = mat.permute(1,0)\n    return torch.norm( torch.t(mat)@mat - torch.eye(mat.shape[1]).cuda())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load pre-trained ResNet18 and fine-tune","metadata":{}},{"cell_type":"code","source":"model = torchvision.models.resnet18(pretrained=True)\nnum_ftrs = model.fc.in_features\n# Here the size of each output sample is set to 2.\n# Alternatively, it can be generalized to nn.Linear(num_ftrs, len(class_names)).\nmodel.fc = nn.Linear(num_ftrs, 2)\n\nmodel.to(device)\n\ncriterion = nn.CrossEntropyLoss()\n\n# Observe that all parameters are being optimized\noptimizer = torch.optim.SGD(model.parameters(), lr=0.001, momentum=0.9)\n\n# Decay LR by a factor of 0.1 every 7 epochs\nexp_lr_scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Write a train code and RUN","metadata":{}},{"cell_type":"code","source":"num_epochs = 20\n# Train the model\ntotal_step = len(train_loader)\nfor epoch in range(num_epochs):\n    # Training step\n    for i, (images, labels) in tqdm(enumerate(train_loader)):\n        images = images.to(device)\n        labels = labels.to(device)\n        \n        # Forward pass\n        outputs = model(images)\n        diff = orth_dist(model.layer2[0].downsample[0].weight) + orth_dist(model.layer3[0].downsample[0].weight) + orth_dist(model.layer4[0].downsample[0].weight)\n        diff += deconv_orth_dist(model.layer1[0].conv1.weight, stride=1) + deconv_orth_dist(model.layer1[1].conv1.weight, stride=1)\n      #  diff += deconv_orth_dist(model.layer1[0].conv2.weight, stride=1) + deconv_orth_dist(model.layer1[1].conv2.weight, stride=1)\n        diff += deconv_orth_dist(model.layer2[0].conv1.weight, stride=2) + deconv_orth_dist(model.layer2[1].conv1.weight, stride=1)\n      #  diff += deconv_orth_dist(model.layer2[0].conv2.weight, stride=1) + deconv_orth_dist(model.layer2[1].conv2.weight, stride=1)\n        diff += deconv_orth_dist(model.layer3[0].conv1.weight, stride=2) + deconv_orth_dist(model.layer3[1].conv1.weight, stride=1)\n      #  diff += deconv_orth_dist(model.layer3[0].conv2.weight, stride=2) + deconv_orth_dist(model.layer3[1].conv2.weight, stride=1)\n        diff += deconv_orth_dist(model.layer4[0].conv1.weight, stride=2) + deconv_orth_dist(model.layer4[1].conv1.weight, stride=1)\n      #  diff += deconv_orth_dist(model.layer4[0].conv2.weight, stride=2) + deconv_orth_dist(model.layer4[1].conv2.weight, stride=1)\n        loss = criterion(outputs, labels) + diff * 0.1\n        \n        # Backward and optimize\n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n        \n        if (i+1) % 2000 == 0:\n            \n            print(\"Epoch [{}/{}], Step [{}/{}], Loss: {:.4f}\"\n                   .format(epoch+1, num_epochs, i+1, total_step, loss.item()))\n\n\n    # Validation step\n    correct = 0\n    total = 0  \n    for images, labels in tqdm(val_loader):\n        images = images.to(device)\n        labels = labels.to(device)\n        predictions = model(images)\n        _, predicted = torch.max(predictions, 1)\n        total += labels.size(0)\n        correct += (labels == predicted).sum()\n    print(f'Epoch: {epoch+1}/{num_epochs}, Val_Acc: {100*correct/total}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test model","metadata":{}},{"cell_type":"code","source":"model.eval()\n\ncorrect = 0\ntotal = 0  \nfor images, labels in tqdm(val_loader):\n    images = images.to(device)\n    labels = labels.to(device)\n    predictions = model(images)\n    _, predicted = torch.max(predictions, 1)\n    total += labels.size(0)\n    correct += (labels == predicted).sum()\nprint(f'Val_Acc: {100*correct/total}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(model)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pip install git+git://github.com/raghakot/keras-vis.git --upgrade --no-deps","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(images)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def print_saliency(i):\n    image_test = images[i]\n    # from: https://towardsdatascience.com/saliency-map-using-pytorch-68270fe45e80\n    image_test = image_test.reshape(1, 3, 224, 224)\n    image_test = image_test.to(device)\n    image_test.requires_grad_()\n\n    # Retrieve output from the image\n    output = model(image_test)\n\n    # Catch the output\n    output_idx = output.argmax()\n    output_max = output[0, output_idx]\n\n    # Do backpropagation to get the derivative of the output based on the image\n    output_max.backward()\n    \n    # Retireve the saliency map and also pick the maximum value from channels on each pixel.\n    # In this case, we look at dim=1. Recall the shape (batch_size, channel, width, height)\n    saliency, _ = torch.max(image_test.grad.data.abs(), dim=1) \n    saliency = saliency.reshape(224, 224)\n\n    # Reshape the image\n    image_test = image_test.reshape(-1, 224, 224)\n\n    # Visualize the image and the saliency map\n    fig, ax = plt.subplots(1, 2)\n    ax[0].imshow(image_test.cpu().detach().numpy().transpose(1, 2, 0))\n    ax[0].axis('off')\n    ax[1].imshow(saliency.cpu(), cmap='hot')\n    ax[1].axis('off')\n    plt.tight_layout()\n    fig.suptitle('The Image and Its Saliency Map')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print_saliency(0)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print_saliency(50)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print_saliency(78)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"torch.save(model.state_dict(), \"orthogonal_cnn.pt\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}