{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":117682,"databundleVersionId":14443416,"sourceType":"competition"},{"sourceId":13698445,"sourceType":"datasetVersion","datasetId":8670484},{"sourceId":13874976,"sourceType":"datasetVersion","datasetId":8840061},{"sourceId":14223014,"sourceType":"datasetVersion","datasetId":8790422},{"sourceId":14230131,"sourceType":"datasetVersion","datasetId":9063178},{"sourceId":262076264,"sourceType":"kernelVersion"},{"sourceId":284699010,"sourceType":"kernelVersion"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"pip install tensorly-torch","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:32:59.385044Z","iopub.execute_input":"2025-12-21T04:32:59.385335Z","iopub.status.idle":"2025-12-21T04:33:03.946648Z","shell.execute_reply.started":"2025-12-21T04:32:59.385313Z","shell.execute_reply":"2025-12-21T04:33:03.945978Z"},"_kg_hide-output":true,"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pip install tensorly","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:33:03.948134Z","iopub.execute_input":"2025-12-21T04:33:03.948364Z","iopub.status.idle":"2025-12-21T04:33:07.647886Z","shell.execute_reply.started":"2025-12-21T04:33:03.948341Z","shell.execute_reply":"2025-12-21T04:33:07.647019Z"},"_kg_hide-input":true,"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nVesuvius competition metric.\n\nExpects standard Kaggle paths and Linux in order to manage dependencies.\n\"\"\"\n\nimport glob\nimport importlib\nimport os\nimport subprocess\nimport sys\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image, ImageSequence\nfrom scipy.ndimage import distance_transform_edt, binary_dilation, grey_dilation, gaussian_gradient_magnitude, map_coordinates, gaussian_filter, binary_fill_holes\nimport networkx as nx\nfrom scipy.sparse import lil_matrix\nfrom scipy.sparse.linalg import lsqr\nfrom scipy.sparse import csr_matrix\nfrom scipy.sparse.linalg import spsolve\nfrom pathlib import Path\nfrom scipy.sparse import coo_matrix\nfrom scipy.sparse.linalg import cg\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nimport torch\nimport torch.nn.functional as F","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:33:07.649003Z","iopub.execute_input":"2025-12-21T04:33:07.649578Z","iopub.status.idle":"2025-12-21T04:33:11.614439Z","shell.execute_reply.started":"2025-12-21T04:33:07.649552Z","shell.execute_reply":"2025-12-21T04:33:11.613841Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class ParticipantVisibleError(Exception):\n    pass\n\n\nclass HostVisibleError(Exception):\n    pass\n\n\ndef load_volume(path):\n    im = Image.open(path)\n    slices = []\n    for i, page in enumerate(ImageSequence.Iterator(im)):\n        slice_array = np.array(page)\n        slices.append(slice_array)\n    volume = np.stack(slices, axis=0)\n    return volume\n\n\ndef install_dependencies():\n    \"\"\"On Kaggle, the topometrics library must be installed during the run. This function handles the entire process.\"\"\"\n    try:\n        import topometrics.leaderboard\n\n        return None\n    # The broad exception is necessary as the initial import can fail for multiple reasons.\n    except:\n        pass\n\n    resources_dir = '/kaggle/input/vesuvius-metric-resources'\n    install_dir = '/kaggle/working/topological-metrics-kaggle'\n\n    try:\n        subprocess.run(\n            f'cd {resources_dir} && uv pip install --no-index --find-links=wheels -r topological-metrics-kaggle/requirements.txt',\n            shell=True,\n            check=True,\n        )\n        subprocess.run(f'cd /kaggle/working && cp -r {resources_dir}/topological-metrics-kaggle .', shell=True, check=True)\n        subprocess.run(\n            f'cd {install_dir} && chmod +x scripts/setup_submodules.sh scripts/build_betti.sh && make build-betti',\n            shell=True,\n            check=True,\n        )\n        subprocess.run(\n            f'cd {install_dir} && uv pip install -e . --no-deps --no-index --no-build-isolation -v',\n            shell=True,\n            check=True,\n        )\n        # Add the new library to Python's path and invalidate caches to ensure it's found.\n        sys.path.append('/kaggle/working/topological-metrics-kaggle/src')\n        importlib.invalidate_caches()\n\n    except Exception as err:\n        raise HostVisibleError(f'Failed to install topometrics library: {err}')\n\n\ndef generate_standard_submission(submission_dir: str) -> None:\n    # Dependencies installed here as generate_standard_submission is the first metric function that gets called by the orchestrator.\n    submission_tifs = glob.glob(f'{submission_dir}/**/*.tif', recursive=True)\n    if len(submission_tifs) == 0:\n        submission_tifs = glob.glob('/kaggle/tmp/**/*.tif', recursive=True)\n    if len(submission_tifs) == 0:\n        raise ParticipantVisibleError('No submission files found')\n    df = pd.DataFrame({'tif_paths': submission_tifs})\n    df['id'] = df['tif_paths'].apply(lambda x: x.split('/')[-1].split('.')[0])\n    os.chdir('/kaggle/working')\n    df[['id', 'tif_paths']].to_csv('submission.csv', index=False)\n\n\ndef score_single_tif(\n    gt,\n    pr,\n    surface_tolerance = 2.0,\n    voi_connectivity=26,\n    voi_transform='one_over_one_plus',\n    voi_alpha=0.3,\n    topo_weight=0.3,\n    surface_dice_weight=0.35,\n    voi_weight=0.35,\n):\n    # The import is here to ensure dependencies are loaded first.\n    try:\n        # Use a standard import now that the path is reliably set.\n        import topometrics.leaderboard\n    except Exception as err:\n        raise HostVisibleError(f'Failed to import topometrics after installation: {err}')\n\n    score_report = topometrics.leaderboard.compute_leaderboard_score(\n        predictions=pr,\n        labels=gt,\n        dims=(0, 1, 2),\n        spacing=(1.0, 1.0, 1.0),  # (z, y, x)\n        surface_tolerance=surface_tolerance,  # in spacing units\n        voi_connectivity=voi_connectivity,\n        voi_transform=voi_transform,\n        voi_alpha=voi_alpha,\n        combine_weights=(topo_weight, surface_dice_weight, voi_weight),  # (Topo, SurfaceDice, VOI)\n        fg_threshold=None,  # None => legacy \"!= 0\"; else uses \"x > threshold\"\n        ignore_label=2,  # voxels with this GT label are ignored\n        ignore_mask=None,  # or pass an explicit boolean mask\n    )\n    return score_report\n\ninstall_dependencies()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:33:11.615237Z","iopub.execute_input":"2025-12-21T04:33:11.615708Z","iopub.status.idle":"2025-12-21T04:34:05.100284Z","shell.execute_reply.started":"2025-12-21T04:33:11.615678Z","shell.execute_reply":"2025-12-21T04:34:05.099601Z"},"_kg_hide-input":true,"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data_path = Path('/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1407735.tif')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:05.101968Z","iopub.execute_input":"2025-12-21T04:34:05.102379Z","iopub.status.idle":"2025-12-21T04:34:05.106054Z","shell.execute_reply.started":"2025-12-21T04:34:05.102361Z","shell.execute_reply":"2025-12-21T04:34:05.105331Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mask = load_volume(data_path)\nmask = mask * (mask != 2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:05.107166Z","iopub.execute_input":"2025-12-21T04:34:05.107433Z","iopub.status.idle":"2025-12-21T04:34:05.776510Z","shell.execute_reply.started":"2025-12-21T04:34:05.107409Z","shell.execute_reply":"2025-12-21T04:34:05.775713Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Decomposition A: Thickness, SDF, Normals (dx, dy, dz)","metadata":{}},{"cell_type":"code","source":"class ExactScrollDecomposition:\n    \"\"\"\n    Exact scroll mask decomposition using:\n    - sdf: signed distance field\n    - normals: gradient of sdf\n    - thickness: distance to mask boundary along normal\n    Fully vectorized inverse transform using sparse accumulation.\n    \"\"\"\n\n    def __init__(self, sdf_truncate=20.0, substeps=5):\n        self.sdf_truncate = sdf_truncate\n        self.substeps = substeps  # sub-voxel sampling along normal\n\n    def transform(self, mask):\n        mask = mask.astype(bool)\n        self.mask_shape = mask.shape\n        sdf = self._compute_sdf(mask)\n        normals = self._normals_from_sdf(sdf)\n        thickness = self._compute_thickness(mask)\n        return {\"sdf\": sdf, \"normals\": normals, \"thickness\": thickness}\n\n    def inverse_transform(self, pred):\n        D, H, W = self.mask_shape\n        sdf = pred[\"sdf\"]\n        normals = pred[\"normals\"]\n        thickness = pred[\"thickness\"]\n\n        # --- Center surface coordinates ---\n        center_mask = sdf <= 0\n        coords = np.argwhere(center_mask)\n        if coords.size == 0:\n            return np.zeros((D, H, W), dtype=np.uint8)\n\n        dz = normals[0][center_mask]\n        dy = normals[1][center_mask]\n        dx = normals[2][center_mask]\n        t_half = thickness[center_mask] / 2.0\n\n        # --- Sub-voxel steps along normals ---\n        num_centers = coords.shape[0]\n        steps = np.linspace(-1, 1, self.substeps)[None, :] * t_half[:, None]\n\n        z_coords = coords[:, 0:1] + dz[:, None] * steps\n        y_coords = coords[:, 1:2] + dy[:, None] * steps\n        x_coords = coords[:, 2:3] + dx[:, None] * steps\n\n        # --- Round and clip voxel indices ---\n        z_idx = np.clip(np.round(z_coords).astype(int), 0, D-1).ravel()\n        y_idx = np.clip(np.round(y_coords).astype(int), 0, H-1).ravel()\n        x_idx = np.clip(np.round(x_coords).astype(int), 0, W-1).ravel()\n\n        # --- Sparse accumulation ---\n        flat_idx = z_idx*H*W + y_idx*W + x_idx\n        ones = np.ones_like(flat_idx, dtype=np.uint16)\n        recon_sparse = coo_matrix((ones, (flat_idx, np.zeros_like(flat_idx))),\n                                  shape=(D*H*W, 1))\n\n        # Convert to dense mask\n        recon_flat = (recon_sparse.toarray().ravel() > 0).astype(np.uint8)\n        recon = recon_flat.reshape(D, H, W)\n\n        return recon\n\n    def _compute_sdf(self, mask):\n        dist_out = distance_transform_edt(~mask)\n        dist_in = distance_transform_edt(mask)\n        sdf = dist_out - dist_in\n        return np.clip(sdf, -self.sdf_truncate, self.sdf_truncate).astype(np.float32)\n\n    def _normals_from_sdf(self, sdf):\n        dz, dy, dx = np.gradient(sdf)\n        grad = np.stack([dz, dy, dx], axis=0)\n        norm = np.linalg.norm(grad, axis=0, keepdims=True) + 1e-8\n        return (grad / norm).astype(np.float32)\n\n    def _compute_thickness(self, mask):\n        dist_in = distance_transform_edt(mask)\n        dist_out = distance_transform_edt(~mask)\n        thickness = dist_in + dist_out\n        return thickness.astype(np.float32)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:05.777469Z","iopub.execute_input":"2025-12-21T04:34:05.777887Z","iopub.status.idle":"2025-12-21T04:34:05.788951Z","shell.execute_reply.started":"2025-12-21T04:34:05.777863Z","shell.execute_reply":"2025-12-21T04:34:05.788181Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"decomp = ExactScrollDecomposition()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:05.789688Z","iopub.execute_input":"2025-12-21T04:34:05.789933Z","iopub.status.idle":"2025-12-21T04:34:05.805746Z","shell.execute_reply.started":"2025-12-21T04:34:05.789916Z","shell.execute_reply":"2025-12-21T04:34:05.805065Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pred = decomp.transform(mask)\nrecon = decomp.inverse_transform(pred)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:05.806468Z","iopub.execute_input":"2025-12-21T04:34:05.806742Z","iopub.status.idle":"2025-12-21T04:34:42.584100Z","shell.execute_reply.started":"2025-12-21T04:34:05.806721Z","shell.execute_reply":"2025-12-21T04:34:42.583242Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(recon[100])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:42.584954Z","iopub.execute_input":"2025-12-21T04:34:42.585218Z","iopub.status.idle":"2025-12-21T04:34:42.776731Z","shell.execute_reply.started":"2025-12-21T04:34:42.585197Z","shell.execute_reply":"2025-12-21T04:34:42.776139Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots(ncols = 5, figsize = (24, 24))\nax[0].imshow(pred['thickness'][100])\nax[0].set_title('Thinkness')\nax[1].imshow(pred['sdf'][100])\nax[1].set_title('sign distance function')\nax[2].imshow(pred['normals'][0, 100])\nax[2].set_title('normals dx')\nax[3].imshow(pred['normals'][1, 100])\nax[3].set_title('normals dy')\nax[4].imshow(pred['normals'][2, 100])\nax[4].set_title('normals dz')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:42.777433Z","iopub.execute_input":"2025-12-21T04:34:42.777780Z","iopub.status.idle":"2025-12-21T04:34:43.748393Z","shell.execute_reply.started":"2025-12-21T04:34:42.777762Z","shell.execute_reply":"2025-12-21T04:34:43.747543Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Mask voxels:\", mask.sum())\nprint(\"Recon voxels:\", recon.sum())\nprint(\"Reconstruction error:\", np.sum(recon != mask))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:43.749274Z","iopub.execute_input":"2025-12-21T04:34:43.749493Z","iopub.status.idle":"2025-12-21T04:34:43.837555Z","shell.execute_reply.started":"2025-12-21T04:34:43.749474Z","shell.execute_reply":"2025-12-21T04:34:43.836961Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"score_report = score_single_tif(mask, recon)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:34:43.838229Z","iopub.execute_input":"2025-12-21T04:34:43.838403Z","iopub.status.idle":"2025-12-21T04:36:05.553682Z","shell.execute_reply.started":"2025-12-21T04:34:43.838389Z","shell.execute_reply":"2025-12-21T04:36:05.552895Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\n    \"total_score:\", score_report.score,\n    \"topo_score:\", score_report.topo.toposcore,\n    \"voi_score:\", score_report.voi.voi_score,\n    \"surface_dice:\", score_report.surface_dice)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:36:05.556529Z","iopub.execute_input":"2025-12-21T04:36:05.557167Z","iopub.status.idle":"2025-12-21T04:36:05.561228Z","shell.execute_reply.started":"2025-12-21T04:36:05.557147Z","shell.execute_reply":"2025-12-21T04:36:05.560565Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Decomposition B: Convolutional Dictionary Learning","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import DictionaryLearning, MiniBatchDictionaryLearning\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport random\nfrom tqdm import tqdm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:36:05.561918Z","iopub.execute_input":"2025-12-21T04:36:05.562193Z","iopub.status.idle":"2025-12-21T04:36:06.182613Z","shell.execute_reply.started":"2025-12-21T04:36:05.562168Z","shell.execute_reply":"2025-12-21T04:36:06.181877Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def cdl_sdt_with_betti1_ripser_deterministic(\n    binary_mask_np,\n    device=\"cuda\",\n    K=64,\n    kernel_size=9,\n    sdt_clip=20.0,\n    lambda_sparse=0.01,\n    lambda_betti1=0.1,\n    n_epochs=200,\n    seed=42,\n):\n    \"\"\"\n    Fully deterministic 3D Convolutional Dictionary Learning on SDT\n    with β1 persistent homology loss (CubicalRipser).\n\n    Args:\n        binary_mask_np (np.ndarray): input mask (D,H,W) with 0/1\n        device (str): \"cuda\" or \"cpu\"\n        K (int): number of convolutional atoms\n        kernel_size (int): size of 3D convolutional atoms\n        sdt_clip (float): clip SDT values\n        lambda_sparse (float): weight for L1 sparsity\n        lambda_betti1 (float): weight for PH β1 loss\n        n_epochs (int): number of training epochs\n        seed (int): deterministic seed\n\n    Returns:\n        dict_atoms (torch.Tensor): (K,1,k,k,k)\n        sparse_codes (torch.Tensor): (1,K,D,H,W)\n        reconstructed_mask (torch.Tensor): (D,H,W) binary\n    \"\"\"\n\n    binary_mask_np = binary_mask_np.astype('float32')\n\n    # -------------------------\n    # 0. Set deterministic seeds\n    # -------------------------\n    np.random.seed(seed)\n    random.seed(seed)\n    torch.manual_seed(seed)\n    if device.startswith(\"cuda\"):\n        torch.cuda.manual_seed(seed)\n        torch.cuda.manual_seed_all(seed)\n        torch.backends.cudnn.deterministic = True\n        torch.backends.cudnn.benchmark = False\n\n    # -------------------------\n    # 1. SDT preprocessing\n    # -------------------------\n    sdt = distance_transform_edt(binary_mask_np) - distance_transform_edt(1 - binary_mask_np)\n    sdt = np.clip(sdt, -sdt_clip, sdt_clip)\n    X = torch.tensor(sdt, dtype=torch.float32, device=device)[None, None]  # (1,1,D,H,W)\n\n    # -------------------------\n    # 2. CDL dictionary & sparse codes\n    # -------------------------\n    class ConvDictionary(nn.Module):\n        def __init__(self, K, k, seed):\n            super().__init__()\n            rng = torch.Generator(device=device).manual_seed(seed)\n            self.D = nn.Parameter(torch.randn(K, 1, k, k, k, device=device, generator=rng))\n            self.normalize()\n\n        def normalize(self):\n            with torch.no_grad():\n                self.D /= (torch.norm(self.D.flatten(1), dim=1).view(-1,1,1,1,1) + 1e-8)\n\n    dict_model = ConvDictionary(K, kernel_size, seed)\n\n    # deterministic sparse code initialization\n    rng_z = torch.Generator(device=device).manual_seed(seed)\n    Z = nn.Parameter(torch.randn(1, K, *X.shape[2:], device=device, generator=rng_z))\n\n    opt_D = torch.optim.Adam([dict_model.D], lr=1e-3)\n    opt_Z = torch.optim.Adam([Z], lr=1e-2)\n\n    pad = kernel_size // 2\n\n    # -------------------------\n    # 3. Training loop\n    # -------------------------\n    for epoch in tqdm(range(n_epochs)):\n\n        # --- sparse coding update ---\n        for _ in range(5):\n            opt_Z.zero_grad()\n\n            # grouped conv: Z shape (1,K,D,H,W), D shape (K,1,k,k,k)\n            X_hat = F.conv3d(Z, dict_model.D, padding=pad, groups=K)\n            X_hat_sum = X_hat.sum(dim=1, keepdim=True)  # (1,1,D,H,W)\n\n            loss_recon = 0.5 * ((X_hat_sum - X) ** 2).mean()\n            loss_sparse = lambda_sparse * Z.abs().mean()\n            loss = loss_recon + loss_sparse\n            loss.backward()\n            opt_Z.step()\n\n        # --- dictionary update ---\n        opt_D.zero_grad()\n        X_hat = F.conv3d(Z, dict_model.D, padding=pad, groups=K)\n        X_hat_sum = X_hat.sum(dim=1, keepdim=True)\n        loss_D = 0.5 * ((X_hat_sum - X) ** 2).mean()\n        loss_D.backward()\n        opt_D.step()\n        dict_model.normalize()\n\n        if epoch % 20 == 0:\n            print(f\"[{epoch:03d}] Recon={loss_recon.item():.4f}, Sparse={loss_sparse.item():.4f}\")\n\n    # -------------------------\n    # 4. Reconstruction & threshold\n    # -------------------------\n    with torch.no_grad():\n        X_hat = F.conv3d(Z, dict_model.D, padding=pad, groups=K)\n        X_hat_sum = X_hat.sum(dim=1, keepdim=True)\n        mask_hat = (X_hat_sum[0,0] > 0).float()\n\n    return dict_model.D.detach(), Z.detach(), mask_hat\n\n\ndef compute_sparse_code_from_mask(mask_np, D, n_iters=100, lr=0.1, lambda_sparse=0.01, sdt_clip=20.0, device=\"cuda\"):\n    \"\"\"\n    Compute sparse code Z for a new input mask using fixed dictionary D.\n\n    Args:\n        mask_np (np.ndarray): input binary mask (D,H,W)\n        D (torch.Tensor): fixed dictionary (K,1,k,k,k)\n        n_iters (int): number of optimization steps\n        lr (float): learning rate for Z optimization\n        lambda_sparse (float): weight for sparsity regularization\n        sdt_clip (float): clip SDT values\n        device (str): \"cuda\" or \"cpu\"\n\n    Returns:\n        Z (torch.Tensor): optimized sparse code (1,K,D,H,W)\n    \"\"\"\n    # 1. Compute SDT\n    sdt = distance_transform_edt(mask_np) - distance_transform_edt(1 - mask_np)\n    sdt = np.clip(sdt, -sdt_clip, sdt_clip)\n    X = torch.tensor(sdt, dtype=torch.float32, device=device).unsqueeze(0).unsqueeze(0)  # (1,1,D,H,W)\n\n    # 2. Initialize Z\n    K = D.shape[0]\n    pad = D.shape[2] // 2  # kernel_size//2\n    Z = torch.zeros(1, K, *X.shape[2:], device=device, requires_grad=True)\n\n    optimizer = torch.optim.Adam([Z], lr=lr)\n\n    # 3. Optimize Z\n    for i in tqdm(range(n_iters), desc=\"Sparse coding\"):\n        optimizer.zero_grad()\n        X_hat = F.conv3d(Z, D, padding=pad, groups=K).sum(dim=1, keepdim=True)\n        loss = 0.5 * F.mse_loss(X_hat, X) + lambda_sparse * Z.abs().mean()\n        loss.backward()\n        optimizer.step()\n        if i % 20==0:\n            print(f\"[{i:03d}] Recon={loss.item():.4f}\")\n\n    return Z.detach()\n\n\ndef reconstruct_mask(Z, D):\n    K = D.shape[0]\n    pad = D.shape[2] // 2\n    sdf_hat = F.conv3d(Z, D, padding=pad, groups=K).sum(dim=1, keepdim=True)\n    mask_hat = new_sdf_hat>=0\n    return mask_hat","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:45:48.753921Z","iopub.execute_input":"2025-12-21T04:45:48.754471Z","iopub.status.idle":"2025-12-21T04:45:48.770142Z","shell.execute_reply.started":"2025-12-21T04:45:48.754449Z","shell.execute_reply":"2025-12-21T04:45:48.769364Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"D, Z, mask_hat = cdl_sdt_with_betti1_ripser_deterministic(mask[:64, :64, :64])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:36:06.200436Z","iopub.execute_input":"2025-12-21T04:36:06.200674Z","iopub.status.idle":"2025-12-21T04:42:19.112840Z","shell.execute_reply.started":"2025-12-21T04:36:06.200637Z","shell.execute_reply":"2025-12-21T04:42:19.112189Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"D[D.abs()<1e-8]=0.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:19.113634Z","iopub.execute_input":"2025-12-21T04:42:19.114194Z","iopub.status.idle":"2025-12-21T04:42:19.199216Z","shell.execute_reply.started":"2025-12-21T04:42:19.114174Z","shell.execute_reply":"2025-12-21T04:42:19.198411Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"D.shape, Z.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:19.200008Z","iopub.execute_input":"2025-12-21T04:42:19.200214Z","iopub.status.idle":"2025-12-21T04:42:19.205627Z","shell.execute_reply.started":"2025-12-21T04:42:19.200193Z","shell.execute_reply":"2025-12-21T04:42:19.205070Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(mask[:64, :64, :64][0])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:19.206368Z","iopub.execute_input":"2025-12-21T04:42:19.207034Z","iopub.status.idle":"2025-12-21T04:42:19.359188Z","shell.execute_reply.started":"2025-12-21T04:42:19.207017Z","shell.execute_reply":"2025-12-21T04:42:19.358643Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(mask_hat[0].cpu())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:19.359847Z","iopub.execute_input":"2025-12-21T04:42:19.360108Z","iopub.status.idle":"2025-12-21T04:42:19.500019Z","shell.execute_reply.started":"2025-12-21T04:42:19.360090Z","shell.execute_reply":"2025-12-21T04:42:19.499370Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def visualize_atoms(D, axis=0):\n    \"\"\"\n    D: torch.Tensor (K,1,k,k,k)\n    axis: which axis to project for visualization\n    \"\"\"\n    D = D.detach().cpu().numpy()\n    K = D.shape[0]\n    fig, axs = plt.subplots(1, K, figsize=(K*2,2))\n    for k in range(K):\n        mip = D[k,0].max(axis=axis)  # max projection\n        axs[k].imshow(mip, cmap='bwr')\n        axs[k].axis('off')\n    plt.show()\n\ndef visualize_atom_slices(D):\n    D = D.detach().cpu().numpy()\n    K, _, k, _, _ = D.shape\n    mid = k//2\n    fig, axs = plt.subplots(3, 3, figsize=(9,3*3))\n    for i in range(3):\n        axs[i,0].imshow(D[i,0,mid,:,:], cmap='bwr'); axs[i,0].set_title('Z-slice')\n        axs[i,1].imshow(D[i,0,:,mid,:], cmap='bwr'); axs[i,1].set_title('Y-slice')\n        axs[i,2].imshow(D[i,0,:,:,mid], cmap='bwr'); axs[i,2].set_title('X-slice')\n        for ax in axs[i]: ax.axis('off')\n    plt.show()\n\ndef visualize_Z(Z, axis=0):\n    \"\"\"\n    Z: torch.Tensor (1,K,D,H,W)\n    \"\"\"\n    Z = Z[0].detach().cpu().numpy()  # (K,D,H,W)\n    K = Z.shape[0]\n    fig, axs = plt.subplots(1,K, figsize=(K*2,2))\n    for k in range(K):\n        mip = Z[k].max(axis=axis)\n        axs[k].imshow(mip, cmap='hot')\n        axs[k].axis('off')\n    plt.show()\n\ndef overlay_Z_on_mask(Z, mask, axis=0):\n    Z = Z[0].detach().cpu().numpy()\n    mask = mask.cpu().numpy() if torch.is_tensor(mask) else mask\n    K = Z.shape[0]\n    fig, axs = plt.subplots(1,K, figsize=(K*2,2))\n    for k in range(K):\n        mip = Z[k].max(axis=axis)\n        axs[k].imshow(mask.max(axis=axis), cmap='gray', alpha=0.3)\n        axs[k].imshow(mip, cmap='hot', alpha=0.6)\n        axs[k].axis('off')\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:19.500719Z","iopub.execute_input":"2025-12-21T04:42:19.500916Z","iopub.status.idle":"2025-12-21T04:42:19.510277Z","shell.execute_reply.started":"2025-12-21T04:42:19.500901Z","shell.execute_reply":"2025-12-21T04:42:19.509698Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"visualize_atoms(D)\nvisualize_atom_slices(D)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:19.510902Z","iopub.execute_input":"2025-12-21T04:42:19.511133Z","iopub.status.idle":"2025-12-21T04:42:21.248158Z","shell.execute_reply.started":"2025-12-21T04:42:19.511109Z","shell.execute_reply":"2025-12-21T04:42:21.247549Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"new_mask = load_volume('/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1006462223.tif')\nnew_mask = new_mask * (new_mask != 2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:21.249533Z","iopub.execute_input":"2025-12-21T04:42:21.250168Z","iopub.status.idle":"2025-12-21T04:42:21.814760Z","shell.execute_reply.started":"2025-12-21T04:42:21.250148Z","shell.execute_reply":"2025-12-21T04:42:21.814104Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(new_mask[:64, :64, :64][0])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:42:21.815477Z","iopub.execute_input":"2025-12-21T04:42:21.815720Z","iopub.status.idle":"2025-12-21T04:42:21.967276Z","shell.execute_reply.started":"2025-12-21T04:42:21.815702Z","shell.execute_reply":"2025-12-21T04:42:21.966677Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"new_Z = compute_sparse_code_from_mask(new_mask[:64, :64, :64], D, n_iters=100)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:48:49.336981Z","iopub.execute_input":"2025-12-21T04:48:49.337542Z","iopub.status.idle":"2025-12-21T04:49:06.097241Z","shell.execute_reply.started":"2025-12-21T04:48:49.337517Z","shell.execute_reply":"2025-12-21T04:49:06.096503Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"new_Z[new_Z.abs()<1e-8] = 0.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:50:22.862195Z","iopub.execute_input":"2025-12-21T04:50:22.862723Z","iopub.status.idle":"2025-12-21T04:50:22.866419Z","shell.execute_reply.started":"2025-12-21T04:50:22.862702Z","shell.execute_reply":"2025-12-21T04:50:22.865806Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# (new_Z==0).sum(), new_Z.numel()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:47:34.708876Z","iopub.execute_input":"2025-12-21T04:47:34.709487Z","iopub.status.idle":"2025-12-21T04:47:34.712543Z","shell.execute_reply.started":"2025-12-21T04:47:34.709465Z","shell.execute_reply":"2025-12-21T04:47:34.711839Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"sparse rate: {round(((D==0).sum()/(D.numel())).item(), 3)}%\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:50:27.354810Z","iopub.execute_input":"2025-12-21T04:50:27.355474Z","iopub.status.idle":"2025-12-21T04:50:27.359961Z","shell.execute_reply.started":"2025-12-21T04:50:27.355448Z","shell.execute_reply":"2025-12-21T04:50:27.359223Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"K = D.shape[0]\npad = D.shape[2] // 2\nnew_sdf_hat = F.conv3d(new_Z, D, padding=pad, groups=K).sum(dim=1, keepdim=True)\nnew_mask_hat = new_sdf_hat>=0","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:50:27.765328Z","iopub.execute_input":"2025-12-21T04:50:27.766089Z","iopub.status.idle":"2025-12-21T04:50:27.770192Z","shell.execute_reply.started":"2025-12-21T04:50:27.766069Z","shell.execute_reply":"2025-12-21T04:50:27.769597Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(new_mask_hat[0, 0].cpu()[3])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:50:32.986837Z","iopub.execute_input":"2025-12-21T04:50:32.987100Z","iopub.status.idle":"2025-12-21T04:50:33.130161Z","shell.execute_reply.started":"2025-12-21T04:50:32.987081Z","shell.execute_reply":"2025-12-21T04:50:33.129579Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"score_report = score_single_tif(new_mask_hat[0, 0].cpu(), new_mask[:64, :64, :64])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:50:35.088500Z","iopub.execute_input":"2025-12-21T04:50:35.088821Z","iopub.status.idle":"2025-12-21T04:50:35.453012Z","shell.execute_reply.started":"2025-12-21T04:50:35.088798Z","shell.execute_reply":"2025-12-21T04:50:35.452429Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\n    \"total_score:\", score_report.score,\n    \"topo_score:\", score_report.topo.toposcore,\n    \"voi_score:\", score_report.voi.voi_score,\n    \"surface_dice:\", score_report.surface_dice)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-21T04:50:35.454017Z","iopub.execute_input":"2025-12-21T04:50:35.454290Z","iopub.status.idle":"2025-12-21T04:50:35.458397Z","shell.execute_reply.started":"2025-12-21T04:50:35.454266Z","shell.execute_reply":"2025-12-21T04:50:35.457777Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}