Add functions std and var (#105)

Author: densmirn

dpnp/backend.pxd

@@ -83,12 +83,14 @@ cdef extern from "backend/backend_iface_fptr.hpp" namespace "DPNPFuncName":  # n
         DPNP_FN_SORT
         DPNP_FN_SQRT
         DPNP_FN_SQUARE
+        DPNP_FN_STD
         DPNP_FN_SUBTRACT
         DPNP_FN_SUM
         DPNP_FN_TAN
         DPNP_FN_TANH
         DPNP_FN_TRANSPOSE
         DPNP_FN_TRUNC
+        DPNP_FN_VAR
 
 cdef extern from "backend/backend_iface_fptr.hpp" namespace "DPNPFuncType":  # need this namespace for Enum import
     cdef enum DPNPFuncType "DPNPFuncType":

dpnp/backend/backend_iface.hpp

@@ -341,6 +341,52 @@ INP_DLLEXPORT void custom_argmax_c(void* array, void* result, size_t size);
 template <typename _DataType, typename _idx_DataType>
 INP_DLLEXPORT void custom_argmin_c(void* array, void* result, size_t size);
 
+/**
+ * @ingroup BACKEND_API
+ * @brief MKL implementation of std function
+ *
+ * @param [in]  array   Input array with data.
+ *
+ * @param [out] result  Output array with indeces.
+ * 
+ * @param [in]  shape   Shape of input array.
+ *
+ * @param [in]  ndim    Number of elements in shape.
+ *
+ * @param [in]  axis    Axis.
+ *
+ * @param [in]  naxis   Number of elements in axis.
+ *
+ * @param [in]  ddof    Delta degrees of freedom.
+ *
+ */
+template <typename _DataType, typename _ResultType>
+INP_DLLEXPORT void custom_std_c(void* array, void* result, size_t* shape,
+                                size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+
+/**
+ * @ingroup BACKEND_API
+ * @brief MKL implementation of var function
+ *
+ * @param [in]  array   Input array with data.
+ *
+ * @param [out] result  Output array with indeces.
+ * 
+ * @param [in]  shape   Shape of input array.
+ *
+ * @param [in]  ndim    Number of elements in shape.
+ *
+ * @param [in]  axis    Axis.
+ *
+ * @param [in]  naxis   Number of elements in axis.
+ *
+ * @param [in]  ddof    Delta degrees of freedom.
+ *
+ */
+template <typename _DataType, typename _ResultType>
+INP_DLLEXPORT void custom_var_c(void* array, void* result, size_t* shape,
+                                size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+
 #if 0 // Example for OpenCL kernel
 template <typename _DataType>
 void custom_dgemm_c_opencl(void* array_1, void* array_2, void* result_1, size_t size);

dpnp/backend/backend_iface_fptr.cpp

@@ -610,6 +610,11 @@ static func_map_t func_map_init()
     fmap[DPNPFuncName::DPNP_FN_SQUARE][eft_FLT][eft_FLT] = {eft_FLT, (void*)custom_elemwise_square_c<float>};
     fmap[DPNPFuncName::DPNP_FN_SQUARE][eft_DBL][eft_DBL] = {eft_DBL, (void*)custom_elemwise_square_c<double>};
 
+    fmap[DPNPFuncName::DPNP_FN_STD][eft_INT][eft_INT] = {eft_DBL, (void*)custom_std_c<int, double>};
+    fmap[DPNPFuncName::DPNP_FN_STD][eft_LNG][eft_LNG] = {eft_DBL, (void*)custom_std_c<long, double>};
+    fmap[DPNPFuncName::DPNP_FN_STD][eft_FLT][eft_FLT] = {eft_FLT, (void*)custom_std_c<float, float>};
+    fmap[DPNPFuncName::DPNP_FN_STD][eft_DBL][eft_DBL] = {eft_DBL, (void*)custom_std_c<double, double>};
+
     fmap[DPNPFuncName::DPNP_FN_SUBTRACT][eft_INT][eft_INT] = {eft_INT,
                                                               (void*)custom_elemwise_subtract_c<int, int, int>};
     fmap[DPNPFuncName::DPNP_FN_SUBTRACT][eft_INT][eft_LNG] = {eft_LNG,
@@ -668,5 +673,10 @@ static func_map_t func_map_init()
     fmap[DPNPFuncName::DPNP_FN_TRUNC][eft_FLT][eft_FLT] = {eft_FLT, (void*)custom_elemwise_trunc_c<float, float>};
     fmap[DPNPFuncName::DPNP_FN_TRUNC][eft_DBL][eft_DBL] = {eft_DBL, (void*)custom_elemwise_trunc_c<double, double>};
 
+    fmap[DPNPFuncName::DPNP_FN_VAR][eft_INT][eft_INT] = {eft_DBL, (void*)custom_var_c<int, double>};
+    fmap[DPNPFuncName::DPNP_FN_VAR][eft_LNG][eft_LNG] = {eft_DBL, (void*)custom_var_c<long, double>};
+    fmap[DPNPFuncName::DPNP_FN_VAR][eft_FLT][eft_FLT] = {eft_FLT, (void*)custom_var_c<float, float>};
+    fmap[DPNPFuncName::DPNP_FN_VAR][eft_DBL][eft_DBL] = {eft_DBL, (void*)custom_var_c<double, double>};
+
     return fmap;
 };

dpnp/backend/backend_iface_fptr.hpp

@@ -112,12 +112,14 @@ enum class DPNPFuncName : size_t
     DPNP_FN_SORT,      /**< Used in numpy.sort() implementation  */
     DPNP_FN_SQRT,      /**< Used in numpy.sqrt() implementation  */
     DPNP_FN_SQUARE,    /**< Used in numpy.square() implementation  */
+    DPNP_FN_STD,       /**< Used in numpy.std() implementation  */
     DPNP_FN_SUBTRACT,  /**< Used in numpy.subtract() implementation  */
     DPNP_FN_SUM,       /**< Used in numpy.sum() implementation  */
     DPNP_FN_TAN,       /**< Used in numpy.tan() implementation  */
     DPNP_FN_TANH,      /**< Used in numpy.tanh() implementation  */
     DPNP_FN_TRANSPOSE, /**< Used in numpy.transpose() implementation  */
     DPNP_FN_TRUNC,     /**< Used in numpy.trunc() implementation  */
+    DPNP_FN_VAR,       /**< Used in numpy.var() implementation  */
     DPNP_FN_LAST       /**< The latest element of the enumeration */
 };
 

dpnp/backend/custom_kernels_statistics.cpp

@@ -140,6 +140,7 @@ void custom_cov_c(void* array1_in, void* result1, size_t nrows, size_t ncols)
 
 template void custom_cov_c<double>(void* array1_in, void* result1, size_t nrows, size_t ncols);
 
+
 template <typename _DataType>
 class custom_max_c_kernel;
 
@@ -282,3 +283,91 @@ template void
 template void
     custom_min_c<long>(void* array1_in, void* result1, size_t* shape, size_t ndim, size_t* axis, size_t naxis);
 template void custom_min_c<int>(void* array1_in, void* result1, size_t* shape, size_t ndim, size_t* axis, size_t naxis);
+
+
+template <typename _DataType, typename _ResultType>
+void custom_std_c(void* array1_in, void* result1, size_t* shape, size_t ndim, size_t* axis, size_t naxis, size_t ddof)
+{
+    _DataType* array1 = reinterpret_cast<_DataType*>(array1_in);
+    _ResultType* result = reinterpret_cast<_ResultType*>(result1);
+
+    _ResultType* var = reinterpret_cast<_ResultType*>(dpnp_memory_alloc_c(1 * sizeof(_ResultType)));
+    custom_var_c<_DataType, _ResultType>(array1, var, shape, ndim, axis, naxis, ddof);
+
+    custom_elemwise_sqrt_c<_ResultType, _ResultType>(var, result, 1);
+
+    dpnp_memory_free_c(var);
+
+#if 0
+    std::cout << "std result " << res[0] << "\n";
+#endif
+}
+
+template void custom_std_c<int, double>(void* array1_in, void* result1, size_t* shape,
+                                        size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+template void custom_std_c<long, double>(void* array1_in, void* result1, size_t* shape,
+                                         size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+template void custom_std_c<float, float>(void* array1_in, void* result1, size_t* shape,
+                                         size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+template void custom_std_c<double, double>(void* array1_in, void* result1, size_t* shape,
+                                           size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+
+
+template <typename _DataType, typename _ResultType>
+class custom_var_c_kernel;
+
+template <typename _DataType, typename _ResultType>
+void custom_var_c(void* array1_in, void* result1, size_t* shape, size_t ndim, size_t* axis, size_t naxis, size_t ddof)
+{
+    cl::sycl::event event;
+    _DataType* array1 = reinterpret_cast<_DataType*>(array1_in);
+    _ResultType* result = reinterpret_cast<_ResultType*>(result1);
+
+    _ResultType* mean = reinterpret_cast<_ResultType*>(dpnp_memory_alloc_c(1 * sizeof(_ResultType)));
+    custom_mean_c<_DataType, _ResultType>(array1, mean, shape, ndim, axis, naxis);
+    _ResultType mean_val = mean[0];
+
+    size_t size = 1;
+    for (size_t i = 0; i < ndim; ++i) {
+        size *= shape[i];
+    }
+
+    _ResultType* squared_deviations = reinterpret_cast<_ResultType*>(dpnp_memory_alloc_c(size * sizeof(_ResultType)));
+
+    cl::sycl::range<1> gws(size);
+    event = DPNP_QUEUE.submit([&](cl::sycl::handler& cgh) {
+        cgh.parallel_for<class custom_var_c_kernel<_DataType, _ResultType> >(
+            gws,
+            [=](cl::sycl::id<1> global_id)
+        {
+            size_t i = global_id[0]; /*for (size_t i = 0; i < size; ++i)*/
+            {
+                _ResultType deviation = (_ResultType)array1[i] - mean_val;
+                squared_deviations[i] = deviation * deviation;
+            }
+        }); /* parallel_for */
+    });     /* queue.submit */
+
+    event.wait();
+
+    custom_mean_c<_ResultType, _ResultType>(squared_deviations, mean, shape, ndim, axis, naxis);
+    mean_val = mean[0];
+
+    result[0] = mean_val * size / (size - ddof);
+
+    dpnp_memory_free_c(mean);
+    dpnp_memory_free_c(squared_deviations);
+
+#if 0
+    std::cout << "var result " << res[0] << "\n";
+#endif
+}
+
+template void custom_var_c<int, double>(void* array1_in, void* result1, size_t* shape,
+                                        size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+template void custom_var_c<long, double>(void* array1_in, void* result1, size_t* shape,
+                                         size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+template void custom_var_c<float, float>(void* array1_in, void* result1, size_t* shape,
+                                         size_t ndim, size_t* axis, size_t naxis, size_t ddof);
+template void custom_var_c<double, double>(void* array1_in, void* result1, size_t* shape,
+                                           size_t ndim, size_t* axis, size_t naxis, size_t ddof);

dpnp/backend_statistics.pyx

@@ -46,18 +46,45 @@ __all__ += [
     "dpnp_max",
     "dpnp_mean",
     "dpnp_median",
-    "dpnp_min"
+    "dpnp_min",
+    "dpnp_std",
+    "dpnp_var",
 ]
 
 
 # C function pointer to the C library template functions
 ctypedef void(*fptr_custom_cov_1in_1out_t)(void * , void * , size_t, size_t)
-
+ctypedef void(*fptr_custom_std_var_1in_1out_t)(void * , void * , size_t * , size_t, size_t * , size_t, size_t)
 
 # C function pointer to the C library template functions
 ctypedef void(*custom_statistic_1in_1out_func_ptr_t)(void * , void * , size_t * , size_t, size_t * , size_t)
 
 
+cdef dparray call_fptr_custom_std_var_1in_1out(DPNPFuncName fptr_name, dparray a, ddof):
+
+    """ Convert string type names (dparray.dtype) to C enum DPNPFuncType """
+    cdef DPNPFuncType param_type = dpnp_dtype_to_DPNPFuncType(a.dtype)
+
+    """ get the FPTR data structure """
+    cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(fptr_name, param_type, DPNP_FT_NONE)
+
+    result_type = dpnp_DPNPFuncType_to_dtype( < size_t > kernel_data.return_type)
+    """ Create result array with type given by FPTR data """
+    cdef dparray result = dparray((1,), dtype=result_type)
+
+    cdef fptr_custom_std_var_1in_1out_t func = <fptr_custom_std_var_1in_1out_t > kernel_data.ptr
+
+    # stub for interface support
+    cdef dparray_shape_type axis
+    cdef Py_ssize_t axis_size = 0
+
+    """ Call FPTR function """
+    func(a.get_data(), result.get_data(), < size_t * > a._dparray_shape.data(),
+         a.ndim, < size_t * > axis.data(), axis_size, ddof)
+
+    return result
+
+
 cpdef dpnp_average(dparray x1):
     array_sum = dpnp_sum(x1)
 
@@ -432,3 +459,11 @@ cpdef dparray dpnp_min(dparray input, axis):
     dpnp_array = dpnp.array(result_array, dtype=input.dtype)
     dpnp_result_array = dpnp_array.reshape(output_shape)
     return dpnp_result_array
+
+
+cpdef dparray dpnp_std(dparray a, size_t ddof):
+    return call_fptr_custom_std_var_1in_1out(DPNP_FN_STD, a, ddof)
+
+
+cpdef dparray dpnp_var(dparray a, size_t ddof):
+    return call_fptr_custom_std_var_1in_1out(DPNP_FN_VAR, a, ddof)

dpnp/dparray.pyx

@@ -518,6 +518,15 @@ cdef class dparray:
         """
         return repeat(self, *args, **kwds)
 
+    def std(self, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
+        """ Returns the variance of the array elements, along given axis.
+
+        .. seealso::
+           :func:`dpnp.var` for full documentation,
+
+        """
+        return std(self, axis, dtype, out, ddof, keepdims)
+
     def transpose(self, *axes):
         """ Returns a view of the array with axes permuted.
 
@@ -528,6 +537,15 @@ cdef class dparray:
         """
         return transpose(self, axes)
 
+    def var(self, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
+        """ Returns the standard deviation of the array elements along given axis.
+
+        .. seealso::
+           :func:`dpnp.std` for full documentation,
+
+        """
+        return var(self, axis, dtype, out, ddof, keepdims)
+
     def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
         """ Return our function pointer regusted by `ufunc` parameter
         """