Leverage dpnp.cumlogsumexp through dpctl.tensor implementation

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -105,8 +105,6 @@ enum class DPNPFuncName : size_t
     DPNP_FN_COV,           /**< Used in numpy.cov() impl  */
     DPNP_FN_CROSS,         /**< Used in numpy.cross() impl  */
     DPNP_FN_CUMPROD,       /**< Used in numpy.cumprod() impl  */
-    DPNP_FN_CUMPROD_EXT,   /**< Used in numpy.cumprod() impl, requires extra
-                              parameters */
     DPNP_FN_CUMSUM,        /**< Used in numpy.cumsum() impl  */
     DPNP_FN_DEGREES,       /**< Used in numpy.degrees() impl  */
     DPNP_FN_DEGREES_EXT,   /**< Used in numpy.degrees() impl, requires extra

dpnp/backend/kernels/dpnp_krnl_mathematical.cpp

@@ -364,14 +364,6 @@ template <typename _DataType_input, typename _DataType_output>
 void (*dpnp_cumprod_default_c)(void *, void *, size_t) =
     dpnp_cumprod_c<_DataType_input, _DataType_output>;
 
-template <typename _DataType_input, typename _DataType_output>
-DPCTLSyclEventRef (*dpnp_cumprod_ext_c)(DPCTLSyclQueueRef,
-                                        void *,
-                                        void *,
-                                        size_t,
-                                        const DPCTLEventVectorRef) =
-    dpnp_cumprod_c<_DataType_input, _DataType_output>;
-
 template <typename _KernelNameSpecialization1,
           typename _KernelNameSpecialization2>
 class dpnp_cumsum_c_kernel;
@@ -1153,15 +1145,6 @@ void func_map_init_mathematical(func_map_t &fmap)
     fmap[DPNPFuncName::DPNP_FN_CUMPROD][eft_DBL][eft_DBL] = {
         eft_DBL, (void *)dpnp_cumprod_default_c<double, double>};
 
-    fmap[DPNPFuncName::DPNP_FN_CUMPROD_EXT][eft_INT][eft_INT] = {
-        eft_LNG, (void *)dpnp_cumprod_ext_c<int32_t, int64_t>};
-    fmap[DPNPFuncName::DPNP_FN_CUMPROD_EXT][eft_LNG][eft_LNG] = {
-        eft_LNG, (void *)dpnp_cumprod_ext_c<int64_t, int64_t>};
-    fmap[DPNPFuncName::DPNP_FN_CUMPROD_EXT][eft_FLT][eft_FLT] = {
-        eft_FLT, (void *)dpnp_cumprod_ext_c<float, float>};
-    fmap[DPNPFuncName::DPNP_FN_CUMPROD_EXT][eft_DBL][eft_DBL] = {
-        eft_DBL, (void *)dpnp_cumprod_ext_c<double, double>};
-
     fmap[DPNPFuncName::DPNP_FN_CUMSUM][eft_INT][eft_INT] = {
         eft_LNG, (void *)dpnp_cumsum_default_c<int32_t, int64_t>};
     fmap[DPNPFuncName::DPNP_FN_CUMSUM][eft_LNG][eft_LNG] = {

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -37,7 +37,6 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName":  # need this na
         DPNP_FN_CHOOSE_EXT
         DPNP_FN_COPY_EXT
         DPNP_FN_CORRELATE_EXT
-        DPNP_FN_CUMPROD_EXT
         DPNP_FN_DEGREES_EXT
         DPNP_FN_DIAG_INDICES_EXT
         DPNP_FN_DIAGONAL_EXT
@@ -127,9 +126,6 @@ cdef extern from "dpnp_iface.hpp":
 
 
 # C function pointer to the C library template functions
-ctypedef c_dpctl.DPCTLSyclEventRef(*fptr_1in_1out_t)(c_dpctl.DPCTLSyclQueueRef,
-                                                     void *, void * , size_t,
-                                                     const c_dpctl.DPCTLEventVectorRef)
 ctypedef c_dpctl.DPCTLSyclEventRef(*fptr_1in_1out_strides_t)(c_dpctl.DPCTLSyclQueueRef,
                                                              void *, const size_t, const size_t,
                                                              const shape_elem_type * , const shape_elem_type * ,

dpnp/dpnp_algo/dpnp_algo.pyx

@@ -147,56 +147,6 @@ cdef dpnp_DPNPFuncType_to_dtype(size_t type):
         utils.checker_throw_type_error("dpnp_DPNPFuncType_to_dtype", type)
 
 
-cdef utils.dpnp_descriptor call_fptr_1in_1out(DPNPFuncName fptr_name,
-                                              utils.dpnp_descriptor x1,
-                                              shape_type_c result_shape,
-                                              utils.dpnp_descriptor out=None,
-                                              func_name=None):
-
-    """ Convert type (x1.dtype) to C enum DPNPFuncType """
-    cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(x1.dtype)
-
-    """ get the FPTR data structure """
-    cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(fptr_name, param1_type, param1_type)
-
-    result_type = dpnp_DPNPFuncType_to_dtype( < size_t > kernel_data.return_type)
-
-    cdef utils.dpnp_descriptor result
-
-    if out is None:
-        """ Create result array with type given by FPTR data """
-        x1_obj = x1.get_array()
-        result = utils.create_output_descriptor(result_shape,
-                                                kernel_data.return_type,
-                                                None,
-                                                device=x1_obj.sycl_device,
-                                                usm_type=x1_obj.usm_type,
-                                                sycl_queue=x1_obj.sycl_queue)
-    else:
-        if out.dtype != result_type:
-            utils.checker_throw_value_error(func_name, 'out.dtype', out.dtype, result_type)
-        if out.shape != result_shape:
-            utils.checker_throw_value_error(func_name, 'out.shape', out.shape, result_shape)
-
-        result = out
-
-        utils.get_common_usm_allocation(x1, result)  # check USM allocation is common
-
-    result_sycl_queue = result.get_array().sycl_queue
-
-    cdef c_dpctl.SyclQueue q = <c_dpctl.SyclQueue> result_sycl_queue
-    cdef c_dpctl.DPCTLSyclQueueRef q_ref = q.get_queue_ref()
-
-    cdef fptr_1in_1out_t func = <fptr_1in_1out_t > kernel_data.ptr
-
-    cdef c_dpctl.DPCTLSyclEventRef event_ref = func(q_ref, x1.get_data(), result.get_data(), x1.size, NULL)
-
-    with nogil: c_dpctl.DPCTLEvent_WaitAndThrow(event_ref)
-    c_dpctl.DPCTLEvent_Delete(event_ref)
-
-    return result
-
-
 cdef utils.dpnp_descriptor call_fptr_1in_1out_strides(DPNPFuncName fptr_name,
                                                       utils.dpnp_descriptor x1,
                                                       object dtype=None,

dpnp/dpnp_algo/dpnp_algo_mathematical.pxi

@@ -43,7 +43,6 @@ __all__ += [
     "dpnp_fmax",
     "dpnp_fmin",
     "dpnp_modf",
-    "dpnp_nancumprod",
     "dpnp_trapz",
 ]
 
@@ -56,18 +55,6 @@ ctypedef c_dpctl.DPCTLSyclEventRef(*ftpr_custom_trapz_2in_1out_with_2size_t)(c_d
                                                                              const c_dpctl.DPCTLEventVectorRef)
 
 
-cpdef utils.dpnp_descriptor dpnp_cumprod(utils.dpnp_descriptor x1):
-    # instead of x1.shape, (x1.size, ) is passed to the function
-    # due to the following:
-    # >>> import numpy
-    # >>> a = numpy.array([[1, 2], [2, 3]])
-    # >>> res = numpy.cumprod(a)
-    # >>> res.shape
-    # (4,)
-
-    return call_fptr_1in_1out(DPNP_FN_CUMPROD_EXT, x1, (x1.size,))
-
-
 cpdef utils.dpnp_descriptor dpnp_ediff1d(utils.dpnp_descriptor x1):
 
     if x1.size <= 1:
@@ -226,19 +213,6 @@ cpdef tuple dpnp_modf(utils.dpnp_descriptor x1):
     return (result1.get_pyobj(), result2.get_pyobj())
 
 
-cpdef utils.dpnp_descriptor dpnp_nancumprod(utils.dpnp_descriptor x1):
-    cur_x1 = x1.get_pyobj().copy()
-
-    cur_x1_flatiter = cur_x1.flat
-
-    for i in range(cur_x1.size):
-        if dpnp.isnan(cur_x1_flatiter[i]):
-            cur_x1_flatiter[i] = 1
-
-    x1_desc = dpnp.get_dpnp_descriptor(cur_x1, copy_when_nondefault_queue=False)
-    return dpnp_cumprod(x1_desc)
-
-
 cpdef utils.dpnp_descriptor dpnp_trapz(utils.dpnp_descriptor y1, utils.dpnp_descriptor x1, double dx):
 
     cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(y1.dtype)

dpnp/dpnp_iface_nanfunctions.py

@@ -39,18 +39,8 @@
 
 import warnings
 
-import numpy
-
 import dpnp
 
-# pylint: disable=no-name-in-module
-from .dpnp_algo import (
-    dpnp_nancumprod,
-)
-from .dpnp_utils import (
-    call_origin,
-)
-
 __all__ = [
     "nanargmax",
     "nanargmin",
@@ -249,19 +239,40 @@ def nanargmin(a, axis=None, out=None, *, keepdims=False):
     return dpnp.argmin(a, axis=axis, out=out, keepdims=keepdims)
 
 
-def nancumprod(x1, **kwargs):
+def nancumprod(a, axis=None, dtype=None, out=None):
     """
     Return the cumulative product of array elements over a given axis treating
-    Not a Numbers (NaNs) as one.
+    Not a Numbers (NaNs) as zero. The cumulative product does not change when
+    NaNs are encountered and leading NaNs are replaced by ones.
 
     For full documentation refer to :obj:`numpy.nancumprod`.
 
-    Limitations
-    -----------
-    Parameter `x` is supported as :class:`dpnp.ndarray`.
-    Keyword argument `kwargs` is currently unsupported.
-    Otherwise the function will be executed sequentially on CPU.
-    Input array data types are limited by supported DPNP :ref:`Data types`.
+    Parameters
+    ----------
+    a : {dpnp.ndarray, usm_ndarray}
+        Input array.
+    axis : {None, int}, optional
+        Axis along which the cumulative product is computed. The default
+        (``None``) is to compute the cumulative product over the flattened
+        array.
+    dtype : {None, dtype}, optional
+        Type of the returned array and of the accumulator in which the elements
+        are summed. If `dtype` is not specified, it defaults to the dtype of
+        `a`, unless `a` has an integer dtype with a precision less than that of
+        the default platform integer. In that case, the default platform
+        integer is used.
+    out : {None, dpnp.ndarray, usm_ndarray}, optional
+        Alternative output array in which to place the result. It must have the
+        same shape and buffer length as the expected output but the type will
+        be cast if necessary.
+
+    Returns
+    -------
+    out : dpnp.ndarray
+        A new array holding the result is returned unless `out` is specified as
+        :class:`dpnp.ndarray`, in which case a reference to `out` is returned.
+        The result has the same size as `a`, and the same shape as `a` if `axis`
+        is not ``None`` or `a` is a 1-d array.
 
     See Also
     --------
@@ -271,22 +282,26 @@ def nancumprod(x1, **kwargs):
     Examples
     --------
     >>> import dpnp as np
-    >>> a = np.array([1., np.nan])
-    >>> result = np.nancumprod(a)
-    >>> [x for x in result]
-    [1.0, 1.0]
-    >>> b = np.array([[1., 2., np.nan], [4., np.nan, 6.]])
-    >>> result = np.nancumprod(b)
-    >>> [x for x in result]
-    [1.0, 2.0, 2.0, 8.0, 8.0, 48.0]
+    >>> np.nancumprod(np.array(1))
+    array([1])
+    >>> np.nancumprod(np.array([1]))
+    array([1])
+    >>> np.nancumprod(np.array([1, np.nan]))
+    array([1., 1.])
+    >>> a = np.array([[1, 2], [3, np.nan]])
+    >>> np.nancumprod(a)
+    array([1., 2., 6., 6.])
+    >>> np.nancumprod(a, axis=0)
+    array([[1., 2.],
+           [3., 2.]])
+    >>> np.nancumprod(a, axis=1)
+    array([[1., 2.],
+           [3., 3.]])
 
     """
 
-    x1_desc = dpnp.get_dpnp_descriptor(x1, copy_when_nondefault_queue=False)
-    if x1_desc and not kwargs:
-        return dpnp_nancumprod(x1_desc).get_pyobj()
-
-    return call_origin(numpy.nancumprod, x1, **kwargs)
+    a, _ = _replace_nan(a, 1)
+    return dpnp.cumprod(a, axis=axis, dtype=dtype, out=out)
 
 
 def nancumsum(a, axis=None, dtype=None, out=None):

dpnp/dpnp_iface_trigonometric.py

@@ -71,6 +71,7 @@
     "cbrt",
     "cos",
     "cosh",
+    "cumlogsumexp",
     "deg2rad",
     "degrees",
     "exp",
@@ -665,6 +666,94 @@ def _get_accumulation_res_dt(a, dtype, _out):
 )
 
 
+def cumlogsumexp(
+    x, /, *, axis=None, dtype=None, include_initial=False, out=None
+):
+    """
+    Calculates the cumulative logarithm of the sum of elements in the input
+    array `x`.
+
+    Parameters
+    ----------
+    x : {dpnp.ndarray, usm_ndarray}
+        Input array, expected to have a real-valued data type.
+    axis : {None, int or tuple of ints}, optional
+        Axis or axes along which values must be computed. If a tuple of unique
+        integers, values are computed over multiple axes. If ``None``, the
+        result is computed over the entire array.
+        Default: ``None``.
+    dtype : {None, dtype}, optional
+        Data type of the returned array. If ``None``, the default data type is
+        inferred from the "kind" of the input array data type.
+
+        - If `x` has a real-valued floating-point data type, the returned array
+          will have the same data type as `x`.
+        - If `x` has a boolean or integral data type, the returned array will
+          have the default floating point data type for the device where input
+          array `x` is allocated.
+        - If `x` has a complex-valued floating-point data type, an error is
+          raised.
+
+        If the data type (either specified or resolved) differs from the data
+        type of `x`, the input array elements are cast to the specified data
+        type before computing the result.
+        Default: ``None``.
+    include_initial : {None, bool}, optional
+        A boolean indicating whether to include the initial value (i.e., the
+        additive identity, zero) as the first value along the provided axis in
+        the output.
+        Default: ``False``.
+    out : {None, dpnp.ndarray, usm_ndarray}, optional
+        The array into which the result is written. The data type of `out` must
+        match the expected shape and the expected data type of the result or
+        (if provided) `dtype`. If ``None`` then a new array is returned.
+        Default: ``None``.
+
+    Returns
+    -------
+    out : dpnp.ndarray
+        An array containing the results. If the result was computed over the
+        entire array, a zero-dimensional array is returned. The returned array
+        has the data type as described in the `dtype` parameter description
+        above.
+
+    Note
+    ----
+    This function is equivalent of `numpy.logaddexp.reduce`.
+
+    See Also
+    --------
+    :obj:`dpnp.logsumexp` : Logarithm of the sum of elements of the inputs,
+                            element-wise.
+
+    Examples
+    --------
+    >>> import dpnp as np
+    >>> a = np.ones(10)
+    >>> np.cumlogsumexp(a)
+    array([1.        , 1.69314718, 2.09861229, 2.38629436, 2.60943791,
+           2.79175947, 2.94591015, 3.07944154, 3.19722458, 3.30258509])
+
+    """
+
+    dpnp.check_supported_arrays_type(x)
+    if x.ndim > 1 and axis is None:
+        usm_x = dpnp.ravel(x).get_array()
+    else:
+        usm_x = dpnp.get_usm_ndarray(x)
+
+    return dpnp_wrap_reduction_call(
+        x,
+        out,
+        dpt.cumulative_logsumexp,
+        _get_accumulation_res_dt,
+        usm_x,
+        axis=axis,
+        dtype=dtype,
+        include_initial=include_initial,
+    )
+
+
 def deg2rad(x1):
     """
     Convert angles from degrees to radians.