ENH: random.standard_t impl (#515)

* ENH: random.standard_t impl

Author: samir-nasibli

dpnp/backend/include/dpnp_iface.hpp

@@ -948,6 +948,20 @@ INP_DLLEXPORT void dpnp_rng_standard_gamma_c(void* result, _DataType shape, size
 template <typename _DataType>
 INP_DLLEXPORT void dpnp_rng_standard_normal_c(void* result, size_t size);
 
+/**
+ * @ingroup BACKEND_API
+ * @brief math library implementation of random number generator (standard Student's t distribution)
+ *
+ * @param [in]  size   Number of elements in `result` arrays.
+ *
+ * @param [in]  df     Degrees of freedom.
+ *
+ * @param [out] result Output array.
+ *
+ */
+template <typename _DataType>
+INP_DLLEXPORT void dpnp_rng_standard_t_c(void* result, _DataType df, size_t size);
+
 /**
  * @ingroup BACKEND_API
  * @brief math library implementation of random number generator (uniform distribution)

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -144,6 +144,7 @@ enum class DPNPFuncName : size_t
     DPNP_FN_RNG_STANDARD_EXPONENTIAL, /**< Used in numpy.random.standard_exponential() implementation  */
     DPNP_FN_RNG_STANDARD_GAMMA,       /**< Used in numpy.random.standard_gamma() implementation  */
     DPNP_FN_RNG_STANDARD_NORMAL,      /**< Used in numpy.random.standard_normal() implementation  */
+    DPNP_FN_RNG_STANDARD_T,           /**< Used in numpy.random.standard_t() implementation  */
     DPNP_FN_RNG_UNIFORM,              /**< Used in numpy.random.uniform() implementation  */
     DPNP_FN_RNG_WEIBULL,              /**< Used in numpy.random.weibull() implementation  */
     DPNP_FN_SIGN,                     /**< Used in numpy.sign() implementation  */

dpnp/backend/kernels/dpnp_krnl_random.cpp

@@ -564,6 +564,74 @@ void dpnp_rng_standard_normal_c(void* result, size_t size)
     dpnp_rng_normal_c(result, mean, stddev, size);
 }
 
+template <typename _DataType>
+void dpnp_rng_standard_t_c(void* result, _DataType df, size_t size)
+{
+    if (!size)
+    {
+        return;
+    }
+    cl::sycl::vector_class<cl::sycl::event> no_deps;
+
+    _DataType* result1 = reinterpret_cast<_DataType*>(result);
+    const _DataType d_zero = 0.0, d_one = 1.0;
+    _DataType shape = df/2;
+    _DataType *sn = nullptr;
+
+    if (dpnp_queue_is_cpu_c())
+    {
+        mkl_rng::gamma<_DataType> gamma_distribution(shape, d_zero, 1.0/shape);
+        auto event_out = mkl_rng::generate(gamma_distribution, DPNP_RNG_ENGINE,
+            size, result1);
+        event_out.wait();
+        event_out = mkl_vm::invsqrt(DPNP_QUEUE, size, result1, result1, no_deps,
+            mkl_vm::mode::ha);
+        event_out.wait();
+
+        sn = reinterpret_cast<_DataType*>(dpnp_memory_alloc_c(size * sizeof(_DataType)));
+        if (sn == nullptr)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_standard_t_c() failed.");
+        }
+
+        mkl_rng::gaussian<_DataType> gaussian_distribution(d_zero, d_one);
+        event_out = mkl_rng::generate(gaussian_distribution, DPNP_RNG_ENGINE, size, sn);
+        event_out.wait();
+
+        event_out = mkl_vm::mul(DPNP_QUEUE, size, result1, sn, result1, no_deps,
+            mkl_vm::mode::ha);
+        dpnp_memory_free_c(sn);
+        event_out.wait();
+    }
+    else
+    {
+        int errcode = vdRngGamma(VSL_RNG_METHOD_GAMMA_GNORM_ACCURATE, get_rng_stream(),
+            size, result1, shape, d_zero, 1.0/shape);
+
+        if (errcode != VSL_STATUS_OK)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_standard_t_c() failed.");
+        }
+
+        vmdInvSqrt(size, result1, result1, VML_HA);
+
+        sn = (_DataType *) mkl_malloc(size * sizeof(_DataType), 64);
+        if (sn == nullptr)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_standard_t_c() failed.");
+        }
+
+        errcode = vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, get_rng_stream(), size, sn,
+            d_zero, d_one);
+        if (errcode != VSL_STATUS_OK)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_standard_t_c() failed.");
+        }
+        vmdMul(size, result1, sn, result1, VML_HA);
+        mkl_free(sn);
+    }
+}
+
 template <typename _DataType>
 void dpnp_rng_uniform_c(void* result, long low, long high, size_t size)
 {
@@ -658,6 +726,8 @@ void func_map_init_random(func_map_t& fmap)
 
     fmap[DPNPFuncName::DPNP_FN_RNG_STANDARD_NORMAL][eft_DBL][eft_DBL] = {eft_DBL,
                                                                          (void*)dpnp_rng_standard_normal_c<double>};
+    fmap[DPNPFuncName::DPNP_FN_RNG_STANDARD_T][eft_DBL][eft_DBL] = {eft_DBL,
+                                                                         (void*)dpnp_rng_standard_t_c<double>};
 
     fmap[DPNPFuncName::DPNP_FN_RNG_UNIFORM][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_uniform_c<double>};
     fmap[DPNPFuncName::DPNP_FN_RNG_UNIFORM][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_rng_uniform_c<float>};

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -117,6 +117,7 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName":  # need this na
         DPNP_FN_RNG_STANDARD_EXPONENTIAL
         DPNP_FN_RNG_STANDARD_GAMMA
         DPNP_FN_RNG_STANDARD_NORMAL
+        DPNP_FN_RNG_STANDARD_T
         DPNP_FN_RNG_UNIFORM
         DPNP_FN_RNG_WEIBULL
         DPNP_FN_SIGN

dpnp/random/dpnp_algo_random.pyx

@@ -67,6 +67,7 @@ __all__ = [
     "dpnp_standard_exponential",
     "dpnp_standard_gamma",
     "dpnp_standard_normal",
+    "dpnp_standard_t",
     "dpnp_uniform",
     "dpnp_weibull"
 ]
@@ -100,6 +101,7 @@ ctypedef void(*fptr_dpnp_rng_standard_cauchy_c_1out_t)(void *, size_t) except +
 ctypedef void(*fptr_dpnp_rng_standard_exponential_c_1out_t)(void *, size_t) except +
 ctypedef void(*fptr_dpnp_rng_standard_gamma_c_1out_t)(void *, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_standard_normal_c_1out_t)(void *, size_t) except +
+ctypedef void(*fptr_dpnp_rng_standard_t_c_1out_t)(void *, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_uniform_c_1out_t)(void *, long, long, size_t) except +
 ctypedef void(*fptr_dpnp_rng_weibull_c_1out_t)(void *, double, size_t) except +
 
@@ -867,6 +869,30 @@ cpdef dparray dpnp_standard_normal(size):
 
     return result
 
+cpdef dparray dpnp_standard_t(double df, size):
+    """
+    Returns an array populated with samples from standard t distribution.
+    `dpnp_standard_t` generates a matrix filled with random floats sampled from a
+    univariate standard t distribution for a given number of degrees of freedom.
+
+    """
+
+    # convert string type names (dparray.dtype) to C enum DPNPFuncType
+    cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(numpy.float64)
+
+    # get the FPTR data structure
+    cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_RNG_STANDARD_T, param1_type, param1_type)
+
+    result_type = dpnp_DPNPFuncType_to_dtype( < size_t > kernel_data.return_type)
+    # ceate result array with type given by FPTR data
+    cdef dparray result = dparray(size, dtype=result_type)
+
+    cdef fptr_dpnp_rng_standard_t_c_1out_t func = <fptr_dpnp_rng_standard_t_c_1out_t > kernel_data.ptr
+    # call FPTR function
+    func(result.get_data(), df, result.size)
+
+    return result
+
 
 cpdef dparray dpnp_uniform(long low, long high, size, dtype):
     """

dpnp/random/dpnp_iface_random.py

@@ -1358,20 +1358,38 @@ def standard_normal(size=None):
 
 
 def standard_t(df, size=None):
-    """Power distribution.
+    """Standard Student’s t distribution.
 
-    Draw samples from a standard Student’s t distribution with df degrees
-    of freedom.
+    Draw samples from a standard Student’s t distribution with
+    df degrees of freedom.
 
     For full documentation refer to :obj:`numpy.random.standard_t`.
 
-    Notes
-    -----
-    The function uses `numpy.random.standard_t` on the backend and
-    will be executed on fallback backend.
+    Limitations
+    -----------
+    Parameter ``df`` is supported as a scalar.
+    Otherwise, :obj:`numpy.random.standard_t(df, size)` samples
+    are drawn.
+    Output array data type is :obj:`dpnp.float64`.
+
+    Examples
+    --------
+    Draw samples from the distribution:
+    >>> df = 2.
+    >>> s = dpnp.random.standard_t(df, 1000000)
 
     """
 
+    if not use_origin_backend(df) and dpnp_queue_is_cpu():
+        # TODO:
+        # array_like of floats for `df`
+        if not dpnp.isscalar(df):
+            pass
+        elif df <= 0:
+            pass
+        else:
+            return dpnp_standard_t(df, size)
+    print("here")
     return call_origin(numpy.random.standard_t, df, size)
 
 

tests/test_random.py

@@ -671,6 +671,23 @@ def test_seed(self):
         self.check_seed('standard_normal', {})
 
 
+class TestDistributionsStandardT(TestDistribution):
+
+    def test_moments(self):
+        df = 300.0
+        expected_mean = 0.0
+        expected_var = df / (df - 2)
+        self.check_moments('standard_t', expected_mean,
+                           expected_var, {'df': df})
+
+    def test_invalid_args(self):
+        df = 0.0   # positive `df` is expected
+        self.check_invalid_args('standard_t', {'df': df})
+
+    def test_seed(self):
+        self.check_seed('standard_t', {'df': 10.0})
+
+
 class TestDistributionsUniform(TestDistribution):
 
     def test_extreme_value(self):

tests_external/skipped_tests_numpy.tbl

@@ -1501,6 +1501,7 @@ tests/test_random.py::TestRandomDist::test_standard_exponential
 tests/test_random.py::TestRandomDist::test_standard_gamma
 tests/test_random.py::TestRandomDist::test_standard_gamma_0
 tests/test_random.py::TestRandomDist::test_standard_normal
+tests/test_random.py::TestRandomDist::test_standard_t
 tests/test_random.py::TestRandomDist::test_uniform
 tests/test_random.py::TestRandomDist::test_uniform_range_bounds
 tests/test_random.py::TestRandomDist::test_weibull
@@ -1595,6 +1596,7 @@ tests/test_randomstate.py::TestRandomDist::test_standard_exponential
 tests/test_randomstate.py::TestRandomDist::test_standard_gamma
 tests/test_randomstate.py::TestRandomDist::test_standard_gamma_0
 tests/test_randomstate.py::TestRandomDist::test_standard_normal
+tests/test_randomstate.py::TestRandomDist::test_standard_t
 tests/test_randomstate.py::TestRandomDist::test_tomaxint
 tests/test_randomstate.py::TestRandomDist::test_uniform
 tests/test_randomstate.py::TestRandomDist::test_uniform_range_bounds