ENH: impl random.f (#528)

Author: samir-nasibli

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -126,6 +126,7 @@ enum class DPNPFuncName : size_t
     DPNP_FN_RNG_BINOMIAL,             /**< Used in numpy.random.binomial() implementation  */
     DPNP_FN_RNG_CHISQUARE,            /**< Used in numpy.random.chisquare() implementation  */
     DPNP_FN_RNG_EXPONENTIAL,          /**< Used in numpy.random.exponential() implementation  */
+    DPNP_FN_RNG_F,                    /**< Used in numpy.random.f() implementation  */
     DPNP_FN_RNG_GAMMA,                /**< Used in numpy.random.gamma() implementation  */
     DPNP_FN_RNG_GAUSSIAN,             /**< Used in numpy.random.randn() implementation  */
     DPNP_FN_RNG_GEOMETRIC,            /**< Used in numpy.random.geometric() implementation  */

dpnp/backend/include/dpnp_iface_random.hpp

@@ -98,6 +98,18 @@ INP_DLLEXPORT void dpnp_rng_chi_square_c(void* result, int df, size_t size);
 template <typename _DataType>
 INP_DLLEXPORT void dpnp_rng_exponential_c(void* result, _DataType beta, size_t size);
 
+/**
+ * @ingroup BACKEND_RANDOM_API
+ * @brief math library implementation of random number generator (F distribution)
+ *
+ * @param [in]  size   Number of elements in `result` arrays.
+ * @param [in]  df_num  Degrees of freedom in numerator.
+ * @param [in]  df_den  Degrees of freedom in denominator.
+ * @param [out] result Output array.
+ */
+template <typename _DataType>
+INP_DLLEXPORT void dpnp_rng_f_c(void* result, _DataType df_num, _DataType df_den, size_t size);
+
 /**
  * @ingroup BACKEND_RANDOM_API
  * @brief math library implementation of random number generator (gamma distribution)

dpnp/backend/kernels/dpnp_krnl_random.cpp

@@ -176,6 +176,72 @@ void dpnp_rng_exponential_c(void* result, _DataType beta, size_t size)
     event_out.wait();
 }
 
+template <typename _DataType>
+void dpnp_rng_f_c(void* result, const _DataType df_num, const _DataType df_den, size_t size)
+{
+    if (!size)
+    {
+        return;
+    }
+    cl::sycl::vector_class<cl::sycl::event> no_deps;
+
+    const _DataType d_zero = (_DataType(0.0));
+
+    _DataType shape = 0.5*df_num;
+    _DataType scale = 2.0/df_num;
+    _DataType *den = nullptr;
+
+    _DataType* result1 = reinterpret_cast<_DataType*>(result);
+
+    if (dpnp_queue_is_cpu_c())
+    {
+        mkl_rng::gamma<_DataType> gamma_distribution1(shape, d_zero, scale);
+        auto event_out = mkl_rng::generate(gamma_distribution1, DPNP_RNG_ENGINE, size, result1);
+        event_out.wait();
+
+        den = reinterpret_cast<_DataType*>(dpnp_memory_alloc_c(size * sizeof(_DataType)));
+        if (den == nullptr)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_f_c() failed.");
+        }
+        shape = 0.5*df_den;
+        scale = 2.0/df_den;
+        mkl_rng::gamma<_DataType> gamma_distribution2(shape, d_zero, scale);
+        event_out = mkl_rng::generate(gamma_distribution2, DPNP_RNG_ENGINE, size, den);
+        event_out.wait();
+
+        event_out = mkl_vm::div(DPNP_QUEUE, size, result1, den, result1, no_deps,
+            mkl_vm::mode::ha);
+        event_out.wait();
+
+        dpnp_memory_free_c(den);
+    }
+    else
+    {
+        int errcode = vdRngGamma(VSL_RNG_METHOD_GAMMA_GNORM_ACCURATE, get_rng_stream(),
+            size, result1, shape, d_zero, scale);
+        if (errcode != VSL_STATUS_OK)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_f_c() failed.");
+        }
+        den = (_DataType *) mkl_malloc(size * sizeof(_DataType), 64);
+        if (den == nullptr)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_f_c() failed.");
+        }
+        shape = 0.5*df_den;
+        scale = 2.0/df_den;
+        errcode = vdRngGamma(VSL_RNG_METHOD_GAMMA_GNORM_ACCURATE, get_rng_stream(),
+            size, den, shape, d_zero, scale);
+        if (errcode != VSL_STATUS_OK)
+        {
+            throw std::runtime_error("DPNP RNG Error: dpnp_rng_f_c() failed.");
+        }
+        vmdDiv(size, result1, den, result1, VML_HA);
+        mkl_free(den);
+    }
+}
+
 template <typename _DataType>
 void dpnp_rng_gamma_c(void* result, _DataType shape, _DataType scale, size_t size)
 {
@@ -732,6 +798,8 @@ void func_map_init_random(func_map_t& fmap)
     fmap[DPNPFuncName::DPNP_FN_RNG_EXPONENTIAL][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_exponential_c<double>};
     fmap[DPNPFuncName::DPNP_FN_RNG_EXPONENTIAL][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_rng_exponential_c<float>};
 
+    fmap[DPNPFuncName::DPNP_FN_RNG_F][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_f_c<double>};
+
     fmap[DPNPFuncName::DPNP_FN_RNG_GAMMA][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_gamma_c<double>};
 
     fmap[DPNPFuncName::DPNP_FN_RNG_GAUSSIAN][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_gaussian_c<double>};

dpnp/backend/tests/test_random.cpp

@@ -32,6 +32,33 @@ TEST (TestBackendRandomBeta, test_seed) {
     }
 }
 
+TEST (TestBackendRandomF, test_seed) {
+    const size_t size = 256;
+    size_t seed = 10;
+    double dfnum = 10.4;
+    double dfden = 12.5;
+
+    auto QueueOptionsDevices = std::vector<QueueOptions>{ QueueOptions::CPU_SELECTOR,
+        QueueOptions::GPU_SELECTOR };
+
+    for (auto device_selector :  QueueOptionsDevices) {
+        dpnp_queue_initialize_c(device_selector);
+        double* result1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
+        double* result2 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
+
+        dpnp_rng_srand_c(seed);
+        dpnp_rng_f_c<double>(result1, dfnum, dfden, size);
+
+        dpnp_rng_srand_c(seed);
+        dpnp_rng_f_c<double>(result2, dfnum, dfden, size);
+
+        for (size_t i = 0; i < size; ++i)
+        {
+            EXPECT_NEAR (result1[i], result2[i], 0.004);
+        }
+    }
+}
+
 TEST (TestBackendRandomNormal, test_seed) {
     const size_t size = 256;
     size_t seed = 10;

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -99,6 +99,7 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName":  # need this na
         DPNP_FN_RNG_BINOMIAL
         DPNP_FN_RNG_CHISQUARE
         DPNP_FN_RNG_EXPONENTIAL
+        DPNP_FN_RNG_F
         DPNP_FN_RNG_GAMMA
         DPNP_FN_RNG_GAUSSIAN
         DPNP_FN_RNG_GEOMETRIC

dpnp/random/dpnp_algo_random.pyx

@@ -47,6 +47,7 @@ __all__ = [
     "dpnp_rng_binomial",
     "dpnp_rng_chisquare",
     "dpnp_rng_exponential",
+    "dpnp_rng_f",
     "dpnp_rng_gamma",
     "dpnp_rng_geometric",
     "dpnp_rng_gumbel",
@@ -79,6 +80,7 @@ ctypedef void(*fptr_dpnp_rng_beta_c_1out_t)(void *, double, double, size_t) exce
 ctypedef void(*fptr_dpnp_rng_binomial_c_1out_t)(void *, int, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_chi_square_c_1out_t)(void *, int, size_t) except +
 ctypedef void(*fptr_dpnp_rng_exponential_c_1out_t)(void *, double, size_t) except +
+ctypedef void(*fptr_dpnp_rng_f_c_1out_t)(void *, const double, const double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_gamma_c_1out_t)(void *, double, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_geometric_c_1out_t)(void *, float, size_t) except +
 ctypedef void(*fptr_dpnp_rng_gaussian_c_1out_t)(void *, double, double, size_t) except +
@@ -226,6 +228,31 @@ cpdef dparray dpnp_rng_exponential(double beta, size):
     return result
 
 
+cpdef dparray dpnp_rng_f(double df_num, double df_den, size):
+    """
+    Returns an array populated with samples from F distribution.
+    `dpnp_rng_f` generates a matrix filled with random floats sampled from a
+    univariate F distribution.
+    """
+
+    dtype = numpy.float64
+    # convert string type names (dparray.dtype) to C enum DPNPFuncType
+    cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(dtype)
+
+    # get the FPTR data structure
+    cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_RNG_F, 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=dtype)
+
+    cdef fptr_dpnp_rng_f_c_1out_t func = <fptr_dpnp_rng_f_c_1out_t > kernel_data.ptr
+    # call FPTR function
+    func(result.get_data(), df_num, df_den, result.size)
+
+    return result
+
+
 cpdef dparray dpnp_rng_gamma(double shape, double scale, size):
     """
     Returns an array populated with samples from gamma distribution.

dpnp/random/dpnp_iface_random.py

@@ -315,13 +315,32 @@ def f(dfnum, dfden, size=None):
 
     For full documentation refer to :obj:`numpy.random.f`.
 
-    Notes
-    -----
-    The function uses `numpy.random.f` on the backend and will be
-    executed on fallback backend.
+    Limitations
+    -----------
+    Parameters ``dfnum`` and ``dfden`` are supported as scalar.
+    Otherwise, :obj:`numpy.random.f(dfnum, dfden, size)` samples are drawn.
+    Output array data type is :obj:`dpnp.float64`.
+    Examples
+    --------
+    >>> dfnum, dfden = 3., 2.
+    >>> s = dpnp.random.f(dfnum, dfden, size)
 
     """
 
+    if not use_origin_backend(dfnum):
+        # TODO:
+        # array_like of floats for `dfnum` and `dfden`
+        if not dpnp.isscalar(dfnum):
+            pass
+        elif not dpnp.isscalar(dfden):
+            pass
+        elif dfnum <= 0:
+            pass
+        elif dfden <= 0:
+            pass
+        else:
+            return dpnp_rng_f(dfnum, dfden, size)
+
     return call_origin(numpy.random.f, dfnum, dfden, size)
 
 

tests/test_random.py

@@ -225,6 +225,34 @@ def test_seed(self):
         self.check_seed('exponential', {'scale': scale})
 
 
+class TestDistributionsF(TestDistribution):
+
+    def test_moments(self):
+        dfnum = 12.56
+        dfden = 13.0
+        # for dfden > 2
+        expected_mean = dfden / (dfden - 2)
+        # for dfden > 4
+        expected_var = 2 * (dfden ** 2) * (dfnum + dfden - 2) / \
+                       (dfnum * ((dfden - 2) ** 2) * ((dfden - 4)))
+        self.check_moments('f', expected_mean, expected_var,
+                           {'dfnum': dfnum, 'dfden': dfden})
+
+    def test_invalid_args(self):
+        size = 10
+        dfnum = -1.0   # positive `dfnum` is expected
+        dfden = 1.0    # OK
+        self.check_invalid_args('f', {'dfnum': dfnum, 'dfden': dfden})
+        dfnum = 1.0    # OK
+        dfden = -1.0   # positive `dfden` is expected
+        self.check_invalid_args('f', {'dfnum': dfnum, 'dfden': dfden})
+
+    def test_seed(self):
+        dfnum = 3.56   # `dfden` param for Wald distr
+        dfden = 2.8    # `dfden` param for Wald distr
+        self.check_seed('f', {'dfnum': dfnum, 'dfden': dfden})
+
+
 class TestDistributionsGamma(TestDistribution):
 
     def test_moments(self):

tests_external/skipped_tests_numpy.tbl

@@ -1466,6 +1466,7 @@ tests/test_random.py::TestRandomDist::test_dirichlet_bad_alpha
 tests/test_random.py::TestRandomDist::test_dirichlet_size
 tests/test_random.py::TestRandomDist::test_exponential
 tests/test_random.py::TestRandomDist::test_exponential_0
+tests/test_random.py::TestRandomDist::test_f
 tests/test_random.py::TestRandomDist::test_gamma
 tests/test_random.py::TestRandomDist::test_gamma_0
 tests/test_random.py::TestRandomDist::test_geometric
@@ -1559,6 +1560,7 @@ tests/test_randomstate.py::TestRandomDist::test_dirichlet
 tests/test_randomstate.py::TestRandomDist::test_dirichlet_size
 tests/test_randomstate.py::TestRandomDist::test_exponential
 tests/test_randomstate.py::TestRandomDist::test_exponential_0
+tests/test_randomstate.py::TestRandomDist::test_f
 tests/test_randomstate.py::TestRandomDist::test_gamma
 tests/test_randomstate.py::TestRandomDist::test_gamma_0
 tests/test_randomstate.py::TestRandomDist::test_geometric