{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"notebook for discussion at:  \nhttps://www.kaggle.com/competitions/vesuvius-challenge-ink-detection/discussion/407972","metadata":{}},{"cell_type":"code","source":"my_lib_dir ='/kaggle/input/ink-00/my_lib'\n\nimport sys\nsys.path.append(my_lib_dir)\nsys.path.append('/kaggle/input/pretrainedmodels/pretrainedmodels-0.7.4')\nsys.path.append('/kaggle/input/efficientnet-pytorch/EfficientNet-PyTorch-master')\nsys.path.append('/kaggle/input/timm-pytorch-image-models/pytorch-image-models-master')\nsys.path.append('/kaggle/input/segmentation-models-pytorch/segmentation_models.pytorch-master')\nsys.path.append('/kaggle/input/einops/einops-master')\n\nfrom helper import *\n\nimport numpy as np\nimport pandas as pd\n\nfrom collections import defaultdict\nfrom glob import glob\nimport PIL.Image as Image\nImage.MAX_IMAGE_PIXELS = 10000000000  # Ignore PIL warnings about large images\n\nimport cv2\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nfrom einops import rearrange, reduce, repeat\nimport segmentation_models_pytorch as smp\nfrom segmentation_models_pytorch.decoders.unet.decoder import UnetDecoder\nfrom timm.models.resnet import *\n\n\nimport matplotlib\nimport matplotlib.pyplot as plt\n#matplotlib.use('TkAgg')\n%matplotlib inline \n  \nprint('import ok !!!')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-05-14T05:30:30.604431Z","iopub.execute_input":"2023-05-14T05:30:30.605048Z","iopub.status.idle":"2023-05-14T05:30:38.612180Z","shell.execute_reply.started":"2023-05-14T05:30:30.605014Z","shell.execute_reply":"2023-05-14T05:30:38.610755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Config(object):\n\tmode = [\n\t\t'train', #'test', 'train'\n\t]\n\tcrop_size  = 224\n\tcrop_depth = 8+4\n\tone_depth  = 8  #6+4\n\nCFG = Config()\nCFG.fragment_z0 = 65//2-5-2 #-1\nCFG.fragment_z1 = CFG.fragment_z0+CFG.crop_depth #+2\nCFG.is_tta = True\n\nif 'train' in CFG.mode:\n\tCFG.stride = CFG.crop_size//4\nif 'test' in CFG.mode:\n\tCFG.stride = CFG.crop_size//8\n","metadata":{"execution":{"iopub.status.busy":"2023-05-14T05:30:38.614351Z","iopub.execute_input":"2023-05-14T05:30:38.614700Z","iopub.status.idle":"2023-05-14T05:30:38.621572Z","shell.execute_reply.started":"2023-05-14T05:30:38.614669Z","shell.execute_reply":"2023-05-14T05:30:38.620373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## dataset ##\nif 'train' in CFG.mode:\n\tdata_dir = '/kaggle/input/vesuvius-challenge-ink-detection/train'\n\tvalid_id =[\n\t    '2a',\n\t]\n\nif 'test' in CFG.mode: \n\tdata_dir = '/kaggle/input/vesuvius-challenge-ink-detection/test'\n\n\tvalid_id = glob(f'{data_dir}/*')\n\tvalid_id = sorted(valid_id)\n\tvalid_id = [f.split('/')[-1] for f in valid_id]\n\nprint('data_dir', data_dir)\nprint('valid_id', valid_id)\n\n\n\ndef do_binarise(m, threshold=0.5):\n    m = m-m.min()\n    m = m/(m.max()+1e-7)\n    m = (m>threshold).astype(np.float32)\n    return m\n\ndef read_data(fragment_id, z0=CFG.fragment_z0, z1=CFG.fragment_z1):\n    volume = []\n    start_timer = timer()\n    for i in range(z0,z1):\n        v = np.array(Image.open(f'{data_dir}/{fragment_id}/surface_volume/{i:02d}.tif'), dtype=np.uint16)\n        v = (v >> 8).astype(np.uint8)\n        #v = (v / 65535.0 * 255).astype(np.uint8)\n        volume.append(v)\n        print(f'\\r @ read_data(): volume{fragment_id}  {time_to_str(timer() - start_timer, \"sec\")}', end='', flush=True)\n    #print('')\n    volume = np.stack(volume, -1)\n    height, width, depth = volume.shape\n    #print(f'fragment_id={fragment_id} volume: {volume.shape}')\n\n    #---\n    mask = cv2.imread(f'{data_dir}/{fragment_id}/mask.png',cv2.IMREAD_GRAYSCALE)\n    mask = do_binarise(mask)\n\n    if 'train' in CFG.mode:\n        ir    = cv2.imread(f'{data_dir}/{fragment_id}/ir.png',cv2.IMREAD_GRAYSCALE)\n        label = cv2.imread(f'{data_dir}/{fragment_id}/inklabels.png',cv2.IMREAD_GRAYSCALE)\n        ir    = ir/255\n        label = do_binarise(label)\n\n    if 'test' in CFG.mode:\n        ir = None\n        label = None\n\n    d = dotdict(\n        fragment_id = fragment_id,\n        volume = volume,\n        ir     = ir,\n        label  = label,\n        mask   = mask,\n    )\n    return d\n\ndef read_data1(fragment_id):\n\tif fragment_id=='2a':\n\t\td = read_data('2')\n\t\ty = { #ab split\n\t\t\t'1': 4560,\n\t\t\t'2': 9456,\n\t\t\t'3': 4060,\n\t\t}['2']\n\t\td = dotdict(\n\t\t\tfragment_id='2a',\n\t\t\tvolume  = d.volume[:y],\n\t\t\tir      = d.ir[:y],\n\t\t\tlabel   = d.label[:y],\n\t\t\tmask    = d.mask[:y],\n\t\t)\n\telse:\n\t\td = read_data(fragment_id)\n\treturn d\n\ndef run_check_data():\n    d=read_data1(valid_id[0])#valid_id[0]\n    print('')\n    print('fragment_id:', d.fragment_id)\n    print('volume:', d.volume.shape, d.volume.min(), d.volume.max())\n    print('mask  :', d.mask.shape, d.mask.min(), d.mask.max())\n    if 'train' in CFG.mode:\n        print('ir    :', d.ir.shape, d.ir.min(), d.ir.max())\n        print('label :', d.label.shape, d.label.min(), d.label.max())\n\n#run_check_data()\nprint('data ok !!!')","metadata":{"execution":{"iopub.status.busy":"2023-05-14T05:30:38.627252Z","iopub.execute_input":"2023-05-14T05:30:38.627801Z","iopub.status.idle":"2023-05-14T05:30:38.649668Z","shell.execute_reply.started":"2023-05-14T05:30:38.627767Z","shell.execute_reply":"2023-05-14T05:30:38.648645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## model ##\n\n\nclass SmpUnetDecoder(UnetDecoder):\n\tdef __init__(self, **kwargs):\n\t\tsuper(SmpUnetDecoder, self).__init__(\n\t\t\t**kwargs)\n\n\tdef forward(self, encoder):\n\t\tfeature = encoder[::-1]  # reverse channels to start from head of encoder\n\t\thead = feature[0]\n\t\tskip = feature[1:] + [None]\n\t\td = self.center(head)\n\n\t\tdecoder = []\n\t\tfor i, decoder_block in enumerate(self.blocks):\n\t\t\t# print(i, d.shape, skip[i].shape if skip[i] is not None else 'none')\n\t\t\t# print(decoder_block.conv1[0])\n\t\t\t# print('')\n\t\t\ts = skip[i]\n\t\t\td = decoder_block(d, s)\n\t\t\tdecoder.append(d)\n\n\t\tlast  = d\n\t\treturn last, decoder\n\nclass Net(nn.Module):\n\tdef __init__(self,):\n\t\tsuper().__init__()\n\n\t\tconv_dim=64\n\t\tencoder_dim  = [conv_dim] + [64, 128, 256, 512 ]\n\t\tself.encoder = resnet34d(pretrained=False,in_chans=CFG.one_depth)\n\n\t\tself.decoder = SmpUnetDecoder(\n\t\t\tencoder_channels=[0] + encoder_dim,\n\t\t\tdecoder_channels=[256, 128, 64, 32, 16],\n\t\t\tn_blocks=5,\n\t\t\tuse_batchnorm=True,\n\t\t\tcenter=False,\n\t\t\tattention_type=None,\n\t\t)\n\t\tself.logit = nn.Conv2d(16,1,kernel_size=1)\n\n\t\t#-- pool attention weight\n\t\tself.weight = nn.ModuleList([\n\t\t\tnn.Sequential(\n\t\t\t\tnn.Conv2d(dim, dim, kernel_size=3, padding=1),\n\t\t\t\tnn.ReLU(inplace=True),\n\t\t\t) for dim in encoder_dim\n\t\t])\n\n\tdef forward(self, batch):\n\t\tv = batch['volume']\n\t\tB,C,H,W = v.shape\n\t\tvv = [\n\t\t\tv[:,i:i+CFG.one_depth] for i in [0,2,4]\n\t\t]\n\t\tK = len(vv)\n\t\tx = torch.cat(vv,0)\n\t\t#x = v\n\n\t\t# ----\n\t\tencoder = []\n\t\tx = self.encoder.conv1(x)\n\t\tx = self.encoder.bn1(x)\n\t\tx = self.encoder.act1(x)   ; encoder.append(x)\n\t\tx = F.avg_pool2d(x,kernel_size=2,stride=2)\n\t\tx = self.encoder.layer1(x) ; encoder.append(x)\n\t\tx = self.encoder.layer2(x) ; encoder.append(x)\n\t\tx = self.encoder.layer3(x) ; encoder.append(x)\n\t\tx = self.encoder.layer4(x) ; encoder.append(x)\n\t\t#print('encoder', [f.shape for f in encoder])\n\n\t\t#encode pooling -------\n\t\t#<todo> add positional encode (z slice no.)\n\t\tfor i in range(len(encoder)):\n\t\t\te = encoder[i]\n\t\t\tf = self.weight[i](encoder[i])\n\t\t\t_, c, h, w = f.shape\n\t\t\tf = rearrange(f, '(K B) c h w -> B K c h w', K=K, B=B, h=h, w=w) #f.reshape(B, K, c, h, w)\n\t\t\te = rearrange(e, '(K B) c h w -> B K c h w', K=K, B=B, h=h, w=w) #e.reshape(B, K, c, h, w)\n\t\t\tw = F.softmax(f, 1)\n\t\t\te = (w * e).sum(1)\n\t\t\tencoder[i] = e\n\n\t\t# ---\n\t\tlast, decoder = self.decoder(encoder)\n\t\t#print('decoder',[f.shape for f in decoder])\n\t\t#print('last',last.shape)\n\t\tlogit = self.logit(last)\n\n\t\toutput = {\n\t\t\t'ink' : torch.sigmoid(logit),\n\t\t}\n\t\treturn output\n\n\ndef run_check_net():\n\n    height,width =  CFG.crop_size, CFG.crop_size\n    depth = CFG.crop_depth\n    batch_size = 2\n\n    batch = {\n        'volume' : torch.from_numpy( np.random.choice(256, (batch_size, depth, height, width))).float(),#.cuda()\n    }\n    net = Net()#.cuda()\n\n    with torch.no_grad():\n        with torch.cuda.amp.autocast(enabled=True):\n            output = net(batch)\n\n    #---\n    print('batch')\n    for k, v in batch.items():\n        print(f'{k:>32} : {v.shape} ')\n\n    print('output')\n    for k, v in output.items():\n        print(f'{k:>32} : {v.shape} ')\n\nrun_check_net()\nprint('net ok !!!')","metadata":{"execution":{"iopub.status.busy":"2023-05-14T05:30:38.651174Z","iopub.execute_input":"2023-05-14T05:30:38.651840Z","iopub.status.idle":"2023-05-14T05:30:40.232471Z","shell.execute_reply.started":"2023-05-14T05:30:38.651764Z","shell.execute_reply":"2023-05-14T05:30:40.231359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# infer here !!!!\n#https://gist.github.com/janpaul123/ca3477c1db6de4346affca37e0e3d5b0\ndef mask_to_rle(mask):\n    m = mask.reshape(-1)\n    #m = np.where(mask > threshold, 1, 0).astype(np.uint8)\n\n    s = np.array((m[:-1] == 0) & (m[1:] == 1))\n    e = np.array((m[:-1] == 1) & (m[1:] == 0))\n\n    s_index = np.where(s)[0] + 2\n    e_index = np.where(e)[0] + 2\n    length = e_index - s_index\n    rle = ' '.join(map(str, sum(zip(s_index, length), ())))\n    return rle\n\ndef metric_to_text(ink, label):\n    text = []\n\n    p = ink.reshape(-1)\n    t = label.reshape(-1)\n    pos = np.log(np.clip(p,1e-7,1))\n    neg = np.log(np.clip(1-p,1e-7,1))\n    bce = -(t*pos +(1-t)*neg).mean()\n    text.append(f'bce={bce:0.5f}')\n\n\n    #print(f'{threshold:0.1f}, {precision:0.3f}, {recall:0.3f}, {fpr:0.3f},  {dice:0.3f},  {score:0.3f}')\n    text.append('th   prec   recall   fpr   dice   score')\n    text.append('---------------------------------------')\n    for threshold in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]:\n        p = ink.reshape(-1)\n        t = label.reshape(-1)\n        p = (p > threshold).astype(np.float32)\n        t = (t > 0.5).astype(np.float32)\n\n        tp = p * t\n        precision = tp.sum() / (p.sum() + 0.0001)\n        recall = tp.sum() / t.sum()\n\n        fp = p * (1 - t)\n        fpr = fp.sum() / (1 - t).sum()\n\n        beta = 0.5\n        #  0.2*1/recall + 0.8*1/prec\n        score = beta * beta / (1 + beta * beta) * 1 / recall + 1 / (1 + beta * beta) * 1 / precision\n        score = 1 / score\n\n        dice = 2 * tp.sum() / (p.sum() + t.sum())\n\n        # print(fold, threshold, precision, recall, fpr,  score)\n        text.append( f'{threshold:0.1f}, {precision:0.3f}, {recall:0.3f}, {fpr:0.3f},  {dice:0.3f},  {score:0.3f}')\n    text = '\\n'.join(text)\n    return text\n\ndef infer_one(net, d):\n    num_net = len(net)\n    for i in range(num_net):\n        net[i] = net[i].cuda()\n        net[i] = net[i].eval()\n\n    #get coord\n    size   = CFG.crop_size\n    stride = CFG.stride\n    H,W,D  = d.volume.shape\n\n    x = np.arange(0,W-size+1,stride)\n    y = np.arange(0,H-size+1,stride)\n    x,y = np.meshgrid(x,y)\n    xy  = np.stack([x,y],-1).reshape(-1,2)\n    print('H,W,len(xy)',H,W,len(xy))\n\n    #---\n    probability = np.zeros((H,W))\n    count = np.zeros((H,W))\n\n    start_timer = timer()\n    batch_iter = np.array_split(xy, len(xy)//32)\n    for t, xy0 in enumerate(batch_iter):\n        #print('\\r', t, len(batch_iter), end='')\n        crop_size  = CFG.crop_size\n\n        volume =[]\n        for x0,y0 in xy0 :\n            v = d.volume[y0:y0 + crop_size, x0:x0 + crop_size]\n            volume.append(v)\n        volume = np.stack(volume)\n        volume = np.ascontiguousarray(volume.transpose(0,3,1,2))\n        volume = volume/255\n        volume = torch.from_numpy(volume).float().cuda()\n        ##print(volume.shape)\n\n        batch = { 'volume': volume }\n\n        k = 0\n        c = 0\n        with torch.no_grad():\n            with torch.cuda.amp.autocast(enabled=True):\n                for i in range(num_net):\n                    if 0:\n                        output = net[i](batch)\n                        k += output['ink'].data.cpu().numpy()\n                        c += 1\n\n                    #--\n                    #TTA <todo>\n                    if CFG.is_tta: #tta \n                        v = [\n                            volume,\n                            torch.rot90(volume, k=1, dims=(-2, -1)),\n                            torch.rot90(volume, k=2, dims=(-2, -1)),\n                            torch.rot90(volume, k=3, dims=(-2, -1)),\n                        ]\n                        K=len(v)\n                        batch = {\n                            'volume': torch.cat(v,0)\n                        }\n                        output = net[i](batch)\n                        ink = output['ink']\n\n                        B,_,h,w = volume.shape\n                        ink = ink.reshape(K, B, 1, h, w)\n                        ink = [\n                            ink[0],\n                            torch.rot90(ink[1], k=-1, dims=(-2, -1)),\n                            torch.rot90(ink[2], k=-2, dims=(-2, -1)),\n                            torch.rot90(ink[3], k=-3, dims=(-2, -1)),\n                        ]\n                        ink = torch.stack(ink, dim=0)\n                        ink = ink.mean(0)\n\n                        k += ink.data.cpu().numpy()\n                        c += 1\n                    #--\n        k = k/c\n        ##print(k.shape)\n\n        batch_size = len(k)\n        for b in range(batch_size):\n            x0,y0 = xy0[b]\n            probability[y0:y0 + crop_size, x0:x0 + crop_size] += k[b,0]\n            count[y0:y0 + crop_size, x0:x0 + crop_size] += 1\n        print(f'\\r @infer_one(): {t} / {len(batch_iter)} : {time_to_str(timer() - start_timer, \"sec\")}', end='', flush=True)\n    print('')\n    probability = probability/(count+0.000001)\n    return probability\n\n\n######################################\ncheckpoint=[\n    '/kaggle/input/ink-00-weight/run12-ref-encoder-pool-resnet34d-sample4f-00-fold-2a-00002672.model.pth',\n]\n\n#----\nnet = []\nfor i,f in enumerate(checkpoint):\n    n = Net()\n    f = torch.load(f, map_location=lambda storage, loc: storage)\n    print(n.load_state_dict(f['state_dict'], strict=True))  # True\n    net.append(n)\n\n#----\nsubmission = defaultdict(list)\nfor t,fragment_id in enumerate(valid_id):\n    d = read_data1(fragment_id)\n    \n    print('==================================')\n    print('fragment_id', d.fragment_id)\n    print('\\tmask', d.mask.shape)\n    print('\\tvolume', d.volume.shape)\n    print('CFG.stride', CFG.stride)\n    print('CFG.crop_size', CFG.crop_size)  \n    print('')\n\n    probability = infer_one(net, d)\n    print('probability', probability.shape)\n\n    probability = d.mask*probability\n    predict = (probability>0.5).astype(np.uint8)\n    \n    #----\n    submission['Id'].append(fragment_id)\n    submission['Predicted'].append(mask_to_rle(predict))\n    \n    #----\n    probability8 = (probability * 255).astype(np.uint8)\n    plt.figure(t), plt.imshow(probability8, cmap='gray')\n    #plt.waitforbuttonpress()\n    if 'train' in CFG.mode:\n        text = metric_to_text(probability, d.label)\n        print(text)\n    print('')\n\nprint('')\nprint('CFG.mode', CFG.mode)\nsubmit_df = pd.DataFrame.from_dict(submission)\nsubmit_df.to_csv('submission.csv', index=False)\nprint(submit_df)\nprint('submission.csv ok!!!')\n\n'''\n @ read_data(): volume2   1 min 01 sec==================================\nfragment_id 2a\n\tmask (9456, 9506)\n\tvolume (9456, 9506, 12)\nCFG.stride 56\nCFG.crop_size 224\n\nH,W,len(xy) 9456 9506 27390\n @infer_one(): 426 / 427 :  5 min 09 sec\nprobability (9456, 9506)\nbce=0.20143\nth   prec   recall   fpr   dice   score\n---------------------------------------\n0.1, 0.340, 0.812, 0.182,  0.480,  0.385\n0.2, 0.470, 0.671, 0.088,  0.553,  0.500\n0.3, 0.578, 0.556, 0.047,  0.567,  0.573\n0.4, 0.681, 0.448, 0.024,  0.541,  0.617\n0.5, 0.780, 0.351, 0.011,  0.484,  0.627\n0.6, 0.862, 0.259, 0.005,  0.399,  0.588\n0.7, 0.927, 0.171, 0.002,  0.288,  0.492\n0.8, 0.970, 0.075, 0.000,  0.140,  0.287\n0.9, 0.992, 0.006, 0.000,  0.012,  0.029\n\n\nCFG.mode ['train']\n   Id                                          Predicted\n0  2a  5602838 1 5612341 5 5621847 6 5631352 7 564085...\nsubmission.csv ok!!!\n\n'''","metadata":{"execution":{"iopub.status.busy":"2023-05-14T05:30:40.234343Z","iopub.execute_input":"2023-05-14T05:30:40.235078Z","iopub.status.idle":"2023-05-14T05:42:37.077881Z","shell.execute_reply.started":"2023-05-14T05:30:40.235038Z","shell.execute_reply":"2023-05-14T05:42:37.076887Z"},"trusted":true},"execution_count":null,"outputs":[]}]}