From 1ffaf50cc7155ac117a3ef5c2ac0cde37858a33d Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Fri, 25 Oct 2024 20:11:33 +0200 Subject: [PATCH 01/11] Add _sinc to ufunc extension --- dpnp/backend/extensions/ufunc/CMakeLists.txt | 1 + .../ufunc/elementwise_functions/common.cpp | 2 + .../ufunc/elementwise_functions/sinc.cpp | 125 ++++++++++++++++++ .../ufunc/elementwise_functions/sinc.hpp | 35 +++++ .../kernels/elementwise_functions/sinc.hpp | 56 ++++++++ 5 files changed, 219 insertions(+) create mode 100644 dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp create mode 100644 dpnp/backend/extensions/ufunc/elementwise_functions/sinc.hpp create mode 100644 dpnp/backend/kernels/elementwise_functions/sinc.hpp diff --git a/dpnp/backend/extensions/ufunc/CMakeLists.txt b/dpnp/backend/extensions/ufunc/CMakeLists.txt index 1d7bc2f4ac1b..f02030935e1e 100644 --- a/dpnp/backend/extensions/ufunc/CMakeLists.txt +++ b/dpnp/backend/extensions/ufunc/CMakeLists.txt @@ -38,6 +38,7 @@ set(_elementwise_sources ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/ldexp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/logaddexp2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/radians.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/sinc.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/spacing.cpp ) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp index e29def5df4ac..4ed4cf1f5632 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp @@ -38,6 +38,7 @@ #include "ldexp.hpp" #include "logaddexp2.hpp" #include "radians.hpp" +#include "sinc.hpp" #include "spacing.hpp" namespace py = pybind11; @@ -62,6 +63,7 @@ void init_elementwise_functions(py::module_ m) init_ldexp(m); init_logaddexp2(m); init_radians(m); + init_sinc(m); init_spacing(m); } } // namespace dpnp::extensions::ufunc diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp new file mode 100644 index 000000000000..cc8432f075ba --- /dev/null +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp @@ -0,0 +1,125 @@ +//***************************************************************************** +// Copyright (c) 2024, 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 + +#include "dpctl4pybind11.hpp" + +#include "kernels/elementwise_functions/sinc.hpp" +#include "populate.hpp" +#include "sinc.hpp" + +// include a local copy of elementwise common header from dpctl tensor: +// dpctl/tensor/libtensor/source/elementwise_functions/elementwise_functions.hpp +// TODO: replace by including dpctl header once available +#include "../../elementwise_functions/elementwise_functions.hpp" + +// dpctl tensor headers +#include "kernels/elementwise_functions/common.hpp" +#include "utils/type_dispatch.hpp" + +namespace dpnp::extensions::ufunc +{ +namespace py = pybind11; +namespace py_int = dpnp::extensions::py_internal; + +namespace impl +{ +namespace ew_cmn_ns = dpctl::tensor::kernels::elementwise_common; +namespace td_ns = dpctl::tensor::type_dispatch; + +/** + * @brief A factory to define pairs of supported types for which + * sycl::sinc function is available. + * + * @tparam T Type of input vector `a` and of result vector `y`. + */ +template +struct OutputType +{ + using value_type = + typename std::disjunction, + td_ns::TypeMapResultEntry, + td_ns::TypeMapResultEntry, + td_ns::DefaultResultEntry>::result_type; +}; + +using dpnp::kernels::sinc::SincFunctor; + +template +using ContigFunctor = ew_cmn_ns::UnaryContigFunctor, + vec_sz, + n_vecs, + enable_sg_loadstore>; + +template +using StridedFunctor = ew_cmn_ns:: + UnaryStridedFunctor>; + +using ew_cmn_ns::unary_contig_impl_fn_ptr_t; +using ew_cmn_ns::unary_strided_impl_fn_ptr_t; + +static unary_contig_impl_fn_ptr_t sinc_contig_dispatch_vector[td_ns::num_types]; +static int sinc_output_typeid_vector[td_ns::num_types]; +static unary_strided_impl_fn_ptr_t + sinc_strided_dispatch_vector[td_ns::num_types]; + +MACRO_POPULATE_DISPATCH_VECTORS(sinc); +} // namespace impl + +void init_sinc(py::module_ m) +{ + using arrayT = dpctl::tensor::usm_ndarray; + using event_vecT = std::vector; + { + impl::populate_sinc_dispatch_vectors(); + using impl::sinc_contig_dispatch_vector; + using impl::sinc_output_typeid_vector; + using impl::sinc_strided_dispatch_vector; + + auto sinc_pyapi = [&](const arrayT &src, const arrayT &dst, + sycl::queue &exec_q, + const event_vecT &depends = {}) { + return py_int::py_unary_ufunc( + src, dst, exec_q, depends, sinc_output_typeid_vector, + sinc_contig_dispatch_vector, sinc_strided_dispatch_vector); + }; + m.def("_sinc", sinc_pyapi, "", py::arg("src"), py::arg("dst"), + py::arg("sycl_queue"), py::arg("depends") = py::list()); + + auto sinc_result_type_pyapi = [&](const py::dtype &dtype) { + return py_int::py_unary_ufunc_result_type( + dtype, sinc_output_typeid_vector); + }; + m.def("_sinc_result_type", sinc_result_type_pyapi); + } +} +} // namespace dpnp::extensions::ufunc diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.hpp b/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.hpp new file mode 100644 index 000000000000..b6cd7df6195c --- /dev/null +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.hpp @@ -0,0 +1,35 @@ +//***************************************************************************** +// Copyright (c) 2024, 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 + +namespace py = pybind11; + +namespace dpnp::extensions::ufunc +{ +void init_sinc(py::module_ m); +} // namespace dpnp::extensions::ufunc diff --git a/dpnp/backend/kernels/elementwise_functions/sinc.hpp b/dpnp/backend/kernels/elementwise_functions/sinc.hpp new file mode 100644 index 000000000000..b135567cf96f --- /dev/null +++ b/dpnp/backend/kernels/elementwise_functions/sinc.hpp @@ -0,0 +1,56 @@ +//***************************************************************************** +// Copyright (c) 2024, 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 + +namespace dpnp::kernels::sinc +{ +template +struct SincFunctor +{ + // is function constant for given argT + using is_constant = typename std::false_type; + // constant value, if constant + // constexpr resT constant_value = resT{}; + // is function defined for sycl::vec + using supports_vec = typename std::false_type; + // do both argT and resT support subgroup store/load operation + using supports_sg_loadstore = typename std::true_type; + + resT operator()(const argT &x) const + { + constexpr argT pi = + static_cast(3.1415926535897932384626433832795029L); + const argT y = pi * x; + + if (y == argT(0)) { + return argT(1); + } + return sycl::sinpi(x) / y; + } +}; +} // namespace dpnp::kernels::sinc From 7a33fc9bbc77b53e2ab61fff0239a888ccda4da6 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Sun, 27 Oct 2024 19:35:29 +0100 Subject: [PATCH 02/11] Add dpnp.sinc function to the interface --- dpnp/dpnp_iface_mathematical.py | 70 +++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index ecbe35c741cf..9984b48958a4 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -130,6 +130,7 @@ "round", "sign", "signbit", + "sinc", "spacing", "subtract", "sum", @@ -3750,6 +3751,75 @@ def real_if_close(a, tol=100): ) +_SINC_DOCSTRING = r""" +Return the normalized sinc function. + +The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument +:math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not +only everywhere continuous but also infinitely differentiable. + +For full documentation refer to :obj:`numpy.sinc`. + +Parameters +---------- +x : {dpnp.ndarray, usm_ndarray} + Input array, expected to have floating-point type. +out : {None, dpnp.ndarray, usm_ndarray}, optional + Output array to populate. + Array must have the correct shape and the expected data type. + Default: ``None``. +order : {"C", "F", "A", "K"}, optional + Memory layout of the newly output array, if parameter `out` is ``None``. + Default: ``"K"``. + +Returns +------- +out : dpnp.ndarray + An array containing the element-wise result of the ``sinc(x)`` function. + The data type of the returned array is determined by the Type Promotion + Rules. + +Limitations +----------- +Parameters `where` and `subok` are supported with their default values. +Keyword argument `kwargs` is currently unsupported. +Otherwise ``NotImplementedError`` exception will be raised. + +Notes +----- +The name sinc is short for "sine cardinal" or "sinus cardinalis". + +The sinc function is used in various signal processing applications, including +in anti-aliasing, in the construction of a Lanczos resampling filter, and in +interpolation. + +For bandlimited interpolation of discrete-time signals, the ideal interpolation +kernel is proportional to the sinc function. + +Examples +-------- +>>> import dpnp as np +>>> x = np.linspace(-4, 4, 41, dtype=dpnp.float64) +>>> np.sinc(x) # result may vary + array([ 0. , -0.04923628, -0.08409186, -0.08903844, -0.05846808, + 0. , 0.06682066, 0.11643488, 0.12613779, 0.08504448, + 0. , -0.10394325, -0.18920668, -0.21623621, -0.15591488, + 0. , 0.23387232, 0.50455115, 0.75682673, 0.93548928, + 1. , 0.93548928, 0.75682673, 0.50455115, 0.23387232, + 0. , -0.15591488, -0.21623621, -0.18920668, -0.10394325, + 0. , 0.08504448, 0.12613779, 0.11643488, 0.06682066, + 0. , -0.05846808, -0.08903844, -0.08409186, -0.04923628, + 0. ]) +""" + +sinc = DPNPUnaryFunc( + "sinc", + ufi._sinc_result_type, + ufi._sinc, + _SINC_DOCSTRING, +) + + _SPACING_DOCSTRING = """ Return the distance between `x` and the nearest adjacent number. From 85807880b158919863f3b3d73b2487c7f49e2428 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Mon, 28 Oct 2024 12:11:03 +0100 Subject: [PATCH 03/11] Add third party tests --- dpnp/dpnp_iface_mathematical.py | 2 +- .../cupy/math_tests/test_special.py | 33 +++++++++++++++++++ 2 files changed, 34 insertions(+), 1 deletion(-) create mode 100644 tests/third_party/cupy/math_tests/test_special.py diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 9984b48958a4..c0b4b2aa8470 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -3763,7 +3763,7 @@ def real_if_close(a, tol=100): Parameters ---------- x : {dpnp.ndarray, usm_ndarray} - Input array, expected to have floating-point type. + Input array, expected to have floating-point data type. out : {None, dpnp.ndarray, usm_ndarray}, optional Output array to populate. Array must have the correct shape and the expected data type. diff --git a/tests/third_party/cupy/math_tests/test_special.py b/tests/third_party/cupy/math_tests/test_special.py new file mode 100644 index 000000000000..d257e2f98e4d --- /dev/null +++ b/tests/third_party/cupy/math_tests/test_special.py @@ -0,0 +1,33 @@ +import unittest + +import pytest + +import dpnp as cupy +from tests.third_party.cupy import testing + + +class TestSpecial(unittest.TestCase): + + @testing.for_dtypes(["e", "f", "d"]) + @testing.numpy_cupy_allclose(rtol=1e-3) + def test_i0(self, xp, dtype): + a = testing.shaped_random((2, 3), xp, dtype) + return xp.i0(a) + + @testing.for_dtypes(["e", "f", "d", "F", "D"]) + @testing.numpy_cupy_allclose(atol=1e-3) + def test_sinc(self, xp, dtype): + + if dtype in [cupy.float16, cupy.float32, cupy.complex64]: + pytest.xfail( + reason="XXX: np2.0: numpy 1.26 uses a wrong " + "promotion; numpy 2.0 is OK" + ) + + a = testing.shaped_random((2, 3), xp, dtype, scale=1) + return xp.sinc(a) + + @testing.for_dtypes(["e", "f", "d", "F", "D"]) + def test_sinc_zero(self, dtype): + a = cupy.sinc(cupy.zeros(1, dtype=dtype)) + testing.assert_array_equal(a, cupy.ones(1, dtype=dtype)) From f5307f6abbf5d57157ea45b3ee034fc36e674637 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Mon, 28 Oct 2024 13:31:39 +0100 Subject: [PATCH 04/11] Add complex type handling --- .../ufunc/elementwise_functions/sinc.cpp | 12 +- .../kernels/elementwise_functions/sinc.hpp | 158 +++++++++++++++++- 2 files changed, 158 insertions(+), 12 deletions(-) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp index cc8432f075ba..961e8f2fcd8e 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp @@ -59,11 +59,13 @@ namespace td_ns = dpctl::tensor::type_dispatch; template struct OutputType { - using value_type = - typename std::disjunction, - td_ns::TypeMapResultEntry, - td_ns::TypeMapResultEntry, - td_ns::DefaultResultEntry>::result_type; + using value_type = typename std::disjunction< + td_ns::TypeMapResultEntry, + td_ns::TypeMapResultEntry, + td_ns::TypeMapResultEntry, + td_ns::TypeMapResultEntry>, + td_ns::TypeMapResultEntry>, + td_ns::DefaultResultEntry>::result_type; }; using dpnp::kernels::sinc::SincFunctor; diff --git a/dpnp/backend/kernels/elementwise_functions/sinc.hpp b/dpnp/backend/kernels/elementwise_functions/sinc.hpp index b135567cf96f..02b534ee35d5 100644 --- a/dpnp/backend/kernels/elementwise_functions/sinc.hpp +++ b/dpnp/backend/kernels/elementwise_functions/sinc.hpp @@ -25,32 +25,176 @@ #pragma once +#define SYCL_EXT_ONEAPI_COMPLEX +#if __has_include() +#include +#else +#include +#endif + #include +// dpctl tensor headers +#include "utils/type_utils.hpp" + namespace dpnp::kernels::sinc { -template +namespace tu_ns = dpctl::tensor::type_utils; + +namespace impl +{ +namespace exprm_ns = sycl::ext::oneapi::experimental; + +template +inline Tp sin(const Tp &in) +{ + if constexpr (tu_ns::is_complex::value) { + using realTp = typename Tp::value_type; + + constexpr realTp q_nan = std::numeric_limits::quiet_NaN(); + + realTp const &in_re = std::real(in); + realTp const &in_im = std::imag(in); + + const bool in_re_finite = sycl::isfinite(in_re); + const bool in_im_finite = sycl::isfinite(in_im); + /* + * Handle the nearly-non-exceptional cases where + * real and imaginary parts of input are finite. + */ + if (in_re_finite && in_im_finite) { + Tp res = exprm_ns::sin(exprm_ns::complex(in)); // sin(in); + if (in_re == realTp(0)) { + res.real(sycl::copysign(realTp(0), in_re)); + } + return res; + } + + /* + * since sin(in) = -I * sinh(I * in), for special cases, + * we calculate real and imaginary parts of z = sinh(I * in) and + * then return { imag(z) , -real(z) } which is sin(in). + */ + const realTp x = -in_im; + const realTp y = in_re; + const bool xfinite = in_im_finite; + const bool yfinite = in_re_finite; + /* + * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN. + * The sign of 0 in the result is unspecified. Choice = normally + * the same as dNaN. + * + * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN). + * The sign of 0 in the result is unspecified. Choice = normally + * the same as d(NaN). + */ + if (x == realTp(0) && !yfinite) { + const realTp sinh_im = q_nan; + const realTp sinh_re = sycl::copysign(realTp(0), x * sinh_im); + return Tp{sinh_im, -sinh_re}; + } + + /* + * sinh(+-Inf +- I 0) = +-Inf + I +-0. + * + * sinh(NaN +- I 0) = d(NaN) + I +-0. + */ + if (y == realTp(0) && !xfinite) { + if (std::isnan(x)) { + const realTp sinh_re = x; + const realTp sinh_im = y; + return Tp{sinh_im, -sinh_re}; + } + const realTp sinh_re = x; + const realTp sinh_im = sycl::copysign(realTp(0), y); + return Tp{sinh_im, -sinh_re}; + } + + /* + * sinh(x +- I Inf) = dNaN + I dNaN. + * + * sinh(x + I NaN) = d(NaN) + I d(NaN). + */ + if (xfinite && !yfinite) { + const realTp sinh_re = q_nan; + const realTp sinh_im = x * sinh_re; + return Tp{sinh_im, -sinh_re}; + } + + /* + * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). + * The sign of Inf in the result is unspecified. Choice = normally + * the same as d(NaN). + * + * sinh(+-Inf +- I Inf) = +Inf + I dNaN. + * The sign of Inf in the result is unspecified. + * Choice = always - here for sinh to have positive result for + * imaginary part of sin. + * + * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y) + */ + if (std::isinf(x)) { + if (!yfinite) { + const realTp sinh_re = -x * x; + const realTp sinh_im = x * (y - y); + return Tp{sinh_im, -sinh_re}; + } + const realTp sinh_re = x * sycl::cos(y); + const realTp sinh_im = + std::numeric_limits::infinity() * sycl::sin(y); + return Tp{sinh_im, -sinh_re}; + } + + /* + * sinh(NaN + I NaN) = d(NaN) + I d(NaN). + * + * sinh(NaN +- I Inf) = d(NaN) + I d(NaN). + * + * sinh(NaN + I y) = d(NaN) + I d(NaN). + */ + const realTp y_m_y = (y - y); + const realTp sinh_re = (x * x) * y_m_y; + const realTp sinh_im = (x + x) * y_m_y; + return Tp{sinh_im, -sinh_re}; + } + else { + if (in == Tp(0)) { + return in; + } + return sycl::sin(in); + } +} +} // namespace impl + +template struct SincFunctor { // is function constant for given argT using is_constant = typename std::false_type; // constant value, if constant - // constexpr resT constant_value = resT{}; + // constexpr Tp constant_value = Tp{}; // is function defined for sycl::vec using supports_vec = typename std::false_type; - // do both argT and resT support subgroup store/load operation - using supports_sg_loadstore = typename std::true_type; + // do both argT and Tp support subgroup store/load operation + using supports_sg_loadstore = typename std::negation< + std::disjunction, tu_ns::is_complex>>; - resT operator()(const argT &x) const + Tp operator()(const argT &x) const { constexpr argT pi = static_cast(3.1415926535897932384626433832795029L); const argT y = pi * x; if (y == argT(0)) { - return argT(1); + return Tp(1); + } + + if constexpr (tu_ns::is_complex::value) { + return impl::sin(y) / y; + } + else { + return sycl::sinpi(x) / y; } - return sycl::sinpi(x) / y; } }; } // namespace dpnp::kernels::sinc From 673b22a8c8204f9a22df8d0e42a7ecd3a9313af1 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Mon, 28 Oct 2024 13:39:26 +0100 Subject: [PATCH 05/11] Mute tests for dono.i0 --- tests/third_party/cupy/math_tests/test_special.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/third_party/cupy/math_tests/test_special.py b/tests/third_party/cupy/math_tests/test_special.py index d257e2f98e4d..5b3b84c29a5e 100644 --- a/tests/third_party/cupy/math_tests/test_special.py +++ b/tests/third_party/cupy/math_tests/test_special.py @@ -8,6 +8,7 @@ class TestSpecial(unittest.TestCase): + @pytest.mark.skip("i0 is not implemented") @testing.for_dtypes(["e", "f", "d"]) @testing.numpy_cupy_allclose(rtol=1e-3) def test_i0(self, xp, dtype): From 4bd7f3a616696b9e45a46481c230c354793f7b48 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Mon, 28 Oct 2024 15:50:12 +0100 Subject: [PATCH 06/11] Add CFD tests --- tests/test_sycl_queue.py | 1 + tests/test_usm_type.py | 1 + 2 files changed, 2 insertions(+) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index e8c0969f8d0c..2dfdd3c1dd1a 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -502,6 +502,7 @@ def test_meshgrid(device): pytest.param( "sin", [-dpnp.pi / 2, -dpnp.pi / 4, 0.0, dpnp.pi / 4, dpnp.pi / 2] ), + pytest.param("sinc", [-5.0, -3.5, 0.0, 2.5, 4.3]), pytest.param("sinh", [-5.0, -3.5, 0.0, 3.5, 5.0]), pytest.param("sort", [2.0, 1.0, 7.0, 4.0]), pytest.param("sort_complex", [1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j]), diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 13752a377dbe..33d0396541b3 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -641,6 +641,7 @@ def test_norm(usm_type, ord, axis): pytest.param( "sin", [-dp.pi / 2, -dp.pi / 4, 0.0, dp.pi / 4, dp.pi / 2] ), + pytest.param("sinc", [-5.0, -3.5, 0.0, 2.5, 4.3]), pytest.param("sinh", [-5.0, -3.5, 0.0, 3.5, 5.0]), pytest.param("sort", [2.0, 1.0, 7.0, 4.0]), pytest.param("sort_complex", [1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j]), From 1f89f5a816d484ae6c532599e89a0070119f82d4 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Mon, 28 Oct 2024 19:25:20 +0100 Subject: [PATCH 07/11] Add more tests to cover different use cases --- tests/test_mathematical.py | 87 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/tests/test_mathematical.py b/tests/test_mathematical.py index 01555df2718e..49e6332f38c2 100644 --- a/tests/test_mathematical.py +++ b/tests/test_mathematical.py @@ -1897,6 +1897,93 @@ def test_wrong_tol_type(self, xp, tol_val): assert_raises(TypeError, xp.real_if_close, a, tol=tol_val) +class TestSinc: + @pytest.mark.parametrize( + "dt", get_all_dtypes(no_none=True, no_bool=True, no_float16=False) + ) + def test_basic(self, dt): + a = numpy.linspace(-1, 1, 100, dtype=dt) + ia = dpnp.array(a) + + result = dpnp.sinc(ia) + expected = numpy.sinc(a) + assert_dtype_allclose(result, expected) + + def test_bool(self): + a = numpy.array([True, False, True]) + ia = dpnp.array(a) + + result = dpnp.sinc(ia) + expected = numpy.sinc(a) + # numpy 1.26 promotes result to float64 dtype, but expected float16 + assert_dtype_allclose( + result, + expected, + check_only_type_kind=( + numpy.lib.NumpyVersion(numpy.__version__) < "2.0.0" + ), + ) + + @pytest.mark.parametrize("dt", get_all_dtypes(no_none=True, no_bool=True)) + def test_zero(self, dt): + if ( + dt == numpy.float16 + and numpy.lib.NumpyVersion(numpy.__version__) < "2.0.0" + ): + pytest.skip("numpy.sinc return NaN") + + a = numpy.array([0.0], dtype=dt) + ia = dpnp.array(a) + + result = dpnp.sinc(ia) + expected = numpy.sinc(a) + assert_dtype_allclose(result, expected) + + @testing.with_requires("numpy>=2.0.0") + def test_zero_fp16(self): + a = numpy.array([0.0], dtype=numpy.float16) + ia = dpnp.array(a) + + result = dpnp.sinc(ia) + expected = numpy.sinc(a) + assert_dtype_allclose(result, expected) + + @pytest.mark.usefixtures("suppress_invalid_numpy_warnings") + def test_nan_infs(self): + a = numpy.array([numpy.inf, -numpy.inf, numpy.nan]) + ia = dpnp.array(a) + + result = dpnp.sinc(ia) + expected = numpy.sinc(a) + assert_equal(result, expected) + + @pytest.mark.usefixtures("suppress_invalid_numpy_warnings") + def test_nan_infs_complex(self): + a = numpy.array( + [ + numpy.inf, + -numpy.inf, + numpy.nan, + complex(numpy.nan), + complex(numpy.nan, numpy.nan), + complex(0, numpy.nan), + complex(numpy.inf, numpy.nan), + complex(numpy.nan, numpy.inf), + complex(-numpy.inf, numpy.nan), + complex(numpy.nan, -numpy.inf), + complex(numpy.inf, numpy.inf), + complex(numpy.inf, -numpy.inf), + complex(-numpy.inf, numpy.inf), + complex(-numpy.inf, -numpy.inf), + ] + ) + ia = dpnp.array(a) + + result = dpnp.sinc(ia) + expected = numpy.sinc(a) + assert_equal(result, expected) + + class TestSpacing: @pytest.mark.parametrize("sign", [1, -1]) @pytest.mark.parametrize("dt", get_float_dtypes()) From 8f54def424a0f803d6fd7bb7a7e95c64c534a0cc Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Tue, 29 Oct 2024 11:09:07 +0100 Subject: [PATCH 08/11] Add a proper handling of corner cases --- tests/test_mathematical.py | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/tests/test_mathematical.py b/tests/test_mathematical.py index 49e6332f38c2..ba74c8db1e30 100644 --- a/tests/test_mathematical.py +++ b/tests/test_mathematical.py @@ -1915,23 +1915,11 @@ def test_bool(self): result = dpnp.sinc(ia) expected = numpy.sinc(a) - # numpy 1.26 promotes result to float64 dtype, but expected float16 - assert_dtype_allclose( - result, - expected, - check_only_type_kind=( - numpy.lib.NumpyVersion(numpy.__version__) < "2.0.0" - ), - ) + # numpy promotes result to float64 dtype, but expected float16 + assert_dtype_allclose(result, expected, check_only_type_kind=True) @pytest.mark.parametrize("dt", get_all_dtypes(no_none=True, no_bool=True)) def test_zero(self, dt): - if ( - dt == numpy.float16 - and numpy.lib.NumpyVersion(numpy.__version__) < "2.0.0" - ): - pytest.skip("numpy.sinc return NaN") - a = numpy.array([0.0], dtype=dt) ia = dpnp.array(a) @@ -1939,13 +1927,15 @@ def test_zero(self, dt): expected = numpy.sinc(a) assert_dtype_allclose(result, expected) + # TODO: add a proper NumPY version once resolved @testing.with_requires("numpy>=2.0.0") def test_zero_fp16(self): a = numpy.array([0.0], dtype=numpy.float16) ia = dpnp.array(a) result = dpnp.sinc(ia) - expected = numpy.sinc(a) + # expected = numpy.sinc(a) # numpy returns NaN, but expected 1.0 + expected = numpy.ones_like(a) assert_dtype_allclose(result, expected) @pytest.mark.usefixtures("suppress_invalid_numpy_warnings") From 7a6501fd6ec2b652cdd3574507aae226051fa896 Mon Sep 17 00:00:00 2001 From: Anton <100830759+antonwolfy@users.noreply.github.com> Date: Fri, 1 Nov 2024 11:39:22 +0100 Subject: [PATCH 09/11] Update dpnp/dpnp_iface_mathematical.py Co-authored-by: vtavana <120411540+vtavana@users.noreply.github.com> --- dpnp/dpnp_iface_mathematical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index c0b4b2aa8470..1605b886c937 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -3799,7 +3799,7 @@ def real_if_close(a, tol=100): Examples -------- >>> import dpnp as np ->>> x = np.linspace(-4, 4, 41, dtype=dpnp.float64) +>>> x = np.linspace(-4, 4, 41, dtype=np.float64) >>> np.sinc(x) # result may vary array([ 0. , -0.04923628, -0.08409186, -0.08903844, -0.05846808, 0. , 0.06682066, 0.11643488, 0.12613779, 0.08504448, From 6f028ea3532904ac590a98ccff1335f39aadbe09 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Fri, 1 Nov 2024 11:53:49 +0100 Subject: [PATCH 10/11] Removed limitation block --- dpnp/dpnp_iface_mathematical.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 1605b886c937..47c909217d3f 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -3779,12 +3779,6 @@ def real_if_close(a, tol=100): The data type of the returned array is determined by the Type Promotion Rules. -Limitations ------------ -Parameters `where` and `subok` are supported with their default values. -Keyword argument `kwargs` is currently unsupported. -Otherwise ``NotImplementedError`` exception will be raised. - Notes ----- The name sinc is short for "sine cardinal" or "sinus cardinalis". From a1ad4d9e9e3ee5fae19b3359404d6f6981a6f697 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Mon, 4 Nov 2024 11:38:08 +0100 Subject: [PATCH 11/11] Use dedicated DPNPSinc class for dpnp.sinc implementation --- dpnp/dpnp_algo/dpnp_elementwise_common.py | 22 ++++++++++++++++++++++ dpnp/dpnp_iface_mathematical.py | 3 ++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/dpnp/dpnp_algo/dpnp_elementwise_common.py b/dpnp/dpnp_algo/dpnp_elementwise_common.py index 0f2b383751ff..7c627867b694 100644 --- a/dpnp/dpnp_algo/dpnp_elementwise_common.py +++ b/dpnp/dpnp_algo/dpnp_elementwise_common.py @@ -48,6 +48,7 @@ "DPNPImag", "DPNPReal", "DPNPRound", + "DPNPSinc", "DPNPUnaryFunc", "resolve_weak_types_2nd_arg_int", ] @@ -584,6 +585,27 @@ def __call__(self, x, decimals=0, out=None, dtype=None): return super().__call__(x, out=out, dtype=dtype) +class DPNPSinc(DPNPUnaryFunc): + """Class that implements dpnp.sinc unary element-wise functions.""" + + def __init__( + self, + name, + result_type_resolver_fn, + unary_dp_impl_fn, + docs, + ): + super().__init__( + name, + result_type_resolver_fn, + unary_dp_impl_fn, + docs, + ) + + def __call__(self, x, out=None, order="K"): + return super().__call__(x, out=out, order=order) + + def acceptance_fn_gcd_lcm( arg1_dtype, arg2_dtype, buf1_dt, buf2_dt, res_dt, sycl_dev ): diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 47c909217d3f..9318f0af4dc3 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -65,6 +65,7 @@ DPNPImag, DPNPReal, DPNPRound, + DPNPSinc, DPNPUnaryFunc, acceptance_fn_gcd_lcm, acceptance_fn_negative, @@ -3806,7 +3807,7 @@ def real_if_close(a, tol=100): 0. ]) """ -sinc = DPNPUnaryFunc( +sinc = DPNPSinc( "sinc", ufi._sinc_result_type, ufi._sinc,