{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":84969,"databundleVersionId":10033515,"sourceType":"competition"},{"sourceId":9862305,"sourceType":"datasetVersion","datasetId":6052780},{"sourceId":9867543,"sourceType":"datasetVersion","datasetId":6040935},{"sourceId":10445850,"sourceType":"datasetVersion","datasetId":6465904},{"sourceId":10471985,"sourceType":"datasetVersion","datasetId":6484063},{"sourceId":10560529,"sourceType":"datasetVersion","datasetId":6533803},{"sourceId":206640467,"sourceType":"kernelVersion"},{"sourceId":211097053,"sourceType":"kernelVersion"}],"dockerImageVersionId":30823,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# **《《《 YOLO 》》》**","metadata":{}},{"cell_type":"code","source":"from IPython.display import clear_output\n!tar xfvz /kaggle/input/ultralytics-for-offline-install/archive.tar.gz\n!pip install --no-index --find-links=./packages ultralytics\n!rm -rf ./packages\ntry:\n    import zarr\nexcept: \n    !cp -r '/kaggle/input/hengck-czii-cryo-et-01/wheel_file' '/kaggle/working/'\n    !pip install /kaggle/working/wheel_file/asciitree-0.3.3/asciitree-0.3.3\n    !pip install --no-index --find-links=/kaggle/working/wheel_file zarr\n    !pip install --no-index --find-links=/kaggle/working/wheel_file connected-components-3d\nfrom typing import List, Tuple, Union\ndeps_path = '/kaggle/input/czii-cryoet-dependencies'\n! pip install -q --no-index --find-links {deps_path} --requirement {deps_path}/requirements.txt\nimport lightning.pytorch as pl\nfrom datetime import datetime\nimport pytz\nimport sys\nsys.path.append('/kaggle/input/hengck-czii-cryo-et-01')\nfrom czii_helper import *\nfrom dataset import *\nfrom model2 import *\nclear_output()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:46:32.069369Z","iopub.execute_input":"2025-02-05T06:46:32.069642Z","iopub.status.idle":"2025-02-05T06:47:58.063027Z","shell.execute_reply.started":"2025-02-05T06:46:32.069620Z","shell.execute_reply":"2025-02-05T06:47:58.062280Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport glob\nimport time\nimport sys\nimport warnings\nimport math\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport cv2\nimport torch\nfrom tqdm import tqdm\nfrom ultralytics import YOLO\nimport zarr\nfrom scipy.spatial import cKDTree\nfrom collections import defaultdict\nfrom sklearn.ensemble import IsolationForest","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:47:58.063882Z","iopub.execute_input":"2025-02-05T06:47:58.064213Z","iopub.status.idle":"2025-02-05T06:48:00.826521Z","shell.execute_reply.started":"2025-02-05T06:47:58.064178Z","shell.execute_reply":"2025-02-05T06:48:00.825813Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pytorch_lightning as pl\n\nfrom monai.networks.nets import UNet\nfrom monai.metrics import DiceMetric\nfrom torch.optim.lr_scheduler import (\n    CosineAnnealingWarmRestarts,\n    OneCycleLR,\n    ReduceLROnPlateau\n)\nimport time","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:00.827489Z","iopub.execute_input":"2025-02-05T06:48:00.827798Z","iopub.status.idle":"2025-02-05T06:48:12.970704Z","shell.execute_reply.started":"2025-02-05T06:48:00.827768Z","shell.execute_reply":"2025-02-05T06:48:12.969769Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"total_start = time.time()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:12.971612Z","iopub.execute_input":"2025-02-05T06:48:12.972386Z","iopub.status.idle":"2025-02-05T06:48:12.976008Z","shell.execute_reply.started":"2025-02-05T06:48:12.972357Z","shell.execute_reply":"2025-02-05T06:48:12.975165Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model_path = '/kaggle/input/czii-yolo-l-trained-with-synthetic-data/best_synthetic.pt'\nmodel = YOLO(model_path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:12.977106Z","iopub.execute_input":"2025-02-05T06:48:12.977492Z","iopub.status.idle":"2025-02-05T06:48:14.015821Z","shell.execute_reply.started":"2025-02-05T06:48:12.977456Z","shell.execute_reply":"2025-02-05T06:48:14.015087Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"runs = sorted(glob.glob('/kaggle/input/czii-cryo-et-object-identification/test/static/ExperimentRuns/*'))\nruns = [os.path.basename(x) for x in runs]\nruns[:5]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:14.018457Z","iopub.execute_input":"2025-02-05T06:48:14.018683Z","iopub.status.idle":"2025-02-05T06:48:14.029107Z","shell.execute_reply.started":"2025-02-05T06:48:14.018665Z","shell.execute_reply":"2025-02-05T06:48:14.028499Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"particle_names = [\n    'apo-ferritin',\n    'beta-amylase',\n    'beta-galactosidase',\n    'ribosome',\n    'thyroglobulin',\n    'virus-like-particle'\n]\n\nparticle_to_index = {\n    'apo-ferritin': 0,\n    'beta-amylase': 1,\n    'beta-galactosidase': 2,\n    'ribosome': 3,\n    'thyroglobulin': 4,\n    'virus-like-particle': 5\n}\n\nindex_to_particle = {index: name for name, index in particle_to_index.items()}\n\nparticle_radius = {\n    'apo-ferritin': 60,\n    'beta-amylase': 65,\n    'beta-galactosidase': 90,\n    'ribosome': 150,\n    'thyroglobulin': 130,\n    'virus-like-particle': 135,\n}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:14.030311Z","iopub.execute_input":"2025-02-05T06:48:14.030550Z","iopub.status.idle":"2025-02-05T06:48:14.046161Z","shell.execute_reply.started":"2025-02-05T06:48:14.030530Z","shell.execute_reply":"2025-02-05T06:48:14.045387Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ------------------- UnionFind 클래스 -------------------\nclass UnionFind:\n    def __init__(self, size):\n        self.parent = np.arange(size)\n        self.rank = np.zeros(size, dtype=int)\n\n    def find(self, u):\n        if self.parent[u] != u:\n            self.parent[u] = self.find(self.parent[u])\n        return self.parent[u]\n\n    def union(self, u, v):\n        u_root = self.find(u)\n        v_root = self.find(v)\n        if u_root == v_root:\n            return\n        if self.rank[u_root] < self.rank[v_root]:\n            self.parent[u_root] = v_root\n        else:\n            self.parent[v_root] = u_root\n            if self.rank[u_root] == self.rank[v_root]:\n                self.rank[u_root] += 1\n\n# ------------------- PredictionAggregator 클래스 -------------------\nclass PredictionAggregator:\n    def __init__(self, first_conf=0.2, conf_coef=0.75):\n        self.first_conf = first_conf\n        self.conf_coef = conf_coef\n        # 각 입자별 최소 aggregated confidence 기준 (튜닝 값)\n        self.particle_confs = np.array([0.5, 0.0, 0.2, 0.5, 0.2, 0.5])\n        \n    def convert_to_8bit(self, volume):\n        lower, upper = np.percentile(volume, (0.5, 99.5))\n        clipped = np.clip(volume, lower, upper)\n        scaled = ((clipped - lower) / (upper - lower + 1e-12) * 255).astype(np.uint8)\n        return scaled\n\n    def make_predictions(self, run_id, model, device_no):\n        volume_path = f'/kaggle/input/czii-cryo-et-object-identification/test/static/ExperimentRuns/{run_id}/VoxelSpacing10.000/denoised.zarr'\n        volume = zarr.open(volume_path, mode='r')[0]\n        volume_8bit = self.convert_to_8bit(volume)\n        num_slices = volume_8bit.shape[0]\n\n        detections = {\n            'particle_type': [],\n            'confidence': [],\n            'x': [],\n            'y': [],\n            'z': []\n        }\n\n        # 슬라이스별로 모델 예측 수행\n        for slice_idx in range(num_slices):\n            img = volume_8bit[slice_idx]\n            input_image = cv2.resize(np.stack([img]*3, axis=-1), (640, 640))\n\n            results = model.predict(\n                input_image,\n                save=False,\n                imgsz=640,\n                conf=self.first_conf,\n                device=device_no,\n                batch=8,\n                verbose=False,\n            )\n\n            for result in results:\n                boxes = result.boxes\n                if boxes is None:\n                    continue\n                cls = boxes.cls.cpu().numpy().astype(int)\n                conf = boxes.conf.cpu().numpy()\n                xyxy = boxes.xyxy.cpu().numpy()\n\n                # 중심 좌표 계산 (resize 보정: 63/64, scaling factor: 10)\n                xc = ((xyxy[:, 0] + xyxy[:, 2]) / 2.0) * 10 * (63/64)\n                yc = ((xyxy[:, 1] + xyxy[:, 3]) / 2.0) * 10 * (63/64)\n                zc = np.full(xc.shape, slice_idx * 10 + 5)\n\n                particle_types = [index_to_particle[c] for c in cls]\n\n                detections['particle_type'].extend(particle_types)\n                detections['confidence'].extend(conf)\n                detections['x'].extend(xc)\n                detections['y'].extend(yc)\n                detections['z'].extend(zc)\n\n        if not detections['particle_type']:\n            return pd.DataFrame()  \n\n        particle_types = np.array(detections['particle_type'])\n        confidences = np.array(detections['confidence'])\n        xs = np.array(detections['x'])\n        ys = np.array(detections['y'])\n        zs = np.array(detections['z'])\n\n        aggregated_data = []\n\n        # 각 particle type 별로 후처리 수행 (beta-amylase는 스킵)\n        for idx, particle in enumerate(particle_names):\n            if particle == 'beta-amylase':\n                continue\n\n            mask = (particle_types == particle)\n            if not np.any(mask):\n                continue\n\n            particle_confidences = confidences[mask]\n            particle_xs = xs[mask]\n            particle_ys = ys[mask]\n            particle_zs = zs[mask]\n            coords = np.vstack((particle_xs, particle_ys, particle_zs)).T  # (N, 3) [x,y,z]\n\n            # 클러스터링 진행 (고정 임계값 사용)\n            z_distance = 20  # z축 차이 허용 범위\n            xy_distance = 20 # xy 평면 차이 허용 범위\n            max_distance = math.sqrt(z_distance**2 + xy_distance**2)\n            tree = cKDTree(coords)\n            pairs = tree.query_pairs(r=max_distance, p=2)\n            \n            uf = UnionFind(len(coords))\n            coords_xy = coords[:, :2]\n            coords_z = coords[:, 2]\n            for u, v in pairs:\n                if abs(coords_z[u] - coords_z[v]) > z_distance:\n                    continue\n                if np.linalg.norm(coords_xy[u] - coords_xy[v]) > xy_distance:\n                    continue\n                uf.union(u, v)\n            \n            # 클러스터링 결과 계산\n            roots = np.array([uf.find(i) for i in range(len(coords))])\n            unique_roots, inverse_indices, counts = np.unique(roots, return_inverse=True, return_counts=True)\n            conf_sums = np.bincount(inverse_indices, weights=particle_confidences)\n            aggregated_confidences = conf_sums / (counts ** self.conf_coef)\n            \n            # 각 클러스터 내부의 좌표 분산(표준편차) 계산: 각 클러스터에 대해 평균 표준편차를 구함\n            n_clusters = len(unique_roots)\n            std_features = []\n            for i in range(n_clusters):\n                indices = np.where(inverse_indices == i)[0]\n                cluster_coords = coords[indices]\n                std_vals = np.std(cluster_coords, axis=0)  # 각 축별 표준편차\n                avg_std = np.mean(std_vals)\n                std_features.append(avg_std)\n            std_features = np.array(std_features)\n            \n            # 각 클러스터에 대한 feature vector 구성: [클러스터 내 점 개수, aggregated confidence, 평균 분산]\n            features = np.vstack((counts, aggregated_confidences, std_features)).T  # shape: (n_clusters, 3)\n            \n            # Isolation Forest를 적용하여 outlier 클러스터 탐지 (5% outlier 비율)\n            iso_forest = IsolationForest(contamination=0.05, random_state=42)\n            outlier_labels = iso_forest.fit_predict(features)  # 1: inlier, -1: outlier\n            inlier_mask = (outlier_labels == 1)\n            \n            # 기존 valid_clusters 조건: 최소 클러스터 점 수 및 aggregated confidence 기준\n            cluster_per_particle = [4, 1, 2, 9, 4, 8]\n            valid_clusters = (counts >= cluster_per_particle[idx]) & (aggregated_confidences > self.particle_confs[idx])\n            \n            # 최종 valid 클러스터: 기존 조건과 Isolation Forest의 inlier 결과를 모두 만족\n            final_valid_clusters = valid_clusters & inlier_mask\n            \n            if not np.any(final_valid_clusters):\n                continue\n            \n            # 클러스터 중심 좌표 계산 (각 클러스터의 평균 좌표)\n            centers_x = np.bincount(inverse_indices, weights=particle_xs) / counts\n            centers_y = np.bincount(inverse_indices, weights=particle_ys) / counts\n            centers_z = np.bincount(inverse_indices, weights=particle_zs) / counts\n            \n            centers_x = centers_x[final_valid_clusters]\n            centers_y = centers_y[final_valid_clusters]\n            centers_z = centers_z[final_valid_clusters]\n            \n            # ------------------- 2차 후처리: Heuristic Scoring -------------------\n            # 각 클러스터의 특징: 점 개수, aggregated confidence, 평균 분산\n            final_counts = counts[final_valid_clusters]\n            final_conf = aggregated_confidences[final_valid_clusters]\n            final_std = std_features[final_valid_clusters]\n            epsilon = 1e-6  # 0으로 나누는 경우 방지\n            # heuristic scoring: 높은 aggregated confidence와 많은 점 수, 낮은 분산일수록 높은 점수를 부여\n            scores = final_conf * (final_counts / (final_std + epsilon))\n            threshold_score = 1.0  # 임계값 (데이터에 따라 조정)\n            final_mask = scores >= threshold_score\n            \n            if not np.any(final_mask):\n                continue\n            \n            final_centers_x = centers_x[final_mask]\n            final_centers_y = centers_y[final_mask]\n            final_centers_z = centers_z[final_mask]\n            \n            aggregated_df = pd.DataFrame({\n                'experiment': [run_id] * len(final_centers_x),\n                'particle_type': [particle] * len(final_centers_x),\n                'x': final_centers_x,\n                'y': final_centers_y,\n                'z': final_centers_z\n            })\n            \n            aggregated_data.append(aggregated_df)\n\n        if aggregated_data:\n            return pd.concat(aggregated_data, ignore_index=True)\n        else:\n            return pd.DataFrame()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:14.047004Z","iopub.execute_input":"2025-02-05T06:48:14.047325Z","iopub.status.idle":"2025-02-05T06:48:14.070295Z","shell.execute_reply.started":"2025-02-05T06:48:14.047287Z","shell.execute_reply":"2025-02-05T06:48:14.069636Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"agent = PredictionAggregator(first_conf=0.5,  conf_coef=0.5)\nsubs = []\n\ntick = time.time()\nfor r in tqdm(runs, total=len(runs)):\n    df = agent.make_predictions(r, model, \"0\")\n    subs.append(df)\ntock = time.time()\nsubmission_ = pd.concat(subs).reset_index(drop=True)\nsubmission_.insert(0, 'id', range(len(submission_)))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:14.070921Z","iopub.execute_input":"2025-02-05T06:48:14.071099Z","iopub.status.idle":"2025-02-05T06:48:51.413823Z","shell.execute_reply.started":"2025-02-05T06:48:14.071083Z","shell.execute_reply":"2025-02-05T06:48:51.413048Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# agent = PredictionAggregator(first_conf=0.5, conf_coef=0.5)\n# subs = []\n\n# tick = time.time()\n# for r in tqdm(runs, total=len(runs)):\n#     df = agent.make_predictions(r, model, \"0\")\n#     subs.append(df)\n# tock = time.time()\n\n# submission_ = pd.concat(subs).reset_index(drop=True)\n# submission_.insert(0, 'id', range(len(submission_)))\n\n# # 'beta-galactosidase'와 'thyroglobulin' 행 제거\n# submission_ = submission_[~submission_['particle_type'].isin(['beta-galactosidase', 'thyroglobulin'])].reset_index(drop=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:51.414669Z","iopub.execute_input":"2025-02-05T06:48:51.415003Z","iopub.status.idle":"2025-02-05T06:48:51.418338Z","shell.execute_reply.started":"2025-02-05T06:48:51.414976Z","shell.execute_reply":"2025-02-05T06:48:51.417491Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# #change by @minfuka\n# submission0 = pd.concat(results[0])\n# submission1 = pd.concat(results[1])\n# submission_ = pd.concat([submission0, submission1]).reset_index(drop=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:51.419233Z","iopub.execute_input":"2025-02-05T06:48:51.419592Z","iopub.status.idle":"2025-02-05T06:48:51.436738Z","shell.execute_reply.started":"2025-02-05T06:48:51.419561Z","shell.execute_reply":"2025-02-05T06:48:51.435824Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# submission_.insert(0, 'id', range(len(submission_)))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:51.437533Z","iopub.execute_input":"2025-02-05T06:48:51.437783Z","iopub.status.idle":"2025-02-05T06:48:51.450212Z","shell.execute_reply.started":"2025-02-05T06:48:51.437751Z","shell.execute_reply":"2025-02-05T06:48:51.449360Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# **《《《 Unet3D(Monai) 》》》**","metadata":{}},{"cell_type":"code","source":"class Model(pl.LightningModule):\n    def __init__(self, spatial_dims=3, in_channels=1, out_channels=7,\n                 channels=(48, 64, 80, 80), strides=(2, 2, 1),\n                 num_res_units=1, lr=1e-3,\n                 scheduler_type='one_cycle'):\n            super().__init__()\n            self.save_hyperparameters()\n\n            # Model\n            self.model = UNet(\n              spatial_dims=self.hparams.spatial_dims,\n              in_channels=self.hparams.in_channels,\n              out_channels=self.hparams.out_channels,\n              channels=self.hparams.channels,\n              strides=self.hparams.strides,\n              num_res_units=self.hparams.num_res_units,\n              norm='batch',  # BatchNorm3d\n              dropout=0.2,\n          )\n\n            # Loss function\n            self.loss_fn = TverskyLoss(\n                include_background=True,\n                to_onehot_y=True,\n                softmax=True,\n                alpha=0.5,\n                beta=0.95\n            )\n\n            # Metric\n            self.metric_fn = DiceMetric(\n                include_background=False,\n                reduction=\"mean\",\n                get_not_nans=False\n            )\n\n            # Learning rate와 scheduler 설정\n            self.lr = lr\n            self.scheduler_type = scheduler_type\n\n            # 결과 저장용 리스트\n            self.training_step_outputs = []\n            self.validation_step_outputs = []\n\n            # Class weights 정의\n            self.class_weights = torch.tensor([1.0, 1.0, 1.0, 2.0, 2.0, 0.0])\n\n            # Storage for validation outputs\n            self.validation_outputs = []\n\n    def forward(self, x):\n        return self.model(x)\n\n    def validation_step(self, batch, batch_idx):\n        x, y = batch['image'], batch['label']\n        y_hat = self(x)\n        val_loss = self.loss_fn(y_hat, y)\n\n        metric_val_outputs = [AsDiscrete(argmax=True, to_onehot=self.hparams.out_channels)(i)\n                             for i in decollate_batch(y_hat)]\n        metric_val_labels = [AsDiscrete(to_onehot=self.hparams.out_channels)(i)\n                            for i in decollate_batch(y)]\n\n        self.metric_fn(y_pred=metric_val_outputs, y=metric_val_labels)\n        metrics = self.metric_fn.aggregate(reduction=\"mean_batch\")\n\n        # 클래스별 가중치를 device로 이동\n        class_weights = self.class_weights.to(metrics.device)\n\n        # 가중치가 적용된 전체 메트릭\n        weighted_metric = (metrics * class_weights).sum() / class_weights.sum()\n\n        # 로깅\n        self.log('val_loss', val_loss, on_step=False, on_epoch=True)\n        self.log('val_metric', weighted_metric, on_step=False, on_epoch=True)\n\n        output = {\n            'val_loss': val_loss.detach(),\n            'val_metric': weighted_metric.detach(),\n            'class_metrics': metrics.detach()\n        }\n\n        self.validation_outputs.append(output)\n\n        # 출력\n        print(f\"\\nEpoch {self.current_epoch}, Validation batch {batch_idx}\")\n        print(f\"Loss: {val_loss:.4f}, Metric: {weighted_metric:.4f}\")\n\n        return output\n\n    def training_step(self, batch, batch_idx):\n        x, y = batch['image'], batch['label']\n        y_hat = self(x)\n        loss = self.loss_fn(y_hat, y)\n\n        # 메트릭 로깅 추가\n        self.log(\"train_loss\", loss, on_step=False, on_epoch=True)  # 여기를 추가\n\n        print(f\"Epoch {self.current_epoch}, Training batch {batch_idx}, Loss: {loss.item():.4f}\")\n\n        self.training_step_outputs.append(loss)\n        return loss\n\n    def on_train_epoch_end(self):\n        epoch_mean = torch.stack(self.training_step_outputs).mean()\n        print(f\"\\n{'='*40}\")\n        print(f\"Epoch {self.current_epoch} Training completed\")\n        print(f\"Average training loss: {epoch_mean:.4f}\")\n        print(f\"{'='*40}\\n\")\n        self.training_step_outputs.clear()\n\n    def on_validation_epoch_start(self):\n        self.validation_outputs = []\n\n    def on_validation_epoch_end(self):\n        if not self.validation_outputs:\n            print(\"No validation outputs found!\")\n            return\n\n        try:\n            # 평균 계산\n            avg_loss = torch.stack([x['val_loss'] for x in self.validation_outputs]).mean()\n            avg_metric = torch.stack([x['val_metric'] for x in self.validation_outputs]).mean()\n            class_metrics = torch.stack([x['class_metrics'] for x in self.validation_outputs]).mean(dim=0)\n\n            print(f\"\\n{'='*70}\")\n            print(f\"Validation Epoch {self.current_epoch} Summary\")\n            print(f\"{'='*70}\")\n            print(f\"Average Loss: {avg_loss:.4f}\")\n            print(f\"Average Weighted Metric: {avg_metric:.4f}\")\n            print(\"\\nClass-wise Performance:\")\n            class_names = ['Ribosome', 'Virus-like', 'Apo-ferritin',\n                          'Thyroglobulin (Hard)', 'β-galactosidase (Hard)', 'Beta-amylase (Not evaluated)']\n\n            for i, (name, metric) in enumerate(zip(class_names, class_metrics)):\n                print(f\"  {name:<20} {metric:.4f}\")\n            print(f\"{'='*70}\\n\")\n\n        except Exception as e:\n            print(f\"Error in validation epoch end: {str(e)}\")\n        finally:\n            # 메트릭 리셋\n            self.metric_fn.reset()\n            self.validation_outputs = []\n\n    def configure_optimizers(self):\n      optimizer = torch.optim.AdamW(self.parameters(), lr=self.lr)\n\n      if self.scheduler_type == 'one_cycle':\n          # dataloader가 설정되기 전에는 cosine scheduler를 사용\n          if not hasattr(self.trainer, 'train_dataloader') or self.trainer.train_dataloader is None:\n              print(\"Warning: train_dataloader not set, switching to cosine scheduler\")\n              scheduler = CosineAnnealingWarmRestarts(\n                  optimizer,\n                  T_0=10,\n                  T_mult=2,\n                  eta_min=1e-6\n              )\n              scheduler_config = {\n                  \"scheduler\": scheduler,\n                  \"interval\": \"epoch\",\n                  \"frequency\": 1\n              }\n          else:\n              steps_per_epoch = len(self.trainer.train_dataloader())\n              total_steps = steps_per_epoch * self.trainer.max_epochs\n\n              scheduler = OneCycleLR(\n                  optimizer,\n                  max_lr=self.lr,\n                  total_steps=total_steps,\n                  pct_start=0.3,\n                  div_factor=25.0,\n                  final_div_factor=1e4\n              )\n              scheduler_config = {\n                  \"scheduler\": scheduler,\n                  \"interval\": \"step\",\n                  \"frequency\": 1\n              }\n\n      elif self.scheduler_type == 'plateau':\n          scheduler = ReduceLROnPlateau(\n              optimizer,\n              mode='max',\n              factor=0.5,\n              patience=100,\n              min_lr=1e-6,\n              verbose=True\n          )\n          scheduler_config = {\n              \"scheduler\": scheduler,\n              \"interval\": \"epoch\",\n              \"monitor\": 'val_metric',\n              \"frequency\": 1\n          }\n\n      else:  # cosine as default\n          scheduler = CosineAnnealingWarmRestarts(\n              optimizer,\n              T_0=10,\n              T_mult=2,\n              eta_min=1e-6\n          )\n          scheduler_config = {\n              \"scheduler\": scheduler,\n              \"interval\": \"epoch\",\n              \"frequency\": 1\n          }\n\n      return {\"optimizer\": optimizer, \"lr_scheduler\": scheduler_config}\n\nchannels = (64, 128, 256, 256)\nstrides_pattern = (2, 2, 1)\nnum_res_units = 1\nlearning_rate = 1e-3\nnum_epochs = 1000\ndef extract_3d_patches_minimal_overlap(arrays: List[np.ndarray], patch_size: int) -> Tuple[List[np.ndarray], List[Tuple[int, int, int]]]:\n    if not arrays or not isinstance(arrays, list):\n        raise ValueError(\"Input must be a non-empty list of arrays\")\n    \n    # Verify all arrays have the same shape\n    shape = arrays[0].shape\n    if not all(arr.shape == shape for arr in arrays):\n        raise ValueError(\"All input arrays must have the same shape\")\n    \n    if patch_size > min(shape):\n        raise ValueError(f\"patch_size ({patch_size}) must be smaller than smallest dimension {min(shape)}\")\n    \n    m, n, l = shape\n    patches = []\n    coordinates = []\n    \n    # Calculate starting positions for each dimension\n    x_starts = calculate_patch_starts(m, patch_size)\n    y_starts = calculate_patch_starts(n, patch_size)\n    z_starts = calculate_patch_starts(l, patch_size)\n    \n    # Extract patches from each array\n    for arr in arrays:\n        for x in x_starts:\n            for y in y_starts:\n                for z in z_starts:\n                    patch = arr[\n                        x:x + patch_size,\n                        y:y + patch_size,\n                        z:z + patch_size\n                    ]\n                    patches.append(patch)\n                    coordinates.append((x, y, z))\n    \n    return patches, coordinates\ndef reconstruct_array(patches: List[np.ndarray], \n                     coordinates: List[Tuple[int, int, int]], \n                     original_shape: Tuple[int, int, int]) -> np.ndarray:\n    reconstructed = np.zeros(original_shape, dtype=np.int64)  # To track overlapping regions\n    \n    patch_size = patches[0].shape[0]\n    \n    for patch, (x, y, z) in zip(patches, coordinates):\n        reconstructed[\n            x:x + patch_size,\n            y:y + patch_size,\n            z:z + patch_size\n        ] = patch\n        \n    \n    return reconstructed\ndef calculate_patch_starts(dimension_size: int, patch_size: int) -> List[int]:\n    if dimension_size <= patch_size:\n        return [0]\n        \n    # Calculate number of patches needed\n    n_patches = np.ceil(dimension_size / patch_size)\n    \n    if n_patches == 1:\n        return [0]\n    \n    # Calculate overlap\n    total_overlap = (n_patches * patch_size - dimension_size) / (n_patches - 1)\n    \n    # Generate starting positions\n    positions = []\n    for i in range(int(n_patches)):\n        pos = int(i * (patch_size - total_overlap))\n        if pos + patch_size > dimension_size:\n            pos = dimension_size - patch_size\n        if pos not in positions:  # Avoid duplicates\n            positions.append(pos)\n    \n    return positions\nimport pandas as pd\n\ndef dict_to_df(coord_dict, experiment_name):\n    # Create lists to store data\n    all_coords = []\n    all_labels = []\n    \n    # Process each label and its coordinates\n    for label, coords in coord_dict.items():\n        all_coords.append(coords)\n        all_labels.extend([label] * len(coords))\n    \n    # Concatenate all coordinates\n    all_coords = np.vstack(all_coords)\n    \n    df = pd.DataFrame({\n        'experiment': experiment_name,\n        'particle_type': all_labels,\n        'x': all_coords[:, 0],\n        'y': all_coords[:, 1],\n        'z': all_coords[:, 2]\n    })\n\n    \n    return df\nfrom typing import List, Tuple, Union\nimport numpy as np\nimport torch\nfrom monai.data import DataLoader, Dataset, CacheDataset, decollate_batch\nfrom monai.transforms import (\n    Compose, \n    EnsureChannelFirstd, \n    Orientationd,  \n    AsDiscrete,  \n    RandFlipd, \n    RandRotate90d, \n    NormalizeIntensityd,\n    RandCropByLabelClassesd,\n)\nTRAIN_DATA_DIR = \"/kaggle/input/create-numpy-dataset-exp-name\"\nimport json\ncopick_config_path = TRAIN_DATA_DIR + \"/copick.config\"\n\nwith open(copick_config_path) as f:\n    copick_config = json.load(f)\n\ncopick_config['static_root'] = '/kaggle/input/czii-cryo-et-object-identification/test/static'\n\ncopick_test_config_path = 'copick_test.config'\n\nwith open(copick_test_config_path, 'w') as outfile:\n    json.dump(copick_config, outfile)\nimport copick\n\nroot = copick.from_file(copick_test_config_path)\n\ncopick_user_name = \"copickUtils\"\ncopick_segmentation_name = \"paintedPicks\"\nvoxel_size = 10\ntomo_type = \"denoised\"\ninference_transforms = Compose([\n    EnsureChannelFirstd(keys=[\"image\"], channel_dim=\"no_channel\"),\n    NormalizeIntensityd(keys=\"image\"),\n    Orientationd(keys=[\"image\"], axcodes=\"RAS\")\n])\nimport cc3d\n\nid_to_name = {1: \"apo-ferritin\", \n              2: \"beta-amylase\",\n              3: \"beta-galactosidase\", \n              4: \"ribosome\", \n              5: \"thyroglobulin\", \n              6: \"virus-like-particle\"}\nBLOB_THRESHOLD = 250\nCERTAINTY_THRESHOLD = 0.05\n\nclasses = [1, 2, 3, 4, 5, 6]\nimport torch\nimport numpy as np\nimport pandas as pd\nimport cc3d\nfrom monai.data import CacheDataset\nfrom monai.transforms import Compose, EnsureType\nfrom torch import nn\nfrom tqdm import tqdm\nfrom monai.networks.nets import UNet\nfrom monai.losses import TverskyLoss\nfrom monai.metrics import DiceMetric\n\ndef load_models(model_paths):\n    models = []\n    for model_path in model_paths:\n        channels = (64, 128, 256, 256)\n        strides_pattern = (2, 2, 1)\n        num_res_units = 1\n        learning_rate = 1e-3\n        num_epochs = 1000\n        model = Model(channels=channels, strides=strides_pattern, num_res_units=num_res_units, lr=learning_rate)\n        \n        weights =torch.load(model_path)['state_dict']\n        model.load_state_dict(weights)\n        model.to('cuda')\n        model.eval()\n        models.append(model)\n    return models\n\n\n# model_paths = [\n#     '/kaggle/input/20250122-v34-7fold/fold0.ckpt',\n# ]\n\n\n# models = load_models(model_paths)\n# def ensemble_prediction_tta(models, input_tensor, threshold=0.5):\n#     probs_list = []\n#     data_copy0 = input_tensor.clone()\n#     data_copy0=torch.flip(data_copy0, dims=[2])\n#     data_copy1 = input_tensor.clone()\n#     data_copy1=torch.flip(data_copy1, dims=[3])\n#     data_copy2 = input_tensor.clone()\n#     data_copy2=torch.flip(data_copy2, dims=[4])\n#     data_copy3 = input_tensor.clone()\n#     data_copy3 = data_copy3.rot90(1, dims=[3, 4])\n#     with torch.no_grad():\n#         model_output0 = model(input_tensor)\n#         model_output1 = model(data_copy0)\n#         model_output1=torch.flip(model_output1, dims=[2])\n#         model_output2 = model(data_copy1)\n#         model_output2=torch.flip(model_output2, dims=[3])\n#         model_output3 = model(data_copy2)\n#         model_output3=torch.flip(model_output3, dims=[4])\n#         probs0 = torch.softmax(model_output0[0], dim=0)\n#         probs1 = torch.softmax(model_output1[0], dim=0)\n#         probs2 = torch.softmax(model_output2[0], dim=0)\n#         probs3 = torch.softmax(model_output3[0], dim=0)\n#         probs_list.append(probs0)\n#         probs_list.append(probs1)\n#         probs_list.append(probs2)\n#         probs_list.append(probs3)\n#     avg_probs = torch.mean(torch.stack(probs_list), dim=0)\n#     thresh_probs = avg_probs > threshold\n#     _, max_classes = thresh_probs.max(dim=0)\n#     return max_classes\n# sub=[]\n# for model in models:\n#     with torch.no_grad():\n#         location_df = []\n#         for run in root.runs:\n#             tomo = run.get_voxel_spacing(10)\n#             tomo = tomo.get_tomogram(tomo_type).numpy()\n#             tomo_patches, coordinates = extract_3d_patches_minimal_overlap([tomo], 128)\n#             tomo_patched_data = [{\"image\": img} for img in tomo_patches]\n#             tomo_ds = CacheDataset(data=tomo_patched_data, transform=inference_transforms, cache_rate=1.0)\n#             pred_masks = []\n#             for i in tqdm(range(len(tomo_ds))):\n#                 input_tensor = tomo_ds[i]['image'].unsqueeze(0).to(\"cuda\")\n#                 max_classes = ensemble_prediction_tta(models, input_tensor, threshold=CERTAINTY_THRESHOLD)\n#                 pred_masks.append(max_classes.cpu().numpy())\n#             reconstructed_mask = reconstruct_array(pred_masks, coordinates, tomo.shape)\n#             location = {}\n#             for c in classes:\n#                 cc = cc3d.connected_components(reconstructed_mask == c)\n#                 stats = cc3d.statistics(cc)\n#                 zyx = stats['centroids'][1:] * 10.012444  # 转换单位\n#                 zyx_large = zyx[stats['voxel_counts'][1:] > BLOB_THRESHOLD]\n#                 xyz = np.ascontiguousarray(zyx_large[:, ::-1])\n#                 location[id_to_name[c]] = xyz\n#             df = dict_to_df(location, run.name)\n#             location_df.append(df)\n#         location_df = pd.concat(location_df)\n#         location_df.insert(loc=0, column='id', value=np.arange(len(location_df)))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:51.451134Z","iopub.execute_input":"2025-02-05T06:48:51.451449Z","iopub.status.idle":"2025-02-05T06:48:52.177002Z","shell.execute_reply.started":"2025-02-05T06:48:51.451404Z","shell.execute_reply":"2025-02-05T06:48:52.176030Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ckpt_paths = [\n    \"/kaggle/input/20250122-v34-7fold/fold0.ckpt\",\n    \"/kaggle/input/20250122-v34-7fold/fold1.ckpt\",\n    # \"/kaggle/input/20250122-v34-7fold/fold2.ckpt\",\n]\nmodel = Model(channels=channels, strides=strides_pattern, num_res_units=num_res_units, lr=learning_rate)\n# 여러 모델을 불러와 리스트에 저장\nmodels_ensemble = []\nfor cp in ckpt_paths:\n    m = Model.load_from_checkpoint(cp)\n    m.eval()\n    m.to(\"cuda\")\n    models_ensemble.append(m)\n\n# Non-random transforms to be cached\ninference_transforms = Compose([\n    EnsureChannelFirstd(keys=[\"image\"], channel_dim=\"no_channel\"),\n    NormalizeIntensityd(keys=\"image\"),\n    Orientationd(keys=[\"image\"], axcodes=\"RAS\")\n])\n\nimport numpy as np\nimport torch\nimport cc3d\nimport time  # 시간 측정용\n\nfrom monai.inferers import sliding_window_inference\nfrom monai.data import CacheDataset\n\n# 간단한 3D flip 함수들\ndef flip_x_3d(tensor: torch.Tensor) -> torch.Tensor:\n    return torch.flip(tensor, dims=[2])  # x축 뒤집기\n\ndef flip_y_3d(tensor: torch.Tensor) -> torch.Tensor:\n    return torch.flip(tensor, dims=[3])  # y축 뒤집기\n\ndef flip_z_3d(tensor: torch.Tensor) -> torch.Tensor:\n    return torch.flip(tensor, dims=[4])  # z축 뒤집기\n\ndef identity_3d(tensor: torch.Tensor) -> torch.Tensor:\n    return tensor\n\n# (forward_transform, inverse_transform) 쌍 목록\ntta_ops = [\n    (identity_3d, identity_3d),\n    # (flip_x_3d, flip_x_3d),\n    (flip_y_3d, flip_y_3d),\n    (flip_z_3d, flip_z_3d),\n]\n\ndef ensemble_tta_predictor(sub_volume: torch.Tensor) -> torch.Tensor:\n    \"\"\"\n    sub_volume: (B=1, C=1, D, H, W) 형태의 3D sub-volume Tensor.\n    여러 모델(models_ensemble) 및 TTA augmentation (예: identity, flip_y, flip_z)을 적용하여,\n    각 예측 결과를 confidence weighted voting 방식으로 재조합한 후 최종 logits를 반환합니다.\n    \n    반환: (B=1, out_channels=7, D, H, W)\n    \"\"\"\n    all_probs = []\n    all_weights = []\n    \n    with torch.no_grad():\n        # 각 모델과 TTA augmentation에 대해 예측 수행\n        for model in models_ensemble:\n            for fwd_op, inv_op in tta_ops:\n                # 1) 전처리: TTA 변환\n                vol_t = fwd_op(sub_volume)  # (1,1,D,H,W)\n                # 2) 모델 추론: logits 산출\n                logits = model(vol_t)       # (1,7,D,H,W)\n                # 3) 후처리: TTA inverse 변환\n                logits_inv = inv_op(logits) # (1,7,D,H,W)\n                # 4) 확률 맵 계산\n                probs = torch.softmax(logits_inv, dim=1)  # (1,7,D,H,W)\n                # 각 voxel에서 최대 확률값을 confidence로 사용 (채널 차원에서)\n                weight, _ = torch.max(probs, dim=1, keepdim=True)  # (1,1,D,H,W)\n                \n                all_probs.append(probs)\n                all_weights.append(weight)\n    \n    # 리스트를 tensor로 스택 (shape: (N, B, C, D, H, W) 및 (N, B, 1, D, H, W))\n    all_probs = torch.stack(all_probs, dim=0)\n    all_weights = torch.stack(all_weights, dim=0)\n    \n    # 가중 평균 확률 계산: 각 예측의 확률에 그 confidence weight를 곱한 후 합산하고, 가중치 합으로 나눔.\n    weighted_sum = torch.sum(all_probs * all_weights, dim=0)  # (B, C, D, H, W)\n    sum_weights = torch.sum(all_weights, dim=0)               # (B, 1, D, H, W)\n    final_probs = weighted_sum / (sum_weights + 1e-8)           # (B, C, D, H, W)\n    \n    # 최종 확률 맵을 logits로 변환 (후속 처리에서 torch.softmax을 적용하기 전에, logit 형태로 제공)\n    final_logits = torch.log(final_probs + 1e-8)\n    return final_logits\n\n\n# -----------------------------------------------------------\n# 슬라이딩 윈도우 + 앙상블 + TTA 인퍼런스 (후처리 없음) 예시\n# -----------------------------------------------------------\nBLOB_THRESHOLD = 250\nclasses = [1, 2, 3, 4, 5, 6]\n\nwith torch.no_grad():\n    location_df2 = []\n\n    for run in root.runs:\n        run_start = time.time()  # 한 번의 run 시작 시간\n\n        print(run)\n\n        # 1) 볼륨(10Å voxel) 로드\n        load_start = time.time()\n        tomo = run.get_voxel_spacing(10)\n        tomo_arr = tomo.get_tomogram(tomo_type).numpy()  # shape: (X, Y, Z)\n        load_end = time.time()\n        print(f\"[Timer] Volume load time: {load_end - load_start:.3f} sec\")\n\n        # 2) Dataset 로드(전처리)\n        prep_start = time.time()\n        data_dict = [{\"image\": tomo_arr}]\n        tomo_ds = CacheDataset(data=data_dict, transform=inference_transforms, cache_rate=1.0)\n        volume_tensor = tomo_ds[0][\"image\"].unsqueeze(0).to(\"cuda\")  # (1,1,X,Y,Z)\n        prep_end = time.time()\n        print(f\"[Timer] Dataset prep time: {prep_end - prep_start:.3f} sec\")\n\n        # 3) Sliding Window Inference\n        #    -> predictor=ensemble_tta_predictor 로 교체\n        infer_start = time.time()\n        out_logits = sliding_window_inference(\n            inputs=volume_tensor,\n            roi_size=(128, 128, 128),\n            sw_batch_size=4,\n            predictor=ensemble_tta_predictor,  # <-- 여기서 앙상블+TTA 진행\n            overlap=0.25,\n            mode=\"gaussian\"\n        )\n        infer_end = time.time()\n        print(f\"[Timer] SW Inference(Ensemble+TTA) time: {infer_end - infer_start:.3f} sec\")\n\n        # 4) Softmax 후 argmax\n        post_start = time.time()\n        out_probs = torch.softmax(out_logits, dim=1)  # (1,7,X,Y,Z)\n        out_probs_np = out_probs[0].cpu().numpy()     # (7, X, Y, Z)\n        reconstructed_mask = np.argmax(out_probs_np, axis=0)  # (X, Y, Z)\n        post_end = time.time()\n        print(f\"[Timer] Postprocess(softmax+argmax) time: {post_end - post_start:.3f} sec\")\n\n        # 클래스별 aspect ratio 임계값 설정 (실제 모양 정보를 반영)\n        aspect_ratio_thresholds = {\n            \"apo-ferritin\": 0.8,\n            \"beta-galactosidase\": 0.6,\n            \"ribosome\": 0.65,\n            \"thyroglobulin\": 0.6,\n            \"virus-like-particle\": 0.8\n        }\n        \n        cc_start = time.time()\n        location = {}\n        \n        # classes는 [1, 2, 3, 4, 5, 6]와 같이 번호가 할당되어 있다고 가정하고,\n        # id_to_name을 통해 실제 클래스 이름으로 매핑합니다.\n        for c in classes:\n            # reconstructed_mask == c 인 영역에 대해 연결 요소 분석 실행\n            cc = cc3d.connected_components(reconstructed_mask == c)\n            stats = cc3d.statistics(cc)\n            \n            # 배경 label(0)은 제외하고 실제 객체는 인덱스 1부터 사용\n            centroids = np.array(stats[\"centroids\"][1:])  # 보통 [z, y, x] 순서로 반환됨\n            bbox_list = stats[\"bounding_boxes\"][1:]         # 각 객체의 bounding box, 형식: (slice_z, slice_y, slice_x)\n            voxel_counts = np.array(stats[\"voxel_counts\"][1:])\n            \n            valid_indices = []\n            # 해당 클래스의 이름과 임계값을 설정 (없으면 기본값 0.7 사용)\n            class_name = id_to_name[c]\n            ar_threshold = aspect_ratio_thresholds.get(class_name, 0.7)\n            \n            for i, bbox in enumerate(bbox_list):\n                # 각 bounding box는 (slice_z, slice_y, slice_x) 형태이며,\n                # 각 slice 객체에서 start와 stop 값을 사용하여 길이를 계산합니다.\n                slice_z, slice_y, slice_x = bbox  # 각각의 축에 대한 slice 객체\n                depth  = slice_z.stop - slice_z.start\n                height = slice_y.stop - slice_y.start\n                width  = slice_x.stop - slice_x.start\n                \n                # Aspect ratio: 최소 길이를 최대 길이로 나눔 (1에 가까울수록 구형)\n                aspect_ratio = min(width, height, depth) / max(width, height, depth)\n                \n                # 조건: voxel 수가 BLOB_THRESHOLD 이상이며, aspect ratio가 클래스에 맞는 임계값 이상일 경우\n                if voxel_counts[i] > BLOB_THRESHOLD and aspect_ratio >= ar_threshold:\n                    valid_indices.append(i)\n            \n            if len(valid_indices) > 0:\n                # valid한 객체들의 중심 좌표를 선택합니다.\n                valid_centroids = centroids[valid_indices]\n                # voxel 크기 보정을 위해 (예: 10.012444)를 곱합니다.\n                valid_centroids = valid_centroids * 10.012444  \n                # cc3d 결과의 좌표 순서가 [z, y, x]이면, 이를 [x, y, z]로 변경합니다.\n                valid_xyz = np.ascontiguousarray(valid_centroids[:, ::-1])\n                location[class_name] = valid_xyz\n        \n        cc_end = time.time()\n        print(f\"[Timer] Connected components + centroids + shape filtering time: {cc_end - cc_start:.3f} sec\")\n        \n        # 이후 dict_to_df 함수를 이용해 DataFrame으로 변환하여 location_df2에 저장합니다.\n        df = dict_to_df(location, run.name)\n        location_df2.append(df)\n\n        run_end = time.time()\n        print(f\"[Timer] Single run total time: {run_end - run_start:.3f} sec\\n\")\n\n    # 모든 run 결과 결합\n    location_df2 = pd.concat(location_df2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:48:52.177848Z","iopub.execute_input":"2025-02-05T06:48:52.178542Z","iopub.status.idle":"2025-02-05T06:51:57.570322Z","shell.execute_reply.started":"2025-02-05T06:48:52.178517Z","shell.execute_reply":"2025-02-05T06:51:57.569537Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# **《《《 Finaly Blend 》》》**","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom sklearn.cluster import DBSCAN\n\n# submission_ 또는 다른 DataFrame과 location_df2를 결합하는 예시\n# df = pd.concat([submission_, location_df, location_df2], ignore_index=True)\ndf = pd.concat([submission_, location_df2], ignore_index=True)\n# df = pd.concat([location_df, location_df2], ignore_index=True)\n# df = location_df.copy()\n\nparticle_names = [\n    'apo-ferritin', \n    'beta-amylase', \n    'beta-galactosidase', \n    'ribosome', \n    'thyroglobulin', \n    'virus-like-particle'\n]\nparticle_radius = {\n    'apo-ferritin': 60,\n    'beta-amylase': 65,\n    'beta-galactosidase': 90,\n    'ribosome': 150,\n    'thyroglobulin': 130,\n    'virus-like-particle': 135,\n}\n\nfinal = []\n\n# 각 particle_type 별로 처리\nfor pidx, p in enumerate(particle_names):\n    pdf = df[df['particle_type'] == p].reset_index(drop=True)\n    p_rad = particle_radius[p]\n    \n    # experiment별로 그룹화\n    grouped = pdf.groupby(['experiment'])\n    \n    for exp, group in grouped:\n        group = group.reset_index(drop=True)\n        \n        coords = group[['x', 'y', 'z']].values\n        # DBSCAN 클러스터링 진행 (eps와 min_samples는 필요에 따라 조정)\n        db = DBSCAN(eps=p_rad, min_samples=2, metric='euclidean').fit(coords)\n        labels = db.labels_\n        group['cluster'] = labels\n        \n        # 각 클러스터에 대해 후처리 진행\n        for cluster_id in np.unique(labels):\n            if cluster_id == -1:\n                continue  # noise는 건너뜁니다.\n            \n            cluster_points = group[group['cluster'] == cluster_id]\n            \n            # 클러스터 내 좌표 평균 및 표준편차 계산\n            avg_coords = cluster_points[['x', 'y', 'z']].mean().values\n            std_coords = cluster_points[['x', 'y', 'z']].std().values\n            \n            # 각 좌표의 표준편차가 particle_radius의 30% 이상이면 해당 클러스터 무시\n            if np.any(std_coords > (p_rad * 0.4)):\n                # 해당 클러스터는 분산이 너무 크므로 noise로 재설정\n                group.loc[group['cluster'] == cluster_id, 'cluster'] = -1\n                continue\n            \n            # 클러스터 내 모든 좌표를 평균 좌표로 대체 후 중복 제거\n            group.loc[group['cluster'] == cluster_id, ['x', 'y', 'z']] = avg_coords\n            group = group.drop_duplicates(subset=['x', 'y', 'z'])\n            \n        final.append(group)\n\n# 모든 그룹을 합치고 cluster 열 제거\ndf_save = pd.concat(final, ignore_index=True)\ndf_save = df_save.drop(columns=['cluster'])\ndf_save = df_save.sort_values(by=['experiment', 'particle_type']).reset_index(drop=True)\ndf_save['id'] = np.arange(0, len(df_save))\n\ndf_save.to_csv('submission.csv', index=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:51:57.571377Z","iopub.execute_input":"2025-02-05T06:51:57.571683Z","iopub.status.idle":"2025-02-05T06:51:58.587800Z","shell.execute_reply.started":"2025-02-05T06:51:57.571659Z","shell.execute_reply":"2025-02-05T06:51:58.586858Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"total_end = time.time()\nprint(f\"[Total] 전체 파이프라인 수행 시간: {total_end - total_start:.3f} 초\")\nprint(f'estimated predict time is {(total_end - total_start)/3*500:.4f} seconds')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:51:58.588720Z","iopub.execute_input":"2025-02-05T06:51:58.588966Z","iopub.status.idle":"2025-02-05T06:51:58.594241Z","shell.execute_reply.started":"2025-02-05T06:51:58.588945Z","shell.execute_reply":"2025-02-05T06:51:58.593449Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!cp -r /kaggle/input/hengck-czii-cryo-et-01/* .","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:51:58.595065Z","iopub.execute_input":"2025-02-05T06:51:58.595285Z","iopub.status.idle":"2025-02-05T06:51:58.923797Z","shell.execute_reply.started":"2025-02-05T06:51:58.595258Z","shell.execute_reply":"2025-02-05T06:51:58.922584Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from czii_helper import *\nfrom dataset import *\nfrom scipy.optimize import linear_sum_assignment\nimport matplotlib.pyplot as plt","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:51:58.925190Z","iopub.execute_input":"2025-02-05T06:51:58.925601Z","iopub.status.idle":"2025-02-05T06:51:58.930101Z","shell.execute_reply.started":"2025-02-05T06:51:58.925567Z","shell.execute_reply":"2025-02-05T06:51:58.929361Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    MODE = 'submit'\nelse:\n    MODE = 'local'\n\n\n\n\n\n\n\nvalid_dir ='/kaggle/input/czii-cryo-et-object-identification/train'\nvalid_id = ['TS_5_4', 'TS_69_2', 'TS_6_4']\n\ndef do_one_eval(truth, predict, threshold):\n    P=len(predict)\n    T=len(truth)\n\n    if P==0:\n        hit=[[],[]]\n        miss=np.arange(T).tolist()\n        fp=[]\n        metric = [P,T,len(hit[0]),len(miss),len(fp)]\n        return hit, fp, miss, metric\n\n    if T==0:\n        hit=[[],[]]\n        fp=np.arange(P).tolist()\n        miss=[]\n        metric = [P,T,len(hit[0]),len(miss),len(fp)]\n        return hit, fp, miss, metric\n\n    #---\n    distance = predict.reshape(P,1,3)-truth.reshape(1,T,3)\n    distance = distance**2\n    distance = distance.sum(axis=2)\n    distance = np.sqrt(distance)\n    p_index, t_index = linear_sum_assignment(distance)\n\n    valid = distance[p_index, t_index] <= threshold\n    p_index = p_index[valid]\n    t_index = t_index[valid]\n    hit = [p_index.tolist(), t_index.tolist()]\n    miss = np.arange(T)\n    miss = miss[~np.isin(miss,t_index)].tolist()\n    fp = np.arange(P)\n    fp = fp[~np.isin(fp,p_index)].tolist()\n\n    metric = [P,T,len(hit[0]),len(miss),len(fp)] #for lb metric F-beta copmutation\n    return hit, fp, miss, metric\n\n\ndef compute_lb(submit_df, overlay_dir):\n    valid_id = list(submit_df['experiment'].unique())\n    print(valid_id)\n\n    eval_df = []\n    for id in valid_id:\n        truth = read_one_truth(id, overlay_dir) #=f'{valid_dir}/overlay/ExperimentRuns')\n        id_df = submit_df[submit_df['experiment'] == id]\n        for p in PARTICLE:\n            p = dotdict(p)\n            print('\\r', id, p.name, end='', flush=True)\n            xyz_truth = truth[p.name]\n            xyz_predict = id_df[id_df['particle_type'] == p.name][['x', 'y', 'z']].values\n            hit, fp, miss, metric = do_one_eval(xyz_truth, xyz_predict, p.radius* 0.5)\n            eval_df.append(dotdict(\n                id=id, particle_type=p.name,\n                P=metric[0], T=metric[1], hit=metric[2], miss=metric[3], fp=metric[4],\n            ))\n    print('')\n    eval_df = pd.DataFrame(eval_df)\n    gb = eval_df.groupby('particle_type').agg('sum').drop(columns=['id'])\n    gb.loc[:, 'precision'] = gb['hit'] / gb['P']\n    gb.loc[:, 'precision'] = gb['precision'].fillna(0)\n    gb.loc[:, 'recall'] = gb['hit'] / gb['T']\n    gb.loc[:, 'recall'] = gb['recall'].fillna(0)\n    gb.loc[:, 'f-beta4'] = 17 * gb['precision'] * gb['recall'] / (16 * gb['precision'] + gb['recall'])\n    gb.loc[:, 'f-beta4'] = gb['f-beta4'].fillna(0)\n\n    gb = gb.sort_values('particle_type').reset_index(drop=False)\n    # https://www.kaggle.com/competitions/czii-cryo-et-object-identification/discussion/544895\n    gb.loc[:, 'weight'] = [1, 0, 2, 1, 2, 1]\n    lb_score = (gb['f-beta4'] * gb['weight']).sum() / gb['weight'].sum()\n    return gb, lb_score\n\n\n#debug\nif 1:\n    if MODE=='local':\n    #if 1:\n        submit_df=pd.read_csv(\n           'submission.csv'\n            # '/kaggle/input/hengck-czii-cryo-et-weights-01/submission.csv'\n        )\n        gb, lb_score = compute_lb(submit_df, f'{valid_dir}/overlay/ExperimentRuns')\n        print(gb)\n        print('lb_score:',lb_score)\n        print('')\n\n\n        #show one ----------------------------------\n        fig = plt.figure(figsize=(18, 8))\n\n        # debug 시각화: valid_id에 있는 모든 experiment에 대해 시각화하기\n        for exp_id in valid_id:\n            # 해당 experiment의 정답 데이터를 읽어옵니다.\n            truth = read_one_truth(exp_id, overlay_dir=f'{valid_dir}/overlay/ExperimentRuns')\n            # submission 데이터에서 해당 experiment에 해당하는 부분만 선택합니다.\n            id_df = submit_df[submit_df['experiment'] == exp_id]\n            \n            # experiment 별 figure 생성 (여기서는 2행 3열의 subplot 구성)\n            fig = plt.figure(figsize=(18, 8))\n            for p in PARTICLE:\n                p = dotdict(p)\n                # 해당 particle type에 대한 정답과 예측 데이터를 추출합니다.\n                xyz_truth = truth[p.name]\n                xyz_predict = id_df[id_df['particle_type'] == p.name][['x', 'y', 'z']].values\n                \n                # 평가 함수 호출\n                hit, fp, miss, _ = do_one_eval(xyz_truth, xyz_predict, p.radius)\n                print(exp_id, p.name)\n                print('\\t num truth   :', len(xyz_truth))\n                print('\\t num predict :', len(xyz_predict))\n                print('\\t num hit     :', len(hit[0]))\n                print('\\t num fp      :', len(fp))\n                print('\\t num miss    :', len(miss))\n            \n                # subplot에 3D scatter plot 그리기\n                ax = fig.add_subplot(2, 3, p.label, projection='3d')\n                if hit[0]:\n                    # 예측 중 맞춘 점들 (빨간색 점)와 대응하는 정답 (빨간 테두리 원)\n                    pt_pred = xyz_predict[hit[0]]\n                    ax.scatter(pt_pred[:, 0], pt_pred[:, 1], pt_pred[:, 2], alpha=0.5, color='r')\n                    pt_truth = xyz_truth[hit[1]]\n                    ax.scatter(pt_truth[:, 0], pt_truth[:, 1], pt_truth[:, 2], s=80, facecolors='none', edgecolors='r')\n                if fp:\n                    # false positive (검은색 점)\n                    pt_fp = xyz_predict[fp]\n                    ax.scatter(pt_fp[:, 0], pt_fp[:, 1], pt_fp[:, 2], alpha=1, color='k')\n                if miss:\n                    # miss (검은 테두리 원)\n                    pt_miss = xyz_truth[miss]\n                    ax.scatter(pt_miss[:, 0], pt_miss[:, 1], pt_miss[:, 2], s=160, alpha=1, facecolors='none', edgecolors='k')\n            \n                ax.set_title(f'{p.name} ({p.difficulty})')\n            \n            plt.suptitle(f'Experiment: {exp_id}', fontsize=16)\n            plt.tight_layout()\n            plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-05T06:51:58.931197Z","iopub.execute_input":"2025-02-05T06:51:58.931541Z","iopub.status.idle":"2025-02-05T06:52:02.621724Z","shell.execute_reply.started":"2025-02-05T06:51:58.931510Z","shell.execute_reply":"2025-02-05T06:52:02.620749Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}