1{
2  "op": "replace",
3  "filename": "scikit-learn/bench/custom_dimred.py",
4  "start_line": 15,
5  "end_line": 59,
6  "content": "class CustomDimReduction:\n    \"\"\"Triplet-based Multiscale Embedding (TriME).\n\n    Combines multiple neighborhood scales to preserve both local and global\n    structure. For each point, we consider neighbors at different radii,\n    and optimize a hinge loss that enforces relative distance ordering.\n    \"\"\"\n\n    def __init__(self, n_components: int = 2, random_state: int | None = None):\n        self.n_components = n_components\n        self.random_state = random_state\n        # Hyperparameters\n        self.n_neighbors_list = [5, 15, 50]  # multiple scales\n        self.n_negatives = 5  # number of negative samples per positive\n        self.lr = 0.1  # learning rate\n        self.n_epochs = 200\n        self.momentum = 0.9\n        self.hinge_margin = 1.0\n\n    def _compute_neighbors(self, X, n_neighbors):\n        \"\"\"Compute k-nearest neighbors using Euclidean distance.\"\"\"\n        from scipy.spatial.distance import cdist\n        D = cdist(X, X, metric='euclidean')\n        n = X.shape[0]\n        # For each point, get indices of nearest neighbors (excluding self)\n        neighbors = np.zeros((n, n_neighbors), dtype=int)\n        for i in range(n):\n            idx = np.argsort(D[i])[1:n_neighbors+1]  # skip self\n            neighbors[i] = idx\n        return neighbors, D\n\n    def _sample_triplets(self, X, neighbors_list, D):\n        \"\"\"Generate triplets (i, j, k) for each scale.\n        \n        Returns list of triplets per scale, each as (i, j, k) arrays.\n        \"\"\"\n        n = X.shape[0]\n        triplets = []\n        for scale_idx, neighbors in enumerate(neighbors_list):\n            n_neighbors = neighbors.shape[1]\n            # For each point i, each neighbor j, sample negatives k\n            scale_triplets = []\n            for i in range(n):\n                for j_idx in range(n_neighbors):\n                    j = neighbors[i, j_idx]\n                    # Sample negatives that are not neighbors at this scale\n                    # We'll consider points not in the neighbor set\n                    possible_negatives = [k for k in range(n) \n                                          if k != i and k not in neighbors[i]]\n                    if len(possible_negatives) == 0:\n                        continue\n                    # Sample self.n_negatives negatives\n                    negs = np.random.choice(possible_negatives, \n                                            size=min(self.n_negatives, len(possible_negatives)), \n                                            replace=False)\n                    for k in negs:\n                        scale_triplets.append((i, j, k))\n            # Convert to arrays for efficiency\n            if scale_triplets:\n                i_arr, j_arr, k_arr = zip(*scale_triplets)\n                triplets.append((np.array(i_arr), np.array(j_arr), np.array(k_arr)))\n        return triplets\n\n    def _loss_gradient(self, Y, triplets, scale_weights):\n        \"\"\"Compute loss and gradient for current embedding Y.\"\"\"\n        n, dim = Y.shape\n        grad = np.zeros_like(Y)\n        loss = 0.0\n        for scale_idx, (i_arr, j_arr, k_arr) in enumerate(triplets):\n            w = scale_weights[scale_idx]\n            # Compute squared distances in embedding\n            D_ij = np.sum((Y[i_arr] - Y[j_arr]) ** 2, axis=1)\n            D_ik = np.sum((Y[i_arr] - Y[k_arr]) ** 2, axis=1)\n            # Hinge loss: max(0, margin + D_ij - D_ik)\n            violations = self.hinge_margin + D_ij - D_ik\n            mask = violations > 0\n            if np.any(mask):\n                loss += w * np.sum(violations[mask])\n                # Gradients\n                # dL/dY_i = 2 * w * (Y_i - Y_j) - 2 * w * (Y_i - Y_k) = 2w(Y_k - Y_j)\n                # dL/dY_j = -2w(Y_i - Y_j)\n                # dL/dY_k = 2w(Y_i - Y_k)\n                factor = 2 * w\n                # For each triplet with violation\n                for idx in np.where(mask)[0]:\n                    i = i_arr[idx]\n                    j = j_arr[idx]\n                    k = k_arr[idx]\n                    grad[i] += factor * (Y[k] - Y[j])\n                    grad[j] += factor * (Y[i] - Y[j])\n                    grad[k] += factor * (Y[i] - Y[k])\n        return loss, grad\n\n    def fit_transform(self, X: NDArray[np.float64]) -> NDArray[np.float64]:\n        \"\"\"Reduce dimensionality of X.\"\"\"\n        rng = np.random.RandomState(self.random_state)\n        np.random.seed(self.random_state if self.random_state is not None else 42)\n        \n        n_samples, n_features = X.shape\n        n_components = self.n_components\n        \n        # Step 1: Compute neighbors at multiple scales\n        neighbors_list = []\n        dists_list = []\n        for k in self.n_neighbors_list:\n            neighbors, D = self._compute_neighbors(X, k)\n            neighbors_list.append(neighbors)\n            dists_list.append(D)\n        \n        # Step 2: Generate triplets for each scale\n        triplets = self._sample_triplets(X, neighbors_list, dists_list[0])\n        \n        # Step 3: Initialize embedding using PCA\n        from sklearn.decomposition import PCA\n        pca = PCA(n_components=n_components, random_state=self.random_state)\n        Y = pca.fit_transform(X).astype(np.float64)\n        \n        # Step 4: Scale weights inversely with neighborhood size (more weight to local)\n        scale_weights = [1.0 / (k ** 0.5) for k in self.n_neighbors_list]\n        # Normalize\n        total = sum(scale_weights)\n        scale_weights = [w / total for w in scale_weights]\n        \n        # Step 5: Stochastic gradient descent with momentum\n        velocity = np.zeros_like(Y)\n        for epoch in range(self.n_epochs):\n            loss, grad = self._loss_gradient(Y, triplets, scale_weights)\n            # Update with momentum\n            velocity = self.momentum * velocity - self.lr * grad / n_samples\n            Y += velocity\n            # Optional: early stopping\n            if epoch % 50 == 0:\n                # Perturb to escape saddle points\n                Y += rng.randn(*Y.shape) * 0.01\n        \n        return Y"
7}

1{
2  "op": "replace",
3  "filename": "scikit-learn/bench/custom_dimred.py",
4  "start_line": 15,
5  "end_line": 150,
6  "content": "class CustomDimReduction:\n    \"\"\"Multiscale Laplacian Eigenmaps with adaptive bandwidth (MLEM).\n\n    Combines multiple Gaussian affinity scales using per-point bandwidths\n    derived from local neighborhoods. The combined affinity matrix is used\n    to compute a spectral embedding that preserves both local and global\n    structure.\n    \"\"\"\n\n    def __init__(self, n_components: int = 2, random_state: int | None = None):\n        self.n_components = n_components\n        self.random_state = random_state\n        # Hyperparameters\n        self.n_neighbors = 15  # for local bandwidth estimation\n        self.scale_factors = [0.5, 1.0, 2.0, 4.0]  # multipliers for bandwidth\n        self.alpha_power = 1.0  # weight for each scale: alpha_s = 1/(scale_factor^alpha_power)\n        self.knn_affinity = 50  # number of neighbors to keep for sparse affinity\n        self.eigen_tol = 1e-4  # tolerance for eigensolver\n\n    def _compute_pairwise_distances(self, X):\n        \"\"\"Compute squared Euclidean distances between all pairs.\"\"\"\n        from scipy.spatial.distance import cdist\n        D = cdist(X, X, metric='sqeuclidean')\n        return D\n\n    def _local_bandwidths(self, D):\n        \"\"\"Compute per-point bandwidth as median distance to k nearest neighbors.\"\"\"\n        n = D.shape[0]\n        k = min(self.n_neighbors, n - 1)\n        # For each point, get distances to k nearest neighbors (excluding self)\n        neighbor_dists = np.zeros((n, k))\n        for i in range(n):\n            dists = D[i]\n            # skip self\n            idx = np.argsort(dists)[1:k+1]\n            neighbor_dists[i] = dists[idx]\n        sigma_sq = np.median(neighbor_dists, axis=1)  # median squared distance\n        # Ensure non-zero\n        sigma_sq = np.maximum(sigma_sq, 1e-8)\n        return sigma_sq\n\n    def _sparse_affinity_matrix(self, D, sigma_sq_i, sigma_sq_j, scale_factor):\n        \"\"\"Compute sparse affinity matrix for a given scale factor.\n        \n        Returns a sparse matrix W where W[i,j] = exp(-D[i,j] / (sigma_i * sigma_j * scale_factor)).\n        Only keep top knn_affinity neighbors per point.\n        \"\"\"\n        from scipy.sparse import csr_matrix, lil_matrix\n        n = D.shape[0]\n        # Compute combined bandwidth matrix\n        # sigma_i * sigma_j = sqrt(sigma_sq_i * sigma_sq_j)\n        sigma_prod = np.sqrt(sigma_sq_i[:, None] * sigma_sq_j[None, :])\n        # Scale by factor\n        scale = scale_factor ** 2  # because D is squared Euclidean\n        # Affinity = exp(-D / (sigma_prod * scale))\n        # For numerical stability, compute only for nearest neighbors\n        W = lil_matrix((n, n), dtype=np.float64)\n        k_keep = min(self.knn_affinity, n - 1)\n        for i in range(n):\n            # Compute affinities to all points\n            aff = np.exp(-D[i] / (sigma_prod[i] * scale))\n            # Set self-affinity to zero\n            aff[i] = 0\n            # Keep top k_keep neighbors\n            idx = np.argsort(aff)[::-1][:k_keep]\n            for j in idx:\n                if aff[j] > 1e-12:  # threshold\n                    W[i, j] = aff[j]\n        # Make symmetric (average)\n        W = (W + W.T) / 2\n        return W.tocsr()\n\n    def fit_transform(self, X: NDArray[np.float64]) -> NDArray[np.float64]:\n        \"\"\"Reduce dimensionality of X.\"\"\"\n        rng = np.random.RandomState(self.random_state)\n        np.random.seed(self.random_state if self.random_state is not None else 42)\n        \n        n_samples, n_features = X.shape\n        n_components = self.n_components\n        \n        # Step 1: Compute pairwise squared Euclidean distances\n        D = self._compute_pairwise_distances(X)\n        \n        # Step 2: Compute per-point bandwidths (local scale)\n        sigma_sq = self._local_bandwidths(D)\n        \n        # Step 3: Combine affinity matrices across scales\n        from scipy.sparse import csr_matrix\n        W_combined = csr_matrix((n_samples, n_samples), dtype=np.float64)\n        \n        for s, scale_factor in enumerate(self.scale_factors):\n            # Compute affinity for this scale\n            W_s = self._sparse_affinity_matrix(D, sigma_sq, sigma_sq, scale_factor)\n            # Weight for this scale\n            weight = 1.0 / (scale_factor ** self.alpha_power)\n            W_combined += weight * W_s\n        \n        # Step 4: Symmetrize and ensure diagonal zero\n        W_combined = (W_combined + W_combined.T) / 2\n        W_combined.setdiag(0)\n        \n        # Step 5: Compute normalized Laplacian L = I - D^{-1/2} W D^{-1/2}\n        # Actually we want generalized eigenvalue problem (D - W) v = λ D v\n        # Equivalent to solving for smallest eigenvalues of L_sym = I - D^{-1/2} W D^{-1/2}\n        from scipy.sparse import diags\n        d = np.array(W_combined.sum(axis=1)).flatten()\n        # Avoid division by zero\n        d[d == 0] = 1e-12\n        D_sqrt_inv = diags(1.0 / np.sqrt(d))\n        # Compute symmetric Laplacian\n        I = diags(np.ones(n_samples))\n        L_sym = I - D_sqrt_inv @ W_combined @ D_sqrt_inv\n        \n        # Step 6: Compute smallest eigenvectors (excluding constant vector)\n        from scipy.sparse.linalg import eigsh\n        # Compute eigenvalues/vectors for smallest magnitude\n        # We need n_components+1 eigenvectors because first eigenvalue is 0\n        n_eig = n_components + 1\n        try:\n            eigvals, eigvecs = eigsh(L_sym, k=n_eig, which='SM', tol=self.eigen_tol)\n        except:\n            # Fallback: use dense eigensolver\n            L_sym_dense = L_sym.toarray()\n            eigvals, eigvecs = np.linalg.eigh(L_sym_dense)\n            # Sort by eigenvalue ascending\n            idx = np.argsort(eigvals)[:n_eig]\n            eigvals = eigvals[idx]\n            eigvecs = eigvecs[:, idx]\n        \n        # Sort eigenvectors by eigenvalue (ascending)\n        idx = np.argsort(eigvals)\n        eigvecs = eigvecs[:, idx]\n        # Skip the first eigenvector (constant) corresponding to eigenvalue ~0\n        embedding = eigvecs[:, 1:n_components+1].real\n        \n        # Step 7: Orthogonalize using PCA (optional)\n        from sklearn.decomposition import PCA\n        pca = PCA(n_components=n_components, random_state=self.random_state)\n        embedding = pca.fit_transform(embedding)\n        \n        return embedding.astype(np.float64)"
7}

1{
2  "op": "replace",
3  "filename": "scikit-learn/bench/custom_dimred.py",
4  "start_line": 15,
5  "end_line": 155,
6  "content": "class CustomDimReduction:\n    \"\"\"Multiscale Spectral Embedding with Force Refinement (MSE-FR).\n\n    Combines multi-scale affinity spectral embedding with a force-directed\n    refinement that adds repulsion between dissimilar points to improve\n    global structure preservation.\n    \"\"\"\n\n    def __init__(self, n_components: int = 2, random_state: int | None = None):\n        self.n_components = n_components\n        self.random_state = random_state\n        # Hyperparameters for spectral embedding\n        self.n_neighbors = 15\n        self.scale_factors = [0.5, 1.0, 2.0, 4.0]\n        self.knn_affinity = 100  # keep more neighbors for global structure\n        # Hyperparameters for force refinement\n        self.refine_epochs = 100\n        self.lr = 0.1\n        self.attraction_weight = 1.0\n        self.repulsion_weight = 0.1\n        self.repulsion_neighbors = 500  # number of repulsive forces per point\n        self.early_exaggeration = 4.0  # first few epochs exaggeration\n\n    def _compute_pairwise_distances(self, X):\n        from scipy.spatial.distance import cdist\n        D = cdist(X, X, metric='sqeuclidean')\n        return D\n\n    def _local_bandwidths(self, D):\n        n = D.shape[0]\n        k = min(self.n_neighbors, n - 1)\n        neighbor_dists = np.zeros((n, k))\n        for i in range(n):\n            dists = D[i]\n            idx = np.argsort(dists)[1:k+1]\n            neighbor_dists[i] = dists[idx]\n        sigma_sq = np.median(neighbor_dists, axis=1)\n        sigma_sq = np.maximum(sigma_sq, 1e-8)\n        return sigma_sq\n\n    def _sparse_affinity_matrix(self, D, sigma_sq_i, sigma_sq_j, scale_factor):\n        from scipy.sparse import csr_matrix, lil_matrix\n        n = D.shape[0]\n        sigma_prod = np.sqrt(sigma_sq_i[:, None] * sigma_sq_j[None, :])\n        scale = scale_factor ** 2\n        W = lil_matrix((n, n), dtype=np.float64)\n        k_keep = min(self.knn_affinity, n - 1)\n        for i in range(n):\n            aff = np.exp(-D[i] / (sigma_prod[i] * scale))\n            aff[i] = 0\n            idx = np.argsort(aff)[::-1][:k_keep]\n            for j in idx:\n                if aff[j] > 1e-12:\n                    W[i, j] = aff[j]\n        W = (W + W.T) / 2\n        return W.tocsr()\n\n    def _spectral_embedding(self, X):\n        \"\"\"Compute multi-scale spectral embedding.\"\"\"\n        D = self._compute_pairwise_distances(X)\n        sigma_sq = self._local_bandwidths(D)\n        \n        from scipy.sparse import csr_matrix, diags\n        n_samples = X.shape[0]\n        W_combined = csr_matrix((n_samples, n_samples), dtype=np.float64)\n        \n        for scale_factor in self.scale_factors:\n            W_s = self._sparse_affinity_matrix(D, sigma_sq, sigma_sq, scale_factor)\n            weight = 1.0 / scale_factor  # favor finer scales\n            W_combined += weight * W_s\n        \n        W_combined = (W_combined + W_combined.T) / 2\n        W_combined.setdiag(0)\n        \n        d = np.array(W_combined.sum(axis=1)).flatten()\n        d[d == 0] = 1e-12\n        D_sqrt_inv = diags(1.0 / np.sqrt(d))\n        I = diags(np.ones(n_samples))\n        L_sym = I - D_sqrt_inv @ W_combined @ D_sqrt_inv\n        \n        from scipy.sparse.linalg import eigsh\n        n_eig = self.n_components + 1\n        try:\n            eigvals, eigvecs = eigsh(L_sym, k=n_eig, which='SM', tol=1e-4)\n        except:\n            L_sym_dense = L_sym.toarray()\n            eigvals, eigvecs = np.linalg.eigh(L_sym_dense)\n            idx = np.argsort(eigvals)[:n_eig]\n            eigvecs = eigvecs[:, idx]\n        \n        idx = np.argsort(eigvals)\n        eigvecs = eigvecs[:, idx]\n        embedding = eigvecs[:, 1:self.n_components+1].real\n        \n        # Normalize\n        from sklearn.preprocessing import StandardScaler\n        embedding = StandardScaler().fit_transform(embedding)\n        return embedding, W_combined\n\n    def _force_refinement(self, Y, W, D_original):\n        \"\"\"Refine embedding using attractive and repulsive forces.\"\"\"\n        n, dim = Y.shape\n        rng = np.random.RandomState(self.random_state)\n        \n        # Precompute repulsion candidate indices (far points in original space)\n        # For each point, select self.repulsion_neighbors farthest points\n        repulse_indices = []\n        for i in range(n):\n            # Sort by original distance descending\n            far_idx = np.argsort(D_original[i])[::-1][:self.repulsion_neighbors]\n            repulse_indices.append(far_idx)\n        \n        # SGD with momentum\n        velocity = np.zeros_like(Y)\n        momentum = 0.9\n        \n        for epoch in range(self.refine_epochs):\n            exaggeration = self.early_exaggeration if epoch < 20 else 1.0\n            grad = np.zeros_like(Y)\n            \n            # Attractive forces from graph edges\n            # For each edge (i,j) with weight W[i,j]\n            rows, cols = W.nonzero()\n            weights = W.data\n            for idx in range(len(rows)):\n                i = rows[idx]\n                j = cols[idx]\n                w = weights[idx] * exaggeration\n                # Attractive force: gradient of (Y_i - Y_j)^2 * w\n                diff = Y[i] - Y[j]\n                grad[i] -= w * diff\n                grad[j] += w * diff\n            \n            # Repulsive forces from far points\n            for i in range(n):\n                far = repulse_indices[i]\n                # Compute repulsion: push away from far points\n                for j in far:\n                    if i == j:\n                        continue\n                    diff = Y[i] - Y[j]\n                    dist_sq = np.sum(diff**2) + 1e-8\n                    # Inverse square law repulsion\n                    strength = self.repulsion_weight / dist_sq\n                    grad[i] += strength * diff\n                    grad[j] -= strength * diff\n            \n            # Update with momentum\n            velocity = momentum * velocity - self.lr * grad / n\n            Y += velocity\n            \n            # Optional: early stopping criterion\n            if epoch % 20 == 0:\n                # Add small noise to escape local minima\n                Y += rng.randn(*Y.shape) * 0.01\n        \n        return Y\n\n    def fit_transform(self, X: NDArray[np.float64]) -> NDArray[np.float64]:\n        \"\"\"Reduce dimensionality of X.\"\"\"\n        rng = np.random.RandomState(self.random_state)\n        np.random.seed(self.random_state if self.random_state is not None else 42)\n        \n        # Step 1: Compute spectral embedding\n        Y_spectral, W = self._spectral_embedding(X)\n        \n        # Step 2: Compute original distances for repulsion\n        D_original = self._compute_pairwise_distances(X)\n        \n        # Step 3: Force refinement\n        Y = self._force_refinement(Y_spectral, W, D_original)\n        \n        # Final PCA orthogonalization\n        from sklearn.decomposition import PCA\n        pca = PCA(n_components=self.n_components, random_state=self.random_state)\n        Y = pca.fit_transform(Y)\n        \n        return Y.astype(np.float64)"
7}

1{
2  "n": 2
3}