gicdm_core.py

"""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))