Implementation of corrcoef

Author: AlexanderKalistratov

dpnp/dpnp_iface_statistics.py

@@ -37,6 +37,8 @@
 
 """
 
+import warnings
+
 import dpctl.tensor as dpt
 import numpy
 from dpctl.tensor._numpy_helper import (
@@ -65,6 +67,7 @@
     "amin",
     "average",
     "bincount",
+    "corrcoef",
     "correlate",
     "cov",
     "max",
@@ -403,6 +406,69 @@ def correlate(x1, x2, mode="valid"):
     return call_origin(numpy.correlate, x1, x2, mode=mode)
 
 
+def corrcoef(x, y=None, rowvar=True, bias=None, ddof=None, *, dtype=None):
+    """
+    Return Pearson product-moment correlation coefficients.
+
+    For full documentation refer to :obj:`numpy.corrcoef`.
+
+    Parameters
+    ----------
+    x : {dpnp.ndarray, usm_ndarray}
+        A 1-D or 2-D array containing multiple variables and observations.
+        Each row of `x` represents a variable, and each column a single
+        observation of all those variables. Also see `rowvar` below.
+    y : {dpnp.ndarray, usm_ndarray}, optional
+        An additional set of variables and observations. `y` has the same
+        shape as `x`.
+    rowvar : {bool}, optional
+        If `rowvar` is True (default), then each row represents a
+        variable, with observations in the columns. Otherwise, the relationship
+        is transposed: each column represents a variable, while the rows
+        contain observations.
+    bias : {None}, optional
+        Has no effect, do not use.
+    ddof : {None}, optional
+        Has no effect, do not use.
+    dtype : data-type, optional
+        Data-type of the result.
+
+    Returns
+    -------
+    R : dpnp.ndarray
+        The correlation coefficient matrix of the variables.
+
+    See Also
+    --------
+    :obj:`dpnp.cov` : Covariance matrix.
+    """
+    if bias is not None or ddof is not None:
+        warnings.warn(
+            "bias and ddof have no effect and are deprecated",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+
+    out = dpnp.cov(x, y, rowvar, dtype=dtype)
+    try:
+        d = dpnp.diag(out)
+    except ValueError:
+        return out / out
+
+    stddev = dpnp.sqrt(d.real)
+    out /= stddev[:, None]
+    out /= stddev[None, :]
+
+    # Clip real and imaginary parts to [-1, 1].  This does not guarantee
+    # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
+    # excessive work.
+    dpnp.clip(out.real, -1, 1, out=out.real)
+    if dpnp.iscomplexobj(out):
+        dpnp.clip(out.imag, -1, 1, out=out.imag)
+
+    return out
+
+
 def cov(
     m,
     y=None,

tests/skipped_tests_gpu.tbl

@@ -322,11 +322,6 @@ tests/third_party/cupy/random_tests/test_sample.py::TestRandomIntegers2::test_bo
 tests/third_party/cupy/random_tests/test_sample.py::TestRandomIntegers2::test_goodness_of_fit
 tests/third_party/cupy/random_tests/test_sample.py::TestRandomIntegers2::test_goodness_of_fit_2
 
-tests/third_party/cupy/statistics_tests/test_correlation.py::TestCorrcoef::test_corrcoef
-tests/third_party/cupy/statistics_tests/test_correlation.py::TestCorrcoef::test_corrcoef_diag_exception
-tests/third_party/cupy/statistics_tests/test_correlation.py::TestCorrcoef::test_corrcoef_rowvar
-tests/third_party/cupy/statistics_tests/test_correlation.py::TestCorrcoef::test_corrcoef_y
-
 tests/third_party/cupy/statistics_tests/test_order.py::TestOrder::test_percentile_defaults[linear]
 tests/third_party/cupy/statistics_tests/test_order.py::TestOrder::test_percentile_defaults[lower]
 tests/third_party/cupy/statistics_tests/test_order.py::TestOrder::test_percentile_defaults[higher]

tests/test_statistics.py

@@ -573,6 +573,51 @@ def test_std_error(self):
             dpnp.std(ia, ddof="1")
 
 
+class TestCorrcoef:
+    @pytest.mark.usefixtures("suppress_divide_invalid_numpy_warnings")
+    @pytest.mark.parametrize("dtype", get_all_dtypes())
+    @pytest.mark.parametrize("rowvar", [True, False])
+    def test_corrcoef(self, dtype, rowvar):
+        dp_array = dpnp.array([[0, 1, 2], [3, 4, 0]], dtype=dtype)
+        np_array = dpnp.asnumpy(dp_array)
+
+        expected = numpy.corrcoef(np_array, rowvar=rowvar)
+        result = dpnp.corrcoef(dp_array, rowvar=rowvar)
+
+        assert_dtype_allclose(result, expected)
+
+    @pytest.mark.usefixtures("suppress_divide_invalid_numpy_warnings")
+    @pytest.mark.parametrize("shape", [(2, 0), (0, 2)])
+    def test_corrcoef_empty(self, shape):
+        dp_array = dpnp.empty(shape, dtype=dpnp.int64)
+        np_array = dpnp.asnumpy(dp_array)
+
+        result = dpnp.corrcoef(dp_array)
+        expected = numpy.corrcoef(np_array)
+        assert_dtype_allclose(result, expected)
+
+    @pytest.mark.usefixtures("suppress_complex_warning")
+    @pytest.mark.parametrize("dt_in", get_all_dtypes(no_bool=True))
+    @pytest.mark.parametrize("dt_out", get_float_complex_dtypes())
+    def test_corrcoef_dtype(self, dt_in, dt_out):
+        dp_array = dpnp.array([[0, 1, 2], [3, 4, 0]], dtype=dt_in)
+        np_array = dpnp.asnumpy(dp_array)
+
+        expected = numpy.corrcoef(np_array, dtype=dt_out)
+        result = dpnp.corrcoef(dp_array, dtype=dt_out)
+        assert expected.dtype == result.dtype
+        assert_allclose(result, expected, rtol=1e-6)
+
+    @pytest.mark.usefixtures("suppress_divide_invalid_numpy_warnings")
+    def test_corrcoef_scalar(self):
+        dp_array = dpnp.array(5)
+        np_array = dpnp.asnumpy(dp_array)
+
+        result = dpnp.corrcoef(dp_array)
+        expected = numpy.corrcoef(np_array)
+        assert_dtype_allclose(result, expected)
+
+
 @pytest.mark.usefixtures("allow_fall_back_on_numpy")
 class TestBincount:
     @pytest.mark.parametrize(

tests/third_party/cupy/statistics_tests/test_correlation.py

@@ -12,26 +12,26 @@
 
 class TestCorrcoef(unittest.TestCase):
     @testing.for_all_dtypes()
-    @testing.numpy_cupy_allclose()
+    @testing.numpy_cupy_allclose(type_check=False)
     def test_corrcoef(self, xp, dtype):
         a = testing.shaped_arange((2, 3), xp, dtype)
         return xp.corrcoef(a)
 
     @testing.for_all_dtypes()
-    @testing.numpy_cupy_allclose()
+    @testing.numpy_cupy_allclose(type_check=False)
     def test_corrcoef_diag_exception(self, xp, dtype):
         a = testing.shaped_arange((1, 3), xp, dtype)
         return xp.corrcoef(a)
 
     @testing.for_all_dtypes()
-    @testing.numpy_cupy_allclose()
+    @testing.numpy_cupy_allclose(type_check=False)
     def test_corrcoef_y(self, xp, dtype):
         a = testing.shaped_arange((2, 3), xp, dtype)
         y = testing.shaped_arange((2, 3), xp, dtype)
         return xp.corrcoef(a, y=y)
 
     @testing.for_all_dtypes()
-    @testing.numpy_cupy_allclose()
+    @testing.numpy_cupy_allclose(type_check=False)
     def test_corrcoef_rowvar(self, xp, dtype):
         a = testing.shaped_arange((2, 3), xp, dtype)
         y = testing.shaped_arange((2, 3), xp, dtype)