Add Multivariate Distance Covariance metric and the corresponding test of independence

dcor/__init__.py

@@ -36,10 +36,13 @@
 )
 from ._hypothesis import HypothesisTest as HypothesisTest
 from ._partial_dcor import (
-    partial_distance_correlation,
+    partial_distance_correlation as partial_distance_correlation,
     partial_distance_covariance as partial_distance_covariance,
 )
-from ._rowwise import RowwiseMode as RowwiseMode, rowwise as rowwise
+from ._rowwise import RowwiseMode as RowwiseMode, rowwise as rowwise, _initialize_rowwise_functions
 from ._utils import CompileMode as CompileMode
 
-__version__ = "0.8.dev0"
+__version__ = "0.6"
+
+# Initialize rowwise functions after all modules are loaded
+_initialize_rowwise_functions()

dcor/_dcor.py

@@ -41,6 +41,11 @@
     array_namespace,
     numpy_namespace,
 )
+##Additional module for Multivariate dcov test-----------------------------------------
+from scipy.special import gammaln
+import math
+from dcor._rowwise import rowwise
+##-------------------------------------------------------------------------------------
 
 Array = TypeVar("Array", bound=ArrayType)
 
@@ -1169,3 +1174,123 @@ def distance_correlation_af_inv(
             compile_mode=compile_mode,
         ),
     )
+
+
+
+
+
+def gamma_ratio(p):
+    """
+    Parameters
+    ----------
+    p : is the dimension of the data
+
+    Returns
+    -------
+    This function evaluates the gamma ratio, which is
+    required to calculate the constants C_p and C_q (in function u_dist_cov_sqr_mv())
+
+   """
+
+    return np.exp(gammaln((p + 1) / 2) - gammaln(p / 2))
+
+
+
+def rndm_projection(X, p):
+    """
+    Parameters
+    ----------
+    X : N x p, array of arrays
+    where, p: number of dimensions (p >= 1) and N: number of samples
+    p : number of dimensions (p >= 1)
+
+    Returns
+    -------
+    X_new : an array of size N
+    DESCRIPTION: Random projection of multivariate array
+    """
+
+    # X_std = multivariate_normal.rvs( np.zeros(p), np.identity(p), size = 1)
+    X_standard = np.random.standard_normal(p)
+
+    X_norm = np.linalg.norm(X_standard)
+    U_sphere = np.array(X_standard) / X_norm  # Normalize X_std
+
+    if p > 1:
+        X_new = U_sphere @ X.T
+    else:
+        X_new = U_sphere * X
+    return X_new
+
+
+def u_dist_cov_sqr_mv(X, Y, n_projs = 500, method = "mergesort"):
+    r"""
+    Numerically Efficient Multivariate Distance Covariance Based on
+         Random Projections, the computation complexity is
+    :math:`O(kn\log n)`, where k is the number of projections (n_projs)
+
+    For more details see,
+    :footcite:`b-dcov_random_projection`.
+
+
+    Parameters
+    ----------
+    X : N x p, array of arrays, where p > 1
+    Y : N x q, array of arrays, where q >= 1
+    where p and q: number of dimensions of variable X and Y, respectively and N: number of samples
+
+    n_projs : Number of projections (integer type), optional
+        DESCRIPTION. The default is 500.(paper suggests: n_projs < N/logN, larger n_projs provides better results)
+    method : fast computation method either "mergesort" or "avl", optional
+        DESCRIPTION. The default is "mergesort".
+
+    Returns
+    -------
+    omega_bar : Float type
+        DESCRIPTION: Produce fastly computed unbiased distance covariance between X and Y
+
+
+
+
+    Examples:
+        >>> import numpy as np
+        >>> import dcor
+        >>> from scipy.stats import multivariate_normal
+        >>> mean_vector = [2, 3, 5, 3, 2, 1]
+        >>> matrix_size = 6
+        >>> np.random.seed(123)  # in order to achieve reproducible results
+        >>> A = 0.5 * np.random.rand(matrix_size, matrix_size)
+        >>> B = np.dot(A, A.transpose())
+        >>> n_samples = 3000
+        >>> mv = multivariate_normal( mean = mean_vector, cov = B)
+        >>> X = mv.rvs(size = n_samples, random_state = 123)
+        >>> X1 = X.T[:4]
+        >>> X2 = X.T[4:]
+        >>> print(f"Computing fast distance covariance = {u_dist_cov_sqr_mv(X1.T, X2.T)}")
+    """
+
+    n_samples = np.shape(X)[0]
+    p = np.shape(X)[1]
+    if Y.T.ndim == 1:
+        q = 1
+    else:
+        q = np.shape(Y)[1]
+
+    sqrt_pi_value = math.sqrt(math.pi)
+    C_p = sqrt_pi_value * gamma_ratio(p)
+    C_q = sqrt_pi_value * gamma_ratio(q)
+
+
+    X_proj = np.empty(( n_projs, n_samples))
+    Y_proj = np.empty(( n_projs, n_samples))
+
+    for i in range(n_projs):
+        X_proj[i, :] = rndm_projection(X, p)
+        Y_proj[i, :] = rndm_projection(Y, q)
+        pass
+
+    omega_ = rowwise(u_distance_covariance_sqr,
+                     X_proj, Y_proj, rowwise_mode = method)
+    omega_bar = C_p * C_q * np.mean(omega_)
+
+    return omega_bar

dcor/_rowwise.py

@@ -17,9 +17,12 @@
 def _generate_rowwise_distance_covariance_sqr(unbiased):
     def rowwise_distance_covariance_sqr(
             x, y, exponent=1, *,
-            method=_dcor.DistanceCovarianceMethod.AUTO,
+            method=None,
             **kwargs):
 
+        if method is None:
+            method = _dcor.DistanceCovarianceMethod.AUTO
+
         if not _dcor._can_use_fast_algorithm(x[0], y[0],
                                              exponent=exponent):
             return NotImplemented
@@ -34,11 +37,13 @@ def rowwise_distance_covariance_sqr(
     return rowwise_distance_covariance_sqr
 
 
-_dcor.distance_covariance_sqr.rowwise_function = (
-    _generate_rowwise_distance_covariance_sqr(unbiased=False))
+def _initialize_rowwise_functions():
+    """Initialize rowwise functions after _dcor module is fully loaded."""
+    _dcor.distance_covariance_sqr.rowwise_function = (
+        _generate_rowwise_distance_covariance_sqr(unbiased=False))
 
-_dcor.u_distance_covariance_sqr.rowwise_function = (
-    _generate_rowwise_distance_covariance_sqr(unbiased=True))
+    _dcor.u_distance_covariance_sqr.rowwise_function = (
+        _generate_rowwise_distance_covariance_sqr(unbiased=True))
 
 
 def _rowwise_distance_covariance(*args, **kwargs):
@@ -50,7 +55,15 @@ def _rowwise_distance_covariance(*args, **kwargs):
     return _sqrt(res_covs)
 
 
-_dcor.distance_covariance.rowwise_function = _rowwise_distance_covariance
+    _dcor.distance_covariance.rowwise_function = _rowwise_distance_covariance
+
+    _dcor.distance_correlation_sqr.rowwise_function = (
+        _generate_rowwise_distance_correlation_sqr(unbiased=False))
+
+    _dcor.u_distance_correlation_sqr.rowwise_function = (
+        _generate_rowwise_distance_correlation_sqr(unbiased=True))
+
+    _dcor.distance_correlation.rowwise_function = _rowwise_distance_correlation
 
 
 def _generate_rowwise_distance_correlation_sqr(unbiased):
@@ -82,13 +95,6 @@ def rowwise_distance_correlation_sqr(x, y, **kwargs):
     return rowwise_distance_correlation_sqr
 
 
-_dcor.distance_correlation_sqr.rowwise_function = (
-    _generate_rowwise_distance_correlation_sqr(unbiased=False))
-
-_dcor.u_distance_correlation_sqr.rowwise_function = (
-    _generate_rowwise_distance_correlation_sqr(unbiased=True))
-
-
 def _rowwise_distance_correlation(*args, **kwargs):
 
     res_corrs = _dcor.distance_correlation_sqr.rowwise_function(
@@ -99,9 +105,6 @@ def _rowwise_distance_correlation(*args, **kwargs):
     return _sqrt(res_corrs)
 
 
-_dcor.distance_correlation.rowwise_function = _rowwise_distance_correlation
-
-
 def rowwise(
     function: Callable[..., Array],
     x: Array,

dcor/distances.py

@@ -15,6 +15,8 @@
 
 from dcor._utils import ArrayType, _sqrt, _transform_to_2d, array_namespace
 
+from numba import njit, prange # Additional repo for dist_sum
+
 from ._utils import _can_be_numpy_double
 
 Array = TypeVar("Array", bound=ArrayType)
@@ -159,3 +161,34 @@ def pairwise_distances(
 
     x, y = _transform_to_2d(x, y)
     return _cdist(x, y, exponent=exponent)
+
+
+
+
+
+@njit(fastmath=True, parallel=True, cache=True)
+def dist_sum(X): 
+    """
+    Parameters
+    ----------
+    X : 1D array.
+
+    Returns
+    -------
+    res : sum of distinct Euclidean distances corresponding to the elements of X.
+
+    Note: To implement numba, one needs to consider "numpy==1.25" or less. 
+    """
+    res = 0
+    for i in prange(len(X)):
+        for j in prange(len(X)):
+            if i < j:
+                res += np.abs(X[i] - X[j])
+                pass
+            pass
+        pass
+    return res
+
+
+
+