{"cells":[{"metadata":{"_uuid":"d417927fe6d49b8bf27ecd90669bc6ab400a6154"},"cell_type":"markdown","source":"### Overview\nThe goal of this competition is classification of mixed protein patterns. However, unlike most image labeling tasks, where binary or multiclass labeling is considered, in this competition each image can have multiple labels. Multiclass multilabel task has its own specific affecting the design of the model and the loss function. Moreover, the classified images are quite different from ImageNet; therefore, despite usage of a pretrained model is quite helpful, a substantial retraining of entire model is needed. An additional challenge is 4-chanel input to the model (RGBY), which is different from ones used in most of pretrained models (RGB input).\n\nIn this kernel I will show how to handle the above challenges and get started with this competition. I will begin with using a light ResNet34 model and low-resolution images to have a baseline that can be used later to select higher end models and explore the effect of image resolution on the prediction accuracy. **The validation F1 score of the model is ~0.65-0.7**, and I was able to get 0.460 public LB score in V11 of the kernel. Though reuslts are slightly different from one run to another because F1 macro metric is unstable, and sevral items of rear classes contribute in the same way as as thousands items of common classes, 1/28.\n\nThe problem of low public LB score of the model, which mentioned in the first versions of the kernel, is resulted by a bug in the evaluation metric (https://www.kaggle.com/c/human-protein-atlas-image-classification/discussion/69366#409041) that relies on the order of records in the submission file rather than IDs."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"from fastai.conv_learner import *\nfrom fastai.dataset import *\n\nimport pandas as pd\nimport numpy as np\nimport os\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import f1_score\nimport scipy.optimize as opt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1dcb0ed47f29e7a457cc54500b97ab4e96405e8e"},"cell_type":"code","source":"PATH = './'\nTRAIN = '../input/train/'\nTEST = '../input/test/'\nLABELS = '../input/train.csv'\nSAMPLE = '../input/sample_submission.csv'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"87ad005d93033f87ab77acf72d2d6f6cb4a33193"},"cell_type":"code","source":"name_label_dict = {\n0:  'Nucleoplasm',\n1:  'Nuclear membrane',\n2:  'Nucleoli',   \n3:  'Nucleoli fibrillar center',\n4:  'Nuclear speckles',\n5:  'Nuclear bodies',\n6:  'Endoplasmic reticulum',   \n7:  'Golgi apparatus',\n8:  'Peroxisomes',\n9:  'Endosomes',\n10:  'Lysosomes',\n11:  'Intermediate filaments',\n12:  'Actin filaments',\n13:  'Focal adhesion sites',   \n14:  'Microtubules',\n15:  'Microtubule ends',  \n16:  'Cytokinetic bridge',   \n17:  'Mitotic spindle',\n18:  'Microtubule organizing center',  \n19:  'Centrosome',\n20:  'Lipid droplets',\n21:  'Plasma membrane',   \n22:  'Cell junctions', \n23:  'Mitochondria',\n24:  'Aggresome',\n25:  'Cytosol',\n26:  'Cytoplasmic bodies',   \n27:  'Rods & rings' }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8147c0b2f23cf4b646e83f3026698f6df243e468"},"cell_type":"code","source":"nw = 2   #number of workers for data loader\narch = resnet34 #specify target architecture","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6053d1ac2002b769636d1857eda3eb12dd43f06e"},"cell_type":"markdown","source":"### Data"},{"metadata":{"trusted":true,"_uuid":"13dd9234ce2f843d5d7169d5ebef515543a632bd"},"cell_type":"code","source":"train_names = list({f[:36] for f in os.listdir(TRAIN)})\ntest_names = list({f[:36] for f in os.listdir(TEST)})\ntr_n, val_n = train_test_split(train_names, test_size=0.1, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b5ae50160b85e023518ee17d62e4ff216eb5b9fb"},"cell_type":"code","source":"def open_rgby(path,id): #a function that reads RGBY image\n    colors = ['red','green','blue','yellow']\n    flags = cv2.IMREAD_GRAYSCALE\n    img = [cv2.imread(os.path.join(path, id+'_'+color+'.png'), flags).astype(np.float32)/255\n           for color in colors]\n    return np.stack(img, axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"09ae7643632127f1097be17609315b967ebbdc0e"},"cell_type":"markdown","source":"Since a multiclass multilabel task is considered, there are several things about the model that should be pointed out. First, the SOFTMAX MUST NOT BE USED as an output layer because it encourages a single label prediction. The common output function for multilabel tasks is sigmoid. However, combining the sigmoid with the loss function (like in BCE with logits loss or in Focal loss used in this kernel) allows log(sigmoid) optimization of the numerical stability of the loss function. Therefore, sigmoid is also removed."},{"metadata":{"trusted":true,"_uuid":"0dd72024e421ac2ced8a21389e15e914fae8d6a8"},"cell_type":"code","source":"class pdFilesDataset(FilesDataset):\n    def __init__(self, fnames, path, transform):\n        self.labels = pd.read_csv(LABELS).set_index('Id')\n        self.labels['Target'] = [[int(i) for i in s.split()] for s in self.labels['Target']]\n        super().__init__(fnames, transform, path)\n    \n    def get_x(self, i):\n        img = open_rgby(self.path,self.fnames[i])\n        if self.sz == 512: return img \n        else: return cv2.resize(img, (self.sz, self.sz),cv2.INTER_AREA)\n    \n    def get_y(self, i):\n        if(self.path == TEST): return np.zeros(len(name_label_dict),dtype=np.int)\n        else:\n            labels = self.labels.loc[self.fnames[i]]['Target']\n            return np.eye(len(name_label_dict),dtype=np.float)[labels].sum(axis=0)\n        \n    @property\n    def is_multi(self): return True\n    @property\n    def is_reg(self):return True\n    #this flag is set to remove the output sigmoid that allows log(sigmoid) optimization\n    #of the numerical stability of the loss function\n    \n    def get_c(self): return len(name_label_dict) #number of classes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2279cbaea32d80e9bf030b4a3e74ade811a52b3a"},"cell_type":"code","source":"def get_data(sz,bs):\n    #data augmentation\n    aug_tfms = [RandomRotate(30, tfm_y=TfmType.NO),\n                RandomDihedral(tfm_y=TfmType.NO),\n                RandomLighting(0.05, 0.05, tfm_y=TfmType.NO)]\n    #mean and std in of each channel in the train set\n    stats = A([0.08069, 0.05258, 0.05487, 0.08282], [0.13704, 0.10145, 0.15313, 0.13814])\n    tfms = tfms_from_stats(stats, sz, crop_type=CropType.NO, tfm_y=TfmType.NO, \n                aug_tfms=aug_tfms)\n    ds = ImageData.get_ds(pdFilesDataset, (tr_n[:-(len(tr_n)%bs)],TRAIN), \n                (val_n,TRAIN), tfms, test=(test_names,TEST))\n    md = ImageData(PATH, ds, bs, num_workers=nw, classes=None)\n    return md","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"24adecf3a1b558b05026ec43e64b8848ac0b05d7"},"cell_type":"code","source":"bs = 16\nsz = 256\nmd = get_data(sz,bs)\n\nx,y = next(iter(md.trn_dl))\nx.shape, y.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"75d66d8969b52d833f2b3268927f8455d20bed03"},"cell_type":"markdown","source":"Plot several examples of input images."},{"metadata":{"trusted":true,"_uuid":"c288311b2416d46e5ec0126c92ccfed170ded223"},"cell_type":"code","source":"def display_imgs(x):\n    columns = 4\n    bs = x.shape[0]\n    rows = min((bs+3)//4,4)\n    fig=plt.figure(figsize=(columns*4, rows*4))\n    for i in range(rows):\n        for j in range(columns):\n            idx = i+j*columns\n            fig.add_subplot(rows, columns, idx+1)\n            plt.axis('off')\n            plt.imshow((x[idx,:,:,:3]*255).astype(np.int))\n    plt.show()\n    \ndisplay_imgs(np.asarray(md.trn_ds.denorm(x)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"55ef5c338f702c3635080e78fe9b18ca4053c90f"},"cell_type":"markdown","source":"Computing image statistics in the train set. The values listed below are computed without augmentation, therefore the result of the next cell may be a little bit different.\n1. train: (array([0.08069, 0.05258, 0.05487, 0.08282]), array([0.13704, 0.10145, 0.15313, 0.13814]))\n1. test: (array([0.05913, 0.0454 , 0.04066, 0.05928]),  array([0.11734, 0.09503, 0.129  , 0.11528]))"},{"metadata":{"trusted":true,"_uuid":"a907d00a3863cc5f92a9a204cb52188466d363f3"},"cell_type":"code","source":"x_tot = np.zeros(4)\nx2_tot = np.zeros(4)\nfor x,y in iter(md.trn_dl):\n    tmp =  md.trn_ds.denorm(x).reshape(16,-1)\n    x = md.trn_ds.denorm(x).reshape(-1,4)\n    x_tot += x.mean(axis=0)\n    x2_tot += (x**2).mean(axis=0)\n\nchannel_avr = x_tot/len(md.trn_dl)\nchannel_std = np.sqrt(x2_tot/len(md.trn_dl) - channel_avr**2)\nchannel_avr,channel_std","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b947af5caaee24f23249fab49e59bc1532ed1711"},"cell_type":"markdown","source":"### Loss function and metrics"},{"metadata":{"_uuid":"4a8b03adc4418a694339dd4e48f494dc4ec86010"},"cell_type":"markdown","source":"One of the challenges of this competition is strong data imbalance. Some classes, like \"Nucleoplasm\", are very common, while there is a number of rare classes, like \"Endosomes\", \"Lysosomes\", and \"Rods & rings\". In addition, in tasks of multiclass, and especially in multilabel, classification there is always an issue with data imbalance: if you predict 1 class out of 10, given the same number of examples per each class, you have 1 positive vs. 9 negative examples. So, it is crucial to use a loss function that accounts for it. Recently proposed focal loss (https://arxiv.org/pdf/1708.02002.pdf) has revolutionized one stage object localization method in 2017. It is design to address the issue of strong data imbalance, demonstrating amazing results on datasets with imbalance level 1:10-1000. In particular, it works quite well for image segmentation task in \"Airbus Ship Detection Challenge\": https://www.kaggle.com/iafoss/unet34-dice-0-87 .  The implementation of focal loss is borrowed from https://becominghuman.ai/investigating-focal-and-dice-loss-for-the-kaggle-2018-data-science-bowl-65fb9af4f36c ."},{"metadata":{"trusted":true,"_uuid":"b3f8762045c4087ae9aecd22f8c13a075f567d7a"},"cell_type":"code","source":"class FocalLoss(nn.Module):\n    def __init__(self, gamma=2):\n        super().__init__()\n        self.gamma = gamma\n        \n    def forward(self, input, target):\n        if not (target.size() == input.size()):\n            raise ValueError(\"Target size ({}) must be the same as input size ({})\"\n                             .format(target.size(), input.size()))\n\n        max_val = (-input).clamp(min=0)\n        loss = input - input * target + max_val + \\\n            ((-max_val).exp() + (-input - max_val).exp()).log()\n\n        invprobs = F.logsigmoid(-input * (target * 2.0 - 1.0))\n        loss = (invprobs * self.gamma).exp() * loss\n        \n        return loss.sum(dim=1).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"06c8f2af90481b71e1f174c584b7d1d2a0697363"},"cell_type":"code","source":"def acc(preds,targs,th=0.0):\n    preds = (preds > th).int()\n    targs = targs.int()\n    return (preds==targs).float().mean()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"584f2534a3139c6210a62f1bad6b0e7197c50e06"},"cell_type":"markdown","source":"### Training"},{"metadata":{"_uuid":"e42b34e9ac57f99e8c55814f3876f20cb9d93fe6"},"cell_type":"markdown","source":"One of the challenges in this competition is 4-chanel input (RGBY) that limits usage of ImageNet pretrained models taking RGB input. However, the input dataset is too tiny to train even a low capacity model like ResNet34 from scratch. I propose to use the following way to walk around this limitation. The most common way to convert RGBY to RGB is just dropping Y channel while keeping RGB without modification. Therefore, I replace the first convolution layer from 7x7 3->64 to 7x7 **4->64** while keeping weights from 3->64 and setting new initial weighs for Y channel to be zero. It allows using the original weights to initialize the network while giving the opportunity to the model to incorporate Y channel into prediction during the following training (when the first layers of the model are unfreezed). In the following hiden cell I put a code from fast.ai library with adding several lines for the replacment of the first convolutional layer."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"96778d73cae9d6a9a3f1acc2fecd1a7f5363e29b"},"cell_type":"code","source":"class ConvnetBuilder_custom():\n    def __init__(self, f, c, is_multi, is_reg, ps=None, xtra_fc=None, xtra_cut=0, \n                 custom_head=None, pretrained=True):\n        self.f,self.c,self.is_multi,self.is_reg,self.xtra_cut = f,c,is_multi,is_reg,xtra_cut\n        if xtra_fc is None: xtra_fc = [512]\n        if ps is None: ps = [0.25]*len(xtra_fc) + [0.5]\n        self.ps,self.xtra_fc = ps,xtra_fc\n\n        if f in model_meta: cut,self.lr_cut = model_meta[f]\n        else: cut,self.lr_cut = 0,0\n        cut-=xtra_cut\n        layers = cut_model(f(pretrained), cut)\n        \n        #replace first convolutional layer by 4->64 while keeping corresponding weights\n        #and initializing new weights with zeros\n        w = layers[0].weight\n        layers[0] = nn.Conv2d(4,64,kernel_size=(7,7),stride=(2,2),padding=(3, 3), bias=False)\n        layers[0].weight = torch.nn.Parameter(torch.cat((w,torch.zeros(64,1,7,7)),dim=1))\n        \n        self.nf = model_features[f] if f in model_features else (num_features(layers)*2)\n        if not custom_head: layers += [AdaptiveConcatPool2d(), Flatten()]\n        self.top_model = nn.Sequential(*layers)\n\n        n_fc = len(self.xtra_fc)+1\n        if not isinstance(self.ps, list): self.ps = [self.ps]*n_fc\n\n        if custom_head: fc_layers = [custom_head]\n        else: fc_layers = self.get_fc_layers()\n        self.n_fc = len(fc_layers)\n        self.fc_model = to_gpu(nn.Sequential(*fc_layers))\n        if not custom_head: apply_init(self.fc_model, kaiming_normal)\n        self.model = to_gpu(nn.Sequential(*(layers+fc_layers)))\n\n    @property\n    def name(self): return f'{self.f.__name__}_{self.xtra_cut}'\n\n    def create_fc_layer(self, ni, nf, p, actn=None):\n        res=[nn.BatchNorm1d(num_features=ni)]\n        if p: res.append(nn.Dropout(p=p))\n        res.append(nn.Linear(in_features=ni, out_features=nf))\n        if actn: res.append(actn)\n        return res\n\n    def get_fc_layers(self):\n        res=[]\n        ni=self.nf\n        for i,nf in enumerate(self.xtra_fc):\n            res += self.create_fc_layer(ni, nf, p=self.ps[i], actn=nn.ReLU())\n            ni=nf\n        final_actn = nn.Sigmoid() if self.is_multi else nn.LogSoftmax()\n        if self.is_reg: final_actn = None\n        res += self.create_fc_layer(ni, self.c, p=self.ps[-1], actn=final_actn)\n        return res\n\n    def get_layer_groups(self, do_fc=False):\n        if do_fc:\n            return [self.fc_model]\n        idxs = [self.lr_cut]\n        c = children(self.top_model)\n        if len(c)==3: c = children(c[0])+c[1:]\n        lgs = list(split_by_idxs(c,idxs))\n        return lgs+[self.fc_model]\n    \nclass ConvLearner(Learner):\n    def __init__(self, data, models, precompute=False, **kwargs):\n        self.precompute = False\n        super().__init__(data, models, **kwargs)\n        if hasattr(data, 'is_multi') and not data.is_reg and self.metrics is None:\n            self.metrics = [accuracy_thresh(0.5)] if self.data.is_multi else [accuracy]\n        if precompute: self.save_fc1()\n        self.freeze()\n        self.precompute = precompute\n\n    def _get_crit(self, data):\n        if not hasattr(data, 'is_multi'): return super()._get_crit(data)\n\n        return F.l1_loss if data.is_reg else F.binary_cross_entropy if data.is_multi else F.nll_loss\n\n    @classmethod\n    def pretrained(cls, f, data, ps=None, xtra_fc=None, xtra_cut=0, custom_head=None, precompute=False,\n                   pretrained=True, **kwargs):\n        models = ConvnetBuilder_custom(f, data.c, data.is_multi, data.is_reg,\n            ps=ps, xtra_fc=xtra_fc, xtra_cut=xtra_cut, custom_head=custom_head, pretrained=pretrained)\n        return cls(data, models, precompute, **kwargs)\n\n    @classmethod\n    def lsuv_learner(cls, f, data, ps=None, xtra_fc=None, xtra_cut=0, custom_head=None, precompute=False,\n                  needed_std=1.0, std_tol=0.1, max_attempts=10, do_orthonorm=False, **kwargs):\n        models = ConvnetBuilder(f, data.c, data.is_multi, data.is_reg,\n            ps=ps, xtra_fc=xtra_fc, xtra_cut=xtra_cut, custom_head=custom_head, pretrained=False)\n        convlearn=cls(data, models, precompute, **kwargs)\n        convlearn.lsuv_init()\n        return convlearn\n    \n    @property\n    def model(self): return self.models.fc_model if self.precompute else self.models.model\n    \n    def half(self):\n        if self.fp16: return\n        self.fp16 = True\n        if type(self.model) != FP16: self.models.model = FP16(self.model)\n        if not isinstance(self.models.fc_model, FP16): self.models.fc_model = FP16(self.models.fc_model)\n    def float(self):\n        if not self.fp16: return\n        self.fp16 = False\n        if type(self.models.model) == FP16: self.models.model = self.model.module.float()\n        if type(self.models.fc_model) == FP16: self.models.fc_model = self.models.fc_model.module.float()\n\n    @property\n    def data(self): return self.fc_data if self.precompute else self.data_\n\n    def create_empty_bcolz(self, n, name):\n        return bcolz.carray(np.zeros((0,n), np.float32), chunklen=1, mode='w', rootdir=name)\n\n    def set_data(self, data, precompute=False):\n        super().set_data(data)\n        if precompute:\n            self.unfreeze()\n            self.save_fc1()\n            self.freeze()\n            self.precompute = True\n        else:\n            self.freeze()\n\n    def get_layer_groups(self):\n        return self.models.get_layer_groups(self.precompute)\n\n    def summary(self):\n        precompute = self.precompute\n        self.precompute = False\n        res = super().summary()\n        self.precompute = precompute\n        return res\n\n    def get_activations(self, force=False):\n        tmpl = f'_{self.models.name}_{self.data.sz}.bc'\n        # TODO: Somehow check that directory names haven't changed (e.g. added test set)\n        names = [os.path.join(self.tmp_path, p+tmpl) for p in ('x_act', 'x_act_val', 'x_act_test')]\n        if os.path.exists(names[0]) and not force:\n            self.activations = [bcolz.open(p) for p in names]\n        else:\n            self.activations = [self.create_empty_bcolz(self.models.nf,n) for n in names]\n\n    def save_fc1(self):\n        self.get_activations()\n        act, val_act, test_act = self.activations\n        m=self.models.top_model\n        if len(self.activations[0])!=len(self.data.trn_ds):\n            predict_to_bcolz(m, self.data.fix_dl, act)\n        if len(self.activations[1])!=len(self.data.val_ds):\n            predict_to_bcolz(m, self.data.val_dl, val_act)\n        if self.data.test_dl and (len(self.activations[2])!=len(self.data.test_ds)):\n            if self.data.test_dl: predict_to_bcolz(m, self.data.test_dl, test_act)\n\n        self.fc_data = ImageClassifierData.from_arrays(self.data.path,\n                (act, self.data.trn_y), (val_act, self.data.val_y), self.data.bs, classes=self.data.classes,\n                test = test_act if self.data.test_dl else None, num_workers=8)\n\n    def freeze(self):\n        self.freeze_to(-1)\n\n    def unfreeze(self):\n        self.freeze_to(0)\n        self.precompute = False\n\n    def predict_array(self, arr):\n        precompute = self.precompute\n        self.precompute = False\n        pred = super().predict_array(arr)\n        self.precompute = precompute\n        return pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"17c6374e2043d19b460ae425eb4009855423bc27"},"cell_type":"code","source":"sz = 256 #image size\nbs = 64  #batch size\n\nmd = get_data(sz,bs)\nlearner = ConvLearner.pretrained(arch, md, ps=0.5) #dropout 50%\nlearner.opt_fn = optim.Adam\nlearner.clip = 1.0 #gradient clipping\nlearner.crit = FocalLoss()\nlearner.metrics = [acc]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5ecad6cfc5289308523f2d7b946c1926c9b7a1c9","scrolled":true},"cell_type":"code","source":"learner.summary","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b695ea17093722c7dc68dc6fb52b2c3f584c7820"},"cell_type":"markdown","source":"I begin with finding the optimal learning rate. The following function runs training with different lr and records the loss. Increase of the loss indicates onset of divergence of training. The optimal lr lies in the vicinity of the minimum of the curve but before the onset of divergence. Based on the following plot, for the current setup the divergence starts at ~0.1, and the recommended learning rate is ~0.02."},{"metadata":{"trusted":true,"_uuid":"61cb0f0ceb2141d4bd3d2f2b181a207f7b503dd2"},"cell_type":"code","source":"learner.lr_find()\nlearner.sched.plot()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aa1fb7a567976a88370f4aeff9f2333b5ab452e4"},"cell_type":"markdown","source":"First, I train only the head of the model while keeping the rest frozen. It allows to avoid corruption of the pretrained weights at the initial stage of training due to random initialization of the head layers. So the power of transfer learning is fully utilized when the training is continued."},{"metadata":{"trusted":true,"_uuid":"fa183a7b7809e4cbd6aa16916d9915a24826d58b"},"cell_type":"code","source":"lr = 2e-2\nlearner.fit(lr,1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a1bb44fad207cbde1f1459fccfcae54878d856a5"},"cell_type":"markdown","source":"Next, I unfreeze all weights and allow training of entire model. One trick that I use is differential learning rate: the lr of the head part is still lr, while the middle layers of the model are trained with lr/3, and the base is trained with even smaller lr/10. Despite the low-level detectors do not vary much from one image data set to another much, the yellow channel should be trained, and also the images are quite different from ImageNet; therefore, the I decrease the learning rate for first layers only by 10 times. If there was no necessity to train an additional channel and the images were more similar to ImageNet, the learning rates could be [lr/100,lr/10,lr]. Another trick is learning rate annealing. Periodic lr increase followed by slow decrease drives the system out of steep minima (when lr is high) towards broader ones (which are explored when lr decreases) that enhances the ability of the model to generalize and reduces overfitting. The length of the cycles gradually increases during training."},{"metadata":{"trusted":true,"_uuid":"032d2c930d002704e5690f095fc0af043c739982"},"cell_type":"code","source":"learner.unfreeze()\nlrs=np.array([lr/10,lr/3,lr])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7b01fe75583db3453bc3712e4c7d9e36f2e5bdb0"},"cell_type":"code","source":"learner.fit(lrs/4,4,cycle_len=2,use_clr=(10,20))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f36c0b2468425ee0f064decbc1ee40777c84eb52"},"cell_type":"code","source":"learner.fit(lrs/4,2,cycle_len=4,use_clr=(10,20))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3753ac42996bbdff557eddb8778ae5b6e3d613af"},"cell_type":"code","source":"learner.sched.plot_lr()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"775ed9f50c239c18f063ae9d385ad0f951ac68fa"},"cell_type":"code","source":"learner.fit(lrs/16,1,cycle_len=8,use_clr=(5,20))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"39aa0f6326e47f224bfc69dd466e0774ebe6a38d"},"cell_type":"markdown","source":"Save the model for further use or training on higher resolution images."},{"metadata":{"trusted":true,"_uuid":"032a66e6f9db41870dbbbfa5ce7a106363f8b8e0"},"cell_type":"code","source":"learner.save('ResNet34_256_1')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0e19235b01078762646c74bfe41e0862305d8740"},"cell_type":"markdown","source":"### Validation score"},{"metadata":{"_uuid":"a37b5b3730cbf98e41571e3392362cd4115d2e00"},"cell_type":"markdown","source":"Evaluate the score with using TTA (test time augmentation)."},{"metadata":{"trusted":true,"_uuid":"e5f8481c905db44b4132967adec66d70521aa06b"},"cell_type":"code","source":"def sigmoid_np(x):\n    return 1.0/(1.0 + np.exp(-x))\n\npreds,y = learner.TTA(n_aug=16)\npreds = np.stack(preds, axis=-1)\npreds = sigmoid_np(preds)\npred = preds.max(axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1381e2d9055e3d15e8914a4a6b80d5cc0c40f8b3"},"cell_type":"markdown","source":"Instead of 0.5, one can adjust the values of the threshold for each class individually to boost the score. The code below does it automatically."},{"metadata":{"trusted":true,"_uuid":"fb4c482fc261eabd13e367305f7d6583d1e4d45f"},"cell_type":"code","source":"def F1_soft(preds,targs,th=0.5,d=50.0):\n    preds = sigmoid_np(d*(preds - th))\n    targs = targs.astype(np.float)\n    score = 2.0*(preds*targs).sum(axis=0)/((preds+targs).sum(axis=0) + 1e-6)\n    return score\n\ndef fit_val(x,y):\n    params = 0.5*np.ones(len(name_label_dict))\n    wd = 1e-5\n    error = lambda p: np.concatenate((F1_soft(x,y,p) - 1.0,\n                                      wd*(p - 0.5)), axis=None)\n    p, success = opt.leastsq(error, params)\n    return p","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"98cadeb22bcdbf5ddde44a98fa9e90d68249741a"},"cell_type":"code","source":"th = fit_val(pred,y)\nth[th<0.1] = 0.1\nprint('Thresholds: ',th)\nprint('F1 macro: ',f1_score(y, pred>th, average='macro'))\nprint('F1 macro (th = 0.5): ',f1_score(y, pred>0.5, average='macro'))\nprint('F1 micro: ',f1_score(y, pred>th, average='micro'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8ca60fed703d8ce6fe9d524359bfae7f4feb7882"},"cell_type":"code","source":"print('Fractions: ',(pred > th).mean(axis=0))\nprint('Fractions (true): ',(y > th).mean(axis=0))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"83132c40a82e38beb33ba62fa19e5a0233a9cbbb"},"cell_type":"markdown","source":"### Submission"},{"metadata":{"trusted":true,"_uuid":"cedaa4bf73c9b3381ac3e6da32b040e05ce24b91"},"cell_type":"code","source":"preds_t,y_t = learner.TTA(n_aug=16,is_test=True)\npreds_t = np.stack(preds_t, axis=-1)\npreds_t = sigmoid_np(preds_t)\npred_t = preds_t.max(axis=-1) #max works better for F1 macro score","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7229f689dae9d172d3f233d31a8d9f7560256b9f"},"cell_type":"markdown","source":"**It is very important to keep the same order of ids as in the sample submission** https://www.kaggle.com/c/human-protein-atlas-image-classification/discussion/69366#409041 since the competition metric relies only on the order of recods ignoring IDs."},{"metadata":{"trusted":true,"_uuid":"2cd07338fab7fcd9b883e01ea44ec8a55d42d60c"},"cell_type":"code","source":"def save_pred(pred, th=0.5, fname='protein_classification.csv'):\n    pred_list = []\n    for line in pred:\n        s = ' '.join(list([str(i) for i in np.nonzero(line>th)[0]]))\n        pred_list.append(s)\n        \n    sample_df = pd.read_csv(SAMPLE)\n    sample_list = list(sample_df.Id)\n    pred_dic = dict((key, value) for (key, value) \n                in zip(learner.data.test_ds.fnames,pred_list))\n    pred_list_cor = [pred_dic[id] for id in sample_list]\n    df = pd.DataFrame({'Id':sample_list,'Predicted':pred_list_cor})\n    df.to_csv(fname, header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3306ba0f7666d604087e82a484c948d76480aad8"},"cell_type":"markdown","source":"Similar to validation, additional adjustment may be done based on the public LB probing results (https://www.kaggle.com/c/human-protein-atlas-image-classification/discussion/68678) to predict approximately the same fraction of images of a particular class as expected from the public LB."},{"metadata":{"_uuid":"d57ff6d1201bfbfb5aae4be4c10c0ff1cce11faa"},"cell_type":"markdown","source":"Somehow the thresholds that I found manually for one of the models are working the best."},{"metadata":{"trusted":true,"_uuid":"a1c388b048f1d3f613c1d93d2ea8259be25c54ec"},"cell_type":"code","source":"th_t = np.array([0.565,0.39,0.55,0.345,0.33,0.39,0.33,0.45,0.38,0.39,\n               0.34,0.42,0.31,0.38,0.49,0.50,0.38,0.43,0.46,0.40,\n               0.39,0.505,0.37,0.47,0.41,0.545,0.32,0.1])\nprint('Fractions: ',(pred_t > th_t).mean(axis=0))\nsave_pred(pred_t,th_t)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ea13d711304554af07359f162e9875631c35894d"},"cell_type":"markdown","source":"Automatic fitting the thresholds based on the public LB statistics."},{"metadata":{"trusted":true,"_uuid":"546164c9576241139af6de50e5c3b585fcee2bbd"},"cell_type":"code","source":"lb_prob = [\n 0.362397820,0.043841336,0.075268817,0.059322034,0.075268817,\n 0.075268817,0.043841336,0.075268817,0.010000000,0.010000000,\n 0.010000000,0.043841336,0.043841336,0.014198783,0.043841336,\n 0.010000000,0.028806584,0.014198783,0.028806584,0.059322034,\n 0.010000000,0.126126126,0.028806584,0.075268817,0.010000000,\n 0.222493880,0.028806584,0.010000000]\n# I replaced 0 by 0.01 since there may be a rounding error leading to 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"966d9e3a03f689369df5ce3e3c74c3d6bda28c29"},"cell_type":"code","source":"def Count_soft(preds,th=0.5,d=50.0):\n    preds = sigmoid_np(d*(preds - th))\n    return preds.mean(axis=0)\n\ndef fit_test(x,y):\n    params = 0.5*np.ones(len(name_label_dict))\n    wd = 1e-5\n    error = lambda p: np.concatenate((Count_soft(x,p) - y,\n                                      wd*(p - 0.5)), axis=None)\n    p, success = opt.leastsq(error, params)\n    return p","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"25b3c95265890a290fb640ad4b4498e13d81881d"},"cell_type":"code","source":"th_t = fit_test(pred_t,lb_prob)\nth_t[th_t<0.1] = 0.1\nprint('Thresholds: ',th_t)\nprint('Fractions: ',(pred_t > th_t).mean(axis=0))\nprint('Fractions (th = 0.5): ',(pred_t > 0.5).mean(axis=0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e613772527bca8a6aa9f4bd9f7b89575966b2947"},"cell_type":"code","source":"save_pred(pred_t,th_t,'protein_classification_f.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7079a7e1a75aba0f02b17df61b699980108e5cba"},"cell_type":"markdown","source":"Save also predictions for a threshold calculated based on the validation set and constant value 0.5:"},{"metadata":{"trusted":true,"_uuid":"1f822ae26a507eb79986de38234ed10f10f9190f"},"cell_type":"code","source":"save_pred(pred_t,th,'protein_classification_v.csv')\nsave_pred(pred_t,0.5,'protein_classification_05.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ea38f2fdcc64a332e795e0402947a42d52a5b50b"},"cell_type":"markdown","source":"Try using the threshold from validation set for classes not present in the public LB:"},{"metadata":{"trusted":true,"_uuid":"6371d4a8ad55d3b0b0d757ebd42b0ebba7efefc4"},"cell_type":"code","source":"class_list = [8,9,10,15,20,24,27]\nfor i in class_list:\n    th_t[i] = th[i]\nsave_pred(pred_t,th_t,'protein_classification_c.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c5d19b4b9434a443b19a23f7b0839b60a6d3c6af"},"cell_type":"markdown","source":"Try fitting thresholds based on the frequency of classes in the train dataset:"},{"metadata":{"trusted":true,"_uuid":"effda6d101ee0ce6c160da72be7f869df4f09d07"},"cell_type":"code","source":"labels = pd.read_csv(LABELS).set_index('Id')\nlabel_count = np.zeros(len(name_label_dict))\nfor label in labels['Target']:\n    l = [int(i) for i in label.split()]\n    label_count += np.eye(len(name_label_dict))[l].sum(axis=0)\nlabel_fraction = label_count.astype(np.float)/len(labels)\nlabel_count, label_fraction","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"083b552fae17124170dbb7dcc673e89b805f97fb"},"cell_type":"code","source":"th_t = fit_test(pred_t,label_fraction)\nth_t[th_t<0.05] = 0.05\nprint('Thresholds: ',th_t)\nprint('Fractions: ',(pred_t > th_t).mean(axis=0))\nsave_pred(pred_t,th_t,'protein_classification_t.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6b932d216fddcddf9fe953bd492319327728772c"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d2cfeff970696dbce2f6c8b74a2dc1770c6f7a73"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"adf36545450adbc82e6a0a483ac0ccc26c75be3a"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"981df8acdea146109b85eddd599b1f6920b97faf"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1e4b77b81ca029ccb4a33562b3e5a66d25572e3f"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6346b7e02d1102fda17c1f38eda9689b19155d66"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c68069b4f3c841fa11a33ffae87327382160573e"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"82257432429cd11de75c9c989dfb0f3664d224b1"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2b98f402ef8bdd6477541697e53d899c5e7b9eeb"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c6a6591f2712e75e3230f77c4277afe5b6396821"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f2aba6fed3d1767b2fb4efc2fdaeca22926b5e22"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"82a0b70e7bc961b7516b137837f642a0d65caf8d"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7092a3f61503d07bb38b28bfd293bc3146bcf382"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"47e6a7bf79ecb5f3c1bd0893f6580d6c66fa55bd"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"27b16167d1a8efaa601422e46578478468524d7f"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}