def load_dexter():
def check_is_nD_array(arr:np.ndarray, n:int, arr_type=''):
def check_distance_matrix_shape(D:np.ndarray):
def check_distance_matrix_shape_fits_vectors(D:np.ndarray, vectors:np.ndarray):
def check_distance_matrix_shape_fits_labels(D:np.ndarray, classes:np.ndarray):
def check_vector_matrix_shape_fits_labels(X:np.ndarray, classes:np.ndarray):
def check_sample_shape_fits(D:np.ndarray, idx:np.ndarray):
def check_valid_metric_parameter(metric:str):
def matrix_split(rows, cols, elem_size=8, nr_matrices=4): # pragma: no cover
def random_sparse_matrix(size, density=0.05):
def save_csr_matrix(file, matrix):
def load_csr_matrix(file):
    def __init__(self, unit='kB'):
    def _factor(self):
    def total(self):
    def used(self):
    def used_real(self):
    def shared(self):
    def buffers(self):
    def cached(self):
    def user_free(self):
    def swap(self):
    def swap_free(self):
    def swap_used(self):
    D = cosine_distance(vectors)
    S = sparse.rand(size, size, density, 'csr')
class FreeMemLinux(object): # pragma: no cover def centering(X:np.ndarray, metric:str='vector', test_set_mask:np.ndarray=None):
    """
    Perform  centering, i.e. shift the origin to the data centroid.

    Centering of vector data `X` with ``n`` objects in an ``m``-dimensional
    feature space.
    The mean of each feature is calculated and subtracted from each
    point [1]_. In distance based mode, it must be checked upstream, that
    the distance matrix is a gram matrix as described below!

    Parameters
    ----------
    X : ndarray
        - An ``(n x m)`` vector data matrix with ``n`` objects in an
          ``m``-dimensional feature space
        - An ``(n x n)`` distance matrix
          of form ``K = X(X.T)``, if `X` is an ``(n x m)`` matrix;
          and of form ``K = (X.T)X``, if `X` is an ``(m x n)`` matrix,
          where ``X.T`` denotes the transpose of `X`.

        NOTE: The type must be defined via parameter 'metric'!

    metric : {'vector', 'inner'}, optional (Default: 'vector')
        Define, whether `X` is vector data or a Gram matrix of inner product
        similarities as described above.

    test_set_mask : ndarray, optional (default: None)
        Hold back data as a test set and perform centering on the remaining 
        data (training set).

    Returns
    ------- 
    X_cent : ndarray

        Centered vectors with shape (n, m), if given vector data.

    K_cent : ndarray

        Centered inner product similarities with shape (n, n), if given Gram matrix.
        
    References
    ----------
    .. [1] Suzuki, I., Hara, K., Shimbo, M., Saerens, M., & Fukumizu, K. (2013). 
           Centering similarity measures to reduce hubs. In Proceedings of the 
           2013 Conference on Empirical Methods in Natural Language Processing 
           (pp 613–623). 
           Retrieved from https://www.aclweb.org/anthology/D/D13/D13-1058.pdf
    """
    # Kernel based centering requires inner product similarities, NOT distances.
    # Since the parameter was previously erroneously called 'distance',
    # this is kept for compatibility reasons.
    if metric in ('similarity', 'distance', 'inner', 'inner_product'):
        if test_set_mask is not None:
            raise NotImplementedError("Kernel based centering does not "
                                      "support train/test splits so far.")
        io.check_distance_matrix_shape(X)
        n = X.shape[0]
        H = np.identity(n) - np.ones((n, n)) / n
        # K = X.T.X must be provided upstream
        return H.dot(X).dot(H)
    elif metric == 'vector':
        n = X.shape[0]
        if test_set_mask is None:
            # center among all data
            return X - np.mean(X, axis=0)
        else: 
            # center among training data
            train_ind = np.setdiff1d(np.arange(n), test_set_mask)
            return X - np.mean(X[train_ind], axis=0)
    else:
        raise ValueError("Parameter 'metric' must be 'inner' or 'vector'.") def weighted_centering(X:np.ndarray, metric:str='cosine', gamma:float=1., 
                       test_set_mask:np.ndarray=None):
    """
    Perform  weighted centering: shift origin to the weighted data mean
    
    Move the origin more actively towards hub objects in the dataset, 
    rather than towards the data centroid [1]_.
    
    Parameters
    ----------
    X : ndarray
        An ``m x n`` vector data matrix with ``n`` objects in an 
        ``m`` dimensional feature space 
    
    metric : {'cosine', 'euclidean'}, optional (default: 'cosine')
        Distance measure used to place more weight on objects that are more 
        likely to become hubs. (Defined for 'cosine' in [1]_, 'euclidean' does 
        not make much sense and might be removed in the future).
    
    gamma : float, optional (default: 1.0)
        Controls how much we emphasize the weighting effect
        
        - ``gamma=0`` : equivalent to normal centering
        - ``gamma>0`` : move origin closer to objects with larger similarity 
          to other objects
    
    test_set_mask : ndarray, optional (default: None)
        Hold back data as a test set and perform centering on the remaining 
        data (training set).
    
    Returns
    ------- 
    X_wcent : ndarray
        Weighted centered vectors.
        
    References
    ----------
    .. [1] Suzuki, I., Hara, K., Shimbo, M., Saerens, M., & Fukumizu, K. (2013). 
           Centering similarity measures to reduce hubs. In Proceedings of the 
           2013 Conference on Empirical Methods in Natural Language Processing 
           (pp 613–623). 
           Retrieved from https://www.aclweb.org/anthology/D/D13/D13-1058.pdf
    """
    n = X.shape[0]
                   
    # Indices of training examples
    if test_set_mask is not None:
        train_set_mask = np.setdiff1d(np.arange(n), test_set_mask)
    else:
        train_set_mask = slice(0, n)
    
    n_train = X[train_set_mask].shape[0]
    d = np.zeros(n)
    
    if metric == 'cosine':
        vectors_sum = X[train_set_mask].sum(axis=0)
        for i in np.arange(n):
            d[i] = n_train * cos(np.array([X[i], vectors_sum/n_train]))[0, 1]
    # Using euclidean distances does not really make sense
    elif metric == 'euclidean':
        for i in range(n):
            displ_v = X[train_set_mask] - d[i]
            d[i] = np.sum(np.sqrt(displ_v * displ_v))
    else:
        raise ValueError("Weighted centering only supports cosine distances.")
    d_sum = np.sum(d ** gamma)
    w = (d ** gamma) / d_sum
    vectors_mean_weighted = np.sum(w.reshape(n, 1) * X, axis=0)
    X_wcent = X - vectors_mean_weighted
    return X_wcent def localized_centering(X:np.ndarray, Y:np.ndarray=None,
                        kappa:int=40, gamma:float=1., n_jobs:int=1):
    """
    Perform localized centering.
    
    Reduce hubness in datasets according to the method proposed in [2]_.
    
    Parameters
    ----------
    X : ndarray
        An ``n x m`` vector data matrix with ``n`` objects in an 
        ``m`` dimensional feature space

    Y : ndarray, optional
        If Y is provided, calculate similarities between all test data in `X`
        versus all training data in `Y`.
    
    kappa : int, optional (default: 40)
        Local segment size, determines the size of the local neighborhood for 
        calculating the local affinity. When ``kappa=n`` localized centering 
        reduces to standard centering.
        "select κ depending on the dataset, so that the correlation between
        Nk(x) and the local affinity <x, cκ(x)> is maximized" [2]_
        
    gamma : float, optional (default: 1.0)
        Control the degree of penalty, so that used the similarity score 
        is smaller depending on how likely a point is to become a hub.
        "Parameter γ can be tuned so as to maximally reduce the skewness 
        of the Nk distribution" [2]_.

    n_jobs : int, optional
        Parallel execution

    Returns
    ------- 
    S_lcent : ndarray
        Secondary similarity (localized centering) matrix.
        
    References
    ----------
    .. [1] Suzuki, I., Hara, K., Shimbo, M., Saerens, M., & Fukumizu, K. (2013). 
           Centering similarity measures to reduce hubs. In Proceedings of the 
           2013 Conference on Empirical Methods in Natural Language Processing 
           (pp 613–623). 
           Retrieved from https://www.aclweb.org/anthology/D/D13/D13-1058.pdf
    
    .. [2] Hara, K., Suzuki, I., Shimbo, M., Kobayashi, K., Fukumizu, K., & 
           Radovanović, M. (2015). Localized centering: Reducing hubness in 
           large-sample data hubness in high-dimensional data. In AAAI ’15: 
           Proceedings of the 29th AAAI Conference on Artificial Intelligence 
           (pp. 2645–2651).
    """
    if n_jobs == -1:
        n_jobs = cpu_count()
    # Rescale vectors to unit length
    div_ = np.sqrt((X ** 2).sum(axis=-1))[..., np.newaxis]
    div_[div_ == 0] = 1e-7
    v = X / div_
    if Y is None: # calc all-against-all in X
        w = v
        n, _ = X.shape
        sim = v.dot(w.T)
        sim_train = sim
    else: # calc sim from test data in X against train data in Y
        div_ = np.sqrt((Y ** 2).sum(axis=-1))[..., np.newaxis]
        div_[div_ == 0] = 1e-7
        w = Y / div_
        n, _ = Y.shape
        sim = v.dot(w.T)
        sim_train = w.dot(w.T)   

    if n_jobs > 1:
        local_affinity_ctype = RawArray(ctypes.c_double, n)
        local_affinity = np.frombuffer(local_affinity_ctype, dtype=np.float64)
        with Pool(processes=n_jobs,
                  initializer=_lcent_load_shared_data,
                  initargs=(w, sim_train, local_affinity)) as pool:
            for _ in pool.imap(
                func=partial(_lcent_calculate_loc_af, kappa=kappa), 
                iterable=range(n)):
                pass # local_affinity is handled within func
    else:
        local_affinity = np.zeros(n)
        for i in range(n):
            # Get the kappa nearest neighbors (highest similarity)
            nn = np.argpartition(sim_train[i, :], kth=-kappa)[-1:-kappa-1:-1]
            # Local centroid
            c_kappa_x = w[nn, :].mean(axis=0)
            local_affinity[i] = np.inner(w[i, :], c_kappa_x)
    # Only change penalty, if all values are positive 
    if gamma != 1 and (local_affinity < 0).sum() == 0:
        local_affinity **= gamma
    sim -= local_affinity
    return sim def dis_sim_global(X:np.ndarray, Y:np.ndarray=None):
    """
    Calculate dissimilarity based on global 'sample-wise centrality' [1]_.
    
    Parameters
    ----------
    X : ndarray
        An ``n x m`` vector data matrix with ``n`` objects in an 
        ``m`` dimensional feature space

    Y : ndarray, optional
        If Y is provided, calculate dissimilarities between all test data
        in `X` and all training data in `Y`.

    Returns
    -------
    D_dsg : ndarray
        Secondary dissimilarity (DisSimGlobal) matrix.
        
    References
    ----------
    .. [1] Hara, K., Suzuki, I., Kobayashi, K., Fukumizu, K., & 
           Radovanović, M. (2016). Flattening the density gradient for 
           eliminating spatial centrality to reduce hubness. Proceedings of 
           the Thirtieth AAAI Conference on Artificial Intelligence (AAAI ’16), 
           1659–1665. Retrieved from http://www.aaai.org/ocs/index.php/AAAI/
           AAAI16/paper/download/12055/11787
    """
    if Y is None:
        Y = X
    if X.shape[1] != Y.shape[1]:
        raise ValueError("X and Y must have same number of features.")
    c = Y.mean(0)
    x_c = euclidean_distances(Y, c[np.newaxis, :], squared=True)
    if id(X) == id(Y): # i.e. no Y was provided
        q_c = x_c
    else: # avoid duplicate calculations
        q_c = euclidean_distances(X, c[np.newaxis, :], squared=True)
    D_xq = euclidean_distances(X, Y, squared=True)
    D_xq -= x_c.T
    D_xq -= q_c
    return D_xq def dis_sim_local(X:np.ndarray, Y:np.ndarray=None, k:int=10, n_jobs:int=1):
    """Calculate dissimilarity based on local 'sample-wise centrality' [1]_.
    
    Parameters
    ----------
    X : ndarray
        An ``n x m`` vector data matrix with ``n`` objects in an 
        ``m`` dimensional feature space.

    Y : ndarray, optional
        If Y is provided, calculate dissimilarities between all test data
        in `X` and all training data in `Y`.

    k : int, optional (default: 10)
        Neighborhood size used for determining the local centroids.
        Can be optimized as to maximally reduce hubness [1]_.

    n_jobs : int, optional
        Parallel execution

    Returns
    -------
    D_dsl : ndarray
        Secondary dissimiliarity (DisSimLocal) matrix.
        
    References
    ----------
    .. [1] Hara, K., Suzuki, I., Kobayashi, K., Fukumizu, K., & 
           Radovanović, M. (2016). Flattening the density gradient for 
           eliminating spatial centrality to reduce hubness. Proceedings of 
           the Thirtieth AAAI Conference on Artificial Intelligence (AAAI ’16), 
           1659–1665. Retrieved from http://www.aaai.org/ocs/index.php/AAAI/
           AAAI16/paper/download/12055/11787
    """
    X = X.copy()
    # all-against-all dissimilarities?
    if Y is None:
        Y = X
    else:
        Y = Y.copy()

    # dataset size and dimensionality
    n_test, m_test = X.shape
    n_train, m_train = Y.shape
    if m_test != m_train:
        raise ValueError("X and Y must have same number of features.")

    # Calc euclidean distances to find nearest neighbors among training data
    D_train = euclidean_distances(Y, squared=True)
    if id(Y) == id(X):
        # Exclude self distances from kNN lists:
        np.fill_diagonal(D_train, np.inf)
        D_test = D_train
    else:
        # ... and between test and training data
        D_test = euclidean_distances(X, Y, squared=True)

    # Local centroid for each point among its k-nearest training neighbors
    if n_jobs > 1:
        c_k_X_ctype = RawArray(ctypes.c_double, X.size)
        c_k_X = np.frombuffer(c_k_X_ctype, dtype=np.float64).reshape(X.shape)
        with Pool(processes=n_jobs,
                  initializer=_dsl_init,
                  initargs=(c_k_X, D_test, Y)) as pool:
            for _ in pool.imap(
                func=partial(_dsl_local_centroids, k=k), 
                iterable=range(n_test)):
                pass # handling inside function
    else:
        c_k_X = np.zeros_like(X)
        for i in range(n_test):
            knn_idx = np.argpartition(D_test[i, :], kth=k)[:k]
            c_k_X[i] = Y[knn_idx].mean(axis=0)
    X -= c_k_X
    X **= 2
    x_c_k = X.sum(axis=1)
    if id(Y) == id(X):
        c_k_Y = c_k_X
        y_c_k = x_c_k
    else:
        if n_jobs > 1:
            c_k_Y_ctype = RawArray(ctypes.c_double, Y.size)
            c_k_Y = np.frombuffer(c_k_Y_ctype, dtype=np.float64).reshape(Y.shape)
            with Pool(processes=n_jobs,
                      initializer=_dsl_init,
                      initargs=(c_k_Y, D_train, Y)) as pool:
                for _ in pool.imap(
                    func=partial(_dsl_local_centroids, k=k),
                    iterable=range(n_train)):
                    pass # handling within function
        else:
            c_k_Y = np.zeros_like(Y)
            for i in range(n_train):
                knn_idx = np.argpartition(D_train[i, :], kth=k)[:k]
                c_k_Y[i] = Y[knn_idx].mean(axis=0)
        Y -= c_k_Y
        Y **= 2
        y_c_k = Y.sum(axis=1)
    # DisSimLocal
    x_y = D_test
    x_y -= x_c_k[:, np.newaxis]
    x_y -= y_c_k
    if id(Y) == id(X):
        np.fill_diagonal(x_y, -np.inf)
    return x_y def cosine_distance(X):
    """Calculate the cosine distance between all pairs of vectors in `X`."""
    xn = np.sqrt(np.sum(X**2, 1))
    Y = X / xn[:, np.newaxis]
    del xn
    D = 1. - Y.dot(Y.T)
    del Y
    D[D < 0] = 0
    D = np.triu(D, 1) + np.triu(D, 1).T
    return D def mutual_proximity_empiric(D:np.ndarray, metric:str='distance',
                             test_ind:np.ndarray=None, verbose:int=0,
                             sample_ind:np.ndarray=None, n_jobs=None,
                             min_nnz:int=0):
    """Transform a distance matrix with Mutual Proximity (empiric distribution).

    Applies Mutual Proximity (MP) [1]_ on a distance/similarity matrix using
    the empiric data distribution (EXACT, rather SLOW). The resulting
    secondary distance/similarity matrix should show lower hubness.

    Parameters
    ----------
    D : ndarray
        Distance or similarity matrix.

        - Shape ``n x n`` matrix for scaling the complete matrix.
        - Shape ``n x s`` matrix (where ``n`` and ``s`` are the dataset and
        sample size, respectively) for scaling only the distances to
        ``s`` sample instances, whose indices must be specified with
        parameter `sample_ind`.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether matrix `D` is a distance or similarity matrix.

    test_ind : ndarray, optional (default: None)
        Define data points to be hold out as part of a test set. Can be:

        - None : Rescale all distances
        - ndarray : Hold out points indexed in this array as test set.

    verbose : int, optional (default: 0)
        Increasing level of output (progress report).

    sample_ind : ndarray
        The index array that determines, which data points the columns in
        `D` correspond to (indices of training data).

    Returns
    -------
    D_mp : ndarray
        Secondary distance MP empiric matrix.

    References
    ----------
    .. [1] Schnitzer, D., Flexer, A., Schedl, M., & Widmer, G. (2012).
           Local and global scaling reduce hubs in space. The Journal of Machine
           Learning Research, 13(1), 2871–2902.
    """
    if sample_ind is None:
        return _mutual_proximity_empiric_full(D=D, metric=metric,
                                              test_set_ind=test_ind,
                                              min_nnz=min_nnz,
                                              verbose=verbose,
                                              n_jobs=n_jobs)
    else:
        return _mutual_proximity_empiric_sample(D=D, idx=sample_ind,
                                                metric=metric,
                                                test_set_ind=test_ind,
                                                verbose=verbose,
                                                n_jobs=n_jobs) def mutual_proximity_gammai(D:np.ndarray, metric:str='distance',
                            min_nnz:int=30, test_set_ind:np.ndarray=None,
                            verbose:int=0):
    """Transform a distance matrix with Mutual Proximity (indep. Gamma distr.).
    
    Applies Mutual Proximity (MP) [1]_ on a distance/similarity matrix. Gammai 
    variant assumes independent Gamma distributed distances (FAST).
    The resulting second. distance/similarity matrix should show lower hubness.

    Parameters
    ----------
    D : ndarray or csr_matrix
        - ndarray: The ``n x n`` symmetric distance or similarity matrix.
        - csr_matrix: The ``n x n`` symmetric similarity matrix.
        
        NOTE: In case of sparse `D`, zeros are interpreted as missing values 
        and ignored during calculations. Thus, results may differ 
        from using a dense version.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether matrix `D` is a distance or similarity matrix.

        NOTE: In case of sparse `D`, only 'similarity' is supported.

    min_nnz : int, optional, default: 30
        Calculate MP between two objects `i` and `j`, iff at least ``min_nnz``
        values are present in both row ``i`` and ``j``.
        Otherwise, return the original similarity.
        Ignored, if `metric` is 'distance'.

    test_set_ind : ndarray, optional (default: None)
        Define data points to be hold out as part of a test set. Can be:

        - None : Rescale all distances
        - ndarray : Hold out points indexed in this array as test set. 
 
    verbose : int, optional (default: 0)
        Increasing level of output (progress report).

    Returns
    -------
    D_mp : ndarray
        Secondary distance MP gammai matrix.

    References
    ----------
    .. [1] Schnitzer, D., Flexer, A., Schedl, M., & Widmer, G. (2012). 
           Local and global scaling reduce hubs in space. 
           The Journal of Machine Learning Research, 13(1), 2871–2902.
    """
    # Initialization
    n = D.shape[0]
    log = ConsoleLogging()

    # Checking input
    io.check_distance_matrix_shape(D)
    io.check_valid_metric_parameter(metric)
    if metric == 'similarity':
        self_value = 1
    else: # metric == 'distance':
        self_value = 0
    if test_set_ind is None:
        train_set_ind = slice(0, n)
    else:
        train_set_ind = np.setdiff1d(np.arange(n), test_set_ind)

    # Start MP
    if verbose:
        log.message('Mutual proximity Gammai rescaling started.', flush=True)
    D = D.copy()
    
    if issparse(D):
        return _mutual_proximity_gammai_sparse(D, min_nnz, test_set_ind,
                                               verbose, log)

    np.fill_diagonal(D, np.nan)

    mu = np.nanmean(D[train_set_ind], 0)
    va = np.nanvar(D[train_set_ind], 0, ddof=1)
    # Avoid downstream div/0 errors
    va[va == 0] = 1e-7
    A = (mu**2) / va
    B = va / mu

    D_mp = np.zeros_like(D)

    # MP gammai
    for i in range(n):
        if verbose and ((i+1)%1000 == 0 or i+1 == n):
            log.message("MP_gammai: {} of {}".format(i+1, n), flush=True)
        j_idx = slice(i+1, n)
        
        if metric == 'similarity':
            p1 = _local_gamcdf(D[i, j_idx], A[i], B[i])
            p2 = _local_gamcdf(D[j_idx, i], A[j_idx], B[j_idx])
            D_mp[i, j_idx] = (p1 * p2).ravel()
        else: # distance
            p1 = 1 - _local_gamcdf(D[i, j_idx], A[i], B[i])
            p2 = 1 - _local_gamcdf(D[j_idx, i], A[j_idx], B[j_idx])
            D_mp[i, j_idx] = (1 - p1 * p2).ravel()

    # Mirroring the matrix
    D_mp += D_mp.T
    # set correct self dist/sim
    np.fill_diagonal(D_mp, self_value)

    return D_mp def mutual_proximity_gaussi(D:np.ndarray, metric:str='distance',
                            sample_size:int=0, min_nnz:int=30,
                            test_set_ind:np.ndarray=None,
                            verbose:int=0, idx:np.ndarray=None,):
    """Transform distances with Mutual Proximity (indep. normal distributions).
    
    Applies Mutual Proximity (MP) [1]_ on a distance/similarity matrix. Gaussi 
    variant assumes independent normal distributions (FAST).
    The resulting second. distance/similarity matrix should show lower hubness.
    
    Parameters
    ----------
    D : ndarray or csr_matrix
        - ndarray: The ``n x n`` symmetric distance or similarity matrix.
        - csr_matrix: The ``n x n`` symmetric similarity matrix.
        
        NOTE: In case of sparse `D`, zeros are interpreted as missing values 
        and ignored during calculations. Thus, results may differ 
        from using a dense version.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether matrix `D` is a distance or similarity matrix.
        
        NOTE: In case of sparse `D`, only 'similarity' is supported.
        
    sample_size : int, optional (default: 0)
        Define sample size from which Gauss parameters are estimated.
        Use all data when set to ``0``.
        Ignored in case of SampleMP (i.e. if provided `idx`).

    min_nnz : int, optional, default: 30
        Calculate MP between two objects `i` and `j`, iff at least ``min_nnz``
        values are present in both row ``i`` and ``j``.
        Otherwise, return the original similarity.
        Ignored, if `metric` is 'distance'.

    test_set_ind : ndarray, optional (default: None)
        Define data points to be hold out as part of a test set. Can be:
        
        - None : Rescale all distances
        - ndarray : Hold out points indexed in this array as test set.

        Ignored in case of SampleMP (i.e. if provided `idx`).

    verbose : int, optional (default: 0)
        Increasing level of output (progress report).

    idx : ndarray, optional (default: None)
        The index array that determines to which data points the columns in
        `D` correspond. Only required for SampleMP.

    Returns
    -------
    D_mp : ndarray
        Secondary distance MP gaussi matrix.
    
    References
    ----------
    .. [1] Schnitzer, D., Flexer, A., Schedl, M., & Widmer, G. (2012). 
           Local and global scaling reduce hubs in space. The Journal of Machine 
           Learning Research, 13(1), 2871–2902.
    """    
    # Initialization   
    log = ConsoleLogging()

    # Checking input
    if idx is None:
        io.check_distance_matrix_shape(D)
    else:
        io.check_sample_shape_fits(D, idx)
    io.check_valid_metric_parameter(metric)
    n = D.shape[0]
    s = D.shape[1]

    if metric == 'similarity':
        self_value = 1
    else: # metric == 'distance':
        self_value = 0  
    if test_set_ind is None:
        train_set_ind = slice(0, n)
    else:
        train_set_ind = np.setdiff1d(np.arange(n), test_set_ind)

    # Start MP Gaussi    
    if verbose:
        log.message('Mutual Proximity Gaussi rescaling started.', flush=True)
    D = D.copy()

    if issparse(D):
        return _mutual_proximity_gaussi_sparse(D, sample_size, min_nnz,
                                               test_set_ind, verbose, log)

    # ignore self dist/sim for parameter estimation
    if idx is None:
        np.fill_diagonal(D, np.nan)
    else:
        for j, i in enumerate(idx):
            D[i, j] = np.nan

    # Calculate mean and std
    if idx is None:
        if sample_size == 0:
            mu = np.nanmean(D[train_set_ind], 0)
            sd = np.nanstd(D[train_set_ind], 0, ddof=0)
        else:
            samples = np.random.shuffle(train_set_ind)[0:sample_size]
            mu = np.nanmean(D[samples], 0)
            sd = np.nanstd(D[samples], 0, ddof=0)
    else:
        mu = np.nanmean(D, 1)
        sd = np.nanstd(D, 1, ddof=0)
    # Avoid downstream div/0 errors
    sd[sd == 0] = 1e-7
    # set self dist/sim back to self_value to avoid scipy warnings
    if idx is None:
        np.fill_diagonal(D, self_value)
    else:
        for j, i in enumerate(idx):
            D[i, j] = self_value

    # MP Gaussi
    D_mp = np.zeros_like(D)
    for i in range(n):
        if verbose and ((i+1)%1000 == 0 or i+1 == n):
            log.message("MP_gaussi: {} of {}.".format(i+1, n), flush=True)
        if idx is None:
            j = slice(i+1, n)
            j_mom = j
        else:
            j = slice(0, s)
            j_mom = idx[j]
        
        if metric == 'similarity':
            p1 = norm.cdf(D[i, j], mu[i], sd[i])
            p2 = norm.cdf(D[i, j], mu[j_mom], sd[j_mom])
            D_mp[i, j] = (p1 * p2).ravel()
        else:
            # sf(.) := 1 - cdf(.)
            p1 = norm.sf(D[i, j], mu[i], sd[i])
            p2 = norm.sf(D[i, j], mu[j_mom], sd[j_mom])
            D_mp[i, j] = (1 - p1 * p2).ravel()

    if idx is None:
        D_mp += D_mp.T        
        np.fill_diagonal(D_mp, self_value)
    else:
        # Ensure correct self distances
        for j, sample in enumerate(idx):
            D_mp[sample, j] = self_value
    return D_mp def goodman_kruskal_index(D:np.ndarray, classes:np.ndarray,
                          metric:str='distance') -> float:
    """Calculate the Goodman-Kruskal clustering index.
        
    Parameters
    ----------
    D : ndarray
        The ``n x n`` symmetric distance (similarity) matrix.
    
    classes : ndarray
        The ``1 x n`` vector of class labels for each point.
    
    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether the matrix `D` is a distance or similarity matrix
    
    Returns
    -------
    gamma : float
        Goodman-Kruskal index in ``[-1, 1]`` (higher=better)
    
    Notes
    -----
    This clustering quality measure relates the number of concordant (`Q_c`) 
    and discordant (`Q_d`) quadruples (`d_ij`, `d_kl`) of a distance matrix.
    We only consider tuples, so that `i`, `j` are from the same class 
    and `k`, `l` are from different classes. Then a quadruple is...
    concordant, if 
    
    .. math:: 
        d_{i,j} < d_{k,l}
    
    discordant, if 
    
    .. math:: 
        d_{i,j} > d_{k,l}
    
    and not counted, otherwise.
    
    The Goodman-Kruskal index gamma is then defined as: 
    
    .. math:: 
        gamma = \\frac{Q_c - Q_d}{Q_c + Q_d}
        
    `gamma` is bounded to ``[-1, 1]``, where larger values indicate better 
    clustering.
    """
    
    # Checking input
    io.check_distance_matrix_shape(D)
    io.check_distance_matrix_shape_fits_labels(D, classes)
    io.check_valid_metric_parameter(metric)
    
    # Calculations
    Q_c = 0.0
    Q_d = 0.0
    cls = np.unique(classes)
    
    # D_kl pairs in different classes
    other = classes[:, np.newaxis] != classes[np.newaxis, :]
    D_other = D[np.triu(other, 1)]
    
    for c in cls:
        sel = classes == c
        if np.sum(sel) > 1:
            sel = sel[:, np.newaxis].astype(bool)
            selD = np.logical_and(sel, sel.T)
            # D_ij pairs within same class
            D_self = D[np.triu(selD, 1).astype(bool).T].T
        else:
            # skip if there is only one item per class
            continue
        # D_kl pairs in different classes (D_other) are computed once for all c
        D_full = np.append(D_self, D_other)

        self_size = np.max(np.shape(D_self))
        other_size = np.max(np.shape(D_other))
        # Sort algorithm must be stable!
        full_idx = np.argsort(D_full, kind='mergesort')[::-1]
        
        # Calc number of quadruples with equal distance
        n_equidistant = 0
        sdf = np.sort(D_full, axis=None)
        equi_mask = np.zeros(sdf.size, dtype=bool)
        # Positions with repeated values
        equi_mask[1:] = sdf[1:] == sdf[:-1]
        equi_dist = sdf[equi_mask]
        # How often does each value occur in self/other:
        for dist in np.unique(equi_dist):
            equi_arg = np.where(D_full == dist)[0]
            self_equi = (equi_arg < self_size).sum()
            other_equi = len(equi_arg) - self_equi
            # Number of dc that are actually equal
            n_equidistant += self_equi * other_equi

        # Calc number of concordant quadruples
        cc = 0
        ccsize = other_size
        for idx in full_idx:
            if idx < self_size:
                cc += ccsize
            else:
                ccsize -= 1

        # Calc number of discordant quadruples
        dc = self_size * other_size - cc - n_equidistant

        Q_c += cc
        Q_d += dc
    
    # Calc Goodman-Kruskal's gamma
    if Q_c + Q_d == 0:
        gamma = 0.0
    else:
        if metric == 'similarity':
            gamma = (Q_c - Q_d) / (Q_c + Q_d)
        else:
            gamma = (Q_d - Q_c) / (Q_c + Q_d)

    return gamma  def hubness(D:np.ndarray, k:int=5, metric='distance',
            verbose:int=0, n_jobs:int=1,
            random_state=None, shuffle_equal=True):
    """Compute hubness of a distance matrix.

    Hubness [1]_ is the skewness of the `k`-occurrence histogram (reverse
    nearest neighbor count, i.e. how often does a point occur in the
    `k`-nearest neighbor lists of other points).

    Parameters
    ----------
    D : ndarray
        The ``n x n`` symmetric distance (similarity) matrix or
        an ``n x m`` partial distances matrix (e.g. for train/test splits,
        with test objects in rows, train objects in column)
        
        NOTE: Partial distance matrices MUST NOT contain self distances.

    k : int, optional (default: 5)
        Neighborhood size for `k`-occurrence.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether matrix `D` is a distance or similarity matrix

    verbose : int, optional (default: 0)
        Increasing level of output (progress report).

    n_jobs : int, optional (default: 1)
        Number of parallel processes spawned for hubness calculation.
        Value 1 (default): One process (not using multiprocessing)
        Value (-1): As many processes as number of available CPUs.

    random_state : int, optional
        Seed the RNG for reproducible results.
        
        NOTE: Currently only compatible with `n_jobs`=1

    shuffle_equal : bool, optional
        If true, shuffle neighbors with identical distances to avoid
        artifact hubness.
        NOTE: This is especially useful for secondary distance measures
        with a restricted number of possible values, e.g. SNN or MP empiric.

    Returns
    -------
    S_k : float
        Hubness (skewness of `k`-occurrence distribution)
    D_k : ndarray
        `k`-nearest neighbor lists
    N_k : ndarray
        `k`-occurrence list

    References
    ----------
    .. [1] Radovanović, M., Nanopoulos, A., & Ivanović, M. (2010).
           Hubs in Space : Popular Nearest Neighbors in High-Dimensional Data.
           Journal of Machine Learning Research, 11, 2487–2531. Retrieved from
           http://jmlr.csail.mit.edu/papers/volume11/radovanovic10a/
           radovanovic10a.pdf
    """
    # Don't use multiprocessing environment when using only one job
    if n_jobs == 1:
        return _hubness_no_multiprocessing(D=D,
                                           k=k,
                                           metric=metric,
                                           verbose=verbose,
                                           random_state=random_state,
                                           shuffle_equal=shuffle_equal)
    if random_state is not None:
        raise ValueError("Seeding the RNG is not compatible with using n_jobs > 1.")
    log = ConsoleLogging()
    io.check_is_nD_array(arr=D, n=2, arr_type='Distance')
    io.check_valid_metric_parameter(metric)
    n, m = D.shape
    if k >= m:
        k_old = k
        k = m - 1
        log.warning("Reducing k from {} to {}, so that it is less than "
                    "the total number of neighbors.".format(k_old, k))
    if metric == 'distance':
        d_self = np.inf
        sort_order = 1
        kth = np.arange(k)
    if metric == 'similarity':
        d_self = -np.inf
        sort_order = -1
        kth = np.arange(m - k, m)

    if verbose:
        log.message("Hubness calculation (skewness of {}-occurrence)".format(k))

    # Initialization
    D = D.copy()
    D_k = np.zeros((n, k), dtype=np.float64)

    if issparse(D):
        pass # correct self-distance must be ensured upstream for sparse
    else:
        if n == m:
            # Set self dist to inf
            np.fill_diagonal(D, d_self)
        else:
            pass # Partial distance matrices MUST NOT contain self distances
        # make non-finite (NaN, Inf) appear at the end of the sorted list
        D[~np.isfinite(D)] = d_self

    # Parallelization
    if n_jobs == -1: # take all cpus
        NUMBER_OF_PROCESSES = mp.cpu_count() # @UndefinedVariable
    else:
        NUMBER_OF_PROCESSES = n_jobs
    D_k_ctype = RawArray(ctypes.c_int32, n*k)
    D_k = np.frombuffer(D_k_ctype, dtype=np.int32).reshape((n, k))
    with Pool(processes=NUMBER_OF_PROCESSES,
              initializer=_hubness_load_shared_data,
              initargs=(D, D_k, )) as pool:
        for _ in pool.imap(
            func=partial(_hubness_nearest_neighbors, n=n, m=m, 
                         d_self=d_self, metric=metric, kth=kth, 
                         sort_order=sort_order, log=log, verbose=verbose,
                         shuffle_equal=shuffle_equal),
            #chunksize=int(1e2),
            iterable=range(n)):
            pass # results handled within func

    # k-occurrence
    N_k = np.bincount(D_k.astype(int).ravel(), minlength=m)
    # Hubness
    S_k = stats.skew(N_k)

    if verbose:
        log.message("Hubness calculation done.", flush=True)

    # return hubness, k-nearest neighbors, N occurence
    return S_k, D_k, N_k def intrinsic_dimension(X:np.ndarray, k1:int=6, k2:int=12,
                        estimator:str='mackay', metric:str='vector',
                        trafo:str=None, mem_threshold:int=5000):
    """Calculate intrinsic dimension based on the MLE by Levina and Bickel [1]_.

    Parameters
    ----------
    X : ndarray
        - An ``m x n`` vector data matrix with ``n`` objects in an
          ``m`` dimensional feature space
        - An ``n x n`` distance matrix.

        NOTE: The type must be defined via parameter `metric`!

    k1 : int, optional (default: 6)
        Start of neighborhood range to search in.

    k2 : int, optional (default: 12)
        End of neighborhood range to search in.

    estimator : {'levina', 'mackay'}, optional (default: 'mackay')
        Determine the summation strategy: see [2]_.

    metric : {'vector', 'distance'}, optional (default: 'vector')
        Determine data type of `X`.

        NOTE: the MLE was derived for euclidean distances. Using
        other dissimilarity measures may lead to undefined results.

    trafo : {None, 'std', 'var'}, optional (default: None)
        Transform vector data. 

        - None: no transformation
        - 'std': standardization
        - 'var': subtract mean, divide by variance (default behavior of
          Laurens van der Maaten's DR toolbox; most likely for other
          ID/DR techniques).

    mem_treshold : int, optional, default: 5000
        Controls speed-memory usage trade-off: If number of points is higher
        than the given value, don't calculate complete distance matrix at
        once (fast, high memory), but per row (slower, less memory).

    Returns
    -------
    d_mle : int
        Intrinsic dimension estimate (rounded to next integer)

    References
    ----------
    .. [1] Levina, E., & Bickel, P. (2004). Maximum likelihood estimation of
           intrinsic dimension. Advances in Neural Information …, 17, 777–784.
           http://doi.org/10.2307/2335172
    .. [2] http://www.inference.phy.cam.ac.uk/mackay/dimension/
    """
    n = X.shape[0]
    if estimator not in ['levina', 'mackay']:
        raise ValueError("Parameter 'estimator' must be 'levina' or 'mackay'.")
    if k1 < 1 or k2 < k1 or k2 >= n:
        raise ValueError("Invalid neighborhood: Please make sure that "
                         "0 < k1 <= k2 < n. (Got k1={} and k2={}).".
                         format(k1, k2))
    X = X.copy().astype(float)

    if metric == 'vector':
        # New array with unique rows
        X = X[np.lexsort(np.fliplr(X).T)]
        
        if trafo is None:
            pass
        elif trafo == 'var':
            X -= X.mean(axis=0) # broadcast
            X /= X.var(axis=0) + 1e-7 # broadcast
        elif trafo == 'std':
            # Standardization
            X -= X.mean(axis=0) # broadcast
            X /= X.std(axis=0) + 1e-7 # broadcast
        else:
            raise ValueError("Transformation must be None, 'std', or 'var'.")

        # Compute matrix of log nearest neighbor distances
        X2 = (X**2).sum(1)

        if n <= mem_threshold: # speed-memory trade-off
            distance = X2.reshape(-1, 1) + X2 - 2*np.dot(X, X.T) #2x br.cast
            distance.sort(1)
            # Replace invalid values with a small number
            distance[distance<=0] = 1e-7
            knnmatrix = .5 * np.log(distance[:, 1:k2+1])
        else:
            knnmatrix = np.zeros((n, k2))
            for i in range(n):
                distance = np.sort(X2[i] + X2 - 2 * np.dot(X, X[i, :]))
                # Replace invalid values with a small number
                distance[distance <= 0] = 1e-7
                knnmatrix[i, :] = .5 * np.log(distance[1:k2+1])
    elif metric == 'distance':
        raise NotImplementedError("ID currently only supports vector data.")
        # XXX perhaps map to sufficiently high dim with MDS, then calc ID??
        #=======================================================================
        # # TODO calculation WRONG
        # X.sort(1)
        # X[X < 0] = 1e-7
        # knnmatrix = np.log(X[:, 1:k2+1])
        #=======================================================================
    elif metric == 'similarity':
        raise NotImplementedError("ID currently only supports vector data.")
        #=======================================================================
        # # TODO calculation WRONG
        # print("WARNING: using similarity data may return "
        #       "undefined results.", file=sys.stderr)
        # X[X < 0] = 0
        # distance = 1 - (X / X.max())
        # knnmatrix = np.log(distance[:, 1:k2+1])
        #=======================================================================
    else:
        raise ValueError("Parameter `metric` must be 'vector'.")

    # Compute the ML estimate
    S = np.cumsum(knnmatrix, 1)
    indexk = np.arange(k1, k2+1) # broadcasted afterwards
    dhat = -(indexk - 2) / (S[:, k1-1:k2] - knnmatrix[:, k1-1:k2] * indexk)
    if estimator == 'levina':
        # Average over estimates and over values of k
        no_dims = dhat.mean()
    if estimator == 'mackay':
        # Average over inverses
        dhat **= -1
        dhat_k = dhat.mean(0)
        no_dims = (dhat_k ** -1).mean()
    return no_dims def score(D:np.ndarray, target:np.ndarray, k=5,
          metric:str='distance', test_set_ind:np.ndarray=None, verbose:int=0,
          sample_idx=None, filter_self=True):
    """Perform `k`-nearest neighbor classification.

    Use the ``n x n`` symmetric distance matrix `D` and target class
    labels `target` to perform a `k`-NN experiment (leave-one-out
    cross-validation or evaluation of test set; see parameter `test_set_ind`).
    Ties are broken by the nearest neighbor.

    Parameters
    ----------
    D : ndarray
        The ``n x n`` symmetric distance (similarity) matrix.

    target : ndarray (of dtype=int)
        The ``n x 1`` target class labels (ground truth).

    k : int or array_like (of dtype=int), optional (default: 5)
        Neighborhood size for `k`-NN classification.
        For each value in `k`, one `k`-NN experiment is performed.

        HINT: Providing more than one value for `k` is a cheap means to perform
        multiple `k`-NN experiments at once. Try e.g. ``k=[1, 5, 20]``.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether matrix `D` is a distance or similarity matrix

    test_sed_ind : ndarray, optional (default: None)
        Define data points to be hold out as part of a test set. Can be:

        - None : Perform a LOO-CV experiment
        - ndarray : Hold out points indexed in this array as test set. Fit
          model to remaining data. Evaluate model on test set.

    verbose : int, optional (default: 0)
        Increasing level of output (progress report).

    sample_idx : ...
        TODO add description

    filter_self : bool, optional, default: True
        Remove self similarities from sparse ``D``.
        This assumes that the highest similarity per row is the self
        similarity.
        
        NOTE: Quadratic dense matrices are always filtered for self
        distances/similarities, even if `filter_self` is set t0 `False`.
        
    Returns
    -------
    acc : ndarray (shape=(n_k x 1), dtype=float)
        Classification accuracy (`n_k`... number of items in parameter `k`)

        HINT: Refering to the above example... 
        ... ``acc[0]`` gives the accuracy of the ``k=1`` experiment.
    corr : ndarray (shape=(n_k x n), dtype=int)
        Raw vectors of correctly classified items

        HINT: ... ``corr[1, :]`` gives these items for the ``k=5`` experiment.
    cmat : ndarray (shape=(n_k x n_t x n_t), dtype=int) 
        Confusion matrix (``n_t`` number of unique items in parameter target)

        HINT: ... ``cmat[2, :, :]`` gives the confusion matrix of
        the ``k=20`` experiment.
    """

    # Check input sanity
    log = ConsoleLogging()
    if sample_idx is None:
        io.check_distance_matrix_shape(D)
    else:
        io.check_sample_shape_fits(D, sample_idx)
    io.check_distance_matrix_shape_fits_labels(D, target)
    io.check_valid_metric_parameter(metric)
    if metric == 'distance':
        d_self = np.inf
        sort_order = 1
    if metric == 'similarity':
        d_self = -np.inf
        sort_order = -1

    # Copy, because data is changed
    D = D.copy()
    target = target.astype(int)
    D_is_sparse = issparse(D)

    if verbose:
        log.message("Start k-NN experiment.")
    # Handle LOO-CV vs. test set mode
    if test_set_ind is None:
        n = D.shape[0]
        test_set_ind = range(n)    # dummy 
        train_set_ind = n   # dummy
    else:
        # number of points to be classified
        n = test_set_ind.size
        # Indices of training examples
        train_set_ind = np.setdiff1d(np.arange(n), test_set_ind)
        if sample_idx is not None:
            raise NotImplementedError("Sample k-NN does not support train/"
                                      "test splits at the moment.")
    # Number of k-NN parameters
    try:
        k_length = k.size
    except AttributeError as e:
        if isinstance(k, int):
            k = np.array([k])
            k_length = k.size
        elif isinstance(k, list):
            k = np.array(k)
            k_length = k.size
        else:
            raise e

    acc = np.zeros((k_length, 1))
    corr = np.zeros((k_length, D.shape[0]))

    cl = np.sort(np.unique(target))
    if D_is_sparse:
        # Add a label for unknown class (object w/o nonzero sim to any others)
        cl = np.append(cl, cl.max()+1)
        n_classes = len(cl) + 1
    else:
        n_classes = len(cl)
    cmat = np.zeros((k_length, n_classes, n_classes))

    classes = target.copy()
    for idx, cur_class in enumerate(cl):
        # change labels to 0, 1, ..., len(cl)-1
        classes[target == cur_class] = idx
    if sample_idx is not None:
        sample_classes = classes[sample_idx]
        j = np.ones(n, int)
        j *= (n+1) # illegal indices will throw index out of bounds error
        j[sample_idx] = np.arange(len(sample_idx))
        for j, sample in enumerate(sample_idx):
            D[sample, j] = d_self
    cl = range(len(cl))

    rnd_classif = np.zeros(k_length)
    # Classify each point in test set
    for i in test_set_ind:
        if verbose and ((i+1)%1000==0 or i+1==n):
            log.message("Prediction: {} of {}.".format(i+1, n), flush=True)

        seed_class = classes[i]

        if D_is_sparse:
            row = D.getrow(i)
        else:
            row = D[i, :]
            if sample_idx is None:
                row[i] = d_self

        # Sort points in training set according to distance
        # Randomize, in case there are several points of same distance
        # (this is especially relevant for SNN rescaling)
        if sample_idx is None:
            rp = train_set_ind
        else:
            rp = np.arange(len(sample_idx))
        if D_is_sparse:
            nnz = row.nnz
            rp = np.random.permutation(nnz)
            d2 = row.data[rp]
            # Partition for each k value
            kth = nnz - k - 1
            # sort the two highest similarities to end
            kth = np.append(kth, [nnz-2, nnz-1])
            # Clip negative indices (nnz < k)
            np.clip(kth, a_min=0, a_max=nnz-1, out=kth)
            # Remove duplicate k values and sort
            kth = np.unique(kth)
            d2idx = np.argpartition(d2, kth=kth)
            d2idx = d2idx[~np.isnan(d2[d2idx])][::-1]
            idx = row.nonzero()[1][rp[d2idx]]
            idx = idx[1:] # rem self sim
        else:
            rp = np.random.permutation(rp)
            d2 = row[rp]
            d2idx = np.argsort(d2, axis=0)[::sort_order]
            d2idx = d2idx[~np.isnan(d2[d2idx])] # filter NaN values
            idx = rp[d2idx]

        # More than one k is useful for cheap multiple k-NN experiments at once
        for j in range(k_length):
            # Make sure no inf/-inf/nan values are used for classification
            if D_is_sparse:
                #print(row[0, idx[0:k[j]]].toarray())
                finite_val = np.isfinite(row[0, idx[0:k[j]]].toarray().ravel())
                #print(finite_val)
            else:
                finite_val = np.isfinite(row[idx[0:k[j]]])
            # However, if no values are finite, classify randomly
            if finite_val.sum() == 0:
                idx = np.random.permutation(idx)
                finite_val = np.ones_like(finite_val)
                rnd_classif[j] += 1
            if sample_idx is None:
                nn_class = classes[idx[0:k[j]]][finite_val]
            else:
                #finite_val = np.isfinite(sample_row[idx[0:k[j]]])
                nn_class = sample_classes[idx[0:k[j]]][finite_val]
            cs = np.bincount(nn_class.astype(int))
            if cs.size > 0:
                max_cs = np.where(cs == np.max(cs))[0]
            else:
                max_cs = np.array([len(cl) - 1]) # misclassification label

            # "tie": use nearest neighbor
            if len(max_cs) > 1:
                if seed_class == nn_class[0]:
                    acc[j] += 1/n 
                    corr[j, i] = 1
                cmat[j, seed_class, nn_class[0]] += 1
            # majority vote
            else:
                if cl[max_cs[0]] == seed_class:
                    acc[j] += 1/n
                    corr[j, i] = 1
                cmat[j, seed_class, cl[max_cs[0]]] += 1

    if np.any(rnd_classif):
        for x in rnd_classif:
            log.warning(("{} queries were classified randomly, because all "
                        "distances were non-finite numbers.").format(x))
    if verbose:
        log.message("Finished k-NN experiment.")

    return acc, corr, cmat def nicdm(D:np.ndarray, k:int=7, metric:str='distance',
          test_ind:np.ndarray=None, n_jobs:int=1):
    """Transform a distance matrix with local scaling variant NICDM.

    Transforms the given distance matrix into new one using NICDM [1]_
    with the given neighborhood radius `k` (average). There are two types of
    local scaling methods implemented. The original one and the non-iterative
    contextual dissimilarity measure, both reduce hubness in distance spaces,
    similarly to Mutual Proximity.

    Parameters
    ----------
    D : ndarray
        The ``n x n`` symmetric distance (similarity) matrix.

    k : int, optional (default: 7)
        Neighborhood radius for local scaling.

    metric : {'distance'}, optional (default: 'distance')
        Currently, only distance matrices are supported.

    test_ind : ndarray, optional (default: None)
        Define data points to be hold out as part of a test set. Can be:

        - None : Rescale all distances
        - ndarray : Hold out points indexed in this array as test set.

    n_jobs : int, optional, default: 1
        Number of processes for parallel computations.

        - `1`: Don't use multiprocessing.
        - `-1`: Use all CPUs

    Returns
    -------
    D_nicdm : ndarray
        Secondary distance NICDM matrix.

    References
    ----------
    .. [1] Schnitzer, D., Flexer, A., Schedl, M., & Widmer, G. (2012).
           Local and global scaling reduce hubs in space. The Journal of Machine
           Learning Research, 13(1), 2871–2902.
    """
    #log = logging.ConsoleLogging()
    # Checking input
    io.check_distance_matrix_shape(D)
    io.check_valid_metric_parameter(metric)
    if metric == 'similarity':
        raise NotImplementedError("NICDM does not support similarity matrices "
                                  "at the moment.")
    else:
        kth = np.arange(k)
        exclude = np.inf
    if n_jobs == -1:
        n_jobs = cpu_count()
    D = np.copy(D)
    n = D.shape[0]

    if test_ind is None:
        train_ind = slice(0, n)
    else:
        train_ind = np.setdiff1d(np.arange(n), test_ind)

    np.fill_diagonal(D, exclude)
    if n_jobs > 1:
        r_ctype = RawArray(ctypes.c_double, n)
        r = np.frombuffer(r_ctype, dtype=np.float64)
        with Pool(processes=n_jobs,
                  initializer=_nicdm_load_shared_data,
                  initargs=(D, train_ind, r, r_ctype)) as pool:
            for i, knn in enumerate(pool.imap(
                func=partial(_nicdm_calculate_r, kth=kth, k=k),
                iterable=range(n))):
                pass # r is handled within func
            r_geom = _local_geomean(r)
        D_nicdm_ctype = RawArray(ctypes.c_double, D.size)
        D_nicdm = np.frombuffer(D_nicdm_ctype, dtype=np.float64).reshape(D.shape)
        with Pool(processes=n_jobs,
                  initializer=_nicdm_load_shared_data,
                  initargs=(D, train_ind, r, r_ctype, D_nicdm, D_nicdm_ctype)) as pool:
            for _ in pool.imap(
                func=partial(_nicdm_calculate_sec_dist, r_geom=r_geom, n=n, metric=metric),
                iterable=range(n)):
                pass # results handled within func
    else: # no multiprocessing
        knn = np.partition(D[:, train_ind], kth=kth, axis=1)[:, :k]
        r = knn.mean(axis=1)
        r_geom = _local_geomean(r)

        D_nicdm = np.zeros_like(D)
        for i in range(n):
            # vectorized inner loop for 100x speed-up (using broadcasting)
            #D_nicdm[i, i+1:] = ((r_geom**2) * D[i, i+1:]) / (r[i] * r[i+1:])
            D_nicdm[i, i+1:] = (r_geom * D[i, i+1:]) / np.sqrt(r[i] * r[i+1:])
        D_nicdm += D_nicdm.T

    return D_nicdm def local_scaling(D:np.ndarray, k:int=7, metric:str='distance',
                  test_ind:np.ndarray=None, n_jobs:int=1):
    """Transform a distance matrix with Local Scaling.

    Transforms the given distance matrix into new one using local scaling [1]_
    with the given `k`-th nearest neighbor. There are two types of local
    scaling methods implemented. The original one and NICDM, both reduce
    hubness in distance spaces, similarly to Mutual Proximity.

    Parameters
    ----------
    D : ndarray or csr_matrix
        The ``n x n`` symmetric distance (similarity) matrix.

    k : int, optional (default: 7)
        Neighborhood radius for local scaling.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether matrix `D` is a distance or similarity matrix.

        NOTE: self similarities in sparse `D_ls` are set to ``np.inf``

    test_ind : ndarray, optional (default: None)
        Define data points to be hold out as part of a test set. Can be:

        - None : Rescale all distances
        - ndarray : Hold out points indexed in this array as test set.

    n_jobs : int, optional, default: 1
        Number of processes for parallel computations.

        - `1`: Don't use multiprocessing.
        - `-1`: Use all CPUs

    Returns
    -------
    D_ls : ndarray
        Secondary distance LocalScaling matrix.

    References
    ----------
    .. [1] Schnitzer, D., Flexer, A., Schedl, M., & Widmer, G. (2012).
           Local and global scaling reduce hubs in space. The Journal of Machine
           Learning Research, 13(1), 2871–2902.
    """
    log = ConsoleLogging()
    # Checking input
    io.check_distance_matrix_shape(D)
    io.check_valid_metric_parameter(metric)
    sparse = issparse(D)
    n = D.shape[0]
    if n_jobs == -1:
        n_jobs = cpu_count()
    if metric == 'similarity':
        kth = n - k
        exclude = -np.inf
        self_tmp_value = np.inf
        self_value = 1.
        log.warning("Similarity matrix support for LS is experimental.")
        if sparse and n_jobs != 1:
            log.warning("Parallel processing not implemented for sparse "
                        "matrices. Using single process instead.")
            n_jobs = 1
    else: # metric == 'distance':
        kth = k - 1
        exclude = np.inf
        self_value = 0
        self_tmp_value = self_value
        if sparse:
            log.error("Sparse distance matrices are not supported.")
            raise NotImplementedError(
                "Sparse distance matrices are not supported.")
    D = np.copy(D)

    if test_ind is None:
        train_ind = slice(0, n) #take all        
    else:
        train_ind = np.setdiff1d(np.arange(n), test_ind)
    if sparse:
        r = np.zeros(n)
        for i in range(n):
            di = D[i, train_ind].toarray()
            di[i] = exclude
            r[i] = np.partition(di, kth=kth)[kth]
        D_ls = lil_matrix(D.shape)
        # Number of nonzero cells per row
        nnz = D.getnnz(axis=1)
    else:
        np.fill_diagonal(D, exclude)
        if n_jobs > 1:
            r_ctype = RawArray(ctypes.c_double, n)
            r = np.frombuffer(r_ctype, dtype=np.float64)
            with Pool(processes=n_jobs,
                      initializer=_ls_load_shared_data,
                      initargs=(D, train_ind, r, r_ctype)) as pool:
                for _ in pool.imap(func=partial(_ls_calculate_r, kth=kth),
                                   iterable=range(n)):
                    pass # results handled within func
        else:
            r = np.partition(D[:, train_ind], kth=kth)[:, kth]

    if sparse or n_jobs == 1:
        D_ls = np.zeros_like(D)
        for i in range(n):
            # vectorized inner loop: calc only triu part
            tmp = np.empty(n-i)
            tmp[0] = self_tmp_value
            if metric == 'similarity':
                if sparse and nnz[i] <= k:  # Don't rescale if there are
                    tmp[1:] = np.nan        # too few neighbors in row
                else:
                    tmp[1:] = np.exp(-1 * D[i, i+1:]**2 / (r[i] * r[i+1:]))
            else:
                tmp[1:] = 1 - np.exp(-1 * D[i, i+1:]**2 / (r[i] * r[i+1:]))
            D_ls[i, i:] = tmp
        # copy triu to tril -> symmetric matrix (diag=zeros)
        # NOTE: does not affect self values, since inf+inf=inf and 0+0=0
        D_ls += D_ls.T
    else:
        D_ls_ctype = RawArray(ctypes.c_double, D.size)
        D_ls = np.frombuffer(D_ls_ctype, dtype=np.float64).reshape(D.shape)
        with Pool(processes=n_jobs,
                  initializer=_ls_load_shared_data,
                  initargs=(D, train_ind, r, r_ctype, D_ls, D_ls_ctype)) as pool:
            for _ in pool.imap(func=partial(_ls_calculate_sec_dist,
                                  n=n, metric=metric,
                                  self_tmp_value=self_tmp_value),
                               iterable=range(n)):
                pass # results handled within func
        # triu is copied to tril within func
    if sparse:
        for i, nz in enumerate(nnz):
            if nz <= k: # too few neighbors
                D_ls[i, :] = D[i, :]
        return D_ls.tocsr()
    else:
        np.fill_diagonal(D_ls, self_value)
        return D_ls def shared_nearest_neighbors(D:np.ndarray, k:int=10, metric='distance',
                             n_jobs:int=1):
    """Transform distance matrix using shared nearest neighbors [1]_.

    SNN similarity is based on computing the overlap between the `k` nearest
    neighbors of two objects. SNN approaches try to symmetrize nearest neighbor
    relations using only rank and not distance information [2]_.

    Parameters
    ----------
    D : np.ndarray
        The ``n x n`` symmetric distance (similarity) matrix.

    k : int, optional (default: 10)
        Neighborhood radius: The `k` nearest neighbors are used to calculate SNN.

    metric : {'distance', 'similarity'}, optional (default: 'distance')
        Define, whether the matrix `D` is a distance or similarity matrix

    n_jobs : int, optional, default: 1
        Number of processes for parallel computations.

        - `1`: Don't use multiprocessing.
        - `-1`: Use all CPUs

    Returns
    -------
    D_snn : ndarray
        Secondary distance SNN matrix

    References
    ---------- 
    .. [1] R. Jarvis and E. A. Patrick, “Clustering using a similarity measure
           based on shared near neighbors,” IEEE Transactions on Computers,
           vol. 22, pp. 1025–1034, 1973.

    .. [2] Flexer, A., & Schnitzer, D. (2013). Can Shared Nearest Neighbors
           Reduce Hubness in High-Dimensional Spaces? 2013 IEEE 13th
           International Conference on Data Mining Workshops, 460–467.
           http://doi.org/10.1109/ICDMW.2013.101
    """
    io.check_distance_matrix_shape(D)
    io.check_valid_metric_parameter(metric)
    n = D.shape[0]
    if metric == 'distance':
        self_value = 0.
        sort_order = 1
        exclude = np.inf
        kth = k
    if metric == 'similarity':
        self_value = 1.
        sort_order = -1
        exclude = -np.inf
        kth = n - k
    distance = D.copy()
    np.fill_diagonal(distance, exclude)
    
    if n_jobs == -1:
        n_jobs = cpu_count()
    if n_jobs > 1:
        knn_ctype = RawArray(ctypes.c_bool, D.size)
        knn = np.frombuffer(knn_ctype, dtype=bool).reshape(D.shape)
        D_snn_ctype = RawArray(ctypes.c_double, D.size)
        D_snn = np.frombuffer(D_snn_ctype, dtype=np.float64).reshape(D.shape)
        with Pool(processes=n_jobs,
                  initializer=_snn_init,
                  initargs=(distance, knn, D_snn)) as pool:
            for _ in pool.imap(
                func=partial(_snn_my_hood, k=k, kth=kth, sort_order=sort_order),
                iterable=range(n)):
                pass
            for _ in pool.imap(
                func=partial(_snn_our_hood, k=k, metric=metric),
                iterable=range(n)):
                pass
    else:
        knn = np.zeros_like(distance, bool)
        # find nearest neighbors for each point
        for i in range(n):
            di = distance[i, :]
            nn = np.argpartition(di, kth=kth)[::sort_order]
            knn[i, nn[0:k]] = True
        D_snn = np.zeros_like(distance)
        for i in range(n):
            knn_i = knn[i, :]
            j_idx = slice(i+1, n)
    
            # using broadcasting
            Dij = np.sum(np.logical_and(knn_i, knn[j_idx, :]), 1)
            if metric == 'distance':
                D_snn[i, j_idx] = 1. - Dij / k
            else: # metric == 'similarity':
                D_snn[i, j_idx] = Dij / k

    D_snn += D_snn.T
    np.fill_diagonal(D_snn, self_value)
    return D_snn
from inspect import signature
from hub_toolbox import io
from hub_toolbox.centering import centering, weighted_centering, \
    localized_centering, dis_sim_global, dis_sim_local
from hub_toolbox.distances import cosine_distance
from hub_toolbox.global_scaling import mutual_proximity_empiric, \
    mutual_proximity_gammai, mutual_proximity_gaussi
from hub_toolbox.goodman_kruskal import goodman_kruskal_index
from hub_toolbox.hubness import hubness
from hub_toolbox.intrinsic_dimension import intrinsic_dimension
from hub_toolbox.knn_classification import score
from hub_toolbox.local_scaling import nicdm, local_scaling
from hub_toolbox.shared_neighbors import shared_nearest_neighbors
import numpy as np
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
This file is part of the HUB TOOLBOX available at
https://github.com/OFAI/hub-toolbox-python3/
The HUB TOOLBOX is licensed under the terms of the GNU GPLv3.

(c) 2011-2018, Dominik Schnitzer, Roman Feldbauer
Austrian Research Institute for Artificial Intelligence (OFAI)
Contact: <roman.feldbauer@ofai.at>
"""

__all__ = ['HubnessAnalysis']

CITATION = \
"""
R. Feldbauer, M. Leodolter, C. Plant and A. Flexer,
"Fast Approximate Hubness Reduction for Large High-Dimensional Data",
2018 IEEE International Conference on Big Knowledge (ICBK), Singapore, 2018,
pp. 358-367. doi: 10.1109/ICBK.2018.00055
(tech report available at http://www.ofai.at/cgi-bin/tr-online?number+2018-02)

or

R. Feldbauer, A. Flexer, "A comprehensive empirical comparison of hubness reduction in high-dimensional spaces"
Knowledge and Information Systems, 2018, https://doi.org/10.1007/s10115-018-1205-y
"""


def _primary_distance(D: np.ndarray, metric):
    """Return `D`, identical. (Dummy function.)"""
    return D


# New types of hubness reduction methods must be added here
SEC_DIST = {'mp': mutual_proximity_empiric,
            'mp_gaussi': mutual_proximity_gaussi,
            'mp_gammai': mutual_proximity_gammai,
            'ls': local_scaling,
            'nicdm': nicdm,

            'snn': shared_nearest_neighbors,