ENH: impl random.power (#520)

* ENH: impl random.power

* disabled incorrect test checks for dpnp.random.power

Author: samir-nasibli

dpnp/backend/include/dpnp_iface.hpp

@@ -883,6 +883,20 @@ INP_DLLEXPORT void dpnp_rng_pareto_c(void* result, double alpha, size_t size);
 template <typename _DataType>
 INP_DLLEXPORT void dpnp_rng_poisson_c(void* result, double lambda, size_t size);
 
+/**
+ * @ingroup BACKEND_API
+ * @brief math library implementation of random number generator (power 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_power_c(void* result, double alpha, size_t size);
+
 /**
  * @ingroup BACKEND_API
  * @brief math library implementation of random number generator (rayleigh distribution)

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -139,6 +139,7 @@ enum class DPNPFuncName : size_t
     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_POWER,                /**< Used in numpy.random.power() implementation  */
     DPNP_FN_RNG_RAYLEIGH,             /**< Used in numpy.random.rayleigh() implementation  */
     DPNP_FN_RNG_STANDARD_CAUCHY,      /**< Used in numpy.random.standard_cauchy() implementation  */
     DPNP_FN_RNG_STANDARD_EXPONENTIAL, /**< Used in numpy.random.standard_exponential() implementation  */

dpnp/backend/kernels/dpnp_krnl_random.cpp

@@ -477,6 +477,30 @@ void dpnp_rng_poisson_c(void* result, double lambda, size_t size)
     event_out.wait();
 }
 
+template <typename _DataType>
+void dpnp_rng_power_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_rayleigh_c(void* result, _DataType scale, size_t size)
 {
@@ -713,6 +737,8 @@ void func_map_init_random(func_map_t& fmap)
 
     fmap[DPNPFuncName::DPNP_FN_RNG_POISSON][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_rng_poisson_c<int>};
 
+    fmap[DPNPFuncName::DPNP_FN_RNG_POWER][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_power_c<double>};
+
     fmap[DPNPFuncName::DPNP_FN_RNG_RAYLEIGH][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_rayleigh_c<double>};
 
     fmap[DPNPFuncName::DPNP_FN_RNG_STANDARD_CAUCHY][eft_DBL][eft_DBL] = {eft_DBL,

dpnp/dpnp_algo/dpnp_algo.pxd

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

dpnp/random/dpnp_algo_random.pyx

@@ -59,6 +59,7 @@ __all__ = [
     "dpnp_normal",
     "dpnp_pareto",
     "dpnp_poisson",
+    "dpnp_rng_power",
     "dpnp_randn",
     "dpnp_random",
     "dpnp_rayleigh",
@@ -96,6 +97,7 @@ ctypedef void(*fptr_dpnp_rng_negative_binomial_c_1out_t)(void *, double, double,
 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_power_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 +
 ctypedef void(*fptr_dpnp_rng_standard_exponential_c_1out_t)(void *, size_t) except +
@@ -661,6 +663,31 @@ cpdef dparray dpnp_poisson(double lam, size):
     return result
 
 
+cpdef dparray dpnp_rng_power(double alpha, size):
+    """
+    Returns an array populated with samples from power distribution.
+    `dpnp_power` generates a matrix filled with random floats sampled from a
+    univariate power 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_POWER, 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_power_c_1out_t func = <fptr_dpnp_rng_power_c_1out_t > kernel_data.ptr
+    # call FPTR function
+    func(result.get_data(), alpha, result.size)
+
+    return result
+
+
 cpdef dparray dpnp_randn(dims):
     """
     Returns an array populated with samples from standard normal distribution.

dpnp/random/dpnp_iface_random.py

@@ -889,18 +889,35 @@ def poisson(lam=1.0, size=None):
 def power(a, size=None):
     """Power distribution.
 
-    Draws samples in [0, 1] from a power distribution with positive exponent
-    a - 1.
+    Draws samples in [0, 1] from a power distribution with positive
+    exponent a - 1.
 
     For full documentation refer to :obj:`numpy.random.power`.
 
-    Notes
-    -----
-    The function uses `numpy.random.power` on the backend and
-    will be executed on fallback backend.
+    Limitations
+    -----------
+    Parameter ``a`` is supported as a scalar.
+    Otherwise, :obj:`numpy.random.power(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.power(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_rng_power(a, size)
+
     return call_origin(numpy.random.power, a, size)
 
 

tests/test_random.py

@@ -598,6 +598,26 @@ def test_seed(self):
         self.check_seed('poisson', {'lam': lam})
 
 
+class TestDistributionsPower(TestDistribution):
+
+    def test_moments(self):
+        a = 30.0
+        neg_a = -a
+        expected_mean = neg_a / (neg_a - 1)
+        expected_var = neg_a / (((neg_a - 1)**2) * (neg_a - 2))
+        self.check_moments('power', expected_mean,
+                           expected_var, {'a': a})
+
+    def test_invalid_args(self):
+        size = 10
+        a = -1.0  # positive `a` is expected
+        self.check_invalid_args('power', {'a': a})
+
+    def test_seed(self):
+        a = 3.0  # a param for pareto distr
+        self.check_seed('power', {'a': a})
+
+
 class TestDistributionsRayleigh(TestDistribution):
 
     def test_extreme_value(self):

tests_external/skipped_tests_numpy.tbl

@@ -1493,6 +1493,7 @@ tests/test_random.py::TestRandomDist::test_random_integers_deprecated
 tests/test_random.py::TestRandomDist::test_random_integers_max_int
 tests/test_random.py::TestRandomDist::test_rayleigh
 tests/test_random.py::TestRandomDist::test_rayleigh_0
+tests/test_random.py::TestRandomDist::test_power
 tests/test_random.py::TestRandomDist::test_scalar_exception_propagation
 tests/test_random.py::TestRandomDist::test_shuffle
 tests/test_random.py::TestRandomDist::test_shuffle_masked
@@ -1588,6 +1589,7 @@ tests/test_randomstate.py::TestRandomDist::test_random_sample
 tests/test_randomstate.py::TestRandomDist::test_rand_singleton
 tests/test_randomstate.py::TestRandomDist::test_rayleigh
 tests/test_randomstate.py::TestRandomDist::test_rayleigh_0
+tests/test_randomstate.py::TestRandomDist::test_power
 tests/test_randomstate.py::TestRandomDist::test_scalar_exception_propagation
 tests/test_randomstate.py::TestRandomDist::test_shuffle
 tests/test_randomstate.py::TestRandomDist::test_shuffle_masked