{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!cp /kaggle/input/gdcm-conda-install/gdcm.tar .\n!tar  -xzf gdcm.tar\n!conda install -q --offline ./gdcm/gdcm-2.8.9-py37h71b2a6d_0.tar.bz2","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"from fastai.vision.all import *\nfrom fastai.medical.imaging import *","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"datapath = Path(\"/kaggle/input/rsna-str-pulmonary-embolism-detection/\")\ntestdatapath = datapath/'test'\n\ncnn2dmodelpath = Path(\"/kaggle/input/rsnape2dmodels/\")\ncnn3dmodelpath = Path(\"/kaggle/input/rsnape3dmodels/\")\nseqmodelpath = Path(\"/kaggle/input/rsnapeseqmodels/\")\n\ntest_df = pd.read_csv(datapath/'test.csv')\nsub_df = pd.read_csv(datapath/'sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"[o for o in sub_df['id'].values if \"df06fad17bc3\" in o]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"device = default_device()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"RAM Until Here: 1GB"},{"metadata":{},"cell_type":"markdown","source":"### Predict Study"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_study_dirnames = [datapath/'test'/o for o in test_df['StudyInstanceUID'].unique()]\nstudy_dirname = test_study_dirnames[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# RGB windows\nlung_window = (1500, -600)\npe_window = (700, 100)\nmediastinal_window = (400, 40)\nwindows = (lung_window, pe_window, mediastinal_window)\n\ndef read_dcm_img(dcm, windows=windows):\n    \"Read single slice in RGB\"\n    return torch.stack([dcm.windowed(*w) for w in windows])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### CNN Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load CNN model\ndef get_dls(tensors, size=256, bs=128):\n    \"Get study dataloader\"\n    tfms = [[RandomResizedCropGPU(size, min_scale=0.9)], []]\n\n    dsets = Datasets(tensors, tfms=tfms, splits=([0,1], [2,3]))\n\n#     batch_tfms = [Normalize.from_stats(*imagenet_stats)]\n    batch_tfms = []\n    dls = dsets.dataloaders(bs=bs, after_batch=batch_tfms, num_workers=2)\n    return dls\n\ndls = get_dls(torch.zeros(4, 3, 512, 512), size=480, bs=32)\ndls.c = 2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learn2d_fold0 = cnn_learner(dls, xresnet34, pretrained=False, loss_func=nn.CrossEntropyLoss(), model_dir='.')\nlearn2d_fold0.path = cnn2dmodelpath\nlearn2d_fold0.load('xresnet34-512-PR-fold0');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learn2d_fold1 = cnn_learner(dls, xresnet34, pretrained=False, loss_func=nn.CrossEntropyLoss(), model_dir='.')\nlearn2d_fold1.path = cnn2dmodelpath\nlearn2d_fold1.load('xresnet34-512-PR-fold1');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learn2d_fold2 = cnn_learner(dls, xresnet34, pretrained=False, loss_func=nn.CrossEntropyLoss(), model_dir='.')\nlearn2d_fold2.path = cnn2dmodelpath\nlearn2d_fold2.load('xresnet34-512-PR-fold2');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learn2d_fold3 = cnn_learner(dls, xresnet34, pretrained=False, loss_func=nn.CrossEntropyLoss(), model_dir='.')\nlearn2d_fold3.path = cnn2dmodelpath\nlearn2d_fold3.load('xresnet34-512-PR-fold3');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learn2d_fold4 = cnn_learner(dls, xresnet34, pretrained=False, loss_func=nn.CrossEntropyLoss(), model_dir='.')\nlearn2d_fold4.path = cnn2dmodelpath\nlearn2d_fold4.load('xresnet34-512-PR-fold4');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model0 = learn2d_fold0.model.eval().to(device)\nmodel1 = learn2d_fold1.model.eval().to(device)\nmodel2 = learn2d_fold2.model.eval().to(device)\nmodel3 = learn2d_fold3.model.eval().to(device)\nmodel4 = learn2d_fold4.model.eval().to(device)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"RAM Until Here: 3.2 GB"},{"metadata":{},"cell_type":"markdown","source":"#### Sequence Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"from fastai.text.all import *","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_pad_idx = None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AWD_LSTM(Module):\n    \"AWD-LSTM inspired by https://arxiv.org/abs/1708.02182\"\n    initrange=0.1\n\n    def __init__(self, emb_sz,n_hid, n_layers, hidden_p=0.2, input_p=0.6, weight_p=0.5, bidir=False):\n        store_attr('emb_sz,n_hid,n_layers')\n        self.bs = 1\n        self.n_dir = 2 if bidir else 1\n        \n        self.rnns = nn.ModuleList([self._one_rnn(emb_sz if l == 0 else n_hid, (n_hid)//self.n_dir, bidir, weight_p, l) for l in range(n_layers)])\n\n        self.input_dp = RNNDropout(input_p)\n        self.hidden_dps = nn.ModuleList([RNNDropout(hidden_p) for l in range(n_layers)])\n        self.reset()\n\n    def forward(self, x, from_embeds=False):\n        \n        if from_embeds: inp = x\n        else: inp = combined_embeddings[x].to(device)\n        bs,sl = inp.shape[:2]\n        if bs!=self.bs: self._change_hidden(bs)\n\n        output = self.input_dp(inp)\n        new_hidden = []\n        for l, (rnn,hid_dp) in enumerate(zip(self.rnns, self.hidden_dps)):\n            output, new_h = rnn(output, self.hidden[l])\n            new_hidden.append(new_h)\n            if l != self.n_layers - 1: output = hid_dp(output)\n        self.hidden = to_detach(new_hidden, cpu=False, gather=False)\n        return output\n\n    def _change_hidden(self, bs):\n        self.hidden = [self._change_one_hidden(l, bs) for l in range(self.n_layers)]\n        self.bs = bs\n\n    def _one_rnn(self, n_in, n_out, bidir, weight_p, l):\n        \"Return one of the inner rnn\"\n        rnn = nn.LSTM(n_in, n_out, 1, batch_first=True, bidirectional=bidir, bias=False)\n        return WeightDropout(rnn, weight_p)\n\n    def _one_hidden(self, l):\n        \"Return one hidden state\"\n        nh = (self.n_hid) // self.n_dir\n        return (one_param(self).new_zeros(self.n_dir, self.bs, nh).to(device), one_param(self).new_zeros(self.n_dir, self.bs, nh).to(device))\n\n    def _change_one_hidden(self, l, bs):\n        if self.bs < bs:\n            nh = (self.n_hid) // self.n_dir\n            return tuple(torch.cat([h, h.new_zeros(self.n_dir, bs-self.bs, nh)], dim=1) for h in self.hidden[l])\n        if self.bs > bs: return (self.hidden[l][0][:,:bs].contiguous(), self.hidden[l][1][:,:bs].contiguous())\n        return self.hidden[l]\n\n    def reset(self):\n        \"Reset the hidden states\"\n        [r.reset() for r in self.rnns if hasattr(r, 'reset')]\n        self.hidden = [self._one_hidden(l) for l in range(self.n_layers)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lstm_width = 512 \nlayers = [lstm_width * 3] + [lstm_width] + [9]\n\nclass MultiHeadedSequenceClassifier(Module):\n    \"dim: input sequence feature dim\"\n    def __init__(self, bptt=72, input_pad_idx=input_pad_idx, n_meta=1, dim=1024):\n        \n        store_attr('input_pad_idx')\n        self.awd_lstm = AWD_LSTM(dim+n_meta, lstm_width, 2, bidir=True)\n        self.encoder = SentenceEncoder(bptt=bptt, module=self.awd_lstm, pad_idx=input_pad_idx)\n        \n        # image level preds\n        self.seq_head = LinearDecoder(1, lstm_width, bias=True)\n \n        # exam level preds\n        self.exam_head = PoolingLinearClassifier(layers, ps=[0.4, 0.1], bptt=bptt)\n        \n    \n    def forward(self, x):\n        out, mask = self.encoder(x) \n       \n        # img level out\n        seq_cls_out,_,_ = self.seq_head(out)\n        seq_cls_out = seq_cls_out.squeeze(-1)\n              \n        # exam level out\n        exam_out,_,_ = self.exam_head((out,mask))\n\n        return (seq_cls_out, exam_out)\n    \n    def predict(self, x):\n        out = self.awd_lstm(x, from_embeds=True)\n        \n        # img level out\n        seq_cls_out,_,_ = self.seq_head(out)\n        seq_cls_out = seq_cls_out.squeeze(-1)\n        \n       \n        # exam level out\n        mask = torch.zeros(x.shape[:-1]).bool().to(device)\n        exam_out,_,_ = self.exam_head((out, mask))\n        return (seq_cls_out, exam_out)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"layers1 = [512 * 3] + [512] + [12]\n\nclass MultiHeadedSoftmaxSequenceClassifier(Module):\n    \"dim: input sequence feature dim\"\n    def __init__(self, bptt=72, input_pad_idx=input_pad_idx, n_meta=2, dim=1024):\n        \n        store_attr('input_pad_idx')\n        self.awd_lstm = AWD_LSTM(dim+n_meta, 512, 2, bidir=True)\n        self.encoder = SentenceEncoder(bptt=bptt, module=self.awd_lstm, pad_idx=input_pad_idx)\n        \n        # image level preds\n        self.seq_head = LinearDecoder(1, 512, bias=True)\n \n        # exam level preds\n        self.exam_head = PoolingLinearClassifier(layers1, ps=[0.4, 0.1], bptt=bptt)\n\n#         self.posnegind_head = PoolingLinearClassifier(layers1, ps=[0.4, 0.1], bptt=bptt)\n#         self.rvlv_head = PoolingLinearClassifier(layers2, ps=[0.4, 0.1], bptt=bptt)\n#         self.lrc_head = PoolingLinearClassifier(layers3, ps=[0.4, 0.1], bptt=bptt)\n#         self.chroacute_head = PoolingLinearClassifier(layers4, ps=[0.4, 0.1], bptt=bptt)\n        \n    \n    def forward(self, x):\n        out, mask = self.encoder(x) \n       \n        # img level out\n        seq_cls_out,_,_ = self.seq_head(out)\n        seq_cls_out = seq_cls_out.squeeze(-1)\n              \n        # exam level out\n        exam_out,_,_ = self.exam_head((out,mask))\n        posneg_out, rvlv_out, lrc_out, chroacute_out = (exam_out[:,:3], \n                                                        exam_out[:,3:6], \n                                                        exam_out[:,6:9], \n                                                        exam_out[:,9:])\n\n        return (seq_cls_out, posneg_out, rvlv_out, lrc_out, chroacute_out)\n    \n    \n    def predict(self, x):\n        out = self.awd_lstm(x, from_embeds=True)\n        \n        # img level out\n        seq_cls_out,_,_ = self.seq_head(out)\n        seq_cls_out = seq_cls_out.squeeze(-1)\n        \n       \n        # exam level out\n        mask = torch.zeros(x.shape[:-1]).bool().to(device)\n        exam_out,_,_ = self.exam_head((out, mask))\n        posneg_out, rvlv_out, lrc_out, chroacute_out = (exam_out[:,:3], \n                                                        exam_out[:,3:6], \n                                                        exam_out[:,6:9], \n                                                        exam_out[:,9:])\n\n        return (seq_cls_out, posneg_out, rvlv_out, lrc_out, chroacute_out)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seq_model0 = SequentialRNN(MultiHeadedSoftmaxSequenceClassifier(bptt=256, dim=1024, n_meta=2))\nseq_model0.load_state_dict(torch.load(seqmodelpath/\"nometa_sequence_softmax_with_preds_fulldata_fold0.pth\"));\n\nseq_model1 = SequentialRNN(MultiHeadedSoftmaxSequenceClassifier(bptt=256, dim=1024, n_meta=2))\nseq_model1.load_state_dict(torch.load(seqmodelpath/\"nometa_sequence_softmax_with_preds_fulldata_fold1.pth\"));\n\nseq_model2 = SequentialRNN(MultiHeadedSoftmaxSequenceClassifier(bptt=256, dim=1024, n_meta=2))\nseq_model2.load_state_dict(torch.load(seqmodelpath/\"nometa_sequence_softmax_with_preds_fulldata_fold2.pth\"));\n\nseq_model3 = SequentialRNN(MultiHeadedSoftmaxSequenceClassifier(bptt=256, dim=1024, n_meta=2))\nseq_model3.load_state_dict(torch.load(seqmodelpath/\"nometa_sequence_softmax_with_preds_fulldata_fold3.pth\"));\n\nseq_model4 = SequentialRNN(MultiHeadedSoftmaxSequenceClassifier(bptt=256, dim=1024, n_meta=2))\nseq_model4.load_state_dict(torch.load(seqmodelpath/\"nometa_sequence_softmax_with_preds_fulldata_fold4.pth\"));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seq_model0 = seq_model0[0].to(device).eval()\nseq_model1 = seq_model1[0].to(device).eval()\nseq_model2 = seq_model2[0].to(device).eval()\nseq_model3 = seq_model3[0].to(device).eval()\nseq_model4 = seq_model4[0].to(device).eval()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"RAM Until Here: 3.4 GB"},{"metadata":{},"cell_type":"markdown","source":"### Metadata"},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_std_dict = {\n#             'img_min': [-1542.35551498553, 849.3331965009891],\n#             'img_max': [3209.4925326455914, 1138.112174280331],\n#             'img_mean': [165.7994337836255, 278.9659609535833],\n#             'img_std': [994.3087633304141, 293.05859626364196],\n#             'img_pct_window': [0.43983955119726603, 0.11747102851831802],\n            'scaled_position': [0.5078721739409284, 0.29139548181397823]\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Predict"},{"metadata":{"trusted":true},"cell_type":"code","source":"from fastai.medical.imaging import *\nimport pydicom","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class EmbeddingHook:\n    def __init__(self, m):\n        self.embeddings, self.m = tensor([]).to(device), m\n        if len(m._forward_hooks) > 0: self.reset()\n        self.hook = Hook(m, self.hook_fn, cpu=False)\n       \n    def hook_fn(self, m, inp, out): \n        \"Stack and save computed embeddings\"\n        self.embeddings = torch.cat([self.embeddings, out])\n    \n    def reset(self): \n        self.m._forward_hooks = OrderedDict()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# meta_feat_cols = ['img_min', 'img_max', 'img_mean', 'img_std', 'img_pct_window', 'scaled_position']\nmeta_feat_cols = ['scaled_position']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def minmax_scaler(o): return (o - min(o))/(max(o) - min(o))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mean, std = mean_std_dict['scaled_position']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_dcm_img_v2(dcm, windows=windows):\n    \"Read single slice in RGB\"\n    return [dcm.windowed(*w) for w in windows]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def predict_study(study_dirname):   \n    # get metadata\n    study_df = test_df.query(f'StudyInstanceUID == \"{study_dirname.stem}\"')\n    sop_ids = study_df['SOPInstanceUID'].values\n    study_files = str(testdatapath) + \"/\" + study_df['StudyInstanceUID'] + \"/\" + study_df['SeriesInstanceUID'] + '/' + sop_ids + '.dcm'\n    dcm_ds = array([pydicom.read_file(o) for o in study_files])\n\n    z_positions = array([int(o.ImagePositionPatient[-1]) for o in dcm_ds])\n    sortidxs = np.argsort(z_positions)\n\n    sop_ids = sop_ids[sortidxs]\n    dcm_ds = dcm_ds[sortidxs]\n    z_positions = z_positions[sortidxs]\n\n    #     imgs = torch.stack([read_dcm_img(o) for o in dcm_ds])\n    imgs = []\n    for o in dcm_ds: imgs += read_dcm_img_v2(o)\n    h,w = imgs[0].size()\n    imgs = torch.stack(imgs).view(-1,3,h,w)\n\n    meta_embeddings = tensor((minmax_scaler(z_positions) - mean)/std).unsqueeze(1)\n\n    emb_hook0 = EmbeddingHook(model0[1][1])\n    emb_hook1 = EmbeddingHook(model1[1][1])\n    emb_hook2 = EmbeddingHook(model2[1][1])\n    emb_hook3 = EmbeddingHook(model3[1][1])\n    emb_hook4 = EmbeddingHook(model4[1][1])\n\n\n    with torch.no_grad():\n        test_dl = learn2d_fold0.dls.test_dl(imgs.numpy(), bs=32, num_workers=2)\n\n        outs0, outs1, outs2, outs3, outs4, outs5 = [],[],[],[],[],[]\n        for xb in test_dl:\n            out0 = model0(*xb); outs0.append(out0)\n            out1 = model1(*xb); outs1.append(out1)\n            out2 = model2(*xb); outs2.append(out2)\n            out3 = model3(*xb); outs3.append(out3)\n            out4 = model4(*xb); outs4.append(out4)\n\n        outs0 = torch.cat(outs0)[:,1].view(-1,1)\n        outs1 = torch.cat(outs1)[:,1].view(-1,1)\n        outs2 = torch.cat(outs2)[:,1].view(-1,1)\n        outs3 = torch.cat(outs3)[:,1].view(-1,1)\n        outs4 = torch.cat(outs4)[:,1].view(-1,1)\n\n        seq_inp0 = torch.cat([emb_hook0.embeddings, meta_embeddings.to(device), outs0],1)\n        seq_inp1 = torch.cat([emb_hook1.embeddings, meta_embeddings.to(device), outs1],1)\n        seq_inp2 = torch.cat([emb_hook2.embeddings, meta_embeddings.to(device), outs2],1)\n        seq_inp3 = torch.cat([emb_hook3.embeddings, meta_embeddings.to(device), outs3],1)\n        seq_inp4 = torch.cat([emb_hook4.embeddings, meta_embeddings.to(device), outs4],1)\n        \n        # sequence pred\n        seq_img_preds0, posneg_preds0, rvlv_preds0, lrc_preds0, chroacute_preds0 = seq_model0.predict(seq_inp0[None,...])\n        seq_img_preds1, posneg_preds1, rvlv_preds1, lrc_preds1, chroacute_preds1 = seq_model1.predict(seq_inp1[None,...])\n        seq_img_preds2, posneg_preds2, rvlv_preds2, lrc_preds2, chroacute_preds2 = seq_model2.predict(seq_inp2[None,...])\n        seq_img_preds3, posneg_preds3, rvlv_preds3, lrc_preds3, chroacute_preds3 = seq_model3.predict(seq_inp3[None,...])\n        seq_img_preds4, posneg_preds4, rvlv_preds4, lrc_preds4, chroacute_preds4 = seq_model4.predict(seq_inp4[None,...])\n        \n        # tta pred - flip sequence order\n        seq_img_preds0_tta, posneg_preds0_tta, rvlv_preds0_tta, lrc_preds0_tta, chroacute_preds0_tta = seq_model0.predict(seq_inp0.flip(dims=[0])[None,...])\n        seq_img_preds1_tta, posneg_preds1_tta, rvlv_preds1_tta, lrc_preds1_tta, chroacute_preds1_tta = seq_model1.predict(seq_inp1.flip(dims=[0])[None,...])\n        seq_img_preds2_tta, posneg_preds2_tta, rvlv_preds2_tta, lrc_preds2_tta, chroacute_preds2_tta = seq_model2.predict(seq_inp2.flip(dims=[0])[None,...])\n        seq_img_preds3_tta, posneg_preds3_tta, rvlv_preds3_tta, lrc_preds3_tta, chroacute_preds3_tta = seq_model3.predict(seq_inp3.flip(dims=[0])[None,...])\n        seq_img_preds4_tta, posneg_preds4_tta, rvlv_preds4_tta, lrc_preds4_tta, chroacute_preds4_tta = seq_model4.predict(seq_inp4.flip(dims=[0])[None,...])\n            \n        # flip back to align image slices\n        seq_img_preds0_tta = seq_img_preds0_tta.flip(dims=[1])\n        seq_img_preds1_tta = seq_img_preds1_tta.flip(dims=[1])\n        seq_img_preds2_tta = seq_img_preds2_tta.flip(dims=[1])\n        seq_img_preds3_tta = seq_img_preds3_tta.flip(dims=[1])\n        seq_img_preds4_tta = seq_img_preds4_tta.flip(dims=[1])\n        \n        \n        seq_img_preds0,posneg_preds0,rvlv_preds0,lrc_preds0,chroacute_preds0 =  ((seq_img_preds0+seq_img_preds0_tta)/2, \n                                                                                 (posneg_preds0+posneg_preds0_tta)/2, \n                                                                                 (rvlv_preds0+rvlv_preds0_tta)/2, \n                                                                                 (lrc_preds0+lrc_preds0_tta)/2, \n                                                                                 (chroacute_preds0+chroacute_preds0_tta)/2)\n        \n        seq_img_preds1,posneg_preds1,rvlv_preds1,lrc_preds1,chroacute_preds1 =  ((seq_img_preds1+seq_img_preds1_tta)/2, \n                                                                                 (posneg_preds1+posneg_preds1_tta)/2, \n                                                                                 (rvlv_preds1+rvlv_preds1_tta)/2, \n                                                                                 (lrc_preds1+lrc_preds1_tta)/2, \n                                                                                 (chroacute_preds1+chroacute_preds1_tta)/2)\n\n        seq_img_preds2,posneg_preds2,rvlv_preds2,lrc_preds2,chroacute_preds2 =  ((seq_img_preds2+seq_img_preds2_tta)/2, \n                                                                                 (posneg_preds2+posneg_preds2_tta)/2, \n                                                                                 (rvlv_preds2+rvlv_preds2_tta)/2, \n                                                                                 (lrc_preds2+lrc_preds2_tta)/2, \n                                                                                 (chroacute_preds2+chroacute_preds2_tta)/2)\n        \n        seq_img_preds3,posneg_preds3,rvlv_preds3,lrc_preds3,chroacute_preds3 =  ((seq_img_preds3+seq_img_preds3_tta)/2, \n                                                                                 (posneg_preds3+posneg_preds3_tta)/2, \n                                                                                 (rvlv_preds3+rvlv_preds3_tta)/2, \n                                                                                 (lrc_preds3+lrc_preds3_tta)/2, \n                                                                                 (chroacute_preds3+chroacute_preds3_tta)/2)\n        \n        seq_img_preds4,posneg_preds4,rvlv_preds4,lrc_preds4,chroacute_preds4 =  ((seq_img_preds4+seq_img_preds4_tta)/2, \n                                                                                 (posneg_preds4+posneg_preds4_tta)/2, \n                                                                                 (rvlv_preds4+rvlv_preds4_tta)/2, \n                                                                                 (lrc_preds4+lrc_preds4_tta)/2, \n                                                                                 (chroacute_preds4+chroacute_preds4_tta)/2)\n        \n        \n    seq_img_preds0 = seq_img_preds0.sigmoid()\n    seq_img_preds1 = seq_img_preds1.sigmoid()\n    seq_img_preds2 = seq_img_preds2.sigmoid()\n    seq_img_preds3 = seq_img_preds3.sigmoid()\n    seq_img_preds4 = seq_img_preds4.sigmoid()\n\n    posneg_preds = ((\n                      posneg_preds0.softmax(1) \n                    + posneg_preds1.softmax(1)\n                    + posneg_preds2.softmax(1)\n                    + posneg_preds3.softmax(1)\n                    + posneg_preds4.softmax(1)) / 5)[0]\n\n    rvlv_preds = ((\n                      rvlv_preds0.softmax(1) \n                    + rvlv_preds1.softmax(1)\n                    + rvlv_preds2.softmax(1)\n                    + rvlv_preds3.softmax(1)\n                    + rvlv_preds4.softmax(1)) / 5)[0]\n\n    lrc_preds = ((\n                      lrc_preds0.sigmoid() \n                    + lrc_preds1.sigmoid()\n                    + lrc_preds2.sigmoid()\n                    + lrc_preds3.sigmoid()\n                    + lrc_preds4.sigmoid()) / 5)[0]\n\n    chroacute_preds = ((\n                      chroacute_preds0.softmax(1) \n                    + chroacute_preds1.softmax(1)\n                    + chroacute_preds2.softmax(1)\n                    + chroacute_preds3.softmax(1)\n                    + chroacute_preds4.softmax(1)) / 5)[0]\n\n    # pos, neg, ind\n    _, negative, indeterminate = posneg_preds\n    \n    # rvlv >= 1, rvlv < 1, neither\n    rvlv_gte, rvlv_lt, _ = rvlv_preds\n    \n    # left, right, central\n    left_pe, right_pe, central_pe = lrc_preds\n    \n    # chronic, acute_and_chronic, neither\n    chro_pe, acute_chro_pe, _ = chroacute_preds\n    \n    seq_exam_preds = torch.stack([negative, rvlv_gte, rvlv_lt, left_pe, chro_pe, right_pe, acute_chro_pe, central_pe, indeterminate])\n\n    seq_img_preds = ((seq_img_preds0\n                    + seq_img_preds1\n                    + seq_img_preds2\n                    + seq_img_preds3\n                    + seq_img_preds4) / 5)[0]\n    \n    return (sop_ids, to_cpu(seq_img_preds), to_cpu(seq_exam_preds))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"RAM Until Here: 4.3 GB"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_study_res(sid, sop_ids, img_preds, exam_preds):\n    \"Get preds from 2D and 3D cnn models\"\n    sub_res = []\n    for sopid, p in zip(sop_ids, to_np(img_preds)):\n        sub_res.append((sopid, p))\n            \n    # exam probas (same order for 3D and no meta sequence models)\n    target_cols = [\n        'negative_exam_for_pe', # exam level\n        'rv_lv_ratio_gte_1', # exam level\n        'rv_lv_ratio_lt_1', # exam level\n        'leftsided_pe', # exam level\n        'chronic_pe', # exam level\n        'rightsided_pe', # exam level\n        'acute_and_chronic_pe', # exam level\n        'central_pe', # exam level\n        'indeterminate' # exam level\n    ]\n    \n    \n    for tcol, p in zip(target_cols, to_np(exam_preds)):\n        sub_res.append((f\"{sid}_{tcol}\", p))\n    return sub_res","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ndo_full = False\nn = 20\n\nif Path('../input/rsna-str-pulmonary-embolism-detection/train').exists() and not do_full: \n    test_study_dirnames = [datapath/'test'/o for o in test_df['StudyInstanceUID'].unique()]\n    test_study_dirnames = np.random.choice(test_study_dirnames, n, replace=False)\n\nsub_res = []\nfor study_dirname in test_study_dirnames:\n    sop_ids, seq_img_preds, seq_exam_preds = predict_study(study_dirname)\n    study_res = get_study_res(study_dirname.stem, sop_ids, seq_img_preds, seq_exam_preds)\n    sub_res += study_res","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"final_sub_df = pd.DataFrame(sub_res, columns=['id', 'label'])\nfinal_sub_df['label']  = np.clip(final_sub_df['label'], 0.0001, 0.9999)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"final_sub_df.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Consistency Check"},{"metadata":{"trusted":true},"cell_type":"code","source":"def check_consistency(sub, test):\n    \n    '''\n    Checks label consistency and returns the errors\n    \n    Args:\n    sub   = submission dataframe (pandas)\n    test  = test.csv dataframe (pandas)\n    '''\n    \n    # EXAM LEVEL\n    df_tmp = sub.loc[sub.id.str.contains('_', regex = False)].reset_index(drop = True)\n    df_tmp['StudyInstanceUID'] = df_tmp['id'].apply(lambda x:str(x).split('_')[0])\n    df_tmp['label_type'] = df_tmp['id'].apply(lambda x:'_'.join(str(x).split('_')[1:]))\n    df_exam = df_tmp.pivot(index = 'StudyInstanceUID', columns = 'label_type', values = 'label')\n    \n    # IMAGE LEVEL\n    df_image = sub.loc[sub.id.isin(test.SOPInstanceUID)].reset_index(drop = True)\n    df_image = df_image.merge(test, how = 'left', left_on = 'id', right_on = 'SOPInstanceUID')\n    df_image.rename(columns = {\"label\": \"pe_present_on_image\"}, inplace = True)\n    del df_image['id']\n    \n    # MERGER\n    df = df_exam.merge(df_image, how = 'left', on = 'StudyInstanceUID')\n    ids    = ['StudyInstanceUID', 'SeriesInstanceUID', 'SOPInstanceUID']\n    labels = [c for c in df.columns if c not in ids]\n    df = df[ids + labels]\n    \n    # SPLIT NEGATIVE AND POSITIVE EXAMS\n    df['positive_images_in_exam'] = df['StudyInstanceUID'].map(df.groupby(['StudyInstanceUID']).pe_present_on_image.max())\n    df_pos = df.loc[df.positive_images_in_exam >  0.5]\n    df_neg = df.loc[df.positive_images_in_exam <= 0.5]\n    \n    # CHECKING CONSISTENCY OF POSITIVE EXAM LABELS\n    rule1a = df_pos.loc[((df_pos.rv_lv_ratio_lt_1  >  0.5)  & \n                         (df_pos.rv_lv_ratio_gte_1 >  0.5)) | \n                        ((df_pos.rv_lv_ratio_lt_1  <= 0.5)  & \n                         (df_pos.rv_lv_ratio_gte_1 <= 0.5))].reset_index(drop = True)\n    rule1a['broken_rule'] = '1a'\n    rule1b = df_pos.loc[(df_pos.central_pe    <= 0.5) & \n                        (df_pos.rightsided_pe <= 0.5) & \n                        (df_pos.leftsided_pe  <= 0.5)].reset_index(drop = True)\n    rule1b['broken_rule'] = '1b'\n    rule1c = df_pos.loc[(df_pos.acute_and_chronic_pe > 0.5) & \n                        (df_pos.chronic_pe           > 0.5)].reset_index(drop = True)\n    rule1c['broken_rule'] = '1c'\n    rule1d = df_pos.loc[(df_pos.indeterminate        > 0.5) | \n                        (df_pos.negative_exam_for_pe > 0.5)].reset_index(drop = True)\n    rule1d['broken_rule'] = '1d'\n\n    # CHECKING CONSISTENCY OF NEGATIVE EXAM LABELS\n    rule2a = df_neg.loc[((df_neg.indeterminate        >  0.5)  & \n                         (df_neg.negative_exam_for_pe >  0.5)) | \n                        ((df_neg.indeterminate        <= 0.5)  & \n                         (df_neg.negative_exam_for_pe <= 0.5))].reset_index(drop = True)\n    rule2a['broken_rule'] = '2a'\n    rule2b = df_neg.loc[(df_neg.rv_lv_ratio_lt_1     > 0.5) | \n                        (df_neg.rv_lv_ratio_gte_1    > 0.5) |\n                        (df_neg.central_pe           > 0.5) | \n                        (df_neg.rightsided_pe        > 0.5) | \n                        (df_neg.leftsided_pe         > 0.5) |\n                        (df_neg.acute_and_chronic_pe > 0.5) | \n                        (df_neg.chronic_pe           > 0.5)].reset_index(drop = True)\n    rule2b['broken_rule'] = '2b'\n    \n    # MERGING INCONSISTENT PREDICTIONS\n    errors = pd.concat([rule1a, rule1b, rule1c, rule1d, rule2a, rule2b], axis = 0)\n    \n    # OUTPUT\n    print('Found', len(errors), 'inconsistent predictions')\n    if len(errors) > 0:\n        print(errors.broken_rule.value_counts())\n        \n    return errors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nconsistency_df0 = check_consistency(final_sub_df, test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"consistency_df0.StudyInstanceUID.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"consistency_df0.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"consistency_df0['broken_rule'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.special import softmax","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def solve_rule_1a (sub, errors, log=False):\n    '''\n    ((df_pos.rv_lv_ratio_lt_1  >  0.5)  & \n                         (df_pos.rv_lv_ratio_gte_1 >  0.5)) | \n                        ((df_pos.rv_lv_ratio_lt_1  <= 0.5)  & \n                         (df_pos.rv_lv_ratio_gte_1 <= 0.5))\n    '''\n    \n    studies_break_rule = errors.query(\"((rv_lv_ratio_lt_1  >  0.5) & (rv_lv_ratio_gte_1 >  0.5)) | \\\n                                  ((rv_lv_ratio_lt_1  <=  0.5) & (rv_lv_ratio_gte_1 <=  0.5))\")['StudyInstanceUID'].unique()\n\n    if log: \n        print (studies_break_rule,studies_break_rule.shape)\n    \n    for i in studies_break_rule:\n        \n        idx = sub[sub['id'].str.contains(i)].index.values\n        old = sub.loc[idx[[1,2]], 'label'].values\n        if old[0] > 0.5 and old[1] > 0.5:\n            new = softmax (old)\n        elif old[0] == 0.5 and old[1] == 0.5:\n            pass\n        else:\n            new = softmax (old)\n        \n        # softmax collateral effect\n        if new[0]== 0.5 and new[1] == 0.5:\n            pass\n        \n        if log: print (old , new)\n        sub.loc[idx[[1,2]], 'label'] = new","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def solve_rule_1b (sub, errors, log=False):\n    \n    studies_break_rule = errors.query(\"(central_pe <= 0.5) & (rightsided_pe <= 0.5) & (leftsided_pe  <= 0.5)\")['StudyInstanceUID'].unique()\n\n    if log: \n        print (studies_break_rule,studies_break_rule.shape)\n    \n    for i in studies_break_rule:\n        \n        idx = sub[sub['id'].str.contains(i)].index.values\n        #print (i, idx[[3,5,7]])\n        #print (sub.loc[(sub['id'].str.contains(i))])\n        old = sub.loc[idx[[3,5,7]], 'label'].values\n        best = np.argmax (old)\n        new = old.copy()\n        \n        if (old[0] == 0.5) and (old[1] == 0.5) and (old[2] == 0.5):\n            pass\n        else:\n            new[best] = new[best] + 0.5\n            new = np.clip (new , 0, 0.51) # remember that is possible to have all > 0.5\n            \n        if log: \n            print (old, new)\n        sub.loc[idx[[3,5,7]], 'label'] = new","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def solve_rule_1c (sub, errors, log=False):\n    \n    studies_break_rule = errors.query(\"(acute_and_chronic_pe > 0.5) & (chronic_pe > 0.5)\")['StudyInstanceUID'].unique()\n\n    if log: \n        print (studies_break_rule, studies_break_rule.shape)\n    \n    for i in studies_break_rule:\n        \n        idx = sub[sub['id'].str.contains(i)].index.values\n        #print (i, idx[1:3])\n        #print (sub.loc[(sub['id'].str.contains(i))])\n        old = sub.loc[idx[[4,6]], 'label'].values\n        new = softmax (old)\n        if log: \n            print (old , new)\n            \n        sub.loc[idx[[4,6]], 'label'] == new","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def solve_rule_1d(sub, errors, log=False):\n    \n    studies_break_rule = errors.query(\"(indeterminate > 0.5) | (negative_exam_for_pe > 0.5)\")['StudyInstanceUID'].unique()\n\n    if log: \n        print (studies_break_rule, studies_break_rule.shape)\n    \n    for i in studies_break_rule:\n        \n        idx = sub[sub['id'].str.contains(i)].index.values\n        #print (sub.loc[(sub['id'].str.contains(i))])\n        \n        old = sub.loc[idx[[0,-1]], 'label'].values\n        #if (old_neg_ind[0] == old_neg_ind[1]):\n        #    if log: print ('negative = 0.5 = indet')\n        #    new_neg_ind = [0.7 , 0.3]\n        #else:\n        #    new_neg_ind = softmax (old_neg_ind)\n        \n        new = np.clip (old, 0. , 0.4)\n        if log: \n            print (old , new)\n            \n        sub.loc[idx[[0,-1]], 'label'] = new","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def solve_negative (sub, errors):\n    # clip all images to (0, 0.5) that are predicted as negative exam\n    negative_sopids = errors.query(\"negative_exam_for_pe > 0.5\")['SOPInstanceUID'].values\n    neg_clip_idxs = sub[sub['id'].isin(negative_sopids)].index\n    sub.loc[neg_clip_idxs, 'label'] = np.clip(sub.loc[neg_clip_idxs, 'label'], 0, 0.5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsolve_negative (sub=final_sub_df, errors=consistency_df0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nconsistency_df0 = check_consistency(sub=final_sub_df, test=test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsolve_rule_1a (sub=final_sub_df, errors=consistency_df0, log=False)\nsolve_rule_1b (sub=final_sub_df, errors=consistency_df0, log=False)\nsolve_rule_1c (sub=final_sub_df, errors=consistency_df0, log=False)\nsolve_rule_1d (sub=final_sub_df, errors=consistency_df0, log=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"consistency_df_final = check_consistency(final_sub_df, test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submit_fixed = True\nif submit_fixed:\n    final_sub_df.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Submission Stats"},{"metadata":{"trusted":true},"cell_type":"code","source":"final_sub_df['label_prefix'] = final_sub_df.apply(lambda o: \"_\".join(o['id'].split(\"_\")[1:]), 1)\nstats = final_sub_df.groupby(\"label_prefix\").agg(['min', 'max', 'mean', 'median'])\nstats","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"final_sub_df.query(\"label_prefix == ''\")['label'].hist();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}