Fix IndexError in MGCX and Time Series Tests for Zero Variance Slices

.circleci/config.yml

@@ -2,7 +2,6 @@ version: 2.1
 
 orbs:
   python: circleci/python@1.2
-  codecov: codecov/codecov@3.1.1
 
 jobs:
   build-and-test:
@@ -25,8 +24,11 @@ jobs:
           path: test-results
       - store_artifacts:
           path: test-results
-      - codecov/upload:
-          file: "coverage.xml"
+      - run:
+          name: Upload to Codecov
+          command: |
+            pip install codecov-cli
+            codecovcli upload-coverage --file coverage.xml
 
 workflows:
   run-tests:

codecov.yml

@@ -19,9 +19,9 @@ coverage:
 # --------------
 # which folders/files to ignore
 ignore:
-  - *tests*
-  - *__init__.py
-  - *hyppo/_utils.py
+  - "*tests*"
+  - "*__init__.py"
+  - "*hyppo/_utils.py"
 
 # Pull request comments:
 # ----------------------

hyppo/conditional/tests/test_FCIT.py

@@ -24,8 +24,7 @@ def test_linear_oned(self, n, obs_stat, obs_pvalue):
 
     @pytest.mark.parametrize(
         "dim, n, obs_stat, obs_pvalue",
-        # 0.56139, -0.16024
-        [(1, 100000, -0.06757, 0.52599), (2, 100000, -4.59882, 0.99876)],
+        [(1, 100000, -0.16024, 0.56139), (2, 100000, -4.59882, 0.99876)],
     )
     def test_null(self, dim, n, obs_stat, obs_pvalue):
         np.random.seed(12)
@@ -56,9 +55,8 @@ def test_null(self, dim, n, obs_stat, obs_pvalue):
     @pytest.mark.parametrize(
         "dim, n, obs_stat, obs_pvalue",
         [
-            #89.271754, 161.35165
-            (1, 100000, 89.184784, 2.91447597e-12),
-            (2, 100000, 161.35105, 4.63412957e-14),
+            (1, 100000, 89.271754, 2.91447597e-12),
+            (2, 100000, 161.35165, 4.63412957e-14),
         ],
     )
     def test_alternative(self, dim, n, obs_stat, obs_pvalue):

hyppo/discrim/tests/test_discrim_one_samp.py

@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_almost_equal, assert_raises, assert_warns
+from numpy.testing import assert_almost_equal, assert_raises
 
 from .. import DiscrimOneSample
 
@@ -32,6 +32,13 @@ def test_diff_one(self):
         assert_almost_equal(stat, obs_stat, decimal=3)
         assert_almost_equal(p, obs_p, decimal=3)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        x = np.concatenate((np.zeros((50, 2)), np.ones((50, 2))), axis=0).astype(dtype)
+        y = np.concatenate((np.zeros(50), np.ones(50)), axis=0).astype(dtype)
+        stat, p, _ = DiscrimOneSample().test(x, y, reps=10)
+        assert_almost_equal(stat, 1.0, decimal=3)
+
 
 class TestOneSampleWarn:
     """Tests errors and warnings derived from one sample test."""

hyppo/independence/dcorr.py

@@ -1,7 +1,8 @@
 import numpy as np
 from numba import jit
 
-from ..tools import check_perm_blocks_dim, chi2_approx, compute_dist
+from ..tools import check_perm_blocks_dim, chi2_approx, compute_dist, convert_xy_float64
+
 from ._utils import _CheckInputs
 from .base import IndependenceTest, IndependenceTestOutput
 
@@ -138,6 +139,11 @@ def statistic(self, x, y):
         stat : float
             The computed Dcorr statistic.
         """
+        # Only convert dtype when inputs are raw data, not precomputed distance
+        # matrices. Distance matrices produced by compute_dist are already float64;
+        # converting them on every permutation call causes large unnecessary copies.
+        if not (self.is_distance or self.is_fast):
+            x, y = convert_xy_float64(x, y)
         distx = x
         disty = y
 
@@ -302,15 +308,18 @@ def _fast_1d_dcov(x, y, bias=False):  # pragma: no cover
 
     # sort inputs
     x_orig = x.ravel()
+    y_orig = y.ravel()
     x = np.sort(x_orig)
-    y = y[np.argsort(x_orig)]
+    y = y_orig[np.argsort(x_orig)]
     x = x.reshape(-1, 1)  # for numba
+    y_col = y.reshape(-1, 1)
 
     # cumulative sum
     si = _cpu_cumsum(x)
     ax = (np.arange(-(n - 2), n + 1, 2) * x.ravel()).reshape(-1, 1) + (si[-1] - 2 * si)
 
-    v = np.hstack((x, y, x * y))
+    v = np.hstack((x, y_col, x * y_col))
+
     nw = v.shape[1]
 
     idx = np.vstack((np.arange(n), np.zeros(n))).astype(np.int64).T
@@ -370,7 +379,8 @@ def _fast_1d_dcov(x, y, bias=False):  # pragma: no cover
     c4 = np.sum(iv3.T @ x)
     d = 4 * ((c1 + c2) - (c3 + c4)) - 2 * covterm
 
-    y_sorted = y[idx[n::-1, r], :]
+    y_sorted = y_col[idx[n::-1, r], :]
+
     si = _cpu_cumsum(y_sorted)
     by = np.zeros((n, 1))
     by[idx[::-1, r]] = (np.arange(-(n - 2), n + 1, 2) * y_sorted.ravel()).reshape(

hyppo/independence/mgc.py

@@ -4,7 +4,8 @@
 import numpy as np
 from scipy.stats import multiscale_graphcorr
 
-from ..tools import compute_dist
+from ..tools import compute_dist, convert_xy_float64
+
 from ._utils import _CheckInputs
 from .base import IndependenceTest
 
@@ -148,6 +149,11 @@ def statistic(self, x, y):
         stat : float
             The computed MGC statistic.
         """
+        # Only convert dtype when inputs are raw data, not precomputed distance
+        # matrices. Distance matrices from compute_dist are already float64;
+        # converting them on every permutation call causes large unnecessary copies.
+        if not self.is_distance:
+            x, y = convert_xy_float64(x, y)
         distx = x
         disty = y
 
@@ -158,9 +164,13 @@ def statistic(self, x, y):
 
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore")
-            mgc = multiscale_graphcorr(distx, disty, compute_distance=None, reps=0)
+            try:
+                mgc = multiscale_graphcorr(distx, disty, compute_distance=None, reps=0)
+                stat = mgc.stat
+            except IndexError:
+                # scipy mgc throws IndexError if disty is all zeros (variance 0)
+                stat = 0.0
 
-        stat = mgc.stat
         self.stat = stat
 
         return stat

hyppo/independence/tests/test_dcorr.py

@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_almost_equal, assert_raises, assert_warns
+from numpy.testing import assert_almost_equal
 
 from ...tools import linear, power
 from .. import Dcorr
@@ -30,12 +30,21 @@ def test_rep(self, n):
         assert pvalue1 == pvalue2
 
     def test_dcorr_sqrt_bug(self):
-        x = np.array([1,2,3,4,5])
-        y = np.array([1,2,9,4,4])
+        x = np.array([1, 2, 3, 4, 5])
+        y = np.array([1, 2, 9, 4, 4])
         stat = Dcorr(bias=True).test(x, y, reps=0)[0]
 
         assert_almost_equal(stat, 0.762676242417, decimal=2)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        np.random.seed(123456789)
+        x, y = linear(100, 1)
+        x_cast = (x * 100).astype(dtype)
+        y_cast = (y * 100).astype(dtype)
+        stat = Dcorr().statistic(x_cast, y_cast)
+        assert_almost_equal(stat, 1.0, decimal=2)
+
 
 class TestDcorrTypeIError:
     def test_oned(self):

hyppo/independence/tests/test_mgc.py

@@ -3,11 +3,10 @@
 from numpy.testing import (
     assert_almost_equal,
     assert_approx_equal,
-    assert_equal,
     assert_warns,
 )
 
-from ...tools import linear, multimodal_independence, power, spiral
+from ...tools import linear, power, spiral
 from .. import MGC
 
 
@@ -68,6 +67,15 @@ def test_rep(self, sim, obs_stat, obs_pvalue):
         assert stat1 == stat2
         assert pvalue1 == pvalue2
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        np.random.seed(12345678)
+        x, y = linear(n=100, p=1)
+        x_cast = (x * 100).astype(dtype)
+        y_cast = (y * 100).astype(dtype)
+        stat = MGC().statistic(x_cast, y_cast)
+        assert_approx_equal(stat, 0.97, significant=1)
+
 
 class TestMGCTypeIError:
     def test_oned(self):

hyppo/kgof/kernel.py

@@ -86,7 +86,7 @@ def pair_gradY_X(self, X, Y):
 
     @abstractmethod
     def pair_gradXY_sum(self, X, Y):
-        """
+        r"""
         Compute \sum_{i=1}^d \frac{\partial^2 k(X, Y)}{\partial x_i \partial y_i}
         evaluated at each x_i in X, and y_i in Y.
         X: n x d numpy array.
@@ -125,7 +125,7 @@ def gradY_X(self, X, Y, dim):
 
     @abstractmethod
     def gradXY_sum(self, X, Y):
-        """
+        r"""
         Compute \sum_{i=1}^d \frac{\partial^2 k(x, Y)}{\partial x_i \partial y_i}
         evaluated at each x_i in X, and y_i in Y.
         X: nx x d numpy array.
@@ -260,7 +260,7 @@ def gradXY_sum(self, X, Y):
         return G
 
     def pair_gradXY_sum(self, X, Y):
-        """
+        r"""
         Compute \sum_{i=1}^d \frac{\partial^2 k(X, Y)}{\partial x_i \partial y_i}
         evaluated at each x_i in X, and y_i in Y.
         X: n x d numpy array.

hyppo/ksample/mean_embedding.py

@@ -166,8 +166,8 @@ def mean_embed_distance(difference, num_randfeatures):
     mu = np.mean(difference, 0)
 
     if num_randfeatures == 1:
-        stat = float(num_samples * mu**2) / float(sigma)
+        stat = float(num_samples * mu.item() ** 2 / sigma.item())
     else:
-        stat = num_samples * mu.dot(np.linalg.solve(sigma, np.transpose(mu)))
+        stat = float(num_samples * mu.dot(np.linalg.solve(sigma, np.transpose(mu))))
 
     return stat

hyppo/ksample/tests/test_ksamp.py

@@ -45,6 +45,14 @@ def test_rep(self, n, obs_stat, obs_pvalue, indep_test):
         assert stat == stat2
         assert pvalue == pvalue2
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        np.random.seed(123456789)
+        inputs = rot_ksamp("linear", 100, 1, k=2)
+        inputs_cast = [mat.astype(dtype) for mat in inputs]
+        stat, _ = KSample("Dcorr").test(*inputs_cast, reps=0)
+        assert_almost_equal(stat, 0.045646974150778084, decimal=1)
+
 
 class TestKSampleErrorWarn:
     """Tests errors and warnings derived from MGC."""

hyppo/time_series/_utils.py

@@ -89,16 +89,23 @@ def compute_stat(x, y, indep_test, compute_distance, max_lag, **kwargs):
 
     # calculate dep_lag when max_lag is 0
     dep_lag = []
-    indep_test = indep_test(compute_distance=compute_distance, **kwargs)
-    indep_test_stat = indep_test.statistic(x, y)
+    indep_test = indep_test(compute_distance=None, **kwargs)
+    indep_test_stat = indep_test.statistic(distx, disty)
     dep_lag.append(indep_test_stat)
 
     # loop over time points and find max test statistic
     n = distx.shape[0]
     for j in range(1, max_lag + 1):
         slice_distx = distx[j:n, j:n]
         slice_disty = disty[0 : (n - j), 0 : (n - j)]
-        stat = indep_test.statistic(slice_distx, slice_disty)
+        try:
+            stat = indep_test.statistic(slice_distx, slice_disty)
+        except IndexError as e:
+            print("ERROR on slice_distx:")
+            print(slice_distx)
+            print("ERROR on slice_disty:")
+            print(slice_disty)
+            raise e
         dep_lag.append((n - j) * stat / n)
 
     # calculate optimal lag and test statistic
@@ -121,6 +128,10 @@ def compute_scale_at_lag(x, y, opt_lag, compute_distance, **kwargs):
     mgc = MGC()
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore")
-        opt_scale = mgc.test(slice_distx, slice_disty, reps=0)[2]["opt_scale"]
+        try:
+            opt_scale = mgc.test(slice_distx, slice_disty, reps=0)[2]["opt_scale"]
+        except ValueError:
+            # test cannot be run due to zero variance
+            opt_scale = [np.nan, np.nan]
 
     return opt_scale

hyppo/tools/common.py

@@ -171,20 +171,26 @@ def check_ndarray_or_dataframe(data, col_id):
 
 def convert_xy_float64(x, y):
     """Convert x or y to np.float64 (if not already done)"""
-    # convert x and y to floats
-    x = np.asarray(x).astype(np.float64)
-    y = np.asarray(y).astype(np.float64)
-
+    x = np.asarray(x)
+    y = np.asarray(y)
+    if x.dtype != np.float64:
+        x = x.astype(np.float64)
+    if y.dtype != np.float64:
+        y = y.astype(np.float64)
     return x, y
 
 
 def convert_xyz_float64(x, y, z):
     """Convert x or y or z to np.float64 (if not already done)"""
-    # convert x and y to floats
-    x = np.asarray(x).astype(np.float64)
-    y = np.asarray(y).astype(np.float64)
-    z = np.asarray(z).astype(np.float64)
-
+    x = np.asarray(x)
+    y = np.asarray(y)
+    z = np.asarray(z)
+    if x.dtype != np.float64:
+        x = x.astype(np.float64)
+    if y.dtype != np.float64:
+        y = y.astype(np.float64)
+    if z.dtype != np.float64:
+        z = z.astype(np.float64)
     return x, y, z