1{
2  "op": "replace",
3  "filename": "scikit-learn/custom_clustering.py",
4  "start_line": 36,
5  "end_line": 109,
6  "content": "import numpy as np\nfrom sklearn.base import BaseEstimator, ClusterMixin\nfrom sklearn.neighbors import NearestNeighbors\nfrom scipy.sparse import csgraph\nfrom scipy.sparse.linalg import eigs\nfrom sklearn.cluster import KMeans\n\n\nclass CustomClustering(BaseEstimator, ClusterMixin):\n    \"\"\"Adaptive Density-Aware Spectral Clustering (ADASC).\n\n    A novel clustering algorithm that combines density estimation with spectral\n    graph theory to handle diverse cluster geometries: convex, non-convex,\n    varying densities, and high-dimensional data.\n\n    Steps:\n    1. Compute local density for each point using k-NN distances.\n    2. Construct an affinity matrix with adaptive bandwidths based on local density.\n    3. Compute normalized Laplacian and its eigenvectors.\n    4. Determine number of clusters via eigengap heuristic if n_clusters is None.\n    5. Cluster eigenvectors using k-means.\n\n    Args:\n        n_clusters: Number of clusters. If None, determined automatically.\n        random_state: Random seed for reproducibility.\n        k_neighbors: Number of neighbors for local density estimation.\n        gamma: Scaling factor for adaptive bandwidth (default 0.5).\n    \"\"\"\n\n    def __init__(self, n_clusters=None, random_state=42, k_neighbors=10, gamma=0.5):\n        self.n_clusters = n_clusters\n        self.random_state = random_state\n        self.k_neighbors = k_neighbors\n        self.gamma = gamma\n        self.labels_ = None\n        self.affinity_matrix_ = None\n        self.eigenvectors_ = None\n\n    def _compute_local_density(self, X):\n        \"\"\"Compute local density as inverse of mean distance to k nearest neighbors.\"\"\"\n        n_samples = X.shape[0]\n        k = min(self.k_neighbors, n_samples - 1)\n        if k <= 0:\n            return np.ones(n_samples)\n        nn = NearestNeighbors(n_neighbors=k + 1)  # include self\n        nn.fit(X)\n        distances, _ = nn.kneighbors(X)\n        # Exclude self-distance (first column)\n        mean_dist = np.mean(distances[:, 1:], axis=1)\n        # Avoid division by zero\n        mean_dist = np.maximum(mean_dist, 1e-10)\n        density = 1.0 / mean_dist\n        return density\n\n    def _compute_affinity_matrix(self, X, densities):\n        \"\"\"Compute affinity matrix with adaptive bandwidth per point.\"\"\"\n        n_samples = X.shape[0]\n        # Pairwise squared Euclidean distances\n        D = np.sum(X**2, axis=1, keepdims=True)\n        dist_sq = D + D.T - 2 * X @ X.T\n        # Ensure non-negative (numerical errors)\n        dist_sq = np.maximum(dist_sq, 0.0)\n        # Adaptive bandwidth: sigma_i = gamma * (1 / density_i)\n        # Higher density => smaller bandwidth (tighter)\n        sigma = self.gamma / densities\n        # Broadcast to matrix: sigma_ij = sigma_i * sigma_j\n        sigma_mat = sigma[:, np.newaxis] * sigma[np.newaxis, :]\n        # Affinity: exp(-dist^2 / sigma_ij)\n        affinity = np.exp(-dist_sq / sigma_mat)\n        # Zero out diagonal for symmetry (already exp(0)=1)\n        np.fill_diagonal(affinity, 0.0)\n        return affinity\n\n    def _spectral_embedding(self, affinity, n_clusters):\n        \"\"\"Compute eigenvectors of normalized Laplacian.\"\"\"\n        n_samples = affinity.shape[0]\n        # Degree matrix\n        deg = np.sum(affinity, axis=1)\n        # Avoid zero degree\n        deg = np.maximum(deg, 1e-10)\n        # Normalized Laplacian: L = I - D^{-1/2} A D^{-1/2}\n        deg_sqrt_inv = 1.0 / np.sqrt(deg)\n        L = np.eye(n_samples) - deg_sqrt_inv[:, np.newaxis] * affinity * deg_sqrt_inv[np.newaxis, :]\n        # Compute eigenvectors corresponding to smallest eigenvalues\n        # Use eigs for sparse matrix; treat L as dense (since n_samples moderate)\n        # Compute k+1 smallest eigenvalues (first eigenvalue is zero)\n        n_eig = min(n_clusters + 1, n_samples - 1)\n        if n_eig <= 0:\n            # Fallback: return identity\n            return np.ones((n_samples, 1))\n        # Use eigs from scipy.sparse.linalg (works with dense arrays)\n        vals, vecs = eigs(L, k=n_eig, which='SR', tol=1e-4)\n        # Sort by increasing eigenvalue\n        idx = np.argsort(np.real(vals))\n        vecs = vecs[:, idx]\n        # Take eigenvectors corresponding to smallest eigenvalues, excluding the first zero eigenvalue\n        # (the first eigenvector is constant, discard)\n        embedding = np.real(vecs[:, 1:n_clusters+1])\n        # Normalize rows to unit length (optional but standard in spectral clustering)\n        row_norms = np.linalg.norm(embedding, axis=1, keepdims=True)\n        row_norms[row_norms == 0] = 1\n        embedding = embedding / row_norms\n        return embedding\n\n    def _determine_n_clusters(self, affinity):\n        \"\"\"Determine number of clusters via eigengap heuristic.\"\"\"\n        n_samples = affinity.shape[0]\n        # Compute normalized Laplacian\n        deg = np.sum(affinity, axis=1)\n        deg = np.maximum(deg, 1e-10)\n        deg_sqrt_inv = 1.0 / np.sqrt(deg)\n        L = np.eye(n_samples) - deg_sqrt_inv[:, np.newaxis] * affinity * deg_sqrt_inv[np.newaxis, :]\n        # Compute eigenvalues (smallest)\n        max_eig = min(20, n_samples - 1)\n        if max_eig <= 1:\n            return 1\n        vals, _ = eigs(L, k=max_eig, which='SR', tol=1e-4)\n        vals = np.real(vals)\n        vals = np.sort(vals)\n        # Compute gaps between consecutive eigenvalues\n        gaps = vals[1:] - vals[:-1]\n        # Skip the first gap (between zero and second eigenvalue)\n        if len(gaps) > 1:\n            # Find largest gap among the first half (to avoid overestimating)\n            search_len = len(gaps) // 2\n            idx = np.argmax(gaps[:search_len]) + 1  # +1 because we skipped zero eigenvalue\n            n_clusters = max(2, idx + 1)  # ensure at least 2 clusters\n        else:\n            n_clusters = 2\n        # Limit to reasonable number\n        n_clusters = min(n_clusters, n_samples // 2)\n        return n_clusters\n\n    def fit(self, X):\n        \"\"\"Fit the clustering model to data X.\"\"\"\n        n_samples = X.shape[0]\n        # If n_clusters is provided and is 1, trivial clustering\n        if self.n_clusters == 1:\n            self.labels_ = np.zeros(n_samples, dtype=int)\n            return self\n\n        # Compute local density\n        densities = self._compute_local_density(X)\n        # Compute affinity matrix\n        affinity = self._compute_affinity_matrix(X, densities)\n        self.affinity_matrix_ = affinity\n\n        # Determine number of clusters if not provided\n        if self.n_clusters is None:\n            n_clusters = self._determine_n_clusters(affinity)\n            # Ensure at least 2 and at most n_samples // 2\n            n_clusters = max(2, min(n_clusters, n_samples // 2))\n        else:\n            n_clusters = self.n_clusters\n\n        # Spectral embedding\n        embedding = self._spectral_embedding(affinity, n_clusters)\n\n        # Cluster embeddings using k-means\n        kmeans = KMeans(n_clusters=n_clusters, random_state=self.random_state, n_init=10)\n        labels = kmeans.fit_predict(embedding)\n        self.labels_ = labels\n        return self\n\n    def predict(self, X):\n        \"\"\"Predict cluster labels for X using nearest centroid in embedding space.\"\"\"\n        if self.labels_ is None:\n            self.fit(X)\n        # For spectral clustering, out-of-sample prediction is non-trivial.\n        # We'll fallback to refitting (stateless) for simplicity.\n        # In a more advanced version, one could use Nystrom extension.\n        # For now, we refit.\n        self.fit(X)\n        return self.labels_\n\n\ndef custom_distance(x, y):\n    \"\"\"Custom distance metric: Euclidean distance.\"\"\"\n    return np.sqrt(np.sum((x - y) ** 2))"
7}

1{
2  "op": "replace",
3  "filename": "scikit-learn/custom_clustering.py",
4  "start_line": 36,
5  "end_line": 214,
6  "content": "import numpy as np\nfrom sklearn.base import BaseEstimator, ClusterMixin\nfrom sklearn.neighbors import NearestNeighbors\nfrom scipy.sparse import csgraph, csr_matrix\nfrom scipy.sparse.linalg import eigs\nfrom sklearn.cluster import KMeans\nfrom sklearn.decomposition import PCA\n\n\nclass CustomClustering(BaseEstimator, ClusterMixin):\n    \"\"\"Mutual k-NN Spectral Clustering with Adaptive Bandwidth (MkSC-AB).\n\n    A novel algorithm that builds a mutual k-nearest neighbor graph with\n    density-adaptive edge weights, then performs spectral clustering.\n    Handles non-convex shapes, varied densities, and high-dimensional data.\n\n    Steps:\n    1. If n_features > 50, apply PCA to reduce dimensionality.\n    2. Compute local density as inverse of mean k-NN distance.\n    3. Build mutual k-NN graph (edge if i in kNN(j) and j in kNN(i)).\n    4. Edge weight = exp(-dist^2 / (sigma_i * sigma_j)), where sigma_i is\n       adaptive bandwidth proportional to local density.\n    5. Compute normalized Laplacian eigenvectors.\n    6. Determine number of clusters via eigengap heuristic if n_clusters is None.\n    7. Cluster eigenvectors using k-means.\n\n    Args:\n        n_clusters: Number of clusters. If None, determined automatically.\n        random_state: Random seed for reproducibility.\n        k_neighbors: Number of neighbors for graph construction (default 10).\n        gamma: Scaling factor for adaptive bandwidth (default 0.5).\n    \"\"\"\n\n    def __init__(self, n_clusters=None, random_state=42, k_neighbors=10, gamma=0.5):\n        self.n_clusters = n_clusters\n        self.random_state = random_state\n        self.k_neighbors = k_neighbors\n        self.gamma = gamma\n        self.labels_ = None\n        self.affinity_matrix_ = None\n        self.eigenvectors_ = None\n        self.pca_ = None\n\n    def _reduce_dimensionality(self, X):\n        \"\"\"Apply PCA if n_features > 50, preserving 95% variance.\"\"\"\n        n_samples, n_features = X.shape\n        if n_features <= 50:\n            return X\n        # Keep enough components to explain 95% variance\n        pca = PCA(n_components=0.95, random_state=self.random_state)\n        X_reduced = pca.fit_transform(X)\n        self.pca_ = pca\n        return X_reduced\n\n    def _compute_local_density(self, X):\n        \"\"\"Compute local density as inverse of mean distance to k nearest neighbors.\"\"\"\n        n_samples = X.shape[0]\n        k = min(self.k_neighbors, n_samples - 1)\n        if k <= 0:\n            return np.ones(n_samples)\n        nn = NearestNeighbors(n_neighbors=k + 1)  # include self\n        nn.fit(X)\n        distances, _ = nn.kneighbors(X)\n        # Exclude self-distance (first column)\n        mean_dist = np.mean(distances[:, 1:], axis=1)\n        mean_dist = np.maximum(mean_dist, 1e-10)\n        density = 1.0 / mean_dist\n        return density\n\n    def _build_mutual_knn_graph(self, X, densities):\n        \"\"\"Build symmetric mutual k-NN graph with adaptive edge weights.\"\"\"\n        n_samples = X.shape[0]\n        k = min(self.k_neighbors, n_samples - 1)\n        if k <= 0:\n            # Fallback to full affinity\n            D = np.sum(X**2, axis=1, keepdims=True)\n            dist_sq = D + D.T - 2 * X @ X.T\n            dist_sq = np.maximum(dist_sq, 0.0)\n            sigma = self.gamma / densities\n            sigma_mat = sigma[:, np.newaxis] * sigma[np.newaxis, :]\n            affinity = np.exp(-dist_sq / sigma_mat)\n            np.fill_diagonal(affinity, 0.0)\n            return affinity\n\n        # Find k nearest neighbors (including self)\n        nn = NearestNeighbors(n_neighbors=k + 1)\n        nn.fit(X)\n        distances, indices = nn.kneighbors(X)\n        # Exclude self\n        distances = distances[:, 1:]\n        indices = indices[:, 1:]\n\n        # Build adjacency matrix (binary mutual k-NN)\n        adjacency = np.zeros((n_samples, n_samples), dtype=bool)\n        for i in range(n_samples):\n            for j_idx, j in enumerate(indices[i]):\n                # Check mutual: if i is in neighbors of j\n                if i in indices[j]:\n                    adjacency[i, j] = True\n                    adjacency[j, i] = True\n\n        # Edge weights: Gaussian with adaptive bandwidth\n        sigma = self.gamma / densities\n        # Precompute pairwise squared distances for edges only (to save memory)\n        # We'll compute distances on the fly for mutual edges\n        affinity = np.zeros((n_samples, n_samples))\n        # For each mutual edge, compute weight\n        rows, cols = np.where(adjacency)\n        for idx in range(len(rows)):\n            i = rows[idx]\n            j = cols[idx]\n            if i < j:  # compute only once per pair\n                dist_sq = np.sum((X[i] - X[j]) ** 2)\n                weight = np.exp(-dist_sq / (sigma[i] * sigma[j]))\n                affinity[i, j] = weight\n                affinity[j, i] = weight\n        # Ensure diagonal zero\n        np.fill_diagonal(affinity, 0.0)\n        return affinity\n\n    def _spectral_embedding(self, affinity, n_clusters):\n        \"\"\"Compute eigenvectors of normalized Laplacian.\"\"\"\n        n_samples = affinity.shape[0]\n        # Degree matrix\n        deg = np.sum(affinity, axis=1)\n        deg = np.maximum(deg, 1e-10)\n        # Normalized Laplacian: L = I - D^{-1/2} A D^{-1/2}\n        deg_sqrt_inv = 1.0 / np.sqrt(deg)\n        L = np.eye(n_samples) - deg_sqrt_inv[:, np.newaxis] * affinity * deg_sqrt_inv[np.newaxis, :]\n        # Compute eigenvectors corresponding to smallest eigenvalues\n        n_eig = min(n_clusters + 1, n_samples - 1)\n        if n_eig <= 0:\n            return np.ones((n_samples, 1))\n        # Use eigs (dense matrix)\n        vals, vecs = eigs(L, k=n_eig, which='SR', tol=1e-4)\n        idx = np.argsort(np.real(vals))\n        vecs = vecs[:, idx]\n        embedding = np.real(vecs[:, 1:n_clusters+1])\n        # Normalize rows\n        row_norms = np.linalg.norm(embedding, axis=1, keepdims=True)\n        row_norms[row_norms == 0] = 1\n        embedding = embedding / row_norms\n        return embedding\n\n    def _determine_n_clusters(self, affinity):\n        \"\"\"Determine number of clusters via eigengap heuristic.\"\"\"\n        n_samples = affinity.shape[0]\n        deg = np.sum(affinity, axis=1)\n        deg = np.maximum(deg, 1e-10)\n        deg_sqrt_inv = 1.0 / np.sqrt(deg)\n        L = np.eye(n_samples) - deg_sqrt_inv[:, np.newaxis] * affinity * deg_sqrt_inv[np.newaxis, :]\n        max_eig = min(20, n_samples - 1)\n        if max_eig <= 1:\n            return 1\n        vals, _ = eigs(L, k=max_eig, which='SR', tol=1e-4)\n        vals = np.real(vals)\n        vals = np.sort(vals)\n        gaps = vals[1:] - vals[:-1]\n        if len(gaps) > 1:\n            search_len = len(gaps) // 2\n            idx = np.argmax(gaps[:search_len]) + 1\n            n_clusters = max(2, idx + 1)\n        else:\n            n_clusters = 2\n        n_clusters = min(n_clusters, n_samples // 2)\n        return n_clusters\n\n    def fit(self, X):\n        \"\"\"Fit the clustering model to data X.\"\"\"\n        n_samples = X.shape[0]\n        if self.n_clusters == 1:\n            self.labels_ = np.zeros(n_samples, dtype=int)\n            return self\n\n        # Dimensionality reduction for high-dimensional data\n        X_processed = self._reduce_dimensionality(X)\n\n        # Local density\n        densities = self._compute_local_density(X_processed)\n\n        # Build mutual k-NN affinity matrix\n        affinity = self._build_mutual_knn_graph(X_processed, densities)\n        self.affinity_matrix_ = affinity\n\n        # Determine number of clusters\n        if self.n_clusters is None:\n            n_clusters = self._determine_n_clusters(affinity)\n            n_clusters = max(2, min(n_clusters, n_samples // 2))\n        else:\n            n_clusters = self.n_clusters\n\n        # Spectral embedding\n        embedding = self._spectral_embedding(affinity, n_clusters)\n\n        # Cluster embeddings using k-means\n        kmeans = KMeans(n_clusters=n_clusters, random_state=self.random_state, n_init=10)\n        labels = kmeans.fit_predict(embedding)\n        self.labels_ = labels\n        return self\n\n    def predict(self, X):\n        \"\"\"Predict cluster labels for X (refit for simplicity).\"\"\"\n        if self.labels_ is None:\n            self.fit(X)\n        # Out-of-sample prediction not implemented; refit\n        self.fit(X)\n        return self.labels_\n\n\ndef custom_distance(x, y):\n    \"\"\"Custom distance metric: Euclidean distance.\"\"\"\n    return np.sqrt(np.sum((x - y) ** 2))"
7}

1{
2  "op": "replace",
3  "filename": "scikit-learn/custom_clustering.py",
4  "start_line": 36,
5  "end_line": 247,
6  "content": "import numpy as np\nfrom sklearn.base import BaseEstimator, ClusterMixin\nfrom sklearn.cluster import KMeans, HDBSCAN, SpectralClustering\nfrom sklearn.neighbors import NearestNeighbors\nfrom sklearn.decomposition import PCA\nfrom sklearn.metrics import silhouette_score\nfrom scipy.cluster.hierarchy import linkage, fcluster\nfrom scipy.spatial.distance import squareform\n\n\nclass CustomClustering(BaseEstimator, ClusterMixin):\n    \"\"\"Ensemble Clustering with Adaptive Base Algorithms (ECABA).\n\n    Novel ensemble method that combines K‑means, HDBSCAN, and spectral clustering\n    via co‑association matrix fusion. Each base algorithm captures different\n    geometric structures; the ensemble integrates their strengths.\n\n    Steps:\n    1. If n_features > 50, apply PCA to retain 95% variance.\n    2. Generate three base clusterings:\n        - K‑means with given n_clusters (or default 8).\n        - HDBSCAN (automatic cluster count, density‑based).\n        - Spectral clustering with mutual k‑NN affinity.\n    3. Build co‑association matrix: entry (i,j) = fraction of base clusterings\n       that assign i and j to the same cluster.\n    4. Perform hierarchical clustering (average linkage) on the co‑association\n       similarity matrix.\n    5. Cut dendrogram at n_clusters (if provided) or choose optimal number\n       via silhouette maximization.\n\n    Args:\n        n_clusters: Number of clusters. If None, determined automatically.\n        random_state: Random seed for reproducibility.\n        k_neighbors: Number of neighbors for spectral clustering (default 10).\n    \"\"\"\n\n    def __init__(self, n_clusters=None, random_state=42, k_neighbors=10):\n        self.n_clusters = n_clusters\n        self.random_state = random_state\n        self.k_neighbors = k_neighbors\n        self.labels_ = None\n        self.pca_ = None\n\n    def _reduce_dimensionality(self, X):\n        \"\"\"Apply PCA if n_features > 50.\"\"\"\n        n_samples, n_features = X.shape\n        if n_features <= 50:\n            return X\n        pca = PCA(n_components=0.95, random_state=self.random_state)\n        X_reduced = pca.fit_transform(X)\n        self.pca_ = pca\n        return X_reduced\n\n    def _base_kmeans(self, X, n_clusters):\n        \"\"\"K‑means clustering.\"\"\"\n        km = KMeans(n_clusters=n_clusters, random_state=self.random_state, n_init=10)\n        return km.fit_predict(X)\n\n    def _base_hdbscan(self, X):\n        \"\"\"HDBSCAN density‑based clustering.\"\"\"\n        min_cluster_size = max(5, X.shape[0] // 50)\n        hdb = HDBSCAN(\n            min_cluster_size=min_cluster_size,\n            min_samples=5,\n            cluster_selection_method=\"eom\",\n            gen_min_span_tree=False,\n        )\n        labels = hdb.fit_predict(X)\n        # Convert noise (-1) to unique cluster IDs to avoid artificial agreement\n        # Map each noise point to a distinct negative number\n        noise_mask = labels == -1\n        if noise_mask.any():\n            # Assign each noise point a unique label below -1\n            unique_noise_labels = -np.arange(2, np.sum(noise_mask) + 2)\n            labels[noise_mask] = unique_noise_labels\n        return labels\n\n    def _base_spectral(self, X, n_clusters):\n        \"\"\"Spectral clustering with mutual k‑NN affinity.\"\"\"\n        # Build mutual k‑NN graph\n        n_samples = X.shape[0]\n        k = min(self.k_neighbors, n_samples - 1)\n        if k <= 0:\n            # Fallback to standard spectral clustering\n            sc = SpectralClustering(\n                n_clusters=n_clusters,\n                affinity=\"nearest_neighbors\",\n                random_state=self.random_state,\n            )\n            return sc.fit_predict(X)\n\n        nn = NearestNeighbors(n_neighbors=k + 1)\n        nn.fit(X)\n        indices = nn.kneighbors(X, return_distance=False)[:, 1:]  # exclude self\n\n        # Mutual adjacency\n        adjacency = np.zeros((n_samples, n_samples), dtype=bool)\n        for i in range(n_samples):\n            for j in indices[i]:\n                if i in indices[j]:\n                    adjacency[i, j] = True\n                    adjacency[j, i] = True\n\n        # Convert to sparse matrix for spectral clustering\n        from scipy.sparse import csr_matrix\n        affinity = csr_matrix(adjacency, dtype=float)\n        # Spectral clustering with precomputed affinity\n        sc = SpectralClustering(\n            n_clusters=n_clusters,\n            affinity=\"precomputed\",\n            random_state=self.random_state,\n        )\n        return sc.fit_predict(affinity)\n\n    def _coassociation_matrix(self, base_labels):\n        \"\"\"Compute co‑association matrix from base clusterings.\"\"\"\n        n_samples = base_labels.shape[1]\n        n_base = base_labels.shape[0]\n        coassoc = np.zeros((n_samples, n_samples))\n        # For each base clustering\n        for b in range(n_base):\n            labels = base_labels[b]\n            # Create indicator matrix where entry (i,j) = 1 if same cluster\n            # Vectorized using broadcasting\n            same = labels[:, np.newaxis] == labels[np.newaxis, :]\n            coassoc += same.astype(float)\n        coassoc /= n_base\n        return coassoc\n\n    def _determine_n_clusters(self, coassoc):\n        \"\"\"Determine optimal number of clusters via silhouette maximization.\"\"\"\n        n_samples = coassoc.shape[0]\n        # Convert similarity to distance\n        distance = 1.0 - coassoc\n        np.fill_diagonal(distance, 0.0)\n        # Try candidate numbers of clusters\n        min_c = 2\n        max_c = min(int(np.sqrt(n_samples)) + 1, 20)\n        best_score = -1\n        best_n = min_c\n        # Use hierarchical clustering with average linkage\n        condensed = squareform(distance, checks=False)\n        Z = linkage(condensed, method='average')\n        for k in range(min_c, max_c + 1):\n            labels = fcluster(Z, k, criterion='maxclust')\n            # Adjust labels to start from 0\n            labels = labels - 1\n            if len(np.unique(labels)) < 2:\n                continue\n            score = silhouette_score(distance, labels, metric='precomputed')\n            if score > best_score:\n                best_score = score\n                best_n = k\n        return best_n\n\n    def fit(self, X):\n        \"\"\"Fit the ensemble clustering model.\"\"\"\n        n_samples = X.shape[0]\n        if self.n_clusters == 1:\n            self.labels_ = np.zeros(n_samples, dtype=int)\n            return self\n\n        # Dimensionality reduction for high‑dimensional data\n        X_processed = self._reduce_dimensionality(X)\n\n        # Determine n_clusters for base algorithms that need it\n        if self.n_clusters is None:\n            # Use default guess: sqrt(n_samples) capped at 20\n            k_base = min(int(np.sqrt(n_samples)) + 1, 20)\n        else:\n            k_base = self.n_clusters\n\n        # Generate base clusterings\n        base_labels = []\n        # 1. K‑means\n        base_labels.append(self._base_kmeans(X_processed, k_base))\n        # 2. HDBSCAN\n        base_labels.append(self._base_hdbscan(X_processed))\n        # 3. Spectral clustering\n        base_labels.append(self._base_spectral(X_processed, k_base))\n\n        base_labels = np.array(base_labels)  # shape (3, n_samples)\n\n        # Co‑association matrix\n        coassoc = self._coassociation_matrix(base_labels)\n\n        # Determine final number of clusters\n        if self.n_clusters is None:\n            n_clusters = self._determine_n_clusters(coassoc)\n        else:\n            n_clusters = self.n_clusters\n\n        # Hierarchical clustering on co‑association similarity\n        distance = 1.0 - coassoc\n        np.fill_diagonal(distance, 0.0)\n        condensed = squareform(distance, checks=False)\n        Z = linkage(condensed, method='average')\n        labels = fcluster(Z, n_clusters, criterion='maxclust') - 1\n        self.labels_ = labels\n        return self\n\n    def predict(self, X):\n        \"\"\"Predict cluster labels for X (refit for simplicity).\"\"\"\n        if self.labels_ is None:\n            self.fit(X)\n        # Out‑of‑sample prediction not implemented; refit\n        self.fit(X)\n        return self.labels_\n\n\ndef custom_distance(x, y):\n    \"\"\"Custom distance metric: Euclidean distance.\"\"\"\n    return np.sqrt(np.sum((x - y) ** 2))"
7}

1{
2  "op": "replace",
3  "filename": "scikit-learn/custom_clustering.py",
4  "start_line": 36,
5  "end_line": 248,
6  "content": "import numpy as np\nfrom sklearn.base import BaseEstimator, ClusterMixin\nfrom sklearn.cluster import KMeans, HDBSCAN\nfrom sklearn.decomposition import PCA\nfrom sklearn.metrics import silhouette_score\n\n\nclass CustomClustering(BaseEstimator, ClusterMixin):\n    \"\"\"Silhouette‑Optimized Hybrid Clustering (SOHC).\n\n    Novel algorithm that dynamically chooses between density‑based (HDBSCAN)\n    and partition‑based (K‑means) clustering based on the silhouette score\n    of each method. This adapts to the underlying geometry of the data.\n\n    Steps:\n    1. If n_features > 50, apply PCA to retain 95% variance.\n    2. Run HDBSCAN (with automatic cluster count) and K‑means (with given\n       n_clusters or default sqrt(n_samples)).\n    3. Compute silhouette scores for both labelings (using Euclidean distance).\n    4. Select the labeling with the higher silhouette score.\n    5. If the selected labeling has only one cluster (or all noise), fall back\n       to the other labeling.\n\n    Args:\n        n_clusters: Number of clusters for K‑means. If None, a default guess\n                    is used (sqrt(n_samples) capped at 20).\n        random_state: Random seed for reproducibility.\n    \"\"\"\n\n    def __init__(self, n_clusters=None, random_state=42):\n        self.n_clusters = n_clusters\n        self.random_state = random_state\n        self.labels_ = None\n        self.pca_ = None\n        self.chosen_method_ = None\n\n    def _reduce_dimensionality(self, X):\n        \"\"\"Apply PCA if n_features > 50.\"\"\"\n        n_samples, n_features = X.shape\n        if n_features <= 50:\n            return X\n        pca = PCA(n_components=0.95, random_state=self.random_state)\n        X_reduced = pca.fit_transform(X)\n        self.pca_ = pca\n        return X_reduced\n\n    def fit(self, X):\n        \"\"\"Fit the hybrid clustering model.\"\"\"\n        n_samples = X.shape[0]\n        if self.n_clusters == 1:\n            self.labels_ = np.zeros(n_samples, dtype=int)\n            self.chosen_method_ = 'trivial'\n            return self\n\n        # Dimensionality reduction for high‑dimensional data\n        X_processed = self._reduce_dimensionality(X)\n\n        # Determine n_clusters for K‑means\n        if self.n_clusters is None:\n            k = min(int(np.sqrt(n_samples)) + 1, 20)\n        else:\n            k = self.n_clusters\n\n        # Run HDBSCAN\n        min_cluster_size = max(5, n_samples // 50)\n        hdb = HDBSCAN(\n            min_cluster_size=min_cluster_size,\n            min_samples=5,\n            cluster_selection_method=\"eom\",\n            gen_min_span_tree=False,\n        )\n        labels_hdb = hdb.fit_predict(X_processed)\n        # Convert noise (-1) to unique negative labels for silhouette computation\n        noise_mask = labels_hdb == -1\n        if noise_mask.any():\n            unique_noise = -np.arange(2, np.sum(noise_mask) + 2)\n            labels_hdb[noise_mask] = unique_noise\n\n        # Run K‑means\n        km = KMeans(n_clusters=k, random_state=self.random_state, n_init=10)\n        labels_km = km.fit_predict(X_processed)\n\n        # Compute silhouette scores (skip if only one cluster or all noise)\n        def valid_silhouette(labels):\n            unique = np.unique(labels)\n            if len(unique) <= 1:\n                return -1.0\n            # If any label is negative (noise), treat as separate clusters for silhouette\n            # but silhouette expects positive labels; we can map to positive\n            if np.any(labels < 0):\n                # Map each distinct label to a positive integer\n                _, mapped = np.unique(labels, return_inverse=True)\n                labels = mapped\n            # Silhouette with Euclidean distance\n            try:\n                return silhouette_score(X_processed, labels, metric='euclidean')\n            except:\n                return -1.0\n\n        score_hdb = valid_silhouette(labels_hdb)\n        score_km = valid_silhouette(labels_km)\n\n        # Choose the labeling with higher silhouette\n        if score_hdb >= score_km:\n            chosen_labels = labels_hdb\n            self.chosen_method_ = 'HDBSCAN'\n        else:\n            chosen_labels = labels_km\n            self.chosen_method_ = 'KMeans'\n\n        # If chosen labeling has only one effective cluster (or all noise),\n        # fall back to the other labeling.\n        unique_chosen = np.unique(chosen_labels)\n        if len(unique_chosen) <= 1 or (np.all(chosen_labels < 0)):\n            if self.chosen_method_ == 'HDBSCAN':\n                chosen_labels = labels_km\n                self.chosen_method_ = 'KMeans (fallback)'\n            else:\n                chosen_labels = labels_hdb\n                self.chosen_method_ = 'HDBSCAN (fallback)'\n\n        # Ensure labels are non‑negative consecutive integers starting from 0\n        _, mapped = np.unique(chosen_labels, return_inverse=True)\n        self.labels_ = mapped\n        return self\n\n    def predict(self, X):\n        \"\"\"Predict cluster labels for X (refit for simplicity).\"\"\"\n        if self.labels_ is None:\n            self.fit(X)\n        # Out‑of‑sample prediction not implemented; refit\n        self.fit(X)\n        return self.labels_\n\n\ndef custom_distance(x, y):\n    \"\"\"Custom distance metric: Euclidean distance.\"\"\"\n    return np.sqrt(np.sum((x - y) ** 2))"
7}

1{
2  "n": 2
3}