{"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.9.19"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":23823,"databundleVersionId":1920183,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"b522f4c9-8b86-4415-abc4-0a0f4efa7336","cell_type":"markdown","source":"# Trying to Cluster Human Atlas Data","metadata":{}},{"id":"338f85e8-a518-4eb2-9ea0-f2ad5873abed","cell_type":"markdown","source":"# Clustering Failure Analysis on Human Protein Atlas Data\n\n## **Project Goal**\nI attempted to use **ResNet50 feature embeddings** to cluster images from the Human Protein Atlas dataset. The link to the dataset is in section 7 of this file.\n\nThe idea was to:  \n1. **Cluster the images** using **DPGMM (Dirichlet Process Gaussian Mixture Model)** and **K-Means**.  \n2. **Train separate neural networks** for each cluster, specializing in different image types.  \n3. **Use the predicted cluster** at inference time to send each image to its corresponding specialized neural net.  \n\nThe hypothesis was that **specialized networks** would outperform a single model trained on the entire dataset.  \n\n---\n\n## **What Went Wrong?**\n### **1️ Clustering Completely Failed**\n- **DPGMM failed to converge**, with ELBO exploding to massive values.\n- **K-Means produced clusters with near-zero silhouette scores** (no meaningful separations).\n- **DBSCAN classified nearly everything as noise (-1)** (no dense groups detected).  \n\n### **2️ The Root Cause: The Data Was Not Clusterable**\n- **Shapiro-Wilk test showed extreme non-Gaussianity** in feature distributions.  \n- **Skewness & Kurtosis confirmed heavy tails and non-normal behavior.**  \n- **PCA visualizations showed no natural separations.**  \n- **Ultimately, the dataset does not have meaningful clusters for unsupervised learning.**  \n\n---\n\n## **Key Takeaways: How to Detect Clustering Failures Early**\nThrough this process, I developed a **quick pipeline** to detect when clustering is a bad idea:  \n\n### **Fast Clustering Failure Detection Workflow**\n1 **Check Gaussianity First**  \n   - Run **Shapiro-Wilk test**  \n   - Check **Skewness & Kurtosis**  \n   - If heavily non-Gaussian → **Avoid GMM-based methods**  \n\n2️ **Run K-Means as a Baseline**  \n   - Compute **Silhouette Score**  \n   - If score **< 0.2**, data is likely **not separable**  \n   - If **Elbow Method is flat**, there’s no real structure  \n\n3️ **Try DBSCAN for Density-Based Clustering**  \n   - If **everything gets labeled as noise (-1)**, the dataset lacks natural clusters  \n\n---\n\n## **Conclusion: Clustering is NOT a Viable Approach for This Dataset**\n- The **dataset lacks natural separable clusters**, making **unsupervised learning infeasible**.  \n- Instead, **supervised learning with CNNs or Vision Transformers is the correct approach**.  \n- This analysis provides a **reliable workflow for detecting bad clustering problems early**, preventing wasted time on doomed models.  \n\n---\n\n**This project documents a real ML failure and provides a structured approach for debugging similar issues in the future.**  \n\n---","metadata":{}},{"id":"eb231699-36bc-49c2-b9e6-6dd2f18f7d3b","cell_type":"markdown","source":"# Contents:\n\n- Section 1: Preparing the File\n- Section 2: Test Section\n- Section 3: Code for Entire Dataset\n- Section 4: DPGMM Tests\n- Section 5: K-Means\n- Section 6: DBSCAN\n- Section 7: Citations","metadata":{}},{"id":"0203d1a4-e000-40de-8be5-1e9c1c58d1e0","cell_type":"markdown","source":"# Section 1: Preparing the File\n\n## Install all required libraries for the program","metadata":{}},{"id":"c9b3088d-1826-4532-b589-49246f27ba02","cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport cv2\nimport tifffile as tiff  # For loading .tif images\nimport seaborn as sns\nimport torch\nimport torchvision.transforms as transforms\nimport torchvision.models as models\nfrom tqdm import tqdm  # Import progress bar\nfrom sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.mixture import BayesianGaussianMixture\nimport pyro\nimport pyro.distributions as dist\nfrom pyro.infer import SVI, TraceMeanField_ELBO\nfrom pyro.optim import Adam\nimport pyro.optim as optim\nimport time\nfrom scipy.stats import shapiro, skew, kurtosis\nfrom sklearn.cluster import KMeans \nfrom sklearn.metrics import silhouette_score\nimport numpy as np\nfrom sklearn.metrics import silhouette_score  \nfrom sklearn.cluster import DBSCAN  ","metadata":{},"outputs":[],"execution_count":null},{"id":"028391db-2556-41a7-bb78-16db62595939","cell_type":"markdown","source":"## Give Paths","metadata":{}},{"id":"896020d3-d37b-44e6-ba73-3ef4ffd123bb","cell_type":"markdown","source":"All this does is give the paths to the folder containing the training images and the corresponding training labels.","metadata":{}},{"id":"6f68b3cf-453f-4b3c-9d78-c37f401a8014","cell_type":"code","source":"\n# Define main dataset path\nDATASET_PATH = r\"D:\\hpa-single-cell-image-classification\"\n\n# Define sub-paths\nTRAIN_PATH = os.path.join(DATASET_PATH, \"train\")  # Path to training images\nLABELS_PATH = os.path.join(DATASET_PATH, \"train.csv\")  # Path to labels CSV\n\n# Print to verify\nprint(\"Train Path:\", TRAIN_PATH)\nprint(\"Labels Path:\", LABELS_PATH)","metadata":{},"outputs":[],"execution_count":null},{"id":"2e97a846-c916-483d-ae13-177d17884e22","cell_type":"markdown","source":"---\n\n## Preparing the Dataset","metadata":{}},{"id":"2abf34a1-32aa-4f9f-a846-a0f88557be78","cell_type":"markdown","source":"## Load the Labels","metadata":{}},{"id":"7f57b321-45ca-4895-872b-c373b6f003fa","cell_type":"code","source":"df = pd.read_csv(LABELS_PATH)\nprint(df.head())  # Check first few rows","metadata":{},"outputs":[],"execution_count":null},{"id":"92508d34-92e1-4bd6-9192-e336a26be6aa","cell_type":"markdown","source":"## Loading a single image","metadata":{}},{"id":"e4d799cd-e8d6-4a2e-ae00-c43f1af45852","cell_type":"code","source":"# Get first image ID from the CSV\nsample_id = df.iloc[0][\"ID\"]\n\n# Load each channel (CHANGE .tif to .png)\nimg_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{sample_id}_blue.png\"), cv2.IMREAD_UNCHANGED)\nimg_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{sample_id}_green.png\"), cv2.IMREAD_UNCHANGED)  # Protein of interest\nimg_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{sample_id}_red.png\"), cv2.IMREAD_UNCHANGED)\nimg_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{sample_id}_yellow.png\"), cv2.IMREAD_UNCHANGED)\n\n# Convert images to RGB format for correct visualization\nimg_blue = cv2.cvtColor(img_blue, cv2.COLOR_BGR2RGB)\nimg_green = cv2.cvtColor(img_green, cv2.COLOR_BGR2RGB)\nimg_red = cv2.cvtColor(img_red, cv2.COLOR_BGR2RGB)\nimg_yellow = cv2.cvtColor(img_yellow, cv2.COLOR_BGR2RGB)\n\n# Display images\nfig, axes = plt.subplots(1, 4, figsize=(20, 5))\naxes[0].imshow(img_blue); axes[0].set_title(\"Nucleus (Blue)\")\naxes[1].imshow(img_green); axes[1].set_title(\"Protein of Interest (Green)\")\naxes[2].imshow(img_red); axes[2].set_title(\"Microtubules (Red)\")\naxes[3].imshow(img_yellow); axes[3].set_title(\"Endoplasmic Reticulum (Yellow)\")\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"437b43bb-0eb7-4cdc-9f38-d185a0d51d8c","cell_type":"code","source":"print(\"Blue Channel Shape:\", img_blue.shape)\nprint(\"Green Channel Shape:\", img_green.shape)  # This is the protein signal\nprint(\"Red Channel Shape:\", img_red.shape)\nprint(\"Yellow Channel Shape:\", img_yellow.shape)\n\n# Check pixel range\nprint(\"Min Pixel Value:\", np.min(img_green))\nprint(\"Max Pixel Value:\", np.max(img_green))","metadata":{},"outputs":[],"execution_count":null},{"id":"52300981-5261-4301-be7b-37244c430cea","cell_type":"markdown","source":"I will now try converting the images to grayscale and also kleep color images. We will then see which does better. \n\n---","metadata":{}},{"id":"f1691a82-5791-44db-9ead-deacd47e2bdc","cell_type":"markdown","source":"## 1A: Non Grayscale","metadata":{}},{"id":"83909dbd-022d-4c9d-836e-6ac54bab96b8","cell_type":"code","source":"# Stack original RGB images into a single 12-channel tensor\nimg_stacked_rgb = np.concatenate([img_blue, img_green, img_red, img_yellow], axis=-1)\nprint(\"Non-Grayscale Shape:\", img_stacked_rgb.shape)  # Should be (2048, 2048, 12)","metadata":{},"outputs":[],"execution_count":null},{"id":"6b59e439-44ac-4d1d-96f6-45bdf4a0e747","cell_type":"markdown","source":"## 1B: Grayscale","metadata":{}},{"id":"eab9601e-88be-4db3-a178-c14329a5a3ca","cell_type":"code","source":"# Convert each channel to grayscale\nimg_blue_gray = cv2.cvtColor(img_blue, cv2.COLOR_RGB2GRAY)\nimg_green_gray = cv2.cvtColor(img_green, cv2.COLOR_RGB2GRAY)\nimg_red_gray = cv2.cvtColor(img_red, cv2.COLOR_RGB2GRAY)\nimg_yellow_gray = cv2.cvtColor(img_yellow, cv2.COLOR_RGB2GRAY)\n\n# Stack grayscale images into a single 4-channel tensor\nimg_stacked_gray = np.stack([img_blue_gray, img_green_gray, img_red_gray, img_yellow_gray], axis=-1)\nprint(\"Grayscale Shape:\", img_stacked_gray.shape)  # Should be (2048, 2048, 4)","metadata":{},"outputs":[],"execution_count":null},{"id":"c9548245-dd7c-4dad-b8de-22da227c2d95","cell_type":"markdown","source":"I will now use a CNN to get features of the pictures to run DP-GMM on them to cluster them. ","metadata":{}},{"id":"df31144c-631c-49af-814b-0c00d3d78b8c","cell_type":"markdown","source":"---\n# Section 2: Test Section","metadata":{}},{"id":"de029356-6f4a-47bf-92ae-0a4815f0b82d","cell_type":"markdown","source":"Here I tested the code to make sure it works on only 10 images. If it failed it would not have taken an hour to run. I tested for the Green, the Multicolor, and the Grayscale images. Later, this same code would run on the entire dataset.","metadata":{}},{"id":"c7468cfc-5a15-4972-b409-b52e527d653a","cell_type":"markdown","source":"## Green test","metadata":{}},{"id":"be3ad9fa-d97f-43cf-8e08-b1454d27e6e0","cell_type":"code","source":"# Load pretrained ResNet model and move to GPU\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nmodel = models.resnet50(weights=\"IMAGENET1K_V1\").to(device)\nmodel = torch.nn.Sequential(*(list(model.children())[:-1]))  # Remove last FC layer\nmodel.eval()\n\n# Define transform (resize & normalize for CNN)\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((224, 224)),\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.5], std=[0.5])  \n])\n\ndef extract_cnn_features_green(img):\n    \"\"\"Extract features from the Green channel (Protein of Interest) using ResNet.\"\"\"\n    img_green = img[:, :, 1]  # ✅ Extract Green channel (index 1)\n\n    # Convert single-channel grayscale to 3-channel by stacking\n    img_green_3ch = np.stack([img_green, img_green, img_green], axis=-1)  # (H, W, 3)\n\n    img_tensor = transform(img_green_3ch).unsqueeze(0).to(device)  # ✅ Move tensor to GPU\n    with torch.no_grad():\n        features = model(img_tensor)  # ✅ Model is on GPU\n    return features.squeeze().cpu().numpy()  # ✅ Move result back to CPU for NumPy\n\n# Directory where images are stored\nTRAIN_PATH = \"D:/hpa-single-cell-image-classification/train\"\nSAVE_PATH = \"test_green_channel_features.npy\"\n\n# Get the first 10 images for testing\nall_files = [f for f in os.listdir(TRAIN_PATH) if \"_green.png\" in f][:10]\nnum_images = len(all_files)\nprint(f\"Processing {num_images} Green channel images...\")\n\n# Initialize or resume from existing file\nif os.path.exists(SAVE_PATH):\n    saved_features = np.load(SAVE_PATH)\n    start_index = len(saved_features)\n    all_features_green = list(saved_features)  # Convert back to list\n    print(f\"Resuming from index {start_index}/{num_images}...\")\nelse:\n    all_features_green = []\n    start_index = 0\n\n# Loop through the first 10 images\nfor i, filename in enumerate(tqdm(all_files[start_index:], desc=\"Extracting CNN Features\", unit=\"image\")):\n    base_id = filename.replace(\"_green.png\", \"\")\n\n    # Load all four grayscale images (single-channel)\n    img_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_blue.png\"), cv2.IMREAD_GRAYSCALE)\n    img_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_green.png\"), cv2.IMREAD_GRAYSCALE)\n    img_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_red.png\"), cv2.IMREAD_GRAYSCALE)\n    img_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_yellow.png\"), cv2.IMREAD_GRAYSCALE)\n\n    # Stack into (H, W, 4)\n    img_stacked_rgb = np.stack([img_blue, img_green, img_red, img_yellow], axis=-1)\n\n    # Extract features for this image\n    features = extract_cnn_features_green(img_stacked_rgb)\n\n    # Store in the list\n    all_features_green.append(features)\n\n    # Save after every image (since we only have 10, we save frequently)\n    np.save(SAVE_PATH, np.array(all_features_green))\n    print(f\"Progress saved at {i + start_index + 1}/{num_images} images.\")\n\n# Final save\nnp.save(SAVE_PATH, np.array(all_features_green))\nprint(\"\\nFeature extraction complete!\")\nprint(\"Final Shape of Extracted Features (Green-Only):\", np.array(all_features_green).shape)\nprint(f\"Features saved successfully to '{SAVE_PATH}'\")","metadata":{},"outputs":[],"execution_count":null},{"id":"d3221df7-10e3-402a-a663-f9f129e3c817","cell_type":"markdown","source":"Making sure the output is the correct shape. ","metadata":{}},{"id":"7892c6a5-5f53-4b88-a060-b9f6559ba6be","cell_type":"code","source":"# Load the saved features\nfeatures_green = np.load(\"test_green_channel_features.npy\")\n\n# Print shape\nprint(\"Loaded Feature Shape:\", features_green.shape)  # Should be (10, 2048)\n\n# Confirm we have exactly 10 vectors\nif features_green.shape[0] == 10:\n    print(\"✅ The file contains 10 feature vectors.\")\nelse:\n    print(\"❌ Something is wrong! The file has\", features_green.shape[0], \"vectors.\")","metadata":{},"outputs":[],"execution_count":null},{"id":"35285db3-6d2c-41fb-bc7f-be8ef062b9cd","cell_type":"markdown","source":"Making sure the feature vectors are unique. I was checking to see if the model was saving correctly. ","metadata":{}},{"id":"fd545932-5915-4c0b-ad26-b8d7edb737b0","cell_type":"code","source":"# Load features\nfeatures_green = np.load(\"test_green_channel_features.npy\")\n\n# Compute differences between consecutive vectors\ndiffs = np.diff(features_green, axis=0)\n\n# Count how many vectors are identical\nidentical_vectors = np.sum(np.all(diffs == 0, axis=1))\n\nif identical_vectors == 0:\n    print(\"✅ All feature vectors are unique.\")\nelse:\n    print(f\"❌ {identical_vectors} feature vectors are duplicates!\")","metadata":{},"outputs":[],"execution_count":null},{"id":"00c90918-f312-4e8d-b0f6-71c8f50871c8","cell_type":"markdown","source":"A visual inspection of some of the vectors.","metadata":{}},{"id":"8d1782d8-f0b2-493d-b190-b17576234182","cell_type":"code","source":"for i in range(3):  # Print first 3 vectors\n    print(f\"\\nFeature Vector {i + 1} (First 10 values):\", features_green[i][:10])","metadata":{},"outputs":[],"execution_count":null},{"id":"0c23af3b-f0cc-4550-b1d7-323b7a624096","cell_type":"markdown","source":"## Multicolor Test","metadata":{}},{"id":"db5de40e-59ca-48b7-884b-b075c65933b8","cell_type":"code","source":"# Load pretrained ResNet model and move to GPU\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nmodel = models.resnet50(weights=\"IMAGENET1K_V1\").to(device)\nmodel = torch.nn.Sequential(*(list(model.children())[:-1]))  # Remove last FC layer\nmodel.eval()\n\n# Define transform (normalize and resize for CNN)\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((224, 224)),  # Resize to match CNN input size\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.5], std=[0.5])  # Normalize intensity\n])\n\ndef create_pseudo_rgb(img):\n    \"\"\"Convert 4-channel image into a 3-channel (R-G-B) image.\"\"\"\n    r = img[:, :, 2]  # Use Red channel (Microtubules)\n    g = img[:, :, 3]  # Use Green channel (Protein of Interest)\n    b = img[:, :, 0]  # Use Blue channel (Nucleus)\n    \n    pseudo_rgb = np.stack([r, g, b], axis=-1)  # Stack into (H, W, 3)\n    return pseudo_rgb.astype(np.uint8)  # Ensure proper dtype\n\ndef extract_cnn_features_pseudo_rgb(img):\n    \"\"\"Extract features from a Pseudo-RGB image using ResNet.\"\"\"\n    img_tensor = transform(img).unsqueeze(0).to(device)  # ✅ Move tensor to GPU\n    with torch.no_grad():\n        features = model(img_tensor)  # ✅ Model is on GPU\n    return features.squeeze().cpu().numpy()  # ✅ Move result back to CPU for NumPy\n\n# Directory where images are stored\nTRAIN_PATH = \"D:/hpa-single-cell-image-classification/train\"\nSAVE_PATH = \"test_pseudo_rgb_features.npy\"\n\n# Get list of first 10 images for testing\nall_files = [f for f in os.listdir(TRAIN_PATH) if \"_green.png\" in f][:10]\nnum_images = len(all_files)\nprint(f\"Processing {num_images} Pseudo-RGB images...\\n\")\n\n# Initialize or resume from existing file\nif os.path.exists(SAVE_PATH):\n    saved_features = np.load(SAVE_PATH)\n    start_index = len(saved_features)\n    all_features_pseudo_rgb = list(saved_features)  # Convert back to list\n    print(f\"Resuming from index {start_index}/{num_images}...\")\nelse:\n    all_features_pseudo_rgb = []\n    start_index = 0\n\n# Loop through first 10 images\nfor i, filename in enumerate(tqdm(all_files[start_index:], desc=\"Extracting CNN Features (Pseudo-RGB)\", unit=\"image\")):\n    base_id = filename.replace(\"_green.png\", \"\")\n\n    # Load all four images as grayscale (single-channel)\n    img_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_blue.png\"), cv2.IMREAD_GRAYSCALE)\n    img_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_green.png\"), cv2.IMREAD_GRAYSCALE)\n    img_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_red.png\"), cv2.IMREAD_GRAYSCALE)\n    img_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_yellow.png\"), cv2.IMREAD_GRAYSCALE)\n\n    # Stack into (H, W, 4)\n    img_stacked_rgb = np.stack([img_blue, img_green, img_red, img_yellow], axis=-1)\n\n    # Convert to Pseudo-RGB\n    pseudo_rgb_image = create_pseudo_rgb(img_stacked_rgb)\n\n    # Extract features\n    features = extract_cnn_features_pseudo_rgb(pseudo_rgb_image)\n\n    # Store in the list\n    all_features_pseudo_rgb.append(features)\n\n    # Save after every image (for testing)\n    np.save(SAVE_PATH, np.array(all_features_pseudo_rgb))\n    print(f\"Progress saved at {i + start_index + 1}/{num_images} images.\")\n\n# Final save\nnp.save(SAVE_PATH, np.array(all_features_pseudo_rgb))\nprint(\"\\nFeature extraction complete!\")\nprint(\"Final Shape of Extracted Features (Pseudo-RGB):\", np.array(all_features_pseudo_rgb).shape)\nprint(f\"Features saved successfully to '{SAVE_PATH}'\")","metadata":{},"outputs":[],"execution_count":null},{"id":"748c9e04-9aef-4ecb-97f0-2122339f9d8c","cell_type":"markdown","source":"Checking the shape of the vectors again, making sure that they are the correct shape.","metadata":{}},{"id":"fa5206b8-548e-406f-8d28-d69ad307d27d","cell_type":"code","source":"# Load extracted feature files\nfeatures_green = np.load(\"test_green_channel_features.npy\")  # Green-Only features\nfeatures_pseudo_rgb = np.load(\"test_pseudo_rgb_features.npy\")  # Pseudo-RGB features\n\nprint(\"Loaded Green-Only Feature Shape:\", features_green.shape)\nprint(\"Loaded Pseudo-RGB Feature Shape:\", features_pseudo_rgb.shape)","metadata":{},"outputs":[],"execution_count":null},{"id":"46c2f3b2-bfca-4029-89b4-b86a599c04aa","cell_type":"markdown","source":"Checking if there are differences between the absolute features, if there wasn't then there would be no need for different datasets.","metadata":{}},{"id":"498cd6bc-7c46-4ea9-a1d8-2a9ac894c87c","cell_type":"code","source":"diff = np.abs(features_green - features_pseudo_rgb)\nprint(\"Mean Absolute Feature Difference:\", np.mean(diff))\nprint(\"Max Absolute Feature Difference:\", np.max(diff))","metadata":{},"outputs":[],"execution_count":null},{"id":"29059638-30d8-43fd-8638-149c91f752bb","cell_type":"markdown","source":"## Gray Test","metadata":{}},{"id":"f72e4fb5-239d-4d72-b842-5fb0c82bd580","cell_type":"code","source":"# Load pretrained ResNet model and move to GPU\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nmodel = models.resnet50(weights=\"IMAGENET1K_V1\").to(device)\nmodel = torch.nn.Sequential(*(list(model.children())[:-1]))  # Remove last FC layer\nmodel.eval()\n\n# Define transform (normalize and resize for CNN)\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((224, 224)),  # Resize to match CNN input size\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.5], std=[0.5])  # Normalize intensity\n])\n\ndef convert_4ch_to_3ch(img_4ch):\n    \"\"\"Convert 4-channel grayscale image to 3-channel for ResNet.\"\"\"\n    r = img_4ch[:, :, 0]  # Blue (Nucleus)\n    g = img_4ch[:, :, 1]  # Green (Protein of Interest)\n    b = (img_4ch[:, :, 2] + img_4ch[:, :, 3]) / 2  # Merge Red (Microtubules) & Yellow (ER)\n\n    grayscale_3ch = np.stack([r, g, b], axis=-1)  # Shape: (H, W, 3)\n    return np.clip(grayscale_3ch, 0, 255).astype(np.uint8)  # Ensure valid pixel range\n\ndef extract_cnn_features_grayscale(img):\n    \"\"\"Extract features from the 3-channel grayscale image using ResNet.\"\"\"\n    img_tensor = transform(img).unsqueeze(0).to(device)  # Convert to tensor & move to GPU\n    with torch.no_grad():\n        features = model(img_tensor)\n    return features.squeeze().cpu().numpy()  # Move result back to CPU for NumPy\n\n# Directory where images are stored\nTRAIN_PATH = \"D:/hpa-single-cell-image-classification/train\"\nSAVE_PATH = \"test_grayscale_3ch_features.npy\"\n\n# Get first 10 images for testing\nall_files = [f for f in os.listdir(TRAIN_PATH) if \"_green.png\" in f][:10]\nnum_images = len(all_files)\n\nprint(f\"Processing {num_images} Grayscale 3-Channel images...\\n\")\n\n# Initialize or resume from existing file\nif os.path.exists(SAVE_PATH):\n    saved_features = np.load(SAVE_PATH)\n    start_index = len(saved_features)\n    all_features_grayscale = list(saved_features)  # Convert back to list\n    print(f\"Resuming from index {start_index}/{num_images}...\")\nelse:\n    all_features_grayscale = []\n    start_index = 0\n\n# Loop through 10 images with a progress bar\nfor i, filename in enumerate(tqdm(all_files[start_index:], desc=\"Extracting CNN Features (Grayscale 3-CH)\", unit=\"image\")):\n    base_id = filename.replace(\"_green.png\", \"\")\n\n    # Load all four grayscale images (single-channel)\n    img_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_blue.png\"), cv2.IMREAD_GRAYSCALE)\n    img_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_green.png\"), cv2.IMREAD_GRAYSCALE)\n    img_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_red.png\"), cv2.IMREAD_GRAYSCALE)\n    img_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_yellow.png\"), cv2.IMREAD_GRAYSCALE)\n\n    # Stack into (H, W, 4)\n    img_stacked_gray = np.stack([img_blue, img_green, img_red, img_yellow], axis=-1)\n\n    # Convert to 3-channel grayscale\n    img_grayscale_3ch = convert_4ch_to_3ch(img_stacked_gray)\n\n    # Extract features for this image\n    features = extract_cnn_features_grayscale(img_grayscale_3ch)\n\n    # Store feature vector\n    all_features_grayscale.append(features)\n\n    # Save after every image (for testing)\n    np.save(SAVE_PATH, np.array(all_features_grayscale))\n    print(f\"Progress saved at {i + start_index + 1}/{num_images} images.\")\n\n# Final save\nnp.save(SAVE_PATH, np.array(all_features_grayscale))\nprint(\"\\nFeature extraction complete!\")\nprint(\"Final Shape of Extracted Features (Grayscale 3-Channel):\", np.array(all_features_grayscale).shape)\nprint(f\"Features saved successfully to '{SAVE_PATH}'\")","metadata":{},"outputs":[],"execution_count":null},{"id":"8589359c-a4fe-44b6-82b9-ed63b8af9ca0","cell_type":"markdown","source":"Checking the shapes of the vectors to see if they are correct.","metadata":{}},{"id":"4bc6032c-82be-4527-8e3c-8cbd96a3e480","cell_type":"code","source":"# Load extracted feature files\nfeatures_green = np.load(\"test_green_channel_features.npy\")  # Green-Only\nfeatures_pseudo_rgb = np.load(\"test_pseudo_rgb_features.npy\")  # Pseudo-RGB\nfeatures_grayscale = np.load(\"test_grayscale_3ch_features.npy\")  # Grayscale 3-CH\n\nprint(\"Feature Shapes:\")\nprint(\"Green-Only:\", features_green.shape)\nprint(\"Pseudo-RGB:\", features_pseudo_rgb.shape)\nprint(\"Grayscale 3-CH:\", features_grayscale.shape)","metadata":{},"outputs":[],"execution_count":null},{"id":"bbce2169-034f-4637-97b7-0614a8f752fa","cell_type":"markdown","source":"Checking to see if the datasets are different through checking mean features.","metadata":{}},{"id":"d3967ac1-51b6-48a8-8f84-6562269809f7","cell_type":"code","source":"diff_green_gray = np.abs(features_green - features_grayscale)\ndiff_pseudo_gray = np.abs(features_pseudo_rgb - features_grayscale)\n\nprint(\"\\nMean Absolute Feature Difference (Green vs Grayscale):\", np.mean(diff_green_gray))\nprint(\"Max Absolute Feature Difference (Green vs Grayscale):\", np.max(diff_green_gray))\n\nprint(\"\\nMean Absolute Feature Difference (Pseudo-RGB vs Grayscale):\", np.mean(diff_pseudo_gray))\nprint(\"Max Absolute Feature Difference (Pseudo-RGB vs Grayscale):\", np.max(diff_pseudo_gray))","metadata":{},"outputs":[],"execution_count":null},{"id":"53474aef-4f0d-49c9-912a-ed667e06257a","cell_type":"markdown","source":"A simple PCA of these tests, to visualize the datasets.","metadata":{}},{"id":"3de2678c-5809-42bb-877c-6cea7b0e4c22","cell_type":"code","source":"# Reduce to 2D for visualization\npca = PCA(n_components=2)\nfeatures_combined = np.vstack((features_green, features_pseudo_rgb, features_grayscale))  # Stack all feature sets\nfeatures_2d = pca.fit_transform(features_combined)\n\n# Plot Green-Only vs Pseudo-RGB vs Grayscale 3-CH\nplt.figure(figsize=(8,6))\nplt.scatter(features_2d[:10, 0], features_2d[:10, 1], label=\"Green-Only\", alpha=0.7, color=\"blue\")\nplt.scatter(features_2d[10:20, 0], features_2d[10:20, 1], label=\"Pseudo-RGB\", alpha=0.7, color=\"red\")\nplt.scatter(features_2d[20:, 0], features_2d[20:, 1], label=\"Grayscale 3-CH\", alpha=0.7, color=\"gray\")\nplt.xlabel(\"PCA Component 1\")\nplt.ylabel(\"PCA Component 2\")\nplt.title(\"PCA Projection: Green-Only vs Pseudo-RGB vs Grayscale 3-CH Features\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"3b92b5b8-6632-4b86-99ea-c616e224aa83","cell_type":"markdown","source":"## **Next, I will run the program on the entire dataset.**","metadata":{}},{"id":"1aa7d9a1-2f61-4ec8-bcb8-fa41a4700b1e","cell_type":"markdown","source":"---\n# Section 3: Code for Entire Dataset","metadata":{}},{"id":"9272bfde-799e-47d3-89c2-52f0bfdf5783","cell_type":"markdown","source":"## For Green Channel only","metadata":{}},{"id":"f186b2bc-bdc2-4630-9a67-c1de310442a1","cell_type":"code","source":"# Load pretrained ResNet model and move to GPU\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nmodel = models.resnet50(weights=\"IMAGENET1K_V1\").to(device)\nmodel = torch.nn.Sequential(*(list(model.children())[:-1]))  # Remove last FC layer\nmodel.eval()\n\n# Define transform (resize & normalize for CNN)\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((224, 224)),\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.5], std=[0.5])  \n])\n\ndef extract_cnn_features_green(img):\n    \"\"\"Extract features from the Green channel (Protein of Interest) using ResNet.\"\"\"\n    img_green = img[:, :, 1]  # ✅ Extract Green channel (index 1)\n\n    # Convert single-channel grayscale to 3-channel by stacking\n    img_green_3ch = np.stack([img_green, img_green, img_green], axis=-1)  # (H, W, 3)\n\n    img_tensor = transform(img_green_3ch).unsqueeze(0).to(device)  # ✅ Move tensor to GPU\n    with torch.no_grad():\n        features = model(img_tensor)  # ✅ Model is on GPU\n    return features.squeeze().cpu().numpy()  # ✅ Move result back to CPU for NumPy\n\n# Directory where images are stored\nTRAIN_PATH = \"D:/hpa-single-cell-image-classification/train\"\nSAVE_PATH = \"green_channel_features.npy\"\n\n# Get all images\nall_files = [f for f in os.listdir(TRAIN_PATH) if \"_green.png\" in f]\nnum_images = len(all_files)\n\nprint(f\"Processing {num_images} Green channel images...\\n\")\n\n# Initialize or resume from existing file\nif os.path.exists(SAVE_PATH):\n    saved_features = np.load(SAVE_PATH)\n    start_index = len(saved_features)\n    all_features_green = list(saved_features)  # Convert back to list\n    print(f\"Resuming from index {start_index}/{num_images}...\")\nelse:\n    all_features_green = []\n    start_index = 0\n\n# Loop through all images\nfor i, filename in enumerate(tqdm(all_files[start_index:], desc=\"Extracting CNN Features\", unit=\"image\")):\n    base_id = filename.replace(\"_green.png\", \"\")\n\n    # Load all four grayscale images (single-channel)\n    img_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_blue.png\"), cv2.IMREAD_GRAYSCALE)\n    img_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_green.png\"), cv2.IMREAD_GRAYSCALE)\n    img_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_red.png\"), cv2.IMREAD_GRAYSCALE)\n    img_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_yellow.png\"), cv2.IMREAD_GRAYSCALE)\n\n    # Stack into (H, W, 4)\n    img_stacked_rgb = np.stack([img_blue, img_green, img_red, img_yellow], axis=-1)\n\n    # Extract features for this image\n    features = extract_cnn_features_green(img_stacked_rgb)\n\n    # Store in the list\n    all_features_green.append(features)\n\n    # Save every 100 images (prevents data loss)\n    if (i + start_index + 1) % 100 == 0:\n        np.save(SAVE_PATH, np.array(all_features_green))\n        print(f\"Progress saved at {i + start_index + 1}/{num_images} images.\")\n\n# Final save\nnp.save(SAVE_PATH, np.array(all_features_green))\nprint(\"\\nFeature extraction complete!\")\nprint(\"Final Shape of Extracted Features (Green-Only):\", np.array(all_features_green).shape)\nprint(f\"Features saved successfully to '{SAVE_PATH}'\")","metadata":{},"outputs":[],"execution_count":null},{"id":"dffc2660-2fd6-48c0-9c90-150ed57466eb","cell_type":"markdown","source":"## For multicolor","metadata":{}},{"id":"d38c3525-1ca7-42d9-83a3-7427fbe0f083","cell_type":"code","source":"# Load pretrained ResNet model and move to GPU\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nmodel = models.resnet50(weights=\"IMAGENET1K_V1\").to(device)\nmodel = torch.nn.Sequential(*(list(model.children())[:-1]))  # Remove last FC layer\nmodel.eval()\n\n# Define transform (normalize and resize for CNN)\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((224, 224)),  # Resize to match CNN input size\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.5], std=[0.5])  # Normalize intensity\n])\n\ndef create_pseudo_rgb(img):\n    \"\"\"Convert 4-channel image into a 3-channel (R-G-B) image.\"\"\"\n    r = img[:, :, 2]  # Use Red channel (Microtubules)\n    g = img[:, :, 3]  # Use Green channel (Protein of Interest)\n    b = img[:, :, 0]  # Use Blue channel (Nucleus)\n    \n    pseudo_rgb = np.stack([r, g, b], axis=-1)  # Stack into (H, W, 3)\n    return pseudo_rgb.astype(np.uint8)  # Ensure proper dtype\n\ndef extract_cnn_features_pseudo_rgb(img):\n    \"\"\"Extract features from a Pseudo-RGB image using ResNet.\"\"\"\n    img_tensor = transform(img).unsqueeze(0).to(device)  # ✅ Move tensor to GPU\n    with torch.no_grad():\n        features = model(img_tensor)  # ✅ Model is on GPU\n    return features.squeeze().cpu().numpy()  # ✅ Move result back to CPU for NumPy\n\n# Directory where images are stored\nTRAIN_PATH = \"D:/hpa-single-cell-image-classification/train\"\nSAVE_PATH = \"pseudo_rgb_features.npy\"\n\n# Get all images\nall_files = [f for f in os.listdir(TRAIN_PATH) if \"_green.png\" in f]\nnum_images = len(all_files)\n\nprint(f\"Processing {num_images} Pseudo-RGB images...\\n\")\n\n# Initialize or resume from existing file\nif os.path.exists(SAVE_PATH):\n    saved_features = np.load(SAVE_PATH)\n    start_index = len(saved_features)\n    all_features_pseudo_rgb = list(saved_features)  # Convert back to list\n    print(f\"Resuming from index {start_index}/{num_images}...\")\nelse:\n    all_features_pseudo_rgb = []\n    start_index = 0\n\n# Loop through all images\nfor i, filename in enumerate(tqdm(all_files[start_index:], desc=\"Extracting CNN Features (Pseudo-RGB)\", unit=\"image\")):\n    base_id = filename.replace(\"_green.png\", \"\")\n\n    # Load all four images as grayscale (single-channel)\n    img_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_blue.png\"), cv2.IMREAD_GRAYSCALE)\n    img_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_green.png\"), cv2.IMREAD_GRAYSCALE)\n    img_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_red.png\"), cv2.IMREAD_GRAYSCALE)\n    img_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_yellow.png\"), cv2.IMREAD_GRAYSCALE)\n\n    # Stack into (H, W, 4)\n    img_stacked_rgb = np.stack([img_blue, img_green, img_red, img_yellow], axis=-1)\n\n    # Convert to Pseudo-RGB\n    pseudo_rgb_image = create_pseudo_rgb(img_stacked_rgb)\n\n    # Extract features\n    features = extract_cnn_features_pseudo_rgb(pseudo_rgb_image)\n\n    # Store in the list\n    all_features_pseudo_rgb.append(features)\n\n    # Save every 100 images (to prevent data loss)\n    if (i + start_index + 1) % 100 == 0:\n        np.save(SAVE_PATH, np.array(all_features_pseudo_rgb))\n        print(f\"Progress saved at {i + start_index + 1}/{num_images} images.\")\n\n# Final save\nnp.save(SAVE_PATH, np.array(all_features_pseudo_rgb))\nprint(\"\\nFeature extraction complete!\")\nprint(\"Final Shape of Extracted Features (Pseudo-RGB):\", np.array(all_features_pseudo_rgb).shape)\nprint(f\"Features saved successfully to '{SAVE_PATH}'\")","metadata":{},"outputs":[],"execution_count":null},{"id":"b58b29f9-a02e-48f8-9161-2c6b835d20ff","cell_type":"markdown","source":"## For Grayscale made on 1B","metadata":{}},{"id":"932be47c-8355-4d99-b015-e2e2c10ed003","cell_type":"markdown","source":"This averages together Red and Yellow, mictotubules and ER together into one single layer. Then Resnet, which takes 3 layers, can take this four layer data. ","metadata":{}},{"id":"7dde321e-6e05-42b5-bac1-0b590e25545b","cell_type":"code","source":"# Load pretrained ResNet model and move to GPU\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nmodel = models.resnet50(weights=\"IMAGENET1K_V1\").to(device)\nmodel = torch.nn.Sequential(*(list(model.children())[:-1]))  # Remove last FC layer\nmodel.eval()\n\n# Define transform (normalize and resize for CNN)\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((224, 224)),  # Resize to match CNN input size\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.5], std=[0.5])  # Normalize intensity\n])\n\ndef convert_4ch_to_3ch(img_4ch):\n    \"\"\"Convert 4-channel grayscale image to 3-channel for ResNet.\"\"\"\n    r = img_4ch[:, :, 0]  # Blue (Nucleus)\n    g = img_4ch[:, :, 1]  # Green (Protein of Interest)\n    b = (img_4ch[:, :, 2] + img_4ch[:, :, 3]) / 2  # Merge Red (Microtubules) & Yellow (ER)\n\n    grayscale_3ch = np.stack([r, g, b], axis=-1)  # Shape: (H, W, 3)\n    return np.clip(grayscale_3ch, 0, 255).astype(np.uint8)  # Ensure valid pixel range\n\ndef extract_cnn_features_grayscale(img):\n    \"\"\"Extract features from the 3-channel grayscale image using ResNet.\"\"\"\n    img_tensor = transform(img).unsqueeze(0).to(device)  # Convert to tensor & move to GPU\n    with torch.no_grad():\n        features = model(img_tensor)\n    return features.squeeze().cpu().numpy()  # Move result back to CPU for NumPy\n\n# Directory where images are stored\nTRAIN_PATH = \"D:/hpa-single-cell-image-classification/train\"\nSAVE_PATH = \"grayscale_3ch_features.npy\"\n\n# Get all images\nall_files = [f for f in os.listdir(TRAIN_PATH) if \"_green.png\" in f]\nnum_images = len(all_files)\n\nprint(f\"Processing {num_images} Grayscale 3-Channel images...\\n\")\n\n# Initialize or resume from existing file\nif os.path.exists(SAVE_PATH):\n    saved_features = np.load(SAVE_PATH)\n    start_index = len(saved_features)\n    all_features_grayscale = list(saved_features)  # Convert back to list\n    print(f\"Resuming from index {start_index}/{num_images}...\")\nelse:\n    all_features_grayscale = []\n    start_index = 0\n\n# Loop through all images\nfor i, filename in enumerate(tqdm(all_files[start_index:], desc=\"Extracting CNN Features (Grayscale 3-CH)\", unit=\"image\")):\n    base_id = filename.replace(\"_green.png\", \"\")\n\n    # Load all four grayscale images (single-channel)\n    img_blue = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_blue.png\"), cv2.IMREAD_GRAYSCALE)\n    img_green = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_green.png\"), cv2.IMREAD_GRAYSCALE)\n    img_red = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_red.png\"), cv2.IMREAD_GRAYSCALE)\n    img_yellow = cv2.imread(os.path.join(TRAIN_PATH, f\"{base_id}_yellow.png\"), cv2.IMREAD_GRAYSCALE)\n\n    # Stack into (H, W, 4)\n    img_stacked_gray = np.stack([img_blue, img_green, img_red, img_yellow], axis=-1)\n\n    # Convert to 3-channel grayscale\n    img_grayscale_3ch = convert_4ch_to_3ch(img_stacked_gray)\n\n    # Extract features\n    features = extract_cnn_features_grayscale(img_grayscale_3ch)\n\n    # Store feature vector\n    all_features_grayscale.append(features)\n\n    # Save every 100 images (to prevent data loss)\n    if (i + start_index + 1) % 100 == 0:\n        np.save(SAVE_PATH, np.array(all_features_grayscale))\n        print(f\"Progress saved at {i + start_index + 1}/{num_images} images.\")\n\n# Final save\nnp.save(SAVE_PATH, np.array(all_features_grayscale))\nprint(\"\\nFeature extraction complete!\")\nprint(\"Final Shape of Extracted Features (Grayscale 3-Channel):\", np.array(all_features_grayscale).shape)\nprint(f\"Features saved successfully to '{SAVE_PATH}'\")","metadata":{},"outputs":[],"execution_count":null},{"id":"dd77ad5c-9a54-4a24-8de3-1dcb0daaf025","cell_type":"markdown","source":"---\n## PCA on the entire dataset","metadata":{}},{"id":"de1f39a6-7ef5-4044-85c3-ad8827f70ac4","cell_type":"code","source":"features_green = np.load(\"green_channel_features.npy\")  # Entire dataset features\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n\n# Reduce to 2D for visualization\npca = PCA(n_components=2)\nfeatures_2d = pca.fit_transform(features_combined)\nplt.figure(figsize=(8,6))\nnum_images = features_green.shape[0]\n\nplt.scatter(features_2d[:num_images, 0], features_2d[:num_images, 1], label=\"Green-Only\", alpha=0.7, color=\"blue\")\nplt.scatter(features_2d[num_images:2*num_images, 0], features_2d[num_images:2*num_images, 1], label=\"Pseudo-RGB\", alpha=0.7, color=\"red\")\nplt.scatter(features_2d[2*num_images:, 0], features_2d[2*num_images:, 1], label=\"Grayscale 3-CH\", alpha=0.7, color=\"gray\")\n\nplt.xlabel(\"PCA Component 1\")\nplt.ylabel(\"PCA Component 2\")\nplt.title(\"PCA Projection: Green-Only vs Pseudo-RGB vs Grayscale 3-CH Features\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"4201330a-9182-4f5c-a43d-7c17ecdcf249","cell_type":"markdown","source":"As can be seen some components are massive and some are not. The graph only shows the green due to this mismatch. In order to fix this, I will standardize the features.","metadata":{}},{"id":"541238ca-d8da-4f0a-8a6b-4af5b72236fe","cell_type":"code","source":"# Load all feature sets\nfeatures_green = np.load(\"green_channel_features.npy\")  \nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n# Stack all feature sets into one array\nfeatures_combined = np.vstack((features_green, features_pseudo_rgb, features_grayscale))\n\n# Standardize the features to have mean=0 and std=1\nscaler = StandardScaler()\nfeatures_scaled = scaler.fit_transform(features_combined)\n\n# Perform PCA reduction\npca = PCA(n_components=2)\nfeatures_2d = pca.fit_transform(features_scaled)\n\n# Define number of images in each feature set\nnum_images = features_green.shape[0]\n\n# Plot PCA projections\nplt.figure(figsize=(8,6))\nplt.scatter(features_2d[:num_images, 0], features_2d[:num_images, 1], label=\"Green-Only\", alpha=0.7, color=\"blue\")\nplt.scatter(features_2d[num_images:2*num_images, 0], features_2d[num_images:2*num_images, 1], label=\"Pseudo-RGB\", alpha=0.7, color=\"red\")\nplt.scatter(features_2d[2*num_images:, 0], features_2d[2*num_images:, 1], label=\"Grayscale 3-CH\", alpha=0.7, color=\"gray\")\n\nplt.xlabel(\"PCA Component 1\")\nplt.ylabel(\"PCA Component 2\")\nplt.title(\"PCA Projection: Green-Only vs Pseudo-RGB vs Grayscale 3-CH Features\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"d408af4f-1c9d-4855-9df6-4818e28a4c79","cell_type":"markdown","source":"Much better! The PCA now shows all the datasets. 3 Different populations have been created, each with different charecteristics. Now I will try to cluster each population. ","metadata":{}},{"id":"b7f6b65a-2b71-451c-a4e7-79bf79542fab","cell_type":"markdown","source":"---\n# Section 4: DPGMM TESTS","metadata":{}},{"id":"1160d62e-2301-4393-97cd-837c0600bfd7","cell_type":"markdown","source":"DPGMM is a Dirichlet Process GMM. It is a soft clustering algorithm, which means that a datapoint can belong to multiple clusters. GMM fits Gaussians to the datasets as shown in the Stanford CS229 Notes. The following tests, which do not appear in the notes, will show whether the data is Gaussian or non-Gaussian. If it is non-Gaussian the GMM will fail. If it is Gaussian the GMM will work.\n\n1. Shapiro-Wilk p-value\n\n    - The Shapiro-Wilk test checks normality.\n    - Null hypothesis **H0**: The data is normally distributed.\n    - Alternative hypothesis **Ha**: The data is not normally distributed.\n    - If p-value < 0.05, we reject **H0**, meaning the data is not normal.\n\n2. Skewness\n\n    - Gaussian distributions have skewness = 0.\n    - Skewness > 0 → Right-skewed (long right tail).\n    - Skewness < 0 → Left-skewed (long left tail).\n    - If skewness is between -0.5 and 0.5, it's approximately symmetric.\n\n3. Kurtosis\n\n    - Gaussian distributions have kurtosis = 0 (excess kurtosis).\n    - Kurtosis > 0 → Heavy tails (leptokurtic, outlier-prone).\n    - Kurtosis < 0 → Light tails (platykurtic, fewer outliers).\n    - Kurtosis between -1 and 1 is close to normal.","metadata":{}},{"id":"63a8caba-2337-4cf1-8474-612488c7908e","cell_type":"code","source":"features_grayscale = np.load(\"grayscale_3ch_features.npy\") # ✅ Load features \n\n# ✅ Sample a subset (to speed up computation) \nsubset_size = 5000 \nfeatures_grayscale = features_grayscale[:subset_size] \n\n\n# ✅ Pick 5 random features from each dataset \nnp.random.seed(42)\nrandom_features_gray = np.random.choice(features_grayscale.shape[1], size=5, replace=False)  \n\n# ✅ Analyze each feature in Grayscale \nprint(\"\\n📊 **Grayscale Features**\") \nfor feature_idx in random_features_gray:     \n    stat, p_value = shapiro(features_grayscale[:, feature_idx])     \n    skewness = skew(features_grayscale[:, feature_idx])     \n    kurt = kurtosis(features_grayscale[:, feature_idx])      \n    print(f\"🔍 Feature {feature_idx}:\")     \n    print(f\"   - Shapiro-Wilk p-value = {p_value}\")    \n    print(f\"   - Skewness = {skewness}\")    \n    print(f\"   - Kurtosis = {kurt}\\n\") ","metadata":{},"outputs":[],"execution_count":null},{"id":"3df3f103-6c53-4609-8204-e0d3a9511ba7","cell_type":"markdown","source":"## Results: \n\n- Feature 1476: Not Gaussian due to very small p-value.\n- Feature 693: Not Gaussian due to skew and heavy tails.\n- Feature 100: Not Gaussian due to skew and heavy tails.\n- Fature 1421: Not Gaussian due to skew and heavy tails.\n- Feature 1810: Not Gaussian due to skew and heavy tails.","metadata":{}},{"id":"19f95a0a-7c21-49d8-b94d-982e92bcbbd2","cell_type":"code","source":"features_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\") \n \n\n# ✅ Sample a subset (to speed up computation) \nsubset_size = 5000 \nfeatures_pseudo_rgb = features_pseudo_rgb[:subset_size] \nfeatures_grayscale = features_grayscale[:subset_size]  \n\n# ✅ Pick 5 random features from each dataset \nnp.random.seed(42) \nrandom_features_pseudo = np.random.choice(features_pseudo_rgb.shape[1], size=5, replace=False) \nrandom_features_gray = np.random.choice(features_grayscale.shape[1], size=5, replace=False)  \n# ✅ Analyze each feature in Pseudo-RGB \nprint(\"\\n📊 **Pseudo-RGB Features**\") \nfor feature_idx in random_features_pseudo:     \n    stat, p_value = shapiro(features_pseudo_rgb[:, feature_idx])     \n    skewness = skew(features_pseudo_rgb[:, feature_idx])     \n    kurt = kurtosis(features_pseudo_rgb[:, feature_idx])      \n    print(f\"🔍 Feature {feature_idx}:\")     \n    print(f\"   - Shapiro-Wilk p-value = {p_value}\")     \n    print(f\"   - Skewness = {skewness}\")     \n    print(f\"   - Kurtosis = {kurt}\\n\")  ","metadata":{},"outputs":[],"execution_count":null},{"id":"dacfe110-567a-4adb-8e09-0e1cad09b4a5","cell_type":"markdown","source":"## Results:\n\n- None of these features are Gaussian because their p-values are extremely low.\n- Features 693, 100, and 1421 are heavily skewed and have extreme kurtosis, indicating high asymmetry and heavy tails (outliers).\n- Feature 1476 is the closest to Gaussian with moderate skew and kurtosis, but its p-value still rejects normality.\n- Feature 1810 has moderate skew/kurtosis but is still not Gaussian.","metadata":{}},{"id":"240908fd-ca8c-47a6-94dc-bb07d76866d4","cell_type":"code","source":"# ✅ Load features\nfeatures_green = np.load(\"green_channel_features.npy\")\n\n# ✅ Sample a subset\nsubset_size = 5000\nfeatures_green = features_green[:subset_size]\n\n# ✅ Pick 5 random features\nnp.random.seed(42)\nrandom_features = np.random.choice(features_green.shape[1], size=5, replace=False)\n\n# ✅ Analyze each feature\nfor feature_idx in random_features:\n    stat, p_value = shapiro(features_green[:, feature_idx])\n    skewness = skew(features_green[:, feature_idx])\n    kurt = kurtosis(features_green[:, feature_idx])\n\n    print(f\"\\n🔍 Feature {feature_idx}:\")\n    print(f\"   - Shapiro-Wilk p-value = {p_value}\")\n    print(f\"   - Skewness = {skewness}\")\n    print(f\"   - Kurtosis = {kurt}\")","metadata":{},"outputs":[],"execution_count":null},{"id":"5a57b19b-091d-4346-a365-dec0b3af54fd","cell_type":"markdown","source":"## Results: \n- None of these features are Gaussian, as all have very low p-values.\n- Features 693 and 1810 are the most non-Gaussian, with high skewness and extreme kurtosis (long-tailed distributions).\n- Features 1476 and 1421 are the closest to normal but still fail normality tests due to moderate skew and heavy tails.","metadata":{}},{"id":"3277bcf3-0b28-460b-8bdd-fdd09a1533fb","cell_type":"markdown","source":"## Running the DPGMM\n\nGMM algorithms use ELBO to work. A full derivation is available in the Stanford CS229 Notes, Chapter 11.3\n\n$$\n\\sum_{i} \\mathbb{E}_{Q_i(z^{(i)})} \\left[ \\log p(x^{(i)} \\mid z^{(i)}; \\theta) \\right]\n$$\nIs the formula for ELBO as per 11.3 in the Stanford CS229 Notes. \n\n$$\n\\sum_{i} \\sum_{z^{(i)}} Q_i(z^{(i)}) \\log p(x^{(i)} \\mid z^{(i)}; \\theta)\n$$\nIs the ELBO expanded.\n\nA further expansion gives:\n\n$$\n\\sum_{i} \\sum_{z^{(i)}} Q_i(z^{(i)}) \\log \\frac{p(x^{(i)}, z^{(i)}; \\theta)}{Q_i(z^{(i)})}\n$$\n\nHere we see the variance term in the numerator. As we shall see, the variance will explode as the DPGMM attempts to fit Gaussian distributions to non-Gaussian distributions. We should expect massive ELBO values.\n\nA common failure mode of DPGMM is a massive ELBO which can be caused by non-Gaussian data, as the following tests will show.","metadata":{}},{"id":"d2ffc1c1-0ee4-43f8-b330-454fdc0e17ee","cell_type":"code","source":"# ✅ Load features\nfeatures_green = np.load(\"green_channel_features.npy\")\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n# ✅ Check for NaNs\nprint(\"Checking for NaNs...\")\nfeatures_green = np.nan_to_num(features_green)\nfeatures_pseudo_rgb = np.nan_to_num(features_pseudo_rgb)\nfeatures_grayscale = np.nan_to_num(features_grayscale)\n\n# ✅ Sample a smaller subset of data\nsubset_size = 5000\nfeatures_green = features_green[:subset_size]\nfeatures_pseudo_rgb = features_pseudo_rgb[:subset_size]\nfeatures_grayscale = features_grayscale[:subset_size]\n\n# ✅ Standardize features\nscaler = StandardScaler()\nfeatures_green = scaler.fit_transform(features_green)\nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb)\nfeatures_grayscale = scaler.transform(features_grayscale)\n\n# ✅ DP-GMM Model Setup (Fixes Applied)\ndpgmm = BayesianGaussianMixture(\n    n_components=50,  \n    weight_concentration_prior=1e-2,  \n    weight_concentration_prior_type=\"dirichlet_process\",\n    covariance_type=\"full\",\n    reg_covar=1e-4,  \n    init_params=\"kmeans\",  \n    n_init=1,  \n    max_iter=5,  # ⬆ More iterations per fit\n    warm_start=True,  \n    tol=1e-4,  # ⬇ Reduce tolerance\n    verbose=2,\n    random_state=42\n)\n\n# ✅ Track ELBO per iteration\nnum_iterations = 5\nelbo_values = []\n\n# ✅ Check if `warm_start` is actually working\nprint(\"Warm Start Enabled:\", dpgmm.warm_start)\n\nprevious_means = None  # Track mean changes\n\nprint(\"\\nFitting DP-GMM on Green-Only Features with Manual Iteration Tracking...\")\nfor i in range(num_iterations):\n    dpgmm.lower_bound_ = -np.inf  # 🔄 Reset ELBO to force recomputation\n    \n    dpgmm.fit(features_green)  # Run one step of EM\n    elbo_values.append(dpgmm.lower_bound_)\n    \n    # ✅ Check if means are changing\n    if previous_means is not None:\n        print(f\"Iteration {i+1}, Mean Shift: {np.linalg.norm(dpgmm.means_ - previous_means)}\")\n    previous_means = np.copy(dpgmm.means_)\n    \n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\n# ✅ Save ELBO values\nnp.save(\"dpgmm_elbo_green_debug.npy\", elbo_values)\n\n# 📊 Plot ELBO over iterations\nplt.figure(figsize=(8,5))\nplt.plot(range(1, num_iterations + 1), elbo_values, marker='o', linestyle='-')\nplt.xlabel(\"Iteration\")\nplt.ylabel(\"ELBO\")\nplt.title(\"ELBO Over Iterations for Green-Only Features (Debugging)\")\nplt.grid()\nplt.show()\n\n# ✅ Get cluster assignments\nclusters_green = dpgmm.predict(features_green)\nnp.save(\"dpgmm_clusters_greenscidebug.npy\", clusters_green)","metadata":{},"outputs":[],"execution_count":null},{"id":"ebef7027-11a3-4f9c-a399-b38364b0c36c","cell_type":"markdown","source":"ELBO as predicted, is massive. ","metadata":{}},{"id":"9a899c89-d94f-414e-b1a8-8ba14fcac354","cell_type":"code","source":"# Load all feature sets\nfeatures_green = np.load(\"green_channel_features.npy\")      # (21806, 2048)\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")    # (21806, 2048)\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")  # (21806, 2048)\n\n# DP-GMM Model Setup\ndpgmm = BayesianGaussianMixture(\n    n_components=100,  # Max clusters, DP picks the optimal number\n    covariance_type=\"full\",\n    weight_concentration_prior_type=\"dirichlet_process\",\n    verbose = 2,\n    max_iter=5,\n    random_state=42\n)\n\n# Run DP-GMM on Green Features\nprint(\"Fitting DP-GMM on Green-Only Features...\")\ndpgmm.fit(features_green)\nclusters_green = dpgmm.predict(features_green)\nnp.save(\"dpgmm_clusters_greenscimk1.npy\", clusters_green)\n\n# Run DP-GMM on Pseudo-RGB Features\nprint(\"Fitting DP-GMM on Pseudo-RGB Features...\")\ndpgmm.fit(features_pseudo_rgb)\nclusters_pseudo_rgb = dpgmm.predict(features_pseudo_rgb)\nnp.save(\"dpgmm_clusters_pseudo_rgbscimk1.npy\", clusters_pseudo_rgb)\n\n# Run DP-GMM on Grayscale Features\nprint(\"Fitting DP-GMM on Grayscale Features...\")\ndpgmm.fit(features_grayscale)\nclusters_grayscale = dpgmm.predict(features_grayscale)\nnp.save(\"dpgmm_clusters_grayscalescimk1.npy\", clusters_grayscale)\n\nprint(\"\\nClustering complete! Clusters found:\")\nprint(f\"- Green-Only: {len(set(clusters_green))} clusters\")\nprint(f\"- Pseudo-RGB: {len(set(clusters_pseudo_rgb))} clusters\")\nprint(f\"- Grayscale: {len(set(clusters_grayscale))} clusters\")\n\n# 📊 Plot Cluster Distributions\nplt.figure(figsize=(10,5))\nplt.hist(clusters_green, bins=50, alpha=0.6, label=\"Green-Only\", color=\"green\")\nplt.hist(clusters_pseudo_rgb, bins=50, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\")\nplt.hist(clusters_grayscale, bins=50, alpha=0.6, label=\"Grayscale 3-CH\", color=\"gray\")\nplt.xlabel(\"Cluster ID\")\nplt.ylabel(\"Number of Samples\")\nplt.title(\"Cluster Distributions Across Feature Representations\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"d644d433-74f0-4475-8ece-4eada1d99492","cell_type":"markdown","source":"Once again, ELBO is massive. ","metadata":{}},{"id":"ee3d35c4-b03a-4240-a9c5-f5db431863f9","cell_type":"markdown","source":"## K-Means initialization, and sparsity, and preventing singularity and multiple runs\n\nAfter altering the algorithm a bit, we shall see that the problem does not go away since it is based on the inherent non-Gaussanity of the data. Massive ELBO values are still expected. ","metadata":{}},{"id":"d3e79622-8ef0-4f4e-831e-a084c80398e6","cell_type":"code","source":"# ✅ Load features\nfeatures_green = np.load(\"green_channel_features.npy\")\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n# ✅ Check for NaNs\nprint(\"Checking for NaNs...\")\nfeatures_green = np.nan_to_num(features_green)\nfeatures_pseudo_rgb = np.nan_to_num(features_pseudo_rgb)\nfeatures_grayscale = np.nan_to_num(features_grayscale)\n\n# ✅ Standardize features\nscaler = StandardScaler()\nfeatures_green = scaler.fit_transform(features_green)\nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb)\nfeatures_grayscale = scaler.transform(features_grayscale)\n\n# ✅ DP-GMM Model Setup\ndpgmm = BayesianGaussianMixture(\n    n_components=100,  # Maximum possible clusters\n    weight_concentration_prior=1e-2,  # Controls sparsity\n    weight_concentration_prior_type=\"dirichlet_process\",  # ✅ Enables DP behavior\n    covariance_type=\"full\",\n    reg_covar=1e-4,  # ✅ Prevent covariance collapse\n    init_params=\"kmeans\",  # ✅ Better initialization\n    n_init=1,  # Single initialization for tracking\n    max_iter=5,  # Run for multiple steps per update\n    warm_start=True,  # ✅ Keep improving model instead of restarting\n    tol=1e-3,  # ✅ Loosen tolerance for better convergence\n    verbose=2,\n    random_state=42\n)\n# ✅ Track ELBO per iteration\nelbo_values = []\nnum_iterations = 5  # Set how many iterations you want to track manually\n\nprint(\"\\nFitting DP-GMM on Green-Only Features with Manual Iteration Tracking...\")\nfor i in range(num_iterations):\n    dpgmm.fit(features_green)\n    elbo_values.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\n# ✅ Save ELBO values\nnp.save(\"dpgmm_elbo_greenmk2.npy\", elbo_values)\n\n# 📊 Plot ELBO over iterations\nplt.figure(figsize=(8,5))\nplt.plot(range(1, num_iterations + 1), elbo_values, marker='o', linestyle='-')\nplt.xlabel(\"Iteration\")\nplt.ylabel(\"ELBO\")\nplt.title(\"ELBO Over Iterations for Green-Only Features\")\nplt.grid()\nplt.show()\n\n# ✅ Get cluster assignments\nclusters_green = dpgmm.predict(features_green)\nnp.save(\"dpgmm_clusters_greenscimk2.npy\", clusters_green)\n\n# ✅ Repeat for Pseudo-RGB Features\nprint(\"\\nFitting DP-GMM on Pseudo-RGB Features...\")\nelbo_values_pseudo = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_pseudo_rgb)\n    elbo_values_pseudo.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_pseudomk2.npy\", elbo_values_pseudo)\nclusters_pseudo_rgb = dpgmm.predict(features_pseudo_rgb)\nnp.save(\"dpgmm_clusters_pseudo_rgbscimk2.npy\", clusters_pseudo_rgb)\n\n# ✅ Repeat for Grayscale Features\nprint(\"\\nFitting DP-GMM on Grayscale Features...\")\nelbo_values_gray = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_grayscale)\n    elbo_values_gray.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_graymk2.npy\", elbo_values_gray)\nclusters_grayscale = dpgmm.predict(features_grayscale)\nnp.save(\"dpgmm_clusters_grayscalescimk2.npy\", clusters_grayscale)\n\n# 📊 Compare Cluster Distributions\nplt.figure(figsize=(10,5))\nplt.hist(clusters_green, bins=50, alpha=0.6, label=\"Green-Only\", color=\"green\")\nplt.hist(clusters_pseudo_rgb, bins=50, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\")\nplt.hist(clusters_grayscale, bins=50, alpha=0.6, label=\"Grayscale 3-CH\", color=\"gray\")\nplt.xlabel(\"Cluster ID\")\nplt.ylabel(\"Number of Samples\")\nplt.title(\"Cluster Distributions Across Feature Representations\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"093128ad-b374-4614-b135-1235afd080cf","cell_type":"markdown","source":"ELBO is massive as expected since the data is non-Gaussian. ","metadata":{}},{"id":"709bb52c-5f7c-4228-86df-633edec5a76b","cell_type":"code","source":"# ✅ Load features\nfeatures_green = np.load(\"green_channel_features.npy\")\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n# ✅ Check for NaNs\nprint(\"Checking for NaNs...\")\nfeatures_green = np.nan_to_num(features_green)\nfeatures_pseudo_rgb = np.nan_to_num(features_pseudo_rgb)\nfeatures_grayscale = np.nan_to_num(features_grayscale)\n\n# ✅ Standardize features\nscaler = StandardScaler()\nfeatures_green = scaler.fit_transform(features_green)\nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb)\nfeatures_grayscale = scaler.transform(features_grayscale)\n\n# ✅ DP-GMM Model Setup\ndpgmm = BayesianGaussianMixture(\n    n_components=100,  # Maximum possible clusters\n    weight_concentration_prior=1e-2,  # Controls sparsity\n    weight_concentration_prior_type=\"dirichlet_process\",  # ✅ Enables DP behavior\n    covariance_type=\"full\",\n    reg_covar=1e-5,  # Prevents singularities\n    init_params='kmeans',  # Better initialization\n    n_init=1,  # Only one initialization since we're tracking iterations\n    max_iter=1,  # Run only one iteration per step\n    warm_start=True,  # ✅ Keeps updating the model instead of restarting\n    verbose=2,\n    random_state=42\n)\n\n# ✅ Track ELBO per iteration\nelbo_values = []\nnum_iterations = 10  # Set how many iterations you want to track manually\n\nprint(\"\\nFitting DP-GMM on Green-Only Features with Manual Iteration Tracking...\")\nfor i in range(num_iterations):\n    dpgmm.fit(features_green)\n    elbo_values.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\n# ✅ Save ELBO values\nnp.save(\"dpgmm_elbo_greenmk2r.npy\", elbo_values)\n\n# 📊 Plot ELBO over iterations\nplt.figure(figsize=(8,5))\nplt.plot(range(1, num_iterations + 1), elbo_values, marker='o', linestyle='-')\nplt.xlabel(\"Iteration\")\nplt.ylabel(\"ELBO\")\nplt.title(\"ELBO Over Iterations for Green-Only Features\")\nplt.grid()\nplt.show()\n\n# ✅ Get cluster assignments\nclusters_green = dpgmm.predict(features_green)\nnp.save(\"dpgmm_clusters_greenscimk2r.npy\", clusters_green)\n\n# ✅ Repeat for Pseudo-RGB Features\nprint(\"\\nFitting DP-GMM on Pseudo-RGB Features...\")\nelbo_values_pseudo = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_pseudo_rgb)\n    elbo_values_pseudo.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_pseudomk2r.npy\", elbo_values_pseudo)\nclusters_pseudo_rgb = dpgmm.predict(features_pseudo_rgb)\nnp.save(\"dpgmm_clusters_pseudo_rgbscimk2r.npy\", clusters_pseudo_rgb)\n\n# ✅ Repeat for Grayscale Features\nprint(\"\\nFitting DP-GMM on Grayscale Features...\")\nelbo_values_gray = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_grayscale)\n    elbo_values_gray.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_graymk2r.npy\", elbo_values_gray)\nclusters_grayscale = dpgmm.predict(features_grayscale)\nnp.save(\"dpgmm_clusters_grayscalescimk2r.npy\", clusters_grayscale)\n\n# 📊 Compare Cluster Distributions\nplt.figure(figsize=(10,5))\nplt.hist(clusters_green, bins=50, alpha=0.6, label=\"Green-Only\", color=\"green\")\nplt.hist(clusters_pseudo_rgb, bins=50, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\")\nplt.hist(clusters_grayscale, bins=50, alpha=0.6, label=\"Grayscale 3-CH\", color=\"gray\")\nplt.xlabel(\"Cluster ID\")\nplt.ylabel(\"Number of Samples\")\nplt.title(\"Cluster Distributions Across Feature Representations\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"5a081070-072b-4d7a-a79c-5f51ecba50fa","cell_type":"markdown","source":"ELBO is still massive due to non-Gaussian data.\n\n## Results:\n\nELBO is massive due to the non-Gaussian data. Whenever GMMs create massive ELBO values the underlying data should be tested for non-Gaussianity. The tests should be carried out before the GMM is run to see if GMM will cluster it anyways. ","metadata":{}},{"id":"98f2d4ee-2113-4c34-b7be-37f75bb7c7c9","cell_type":"markdown","source":"## No K-Means Initialization","metadata":{}},{"id":"0ad80a83-8e35-4eb7-b87d-ddd94f3b6191","cell_type":"markdown","source":"Without K-Means initialization the problem still does not disappear due to non-Gaussianity of the data. ","metadata":{}},{"id":"123f4deb-3ddf-4831-967a-cc81aca57171","cell_type":"code","source":"# ✅ Load features\nfeatures_green = np.load(\"green_channel_features.npy\")\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n# ✅ Check for NaNs\nprint(\"Checking for NaNs...\")\nfeatures_green = np.nan_to_num(features_green)\nfeatures_pseudo_rgb = np.nan_to_num(features_pseudo_rgb)\nfeatures_grayscale = np.nan_to_num(features_grayscale)\n\n# ✅ Standardize features\nscaler = StandardScaler()\nfeatures_green = scaler.fit_transform(features_green)\nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb)\nfeatures_grayscale = scaler.transform(features_grayscale)\n\n# ✅ DP-GMM Model Setup (No Restart Between Iterations)\ndpgmm = BayesianGaussianMixture(\n    n_components=100,  # Maximum possible clusters\n    weight_concentration_prior=1e-2,  # Controls sparsity\n    weight_concentration_prior_type=\"dirichlet_process\",  # ✅ Enables DP behavior\n    covariance_type=\"full\",\n    reg_covar=1e-5,  # Prevents singularities\n    n_init=1,  # Single initialization for tracking\n    max_iter=1,  # Run one iteration at a time\n   # warm_start=True,  # ✅ Keeps improving model instead of resetting\n    verbose=2,\n    random_state=42\n)\n\n# ✅ Track ELBO per iteration\nnum_iterations = 10  # Number of manual iterations\nelbo_values = []\n\nprint(\"\\nFitting DP-GMM on Green-Only Features with Manual Iteration Tracking...\")\nfor i in range(num_iterations):\n    dpgmm.fit(features_green)\n    elbo_values.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\n# ✅ Save ELBO values\nnp.save(\"dpgmm_elbo_green_mk3.npy\", elbo_values)\n\n# 📊 Plot ELBO over iterations\nplt.figure(figsize=(8,5))\nplt.plot(range(1, num_iterations + 1), elbo_values, marker='o', linestyle='-')\nplt.xlabel(\"Iteration\")\nplt.ylabel(\"ELBO\")\nplt.title(\"ELBO Over Iterations for Green-Only Features\")\nplt.grid()\nplt.show()\n\n# ✅ Get cluster assignments\nclusters_green = dpgmm.predict(features_green)\nnp.save(\"dpgmm_clusters_greenscimk3.npy\", clusters_green)\n\n# ✅ Repeat for Pseudo-RGB Features\nprint(\"\\nFitting DP-GMM on Pseudo-RGB Features...\")\nelbo_values_pseudo = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_pseudo_rgb)\n    elbo_values_pseudo.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_pseudo_mk3.npy\", elbo_values_pseudo)\nclusters_pseudo_rgb = dpgmm.predict(features_pseudo_rgb)\nnp.save(\"dpgmm_clusters_pseudo_rgbscimk3.npy\", clusters_pseudo_rgb)\n\n# ✅ Repeat for Grayscale Features\nprint(\"\\nFitting DP-GMM on Grayscale Features...\")\nelbo_values_gray = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_grayscale)\n    elbo_values_gray.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_gray_mk3.npy\", elbo_values_gray)\nclusters_grayscale = dpgmm.predict(features_grayscale)\nnp.save(\"dpgmm_clusters_grayscalescimk3.npy\", clusters_grayscale)\n\n# 📊 Compare Cluster Distributions\nplt.figure(figsize=(10,5))\nplt.hist(clusters_green, bins=50, alpha=0.6, label=\"Green-Only\", color=\"green\")\nplt.hist(clusters_pseudo_rgb, bins=50, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\")\nplt.hist(clusters_grayscale, bins=50, alpha=0.6, label=\"Grayscale 3-CH\", color=\"gray\")\nplt.xlabel(\"Cluster ID\")\nplt.ylabel(\"Number of Samples\")\nplt.title(\"Cluster Distributions Across Feature Representations\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"4d1572ce-4124-4761-8ea4-bbb6125f23c6","cell_type":"markdown","source":"ELBO is still massive due to non-Gaussian data. ","metadata":{}},{"id":"0e672fe7-4180-4f31-9a40-6a10308df559","cell_type":"code","source":"# ✅ Load features\nfeatures_green = np.load(\"green_channel_features.npy\")\nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\")\nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")\n\n# ✅ Check for NaNs\nprint(\"Checking for NaNs...\")\nfeatures_green = np.nan_to_num(features_green)\nfeatures_pseudo_rgb = np.nan_to_num(features_pseudo_rgb)\nfeatures_grayscale = np.nan_to_num(features_grayscale)\n\n# ✅ Standardize features\nscaler = StandardScaler()\nfeatures_green = scaler.fit_transform(features_green)\nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb)\nfeatures_grayscale = scaler.transform(features_grayscale)\n\n# ✅ DP-GMM Model Setup (No Restart Between Iterations)\ndpgmm = BayesianGaussianMixture(\n    n_components=100,  # Maximum possible clusters\n    weight_concentration_prior=1e-2,  # Controls sparsity\n    weight_concentration_prior_type=\"dirichlet_process\",  # ✅ Enables DP behavior\n    covariance_type=\"full\",\n    reg_covar=1e-5,  # Prevents singularities\n    n_init=1,  # Single initialization for tracking\n    max_iter=1,  # Run one iteration at a time\n    warm_start=True,  # ✅ Keeps improving model instead of resetting\n    verbose=2,\n    random_state=42\n)\n\n# ✅ Track ELBO per iteration\nnum_iterations = 10  # Number of manual iterations\nelbo_values = []\n\nprint(\"\\nFitting DP-GMM on Green-Only Features with Manual Iteration Tracking...\")\nfor i in range(num_iterations):\n    dpgmm.fit(features_green)\n    elbo_values.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\n# ✅ Save ELBO values\nnp.save(\"dpgmm_elbo_green_mk3r.npy\", elbo_values)\n\n# 📊 Plot ELBO over iterations\nplt.figure(figsize=(8,5))\nplt.plot(range(1, num_iterations + 1), elbo_values, marker='o', linestyle='-')\nplt.xlabel(\"Iteration\")\nplt.ylabel(\"ELBO\")\nplt.title(\"ELBO Over Iterations for Green-Only Features\")\nplt.grid()\nplt.show()\n\n# ✅ Get cluster assignments\nclusters_green = dpgmm.predict(features_green)\nnp.save(\"dpgmm_clusters_greenscimk3r.npy\", clusters_green)\n\n# ✅ Repeat for Pseudo-RGB Features\nprint(\"\\nFitting DP-GMM on Pseudo-RGB Features...\")\nelbo_values_pseudo = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_pseudo_rgb)\n    elbo_values_pseudo.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_pseudo_mk3r.npy\", elbo_values_pseudo)\nclusters_pseudo_rgb = dpgmm.predict(features_pseudo_rgb)\nnp.save(\"dpgmm_clusters_pseudo_rgbscimk3r.npy\", clusters_pseudo_rgb)\n\n# ✅ Repeat for Grayscale Features\nprint(\"\\nFitting DP-GMM on Grayscale Features...\")\nelbo_values_gray = []\nfor i in range(num_iterations):\n    dpgmm.fit(features_grayscale)\n    elbo_values_gray.append(dpgmm.lower_bound_)\n    print(f\"Iteration {i+1}, ELBO: {dpgmm.lower_bound_}\")\n\nnp.save(\"dpgmm_elbo_gray_mk3r.npy\", elbo_values_gray)\nclusters_grayscale = dpgmm.predict(features_grayscale)\nnp.save(\"dpgmm_clusters_grayscalescimk3r.npy\", clusters_grayscale)\n\n# 📊 Compare Cluster Distributions\nplt.figure(figsize=(10,5))\nplt.hist(clusters_green, bins=50, alpha=0.6, label=\"Green-Only\", color=\"green\")\nplt.hist(clusters_pseudo_rgb, bins=50, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\")\nplt.hist(clusters_grayscale, bins=50, alpha=0.6, label=\"Grayscale 3-CH\", color=\"gray\")\nplt.xlabel(\"Cluster ID\")\nplt.ylabel(\"Number of Samples\")\nplt.title(\"Cluster Distributions Across Feature Representations\")\nplt.legend()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"66230654-05b9-4e9a-90df-5a47bc9f0301","cell_type":"markdown","source":"ELBO is still massive due to non-Gaussian data.","metadata":{}},{"id":"2566e38f-0ec2-4964-89db-2a478e2b807e","cell_type":"markdown","source":"---\n\n# Section 5: K-Means\n\n\nK-Means requires a k-value, or the number of clusters to be set before running. K-Means is a hard clustering algorithm, which means each datapoint can only be in one cluster. For this reason the GMM was tested beforehand. In Chapter 10 of the Stanford CS229 Notes this algorithm is discussed in great detail. \n\nThe method below is not part of the notes. Here, we run the K-Means algorithm on small subsets of the dataset. We search for the \"Knee of the curve\" to find the optimal k-value rather than guessing several k-values. When the optimal k-value is found we run K-Means on the entire dataset. ","metadata":{}},{"id":"1d2a6549-eb53-4e43-b3cd-f9787a5dfcdc","cell_type":"code","source":"# ✅ Load features (Green Channel) \nfeatures_green = np.load(\"green_channel_features.npy\")  \n# ✅ Sample a smaller subset \nsubset_size = 5000 \nfeatures_green = features_green[:subset_size]  \n# ✅ Standardize the features \nscaler = StandardScaler() \nfeatures_green = scaler.fit_transform(features_green)  \n# ✅ Try different values of k using the Elbow Method \nwcss = []  # Within-cluster sum of squares \nk_values = range(2, 15)  # Test k from 2 to 15  \nfor k in k_values:     \n    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)     \n    kmeans.fit(features_green)     \n    wcss.append(kmeans.inertia_)  # Store the sum of squared distances  \n# 📊 Plot Elbow Method results \nplt.figure(figsize=(8, 5)) \nplt.plot(k_values, wcss, marker='o', linestyle='-') \nplt.xlabel(\"Number of Clusters (k)\") \nplt.ylabel(\"WCSS (Within-Cluster Sum of Squares)\") \nplt.title(\"Elbow Method for Optimal k\") \nplt.grid() \nplt.show() ","metadata":{},"outputs":[],"execution_count":null},{"id":"e8206278-06f7-4532-b2d8-935823652d16","cell_type":"code","source":"# ✅ Load features \nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\") \nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")  \n# ✅ Sample a subset (to speed up computation) \nsubset_size = 5000 \nfeatures_pseudo_rgb = features_pseudo_rgb[:subset_size] \nfeatures_grayscale = features_grayscale[:subset_size]  \n# ✅ Standardize the features \nscaler = StandardScaler() \nfeatures_pseudo_rgb = scaler.fit_transform(features_pseudo_rgb) \nfeatures_grayscale = scaler.transform(features_grayscale)  \n# ✅ Try different values of k using the Elbow Method \nk_values = range(2, 15)  # Testing k from 2 to 15  \n# 📊 Run K-Means for Pseudo-RGB \nwcss_pseudo = [] \nfor k in k_values:     \n    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)     \n    kmeans.fit(features_pseudo_rgb)     \n    wcss_pseudo.append(kmeans.inertia_)  \n# 📊 Run K-Means for Grayscale \nwcss_gray = [] \nfor k in k_values:     \n    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)     \n    kmeans.fit(features_grayscale)     \n    wcss_gray.append(kmeans.inertia_)  \n# ✅ Plot the Elbow Method results for both \nplt.figure(figsize=(10, 5))  \nplt.subplot(1, 2, 1) \nplt.plot(k_values, wcss_pseudo, marker='o', linestyle='-') \nplt.xlabel(\"Number of Clusters (k)\") \nplt.ylabel(\"WCSS (Within-Cluster Sum of Squares)\") \nplt.title(\"Elbow Method for Pseudo-RGB Features\") \nplt.grid()  \nplt.subplot(1, 2, 2) \nplt.plot(k_values, wcss_gray, marker='o', linestyle='-') \nplt.xlabel(\"Number of Clusters (k)\") \nplt.ylabel(\"WCSS (Within-Cluster Sum of Squares)\") \nplt.title(\"Elbow Method for Grayscale Features\") \nplt.grid()  \nplt.tight_layout() \nplt.show() ","metadata":{},"outputs":[],"execution_count":null},{"id":"78c84d8a-94d3-4418-a286-950d12c21599","cell_type":"markdown","source":"Now the cluster k-values are chosen. They are k=7 for Green, k=8 for Pseudo-RGB, and k=9 for Grayscale. Next, the entire K-Means algorithm is run. The silhoutte scores are obtained and discussed below. ","metadata":{}},{"id":"b33327e7-5d89-4ee4-8606-5bfff83dd29f","cell_type":"code","source":"# ✅ Load features \nfeatures_green = np.load(\"green_channel_features.npy\") \nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\") \nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")  \n# ✅ Sample a subset (to speed up computation) \nsubset_size = 5000 \nfeatures_green = features_green[:subset_size] \nfeatures_pseudo_rgb = features_pseudo_rgb[:subset_size] \nfeatures_grayscale = features_grayscale[:subset_size]  \n# ✅ Standardize the features \nscaler = StandardScaler() \nfeatures_green = scaler.fit_transform(features_green) \nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb) \nfeatures_grayscale = scaler.transform(features_grayscale)  \n# ✅ Run K-Means for Green Channel (k=7) \nkmeans_green = KMeans(n_clusters=7, random_state=42, n_init=10) \nclusters_green = kmeans_green.fit_predict(features_green) \nsilhouette_green = silhouette_score(features_green, clusters_green)  \n# ✅ Run K-Means for Pseudo-RGB (k=8) \nkmeans_pseudo = KMeans(n_clusters=8, random_state=42, n_init=10) \nclusters_pseudo = kmeans_pseudo.fit_predict(features_pseudo_rgb) \nsilhouette_pseudo = silhouette_score(features_pseudo_rgb, clusters_pseudo)  \n# ✅ Run K-Means for Grayscale (k=9) \nkmeans_gray = KMeans(n_clusters=9, random_state=42, n_init=10) \nclusters_gray = kmeans_gray.fit_predict(features_grayscale) \nsilhouette_gray = silhouette_score(features_grayscale, clusters_gray)  \n# ✅ Save cluster assignments \nnp.save(\"kmeans_clusters_greentest.npy\", clusters_green) \nnp.save(\"kmeans_clusters_pseudotest.npy\", clusters_pseudo) \nnp.save(\"kmeans_clusters_graytest.npy\", clusters_gray)  \n# ✅ Print Silhouette Scores \nprint(f\"Silhouette Score (Green): {silhouette_green}\") \nprint(f\"Silhouette Score (Pseudo-RGB): {silhouette_pseudo}\") \nprint(f\"Silhouette Score (Grayscale): {silhouette_gray}\")  \n# 📊 Compare Cluster Distributions \nplt.figure(figsize=(10,5)) \nplt.hist(clusters_green, bins=7, alpha=0.6, label=\"Green-Only\", color=\"green\") \nplt.hist(clusters_pseudo, bins=8, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\") \nplt.hist(clusters_gray, bins=9, alpha=0.6, label=\"Grayscale\", color=\"gray\") \nplt.xlabel(\"Cluster ID\") \nplt.ylabel(\"Number of Samples\") \nplt.title(\"K-Means Cluster Distributions Across Feature Representations\") \nplt.legend() \nplt.show() ","metadata":{},"outputs":[],"execution_count":null},{"id":"7fb93872-b45d-4467-b309-a86a55925ab2","cell_type":"markdown","source":"Silhoutte scores vary between -1 and +1. It must be cohesive(or similar) to datapoints within its own cluster and seperated(or different) from the other clusters. \n\n- Greater than 0.7: Strong\n- Over 0.5: Reasonable\n- Over 0.25: Weak\n\nAs can be seen, the above silhoutte scores are all extremely small, hence K-Means did not do a good job of clustering.","metadata":{}},{"id":"b43cbe29-9193-4d84-8440-db2b63919896","cell_type":"markdown","source":"---\n\n# Section 6: DBSCAN","metadata":{}},{"id":"b139cde6-8ca6-4563-b85d-a28925029fe9","cell_type":"code","source":"features_green = np.load(\"green_channel_features.npy\") \nfeatures_pseudo_rgb = np.load(\"pseudo_rgb_features.npy\") \nfeatures_grayscale = np.load(\"grayscale_3ch_features.npy\")  \n# ✅ Sample a subset (to speed up computation) \nsubset_size = 5000 \nfeatures_green = features_green[:subset_size] \nfeatures_pseudo_rgb = features_pseudo_rgb[:subset_size] \nfeatures_grayscale = features_grayscale[:subset_size]  \n# ✅ Standardize the features (DBSCAN is distance-based, so scaling is crucial) \nscaler = StandardScaler() \nfeatures_green = scaler.fit_transform(features_green) \nfeatures_pseudo_rgb = scaler.transform(features_pseudo_rgb) \nfeatures_grayscale = scaler.transform(features_grayscale)  \n# ✅ Set DBSCAN parameters (Will fine-tune later if needed) \neps_value = 0.5  # Controls neighborhood size (will adjust if needed) \nmin_samples_value = 10  # Minimum points per dense region  \n# ✅ Run DBSCAN for Green Channel \ndbscan_green = DBSCAN(eps=eps_value, min_samples=min_samples_value, n_jobs=-1) \nclusters_green = dbscan_green.fit_predict(features_green) \nsilhouette_green = silhouette_score(features_green, clusters_green) if len(set(clusters_green)) > 1 else -1  \n# ✅ Run DBSCAN for Pseudo-RGB \ndbscan_pseudo = DBSCAN(eps=eps_value, min_samples=min_samples_value, n_jobs=-1) \nclusters_pseudo = dbscan_pseudo.fit_predict(features_pseudo_rgb) \nsilhouette_pseudo = silhouette_score(features_pseudo_rgb, clusters_pseudo) if len(set(clusters_pseudo)) > 1 else -1  \n# ✅ Run DBSCAN for Grayscale \ndbscan_gray = DBSCAN(eps=eps_value, min_samples=min_samples_value, n_jobs=-1) \nclusters_gray = dbscan_gray.fit_predict(features_grayscale) \nsilhouette_gray = silhouette_score(features_grayscale, clusters_gray) if len(set(clusters_gray)) > 1 else -1  \n# ✅ Save cluster assignments \nnp.save(\"dbscan_clusters_green.npy\", clusters_green) \nnp.save(\"dbscan_clusters_pseudo.npy\", clusters_pseudo) \nnp.save(\"dbscan_clusters_gray.npy\", clusters_gray)  \n# ✅ Print Silhouette Scores \nprint(f\"DBSCAN Silhouette Score (Green): {silhouette_green}\") \nprint(f\"DBSCAN Silhouette Score (Pseudo-RGB): {silhouette_pseudo}\") \nprint(f\"DBSCAN Silhouette Score (Grayscale): {silhouette_gray}\")  \n# 📊 Compare Cluster Distributions (excluding noise points labeled as -1) \nplt.figure(figsize=(10,5)) \nplt.hist(clusters_green[clusters_green != -1], bins=20, alpha=0.6, label=\"Green-Only\", color=\"green\") \nplt.hist(clusters_pseudo[clusters_pseudo != -1], bins=20, alpha=0.6, label=\"Pseudo-RGB\", color=\"red\") \nplt.hist(clusters_gray[clusters_gray != -1], bins=20, alpha=0.6, label=\"Grayscale\", color=\"gray\") \nplt.xlabel(\"Cluster ID\") \nplt.ylabel(\"Number of Samples\") \nplt.title(\"DBSCAN Cluster Distributions Across Feature Representations\") \nplt.legend() \nplt.show()  ","metadata":{},"outputs":[],"execution_count":null},{"id":"ba74c90a-488d-4837-ab04-752a317f75cf","cell_type":"markdown","source":"As can be seen, DBSCAN clustered everything as noise, meaning that this dataset cannot be clustered by any of these algorithms. Some other algorithm must be used to train on this data, and my original hypothesis will not work. ","metadata":{}},{"id":"cf8ea0f4-34b4-4de6-a054-c707f94f8d82","cell_type":"markdown","source":"---\n# Section 7: Citations\n\n- Stanford CS229 Lecture Notes:\n   https://cs229.stanford.edu/main_notes.pdf\n\n- DPGMM\n    https://mlg.eng.cam.ac.uk/pub/pdf/GoeRas10.pdf\n\n- Silhoutte Scores:\n    https://en.wikipedia.org/wiki/Silhouette_(clustering)\n\n- DBSCAN\n    https://en.wikipedia.org/wiki/DBSCAN\n\n- Human Protein Atlas:\n    https://www.kaggle.com/competitions/hpa-single-cell-image-classification","metadata":{}}]}