{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":18647,"databundleVersionId":1126921,"sourceType":"competition"},{"sourceId":13656575,"sourceType":"datasetVersion","datasetId":8682287}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# HMRF Probabilistic Model applied on PANDA Kaggle dataset\n\n# Importing all the necessary Python libraries\nimport os, sys, numpy as np, pandas as pd\nfrom pathlib import Path\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nimport warnings\nimport openslide\nimport cv2\nimport mahotas\nfrom skimage.color import rgb2gray, rgb2hsv\nfrom skimage.filters import threshold_otsu, gaussian\nfrom skimage.feature import local_binary_pattern\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import log_loss, accuracy_score, confusion_matrix, cohen_kappa_score\nfrom scipy.stats import entropy\nfrom sklearn.cluster import KMeans\nwarnings.filterwarnings('ignore')\n\nos.system('apt-get update -qq')\nos.system('apt-get install -y -qq openslide-tools > /dev/null')\nos.system('pip install -q openslide-python scikit-image mahotas tqdm > /dev/null')\n\n# Loading the dataset\ndata_path = Path(\"/kaggle/input/prostate-cancer-grade-assessment\")\ntrain_csv = data_path / \"train.csv\"\n\n# Handling if the dataset is missing \nif not train_csv.exists():\n    raise FileNotFoundError(f\"train.csv not found at {train_csv}. Make sure dataset is mounted to /kaggle/input/...\")\n\ntrain_df = pd.read_csv(train_csv)\nprint(f\"Loaded {len(train_df)} training samples\")\nprint(\"ISUP grade distribution:\")\nprint(train_df['isup_grade'].value_counts().sort_index())\n\n\ndef safe_array(x):\n    return np.nan_to_num(np.array(x, dtype=float), nan=0.0, posinf=0.0, neginf=0.0)\n\nEPS = 1e-10\n\n#Patch Extraction from the slides\ndef open_slide(path):\n    return openslide.OpenSlide(str(path))\n\n# Function to find the tissue regions from the slides\n    def find_tissue_regions(slide, thumb_size=2048, min_tissue_area=0.01):\n    thumb = np.array(slide.get_thumbnail((thumb_size, thumb_size)))[:,:,:3]\n    gray = rgb2gray(thumb)\n    gray = gaussian(gray, sigma=2)\n    try:\n        thresh = threshold_otsu(gray)\n        mask = gray < thresh\n    except Exception:\n        mask = gray < 0.9\n    from skimage.measure import label, regionprops\n    labeled = label(mask)\n    regions = regionprops(labeled)\n    W, H = slide.dimensions\n    scale_x, scale_y = W / thumb_size, H / thumb_size\n    tissue_coords = []\n    total_area = thumb_size * thumb_size\n    for region in regions:\n        if region.area / total_area > min_tissue_area:\n            cy, cx = region.centroid\n            full_x = int(cx * scale_x)\n            full_y = int(cy * scale_y)\n            tissue_coords.append((full_x, full_y))\n    if not tissue_coords:\n        tissue_coords = [(W//2, H//2)]\n    return tissue_coords\n\n\n# Finction to extract adaptive patches\ndef extract_adaptive_patches(slide, tissue_coords, n_patches=25, patch_size=256):\n    W, H = slide.dimensions\n    patches_list = []\n    coords_list = []\n    grid_size = int(np.sqrt(n_patches))\n    for center_x, center_y in tissue_coords[:3]:\n        stride = 128\n        radius = grid_size // 2\n        for dy in range(-radius, radius + 1):\n            for dx in range(-radius, radius + 1):\n                x = center_x + dx * stride\n                y = center_y + dy * stride\n                if x < 0 or y < 0 or x + patch_size > W or y + patch_size > H:\n                    continue\n                try:\n                    region = slide.read_region((x, y), 0, (patch_size, patch_size))\n                    patch = np.array(region.convert('RGB'))\n                except Exception:\n                    continue\n                if np.mean(patch) < 240:\n                    patches_list.append(patch)\n                    coords_list.append((x, y))\n    return patches_list, coords_list\n\n# Feature Extraction\nfrom skimage.filters import gaussian as sk_gaussian\n\n# Function to extract the enhanced features\ndef extract_enhanced_features(patch):\n    features = {}\n    def safe_stat(arr, stat_func, default=0.0):\n        try:\n            arr = np.array(arr, dtype=float)\n            if arr.size == 0:\n                return default\n            val = stat_func(arr)\n            if np.isnan(val) or np.isinf(val):\n                return default\n            return float(val)\n        except Exception:\n            return default\n\n    for i, channel in enumerate(['r','g','b']):\n        ch = patch[:,:,i].astype(float)\n        features[f'{channel}_mean'] = safe_stat(ch, np.mean, 128.0)\n        features[f'{channel}_std']  = safe_stat(ch, np.std, 10.0)\n        features[f'{channel}_median'] = safe_stat(ch, np.median, 128.0)\n        features[f'{channel}_q25'] = safe_stat(ch, lambda x: np.percentile(x,25), 100.0)\n        features[f'{channel}_q75'] = safe_stat(ch, lambda x: np.percentile(x,75), 150.0)\n\n    try:\n        hsv = rgb2hsv(patch)\n        for i, channel in enumerate(['h','s','v']):\n            ch = hsv[:,:,i]\n            features[f'{channel}_mean'] = safe_stat(ch, np.mean, 0.5)\n            features[f'{channel}_std'] = safe_stat(ch, np.std, 0.1)\n    except Exception:\n        for channel in ['h','s','v']:\n            features[f'{channel}_mean'] = 0.5\n            features[f'{channel}_std']  = 0.1\n\n    try:\n        gray = cv2.cvtColor(patch, cv2.COLOR_RGB2GRAY)\n    except:\n        gray = np.mean(patch, axis=2).astype(np.uint8)\n\n    try:\n        haralick_feats = mahotas.features.haralick(gray).mean(axis=0)\n        for i, val in enumerate(haralick_feats[:10]):\n            features[f'haralick_{i}'] = safe_stat(val, lambda x: x, 0.0)\n    except:\n        for i in range(10):\n            features[f'haralick_{i}'] = 0.0\n\n    try:\n        radius = 3\n        n_points = 8 * radius\n        lbp = local_binary_pattern(gray, n_points, radius, method='uniform')\n        lbp_hist, _ = np.histogram(lbp.ravel(), bins=n_points + 2, range=(0, n_points + 2))\n        lbp_hist = lbp_hist.astype(float) / (lbp_hist.sum() + EPS)\n        for i, val in enumerate(lbp_hist[:10]):\n            features[f'lbp_{i}'] = safe_stat(val, lambda x: x, 0.1)\n    except:\n        for i in range(10):\n            features[f'lbp_{i}'] = 0.1\n\n    # Edges\n    try:\n        edges = cv2.Canny(gray, 50, 150)\n        features['edge_density'] = safe_stat(edges/255.0, np.mean, 0.1)\n    except:\n        features['edge_density'] = 0.1\n\n    # Focus measure\n    try:\n        lap_var = cv2.Laplacian(gray, cv2.CV_64F).var()\n        features['focus_lap'] = safe_stat(lap_var, lambda x: x, 100.0)\n    except:\n        features['focus_lap'] = 100.0\n\n    # Tissue density\n    try:\n        tissue_mask = np.mean(patch, axis=2) < 240\n        features['tissue_density'] = safe_stat(tissue_mask, np.mean, 0.5)\n    except:\n        features['tissue_density'] = 0.5\n\n    # Color contrast\n    try:\n        features['color_contrast'] = safe_stat(patch, lambda x: np.std(x, axis=(0,1)).mean(), 20.0)\n    except:\n        features['color_contrast'] = 20.0\n\n    # Spatial gradient\n    try:\n        gx = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)\n        gy = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=3)\n        grad_mag = np.sqrt(gx**2 + gy**2)\n        features['grad_mean'] = safe_stat(grad_mag, np.mean, 10.0)\n        features['grad_std'] = safe_stat(grad_mag, np.std, 5.0)\n    except:\n        features['grad_mean'] = 10.0\n        features['grad_std'] = 5.0\n\n    return features\n\n\n# Function to build the feature matrix\ndef build_feature_matrix(patches):\n    feature_list = []\n    for patch in tqdm(patches, desc=\"Extracting features\"):\n        features = extract_enhanced_features(patch)\n        feature_list.append(features)\n    df = pd.DataFrame(feature_list)\n\n    # Handling NaN values\n    df = df.replace([np.inf, -np.inf], np.nan)\n    for col in df.columns:\n        col_mean = df[col].mean()\n        if np.isnan(col_mean) or np.isinf(col_mean):\n            df[col] = df[col].fillna(0.0)\n        else:\n            df[col] = df[col].fillna(col_mean)\n\n    df = df.fillna(0.0)\n\n    X_raw = df.values.astype(float)\n    X_raw = np.nan_to_num(X_raw, nan=0.0, posinf=0.0, neginf=0.0)\n\n    scaler = StandardScaler()\n    if X_raw.shape[0] == 0:\n        return np.zeros((0,0)), [], scaler\n    try:\n        X = scaler.fit_transform(X_raw)\n    except Exception:\n        X = X_raw - np.mean(X_raw, axis=0, keepdims=True)\n        X = np.nan_to_num(X, nan=0.0)\n    X = np.nan_to_num(X, nan=0.0, posinf=0.0, neginf=0.0)\n    return X, df.columns.tolist(), scaler\n\n# Building a spatial graph\ndef build_spatial_graph(coords, patch_size=512, connectivity=4):\n    N = len(coords)\n    neighbors = [[] for _ in range(N)]\n    coord_to_idx = {c: i for i, c in enumerate(coords)}\n    for i, (x, y) in enumerate(coords):\n        candidates = [\n            (x + patch_size, y),\n            (x - patch_size, y),\n            (x, y + patch_size),\n            (x, y - patch_size)\n        ]\n        if connectivity == 8:\n            candidates += [\n                (x + patch_size, y + patch_size),\n                (x + patch_size, y - patch_size),\n                (x - patch_size, y + patch_size),\n                (x - patch_size, y - patch_size)\n            ]\n        for coord in candidates:\n            j = coord_to_idx.get(coord, None)\n            if j is not None and j != i:\n                neighbors[i].append(j)\n    return neighbors\n\n# Initialization of Gaussian Log Likelihood\ndef initialize_gmm(X, K, method='kmeans'):\n    N, D = X.shape\n    X = np.nan_to_num(X, nan=0.0, posinf=0.0, neginf=0.0)\n    mus = np.zeros((K, D), dtype=float)\n    Sigmas = np.zeros((K, D, D), dtype=float)\n\n    if method == 'kmeans' and N >= K:\n        kmeans = KMeans(n_clusters=K, random_state=42, n_init=10)\n        labels = kmeans.fit_predict(X)\n        for k in range(K):\n            idx = np.where(labels == k)[0]\n            if len(idx) > 1:\n                mu = X[idx].mean(axis=0)\n                mu = np.nan_to_num(mu, nan=0.0)\n                mus[k] = mu\n                cov = np.cov(X[idx].T)\n                if cov.ndim == 0:\n                    cov = np.eye(D) * max(float(cov), 1e-6)\n                cov = np.nan_to_num(cov, nan=0.0)\n                cov = cov + 1e-4 * np.eye(D)\n                Sigmas[k] = cov\n            elif len(idx) == 1:\n                mus[k] = np.nan_to_num(X[idx[0]], nan=0.0)\n                Sigmas[k] = np.eye(D) * 1e-3\n            else:\n                mus[k] = np.zeros(D)\n                Sigmas[k] = np.eye(D)\n    else:\n        rng = np.random.RandomState(42)\n        if N >= K:\n            indices = rng.choice(N, K, replace=False)\n        else:\n            indices = rng.choice(N, K, replace=True)\n        base_cov = np.cov(X.T) if N > 1 else np.eye(X.shape[1])\n        base_cov = np.nan_to_num(base_cov, nan=0.0)\n        base_cov = base_cov + 1e-4 * np.eye(D)\n        for k, idx in enumerate(indices):\n            mus[k] = np.nan_to_num(X[idx], nan=0.0)\n            Sigmas[k] = base_cov.copy()\n\n    pi = np.ones(K) / K\n    return mus, Sigmas, pi\n\n\ndef gaussian_log_likelihood(x, mu, Sigma, min_loglik=-1e10):\n    x = safe_array(x)\n    mu = safe_array(mu)\n    Sigma = np.array(Sigma, dtype=float)\n    if np.any(np.isnan(x)) or np.any(np.isnan(mu)) or np.any(np.isnan(Sigma)):\n        return min_loglik\n    D = x.shape[0]\n    Sigma = Sigma + 1e-6 * np.eye(D)\n    try:\n        # Use slogdet for determinant and solve for Mahalanobis\n        sign, logdet = np.linalg.slogdet(Sigma)\n        if sign <= 0 or np.isnan(logdet) or np.isinf(logdet):\n            Sigma = Sigma + 1e-3 * np.eye(D)\n            sign, logdet = np.linalg.slogdet(Sigma)\n            if sign <= 0:\n                return min_loglik\n        diff = x - mu\n        sol = np.linalg.solve(Sigma, diff)\n        mahalanobis = float(diff @ sol)\n        if np.isnan(mahalanobis) or np.isinf(mahalanobis):\n            return min_loglik\n        loglik = -0.5 * (mahalanobis + logdet + D * np.log(2 * np.pi))\n        loglik = np.maximum(loglik, min_loglik)\n        return float(loglik)\n    except Exception:\n        return min_loglik\n\n\n# Creating a robust HMRF Model\ndef hmrf_em_icm(X, neighbors, K=6, beta=1.0, max_iter=30, tol=1e-4, verbose=True):\n    N, D = X.shape\n    X = np.nan_to_num(X, nan=0.0, posinf=0.0, neginf=0.0)\n    if N == 0:\n        return np.array([], dtype=int), np.zeros((0, K)), np.zeros((K, D)), np.zeros((K, D, D))\n    mus, Sigmas, pi = initialize_gmm(X, K, method='kmeans')\n    labels = np.random.randint(0, K, size=N)\n    prev_energy = -np.inf\n\n    for it in range(max_iter):\n        for k in range(K):\n            idx = np.where(labels == k)[0]\n            if len(idx) > 1:\n                mu = X[idx].mean(axis=0)\n                mu = np.nan_to_num(mu, nan=0.0)\n                mus[k] = mu\n                cov = np.cov(X[idx].T)\n                if cov.ndim == 0:\n                    cov = np.eye(D) * max(float(cov), 1e-6)\n                cov = np.nan_to_num(cov, nan=0.0)\n                cov = cov + 1e-4 * np.eye(D)\n                Sigmas[k] = cov\n            elif len(idx) == 1:\n                mus[k] = np.nan_to_num(X[idx[0]], nan=0.0)\n                Sigmas[k] = np.eye(D) * 1e-3\n            else:\n                Sigmas[k] = Sigmas[k] + 1e-4 * np.eye(D)\n        changed = 0\n        energy = 0.0\n        order = np.random.permutation(N)\n        for i in order:\n            best_label = labels[i]\n            best_score = -np.inf\n            for k in range(K):\n                data_ll = gaussian_log_likelihood(X[i], mus[k], Sigmas[k])\n                spatial_prior = 0\n                for j in neighbors[i]:\n                    if labels[j] == k:\n                        spatial_prior += 1\n                score = data_ll + beta * spatial_prior\n                if score > best_score:\n                    best_score = score\n                    best_label = k\n            if labels[i] != best_label:\n                labels[i] = best_label\n                changed += 1\n            energy += best_score\n\n        if verbose:\n            print(f\"Iter {it+1}/{max_iter}: changed={changed}, energy={energy:.4f}\")\n        if abs(energy - prev_energy) < tol * (abs(prev_energy) + 1e-10):\n            if verbose:\n                print(f\"Converged at iter {it+1}\")\n            break\n        prev_energy = energy\n        if changed == 0:\n            if verbose:\n                print(\"No label changes, stopping early\")\n            break\n\n    log_post = np.full((N, K), -1e12, dtype=float)\n    for k in range(K):\n        for i in range(N):\n            log_post[i, k] = gaussian_log_likelihood(X[i], mus[k], Sigmas[k], min_loglik=-1e12)\n\n    max_ll = np.max(log_post, axis=1, keepdims=True)\n    max_ll = np.nan_to_num(max_ll, nan=0.0, posinf=0.0, neginf=0.0)\n    exp_ll = np.exp(log_post - max_ll)\n    exp_ll = np.nan_to_num(exp_ll, nan=0.0, posinf=0.0, neginf=0.0)\n\n    row_sums = exp_ll.sum(axis=1, keepdims=True)\n    zero_rows = (row_sums.squeeze() == 0)\n    if np.any(zero_rows):\n        exp_ll[zero_rows, :] = 1.0 / K\n        row_sums = exp_ll.sum(axis=1, keepdims=True)\n\n    posteriors = exp_ll / (row_sums + EPS)\n    posteriors = np.nan_to_num(posteriors, nan=1.0/K, posinf=1.0/K, neginf=0.0)\n\n    return labels, posteriors, mus, Sigmas\n\n# Function to generate the aggregation metrics\ndef aggregate_patch_posteriors(posteriors, method='weighted_mean', weights=None):\n    posteriors = np.nan_to_num(posteriors, nan=0.0)\n    row_sums = posteriors.sum(axis=1, keepdims=True)\n    row_sums = np.where(row_sums == 0, 1.0, row_sums)\n    posteriors = posteriors / row_sums\n\n    K = posteriors.shape[1]\n    if method == 'mean':\n        result = posteriors.mean(axis=0)\n    elif method == 'weighted_mean':\n        if weights is None:\n            weights = posteriors.max(axis=1)\n        weights = np.nan_to_num(weights, nan=1.0)\n        weight_sum = weights.sum()\n        if weight_sum == 0:\n            weight_sum = 1.0\n        weights = weights / weight_sum\n        result = (posteriors.T @ weights)\n    elif method == 'max':\n        votes = posteriors.argmax(axis=1)\n        result = np.bincount(votes, minlength=K)\n        result = result / (result.sum() + EPS)\n    else:\n        result = posteriors.mean(axis=0)\n\n    result = np.nan_to_num(result, nan=1.0/K)\n    result = result / (result.sum() + EPS)\n    return result\n\ndef compute_metrics(y_true, y_pred_probs, verbose=True):\n    y_true = np.array(y_true).astype(int)\n    y_pred_probs = np.nan_to_num(y_pred_probs, nan=0.0)\n    row_sums = y_pred_probs.sum(axis=1, keepdims=True)\n    row_sums = np.where(row_sums == 0, 1.0, row_sums)\n    y_pred_probs = y_pred_probs / row_sums\n\n    K = y_pred_probs.shape[1]\n    y_pred_map = y_pred_probs.argmax(axis=1)\n    valid_mask = (y_true >= 0) & (y_true < K)\n    y_true = y_true[valid_mask]\n    y_pred_probs = y_pred_probs[valid_mask]\n    y_pred_map = y_pred_map[valid_mask]\n\n    if len(y_true) == 0:\n        if verbose:\n            print(\"No valid samples for evaluation\")\n        return None\n\n    acc = accuracy_score(y_true, y_pred_map)\n    qwk = cohen_kappa_score(y_true, y_pred_map, weights='quadratic')\n    epsilon = 1e-12\n    y_pred_probs_clipped = np.clip(y_pred_probs, epsilon, 1 - epsilon)\n    y_pred_probs_clipped = y_pred_probs_clipped / y_pred_probs_clipped.sum(axis=1, keepdims=True)\n    try:\n        ll = log_loss(y_true, y_pred_probs_clipped, labels=list(range(K)))\n    except ValueError as e:\n        print(\"Warning in log_loss:\", e)\n        probs_uniform = np.ones_like(y_pred_probs_clipped) / K\n        ll = log_loss(y_true, probs_uniform, labels=list(range(K)))\n\n    N = len(y_true)\n    onehot = np.zeros((N, K))\n    onehot[np.arange(N), y_true] = 1\n    brier = np.mean(np.sum((y_pred_probs - onehot)**2, axis=1))\n    cm = confusion_matrix(y_true, y_pred_map, labels=list(range(K)))\n\n    # Print all the evaluation metrics\n    if verbose:\n        print(f\"Samples evaluated: {len(y_true)}\")\n        print(f\"Accuracy: {acc:.4f}\")\n        print(f\"Quadratic Weighted Kappa: {qwk:.4f}\")\n        print(f\"Log Loss: {ll:.4f}\")\n        print(f\"Brier Score: {brier:.4f}\")\n        print(f\"\\nConfusion Matrix:\\n{cm}\")\n\n    return {\n        'accuracy': acc,\n        'qwk': qwk,\n        'log_loss': ll,\n        'brier': brier,\n        'confusion_matrix': cm,\n        'n_samples': len(y_true)\n    }\n\n# Process all the slides\ndef process_slide(slide_id, data_path, n_patches=49, patch_size=512):\n    slide_path = data_path / \"train_images\" / f\"{slide_id}.tiff\"\n    if not slide_path.exists():\n        print(f\"Slide {slide_id} not found at {slide_path}\")\n        return None\n    try:\n        slide = open_slide(slide_path)\n        tissue_coords = find_tissue_regions(slide)\n        patches, coords = extract_adaptive_patches(slide, tissue_coords, n_patches=n_patches, patch_size=patch_size)\n        if len(patches) < 3:\n            slide.close()\n            return None\n        X, feature_names, scaler = build_feature_matrix(patches)\n        if np.any(np.isnan(X)):\n            X = np.nan_to_num(X, nan=0.0)\n        neighbors = build_spatial_graph(coords, patch_size=patch_size)\n        slide.close()\n        return {\n            'X': X,\n            'coords': coords,\n            'neighbors': neighbors,\n            'n_patches': len(patches),\n            'scaler': scaler\n        }\n    except Exception as e:\n        print(f\"Error processing {slide_id}: {e}\")\n        return None\n\n# Running the main HMRF Pipeline\ndef run_hmrf_pipeline(n_slides=50, K=6, beta=1.0, aggregation='weighted_mean', verbose=True):\n    print(\"\\n\" + \"=\"*60)\n    print(f\"Running HMRF Pipeline: K={K}, beta={beta}, n_slides={n_slides}\")\n    print(\"=\"*60 + \"\\n\")\n    results = []\n    skipped = 0\n    y_true_list = []\n    y_pred_probs_list = []\n\n    max_idx = min(n_slides, len(train_df))\n    for idx in tqdm(range(max_idx), desc=\"Processing slides\"):\n        row = train_df.iloc[idx]\n        slide_id = row['image_id']\n        true_grade = int(row['isup_grade'])\n        if true_grade < 0 or true_grade >= K:\n            skipped += 1\n            continue\n        slide_data = process_slide(slide_id, data_path)\n        if slide_data is None:\n            skipped += 1\n            continue\n        try:\n            labels, posteriors, mus, Sigmas = hmrf_em_icm(\n                slide_data['X'], slide_data['neighbors'], K=K, beta=beta, max_iter=20, verbose=False)\n\n            posteriors = np.nan_to_num(posteriors, nan=0.0)\n            row_sums = posteriors.sum(axis=1, keepdims=True)\n            row_sums = np.where(row_sums == 0, 1.0, row_sums)\n            posteriors = posteriors / row_sums\n\n            if posteriors.size == 0:\n                print(f\"Warning: empty posteriors for {slide_id}, skipping\")\n                skipped += 1\n                continue\n\n            slide_posterior = aggregate_patch_posteriors(posteriors, method=aggregation)\n            if np.any(np.isnan(slide_posterior)) or slide_posterior.sum() == 0:\n                slide_posterior = np.ones(K) / K\n\n            map_grade = int(np.argmax(slide_posterior))\n            expected_grade = float(np.sum(np.arange(K) * slide_posterior))\n            uncertainty = float(entropy(slide_posterior))\n\n            results.append({\n                'slide_id': slide_id,\n                'true_grade': true_grade,\n                'posterior': slide_posterior,\n                'map_grade': map_grade,\n                'expected_grade': expected_grade,\n                'uncertainty': uncertainty,\n                'n_patches': slide_data['n_patches']\n            })\n\n            # for metrics\n            y_true_list.append(true_grade)\n            y_pred_probs_list.append(slide_posterior)\n\n        except Exception as e:\n            print(f\"Error in HMRF for {slide_id}: {e}\")\n            skipped += 1\n            continue\n\n    if len(y_true_list) > 0:\n        y_pred_probs_arr = np.vstack(y_pred_probs_list)\n        metrics = compute_metrics(y_true_list, y_pred_probs_arr, verbose=verbose)\n    else:\n        metrics = None\n\n    print(f\"\\nProcessed: {len(results)}, Skipped: {skipped}\")\n    return results, metrics\n\n# Performing the grid search\ndef grid_search(n_slides=30, Ks=[4,5,6], betas=[0.5,1.0,1.5], aggregation='weighted_mean'):\n    best_qwk = -np.inf\n    best_params = None\n    grid_results = []\n    for K in Ks:\n        for beta in betas:\n            print(f\"\\n--- Grid try K={K}, beta={beta} ---\")\n            results, metrics = run_hmrf_pipeline(n_slides=n_slides, K=K, beta=beta, aggregation=aggregation, verbose=False)\n            qwk = metrics['qwk'] if metrics is not None else -np.inf\n            grid_results.append({'K':K, 'beta':beta, 'qwk':qwk, 'metrics': metrics})\n            print(f\"Result: qwk={qwk:.4f}\")\n            if qwk > best_qwk:\n                best_qwk = qwk\n                best_params = {'K':K, 'beta':beta}\n    print(f\"\\nBest params: {best_params}, best_qwk={best_qwk:.4f}\")\n    return best_params, grid_results\n\n# Running an example\nif __name__ == \"__main__\":\n    best_params, grid_results = grid_search(\n        n_slides=100,            \n        Ks=[4,5,6,7],            \n        betas=[0.6,0.8,1.0,1.2], \n        aggregation='weighted_mean'\n    )\n    print(\"Quick grid best:\", best_params)\n\n    \n    results, metrics = run_hmrf_pipeline(\n        n_slides=1000,                 \n        K=best_params['K'],\n        beta=best_params['beta'],\n        aggregation='weighted_mean',\n        verbose=True\n    )\n    print(\"\\nSummary metrics:\")\n    print(metrics)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T11:10:30.593215Z","iopub.execute_input":"2025-11-15T11:10:30.593579Z"}},"outputs":[],"execution_count":null}]}