notebooks/resnet18d-inference.ipynb

{"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":"## Introduction\nThis notebook provides the code for inference on the 1st Herculaneum fragment based on a ResNet18 model trained on fragment 2 & 3.\n\nResnet Training notebook:\n1. [Vesuvius Challenge - 3D ResNet Training](https://www.kaggle.com/code/samfc10/vesuvius-challenge-3d-resnet-training)\n","metadata":{}},{"cell_type":"markdown","source":"#### Base parameter:\n1. Do inference on pre-trained ResNet18 model\n2. 1 fold cross validation (Use 2,3 to train and 1 to val)\n3. Inference on 192 x 192 x 16 windows","metadata":{}},{"cell_type":"markdown","source":"## Setup","metadata":{}},{"cell_type":"code","source":"import os,cv2\nimport gc\nimport sys\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n# from torch.cuda import amp\nfrom torch.utils.data import Dataset, DataLoader\nimport PIL.Image as Image\n\nsys.path.append(\"/kaggle/input/resnet3d\")\nfrom resnet3d import generate_model\nimport torch as tc","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:11:54.52944Z","iopub.execute_input":"2023-06-26T11:11:54.53075Z","iopub.status.idle":"2023-06-26T11:11:54.949838Z","shell.execute_reply.started":"2023-06-26T11:11:54.530693Z","shell.execute_reply":"2023-06-26T11:11:54.948751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Configuration","metadata":{}},{"cell_type":"code","source":"class CFG:\n    # ============== paths =============\n    comp_dir_path = '/kaggle/input/'\n    comp_folder_name = 'vesuvius-challenge-ink-detection'\n    comp_dataset_path = f'{comp_dir_path}{comp_folder_name}/'\n\n    # ============== training config ========\n    in_chans = 16   # The number of layers of the papyrus you want to read at both side\n    prd_size= 192   # size of the crops \n    stride = prd_size // 8     # stride = 32\n    batch_size = 24 #32\n    seed = 42\n    num_workers=2","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:12:02.670162Z","iopub.execute_input":"2023-06-26T11:12:02.671195Z","iopub.status.idle":"2023-06-26T11:12:02.679398Z","shell.execute_reply.started":"2023-06-26T11:12:02.671154Z","shell.execute_reply":"2023-06-26T11:12:02.678108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Create test dataset","metadata":{}},{"cell_type":"code","source":"def read_image(fragment_id):\n    \"\"\"\n    return the 16 middle layers as a numpy array\"\"\"\n    images = []\n    mid = 65 // 2\n    start = mid - CFG.in_chans // 2\n    end = mid + CFG.in_chans // 2\n    idxs = range(start, end)\n\n    for i in tqdm(idxs):\n        #the tif files 1 by 1\n        image = cv2.imread(CFG.comp_dataset_path + f\"{mode}/{fragment_id}/surface_volume/{i:02}.tif\", 0)\n\n        pad0 = (CFG.prd_size - image.shape[0] % CFG.prd_size)  # %: caluclates the rest according to a size for padding\n        pad1 = (CFG.prd_size - image.shape[1] % CFG.prd_size)\n\n        image = np.pad(image, [(0, pad0), (0, pad1)], constant_values=0)\n\n        images.append(image)\n    images = np.stack(images, axis=2)  # Stack all the images along new axis \n    \n    return images","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:12:26.881072Z","iopub.execute_input":"2023-06-26T11:12:26.881486Z","iopub.status.idle":"2023-06-26T11:12:26.891801Z","shell.execute_reply.started":"2023-06-26T11:12:26.881444Z","shell.execute_reply":"2023-06-26T11:12:26.890763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class CustomDataset(Dataset):\n    def __init__(self, images, cfg, xys, labels=None):\n        self.images = images\n        self.cfg = cfg\n        self.labels = labels\n        self.xys=xys\n\n    def __len__(self):\n        # return len(self.xyxys)\n        return len(self.images)\n\n    def __getitem__(self, idx):\n        # x1, y1, x2, y2 = self.xyxys[idx]\n        image = self.images[idx]\n        image=tc.from_numpy(image).permute(2,0,1).to(tc.float32)/255\n        image = (image - 0.45)/0.225\n        return image,self.xys[idx]\n","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:12:32.379995Z","iopub.execute_input":"2023-06-26T11:12:32.380612Z","iopub.status.idle":"2023-06-26T11:12:32.388513Z","shell.execute_reply.started":"2023-06-26T11:12:32.380573Z","shell.execute_reply":"2023-06-26T11:12:32.387429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_test_dataset(fragment_id):\n    test_images = read_image(fragment_id)\n    \n    # create lists of pixels for test dataset\n    x1_list = list(range(0, test_images.shape[1]-CFG.prd_size+1, CFG.stride))  # define the stride of the image\n    y1_list = list(range(0, test_images.shape[0]-CFG.prd_size+1, CFG.stride))\n    \n    test_images_list = []\n    # list of all the x and y values of the test dataset\n    xyxys = []\n    for y1 in y1_list:\n        for x1 in x1_list:\n            #define the crops\n            y2 = y1 + CFG.prd_size\n            x2 = x1 + CFG.prd_size\n            if np.all(test_images[y1:y2, x1:x2]==0):\n                # avoid feeding zero images\n                continue\n            test_images_list.append(test_images[y1:y2, x1:x2])\n            xyxys.append((x1, y1, x2, y2))\n    xyxys = np.stack(xyxys)\n            \n    test_dataset = CustomDataset(test_images_list, CFG,xys=xyxys)\n    \n    test_loader = DataLoader(test_dataset,\n                          batch_size=CFG.batch_size,\n                          shuffle=False,\n                          num_workers=CFG.num_workers, pin_memory=True, drop_last=False)\n    \n    return test_loader, xyxys","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:12:32.903442Z","iopub.execute_input":"2023-06-26T11:12:32.903841Z","iopub.status.idle":"2023-06-26T11:12:32.915931Z","shell.execute_reply.started":"2023-06-26T11:12:32.903807Z","shell.execute_reply":"2023-06-26T11:12:32.914703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Build the Resnet18 model","metadata":{}},{"cell_type":"code","source":"class Decoder(nn.Module):\n    def __init__(self, encoder_dims, upscale):\n        super().__init__()\n        self.convs = nn.ModuleList([\n            nn.Sequential(\n                nn.Conv2d(encoder_dims[i]+encoder_dims[i-1], encoder_dims[i-1], 3, 1, 1, bias=False),\n                nn.BatchNorm2d(encoder_dims[i-1]),\n                nn.ReLU(inplace=True)\n            ) for i in range(1, len(encoder_dims))])\n\n        self.logit = nn.Conv2d(encoder_dims[0], 1, 1, 1, 0)\n        self.up = nn.Upsample(scale_factor=upscale, mode=\"bilinear\")\n\n    def forward(self, feature_maps):\n        for i in range(len(feature_maps)-1, 0, -1):\n            f_up = F.interpolate(feature_maps[i], scale_factor=2, mode=\"bilinear\")\n            f = torch.cat([feature_maps[i-1], f_up], dim=1)\n            f_down = self.convs[i-1](f)\n            feature_maps[i-1] = f_down\n\n        x = self.logit(feature_maps[0])\n        mask = self.up(x)\n        return mask\n\n\nclass SegModel(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.encoder = generate_model(model_depth=18, n_input_channels=1)\n        self.decoder = Decoder(encoder_dims=[64, 128, 256, 512], upscale=4)\n        \n    def forward(self, x):\n        if x.ndim==4:\n            x=x[:,None]\n            \n        feat_maps = self.encoder(x)\n        feat_maps_pooled = [torch.mean(f, dim=2) for f in feat_maps]\n        pred_mask = self.decoder(feat_maps_pooled)\n        return pred_mask\n    \n    def load_pretrained_weights(self, state_dict):\n        # Convert 3 channel weights to single channel\n        # ref - https://timm.fast.ai/models#Case-1:-When-the-number-of-input-channels-is-1\n        conv1_weight = state_dict['conv1.weight']\n        state_dict['conv1.weight'] = conv1_weight.sum(dim=1, keepdim=True)\n        print(self.encoder.load_state_dict(state_dict, strict=False))","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:46:02.972274Z","iopub.execute_input":"2023-06-26T11:46:02.973533Z","iopub.status.idle":"2023-06-26T11:46:02.991023Z","shell.execute_reply.started":"2023-06-26T11:46:02.973473Z","shell.execute_reply":"2023-06-26T11:46:02.989996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Load the trained weights","metadata":{}},{"cell_type":"code","source":"in_submission=True\nIS_DEBUG = False # True False\nmode = 'train' if IS_DEBUG else 'test'\nTH = 0.55  # The treshold for ink or no ink\n\n#take the test datasets\nif mode == 'test':\n    fragment_ids = sorted(os.listdir(CFG.comp_dataset_path + mode))\nelse:\n    fragment_ids = [1] \n\nmodel = SegModel()\nmodel = nn.DataParallel(model, device_ids=[0, 1])\nmodel = model.cuda()#.eval() #sending the model to a cuda device\n#model.load_state_dict(tc.load(\"/kaggle/input/3d-resnet-baseline-inference-model-data/resnet3d-34_3d_seg_epoch_14.pth\"))\nmodel.module.load_state_dict(tc.load(\"/kaggle/input/resnet18/resnet18_3d_seg_32_0.58.pt\"))\nmodel.training","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:46:13.269631Z","iopub.execute_input":"2023-06-26T11:46:13.269982Z","iopub.status.idle":"2023-06-26T11:46:14.059429Z","shell.execute_reply.started":"2023-06-26T11:46:13.269953Z","shell.execute_reply":"2023-06-26T11:46:14.058357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### To Do: Add better visualizations about the model","metadata":{}},{"cell_type":"markdown","source":"### Helper functions for inference","metadata":{}},{"cell_type":"markdown","source":"#### 1. Run length encoding of the predictions\nThe rle (run lenght encoding) function is provided to submit the predictions the ofrmat of competition requirements","metadata":{}},{"cell_type":"code","source":"\n# ref.: https://www.kaggle.com/stainsby/fast-tested-rle\ndef rle(img):\n    '''\n    img: numpy array, 1 = mask, 0 = background\n    Returns run length encoding as string formated\n    '''\n    ## DECIDE OPTIMAL threshold\n    #thr = 0.5\n    \n    pixels = img.flatten()\n    #pixels = (pixels >= thr).astype(int)\n    pixels = np.concatenate([[0], pixels, [0]])\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 1  # returns indices where consecutive elements are different\n    print('runs', runs.shape)\n    \n    runs[1::2] -= runs[::2]  #subtract the values of the even indices from the corresponding one of the odd indices and write the result \n    return ' '.join(str(x) for x in runs)","metadata":{"execution":{"iopub.status.busy":"2023-06-26T14:46:59.714563Z","iopub.execute_input":"2023-06-26T14:46:59.715302Z","iopub.status.idle":"2023-06-26T14:46:59.723349Z","shell.execute_reply.started":"2023-06-26T14:46:59.715262Z","shell.execute_reply":"2023-06-26T14:46:59.722184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 2. TTA","metadata":{}},{"cell_type":"code","source":"def TTA(x:tc.Tensor,model:nn.Module):\n    #x.shape=(batch,c,h,w)\n    shape=x.shape\n    #print(\"shape of the tensor that is predicted\" ,shape)\n    # rotate the tensor by (k*90 degrees) from the -2 to -1 axis \n    x=[x,*[tc.rot90(x,k=i,dims=(-2,-1)) for i in range(1,4)]]\n    \n    #print(\"shape of the augmented tensor that is predicted\" ,x[0].shape)\n    x=tc.cat(x,dim=0)  #concatenate in given dimension \n    \n    #print(\"shape after concatenation\" ,x[0].shape)\n    \n    # make a prediction\n    x=model(x)\n    x=torch.sigmoid(x) # create output between 0 and 1\n    #rerotate the augmented data in the original position ? \n    x=x.reshape(4,shape[0],*shape[2:])\n    x=[tc.rot90(x[i],k=-i,dims=(-2,-1)) for i in range(4)]\n    x=tc.stack(x,dim=0)\n    return x.mean(0)","metadata":{"execution":{"iopub.status.busy":"2023-06-26T11:49:21.421939Z","iopub.execute_input":"2023-06-26T11:49:21.422793Z","iopub.status.idle":"2023-06-26T11:49:21.433097Z","shell.execute_reply.started":"2023-06-26T11:49:21.422748Z","shell.execute_reply":"2023-06-26T11:49:21.43175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3. L1/ Hessian Denoising of the predictions\nThis idea was proposed during the competition by Brett Olsen:\\\n[Improving performance with L1/Hessian denoising](https://www.kaggle.com/code/brettolsen/improving-performance-with-l1-hessian-denoising)\\\n\nThe approach is to exploit known properties of the ink distribution, in particular that:\n- The **ink is sparse**, most regions will not contain ink. An **L1 regularization term** will penalyze noisy data\n- The **ink** is continuous, a single pixel is more likely to contain ink when it is newt to another pixel containing ink. A **Hessian matrix term** is used to penalize strongly variable values of ink/no-ink ","metadata":{}},{"cell_type":"code","source":"import cupy as cp\nxp = cp\n\ndelta_lookup = {\n    \"xx\": xp.array([[1, -2, 1]], dtype=float),\n    \"yy\": xp.array([[1], [-2], [1]], dtype=float),\n    \"xy\": xp.array([[1, -1], [-1, 1]], dtype=float),\n}\n\ndef operate_derivative(img_shape, pair):\n    assert len(img_shape) == 2\n    delta = delta_lookup[pair]\n    fft = xp.fft.fftn(delta, img_shape)\n    return fft * xp.conj(fft)\n\ndef soft_threshold(vector, threshold):\n    return xp.sign(vector) * xp.maximum(xp.abs(vector) - threshold, 0)\n\ndef back_diff(input_image, dim):\n    assert dim in (0, 1)\n    r, n = xp.shape(input_image)\n    size = xp.array((r, n))\n    position = xp.zeros(2, dtype=int)\n    temp1 = xp.zeros((r+1, n+1), dtype=float)\n    temp2 = xp.zeros((r+1, n+1), dtype=float)\n    \n    temp1[position[0]:size[0], position[1]:size[1]] = input_image\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    \n    size[dim] += 1\n    position[dim] += 1\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    temp1 -= temp2\n    size[dim] -= 1\n    return temp1[0:size[0], 0:size[1]]\n\ndef forward_diff(input_image, dim):\n    assert dim in (0, 1)\n    r, n = xp.shape(input_image)\n    size = xp.array((r, n))\n    position = xp.zeros(2, dtype=int)\n    temp1 = xp.zeros((r+1, n+1), dtype=float)\n    temp2 = xp.zeros((r+1, n+1), dtype=float)\n        \n    size[dim] += 1\n    position[dim] += 1\n\n    temp1[position[0]:size[0], position[1]:size[1]] = input_image\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    \n    size[dim] -= 1\n    temp2[0:size[0], 0:size[1]] = input_image\n    temp1 -= temp2\n    size[dim] += 1\n    return -temp1[position[0]:size[0], position[1]:size[1]]\n\ndef iter_deriv(input_image, b, scale, mu, dim1, dim2):\n    g = back_diff(forward_diff(input_image, dim1), dim2)\n    d = soft_threshold(g + b, 1 / mu)\n    b = b + (g - d)\n    L = scale * back_diff(forward_diff(d - b, dim2), dim1)\n    return L, b\n\ndef iter_xx(*args):\n    return iter_deriv(*args, dim1=1, dim2=1)\n\ndef iter_yy(*args):\n    return iter_deriv(*args, dim1=0, dim2=0)\n\ndef iter_xy(*args):\n    return iter_deriv(*args, dim1=0, dim2=1)\n\ndef iter_sparse(input_image, bsparse, scale, mu):\n    d = soft_threshold(input_image + bsparse, 1 / mu)\n    bsparse = bsparse + (input_image - d)\n    Lsparse = scale * (d - bsparse)\n    return Lsparse, bsparse\n\ndef denoise_image(input_image, iter_num=100, fidelity=150, sparsity_scale=10, continuity_scale=0.5, mu=1):\n    image_size = xp.shape(input_image)\n    #print(\"Initialize denoising\")\n    norm_array = (\n        operate_derivative(image_size, \"xx\") + \n        operate_derivative(image_size, \"yy\") + \n        2 * operate_derivative(image_size, \"xy\")\n    )\n    norm_array += (fidelity / mu) + sparsity_scale ** 2\n    b_arrays = {\n        \"xx\": xp.zeros(image_size, dtype=float),\n        \"yy\": xp.zeros(image_size, dtype=float),\n        \"xy\": xp.zeros(image_size, dtype=float),\n        \"L1\": xp.zeros(image_size, dtype=float),\n    }\n    g_update = xp.multiply(fidelity / mu, input_image)\n    for i in tqdm(range(iter_num), total=iter_num):\n        #print(f\"Starting iteration {i+1}\")\n        g_update = xp.fft.fftn(g_update)\n        if i == 0:\n            g = xp.fft.ifftn(g_update / (fidelity / mu)).real\n        else:\n            g = xp.fft.ifftn(xp.divide(g_update, norm_array)).real\n        g_update = xp.multiply((fidelity / mu), input_image)\n        \n        #print(\"XX update\")\n        L, b_arrays[\"xx\"] = iter_xx(g, b_arrays[\"xx\"], continuity_scale, mu)\n        g_update += L\n        \n        #print(\"YY update\")\n        L, b_arrays[\"yy\"] = iter_yy(g, b_arrays[\"yy\"], continuity_scale, mu)\n        g_update += L\n        \n        #print(\"XY update\")\n        L, b_arrays[\"xy\"] = iter_xy(g, b_arrays[\"xy\"], 2 * continuity_scale, mu)\n        g_update += L\n        \n        #print(\"L1 update\")\n        L, b_arrays[\"L1\"] = iter_sparse(g, b_arrays[\"L1\"], sparsity_scale, mu)\n        g_update += L\n        \n    g_update = xp.fft.fftn(g_update)\n    g = xp.fft.ifftn(xp.divide(g_update, norm_array)).real\n    \n    g[g < 0] = 0\n    g -= g.min()\n    g /= g.max()\n    return g","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"cell_type":"code","source":"results = []\nfor fragment_id in fragment_ids:\n    print(\"Start Inference on Fragment\", fragment_id)\n    if not in_submission:\n        break\n    \n    print(\"Load test dataset\")\n    test_loader, xyxys = make_test_dataset(fragment_id)\n    print(\"Total number of batches:\",len(test_loader))\n    \n    #mask that says where the fragment is \n    binary_mask = cv2.imread(CFG.comp_dataset_path + f\"{mode}/{fragment_id}/mask.png\", 0)\n    \n    #cv2.imshow(\"binary_mask\", binary_mask)\n    binary_mask = (binary_mask / 255).astype(int) # --> Change to 0 or 1\n    \n    ori_h = binary_mask.shape[0]\n    ori_w = binary_mask.shape[1]\n\n    # add the padding\n    pad0 = (CFG.prd_size - binary_mask.shape[0] % CFG.prd_size)\n    pad1 = (CFG.prd_size - binary_mask.shape[1] % CFG.prd_size)\n\n    # pad the binary mask\n    binary_mask = np.pad(binary_mask, [(0, pad0), (0, pad1)], constant_values=0)\n    \n    # create predictions\n    print(\"start predictions\")\n    mask_pred = np.zeros(binary_mask.shape)\n    mask_count = np.zeros(binary_mask.shape)\n    for step, (images,xys) in tqdm(enumerate(test_loader), total=len(test_loader)):\n        # 1 test loader = batch of 24 image crops --> in total: + than 1000 batches!\n        images = images.cuda()\n        batch_size = images.size(0)\n\n        with torch.no_grad():  \n            # Do the inference\n            y_preds=TTA(images,model)\n            #print(y_preds)\n            \n        \n        for k, (x1, y1, x2, y2) in enumerate(xys):\n            # Add the inference to the mask\n            mask_pred[y1:y2, x1:x2] += y_preds[k].squeeze(0).cpu().numpy()\n            mask_count[y1:y2, x1:x2] += 1\n        \n    mask_pred /= (mask_count+1e-7)\n    \n    ##### NICE PLOTS #####\n\n    fig, axes = plt.subplots(1, 3, figsize=(15, 8))\n    axes[0].imshow(mask_count)\n    axes[1].imshow(mask_pred.copy())\n\n    mask_pred=xp.array(mask_pred)\n    #mask_pred=denoise_image(mask_pred, iter_num=250)\n    mask_pred=mask_pred.get()\n\n    mask_pred = mask_pred[:ori_h, :ori_w]\n    binary_mask = binary_mask[:ori_h, :ori_w]\n\n    mask_pred_1 = (mask_pred >= TH).astype(np.uint8)\n    mask_pred_1 =mask_pred_1.astype(int)\n    mask_pred_1 *= binary_mask\n\n    axes[2].imshow(mask_pred_1)\n    plt.show()\n    \n    inklabels_rle = rle(mask_pred_1)\n    results.append((fragment_id, inklabels_rle))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!apt-get -qq install -y graphviz && pip install -q pydot\n#!pip install torchviz\n#!pip install graphviz","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#make_dot(mask_pred, params = dict())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# clean ntb memory\n\ndel mask_pred, mask_count\ndel test_loader\n    \ngc.collect()\ntorch.cuda.empty_cache()\nplt.clf()\nfig.clear()\nplt.close(fig)","metadata":{"execution":{"iopub.status.busy":"2023-06-26T13:24:23.185111Z","iopub.execute_input":"2023-06-26T13:24:23.186426Z","iopub.status.idle":"2023-06-26T13:24:23.883424Z","shell.execute_reply.started":"2023-06-26T13:24:23.186383Z","shell.execute_reply":"2023-06-26T13:24:23.882314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Add in trainign ntb !","metadata":{}},{"cell_type":"code","source":"# get metrics if in training fragments\ndef metric_to_text(ink, label):\n    text = []\n\n    #calculate bce\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.25, 0.30, 0.35, 0.4,0.45, 0.5,0.55, 0.6,0.65, 0.7]:\n        p = ink.reshape(-1) # reshape predictions\n        t = label.reshape(-1)  #rechape true labels\n        p = (p > threshold).astype(np.float32)  # make predictions as a 0 or 1\n        t = (t > 0.5).astype(np.float32)\n\n        tp = p * t # true positives\n        precision = tp.sum() / (p.sum() + 0.0001)\n        recall = tp.sum() / t.sum()\n\n        fp = p * (1 - t)  # false positives\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.2f}, {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 load_labels(fragment_id):\n    img = Image.open(f\"{CFG.comp_dataset_path}/{mode}/{fragment_id}/inklabels.png\")\n    return np.array(img)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if mode == \"train\":\n    labels = load_labels(\"1\")\n    pad0 = (CFG.prd_size - ori_h % CFG.prd_size)\n    pad1 = (CFG.prd_size - ori_w % CFG.prd_size)\n\n    labels = np.pad(labels, [(0, pad0), (0, pad1)], constant_values=0)\n    print(\"true labels shape\", labels.shape)\n    print(\"predicted mask shape\", mask_pred.shape)\n    \n    text = metric_to_text(mask_pred, labels)\n    print(text)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## submission","metadata":{}},{"cell_type":"code","source":"! cp /kaggle/input/vesuvius-challenge-ink-detection/sample_submission.csv submission.csv\nif in_submission:\n    sub = pd.DataFrame(results, columns=['Id', 'Predicted'])\n    #sub\n    sample_sub = pd.read_csv(CFG.comp_dataset_path + 'sample_submission.csv')\n    sample_sub = pd.merge(sample_sub[['Id']], sub, on='Id', how='left')\n    #sample_sub\n    sample_sub.to_csv(\"submission.csv\", index=False)\n    print(\"ok\")","metadata":{"execution":{"iopub.status.busy":"2023-06-26T14:12:18.949494Z","iopub.execute_input":"2023-06-26T14:12:18.95066Z","iopub.status.idle":"2023-06-26T14:12:20.565755Z","shell.execute_reply.started":"2023-06-26T14:12:18.95061Z","shell.execute_reply":"2023-06-26T14:12:20.564413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}