{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":45867,"databundleVersionId":6924515,"sourceType":"competition"},{"sourceId":1760030,"sourceType":"datasetVersion","datasetId":1046169},{"sourceId":6774553,"sourceType":"datasetVersion","datasetId":3898019},{"sourceId":7486642,"sourceType":"datasetVersion","datasetId":4358578},{"sourceId":7486941,"sourceType":"datasetVersion","datasetId":4358790}],"dockerImageVersionId":30635,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!ls /kaggle/input/pyvips-python-and-deb-package-gpu\n# intall the deb packages\n!yes | dpkg -i --force-depends /kaggle/input/pyvips-python-and-deb-package-gpu/linux_packages/archives/*.deb\n# install the python wrapper\n!pip install pyvips -f /kaggle/input/pyvips-python-and-deb-package-gpu/python_packages/ --no-index\n!pip list | grep pyvips\n!pip install /kaggle/input/einops-030/einops-0.3.0-py2.py3-none-any.whl","metadata":{"execution":{"iopub.status.busy":"2024-01-30T22:48:43.015351Z","iopub.execute_input":"2024-01-30T22:48:43.015847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport random\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\n\nfrom sklearn.metrics import balanced_accuracy_score\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\nimport torch.optim as optim\nimport torch.optim.lr_scheduler as lr_scheduler\n\nlabel_names = ['CC', 'EC', 'HGSC', 'LGSC', 'MC']\nnum_classes = len(label_names)\nlabel_dict = {label_names[i]: i for i in range(len(label_names))}\n\n_CLS_MAPPING = {\n    'CC': 0,\n    'EC': 1, \n    'HGSC': 2, \n    'LGSC': 3, \n    'MC': 4,\n    'Other': 5,\n}\n\n_CLS_INV_MAPPING = {\n    0: 'CC', \n    1: 'EC', \n    2: 'HGSC', \n    3: 'LGSC', \n    4: 'MC',\n    5: 'Other',\n}","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class WSIFeatDataset(Dataset):\n    def __init__(self, data_csv: pd.DataFrame, feature_dir: str, ratio, phase: int):\n        super().__init__()\n        self.data_csv = data_csv\n        self.feature_dir = feature_dir\n        self.ratio = ratio\n        assert phase in [0, 1]\n        self.phase = phase\n\n    def __len__(self):\n        return len(self.data_csv)\n\n    def __getitem__(self, idx):\n        sample = self.data_csv.iloc[idx]\n        file_name = str(sample['image_id'])\n\n        feature_np = np.load(os.path.join(self.feature_dir, file_name + '.npy'))\n        features = torch.from_numpy(feature_np)\n\n        random.shuffle(features)\n        if 0 < self.ratio <= 1:\n            features = features[:int(len(features) * self.ratio)]\n        elif self.ratio > 1:\n            features = features[:min(len(features), self.ratio)]\n\n        if self.phase == 0:\n            label = torch.tensor(label_dict[sample['label']])\n            return file_name, features, label\n        else:\n            return file_name, features","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ABMIL(nn.Module):\n    def __init__(self, in_dim, feat_dim, attn_dim, num_classes):\n        super().__init__()\n        self.downlinear = nn.Sequential(nn.Linear(in_dim, feat_dim), nn.ReLU())\n        self.attention_V = nn.Sequential(nn.Linear(feat_dim, attn_dim), nn.Tanh())\n        self.attention_U = nn.Sequential(nn.Linear(feat_dim, attn_dim), nn.Sigmoid())\n        self.attention_weights = nn.Linear(attn_dim, 1)\n        self.classifier = nn.Linear(feat_dim, num_classes)\n\n        self.apply(self._init_weights)\n\n    def _init_weights(self, m):\n        if isinstance(m, nn.Linear):\n            nn.init.normal_(m.weight, std=.02)\n            if isinstance(m, nn.Linear) and m.bias is not None:\n                nn.init.constant_(m.bias, 0)\n\n    def forward(self, x):\n        x = self.downlinear(x)\n\n        A_V = self.attention_V(x)\n        A_U = self.attention_U(x)\n        A = self.attention_weights(A_V * A_U)\n        A = torch.transpose(A, 1, 0)\n        A = torch.softmax(A, dim=1)\n        x = torch.mm(A, x)\n\n        scores = self.classifier(x)\n\n        return scores","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def set_seed(seed=42):\n    ##### why 42? The Answer to the Ultimate Question of Life, the Universe, and Everything is 42.\n    random.seed(seed) # python\n    np.random.seed(seed) # numpy\n    torch.manual_seed(seed) # pytorch\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = False\n    # Set a fixed value for the hash seed\n    os.environ['PYTHONHASHSEED'] = str(seed) # new add\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Trained = True\nset_seed()\nif(Trained):\n    data_dir = '/kaggle/input/UBC-OCEAN'\n    train_csv = pd.read_csv(os.path.join(data_dir, 'train.csv'))\n    train_data = train_csv.iloc[np.r_[0:100, 200:536]].reset_index(drop=True)\n    val_data = train_csv.iloc[100:200].reset_index(drop=True)\n    label_names = ['CC', 'EC', 'HGSC', 'LGSC', 'MC']\n    num_classes = len(label_names)\n    label_dict = {label_names[i]: i for i in range(len(label_names))}\n    epochs = 20\n    in_dim = 768\n    model_type = 16\n    ratio = 1.0\n    mil_type = 'abmil'\n    accumulate = True\n    test = False\n    feature_dir = '/kaggle/input/ubc-ocean-feature-set/Feature_set'\n    mil_model_name = f'final_wsi_vitp{model_type}_{mil_type}_{ratio}_{epochs}ep.pth'\n    device = torch.device('cuda:0')\n\n    train_dataset = WSIFeatDataset(train_data, feature_dir, ratio, 0) if test else WSIFeatDataset(train_csv, feature_dir, ratio, 0)\n    train_loader = DataLoader(train_dataset, batch_size=1, shuffle=True, num_workers=4, pin_memory=True)\n    val_dataset = WSIFeatDataset(val_data, feature_dir, ratio, 0)\n    val_loader = DataLoader(val_dataset, batch_size=1, shuffle=True, num_workers=4, pin_memory=True)\n\n    if mil_type == 'abmil':\n        model = ABMIL(in_dim, 512, 128, num_classes)\n    elif mil_type == 'dsmil':\n        model = DSMIL(IClassifier(in_dim, num_classes), BClassifier(in_dim, num_classes))\n    elif mil_type == 'transmil':\n        model = TransMIL(in_dim, num_classes)\n\n    optimizer = optim.Adam(model.parameters(), 5e-4, weight_decay=5e-4)\n    scheduler = lr_scheduler.CosineAnnealingLR(optimizer, epochs, 5e-5)\n\n    model = model.to(device)\n\n    max_acc = 0.\n    balanced_acc = 0.\n    for epoch in range(1, epochs + 1):\n        loss_sum = 0.\n        n = 0\n\n        loop = tqdm(train_loader, total=len(train_loader))\n        model.train()\n        for file_name, features, label in loop:\n            label = label.to(device)\n            features = features.squeeze(0).to(device)\n\n            if mil_type == 'abmil':\n                scores = model(features)\n                loss = F.cross_entropy(scores, label)\n            elif mil_type == 'dsmil':\n                classes, bag_prediction, _, _ = model(features)\n                max_prediction, index = torch.max(classes, 0, True)\n                loss_bag = F.cross_entropy(bag_prediction, label)\n                loss_max = F.cross_entropy(max_prediction.view(1, -1), label)\n                loss = 0.5 * loss_bag + 0.5 * loss_max\n            elif mil_type == 'transmil':\n                scores = model(features.unsqueeze(0))\n                loss = F.cross_entropy(scores, label)\n\n            if accumulate:\n                loss = loss / 4\n                loss.backward(retain_graph=True)\n                if (n + 1) % 4 == 0 or (n + 1) == len(train_loader):\n                    optimizer.step()\n                    optimizer.zero_grad()\n            else:\n                optimizer.zero_grad()\n                loss.backward()\n                optimizer.step()\n\n            n += 1\n            loss_sum += loss.item()\n\n            loop.set_description(f'Train [{epoch}/{epochs}]')\n            loop.set_postfix(loss=loss.item(), loss_mean=loss_sum / n)\n\n        if test:\n            submission = pd.DataFrame(columns=['image_id', 'label'])\n            with torch.no_grad():\n                acc = 0\n                loss_val = 0\n                y_true = []\n                y_pred = []\n\n                loop = tqdm(val_loader, total=len(val_loader))\n                model.eval()\n\n                total_length = len(val_loader)\n                for i, (file_name, features, label) in enumerate(loop):\n                    label = label.to(device)\n                    features = features.squeeze(0).to(device)\n\n                    if mil_type == 'abmil':\n                        scores = model(features)\n                        scores = torch.softmax(scores, 1)\n                    elif mil_type == 'dsmil':\n                        classes, bag_prediction, _, _ = model(features)\n                        max_prediction, index = torch.max(classes, 0, True)\n                        scores = 0.5 * torch.softmax(max_prediction, 1) + 0.5 * torch.softmax(bag_prediction, 1)\n                    elif mil_type == 'transmil':\n                        scores = model(features.unsqueeze(0))\n                        scores = torch.softmax(scores, 1)\n\n                    pred = torch.argmax(scores)\n\n                    y_pred.append(pred.item())\n                    y_true.append(label.item())\n\n                    if pred == label.squeeze(0):\n                        acc += 1\n\n                    if epoch == epochs:\n                        submission.loc[len(submission)] = [file_name, _CLS_INV_MAPPING[pred.item()]]\n                    loop.set_description(f'Val [{epoch}/{epochs}]')\n                    loop.set_postfix(acc=acc / len(val_loader), max_acc=max_acc, balanced_acc=balanced_acc)\n\n                if acc / len(val_loader) > max_acc:\n                    max_acc = acc / len(val_loader)\n                y_true = np.array(y_true)\n                y_pred = np.array(y_pred)\n                if balanced_accuracy_score(y_true, y_pred) > balanced_acc:\n                    balanced_acc = balanced_accuracy_score(y_true, y_pred)\n\n            submission.to_csv('validation.csv', index=False)\n\n        scheduler.step()\n\n        if not test:\n            torch.save(model.state_dict(), mil_model_name)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport cv2\nfrom PIL import Image\nMAX_IMAGE_PIXELS = 10_000_000_000 \nImage.MAX_IMAGE_PIXELS = MAX_IMAGE_PIXELS\nfrom pathlib import Path\nfrom typing import List, Tuple, Callable, Optional, Union, Dict, Literal, Any\nfrom matplotlib.colors import rgb_to_hsv\nfrom skimage.filters import threshold_otsu\nfrom transformers import AutoImageProcessor, ViTModel\n\ndef get_metadata(slide_id: int, _DATA_DIR, IS_TRAINING = False):\n    \"\"\"Returns the metadata of the slide of id slide_id contained in the train.csv file.\"\"\"\n    if IS_TRAINING:\n        csv_path = _DATA_DIR / \"train.csv\"\n    else:\n        csv_path = _DATA_DIR / \"test.csv\"\n    df_ = pd.read_csv(csv_path)\n    df_ = df_.set_index(\"image_id\")\n    return df_.loc[slide_id]\n\ndef is_tma(slide_id: int, _DATA_DIR) -> bool:\n    \"\"\"Determines whether a slide is a TMA.\"\"\"\n    slide_metadata = get_metadata(slide_id, _DATA_DIR)\n    return slide_metadata.is_tma\n\ndef get_path_fullres_image(slide_id: int, DATA_DIR, IS_TRAINING = False) -> str:\n    \"\"\"Path to a full-resolution image.\"\"\"\n    if IS_TRAINING:\n        img_dir = Path(DATA_DIR).joinpath('train_images')\n    else:\n        img_dir = Path(DATA_DIR).joinpath('test_images')\n    img_path = img_dir.joinpath(f'{slide_id}.png')\n    return str(img_path)\n\ndef get_image_dimensions(slide_id: int) -> Tuple[int, int]:\n    \"\"\"Reads width and height from a full-resolution image.\"\"\"\n    try:\n        Image.MAX_IMAGE_PIXELS = 10_000_000_000\n        fullres_fpath = get_path_fullres_image(slide_id)\n        img = Image.open(fullres_fpath)\n        fullres_w, fullres_h = img.size\n        return fullres_w, fullres_h\n    except Exception:  # noqa\n        slide_metadata = get_metadata(slide_id)\n        fullres_w, fullres_h = slide_metadata.image_width, slide_metadata.image_height\n        return fullres_w, fullres_h\n\ndef get_path_thumbnail_image(slide_id: int, DATA_DIR, IS_TRAINING = False) -> str:\n    \"\"\"Path to a thumbnail image (if it exists).\"\"\"\n    if IS_TRAINING and is_tma(slide_id, DATA_DIR):\n        img_dir = Path(DATA_DIR).joinpath('train_images')\n        img_path = img_dir.joinpath(f'{slide_id}.png')\n        return str(img_path)\n    elif IS_TRAINING and (not is_tma(slide_id, DATA_DIR)):\n        img_dir = Path(DATA_DIR).joinpath('train_thumbnails')\n        img_path = img_dir.joinpath(f'{slide_id}_thumbnail.png')\n        return str(img_path)\n    elif not IS_TRAINING:\n        # Note: \n        # - if a thumbnail exist, we use it\n        # - otherwise, we use the full resolution image\n        img_dir = Path(DATA_DIR).joinpath('test_thumbnails')\n        expected_img_path = img_dir.joinpath(f'{slide_id}_thumbnail.png')\n        if expected_img_path.is_file():\n            return str(expected_img_path)\n        else:\n            return get_path_fullres_image(slide_id)\n\n\ndef _compute_matter_mask_from_array(arr: np.ndarray) -> np.ndarray:\n    # Convert to HSV\n    _hsv = rgb_to_hsv(arr)\n\n    # Apply matter detection (using Otsu thresholding)\n    threshold_h = threshold_otsu(_hsv[:, :, 0])\n    threshold_s = threshold_otsu(_hsv[:, :, 1])\n    _mask = np.logical_and(_hsv[:, :, 0] > threshold_h, _hsv[:, :, 1] > threshold_s)\n    _kernel = np.ones((1, 1))\n    mask = cv2.dilate(_mask.astype(np.int32), _kernel, iterations=1)\n\n    return mask\n\ndef compute_matter_mask(slide_id: int, DATA_DIR):\n    \"\"\"Applies matter detection to an image (thumbnail).\"\"\"\n    # Load thumbnail\n    path_thumbnail = get_path_thumbnail_image(slide_id, DATA_DIR)\n    thumbnail_img = Image.open(path_thumbnail)\n    if thumbnail_img.im is None:\n        thumbnail_img.load()\n    thumbnail_arr = np.asarray(thumbnail_img)\n        \n    # Compute matter mask\n    mask = _compute_matter_mask_from_array(thumbnail_arr)\n    return mask\n\ndef validate_weights_path(weights_path: Union[str, Path]) -> None:\n    \"\"\"Ensures that model weights exist.\"\"\"\n    if not Path(weights_path).is_file():\n        raise FileNotFoundError(f'Could not load model weights ({weights_path})!')\n\ndef calculate_entropy(prob_dist):\n    # Ensuring the probability distribution sums to 1\n    prob_dist = prob_dist / prob_dist.sum()\n    # Calculating entropy\n    entropy = -torch.sum(prob_dist * torch.log(prob_dist))\n    return entropy\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pyvips\nimport os\nimport cv2\nimport skimage\nfrom skimage.filters import threshold_otsu\nfrom matplotlib.colors import rgb_to_hsv\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nMAX_IMAGE_PIXELS = 10_000_000_000 \nImage.MAX_IMAGE_PIXELS = MAX_IMAGE_PIXELS\nfrom transformers import AutoImageProcessor, ViTModel\n\nfrom tqdm import tqdm\nimport random\nimport io\n\ndef tile_image(image, tile_width, tile_height, image_mask, resize_height_ratio, resize_width_ratio):\n    # 加载大图像，移除 access='sequential'\n\n    # 计算需要多少行列的tiles\n    cols = image.width // tile_width\n    rows = image.height // tile_height\n\n    feats = []\n    # 切割并保存tiles\n    count = 0\n    for row in range(rows):\n        for col in range(cols):\n            # 计算当前tile的位置\n            left = col * tile_width\n            top = row * tile_height\n            \n            x_center_resized = round((row + 0.5) * tile_width * resize_width_ratio)\n            y_center_resized = round((col + 0.5) * tile_height * resize_height_ratio)\n            if(image_mask[x_center_resized,y_center_resized] < 0.6):\n                continue\n            if (count >= 1000):\n                break\n\n            # 切割tile\n            tile = image.crop(left, top, tile_width, tile_height)\n            img_array = np.ndarray(buffer=tile.write_to_memory(),\n                                    dtype=np.uint8,\n                                    shape=[tile.height, tile.width, tile.bands])\n\n            # 将NumPy数组转换为PIL Image对象\n            pil_image = Image.fromarray(img_array)\n\n\n            inputs = image_processor(pil_image, return_tensors=\"pt\")\n            inputs[\"pixel_values\"] = inputs[\"pixel_values\"].cuda()\n\n            # get the features\n            with torch.no_grad():\n                outputs = model(**inputs)\n                features = outputs.last_hidden_state[:, 0, :]  # (1, 768) shape\n                feats.append(features)\n                count += 1\n        if(count >= 30000):\n            break\n    print(len(feats))\n    random.shuffle(feats)\n    # 取前1000个元素，如果列表长度小于1000，则取整个列表\n    feats = feats[: min(1000, len(feats))]\n    return feats\n\nTest = True\n\nif(Test):\n    DATA_DIR = '/kaggle/input/UBC-OCEAN'\n    IS_TRAINING = False\n\n    _CLS_MAPPING = {\n        'CC': 0,\n        'EC': 1, \n        'HGSC': 2, \n        'LGSC': 3, \n        'MC': 4,\n        'Other': 5,\n    }\n\n    _CLS_INV_MAPPING = {\n        0: 'CC', \n        1: 'EC', \n        2: 'HGSC', \n        3: 'LGSC', \n        4: 'MC',\n        5: 'Other',\n    }\n\n    _CLASSES = ['CC', 'EC', 'HGSC', 'LGSC', 'MC']\n    NUM_CLASSES = len(_CLASSES)\n    _DATA_DIR = Path(DATA_DIR)\n    FEATURE_DIR = '/kaggle/input/ubc-ocean-feature-set/Feature_set'\n    if not os.path.exists(FEATURE_DIR):\n        os.mkdir(FEATURE_DIR)\n\n    device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\n    if IS_TRAINING:\n        csv_path = _DATA_DIR / \"train.csv\"\n    else:\n        csv_path = _DATA_DIR / \"test.csv\"\n    df = pd.read_csv(csv_path)\n    imageid_list = df['image_id']\n\n    phikon_path = \"/kaggle/input/phikons-model/phikons\"\n\n    image_processor = AutoImageProcessor.from_pretrained(phikon_path)\n    model = ViTModel.from_pretrained(phikon_path, add_pooling_layer=False)\n    model = model.eval().cuda()\n    \n    submission = pd.DataFrame(columns=['image_id', 'label'])\n    \n    tile_width = 224\n    tile_height = 224\n\n    model_base_dir = '/kaggle/working'\n    model_mil = ABMIL(in_dim, 512, 128, num_classes)\n    model_mil.load_state_dict(torch.load(f'{model_base_dir}/final_wsi_vitp16_abmil_1.0_20ep.pth', map_location='cpu'))\n    model_mil = model_mil.to(device)\n    model_mil.eval()\n\n    for idx, image_id in tqdm(enumerate(df[\"image_id\"]), total=len(df[\"image_id\"])):\n        \n        raw_img_path = get_path_fullres_image(image_id, _DATA_DIR, False)\n        # thumbnail_path = get_path_thumbnail_image(image_id, DATA_DIR, True)\n\n        print(\"image_id :\", image_id)\n        image_mask = compute_matter_mask(image_id, _DATA_DIR)\n        print(\"thumbnail_shape: \",image_mask.shape)\n\n        thumbnail_height, thumbnail_width = image_mask.shape\n\n        raw_image = pyvips.Image.new_from_file(raw_img_path)\n        raw_height = raw_image.height\n        raw_width = raw_image.width\n        print(\"rawimage_shape: \",raw_height, raw_width)\n\n        resize_height_ratio = thumbnail_height / raw_height\n        resize_width_ratio = thumbnail_width / raw_width\n        feats = tile_image(raw_image, tile_width, tile_height, image_mask, resize_height_ratio, resize_width_ratio)\n\n        wsl_feat = torch.cat(feats, dim=0)\n        if wsl_feat.size(0) < 1000:\n            repeat_num = 1000 // wsl_feat.size(0)\n            pick_num = 1000 % wsl_feat.size(0)\n            picks = random.sample(range(0, wsl_feat.size(0)), pick_num)\n            wsl_feat = torch.cat([wsl_feat.repeat(repeat_num, 1), wsl_feat[picks, :]], dim=0)\n\n        with torch.no_grad():\n\n            wsl_feat = wsl_feat.squeeze(0).to(device)\n\n            scores = model_mil(wsl_feat)\n            scores = torch.softmax(scores, 1)\n            print(scores)\n            entropy = calculate_entropy(scores)\n            pred = torch.argmax(scores)\n\n            if(entropy.item() > 1.25):\n                submission.loc[len(submission)] = [image_id, _CLS_INV_MAPPING[5]]\n            else:\n                submission.loc[len(submission)] = [image_id, _CLS_INV_MAPPING[pred.item()]]\n\n        submission.to_csv('submission.csv', index=False)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}