{"cells":[{"metadata":{},"cell_type":"markdown","source":"# RXRX.AI\n\n> disclosure: some lines of code were taken from other kernels available on the Kaggle's competition link https://www.kaggle.com/c/recursion-cellular-image-classification/kernels . Authors will be cited as much as possible where appropriate.   \n\n> disclosure 2: this notebook assumes GPU accessibility."},{"metadata":{"trusted":true},"cell_type":"code","source":"!git clone https://github.com/recursionpharma/rxrx1-utils.git && mv rxrx1-utils rxrxutils","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\nimport os\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport rxrxutils.rxrx.io as rio\nfrom scipy import misc\n\nfrom PIL import Image\n\nimport torch\nimport torch.nn as nn\nimport torch.utils.data as D\nfrom torch.optim.lr_scheduler import ExponentialLR\nimport torch.nn.functional as F\n\nfrom torchvision import models, transforms\n\nfrom ignite.engine import Events, create_supervised_evaluator, create_supervised_trainer\nfrom ignite.metrics import Loss, Accuracy\nfrom ignite.contrib.handlers.tqdm_logger import ProgressBar\nfrom ignite.handlers import  EarlyStopping, ModelCheckpoint\n\nfrom tqdm import tqdm_notebook\n\nfrom sklearn.model_selection import train_test_split\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Data folder overview\n!ls -1 ../input","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading a site and visualizing individual channels\n\nThis exploration is inspired by the competition's [creator notebook](https://colab.research.google.com/github/recursionpharma/rxrx1-utils/blob/master/notebooks/visualization.ipynb) and [utils](https://github.com/recursionpharma/rxrx1-utils).   \n  \nThe input for our model will be a 512x512x6 image tensor representing a site, so we will make sure the utilities provided load the site as a tensor with the proper shape. Here, we request the image in experiment RPE-05 on plate 3 in well D19 at site 2."},{"metadata":{"trusted":true},"cell_type":"code","source":"t = rio.load_site('train', 'RPE-05', 3, 'D19', 2, base_path=\"../input\")\nprint(t.shape)\nt_tensor_default = transforms.ToTensor()(t)\nprint(t_tensor_default.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This seems to work, now let's visualize individual channels."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(2, 3, figsize=(24, 16))\n\nfor i, ax in enumerate(axes.flatten()):\n  ax.axis('off')\n  ax.set_title('channel {}'.format(i + 1))\n  _ = ax.imshow(t[:, :, i], cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With the utils provided, we can also convert a site to a RGB format with the `convert_tensor_to_rgb` method. It associates an RGB color with each channel, then aggregates the color channels across the six cellular channels.   "},{"metadata":{"trusted":true},"cell_type":"code","source":"x = rio.convert_tensor_to_rgb(t)\nprint(x.shape)\n\n# plot RGB Image\nplt.figure(figsize=(8, 8))\nplt.axis('off')\n_ = plt.imshow(x)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The utils also include a wrapper `load_site_as_rgb` combining the last two functions."},{"metadata":{"trusted":true},"cell_type":"code","source":"y = rio.load_site_as_rgb('train', 'HUVEC-07', 4, 'K09', 1)\n\nplt.figure(figsize=(8, 8))\nplt.axis('off')\n\n_ = plt.imshow(y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# convert to Tensor\ny_tensor = transforms.ToTensor()(y)\nprint(y_tensor.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Using 6-channel for the base model\n\nFor training the base model, we will use the 6-channel method `rio.load_site` to convert images from a site to a 6x512x512 pytorch Tensor, because it will take less time to load and transform images. Also, it may be possible that opting for a 6-channel approach prevents losing information about the site. "},{"metadata":{},"cell_type":"markdown","source":"#### Metadata\n\nThe metadata for RxRx1 during the Kaggle competition is broken up into four files: train.csv, train_controls.csv, test.csv and test_controls.csv. It is often more convenient to view all the metadata at once, so we have provided a helper function called combine_metadata for doing just that."},{"metadata":{"trusted":true},"cell_type":"code","source":"md = rio.combine_metadata()\nmd.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading images & base model training\n\nThis part is dedicated to create a base model for comparison with full-blown solutions. The flow and some ideas were drawn from [Michael Diskin kernel](https://www.kaggle.com/yhn112/resnet18-baseline-pytorch-ignite)."},{"metadata":{"trusted":true},"cell_type":"code","source":"path_data = '../input'\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(device)\ntorch.manual_seed(0)\nclasses = 1108\nbatch_size = 32","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Implement dataset class & loaders\n\nAs seen in the _Loading a site and visualizing individual channels_ , useful methods were provided to easily load images to pass to our model later on. We will use `rio.load_site` to load site images and convert them to 6x512x512 tensors in our Pytorch image dataset class (ImageDS)."},{"metadata":{"trusted":true},"cell_type":"code","source":"class ImagesDS(D.Dataset):\n    transform = transforms.Compose([\n        transforms.ToTensor(),\n        transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))\n    ])\n\n    def __init__(self, df, mode='train', site=1, channels=[1,2,3,4,5,6]):\n        self.records = df.to_records(index=False)\n        self.channels = channels\n        self.site = site\n        self.mode = mode\n        self.len = df.shape[0]\n        self.first = None\n        \n    def _get_img(self, index):\n        record = self.records[index]\n        return transforms.ToTensor()(rio.load_site(self.mode, record.experiment, record.plate, record.well, self.site, base_path=path_data)).float().cuda()\n        \n    def __getitem__(self, index):\n        img = self._get_img(index)\n        if self.mode == 'train':\n            return img, int(self.records[index].sirna)\n        else:\n            return img, self.records[index].id_code\n\n    def __len__(self):\n        return self.len","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# dataframes for training, cross-validation, and testing\ndf = pd.read_csv(path_data+'/train.csv')\ndf_train, df_val = train_test_split(df, test_size = 0.025, random_state=42)\ndf_test = pd.read_csv(path_data+'/test.csv')\n\n# pytorch training dataset & loader\nds = ImagesDS(df_train, mode='train')\nloader = D.DataLoader(ds, batch_size=batch_size, shuffle=True, num_workers=0)\n\n# pytorch cross-validation dataset & loader\nds_val = ImagesDS(df_val, mode='train')\nval_loader = D.DataLoader(ds_val, batch_size=batch_size, shuffle=True, num_workers=0)\n\n# pytorch test dataset & loader\nds_test = ImagesDS(df_test, mode='test')\ntloader = D.DataLoader(ds_test, batch_size=batch_size, shuffle=False, num_workers=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Prepare base model   \n\nFor our base model, we will use a simple model architecture that will output 1108 classes, since our goal is to classify images into 1108 different siRNA modifications. The simple model architecture was inspired by [Pytorch's tutorial](https://pytorch.org/tutorials/beginner/blitz/cifar10_tutorial.html) and  adapted for the problem at hand."},{"metadata":{"trusted":true},"cell_type":"code","source":"class Net(nn.Module):\n    def __init__(self):\n        super(Net, self).__init__()\n        self.pool = nn.MaxPool2d(2, 2)\n        self.conv1 = nn.Conv2d(6, 16, kernel_size=7, stride=2, padding=3, bias=False)\n        self.conv2 = nn.Conv2d(16, 64, kernel_size=7, stride=2, padding=3, bias=False)\n        self.fc1 = nn.Linear(64 * 32 * 32, 2216, bias=False)\n        self.fc2 = nn.Linear(2216, 1108, bias=False)\n\n    def forward(self, x):\n        x = self.pool(F.relu(self.conv1(x)))\n        x = self.pool(F.relu(self.conv2(x)))\n        x = x.view(-1, 64 * 32 * 32)\n        x = F.relu(self.fc1(x))\n        x = self.fc2(x)\n        return x\n\nmodel = Net()\nmodel.cuda()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"criterion = nn.CrossEntropyLoss()\ncriterion = criterion.cuda()\noptimizer = torch.optim.Adam(model.parameters(), lr=0.0003)\n\nmetrics = {\n    'loss': Loss(criterion),\n    'accuracy': Accuracy(),\n}\n\ntrainer = create_supervised_trainer(model, optimizer, criterion, device=device)\nval_evaluator = create_supervised_evaluator(model, metrics=metrics, device=device)\n\n@trainer.on(Events.EPOCH_COMPLETED)\ndef compute_and_display_val_metrics(engine):\n    epoch = engine.state.epoch\n    metrics = val_evaluator.run(val_loader).metrics\n    print(\"Validation Results - Epoch: {}  Average Loss: {:.4f} | Accuracy: {:.4f} \"\n          .format(engine.state.epoch, \n                      metrics['loss'], \n                      metrics['accuracy']))\n\ncheckpoints = ModelCheckpoint('models', 'Model', save_interval=1, n_saved=2, create_dir=True)\ntrainer.add_event_handler(Events.EPOCH_COMPLETED, checkpoints, {'BaseModel': model})\n\npbar = ProgressBar(bar_format='')\npbar.attach(trainer, output_transform=lambda x: {'loss': x})\n\ntrainer.run(loader, max_epochs=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.eval()\nwith torch.no_grad():\n    preds = np.empty(0)\n    for x, _ in tqdm_notebook(tloader): \n        x = x.to(device)\n        output = model(x)\n        idx = output.max(dim=-1)[1].cpu().numpy()\n        preds = np.append(preds, idx, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(path_data + '/test.csv')\nsubmission['sirna'] = preds.astype(int)\nsubmission.to_csv('submission.csv', index=False, columns=['id_code','sirna'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a href=\"submission.csv\">Download submission file for Base Model</a>"},{"metadata":{},"cell_type":"markdown","source":"## Conclusion\n\nThis gives us a cross-validation score of `0.0011` (`.1%` accuracy), and a test score of 0.002 (`.2%` accuracy). This score is a bit better than chance since we have 1108 classes. An accuracy reflecting chance would be 1/1108, which is equivalent to ~0.09% accuracy. We will explore how we can improve on this score in a next kernel."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.5"}},"nbformat":4,"nbformat_minor":1}