{"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":"none","dataSources":[{"sourceId":94147,"databundleVersionId":11390004,"sourceType":"competition"}],"dockerImageVersionId":30918,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Prediction using CNN feature extraction\n\nThis notebook attempts to predict Cell Type abundance from the Image data.\n\nIt splits the image into patches, and uses CNN to extract feature from each. \n\nExtracing the CNN features takes a long time - about 40 minutes without GPU. We use pickle files to store it and load it quickly during development.","metadata":{}},{"cell_type":"code","source":"import h5py\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport matplotlib.pyplot as plt\nimport os\nimport pickle\nfrom matplotlib.patches import Rectangle\n\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.base import BaseEstimator, TransformerMixin\n\nfrom tensorflow.keras.applications import ResNet50\nfrom tensorflow.keras.applications.resnet50 import preprocess_input\nfrom tensorflow.keras.models import Model\n\n# New import for image registration:\nfrom skimage.registration import phase_cross_correlation\n\n# -----------------------------\n# 1. Functions for Patch Extraction and CNN Features\n# -----------------------------\ndef extract_patch(image, center, patch_size):\n    \"\"\"\n    Extract a square patch from the image centered at the given coordinate.\n    Assumes image shape is (height, width, channels) and center is (x, y).\n    \"\"\"\n    x, y = int(center[0]), int(center[1])\n    half_size = patch_size // 2\n    # Ensure indices are within bounds\n    y_min = max(y - half_size, 0)\n    y_max = min(y + half_size, image.shape[0])\n    x_min = max(x - half_size, 0)\n    x_max = min(x + half_size, image.shape[1])\n    patch = image[y_min:y_max, x_min:x_max, :]\n    return patch\n\ndef extract_cnn_features(patch, cnn_model):\n    \"\"\"\n    Resize, preprocess, and extract CNN features from a given image patch.\n    \n    Parameters:\n      patch (ndarray): The image patch to process.\n      cnn_model (Model): The pre-trained CNN model for feature extraction.\n      \n    Returns:\n      features (ndarray): Flattened feature vector from the CNN.\n    \"\"\"\n    # Resize patch to the input size expected by ResNet50 (e.g., 224x224)\n    patch_resized = cv2.resize(patch, (224, 224))\n    patch_preprocessed = preprocess_input(np.expand_dims(patch_resized, axis=0))\n    features = cnn_model.predict(patch_preprocessed, verbose=0)\n    return features.flatten()\n\n# -----------------------------\n# 2. Custom Transformer to Extract CNN Features from Patches\n# -----------------------------\nclass PatchFeatureExtractor(BaseEstimator, TransformerMixin):\n    def __init__(self, image, patch_size, cnn_model):\n        \"\"\"\n        Parameters:\n          image (ndarray): The whole-slide HE image as a numpy array.\n          patch_size (int): Size (in pixels) of the square patch to extract.\n          cnn_model (Model): Pre-trained CNN model for feature extraction.\n        \"\"\"\n        self.image = image\n        self.patch_size = patch_size\n        self.cnn_model = cnn_model\n\n    def fit(self, X, y=None):\n        return self\n\n    def transform(self, X):\n        # Collect all resized patches in a list\n        patches = []\n        for coord in X:\n            patch = extract_patch(self.image, coord, self.patch_size)\n            # Resize patch to the input size expected by ResNet50 (224x224)\n            patch_resized = cv2.resize(patch, (224, 224))\n            patches.append(patch_resized)\n        # Convert list to numpy array: shape (n_samples, 224, 224, channels)\n        patches = np.array(patches)\n        # Preprocess all patches at once using the appropriate preprocessing function\n        patches_preprocessed = preprocess_input(patches.astype(np.float32))\n        # Run the CNN model once on the entire batch\n        features = self.cnn_model.predict(patches_preprocessed, verbose=0)\n        # Flatten features if necessary (ResNet50 with pooling='avg' already outputs 2D arrays)\n        features = features.reshape(features.shape[0], -1)\n        return features\n        \n# -----------------------------\n# 3. Pipeline Class for the Elucidata Challenge with Caching and Visualization Options\n# -----------------------------\nclass CellTypePipeline:\n    \"\"\"\n    Pipeline for loading data, extracting image patch features using a CNN,\n    training a multi-output regression model, and generating a submission file.\n    \n    Optionally, CNN features can be cached (saved/loaded) using pickle to speed up re-runs.\n    Additional visualization methods are provided to verify that the spot coordinates \n    and extracted patches align with the HE slide image.\n    \"\"\"\n    \n    def __init__(self, h5_file_path, patch_size=64):\n        self.h5_file_path = h5_file_path\n        self.patch_size = patch_size\n        self.train_spot_tables = {}\n        self.train_images = {}\n        self.cell_type_columns = None\n        self.cnn_model = None  # To be initialized\n        self.feature_extractor_pipeline = None\n\n    def initialize_cnn_model(self):\n        \"\"\"\n        Initialize a pre-trained ResNet50 model (without top layers) for feature extraction.\n        \"\"\"\n        base_model = ResNet50(weights='imagenet', include_top=False, pooling='avg')\n        self.cnn_model = Model(inputs=base_model.input, outputs=base_model.output)\n        print(\"CNN feature extractor initialized.\")\n\n    def load_train_data(self):\n        \"\"\"\n        Load training spot data from the H5 file and store each slide as a DataFrame.\n        \"\"\"\n        with h5py.File(self.h5_file_path, \"r\") as f:\n            train_spots = f[\"spots/Train\"]\n            for slide_name in train_spots.keys():\n                spot_array = np.array(train_spots[slide_name])\n                df = pd.DataFrame(spot_array, columns=[\"x\", \"y\"] + [f\"C{i}\" for i in range(1, 36)])\n                self.train_spot_tables[slide_name] = df\n        print(\"Training spot data loaded successfully.\")\n        \n    def load_train_images(self):\n        \"\"\"\n        Load training HE images from the H5 file.\n        Adjust the key if your H5 file uses a different naming convention.\n        \"\"\"\n        with h5py.File(self.h5_file_path, \"r\") as f:\n            # Adjust key if necessary (e.g., f[\"images/Train\"] if needed)\n            train_imgs = f[\"images/Train\"]\n            for slide_name in train_imgs.keys():\n                image_array = np.array(train_imgs[slide_name])\n                self.train_images[slide_name] = image_array\n        print(\"Training images loaded successfully.\")\n\n    def load_test_data(self, slide_id):\n        \"\"\"\n        Load test spot data for a given slide.\n        \"\"\"\n        with h5py.File(self.h5_file_path, \"r\") as f:\n            test_spots = f[\"spots/Test\"]\n            if slide_id not in test_spots:\n                raise ValueError(f\"Slide {slide_id} not found in test spot data.\")\n            spot_array = np.array(test_spots[slide_id])\n            test_df = pd.DataFrame(spot_array, columns=[\"x\", \"y\"])\n        print(f\"Test spot data for slide {slide_id} loaded successfully.\")\n        return test_df\n\n    def load_test_image(self, slide_id):\n        \"\"\"\n        Load test HE image for a given slide.\n        \"\"\"\n        with h5py.File(self.h5_file_path, \"r\") as f:\n            test_imgs = f[\"images/Test\"]\n            if slide_id not in test_imgs:\n                raise ValueError(f\"Slide {slide_id} not found in test images.\")\n            image_array = np.array(test_imgs[slide_id])\n        print(f\"Test image for slide {slide_id} loaded successfully.\")\n        return image_array\n\n    def prepare_training_set(self, slide_id='S_1', cache_path=None):\n        \"\"\"\n        Prepare training features and targets for a given slide.\n        Uses the HE image to extract patches and then CNN features.\n        \n        If cache_path is provided and exists, the method will load cached features.\n        Otherwise, it will compute the features and then save them to the provided cache_path.\n        \"\"\"\n        if cache_path is not None and os.path.exists(cache_path):\n            print(f\"Loading cached training features from {cache_path} for slide {slide_id} ...\")\n            with open(cache_path, \"rb\") as f:\n                X_features, y = pickle.load(f)\n            return X_features, y\n        \n        if slide_id not in self.train_spot_tables:\n            raise ValueError(f\"Slide {slide_id} not found in training spot data.\")\n        if slide_id not in self.train_images:\n            raise ValueError(f\"Slide {slide_id} image not loaded.\")\n            \n        df = self.train_spot_tables[slide_id]\n        # Assume first two columns are coordinates and the rest are cell type abundances.\n        feature_cols = ['x', 'y']\n        target_cols = [col for col in df.columns if col not in feature_cols]\n        self.cell_type_columns = target_cols  # Assumed consistent across slides\n        \n        # Extract coordinates (for patch extraction)\n        X_coords = df[feature_cols].values.astype(float)\n        # Cell type abundance targets\n        y = df[target_cols].values.astype(float)\n        \n        # Build a feature extractor pipeline for this slide.\n        he_image = self.train_images[slide_id]\n        patch_extractor = PatchFeatureExtractor(he_image, self.patch_size, self.cnn_model)\n        self.feature_extractor_pipeline = Pipeline([\n            ('patch_extractor', patch_extractor),\n            ('scaler', StandardScaler())\n        ])\n        # Extract features for training\n        print(f\"Extracting CNN features for slide {slide_id} ...\")\n        X_features = self.feature_extractor_pipeline.fit_transform(X_coords)\n        \n        if cache_path is not None:\n            print(f\"Saving training features for slide {slide_id} to {cache_path} ...\")\n            with open(cache_path, \"wb\") as f:\n                pickle.dump((X_features, y), f)\n                \n        print(f\"Extracted CNN features for slide {slide_id}.\")\n        return X_features, y\n\n    def prepare_all_training_set(self, cache_dir=None, align_spots=True):\n        \"\"\"\n        Prepare training features and targets for all slides in the training set (S_1 to S_6).\n        If cache_dir is provided, each slide will use a separate cache file named train_features_<slide_id>.pkl.\n        If align_spots is True, the function computes the optimal alignment for each slide and adjusts\n        the spot coordinates accordingly.\n        \n        Returns:\n            X_all (ndarray): Concatenated features for all slides.\n            y_all (ndarray): Concatenated target values for all slides.\n        \"\"\"\n        X_list = []\n        y_list = []\n        # Process slides in sorted order (e.g., S_1, S_2, ... S_6)\n        for slide_id in sorted(self.train_spot_tables.keys()):\n            if align_spots:\n                # Compute the optimal alignment shift for the current slide\n                #optimal_shift, error, diffphase = self.compute_optimal_shift(slide_id, display=False)\n                #print(f\"For slide {slide_id}, applying optimal shift {optimal_shift} (error: {error:.4f}).\")\n                # Adjust the spot coordinates by subtracting the computed shift.\n                # (This assumes that the computed shift indicates how much the spots need to be shifted\n                # to align with the tissue mask.)\n                df = self.train_spot_tables[slide_id]\n                coords = df[['x', 'y']].values.astype(float)\n                #adjusted_coords = coords - optimal_shift  # shift the spots into alignment\n                #df[['x', 'y']] = adjusted_coords\n                self.train_spot_tables[slide_id] = df  # update the stored DataFrame\n    \n            slide_cache_path = os.path.join(cache_dir, f\"train_features_{slide_id}.pkl\") if cache_dir else None\n            X, y = self.prepare_training_set(slide_id=slide_id, cache_path=slide_cache_path)\n            X_list.append(X)\n            y_list.append(y)\n        X_all = np.concatenate(X_list, axis=0)\n        y_all = np.concatenate(y_list, axis=0)\n        print(\"All training features extracted and concatenated.\")\n        return X_all, y_all\n\n    def build_regression_pipeline(self):\n        \"\"\"\n        Build and return a regression pipeline that uses the pre-extracted CNN features.\n        \"\"\"\n        pipeline = Pipeline([\n            ('regressor', MultiOutputRegressor(RandomForestRegressor(n_estimators=100, random_state=42)))\n        ])\n        return pipeline\n\n    def train(self, X, y):\n        \"\"\"\n        Train the regression model on the provided features and targets.\n        \"\"\"\n        reg_pipeline = self.build_regression_pipeline()\n        reg_pipeline.fit(X, y)\n        print(\"Regression model training complete.\")\n        return reg_pipeline\n\n    def predict(self, reg_model, X_test):\n        \"\"\"\n        Predict cell type abundances on test features.\n        \"\"\"\n        predictions = reg_model.predict(X_test)\n        return predictions\n\n    def create_submission(self, test_df, predictions, submission_filename=\"submission.csv\"):\n        \"\"\"\n        Create a submission CSV file with predicted cell type abundances.\n        \"\"\"\n        pred_df = pd.DataFrame(predictions, columns=self.cell_type_columns, index=test_df.index)\n        pred_df.insert(0, 'ID', pred_df.index)\n        pred_df.to_csv(submission_filename, index=False)\n        print(f\"Submission file '{submission_filename}' created!\")\n\n    # -----------------------------\n    # Visualization Methods\n    # -----------------------------\n    def visualize_spot_overlay(self, slide_id, flip_y=False):\n        \"\"\"\n        Visualize the overlay of spot coordinates on the HE slide image.\n        Optionally, flip the y-axis if needed.\n        \"\"\"\n        if slide_id not in self.train_images or slide_id not in self.train_spot_tables:\n            raise ValueError(f\"Slide {slide_id} data not found.\")\n        image = self.train_images[slide_id]\n        df = self.train_spot_tables[slide_id]\n        coords = df[['x', 'y']].values.astype(float)\n        if flip_y:\n            coords[:, 1] = image.shape[0] - coords[:, 1]\n        plt.figure(figsize=(10, 10))\n        plt.imshow(image)\n        plt.scatter(coords[:, 0], coords[:, 1], marker='o', color='red', s=25)\n        plt.title(f\"Overlay of Spot Coordinates for Slide {slide_id}\")\n        plt.show()\n\n    def visualize_extracted_patches(self, slide_id, num_patches=5, flip_y=False):\n        \"\"\"\n        Visualize a few extracted patches from the slide to verify correct extraction.\n        \"\"\"\n        if slide_id not in self.train_images or slide_id not in self.train_spot_tables:\n            raise ValueError(f\"Slide {slide_id} data not found.\")\n        image = self.train_images[slide_id]\n        df = self.train_spot_tables[slide_id]\n        coords = df[['x', 'y']].values.astype(float)\n        if flip_y:\n            coords[:, 1] = image.shape[0] - coords[:, 1]\n        fig, axes = plt.subplots(1, num_patches, figsize=(num_patches * 3, 3))\n        for i in range(num_patches):\n            patch = extract_patch(image, coords[i], self.patch_size)\n            axes[i].imshow(patch)\n            axes[i].set_title(f\"Patch {i}\")\n            axes[i].axis(\"off\")\n        plt.suptitle(f\"Extracted Patches for Slide {slide_id}\")\n        plt.show()\n    \n    def visualize_cnn_input(self, slide_id, index=0, flip_y=False):\n        \"\"\"\n        Visualize the resized patch for a given spot index with an overlay of the spot position,\n        plot the distribution of the 35 cell type abundancies for that spot,\n        print the patch's top-left (X, Y) coordinates, and display the full slide image\n        with a rectangle marking the patch position.\n        \"\"\"\n        if slide_id not in self.train_images or slide_id not in self.train_spot_tables:\n            raise ValueError(f\"Slide {slide_id} data not found.\")\n            \n        # Get image and corresponding DataFrame\n        image = self.train_images[slide_id]\n        df = self.train_spot_tables[slide_id]\n        \n        # Get coordinates from the DataFrame\n        coords = df[['x', 'y']].values.astype(float)\n        if flip_y:\n            coords[:, 1] = image.shape[0] - coords[:, 1]\n        coord = coords[index]\n        \n        # Compute patch boundaries and extract patch\n        half_size = self.patch_size // 2\n        x = int(coord[0])\n        y = int(coord[1])\n        x_min = max(x - half_size, 0)\n        y_min = max(y - half_size, 0)\n        patch = extract_patch(image, coord, self.patch_size)\n        \n        # Resize patch to CNN input size (e.g., 224x224)\n        patch_resized = cv2.resize(patch, (224, 224))\n        \n        # Compute the spot's relative position within the patch\n        rel_x = x - x_min\n        rel_y = y - y_min\n        scale_x = 224 / patch.shape[1]\n        scale_y = 224 / patch.shape[0]\n        spot_resized_x = rel_x * scale_x\n        spot_resized_y = rel_y * scale_y\n    \n        # Retrieve cell type abundancies for this spot.\n        abundances = df.iloc[index][[col for col in df.columns if col not in ['x', 'y']]]\n        \n        # Print the patch's top-left coordinates\n        print(f\"Patch top-left coordinates: (x_min: {x_min}, y_min: {y_min})\")\n        \n        # Create three subplots: left for the full patch, middle for resized patch, right for the bar chart.\n        fig, (ax_full, ax_img, ax_bar) = plt.subplots(1, 3, figsize=(18, 5))\n\n        # Right subplot: Full slide image with a rectangle overlay marking the patch location.\n        ax_full.imshow(image)\n        # Draw a rectangle at (x_min, y_min) with width and height equal to patch_size.\n        rect = Rectangle((x_min, y_min), self.patch_size, self.patch_size, linewidth=2, edgecolor='red', facecolor='none')\n        ax_full.add_patch(rect)\n        ax_full.set_title(\"Full Slide with Patch Overlay\")\n        ax_full.axis(\"off\")\n        \n        # Middle subplot: Resized patch with spot overlay.\n        ax_img.imshow(patch_resized)\n        ax_img.scatter([spot_resized_x], [spot_resized_y], marker='x', color='red', s=50)\n        ax_img.set_title(f\"Resized Patch (Index {index})\")\n        # Annotate the patch with its top-left coordinates.\n        ax_img.text(5, 20, f\"({x_min}, {y_min})\", color='yellow', fontsize=12, \n                    bbox=dict(facecolor='black', alpha=0.5))\n        ax_img.axis(\"off\")\n        \n        # Left subplot: Bar chart of the cell type abundancies.\n        cell_types = abundances.index.tolist()  # e.g., ['C1', 'C2', ..., 'C35']\n        ax_bar.bar(cell_types, abundances.values)\n        ax_bar.set_title(\"Cell Type Abundance Distribution\")\n        ax_bar.set_xticklabels(cell_types, rotation=90)\n        ax_bar.set_ylabel(\"Abundance\")\n        \n        plt.tight_layout()\n        plt.show()\n            \n    def compute_optimal_shift(self, slide_id, flip_y=False, upsample_factor=10, spot_radius=3, display=False):\n        \"\"\"\n        Compute the optimal translational shift to align the spot coordinates with the tissue image.\n        \n        This function creates a binary tissue mask from the HE image (using Otsu thresholding)\n        and a corresponding binary spot mask by drawing filled circles at each spot coordinate.\n        It then computes the optimal (sub-pixel) shift between these masks using phase cross-correlation.\n        \n        Parameters:\n          slide_id (str): Identifier of the slide to process.\n          flip_y (bool): If True, flip the y-axis of spot coordinates.\n          upsample_factor (int): Upsampling factor for sub-pixel precision.\n          spot_radius (int): Radius (in pixels) for drawing spots in the mask.\n          display (bool): If True, display the tissue mask, spot mask, their overlay, and the HE image with corrected spots.\n          \n        Returns:\n          optimal_shift (ndarray): Optimal shift as (x_shift, y_shift).\n          error (float): Registration error.\n          diffphase (float): Diffphase value from phase_cross_correlation.\n        \"\"\"\n        if slide_id not in self.train_images or slide_id not in self.train_spot_tables:\n            raise ValueError(f\"Slide {slide_id} data not found.\")\n        \n        # Retrieve image and spot coordinates\n        image = self.train_images[slide_id]\n        df = self.train_spot_tables[slide_id]\n        coords = df[['x', 'y']].values.astype(float)\n        if flip_y:\n            coords[:, 1] = image.shape[0] - coords[:, 1]\n        \n        # Convert image to grayscale if needed\n        if len(image.shape) == 3 and image.shape[2] == 3:\n            image_gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)\n        else:\n            image_gray = image.copy()\n        \n        # Convert grayscale image to uint8 if needed (required for Otsu thresholding)\n        if image_gray.dtype != np.uint8:\n            image_gray = cv2.normalize(image_gray, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)\n                \n        # Generate tissue mask using Otsu thresholding\n        _, tissue_mask = cv2.threshold(image_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n        \n        # Create a blank mask for the spots (same dimensions as tissue_mask)\n        spot_mask = np.zeros_like(tissue_mask, dtype=np.uint8)\n        for pt in coords:\n            x, y = int(round(pt[0])), int(round(pt[1]))\n            cv2.circle(spot_mask, (x, y), radius=spot_radius, color=255, thickness=-1)\n        \n        # Compute the shift using phase cross-correlation\n        from skimage.registration import phase_cross_correlation\n        shift, error, diffphase = phase_cross_correlation(tissue_mask, spot_mask, upsample_factor=upsample_factor)\n        # shift is returned as (row_shift, col_shift) => (y, x)\n        optimal_shift = np.array([shift[1], shift[0]])  # Convert to (x_shift, y_shift)\n        print(\"Optimal shift (x, y):\", optimal_shift)\n        print(\"Registration error:\", error)\n        \n        if display:\n            # Display tissue mask, spot mask, and overlay\n            plt.figure(figsize=(15, 5))\n            plt.subplot(1, 3, 1)\n            plt.title(\"Tissue Mask\")\n            plt.imshow(tissue_mask, cmap='gray')\n            plt.subplot(1, 3, 2)\n            plt.title(\"Spot Mask\")\n            plt.imshow(spot_mask, cmap='gray')\n            plt.subplot(1, 3, 3)\n            plt.title(\"Overlay of Tissue and Spots\")\n            plt.imshow(tissue_mask, cmap='gray')\n            plt.imshow(spot_mask, cmap='jet', alpha=0.5)\n            plt.tight_layout()\n            plt.show()\n            \n            # Plot the original HE image with corrected spot coordinates\n            adjusted_coords = coords - optimal_shift\n            plt.figure(figsize=(10, 10))\n            plt.imshow(image)\n            plt.scatter(adjusted_coords[:, 0], adjusted_coords[:, 1], marker='o', color='lime', s=25)\n            plt.title(\"Corrected Spots Overlay on Tissue Image\")\n            plt.show()\n        \n        return optimal_shift, error, diffphase\n\n    def manual_shift_alignment(self, slide_id, x_shift, y_shift, flip_y=False, display=True):\n        \"\"\"\n        Manually apply a specified x,y shift to the spot coordinates for a given slide,\n        and plot the image with original and shifted spot overlays.\n        \n        Parameters:\n          slide_id (str): Identifier of the slide.\n          x_shift (float): Shift to apply in the x-direction (positive shifts right).\n          y_shift (float): Shift to apply in the y-direction (positive shifts down).\n          flip_y (bool): If True, flip the y-axis of the spot coordinates.\n          display (bool): If True, display the plots.\n          \n        Returns:\n          original_coords (ndarray): The original spot coordinates.\n          shifted_coords (ndarray): The adjusted spot coordinates.\n        \"\"\"\n        if slide_id not in self.train_images or slide_id not in self.train_spot_tables:\n            raise ValueError(f\"Slide {slide_id} data not found.\")\n        \n        # Retrieve the image and spot DataFrame\n        image = self.train_images[slide_id]\n        df = self.train_spot_tables[slide_id]\n        \n        # Extract original coordinates\n        original_coords = df[['x', 'y']].values.astype(float)\n        if flip_y:\n            original_coords[:, 1] = image.shape[0] - original_coords[:, 1]\n        \n        # Compute the shifted coordinates: here we add the shift.\n        # (Adjust as needed if your coordinate system requires subtracting the shift.)\n        shift_vector = np.array([x_shift, y_shift])\n        shifted_coords = original_coords + shift_vector\n        \n        if display:\n            plt.figure(figsize=(14, 7))\n            \n            # Plot original overlay\n            plt.subplot(1, 2, 1)\n            plt.imshow(image)\n            plt.scatter(original_coords[:, 0], original_coords[:, 1], \n                        marker='o', color='red', s=25, label='Original Spots')\n            plt.title(f\"Slide {slide_id} - Original Spots\")\n            plt.legend()\n            plt.axis(\"off\")\n            \n            # Plot shifted overlay\n            plt.subplot(1, 2, 2)\n            plt.imshow(image)\n            plt.scatter(shifted_coords[:, 0], shifted_coords[:, 1], \n                        marker='o', color='lime', s=25, label='Shifted Spots')\n            plt.title(f\"Slide {slide_id} - Shifted Spots\\n(x_shift: {x_shift}, y_shift: {y_shift})\")\n            plt.legend()\n            plt.axis(\"off\")\n            \n            plt.tight_layout()\n            plt.show()\n            \n        return original_coords, shifted_coords\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:07.811235Z","iopub.execute_input":"2025-04-02T20:18:07.811697Z","iopub.status.idle":"2025-04-02T20:18:07.866331Z","shell.execute_reply.started":"2025-04-02T20:18:07.811658Z","shell.execute_reply":"2025-04-02T20:18:07.865197Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# -----------------------------\n# 4. Example Usage with Caching and Visualization Options\n# -----------------------------\n\n# Path to the provided H5 data file\nh5_file_path = \"/kaggle/input/el-hackathon-2025/elucidata_ai_challenge_data.h5\"\n\n# Optionally specify a directory for caching training features (for slides S_1 to S_6)\ntrain_cache_dir = \"train_features_cache\"\nos.makedirs(train_cache_dir, exist_ok=True)\ntest_cache_path = \"test_features_S_7.pkl\"      # For slide S_7 test features\n\n# Initialize the pipeline with desired patch size (in pixels)\npipeline_obj = CellTypePipeline(h5_file_path, patch_size=64)\n\n# Initialize the CNN feature extractor (ResNet50)\npipeline_obj.initialize_cnn_model()\n\n# Load training spots and images\npipeline_obj.load_train_data()\npipeline_obj.load_train_images()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:07.867987Z","iopub.execute_input":"2025-04-02T20:18:07.868334Z","iopub.status.idle":"2025-04-02T20:18:11.135072Z","shell.execute_reply.started":"2025-04-02T20:18:07.868306Z","shell.execute_reply":"2025-04-02T20:18:11.134025Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Patches in Slide 1","metadata":{}},{"cell_type":"code","source":"\n\n# Optional: Visualize the overlay of spot coordinates on a training slide (e.g., S_1)\npipeline_obj.visualize_spot_overlay(slide_id='S_1', flip_y=False)\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:11.137356Z","iopub.execute_input":"2025-04-02T20:18:11.137756Z","iopub.status.idle":"2025-04-02T20:18:12.427164Z","shell.execute_reply.started":"2025-04-02T20:18:11.137727Z","shell.execute_reply":"2025-04-02T20:18:12.425962Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.compute_optimal_shift(slide_id=\"S_1\", display=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:12.428749Z","iopub.execute_input":"2025-04-02T20:18:12.429230Z","iopub.status.idle":"2025-04-02T20:18:16.132046Z","shell.execute_reply.started":"2025-04-02T20:18:12.429182Z","shell.execute_reply":"2025-04-02T20:18:16.130649Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.manual_shift_alignment(slide_id=\"S_1\", x_shift=-50, y_shift=-50)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:20:16.776146Z","iopub.execute_input":"2025-04-02T20:20:16.776486Z","iopub.status.idle":"2025-04-02T20:20:18.925571Z","shell.execute_reply.started":"2025-04-02T20:20:16.776458Z","shell.execute_reply":"2025-04-02T20:20:18.924041Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# Optional: Visualize a few extracted patches from a training slide (e.g., S_1)\npipeline_obj.visualize_extracted_patches(slide_id='S_1', num_patches=5, flip_y=False)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:18.282609Z","iopub.execute_input":"2025-04-02T20:18:18.282921Z","iopub.status.idle":"2025-04-02T20:18:19.065170Z","shell.execute_reply.started":"2025-04-02T20:18:18.282895Z","shell.execute_reply":"2025-04-02T20:18:19.063822Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Optional: Visualize the CNN input (resized patch) for a specific spot (e.g., index 0 from S_1)\npipeline_obj.visualize_cnn_input(slide_id='S_1', index=0, flip_y=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:19.066262Z","iopub.execute_input":"2025-04-02T20:18:19.066634Z","iopub.status.idle":"2025-04-02T20:18:20.606820Z","shell.execute_reply.started":"2025-04-02T20:18:19.066603Z","shell.execute_reply":"2025-04-02T20:18:20.605622Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Optional: Visualize the CNN input (resized patch) for a specific spot (e.g., index 0 from S_1)\npipeline_obj.visualize_cnn_input(slide_id='S_1', index=1, flip_y=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:20.607880Z","iopub.execute_input":"2025-04-02T20:18:20.608193Z","iopub.status.idle":"2025-04-02T20:18:22.115353Z","shell.execute_reply.started":"2025-04-02T20:18:20.608165Z","shell.execute_reply":"2025-04-02T20:18:22.114399Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.visualize_cnn_input(slide_id='S_1', index=10, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:22.118310Z","iopub.execute_input":"2025-04-02T20:18:22.118636Z","iopub.status.idle":"2025-04-02T20:18:23.659749Z","shell.execute_reply.started":"2025-04-02T20:18:22.118609Z","shell.execute_reply":"2025-04-02T20:18:23.658396Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Patches in Slide 2","metadata":{}},{"cell_type":"code","source":"pipeline_obj.visualize_spot_overlay(slide_id='S_2', flip_y=False)\npipeline_obj.visualize_extracted_patches(slide_id='S_2', num_patches=5, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:23.661927Z","iopub.execute_input":"2025-04-02T20:18:23.662286Z","iopub.status.idle":"2025-04-02T20:18:25.269695Z","shell.execute_reply.started":"2025-04-02T20:18:23.662255Z","shell.execute_reply":"2025-04-02T20:18:25.268556Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.compute_optimal_shift(slide_id=\"S_2\", display=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:25.270932Z","iopub.execute_input":"2025-04-02T20:18:25.271312Z","iopub.status.idle":"2025-04-02T20:18:29.071708Z","shell.execute_reply.started":"2025-04-02T20:18:25.271276Z","shell.execute_reply":"2025-04-02T20:18:29.070540Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.manual_shift_alignment(slide_id=\"S_2\", x_shift=-60, y_shift=-60)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:20:57.790014Z","iopub.execute_input":"2025-04-02T20:20:57.790393Z","iopub.status.idle":"2025-04-02T20:20:59.936935Z","shell.execute_reply.started":"2025-04-02T20:20:57.790360Z","shell.execute_reply":"2025-04-02T20:20:59.935576Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.visualize_cnn_input(slide_id='S_2', index=0, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:29.073072Z","iopub.execute_input":"2025-04-02T20:18:29.073406Z","iopub.status.idle":"2025-04-02T20:18:30.637673Z","shell.execute_reply.started":"2025-04-02T20:18:29.073375Z","shell.execute_reply":"2025-04-02T20:18:30.636136Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.visualize_cnn_input(slide_id='S_3', index=0, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:30.638926Z","iopub.execute_input":"2025-04-02T20:18:30.639278Z","iopub.status.idle":"2025-04-02T20:18:32.188630Z","shell.execute_reply.started":"2025-04-02T20:18:30.639247Z","shell.execute_reply":"2025-04-02T20:18:32.187312Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.visualize_cnn_input(slide_id='S_4', index=0, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:32.189699Z","iopub.execute_input":"2025-04-02T20:18:32.189971Z","iopub.status.idle":"2025-04-02T20:18:33.749209Z","shell.execute_reply.started":"2025-04-02T20:18:32.189944Z","shell.execute_reply":"2025-04-02T20:18:33.747978Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.visualize_cnn_input(slide_id='S_5', index=0, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:33.750363Z","iopub.execute_input":"2025-04-02T20:18:33.750757Z","iopub.status.idle":"2025-04-02T20:18:35.647985Z","shell.execute_reply.started":"2025-04-02T20:18:33.750718Z","shell.execute_reply":"2025-04-02T20:18:35.646567Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Slide S_3","metadata":{}},{"cell_type":"code","source":"pipeline_obj.visualize_spot_overlay(slide_id='S_2', flip_y=False)\npipeline_obj.visualize_extracted_patches(slide_id='S_2', num_patches=5, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:35.649138Z","iopub.execute_input":"2025-04-02T20:18:35.649612Z","iopub.status.idle":"2025-04-02T20:18:37.283630Z","shell.execute_reply.started":"2025-04-02T20:18:35.649555Z","shell.execute_reply":"2025-04-02T20:18:37.282456Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.compute_optimal_shift(slide_id=\"S_3\", display=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:37.284657Z","iopub.execute_input":"2025-04-02T20:18:37.284930Z","iopub.status.idle":"2025-04-02T20:18:41.144601Z","shell.execute_reply.started":"2025-04-02T20:18:37.284906Z","shell.execute_reply":"2025-04-02T20:18:41.143680Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.compute_optimal_shift(slide_id=\"S_4\", display=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:41.145443Z","iopub.execute_input":"2025-04-02T20:18:41.145763Z","iopub.status.idle":"2025-04-02T20:18:45.102328Z","shell.execute_reply.started":"2025-04-02T20:18:41.145735Z","shell.execute_reply":"2025-04-02T20:18:45.101262Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.compute_optimal_shift(slide_id=\"S_5\", display=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:45.103591Z","iopub.execute_input":"2025-04-02T20:18:45.103915Z","iopub.status.idle":"2025-04-02T20:18:48.961874Z","shell.execute_reply.started":"2025-04-02T20:18:45.103889Z","shell.execute_reply":"2025-04-02T20:18:48.960813Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pipeline_obj.compute_optimal_shift(slide_id=\"S_6\", display=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:48.963014Z","iopub.execute_input":"2025-04-02T20:18:48.963296Z","iopub.status.idle":"2025-04-02T20:18:52.682806Z","shell.execute_reply.started":"2025-04-02T20:18:48.963270Z","shell.execute_reply":"2025-04-02T20:18:52.681741Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#pipeline_obj.visualize_cnn_input(slide_id='S_7', index=0, flip_y=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:22:06.568570Z","iopub.execute_input":"2025-04-02T20:22:06.568946Z","iopub.status.idle":"2025-04-02T20:22:06.610568Z","shell.execute_reply.started":"2025-04-02T20:22:06.568917Z","shell.execute_reply":"2025-04-02T20:22:06.608705Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Training","metadata":{}},{"cell_type":"code","source":"skip_training = False\n\nif not skip_training:\n    # Prepare training features and targets from all slides (S_1 to S_6)\n    X_train, y_train = pipeline_obj.prepare_all_training_set(cache_dir=train_cache_dir)\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:52.683864Z","iopub.execute_input":"2025-04-02T20:18:52.684154Z","iopub.status.idle":"2025-04-02T20:18:52.688799Z","shell.execute_reply.started":"2025-04-02T20:18:52.684128Z","shell.execute_reply":"2025-04-02T20:18:52.687648Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if not skip_training:\n    X_train","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:52.690039Z","iopub.execute_input":"2025-04-02T20:18:52.690338Z","iopub.status.idle":"2025-04-02T20:18:52.708013Z","shell.execute_reply.started":"2025-04-02T20:18:52.690314Z","shell.execute_reply":"2025-04-02T20:18:52.706821Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport pickle\nimport numpy as np\nimport cv2\nimport matplotlib.pyplot as plt\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.multioutput import MultiOutputClassifier\n\n\ndef convert_to_rank(arr):\n    \"\"\"\n    Converts each row in the array to a rank order.\n    The highest value gets rank 1, the second highest gets rank 2, and so on.\n    \n    Parameters:\n      arr (numpy.ndarray): 2D array with shape (n_samples, 35) of absolute abundances.\n      \n    Returns:\n      numpy.ndarray: 2D array of the same shape with rank orders.\n    \"\"\"\n    return -1 * np.argsort(-arr, axis=1) + 1\n\n\n\nif not skip_training:\n\n    # Prepare training features and targets (absolute abundances) from slides S_1 to S_6\n    X_train, y_train = pipeline_obj.prepare_all_training_set(cache_dir=train_cache_dir)\n    \n    # Convert the absolute abundance targets to rank order (1 = highest abundance)\n    y_train_rank = convert_to_rank(y_train)\n    \n    # Convert the rank order to quartile bins:\n    # Ranks 1-5   -> ceil(rank/5) = 1, so quartile = 8 - 1 = 7 (top abundant)\n    # Ranks 6-10  -> ceil(rank/5) = 2, so quartile = 8 - 2 = 6\n    # ...\n    # Ranks 31-35 -> ceil(rank/5) = 7, so quartile = 8 - 7 = 1 (least abundant)\n    y_train_quartile = (8 - np.ceil(y_train_rank / 5)).astype(int)\n    \n    # Train a model on the extracted CNN features with quartile targets\n    reg_model = pipeline_obj.train(X_train, y_train_quartile)\n    \n    # Load test data and image for slide S_7 (as per challenge description)\n    test_df = pipeline_obj.load_test_data(slide_id='S_7')\n    test_image = pipeline_obj.load_test_image(slide_id='S_7')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:52.709115Z","iopub.execute_input":"2025-04-02T20:18:52.709430Z","iopub.status.idle":"2025-04-02T20:18:52.726662Z","shell.execute_reply.started":"2025-04-02T20:18:52.709386Z","shell.execute_reply":"2025-04-02T20:18:52.725459Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if not skip_training:\n    # Build a feature extractor for test slide using its HE image\n    test_patch_extractor = PatchFeatureExtractor(test_image, pipeline_obj.patch_size, pipeline_obj.cnn_model)\n    test_feature_pipeline = Pipeline([\n        ('patch_extractor', test_patch_extractor),\n        ('scaler', StandardScaler())\n    ])\n    X_test_coords = test_df[['x', 'y']].values.astype(float)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:52.727741Z","iopub.execute_input":"2025-04-02T20:18:52.728031Z","iopub.status.idle":"2025-04-02T20:18:52.746670Z","shell.execute_reply.started":"2025-04-02T20:18:52.727993Z","shell.execute_reply":"2025-04-02T20:18:52.745588Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if not skip_training:\n    # Check for cached test features\n    if os.path.exists(test_cache_path):\n        print(f\"Loading cached test features from {test_cache_path} ...\")\n        with open(test_cache_path, \"rb\") as f:\n            X_test_features = pickle.load(f)\n    else:\n        print(\"Extracting CNN features for test data ...\")\n        X_test_features = test_feature_pipeline.fit_transform(X_test_coords)\n        print(f\"Saving test features to {test_cache_path} ...\")\n        with open(test_cache_path, \"wb\") as f:\n            pickle.dump(X_test_features, f)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:52.747662Z","iopub.execute_input":"2025-04-02T20:18:52.747986Z","iopub.status.idle":"2025-04-02T20:18:52.764586Z","shell.execute_reply.started":"2025-04-02T20:18:52.747958Z","shell.execute_reply":"2025-04-02T20:18:52.763653Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if not skip_training:\n\n    # Predict cell type abundances for test data\n    predictions = pipeline_obj.predict(reg_model, X_test_features)\n    \n    # Create submission file\n    pipeline_obj.create_submission(test_df, predictions, submission_filename=\"submission.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-02T20:18:52.765664Z","iopub.execute_input":"2025-04-02T20:18:52.765926Z","iopub.status.idle":"2025-04-02T20:18:52.784585Z","shell.execute_reply.started":"2025-04-02T20:18:52.765903Z","shell.execute_reply":"2025-04-02T20:18:52.783627Z"}},"outputs":[],"execution_count":null}]}