"""Core GICDM implementation reproducing the paper arXiv:2602.16449. Implements: - ICDM iterative hubness reduction (Section 3 / Algorithm in paper lines 485-498) - GICDM out-of-sample generated-point scaling (Eq. 1, Algorithm 1) - Hubness statistics h5_1(1%) and A5 (Table 4) - Clipped Density / Clipped Coverage fidelity & coverage metrics - Crossover dimension d* (Proposition 5.3) - Raisa et al. (2025) synthetic benchmark scenarios """ import numpy as np from scipy import stats from scipy.special import ncfdtr, ncfdtri # ---------------------------------------------------------------------------- # Distances & k-nearest-neighbour helpers # ---------------------------------------------------------------------------- def pairwise_sq_dists(X, Y=None): if Y is None: Y = X Xn = np.sum(X ** 2, axis=1, keepdims=True) Yn = np.sum(Y ** 2, axis=1, keepdims=True) D2 = Xn + Yn.T - 2.0 * (X @ Y.T) np.maximum(D2, 0, out=D2) return np.sqrt(D2) def knn_distances(D, k, self_included=True): """Return (N, k) sorted distances to k nearest neighbours. If self_included, the 0-th distance (0) is the point itself. """ N = D.shape[0] if self_included: idx = np.argpartition(D, k, axis=1)[:, :k] else: # exclude self (diagonal) Dc = D + np.eye(N) * 1e18 idx = np.argpartition(Dc, k, axis=1)[:, :k] idx = idx[np.arange(N)[:, None], np.argsort(D[np.arange(N)[:, None], idx], axis=1)] return D[np.arange(N)[:, None], idx] # ---------------------------------------------------------------------------- # ICDM (Iterative Contextual Dissimilarity Measure) # Matches official implementation: metrics/hubness_processor/hubness_reduction_methods.py # ---------------------------------------------------------------------------- def _average_k_dist(D, k): """Average of the k nearest distances (including self at distance 0), i.e. sum of (k+1) smallest distances / k.""" k_dists = np.partition(D, k, axis=1)[:, : k + 1] return k_dists.sum(axis=1) / k def icdm_scaling(D, K, n_iter=10, return_mu=False): """ICDM scaling factors delta_i, matching official `icdm_delta_low_memory`. Iteratively applies NICDM: d_{ij} <- d_{ij} / (sqrt(r_i) sqrt(r_j)) with r_i = average of k nearest dists (k=K) delta_i <- delta_i / sqrt(r_i) Final: secondary dissimilarity d^T_{ij} = d_{ij} delta_i delta_j. Returns delta_i (length N) and optionally final average-neighbour distances. """ N = D.shape[0] d = D.copy() deltas = np.ones(N) for _ in range(n_iter): r = _average_k_dist(d, K) sqrt_r = np.maximum(np.sqrt(r), 1e-12) d = d / (sqrt_r[:, None] * sqrt_r[None, :]) deltas = deltas / sqrt_r if return_mu: mu_final = _average_k_dist(d, K) return deltas, mu_final return deltas # ---------------------------------------------------------------------------- # GICDM (Algorithm 1) # ---------------------------------------------------------------------------- def gicdm(Xr, Xg, K1, K2, q=0.95, n_iter=10, return_info=False): """Generative ICDM (Algorithm 1), matching the official implementation. Xr : (N, d) real points Xg : (M, d) generated points K1, K2 : two filter scales (K2 = 10*K1 in the paper). For each scale k, ICDM is applied with neighbourhood size 2k (paper: GICDM K = 2k). Returns: D_final : (M, N) GICDM real-to-generated dissimilarity matrix keep : (M,) boolean mask of generated points passing multi-scale filter Drr_gicdm : (N, N) GICDM real-to-real dissimilarity matrix """ Drr = pairwise_sq_dists(Xr) Drg = pairwise_sq_dists(Xg, Xr) # (M, N): generated-to-real N, M = Xr.shape[0], Xg.shape[0] keep = np.ones(M, dtype=bool) info = {} for k in (K1, K2): # ICDM with neighbourhood 2k delta_r = icdm_scaling(Drr, K=2 * k, n_iter=n_iter) # --- real filter (Algorithm 1 lines 3-6) --- # k nearest real neighbours (exclude self) nn_r = np.argpartition(Drr, k, axis=1)[:, : k + 1] is_self = nn_r == np.arange(N)[:, None] no_self = ~is_self.any(axis=1) is_self[no_self, -1] = True nn_r = nn_r[~is_self].reshape(N, k) real_avg_delta = delta_r[nn_r].mean(axis=1) r_ri = np.abs(real_avg_delta - delta_r) / real_avg_delta T_k = np.quantile(r_ri, q) # --- generated filter (Algorithm 1 lines 9-14) --- # k+1 nearest real neighbours of each generated point nn_g = np.argpartition(Drg.T, k, axis=1)[:, : k + 1] # r_synthetic = average of (delta_real_nn * d_orig) over k+1 neighbours d_g_nn = Drg[np.arange(M)[:, None], nn_g] delta_r_g_nn = delta_r[nn_g] r_synthetic = (delta_r_g_nn * d_g_nn).sum(axis=1) / (k + 1) delta_g = 1.0 / r_synthetic # Eq. (1) with mu_bar absorbed synth_avg_delta = delta_r_g_nn.mean(axis=1) r_gj = np.abs(synth_avg_delta - delta_g) / synth_avg_delta keep &= (r_gj <= T_k) info[k] = dict(T_k=float(T_k), delta_g=delta_g) # final GICDM dissimilarities (Algorithm 1 line 17) delta_g_K1 = info[K1]['delta_g'] delta_r_K1 = icdm_scaling(Drr, K=2 * K1, n_iter=n_iter) D_final = Drg * (delta_r_K1[None, :] * delta_g_K1[None, :]) Drr_gicdm = Drr * np.outer(delta_r_K1, delta_r_K1) if return_info: return D_final, keep, info, (delta_r_K1, delta_g_K1) return D_final, keep, Drr_gicdm # ---------------------------------------------------------------------------- # Hubness statistics # ---------------------------------------------------------------------------- def k_occurrence(D, k=5): """O_k(x_i) = number of points for which x_i is among their k NN (excl self).""" N = D.shape[0] Dc = D + np.eye(N) * 1e18 knn = np.argsort(Dc, axis=1)[:, :k] occ = np.zeros(N, dtype=int) for j in range(k): occ[knn[:, j]] += 1 return occ def hubness_stats(D, k=5, q=0.01): """Return h5_1(1%) and A5 (proportion of antihubs).""" occ = k_occurrence(D, k) mean_occ = occ.mean() n = len(occ) topq = max(1, int(np.floor(q * n))) top_vals = np.sort(occ)[::-1][:topq] h5 = top_vals.mean() / mean_occ if mean_occ > 0 else np.nan A5 = float(np.mean(occ == 0)) return float(h5), A5 # ---------------------------------------------------------------------------- # Clipped Density / Clipped Coverage (Salvy et al. 2026) # ---------------------------------------------------------------------------- def clipped_density(Xr, Xg, k=5, dissim=None, keep=None): """Clipped Density fidelity metric. For each generated point, distance to its k-th real NN; threshold = distance from each real point to its k-th real NN (clip). Score averages clip term. If dissim ('gicdm') is provided, use GICDM dissimilarities instead of raw dist. keep: boolean mask of generated points to include (filtered-out points get 0). """ Drg = pairwise_sq_dists(Xg, Xr) # (M,N) Drr = pairwise_sq_dists(Xr) if dissim is None: d_rg = Drg d_rr = Drr else: d_rg, d_rr = dissim # k-th NN distance for each real point (in its own set) kth_real = np.sort(d_rr + np.eye(len(Xr)) * 1e18, axis=1)[:, k - 1] kth_gen = np.sort(d_rg, axis=1)[:, k - 1] # clip each generated point's distance at its matched real threshold thresh = kth_real[np.argmin(d_rg, axis=1)] clip = np.clip(kth_gen / thresh, 0, 1) if keep is not None: clip = clip * keep.astype(float) # filtered points -> fidelity 0 return float(np.mean(clip)) def clipped_coverage(Xr, Xg, k=5, dissim=None, keep=None): """Clipped Coverage: for each real point, does a kept generated point fall within its k-th NN threshold?""" Drg = pairwise_sq_dists(Xg, Xr) Drr = pairwise_sq_dists(Xr) if dissim is None: d_rg = Drg d_rr = Drr else: d_rg, d_rr = dissim kth_real = np.sort(d_rr + np.eye(len(Xr)) * 1e18, axis=1)[:, k - 1] min_d = d_rg.min(axis=0) covered = (min_d <= kth_real).astype(float) if keep is not None: # a generated point only contributes if it is kept d_rg_k = d_rg[keep, :] if d_rg_k.shape[0] == 0: return 0.0 min_d2 = d_rg_k.min(axis=0) covered = (min_d2 <= kth_real).astype(float) return float(np.mean(covered))