""" gen_data.py ================ Parameters ---------- n_samples : total number of points (inliers + outliers) n_features : number of features n_clusters : number of clusters outlier_fraction : fraction of points that are outliers [0, 0.5) difficulty : float in [0, 1]. 0 = very easy, 1 = very hard. Controls: cluster separation, density heterogeneity, non-convexity, anisotropy, imbalance, and outlier ambiguity. random_state : int or None return_params : if True, also return the dict of sampled parameters Internal difficulty knobs (can be overridden via kwargs) -------------------------------------------------------- min_sep_factor : minimum inter-center distance as a multiple of the mean cluster sigma. Interpolated [5 → 0.4] with difficulty. imbalance_alpha : Dirichlet concentration. [100 -> 0.3] with difficulty. density_ratio : max/min sigma ratio. [1 -> 8] with difficulty. nonconvex_prob : probability each cluster is non-convex. [0 -> 1]. rotate : whether to apply a random rotation to up to 5 features (n_features>=2). inl_spread : base inlier spread per cluster. Default 0.3. outl_spread_factor: fraction of the inlier range used as an additional boundary margin for global outliers. Default 0.5. collective_frac : fraction of outliers placed as dense groups. Default 0.2. """ from __future__ import annotations import warnings import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from sklearn.decomposition import PCA import os from sklearn.preprocessing import RobustScaler def save_dataset_scatter(X, y_clust, y_bin, out_path): pca = PCA(n_components=2, random_state=0) Z = pca.fit_transform(X) idx = np.random.permutation(1000) Z,y_clust,y_bin = Z[idx,:],y_clust[idx],y_bin[idx] fig, ax = plt.subplots(figsize=(3, 3)) inl = y_bin == 0 ax.scatter(Z[inl, 0], Z[inl, 1], c=y_clust[inl], cmap="tab10", s=8, alpha=0.5, linewidths=0, label="Inliers") out = y_bin == 1 ax.scatter(Z[out, 0], Z[out, 1], c="red", marker="x", s=20, alpha=0.7, linewidths=0.8, label="Outliers") ax.set_xticks([]) ax.set_yticks([]) ax.legend(loc="upper right") ax.set_title("2D PCA (1000 points)", fontsize=12) plt.tight_layout() plt.savefig(out_path, dpi=200, bbox_inches="tight") plt.close() print(f"Plot saved to {out_path}") def _smooth_fold(X, a=1.5, b=0.75): """Sinusoidal fold in 2-D subspace.""" X[:, 1] += a * np.sin(X[:, 0]) X[:, 0] += b * np.sin(X[:, 1]) return X def _smooth_holes(X, k=3, strength=0.4): """Radial oscillation (ring-like distortion) in 2-D subspace.""" r = np.sqrt(X[:, 0] ** 2 + X[:, 1] ** 2) / np.sqrt(2) X[:, 0] += strength * np.cos(k * r) X[:, 1] += strength * np.sin(k * r) return X def _apply_nonconvex(cluster, rng): """ Apply non-convex transforms to multiple random pairs of dimensions. For n_features == 1 the transform is skipped silently. """ n, d = cluster.shape if d < 2: return cluster # Number of pairs to transform: 1 for d<4, up to 3 for d>=8 n_pairs = min(3, max(1, d // 4)) dims = rng.choice(d, size=n_pairs * 2, replace=False) for i in range(n_pairs): i0, i1 = int(dims[2 * i]), int(dims[2 * i + 1]) sub = cluster[:, [i0, i1]].copy() shape = rng.choice(["fold", "holes", "both"]) a = rng.uniform(1.2, 2.0) b = rng.uniform(0.3, 0.8) k = rng.integers(2, 5) st = rng.uniform(0.2, 0.50) if shape in ("fold", "both"): sub = _smooth_fold(sub, a=a, b=b) if shape in ("holes", "both"): sub = _smooth_holes(sub, k=k, strength=st) cluster[:, [i0, i1]] = sub return cluster #def _random_rotation(d, rng): # """Random orthogonal matrix via QR decomposition.""" # A = rng.standard_normal((d, d)) # Q, R = np.linalg.qr(A) # Q *= np.sign(np.diag(R)) # make deterministic sign # return Q.astype(np.float32) def _place_centers(n_clusters, n_features, mean_sigma, min_sep_factor, rng, max_attempts=10000): """ Sample cluster centers with a minimum pairwise distance of `min_sep_factor * mean_sigma`. Falls back to random placement if the constraint cannot be satisfied after `max_attempts`. """ min_dist = min_sep_factor * mean_sigma # Adapt bounding box so clusters don't all crowd the origin scale = max(1.0, np.log10(max(n_clusters, 2))) * max(1.0, n_clusters ** (1.0 / n_features) ) lo, hi = -scale * 3, scale * 3 centers: list[np.ndarray] = [] attempts = 0 while len(centers) < n_clusters and attempts < max_attempts: c = rng.uniform(lo, hi, size=n_features) if all(np.linalg.norm(c - existing) >= min_dist for existing in centers): centers.append(c) attempts += 1 if len(centers) < n_clusters: # Constraint unsatisfiable → fall back to unconstrained warnings.warn(f"Could not place {n_clusters} centers with min_sep_factor={min_sep_factor:.2f}; " "falling back to unconstrained placement.", RuntimeWarning, stacklevel=3 ) extra = rng.uniform(lo, hi, size=(n_clusters - len(centers), n_features)) centers.extend(extra) return np.array(centers, dtype=np.float32) def gen_data(n_samples, n_features, n_clusters, outlier_fraction, difficulty=0.5, *, min_sep_factor=None, imbalance_alpha=None, density_ratio=None, nonconvex_prob=None, rotate=True, inl_spread=0.3, outl_spread_factor=0.5, collective_frac=0.2, random_state=None, return_params=False, control_plot=False): rng = np.random.default_rng(random_state) d = difficulty # ── difficulty-dependent parameters ────────────────────────────── min_sep = min_sep_factor if min_sep_factor is not None else np.interp(d,[0,1],[5.0,0.4]) alpha = imbalance_alpha if imbalance_alpha is not None else np.interp(d,[0,1],[100.0,0.3]) dratio = density_ratio if density_ratio is not None else np.interp(d,[0,1],[1.0,8.0]) nc_prob = nonconvex_prob if nonconvex_prob is not None else np.interp(d,[0,1],[0.0,1.0]) # ── counts ─────────────────────────────────────────────────────── n_in = int(n_samples * (1 - outlier_fraction)) n_out = n_samples - n_in # ── densities ──────────────────────────────────────────────────── if dratio <= 1: sigmas = np.full(n_clusters, inl_spread, dtype=np.float32) else: logh = 0.5*np.log(dratio) sigmas = (inl_spread * np.exp(rng.uniform(-logh, logh, n_clusters))).astype(np.float32) mean_sigma = sigmas.mean() # ── centers ────────────────────────────────────────────────────── centers = _place_centers(n_clusters, n_features, mean_sigma, min_sep, rng).astype(np.float32) # ── sizes ──────────────────────────────────────────────────────── w = rng.dirichlet(np.ones(n_clusters)*alpha) counts = np.floor(w*n_in).astype(int); counts[-1] += n_in - counts.sum() # ── inliers ────────────────────────────────────────────────────── X_in, y_in = [], [] for i in range(n_clusters): n = counts[i] if n == 0: continue center, sigma = centers[i], sigmas[i] # Higher difficulty increases per-feature spread variability. anisotropy = 1 + 4 * d # difficulty-driven base = rng.uniform(0.2, anisotropy, n_features).astype(np.float32) scales = base * sigma cluster = rng.standard_normal((n,n_features), dtype=np.float32) cluster = cluster * scales + center if rng.random()=2: cluster = _apply_nonconvex(cluster, rng) X_in.append(cluster) y_in.append(np.full(n,i,np.int32)) X_in = np.vstack(X_in) y_in = np.concatenate(y_in) # ── OUTLIERS: global, local, collective, and extreme ───────────── if n_out == 0: X_out = np.empty((0,n_features),np.float32) else: lo, hi = X_in.min(0), X_in.max(0) margin = outl_spread_factor*(hi-lo) lo, hi = lo-margin, hi+margin f_global, f_local, f_collective, f_extreme = 0.25, 0.35, collective_frac, 0.15 n_global = int(n_out*f_global) n_local = int(n_out*f_local) n_coll = int(n_out*f_collective) n_ext = n_out - (n_global+n_local+n_coll) parts = [] if n_global>0: parts.append(rng.uniform(lo,hi,(n_global,n_features)).astype(np.float32)) if n_local>0: idx = rng.choice(len(X_in), n_local, replace=False) pts = X_in[idx] tgt = centers[rng.integers(len(centers), size=n_local)] dirs = tgt - pts noise = rng.normal(0,0.2,pts.shape).astype(np.float32) parts.append(pts + 0.5*dirs + noise) created = 0 while created < n_coll: g = min(rng.integers(5,30), n_coll-created) base = centers[rng.integers(len(centers))] shift = rng.normal(0,3*mean_sigma,n_features).astype(np.float32) c = base + shift parts.append(rng.normal(c, mean_sigma*0.1, (g,n_features)).astype(np.float32)) created += g if n_ext>0: Xtmp = rng.normal(0,1,(n_ext,n_features)).astype(np.float32) scale = rng.uniform(6,10) parts.append(Xtmp * scale * mean_sigma) X_out = np.vstack(parts)[:n_out].astype(np.float32) # ── merge ──────────────────────────────────────────────────────── X = np.vstack([X_in,X_out]).astype(np.float32) yc = np.concatenate([y_in, -1*np.ones(len(X_out),np.int32)]) yb = (yc==-1).astype(np.int32) perm = rng.permutation(len(X)) X, yc, yb = X[perm], yc[perm], yb[perm] if rotate and n_features >= 2: k = min(5, n_features) idx = rng.choice(n_features, size=k, replace=False) A = rng.standard_normal((k, k)).astype(np.float32) Q, _ = np.linalg.qr(A) Q = Q.astype(np.float32) X[:, idx] = (X[:, idx] @ Q).astype(np.float32) # control plot if control_plot: img_dir = "control_plots" os.makedirs(img_dir, exist_ok=True) name = f"{n_samples}_{n_features}_{n_clusters}_{outlier_fraction}_{difficulty}_seed{random_state}".replace(".", "-") img_path = f"{img_dir}/cp_{name}.pdf" save_dataset_scatter(X, yc, yb, img_path) scaler = RobustScaler() X_ins = scaler.fit_transform(X_in.astype("float32")) if return_params: return X_ins, y_in, X, yb, yc, dict( difficulty=d, min_sep=min_sep, alpha=alpha, dratio=dratio, nc_prob=nc_prob ) return X_ins, y_in, X, yb, yc if __name__ == "__main__": X_in, y_in, Xo, yb, yc = gen_data( n_samples=2000, n_features=10, n_clusters=5, outlier_fraction=0.1, difficulty=0.2, random_state=42, control_plot = True) X_in, y_in, Xo, yb, yc = gen_data( n_samples=2000, n_features=10, n_clusters=5, outlier_fraction=0.1, difficulty=0.8, random_state=42, control_plot = True)