ENH: impl random.pareto (#518)

Author: samir-nasibli

dpnp/backend/include/dpnp_iface.hpp

@@ -855,6 +855,20 @@ INP_DLLEXPORT void dpnp_rng_negative_binomial_c(void* result, double a, double p
 template <typename _DataType>
 INP_DLLEXPORT void dpnp_rng_normal_c(void* result, _DataType mean, _DataType stddev, size_t size);
 
+/**
+ * @ingroup BACKEND_API
+ * @brief math library implementation of random number generator (Pareto distribution)
+ *
+ * @param [in]  size   Number of elements in `result` arrays.
+ *
+ * @param [in]  alpha  Shape of the distribution, alpha.
+ *
+ * @param [out] result Output array.
+ *
+ */
+template <typename _DataType>
+INP_DLLEXPORT void dpnp_rng_pareto_c(void* result, double alpha, size_t size);
+
 /**
  * @ingroup BACKEND_API
  * @brief math library implementation of random number generator (poisson distribution)

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -137,6 +137,7 @@ enum class DPNPFuncName : size_t
     DPNP_FN_RNG_MULTIVARIATE_NORMAL,  /**< Used in numpy.random.multivariate_normal() implementation  */
     DPNP_FN_RNG_NEGATIVE_BINOMIAL,    /**< Used in numpy.random.negative_binomial() implementation  */
     DPNP_FN_RNG_NORMAL,               /**< Used in numpy.random.normal() implementation  */
+    DPNP_FN_RNG_PARETO,               /**< Used in numpy.random.pareto() implementation  */
     DPNP_FN_RNG_POISSON,              /**< Used in numpy.random.poisson() implementation  */
     DPNP_FN_RNG_RAYLEIGH,             /**< Used in numpy.random.rayleigh() implementation  */
     DPNP_FN_RNG_STANDARD_CAUCHY,      /**< Used in numpy.random.standard_cauchy() implementation  */

dpnp/backend/kernels/dpnp_krnl_random.cpp

@@ -439,6 +439,29 @@ void dpnp_rng_normal_c(void* result, _DataType mean, _DataType stddev, size_t si
     event_out.wait();
 }
 
+template <typename _DataType>
+void dpnp_rng_pareto_c(void* result, double alpha, size_t size)
+{
+    if (!size)
+    {
+        return;
+    }
+    cl::sycl::vector_class<cl::sycl::event> no_deps;
+
+    const _DataType d_zero = _DataType(0.0);
+    const _DataType d_one = _DataType(1.0);
+    _DataType neg_rec_alp = -1.0/alpha;
+
+    _DataType* result1 = reinterpret_cast<_DataType*>(result);
+
+    mkl_rng::uniform<_DataType> distribution(d_zero, d_one);
+    auto event_out = mkl_rng::generate(distribution, DPNP_RNG_ENGINE, size, result1);
+    event_out.wait();
+
+    event_out = mkl_vm::powx(DPNP_QUEUE, size, result1, neg_rec_alp, result1, no_deps, mkl_vm::mode::ha);
+    event_out.wait();
+}
+
 template <typename _DataType>
 void dpnp_rng_poisson_c(void* result, double lambda, size_t size)
 {
@@ -618,6 +641,8 @@ void func_map_init_random(func_map_t& fmap)
 
     fmap[DPNPFuncName::DPNP_FN_RNG_NORMAL][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_normal_c<double>};
 
+    fmap[DPNPFuncName::DPNP_FN_RNG_PARETO][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_pareto_c<double>};
+
     fmap[DPNPFuncName::DPNP_FN_RNG_POISSON][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_rng_poisson_c<int>};
 
     fmap[DPNPFuncName::DPNP_FN_RNG_RAYLEIGH][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_rayleigh_c<double>};

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -110,6 +110,7 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName":  # need this na
         DPNP_FN_RNG_MULTIVARIATE_NORMAL
         DPNP_FN_RNG_NEGATIVE_BINOMIAL
         DPNP_FN_RNG_NORMAL
+        DPNP_FN_RNG_PARETO
         DPNP_FN_RNG_POISSON
         DPNP_FN_RNG_RAYLEIGH
         DPNP_FN_RNG_STANDARD_CAUCHY

dpnp/random/dpnp_algo_random.pyx

@@ -57,6 +57,7 @@ __all__ = [
     "dpnp_multivariate_normal",
     "dpnp_negative_binomial",
     "dpnp_normal",
+    "dpnp_pareto",
     "dpnp_poisson",
     "dpnp_randn",
     "dpnp_random",
@@ -92,6 +93,7 @@ ctypedef void(*fptr_dpnp_rng_multivariate_normal_c_1out_t)(void *,
                                                            size_t) except +
 ctypedef void(*fptr_dpnp_rng_negative_binomial_c_1out_t)(void *, double, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_normal_c_1out_t)(void *, double, double, size_t) except +
+ctypedef void(*fptr_dpnp_rng_pareto_c_1out_t)(void *, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_poisson_c_1out_t)(void *, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_rayleigh_c_1out_t)(void *, double, size_t) except +
 ctypedef void(*fptr_dpnp_rng_standard_cauchy_c_1out_t)(void *, size_t) except +
@@ -595,6 +597,32 @@ cpdef dparray dpnp_normal(double loc, double scale, size):
     return result
 
 
+cpdef dparray dpnp_pareto(double alpha, size):
+    """
+    Returns an array populated with samples from Pareto distribution.
+    `dpnp_pareto` generates a matrix filled with random floats sampled from a
+    univariate Pareto distribution of `alpha`.
+
+    """
+
+    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_PARETO, 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_pareto_c_1out_t func = <fptr_dpnp_rng_pareto_c_1out_t > kernel_data.ptr
+    # call FPTR function
+    func(result.get_data(), alpha, result.size)
+
+    return result
+
+
 cpdef dparray dpnp_poisson(double lam, size):
     """
     Returns an array populated with samples from Poisson distribution.

dpnp/random/dpnp_iface_random.py

@@ -809,13 +809,30 @@ def pareto(a, size=None):
 
     For full documentation refer to :obj:`numpy.random.pareto`.
 
-    Notes
-    -----
-    The function uses `numpy.random.pareto` on the backend and
-    will be executed on fallback backend.
+    Limitations
+    -----------
+    Parameter ``a`` is supported as a scalar.
+    Otherwise, :obj:`numpy.random.pareto(a, size)` samples are drawn.
+    Output array data type is :obj:`dpnp.float64`.
+
+    Examples
+    --------
+    Draw samples from the distribution:
+    >>> a = .5  # alpha
+    >>> s = dpnp.random.pareto(a, 1000)
 
     """
 
+    if not use_origin_backend(a):
+        # TODO:
+        # array_like of floats for `a`
+        if not dpnp.isscalar(a):
+            pass
+        elif a <= 0:
+            pass
+        else:
+            return dpnp_pareto(a, size)
+
     return call_origin(numpy.random.pareto, a, size)
 
 

tests/test_random.py

@@ -557,6 +557,25 @@ def test_seed(self):
         self.check_seed('normal', {'loc': loc, 'scale': scale})
 
 
+class TestDistributionsPareto(TestDistribution):
+
+    def test_moments(self):
+        a = 30.0
+        expected_mean = a / (a - 1)
+        expected_var = a / (((a - 1)**2) * (a - 2))
+        self.check_moments('pareto', expected_mean,
+                           expected_var, {'a': a})
+
+    def test_invalid_args(self):
+        size = 10
+        a = -1.0  # positive `a` is expected
+        self.check_invalid_args('pareto', {'a': a})
+
+    def test_seed(self):
+        a = 3.0  # a param for pareto distr
+        self.check_seed('pareto', {'a': a})
+
+
 class TestDistributionsPoisson(TestDistribution):
 
     def test_extreme_value(self):