implement dpnp.kaiser

Author: Vahid Tavanashad

dpnp/backend/extensions/window/CMakeLists.txt

@@ -26,6 +26,7 @@
 
 set(python_module_name _window_impl)
 set(_module_src
+    ${CMAKE_CURRENT_SOURCE_DIR}/kaiser.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/window_py.cpp
 )
 

dpnp/backend/extensions/window/kaiser.cpp

@@ -0,0 +1,180 @@
+//*****************************************************************************
+// Copyright (c) 2025, Intel Corporation
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+// - Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// - Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+// THE POSSIBILITY OF SUCH DAMAGE.
+//*****************************************************************************
+
+#include "kaiser.hpp"
+#include "utils/output_validation.hpp"
+#include "utils/type_dispatch.hpp"
+#include "utils/type_utils.hpp"
+#include <sycl/sycl.hpp>
+
+/**
+ * Version of SYCL DPC++ 2025.1 compiler where an issue with
+ * sycl::ext::intel::math::cyl_bessel_i0(x) is fully resolved.
+ */
+#ifndef __SYCL_COMPILER_BESSEL_I0_SUPPORT
+#define __SYCL_COMPILER_BESSEL_I0_SUPPORT 20241208L
+#endif
+
+#if __SYCL_COMPILER_VERSION >= __SYCL_COMPILER_BESSEL_I0_SUPPORT
+#include <sycl/ext/intel/math.hpp>
+#endif
+
+#include "../kernels/elementwise_functions/i0.hpp"
+
+namespace dpnp::extensions::window
+{
+namespace dpctl_td_ns = dpctl::tensor::type_dispatch;
+
+typedef sycl::event (*kaiser_fn_ptr_t)(sycl::queue &,
+                                       char *,
+                                       const std::size_t,
+                                       const float,
+                                       const std::vector<sycl::event> &);
+
+static kaiser_fn_ptr_t kaiser_dispatch_vector[dpctl_td_ns::num_types];
+
+template <typename T>
+class KaiserFunctor
+{
+private:
+    T *data = nullptr;
+    const std::size_t N;
+    const float beta;
+
+public:
+    KaiserFunctor(T *data, const std::size_t N, const float beta)
+        : data(data), N(N), beta(beta)
+    {
+    }
+
+    void operator()(sycl::id<1> id) const
+    {
+#if __SYCL_COMPILER_VERSION >= __SYCL_COMPILER_BESSEL_I0_SUPPORT
+        using sycl::ext::intel::math::cyl_bessel_i0;
+#else
+        using dpnp::kernels::i0::impl::cyl_bessel_i0;
+#endif
+
+        const auto i = id.get(0);
+        const T alpha = (N - 1) / T(2);
+        const T tmp = (i - alpha) / alpha;
+        data[i] = cyl_bessel_i0(beta * sycl::sqrt(1 - tmp * tmp)) /
+                  cyl_bessel_i0(beta);
+    }
+};
+
+template <typename T, template <typename> class Functor>
+sycl::event kaiser_impl(sycl::queue &q,
+                        char *result,
+                        const std::size_t nelems,
+                        const float beta,
+                        const std::vector<sycl::event> &depends)
+{
+    dpctl::tensor::type_utils::validate_type_for_device<T>(q);
+
+    T *res = reinterpret_cast<T *>(result);
+
+    sycl::event kaiser_ev = q.submit([&](sycl::handler &cgh) {
+        cgh.depends_on(depends);
+
+        using KaiserKernel = Functor<T>;
+        cgh.parallel_for<KaiserKernel>(sycl::range<1>(nelems),
+                                       KaiserKernel(res, nelems, beta));
+    });
+
+    return kaiser_ev;
+}
+
+template <typename fnT, typename T>
+struct KaiserFactory
+{
+    fnT get()
+    {
+        if constexpr (std::is_floating_point_v<T>) {
+            return kaiser_impl<T, KaiserFunctor>;
+        }
+        else {
+            return nullptr;
+        }
+    }
+};
+
+std::pair<sycl::event, sycl::event>
+    py_kaiser(sycl::queue &exec_q,
+              const float beta,
+              const dpctl::tensor::usm_ndarray &result,
+              const std::vector<sycl::event> &depends)
+{
+    dpctl::tensor::validation::CheckWritable::throw_if_not_writable(result);
+
+    int nd = result.get_ndim();
+    if (nd != 1) {
+        throw py::value_error("Array should be 1d");
+    }
+
+    if (!dpctl::utils::queues_are_compatible(exec_q, {result.get_queue()})) {
+        throw py::value_error(
+            "Execution queue is not compatible with allocation queue.");
+    }
+
+    const bool is_result_c_contig = result.is_c_contiguous();
+    if (!is_result_c_contig) {
+        throw py::value_error("The result input array is not c-contiguous.");
+    }
+
+    size_t nelems = result.get_size();
+    if (nelems == 0) {
+        return std::make_pair(sycl::event{}, sycl::event{});
+    }
+
+    int result_typenum = result.get_typenum();
+    auto array_types = dpctl_td_ns::usm_ndarray_types();
+    int result_type_id = array_types.typenum_to_lookup_id(result_typenum);
+    auto fn = kaiser_dispatch_vector[result_type_id];
+
+    if (fn == nullptr) {
+        throw std::runtime_error("Type of given array is not supported");
+    }
+
+    char *result_typeless_ptr = result.get_data();
+    sycl::event kaiser_ev =
+        fn(exec_q, result_typeless_ptr, nelems, beta, depends);
+    sycl::event args_ev =
+        dpctl::utils::keep_args_alive(exec_q, {result}, {kaiser_ev});
+
+    return std::make_pair(args_ev, kaiser_ev);
+}
+
+void init_kaiser_dispatch_vectors()
+{
+    dpctl_td_ns::DispatchVectorBuilder<kaiser_fn_ptr_t, KaiserFactory,
+                                       dpctl_td_ns::num_types>
+        contig;
+    contig.populate_dispatch_vector(kaiser_dispatch_vector);
+
+    return;
+}
+
+} // namespace dpnp::extensions::window

dpnp/backend/extensions/window/kaiser.hpp

@@ -0,0 +1,41 @@
+//*****************************************************************************
+// Copyright (c) 2025, Intel Corporation
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+// - Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// - Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+// THE POSSIBILITY OF SUCH DAMAGE.
+//*****************************************************************************
+
+#pragma once
+
+#include <dpctl4pybind11.hpp>
+#include <sycl/sycl.hpp>
+
+namespace dpnp::extensions::window
+{
+extern std::pair<sycl::event, sycl::event>
+    py_kaiser(sycl::queue &exec_q,
+              const float beta,
+              const dpctl::tensor::usm_ndarray &result,
+              const std::vector<sycl::event> &depends);
+
+extern void init_kaiser_dispatch_vectors(void);
+
+} // namespace dpnp::extensions::window

dpnp/backend/extensions/window/window_py.cpp

@@ -35,6 +35,7 @@
 #include "common.hpp"
 #include "hamming.hpp"
 #include "hanning.hpp"
+#include "kaiser.hpp"
 
 namespace window_ns = dpnp::extensions::window;
 namespace py = pybind11;
@@ -111,4 +112,12 @@ PYBIND11_MODULE(_window_impl, m)
               py::arg("sycl_queue"), py::arg("result"),
               py::arg("depends") = py::list());
     }
+
+    {
+        window_ns::init_kaiser_dispatch_vectors();
+
+        m.def("_kaiser", window_ns::py_kaiser, "Call Kaiser kernel",
+              py::arg("sycl_queue"), py::arg("beta"), py::arg("result"),
+              py::arg("depends") = py::list());
+    }
 }