{"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 (0.453 after reset of the LB). 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":{"_kg_hide-input":false,"trusted":true,"_uuid":"8ead08116baa59092b769e8f05de7426ef7c0faa","_kg_hide-output":true},"cell_type":"code","source":"!pip install fastai==0.7.0 --no-deps\n!pip install torch==0.4.1 torchvision==0.2.1","execution_count":null,"outputs":[]},{"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\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1dcb0ed47f29e7a457cc54500b97ab4e96405e8e"},"cell_type":"code","source":"PATH = './'\nTRAIN = '../input/human-protein-atlas-image-classification/train/'\nTEST = '../input/human-protein-atlas-image-classification/test/'\nLABELS = '../input/human-protein-atlas-image-classification/train.csv'\nSPLIT = '../input/protein-trainval-split/'\nnw = 2   #number of workers for data loader","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":{"_uuid":"6053d1ac2002b769636d1857eda3eb12dd43f06e"},"cell_type":"markdown","source":"### Data"},{"metadata":{"trusted":true,"_uuid":"d2376dd252c2cd9a88fd60866e0ac4bec8ea1314"},"cell_type":"code","source":"#using a split that includes all classes in val\nwith open(os.path.join(SPLIT,'tr_names.txt'), 'r') as text_file:\n    tr_n = text_file.read().split(',')\nwith open(os.path.join(SPLIT,'val_names.txt'), 'r') as text_file:\n    val_n = text_file.read().split(',')\ntest_names = sorted({f[:36] for f in os.listdir(TEST)})\nprint(len(tr_n),len(val_n))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e14c49dde171f5ffb21cbe67a7d4bb68a88e5084"},"cell_type":"code","source":"#creating duplicates for rare classes in train set\nclass Oversampling:\n    def __init__(self,path):\n        self.train_labels = pd.read_csv(path).set_index('Id')\n        self.train_labels['Target'] = [[int(i) for i in s.split()] \n                                       for s in self.train_labels['Target']]  \n        #set the minimum number of duplicates for each class\n        self.multi = [1,1,1,1,1,1,1,1,\n                      4,4,4,1,1,1,1,4,\n                      1,1,1,1,2,1,1,1,\n                      1,1,1,4]\n\n    def get(self,image_id):\n        labels = self.train_labels.loc[image_id,'Target'] if image_id \\\n          in self.train_labels.index else []\n        m = 1\n        for l in labels:\n            if m < self.multi[l]: m = self.multi[l]\n        return m\n    \ns = Oversampling(os.path.join(PATH,LABELS))\ntr_n = [idx for idx in tr_n for _ in range(s.get(idx))]\nprint(len(tr_n),flush=True)","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        return open_rgby(self.path,self.fnames[i])\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,is_test=False):\n    #data augmentation\n    if is_test:\n        aug_tfms = [RandomRotate(30, tfm_y=TfmType.NO),\n                RandomDihedral(tfm_y=TfmType.NO)]\n    else:\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                Cutout(n_holes=25, length=10*sz//128, 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 = 32\nsz = 256\nmd = get_data(sz,bs,is_test=True)\n\nx,y = next(iter(md.aug_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","scrolled":true},"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)\n#x2_tot = np.zeros(4)\n#for x,y in iter(md.trn_dl):\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\n#channel_avr = x_tot/len(md.trn_dl)\n#channel_std = np.sqrt(x2_tot/len(md.trn_dl) - channel_avr**2)\n#channel_avr,channel_std","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"_kg_hide-output":false,"_uuid":"388ac556ebe644259b7bde0c85f8c2dcd05af655"},"cell_type":"markdown","source":"### Model"},{"metadata":{"_uuid":"890a3d017780d3baf8d09545131a77312be4bf8c"},"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. Therefore, I replace the first convolution layer from 7x7 3->64 to 7x7 **4->64** while keeping weights from 3->64. As suggested by @[Skylum] OleksandrSavsunenko in comments, the weights for the last channel can be initialized with values from other channels (the model requires substantial retraining in any case). Also, as @Artyom Palvelev pinted out, the weights should not be set to zero."},{"metadata":{"trusted":true,"_uuid":"4f82ac4cd3a7aaf5b61a5ffa4e882ee0dcec125c"},"cell_type":"code","source":"class Resnet34_4(nn.Module):\n    def __init__(self, pre=True):\n        super().__init__()\n        encoder = resnet34(pretrained=pre)\n        \n        self.conv1 = nn.Conv2d(4, 64, kernel_size=7, stride=2, padding=3, bias=False)\n        if(pre):\n            w = encoder.conv1.weight\n            self.conv1.weight = nn.Parameter(torch.cat((w,\n                                    0.5*(w[:,:1,:,:]+w[:,2:,:,:])),dim=1))\n        \n        self.bn1 = encoder.bn1\n        self.relu = nn.ReLU(inplace=True) \n        self.maxpool = nn.MaxPool2d(kernel_size=3, stride=2, padding=1)\n        self.layer0 = nn.Sequential(self.conv1,self.relu,self.bn1,self.maxpool)\n        self.layer1 = encoder.layer1\n        self.layer2 = encoder.layer2\n        self.layer3 = encoder.layer3\n        self.layer4 = encoder.layer4\n        #the head will be added automatically by fast.ai\n        \n    def forward(self, x):\n        x = self.layer0(x)\n        x = self.layer1(x)\n        x = self.layer2(x)\n        x = self.layer3(x)\n        x = self.layer4(x)\n        \n        return x","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":"a68caa5f46ccb2ca2fca58da2edc218352f1ed0b"},"cell_type":"markdown","source":"Class for accumulation the statistics on each batch and culculation F1 score for entire val dataset."},{"metadata":{"trusted":true,"_uuid":"30197d3107d0dea338fb88c1b5c50ec356f9e17d"},"cell_type":"code","source":"class F1:\n    __name__ = 'F1 macro'\n    def __init__(self,n=28):\n        self.n = n\n        self.TP = np.zeros(self.n)\n        self.FP = np.zeros(self.n)\n        self.FN = np.zeros(self.n)\n\n    def __call__(self,preds,targs,th=0.0):\n        preds = (preds > th).int()\n        targs = targs.int()\n        self.TP += (preds*targs).float().sum(dim=0)\n        self.FP += (preds > targs).float().sum(dim=0)\n        self.FN += (preds < targs).float().sum(dim=0)\n        score = (2.0*self.TP/(2.0*self.TP + self.FP + self.FN + 1e-6)).mean()\n        return score\n\n    def reset(self):\n        #macro F1 score\n        score = (2.0*self.TP/(2.0*self.TP + self.FP + self.FN + 1e-6))\n        print('F1 macro:',score.mean(),flush=True)\n        #print('F1:',score)\n        self.TP = np.zeros(self.n)\n        self.FP = np.zeros(self.n)\n        self.FN = np.zeros(self.n)\n\nclass F1_callback(Callback):\n    def __init__(self, n=28):\n        self.f1 = F1(n)\n\n    def on_epoch_end(self, metrics):\n        self.f1.reset()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"584f2534a3139c6210a62f1bad6b0e7197c50e06"},"cell_type":"markdown","source":"### Training"},{"metadata":{"trusted":true,"_uuid":"17c6374e2043d19b460ae425eb4009855423bc27","scrolled":true,"_kg_hide-output":true},"cell_type":"code","source":"sz = 256 #image size\nbs = 100  #batch size\n\nmd = get_data(sz,bs)\nlearner = ConvLearner.pretrained(Resnet34_4, md, ps=0.7) #dropout 50%\nlearner.opt_fn = optim.Adam\nlearner.clip = 1.0 #gradient clipping\nlearner.crit = FocalLoss()\nf1_callback = F1_callback()\nlearner.metrics = [acc,f1_callback.f1]\nlearner.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.05, and the recommended learning rate is ~0.005."},{"metadata":{"trusted":true,"_uuid":"61cb0f0ceb2141d4bd3d2f2b181a207f7b503dd2"},"cell_type":"code","source":"with warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    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. **When you check F1 macro metric, look at \"F1 macro: XXX\"** before the epoch number, which shows the result calculated for entire val dataset."},{"metadata":{"trusted":true,"_uuid":"fa183a7b7809e4cbd6aa16916d9915a24826d58b"},"cell_type":"code","source":"lr = 0.5e-3\nwith warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    learner.fit(lr,1,callbacks=[f1_callback])","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":"with warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    learner.fit(lrs/4,4,cycle_len=2,use_clr=(10,20),callbacks=[f1_callback])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f36c0b2468425ee0f064decbc1ee40777c84eb52"},"cell_type":"code","source":"with warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    learner.fit(lrs/16,2,cycle_len=4,use_clr=(10,20),callbacks=[f1_callback])","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":"d20a5444b01a8ab1773ce0beb7e253f90d1ffe72"},"cell_type":"code","source":"with warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    learner.fit(lrs/32,1,cycle_len=8,use_clr=(10,20),callbacks=[f1_callback])","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":"md = get_data(sz,bs,is_test=True)\nlearner.set_data(md)\npreds,y = learner.TTA(n_aug=8)\npreds = np.stack(preds, axis=-1)\npred = preds.mean(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":"c093ae5b146cd2803a74d8fd42580f3913bb3f99"},"cell_type":"code","source":"def sigmoid_np(x):\n    return 1.0/(1.0 + np.exp(-x))\n\ndef F1_soft(preds,targs,th=0.0,d=25.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 = np.zeros(len(name_label_dict))\n    wd = 1e-5\n    error = lambda p: np.concatenate((F1_soft(x,y,p) - 1.0,\n                                      wd*p), axis=None)\n    p, success = opt.leastsq(error, params)\n    return p","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9ef8d61efd6e02b8f54798188cd48a595a0cf41f"},"cell_type":"code","source":"th = fit_val(pred,y)\nprint('Thresholds: ',th)\nprint('F1 macro: ',f1_score(y, pred>th, average='macro'))\nprint('F1 macro (th = 0.0): ',f1_score(y, pred>0.0, average='macro'))\nprint('F1 micro: ',f1_score(y, pred>th, average='micro'))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"31376a63bc7dfd7111f80f99dec92f85625fc377"},"cell_type":"markdown","source":"Using CV to prevent overfitting the thresholds:"},{"metadata":{"trusted":true,"_uuid":"516debe8dafc1fbdc6056b21d31a7c947c423c71"},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nth, score, cv = 0,0,10\nfor i in range(cv):\n    xt,xv,yt,yv = train_test_split(pred,y,test_size=0.5,random_state=i)\n    th_i = fit_val(xt,yt)\n    th += th_i\n    score += f1_score(yv, xv>th_i, average='macro')\nth/=cv\nscore/=cv\nprint('Thresholds: ',th)\nprint('F1 macro avr:',score)\nprint('F1 macro: ',f1_score(y, pred>th, 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 > 0.5).mean(axis=0))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9038a911caa06dab9e363484d3b5db3f828d83ab"},"cell_type":"markdown","source":"Plot the distribution of lables to see how the model performs. Be aware about **log scale**."},{"metadata":{"trusted":true,"_uuid":"2919d16f621ec4ab661e4e91ecd64fabf7637b53","scrolled":true},"cell_type":"code","source":"f1 = f1_score(y, pred>th, average=None)\nfor i in range(len(name_label_dict)):\n    bins = np.linspace(pred[:,i].min(), pred[:,i].max(), 50)\n    plt.hist(pred[y[:,i] == 0][:,i], bins, alpha=0.5, log=True, label='false')\n    plt.hist(pred[y[:,i] == 1][:,i], bins, alpha=0.5, log=True, label='true')\n    plt.legend(loc='upper right')\n    print(name_label_dict[i],i, f1[i], th[i])\n    plt.axvline(x=th[i], color='k', linestyle='--')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"83132c40a82e38beb33ba62fa19e5a0233a9cbbb"},"cell_type":"markdown","source":"### Submission"},{"metadata":{"trusted":true,"_uuid":"cedaa4bf73c9b3381ac3e6da32b040e05ce24b91","scrolled":true},"cell_type":"code","source":"preds_t,y_t = learner.TTA(n_aug=8,is_test=True)\npreds_t = np.stack(preds_t, axis=-1)\npred_t = preds_t.mean(axis=-1)","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.0, 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    df = pd.DataFrame({'Id':learner.data.test_ds.fnames,'Predicted':pred_list})\n    df.sort_values(by='Id').to_csv(fname, header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aa26b302d07d0c19de2e905cbd475485aac4dc7c"},"cell_type":"markdown","source":"Save predictions for thresholds calculated based on the validation set and constant value 0.0:"},{"metadata":{"trusted":true,"_uuid":"44a6cae77d5253416dc3db682c60174977b705cd"},"cell_type":"code","source":"save_pred(pred_t,th,'protein_classification_v.csv')\nsave_pred(pred_t,0.0,'protein_classification_0.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2403e30c65590c5cffa4a25a926e0bb119f74544"},"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 (I replaced 0 by 0.01 since there may be a rounding error leading to 0). Automatic fitting the thresholds based on the public LB statistics:\n"},{"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.0,d=50.0):\n    preds = sigmoid_np(d*(preds - th))\n    return preds.mean(axis=0)\n\ndef fit_test(x,y):\n    params = np.zeros(len(name_label_dict))\n    wd = 1e-5\n    error = lambda p: np.concatenate((Count_soft(x,p) - y,\n                                      wd*p), 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)\nprint('Thresholds: ',th_t)\nprint('Fractions: ',(pred_t > th_t).mean(axis=0))\nprint('Fractions (th = 0.0): ',(pred_t > 0.0).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":"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":"76c2edf326d6b2e6c221dd1d1a6c3f4714850d56"},"cell_type":"markdown","source":"Try fitting thresholds based on the frequency of classes in the train dataset:"},{"metadata":{"trusted":true,"_uuid":"7a65a65b0cb44c4458dbafb45cb6aba78c2aba8f"},"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":"704580dfbfd04d451450dfee336187e52a027206"},"cell_type":"code","source":"th_t = fit_test(pred_t,label_fraction)\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":{"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}