ENH: impl random.triangular (#573)

Author: samir-nasibli

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -152,6 +152,7 @@ enum class DPNPFuncName : size_t
     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_TRIANGULAR,           /**< Used in numpy.random.triangular() 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/include/dpnp_iface_random.hpp

@@ -370,6 +370,20 @@ INP_DLLEXPORT void dpnp_rng_standard_normal_c(void* result, const size_t size);
 template <typename _DataType>
 INP_DLLEXPORT void dpnp_rng_standard_t_c(void* result, const _DataType df, const size_t size);
 
+/**
+ * @ingroup BACKEND_RANDOM_API
+ * @brief math library implementation of random number generator (Triangular distribution)
+ *
+ * @param [in]  size   Number of elements in `result` arrays.
+ * @param [in]  x_min  Lower limit.
+ * @param [in]  x_mode The value where the peak of the distribution occurs.
+ * @param [in]  x_max  Upper limit.
+ * @param [out] result Output array.
+ */
+template <typename _DataType>
+INP_DLLEXPORT void dpnp_rng_triangular_c(
+    void* result, const _DataType x_min, const _DataType x_mode, const _DataType x_max, const size_t size);
+
 /**
  * @ingroup BACKEND_RANDOM_API
  * @brief math library implementation of random number generator (uniform distribution)

dpnp/backend/kernels/dpnp_krnl_random.cpp

@@ -763,6 +763,76 @@ void dpnp_rng_standard_t_c(void* result, const _DataType df, const size_t size)
     }
 }
 
+template <typename _KernelNameSpecialization>
+class dpnp_rng_triangular_ration_acceptance_c_kernel;
+
+template <typename _DataType>
+void dpnp_rng_triangular_c(
+    void* result, const _DataType x_min, const _DataType x_mode, const _DataType x_max, const size_t size)
+{
+    if (!size)
+    {
+        return;
+    }
+    _DataType* result1 = reinterpret_cast<_DataType*>(result);
+    const _DataType d_zero = (_DataType(0));
+    const _DataType d_one = (_DataType(1));
+
+    _DataType ratio, lpr, rpr;
+
+    mkl_rng::uniform<_DataType> uniform_distribution(d_zero, d_one);
+    auto event_out = mkl_rng::generate(uniform_distribution, DPNP_RNG_ENGINE, size, result1);
+    event_out.wait();
+
+    {
+        _DataType wtot, wl, wr;
+
+        wtot = x_max - x_min;
+        wl = x_mode - x_min;
+        wr = x_max - x_mode;
+
+        ratio = wl / wtot;
+        lpr = wl * wtot;
+        rpr = wr * wtot;
+    }
+
+    if (!(0 <= ratio && ratio <= 1))
+    {
+        throw std::runtime_error("DPNP RNG Error: dpnp_rng_triangular_c() failed.");
+    }
+
+    cl::sycl::range<1> gws(size);
+    auto kernel_parallel_for_func = [=](cl::sycl::id<1> global_id) {
+        size_t i = global_id[0];
+        if (ratio <= 0)
+        {
+            result1[i] = x_max - cl::sycl::sqrt(result1[i] * rpr);
+        }
+        else if (ratio >= 1)
+        {
+            result1[i] = x_min + cl::sycl::sqrt(result1[i] * lpr);
+        }
+        else
+        {
+            _DataType ui = result1[i];
+            if (ui > ratio)
+            {
+                result1[i] = x_max - cl::sycl::sqrt((1.0 - ui) * rpr);
+            }
+            else
+            {
+                result1[i] = x_min + cl::sycl::sqrt(ui * lpr);
+            }
+        }
+    };
+    auto kernel_func = [&](cl::sycl::handler& cgh) {
+        cgh.parallel_for<class dpnp_rng_triangular_ration_acceptance_c_kernel<_DataType>>(gws,
+                                                                                          kernel_parallel_for_func);
+    };
+    event_out = DPNP_QUEUE.submit(kernel_func);
+    event_out.wait();
+}
+
 template <typename _DataType>
 void dpnp_rng_uniform_c(void* result, const long low, const long high, const size_t size)
 {
@@ -867,6 +937,8 @@ void func_map_init_random(func_map_t& fmap)
                                                                          (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_TRIANGULAR][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_rng_triangular_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>};
     fmap[DPNPFuncName::DPNP_FN_RNG_UNIFORM][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_rng_uniform_c<int>};

dpnp/dpnp_algo/dpnp_algo.pxd

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

dpnp/random/dpnp_algo_random.pyx

@@ -71,6 +71,7 @@ __all__ = [
     "dpnp_rng_standard_gamma",
     "dpnp_rng_standard_normal",
     "dpnp_rng_standard_t",
+    "dpnp_rng_triangular",
     "dpnp_rng_uniform",
     "dpnp_rng_weibull"
 ]
@@ -113,6 +114,11 @@ ctypedef void(*fptr_dpnp_rng_standard_exponential_c_1out_t)(void * , const size_
 ctypedef void(*fptr_dpnp_rng_standard_gamma_c_1out_t)(void * , const double, const size_t) except +
 ctypedef void(*fptr_dpnp_rng_standard_normal_c_1out_t)(void * , const size_t) except +
 ctypedef void(*fptr_dpnp_rng_standard_t_c_1out_t)(void * , const double, const size_t) except +
+ctypedef void(*fptr_dpnp_rng_triangular_c_1out_t)(void * ,
+                                                  const double,
+                                                  const double,
+                                                  const double,
+                                                  const size_t) except +
 ctypedef void(*fptr_dpnp_rng_uniform_c_1out_t)(void * , const long, const long, const size_t) except +
 ctypedef void(*fptr_dpnp_rng_weibull_c_1out_t)(void * , const double, const size_t) except +
 
@@ -988,6 +994,32 @@ cpdef dparray dpnp_rng_standard_t(double df, size):
     return result
 
 
+cpdef dparray dpnp_rng_triangular(double left, double mode, double right, size):
+    """
+    Returns an array populated with samples from triangular distribution.
+    `dpnp_rng_triangular` generates a matrix filled with random floats sampled from a
+    univariate triangular 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_TRIANGULAR, 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_triangular_c_1out_t func = <fptr_dpnp_rng_triangular_c_1out_t > kernel_data.ptr
+    # call FPTR function
+    func(result.get_data(), left, mode, right, result.size)
+
+    return result
+
+
 cpdef dparray dpnp_rng_uniform(long low, long high, size, dtype):
     """
     Returns an array populated with samples from standard uniform distribution.

dpnp/random/dpnp_iface_random.py

@@ -1446,7 +1446,6 @@ def standard_t(df, size=None):
             pass
         else:
             return dpnp_rng_standard_t(df, size)
-    print("here")
     return call_origin(numpy.random.standard_t, df, size)
 
 
@@ -1458,13 +1457,39 @@ def triangular(left, mode, right, size=None):
 
     For full documentation refer to :obj:`numpy.random.triangular`.
 
-    Notes
-    -----
-    The function uses `numpy.random.triangular` on the backend and
-    will be executed on fallback backend.
+    Limitations
+    -----------
+    Parameter ``left``, ``mode`` and ``right`` are supported as scalar.
+    Otherwise, :obj:`numpy.random.triangular(left, mode, right, size)`
+    samples are drawn.
+    Output array data type is :obj:`dpnp.float64`.
+
+    Examples
+    --------
+    Draw samples from the distribution:
+    >>> df = 2.
+    >>> s = dpnp.random.triangular(-3, 0, 8, 1000000)
 
     """
 
+    if not use_origin_backend(left):
+        # TODO:
+        # array_like of floats for `left`, `mode`, `right`.
+        if not dpnp.isscalar(left):
+            pass
+        elif not dpnp.isscalar(mode):
+            pass
+        elif not dpnp.isscalar(right):
+            pass
+        elif left > mode:
+            pass
+        elif mode > right:
+            pass
+        elif left == right:
+            pass
+        else:
+            return dpnp_rng_triangular(left, mode, right, size)
+
     return call_origin(numpy.random.triangular, left, mode, right, size)
 
 

tests/test_random.py

@@ -757,6 +757,35 @@ def test_seed(self):
         self.check_seed('standard_t', {'df': 10.0})
 
 
+class TestDistributionsTriangular(TestDistribution):
+
+    def test_moments(self):
+        left = 1.0
+        mode = 2.0
+        right = 3.0
+        expected_mean = (left + mode + right) / 3
+        expected_var = (left ** 2 + mode ** 2 + right ** 2 - left * mode -left * right -mode * right) / 18
+        self.check_moments('triangular', expected_mean,
+                           expected_var, {'left': left, 'mode': mode, 'right': right})
+
+    def test_invalid_args(self):
+        left = 2.0   # `left` is expected <= `mode`
+        mode = 1.0   # `mode` is expected > `left`
+        right = 3.0  # OK
+        self.check_invalid_args('triangular', {'left': left, 'mode': mode, 'right': right})
+
+        left = 1.0   # OK
+        mode = 3.0   # `mode` is expected <= `right`
+        right = 2.0  # `right` is expected > `mode`
+        self.check_invalid_args('triangular', {'left': left, 'mode': mode, 'right': right})
+
+    def test_seed(self):
+        left = 1.0
+        mode = 2.0
+        right = 3.0
+        self.check_seed('triangular', {'left': left, 'mode': mode, 'right': right})
+
+
 class TestDistributionsUniform(TestDistribution):
 
     def test_extreme_value(self):

tests_external/skipped_tests_numpy.tbl

@@ -1506,6 +1506,7 @@ 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_triangular
 tests/test_random.py::TestRandomDist::test_uniform
 tests/test_random.py::TestRandomDist::test_uniform_range_bounds
 tests/test_random.py::TestRandomDist::test_weibull
@@ -1605,6 +1606,7 @@ 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_triangular
 tests/test_randomstate.py::TestRandomDist::test_uniform
 tests/test_randomstate.py::TestRandomDist::test_uniform_range_bounds
 tests/test_randomstate.py::TestRandomDist::test_weibull