From a24d367ce24a8e9d762efe10ff8e52ce1ac48acf Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 25 Mar 2025 04:03:34 -0700 Subject: [PATCH 01/26] Initial impl of dpnp.inter() --- dpnp/CMakeLists.txt | 1 + dpnp/backend/extensions/math/CMakeLists.txt | 87 +++++ dpnp/backend/extensions/math/common.cpp | 75 ++++ dpnp/backend/extensions/math/common.hpp | 93 +++++ .../extensions/math/dispatch_table.hpp | 292 +++++++++++++++ dpnp/backend/extensions/math/interpolate.cpp | 338 ++++++++++++++++++ dpnp/backend/extensions/math/interpolate.hpp | 66 ++++ dpnp/backend/extensions/math/math_py.cpp | 37 ++ dpnp/dpnp_iface_statistics.py | 175 +++++++++ 9 files changed, 1164 insertions(+) create mode 100644 dpnp/backend/extensions/math/CMakeLists.txt create mode 100644 dpnp/backend/extensions/math/common.cpp create mode 100644 dpnp/backend/extensions/math/common.hpp create mode 100644 dpnp/backend/extensions/math/dispatch_table.hpp create mode 100644 dpnp/backend/extensions/math/interpolate.cpp create mode 100644 dpnp/backend/extensions/math/interpolate.hpp create mode 100644 dpnp/backend/extensions/math/math_py.cpp diff --git a/dpnp/CMakeLists.txt b/dpnp/CMakeLists.txt index 6be90d849dc4..6c59141bd1d9 100644 --- a/dpnp/CMakeLists.txt +++ b/dpnp/CMakeLists.txt @@ -60,6 +60,7 @@ add_subdirectory(backend/extensions/blas) add_subdirectory(backend/extensions/fft) add_subdirectory(backend/extensions/indexing) add_subdirectory(backend/extensions/lapack) +add_subdirectory(backend/extensions/math) add_subdirectory(backend/extensions/statistics) add_subdirectory(backend/extensions/ufunc) add_subdirectory(backend/extensions/vm) diff --git a/dpnp/backend/extensions/math/CMakeLists.txt b/dpnp/backend/extensions/math/CMakeLists.txt new file mode 100644 index 000000000000..fed91cfd3f9a --- /dev/null +++ b/dpnp/backend/extensions/math/CMakeLists.txt @@ -0,0 +1,87 @@ +# ***************************************************************************** +# Copyright (c) 2016-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. +# ***************************************************************************** + + +set(python_module_name _math_impl) +set(_module_src + ${CMAKE_CURRENT_SOURCE_DIR}/common.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/interpolate.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/math_py.cpp +) + +pybind11_add_module(${python_module_name} MODULE ${_module_src}) +add_sycl_to_target(TARGET ${python_module_name} SOURCES ${_module_src}) + +if(_dpnp_sycl_targets) + # make fat binary + target_compile_options( + ${python_module_name} + PRIVATE + -fsycl-targets=${_dpnp_sycl_targets} + ) + target_link_options( + ${python_module_name} + PRIVATE + -fsycl-targets=${_dpnp_sycl_targets} + ) +endif() + +if (WIN32) + if (${CMAKE_VERSION} VERSION_LESS "3.27") + # this is a work-around for target_link_options inserting option after -link option, cause + # linker to ignore it. + set(CMAKE_CXX_LINK_FLAGS "${CMAKE_CXX_LINK_FLAGS} -fsycl-device-code-split=per_kernel") + endif() +endif() + +set_target_properties(${python_module_name} PROPERTIES CMAKE_POSITION_INDEPENDENT_CODE ON) + +target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../include) +target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../src) + +target_include_directories(${python_module_name} PUBLIC ${Dpctl_INCLUDE_DIR}) +target_include_directories(${python_module_name} PUBLIC ${Dpctl_TENSOR_INCLUDE_DIR}) + +if (WIN32) + target_compile_options(${python_module_name} PRIVATE + /clang:-fno-approx-func + /clang:-fno-finite-math-only + ) +else() + target_compile_options(${python_module_name} PRIVATE + -fno-approx-func + -fno-finite-math-only + ) +endif() + +target_link_options(${python_module_name} PUBLIC -fsycl-device-code-split=per_kernel) + +if (DPNP_GENERATE_COVERAGE) + target_link_options(${python_module_name} PRIVATE -fprofile-instr-generate -fcoverage-mapping) +endif() + +install(TARGETS ${python_module_name} + DESTINATION "dpnp/backend/extensions/math" +) diff --git a/dpnp/backend/extensions/math/common.cpp b/dpnp/backend/extensions/math/common.cpp new file mode 100644 index 000000000000..93723e838b21 --- /dev/null +++ b/dpnp/backend/extensions/math/common.cpp @@ -0,0 +1,75 @@ +//***************************************************************************** +// 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 "common.hpp" +#include "utils/type_dispatch.hpp" +#include + +namespace dpctl_td_ns = dpctl::tensor::type_dispatch; + +namespace math +{ +namespace common +{ +pybind11::dtype dtype_from_typenum(int dst_typenum) +{ + dpctl_td_ns::typenum_t dst_typenum_t = + static_cast(dst_typenum); + switch (dst_typenum_t) { + case dpctl_td_ns::typenum_t::BOOL: + return py::dtype("?"); + case dpctl_td_ns::typenum_t::INT8: + return py::dtype("i1"); + case dpctl_td_ns::typenum_t::UINT8: + return py::dtype("u1"); + case dpctl_td_ns::typenum_t::INT16: + return py::dtype("i2"); + case dpctl_td_ns::typenum_t::UINT16: + return py::dtype("u2"); + case dpctl_td_ns::typenum_t::INT32: + return py::dtype("i4"); + case dpctl_td_ns::typenum_t::UINT32: + return py::dtype("u4"); + case dpctl_td_ns::typenum_t::INT64: + return py::dtype("i8"); + case dpctl_td_ns::typenum_t::UINT64: + return py::dtype("u8"); + case dpctl_td_ns::typenum_t::HALF: + return py::dtype("f2"); + case dpctl_td_ns::typenum_t::FLOAT: + return py::dtype("f4"); + case dpctl_td_ns::typenum_t::DOUBLE: + return py::dtype("f8"); + case dpctl_td_ns::typenum_t::CFLOAT: + return py::dtype("c8"); + case dpctl_td_ns::typenum_t::CDOUBLE: + return py::dtype("c16"); + default: + throw py::value_error("Unrecognized dst_typeid"); + } +} + +} // namespace common +} // namespace math diff --git a/dpnp/backend/extensions/math/common.hpp b/dpnp/backend/extensions/math/common.hpp new file mode 100644 index 000000000000..1d436440c392 --- /dev/null +++ b/dpnp/backend/extensions/math/common.hpp @@ -0,0 +1,93 @@ +//***************************************************************************** +// 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 +#include +#include + +// clang-format off +// math_utils.hpp doesn't include sycl header but uses sycl types +// so sycl.hpp must be included before math_utils.hpp +#include +#include "utils/math_utils.hpp" +// clang-format on + +namespace math +{ +namespace common +{ +template +struct Less +{ + bool operator()(const T &lhs, const T &rhs) const + { + return std::less{}(lhs, rhs); + } +}; + +template +struct Less> +{ + bool operator()(const std::complex &lhs, + const std::complex &rhs) const + { + return dpctl::tensor::math_utils::less_complex(lhs, rhs); + } +}; + +template +struct IsNan +{ + static bool isnan(const T &v) + { + if constexpr (std::is_floating_point_v || + std::is_same_v) { + return sycl::isnan(v); + } + + return false; + } +}; + +template +struct IsNan> +{ + static bool isnan(const std::complex &v) + { + T real1 = std::real(v); + T imag1 = std::imag(v); + return sycl::isnan(real1) || sycl::isnan(imag1); + } +}; + + +// This function is a copy from dpctl because it is not available in the public +// headers of dpctl. +pybind11::dtype dtype_from_typenum(int dst_typenum); + +} // namespace common +} // namespace math diff --git a/dpnp/backend/extensions/math/dispatch_table.hpp b/dpnp/backend/extensions/math/dispatch_table.hpp new file mode 100644 index 000000000000..4cfd3d2a09a4 --- /dev/null +++ b/dpnp/backend/extensions/math/dispatch_table.hpp @@ -0,0 +1,292 @@ +//***************************************************************************** +// 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 +#include + +#include "utils/type_dispatch.hpp" +#include +#include +#include +#include + +#include "common.hpp" + +namespace dpctl_td_ns = dpctl::tensor::type_dispatch; +namespace py = pybind11; + +namespace math +{ +namespace common +{ + +template +struct one_of +{ + static_assert(std::is_same_v>, + "one_of: second parameter cannot be empty std::tuple"); + static_assert(false, "one_of: second parameter must be std::tuple"); +}; + +template +struct one_of> +{ + static constexpr bool value = + std::is_same_v || one_of>::value; +}; + +template +struct one_of> +{ + static constexpr bool value = std::is_same_v; +}; + +template +constexpr bool one_of_v = one_of::value; + +template +using Table = FnT[dpctl_td_ns::num_types]; +template +using Table2 = Table[dpctl_td_ns::num_types]; + +using TypeId = int32_t; +using TypesPair = std::pair; + +struct int_pair_hash +{ + inline size_t operator()(const TypesPair &p) const + { + std::hash hasher; + return hasher(size_t(p.first) << (8 * sizeof(TypeId)) | + size_t(p.second)); + } +}; + +using SupportedTypesList = std::vector; +using SupportedTypesList2 = std::vector; +using SupportedTypesSet = std::unordered_set; +using SupportedTypesSet2 = std::unordered_set; + +using DType = py::dtype; +using DTypePair = std::pair; + +using SupportedDTypeList = std::vector; +using SupportedDTypeList2 = std::vector; + +template + typename Func> +struct TableBuilder2 +{ + template + struct impl + { + static constexpr bool is_defined = + one_of_v, SupportedTypes>; + + _FnT get() + { + if constexpr (is_defined) { + return Func::impl; + } + else { + return nullptr; + } + } + }; + + using type = + dpctl_td_ns::DispatchTableBuilder; +}; + +template +class DispatchTable2 +{ +public: + DispatchTable2(std::string first_name, std::string second_name) + : first_name(first_name), second_name(second_name) + { + } + + template + typename Func> + void populate_dispatch_table() + { + using TBulder = typename TableBuilder2::type; + TBulder builder; + + builder.populate_dispatch_table(table); + populate_supported_types(); + } + + FnT get_unsafe(int first_typenum, int second_typenum) const + { + auto array_types = dpctl_td_ns::usm_ndarray_types(); + const int first_type_id = + array_types.typenum_to_lookup_id(first_typenum); + const int second_type_id = + array_types.typenum_to_lookup_id(second_typenum); + + return table[first_type_id][second_type_id]; + } + + FnT get(int first_typenum, int second_typenum) const + { + auto fn = get_unsafe(first_typenum, second_typenum); + + if (fn == nullptr) { + auto array_types = dpctl_td_ns::usm_ndarray_types(); + const int first_type_id = + array_types.typenum_to_lookup_id(first_typenum); + const int second_type_id = + array_types.typenum_to_lookup_id(second_typenum); + + py::dtype first_dtype = dtype_from_typenum(first_type_id); + auto first_type_pos = + std::find(supported_first_type.begin(), + supported_first_type.end(), first_dtype); + if (first_type_pos == supported_first_type.end()) { + py::str types = py::str(py::cast(supported_first_type)); + py::str dtype = py::str(first_dtype); + + py::str err_msg = + py::str("'" + first_name + "' has unsupported type '") + + dtype + + py::str("'." + " Supported types are: ") + + types; + + throw py::value_error(static_cast(err_msg)); + } + + py::dtype second_dtype = dtype_from_typenum(second_type_id); + auto second_type_pos = + std::find(supported_second_type.begin(), + supported_second_type.end(), second_dtype); + if (second_type_pos == supported_second_type.end()) { + py::str types = py::str(py::cast(supported_second_type)); + py::str dtype = py::str(second_dtype); + + py::str err_msg = + py::str("'" + second_name + "' has unsupported type '") + + dtype + + py::str("'." + " Supported types are: ") + + types; + + throw py::value_error(static_cast(err_msg)); + } + + py::str first_dtype_str = py::str(first_dtype); + py::str second_dtype_str = py::str(second_dtype); + py::str types = py::str(py::cast(all_supported_types)); + + py::str err_msg = + py::str("'" + first_name + "' and '" + second_name + + "' has unsupported types combination: ('") + + first_dtype_str + py::str("', '") + second_dtype_str + + py::str("')." + " Supported types combinations are: ") + + types; + + throw py::value_error(static_cast(err_msg)); + } + + return fn; + } + + const SupportedDTypeList &get_supported_first_type() const + { + return supported_first_type; + } + + const SupportedDTypeList &get_supported_second_type() const + { + return supported_second_type; + } + + const SupportedDTypeList2 &get_all_supported_types() const + { + return all_supported_types; + } + +private: + void populate_supported_types() + { + SupportedTypesSet first_supported_types_set; + SupportedTypesSet second_supported_types_set; + SupportedTypesSet2 all_supported_types_set; + + for (int i = 0; i < dpctl_td_ns::num_types; ++i) { + for (int j = 0; j < dpctl_td_ns::num_types; ++j) { + if (table[i][j] != nullptr) { + all_supported_types_set.emplace(i, j); + first_supported_types_set.emplace(i); + second_supported_types_set.emplace(j); + } + } + } + + auto to_supported_dtype_list = [](const auto &supported_set, + auto &supported_list) { + SupportedTypesList lst(supported_set.begin(), supported_set.end()); + std::sort(lst.begin(), lst.end()); + supported_list.resize(supported_set.size()); + std::transform(lst.begin(), lst.end(), supported_list.begin(), + [](TypeId i) { return dtype_from_typenum(i); }); + }; + + to_supported_dtype_list(first_supported_types_set, + supported_first_type); + to_supported_dtype_list(second_supported_types_set, + supported_second_type); + + SupportedTypesList2 lst(all_supported_types_set.begin(), + all_supported_types_set.end()); + std::sort(lst.begin(), lst.end()); + all_supported_types.resize(all_supported_types_set.size()); + std::transform(lst.begin(), lst.end(), all_supported_types.begin(), + [](TypesPair p) { + return DTypePair(dtype_from_typenum(p.first), + dtype_from_typenum(p.second)); + }); + } + + std::string first_name; + std::string second_name; + + SupportedDTypeList supported_first_type; + SupportedDTypeList supported_second_type; + SupportedDTypeList2 all_supported_types; + + Table2 table; +}; + +} // namespace common +} // namespace math diff --git a/dpnp/backend/extensions/math/interpolate.cpp b/dpnp/backend/extensions/math/interpolate.cpp new file mode 100644 index 000000000000..9b69affc5ea0 --- /dev/null +++ b/dpnp/backend/extensions/math/interpolate.cpp @@ -0,0 +1,338 @@ +//***************************************************************************** +// 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 +#include +#include +#include +#include + +#include +#include + +// dpctl tensor headers +#include "dpctl4pybind11.hpp" +#include "utils/type_dispatch.hpp" + +#include "interpolate.hpp" + +namespace dpctl_td_ns = dpctl::tensor::type_dispatch; +using dpctl::tensor::usm_ndarray; + +using namespace math::interpolate; + +namespace +{ + +template +struct interpolate_kernel +{ + static sycl::event impl(sycl::queue &exec_q, + const void *vx, + const void *vidx, + const void *vxp, + const void *vfp, + void *vout, + const size_t xp_size, + const size_t n, + const std::vector &depends) + { + const T *x = static_cast(vx); + const T *xp = static_cast(vxp); + const T *fp = static_cast(vfp); + const IndexT *idx = static_cast(vidx); + T *out = static_cast(vout); + + // size_t n = x.get_size(); + + // std::vector result(n); + // sycl::event copy_ev = exec_q.memcpy(result.data(), output.get_data(), n * sizeof(float), {ev}); + + // for (size_t i = 0; i< n; i++){ + // std::cout << result[i] << " "; + // } + + return exec_q.submit([&](sycl::handler &cgh) { + cgh.depends_on(depends); + + cgh.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { + // T left = fp[0]; + // T right = fp[xp_size - 1]; + + // // IndexT x_idx = idx[i] - 1; + + // // if (sycl::isnan(x[i])) { + // // out[i] = x[i]; + // // } + // // else if (x_idx < 0) { + // // out[i] = left; + // // } + // // //old version check + // // // else if (x[i] == xp[xp_size - 1]) { + // // // out[i] = right; + // // // } + // // // else if (idx[i] >= xp_size - 1) { + // // // out[i] = right; + // // // } + // // // new version check + // // else if (idx[i] == xp_size) { + // // out[i] = right; + // // } + // // else if (idx[i] == xp_size - 1) { + // // out[i] = fp[x_idx]; + // // } + // // else if (x[i] == xp[x_idx]) { + // // out[i] = fp[x_idx]; + // // } + + // IndexT j = idx[i]; + + // if (sycl::isnan(x[i])) { + // out[i] = x[i]; + // } + // else if (j == 0) { + // out[i] = left; + // } + // else if (j == xp_size) { + // out[i] = right; + // } + // else { + // IndexT x_idx = j - 1; + + // if (x[i] == xp[x_idx]) { + // out[i] = fp[x_idx]; + // } + // else { + // T slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); + // T res = slope * (x[i] - xp[x_idx]) + fp[x_idx]; + // // T res = (x[i] - xp[x_idx]) + fp[x_idx]; + + // if (sycl::isnan(res)) { + // res = slope * (x[i] - xp[x_idx + 1]) + fp[x_idx + 1]; + // if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { + // res = fp[x_idx]; + // } + // } + // out[i] = res; + // } + // } + + T left = fp[0]; + T right = fp[xp_size - 1]; + IndexT x_idx = idx[i] - 1; + + if (sycl::isnan(x[i])) { + out[i] = x[i]; + } + else if (x_idx < 0) { + out[i] = left; + } + else if (x[i] == xp[xp_size - 1]) { + out[i] = right; + } + else if (x_idx >= xp_size - 1) { + out[i] = right; + } + else if (x[i] == xp[x_idx]) { + out[i] = fp[x_idx]; + } + else { + T slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); + T res = slope * (x[i] - xp[x_idx]) + fp[x_idx]; + + if (sycl::isnan(res)) { + res = slope * (x[i] - xp[x_idx + 1]) + fp[x_idx + 1]; + if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { + res = fp[x_idx]; + } + } + + out[i] = res; + } + + // out[i] = x[i]; + }); + }); + } +}; + +// using SupportedTypes = std::tuple, std::tuple>; +// using SupportedTypes = std::tuple, +// std::tuple, +// std::tuple, +// std::tuple, +// std::tuple, +// std::tuple, +// std::tuple>, +// std::tuple>, +// std::tuple>, +// std::tuple>, +// std::tuple, +// std::tuple, +// std::tuple, +// std::tuple, +// std::tuple>, +// std::tuple>, +// std::tuple, int64_t>, +// std::tuple, int64_t>, +// std::tuple, float>, +// std::tuple, double>>; + +using SupportedTypes = std::tuple< + // std::tuple, + // std::tuple, + std::tuple, + std::tuple>; + +} // namespace + +Interpolate::Interpolate() : dispatch_table("x", "idx") +{ + dispatch_table.populate_dispatch_table(); +} + +std::tuple +Interpolate::call(const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + const size_t xp_size, + dpctl::tensor::usm_ndarray &output, + const std::vector &depends) +{ + // validate(x, xp, fp, output); + + if (x.get_size() == 0) { + return {sycl::event(), sycl::event()}; + } + + const int x_typenum = x.get_typenum(); + const int idx_typenum = idx.get_typenum(); + + auto interp_func = dispatch_table.get(x_typenum, idx_typenum); + + auto exec_q = x.get_queue(); + + // size_t n = x.get_size(); + // const size_t m = xp.get_size(); + + // std::vector x_h(n); + // std::vector xp_h(n); + // std::vector fp_h(n); + // std::vector output_h(n); + + // sycl::event copy_1 = exec_q.memcpy(x_h.data(), x.get_data(), n * sizeof(float)); + // sycl::event copy_2 = exec_q.memcpy(xp_h.data(), xp.get_data(), m * sizeof(float)); + // sycl::event copy_3 = exec_q.memcpy(fp_h.data(), fp.get_data(), m * sizeof(float)); + // sycl::event copy_4 = exec_q.memcpy(output_h.data(), output.get_data(), n * sizeof(float)); + + // copy_1.wait(); + // copy_2.wait(); + // copy_3.wait(); + // copy_4.wait(); + + // std::cout << "x: " << std::endl; + // for (size_t i = 0; i< n; i++){ + // std::cout << x_h[i] << " "; + // } + // std::cout << "\n"; + + // std::cout << "xp: " << std::endl; + // for (size_t i = 0; i< m; i++){ + // std::cout << xp_h[i] << " "; + // } + // std::cout << "\n"; + + + // std::cout << "fp: " << std::endl; + // for (size_t i = 0; i< m; i++){ + // std::cout << fp_h[i] << " "; + // } + // std::cout << "\n"; + + + // std::cout << "out: " << std::endl; + // for (size_t i = 0; i< n; i++){ + // std::cout << output_h[i] << " "; + // } + // std::cout << "\n"; + + auto ev = + interp_func(exec_q, x.get_data(), idx.get_data(), xp.get_data(), fp.get_data(), + output.get_data(), xp.get_size(), x.get_size(), depends); + + ev.wait(); + + auto args_ev = dpctl::utils::keep_args_alive( + exec_q, {x, idx, xp, fp, output}, {ev}); + + // size_t n = x.get_size(); + + // std::vector result(n); + // sycl::event copy_ev = exec_q.memcpy(result.data(), output.get_data(), n * sizeof(float), {ev}); + + // copy_ev.wait(); + + // std::cout << "out_host: " << std::endl; + // for (size_t i = 0; i< n; i++){ + // std::cout << result[i] << " "; + // } + + // std::cout << "\n"; + + return {args_ev, ev}; +} + +std::unique_ptr interp; + +void math::interpolate::populate_interpolate(py::module_ m) +{ + using namespace std::placeholders; + + interp.reset(new Interpolate()); + + auto interp_func = + [interpp = interp.get()]( + const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + const size_t xp_size, + dpctl::tensor::usm_ndarray &output, + const std::vector &depends) { + return interpp->call(x, idx, xp, fp, xp_size, output, depends); + }; + + m.def("interpolate", interp_func, "Perform linear interpolation.", + py::arg("x"), py::arg("idx"), py::arg("xp"), py::arg("fp"), + py::arg("xp_size"), py::arg("output"), py::arg("depends") = py::list()); + + auto interpolate_dtypes = [interpp = interp.get()]() { + return interpp->dispatch_table.get_all_supported_types(); + }; + + m.def("interpolate_dtypes", interpolate_dtypes, + "Get the supported data types for interpolation."); +} diff --git a/dpnp/backend/extensions/math/interpolate.hpp b/dpnp/backend/extensions/math/interpolate.hpp new file mode 100644 index 000000000000..95c34d945dad --- /dev/null +++ b/dpnp/backend/extensions/math/interpolate.hpp @@ -0,0 +1,66 @@ +//***************************************************************************** +// 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 + +#include "dispatch_table.hpp" + +namespace dpctl_td_ns = dpctl::tensor::type_dispatch; + +namespace math +{ +namespace interpolate +{ +struct Interpolate +{ + using FnT = sycl::event (*)(sycl::queue &, + const void *, + const void *, + const void *, + const void *, + void *, + const size_t, + const size_t, + const std::vector &); + + common::DispatchTable2 dispatch_table; + + Interpolate(); + + std::tuple + call(const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + const size_t xp_size, + dpctl::tensor::usm_ndarray &output, + const std::vector &depends); +}; + +void populate_interpolate(py::module_ m); +} // namespace interpolate +} // namespace math diff --git a/dpnp/backend/extensions/math/math_py.cpp b/dpnp/backend/extensions/math/math_py.cpp new file mode 100644 index 000000000000..6b8b5b9b5f28 --- /dev/null +++ b/dpnp/backend/extensions/math/math_py.cpp @@ -0,0 +1,37 @@ +//***************************************************************************** +// 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. +//***************************************************************************** +// +// This file defines functions of dpnp.backend._math_impl extensions +// +//***************************************************************************** + +#include + +#include "interpolate.hpp" + +PYBIND11_MODULE(_math_impl, m) +{ + math::interpolate::populate_interpolate(m); +} diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 3958127789ab..08b830b1d360 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -54,6 +54,8 @@ to_supported_dtypes, ) +import dpnp.backend.extensions.math._math_impl as math_ext + from .dpnp_utils import get_usm_allocations from .dpnp_utils.dpnp_utils_reduction import dpnp_wrap_reduction_call from .dpnp_utils.dpnp_utils_statistics import dpnp_cov, dpnp_median @@ -66,6 +68,7 @@ "corrcoef", "correlate", "cov", + "interp", "max", "mean", "median", @@ -1051,6 +1054,178 @@ def cov( ) +def interp(x, xp, fp, left=None, right=None, period=None): + """ + One-dimensional linear interpolation for monotonically increasing sample points. + + Returns the one-dimensional piecewise linear interpolant to a function + with given discrete data points (`xp`, `fp`), evaluated at `x`. + + Parameters + ---------- + x : array_like + The x-coordinates at which to evaluate the interpolated values. + + xp : 1-D sequence of floats + The x-coordinates of the data points, must be increasing if argument + `period` is not specified. Otherwise, `xp` is internally sorted after + normalizing the periodic boundaries with ``xp = xp % period``. + + fp : 1-D sequence of float or complex + The y-coordinates of the data points, same length as `xp`. + + left : optional float or complex corresponding to fp + Value to return for `x < xp[0]`, default is `fp[0]`. + + right : optional float or complex corresponding to fp + Value to return for `x > xp[-1]`, default is `fp[-1]`. + + period : None or float, optional + A period for the x-coordinates. This parameter allows the proper + interpolation of angular x-coordinates. Parameters `left` and `right` + are ignored if `period` is specified. + + Returns + ------- + y : float or complex (corresponding to fp) or ndarray + The interpolated values, same shape as `x`. + + Raises + ------ + ValueError + If `xp` and `fp` have different length + If `xp` or `fp` are not 1-D sequences + If `period == 0` + + See Also + -------- + scipy.interpolate + + Warnings + -------- + The x-coordinate sequence is expected to be increasing, but this is not + explicitly enforced. However, if the sequence `xp` is non-increasing, + interpolation results are meaningless. + + Note that, since NaN is unsortable, `xp` also cannot contain NaNs. + + A simple check for `xp` being strictly increasing is:: + + np.all(np.diff(xp) > 0) + + Examples + -------- + >>> import dpnp as np + >>> xp = np.array([1, 2, 3]) + >>> fp = np.array([3 ,2 ,0]) + >>> x = np.array([2.5]) + >>> np.interp(2.5, xp, fp) + 1.0 + >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp) + array([3. , 3. , 2.5 , 0.56, 0. ]) + >>> UNDEF = -99.0 + >>> np.interp(3.14, xp, fp, right=UNDEF) + -99.0 + + Plot an interpolant to the sine function: + + >>> x = np.linspace(0, 2*np.pi, 10) + >>> y = np.sin(x) + >>> xvals = np.linspace(0, 2*np.pi, 50) + >>> yinterp = np.interp(xvals, x, y) + >>> import matplotlib.pyplot as plt + >>> plt.plot(x, y, 'o') + [] + >>> plt.plot(xvals, yinterp, '-x') + [] + >>> plt.show() + + Interpolation with periodic x-coordinates: + + >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] + >>> xp = [190, -190, 350, -350] + >>> fp = [5, 10, 3, 4] + >>> np.interp(x, xp, fp, period=360) + array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75]) + + Complex interpolation: + + >>> x = [1.5, 4.0] + >>> xp = [2,3,5] + >>> fp = [1.0j, 0, 2+3j] + >>> np.interp(x, xp, fp) + array([0.+1.j , 1.+1.5j]) + + """ + + dpnp.check_supported_arrays_type(x, xp, fp) + + if xp.ndim != 1 or fp.ndim != 1: + raise ValueError('xp and fp must be 1D arrays') + if xp.size != fp.size: + raise ValueError('fp and xp are not of the same length') + if xp.size == 0: + raise ValueError('array of sample points is empty') + if not x.flags.c_contiguous: + raise NotImplementedError('Non-C-contiguous x is currently not ' + 'supported') + x_dtype = dpnp.common_type(x, xp) + if not dpnp.can_cast(x_dtype, dpnp.default_float_type()): + raise TypeError('Cannot cast array data from' + ' {} to {} according to the rule \'safe\'' + .format(x_dtype, dpnp.default_float_type())) + + if period is not None: + # The handling of "period" below is modified from NumPy's + + if period == 0: + raise ValueError("period must be a non-zero value") + period = dpnp.abs(period) + left = None + right = None + + x = x.astype(dpnp.default_float_type()) + xp = xp.astype(dpnp.default_float_type()) + + # normalizing periodic boundaries + x %= period + xp %= period + asort_xp = dpnp.argsort(xp) + xp = xp[asort_xp] + fp = fp[asort_xp] + xp = dpnp.concatenate((xp[-1:]-period, xp, xp[0:1]+period)) + fp = dpnp.concatenate((fp[-1:], fp, fp[0:1])) + assert xp.flags.c_contiguous + assert fp.flags.c_contiguous + + # NumPy always returns float64 or complex128, so we upcast all values + # on the fly in the kernel + out_dtype = 'f8' + output = dpnp.empty(x.shape, dtype=out_dtype) + idx = dpnp.searchsorted(xp, x, side='right') + left = fp[0] if left is None else dpnp.array(left, fp.dtype) + right = fp[-1] if right is None else dpnp.array(right, fp.dtype) + + idx = dpnp.array(idx, dtype='uint64') + + queue = x.sycl_queue + _manager = dpu.SequentialOrderManager[queue] + mem_ev, ht_ev = math_ext.interpolate( + x.get_array(), + idx.get_array(), + xp.get_array(), + fp.get_array(), + xp.size, + # left, + # right, + output.get_array(), + depends=_manager.submitted_events, + ) + _manager.add_event_pair(mem_ev, ht_ev) + + return output + + def max(a, axis=None, out=None, keepdims=False, initial=None, where=True): """ Return the maximum of an array or maximum along an axis. From e2b20b0bd6c2711b47983483f35361d172765730 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 2 Apr 2025 04:31:53 -0700 Subject: [PATCH 02/26] Second impl with dispatch_vector[only floating] --- dpnp/backend/extensions/math/interpolate.cpp | 322 +++--------------- dpnp/backend/extensions/math/interpolate.hpp | 41 +-- .../extensions/math/interpolate_kernel.hpp | 80 +++++ dpnp/backend/extensions/math/math_py.cpp | 2 +- dpnp/dpnp_iface_statistics.py | 6 +- .../third_party/cupy/test_type_routines.py | 4 + 6 files changed, 142 insertions(+), 313 deletions(-) create mode 100644 dpnp/backend/extensions/math/interpolate_kernel.hpp diff --git a/dpnp/backend/extensions/math/interpolate.cpp b/dpnp/backend/extensions/math/interpolate.cpp index 9b69affc5ea0..a620299bfffa 100644 --- a/dpnp/backend/extensions/math/interpolate.cpp +++ b/dpnp/backend/extensions/math/interpolate.cpp @@ -37,302 +37,78 @@ #include "utils/type_dispatch.hpp" #include "interpolate.hpp" +#include "interpolate_kernel.hpp" -namespace dpctl_td_ns = dpctl::tensor::type_dispatch; -using dpctl::tensor::usm_ndarray; +namespace dpnp::extensions::math +{ -using namespace math::interpolate; +namespace py = pybind11; +namespace td_ns = dpctl::tensor::type_dispatch; -namespace -{ +static kernels::interpolate_fn_ptr_t interpolate_dispatch_table[td_ns::num_types]; -template -struct interpolate_kernel +template +struct InterpolateFactory { - static sycl::event impl(sycl::queue &exec_q, - const void *vx, - const void *vidx, - const void *vxp, - const void *vfp, - void *vout, - const size_t xp_size, - const size_t n, - const std::vector &depends) + fnT get() { - const T *x = static_cast(vx); - const T *xp = static_cast(vxp); - const T *fp = static_cast(vfp); - const IndexT *idx = static_cast(vidx); - T *out = static_cast(vout); - - // size_t n = x.get_size(); - - // std::vector result(n); - // sycl::event copy_ev = exec_q.memcpy(result.data(), output.get_data(), n * sizeof(float), {ev}); - - // for (size_t i = 0; i< n; i++){ - // std::cout << result[i] << " "; - // } - - return exec_q.submit([&](sycl::handler &cgh) { - cgh.depends_on(depends); - - cgh.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { - // T left = fp[0]; - // T right = fp[xp_size - 1]; - - // // IndexT x_idx = idx[i] - 1; - - // // if (sycl::isnan(x[i])) { - // // out[i] = x[i]; - // // } - // // else if (x_idx < 0) { - // // out[i] = left; - // // } - // // //old version check - // // // else if (x[i] == xp[xp_size - 1]) { - // // // out[i] = right; - // // // } - // // // else if (idx[i] >= xp_size - 1) { - // // // out[i] = right; - // // // } - // // // new version check - // // else if (idx[i] == xp_size) { - // // out[i] = right; - // // } - // // else if (idx[i] == xp_size - 1) { - // // out[i] = fp[x_idx]; - // // } - // // else if (x[i] == xp[x_idx]) { - // // out[i] = fp[x_idx]; - // // } - - // IndexT j = idx[i]; - - // if (sycl::isnan(x[i])) { - // out[i] = x[i]; - // } - // else if (j == 0) { - // out[i] = left; - // } - // else if (j == xp_size) { - // out[i] = right; - // } - // else { - // IndexT x_idx = j - 1; - - // if (x[i] == xp[x_idx]) { - // out[i] = fp[x_idx]; - // } - // else { - // T slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); - // T res = slope * (x[i] - xp[x_idx]) + fp[x_idx]; - // // T res = (x[i] - xp[x_idx]) + fp[x_idx]; - - // if (sycl::isnan(res)) { - // res = slope * (x[i] - xp[x_idx + 1]) + fp[x_idx + 1]; - // if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { - // res = fp[x_idx]; - // } - // } - // out[i] = res; - // } - // } - - T left = fp[0]; - T right = fp[xp_size - 1]; - IndexT x_idx = idx[i] - 1; - - if (sycl::isnan(x[i])) { - out[i] = x[i]; - } - else if (x_idx < 0) { - out[i] = left; - } - else if (x[i] == xp[xp_size - 1]) { - out[i] = right; - } - else if (x_idx >= xp_size - 1) { - out[i] = right; - } - else if (x[i] == xp[x_idx]) { - out[i] = fp[x_idx]; - } - else { - T slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); - T res = slope * (x[i] - xp[x_idx]) + fp[x_idx]; - - if (sycl::isnan(res)) { - res = slope * (x[i] - xp[x_idx + 1]) + fp[x_idx + 1]; - if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { - res = fp[x_idx]; - } - } - - out[i] = res; - } - - // out[i] = x[i]; - }); - }); + if constexpr (std::is_floating_point_v) { + return kernels::interpolate_impl; + } + else { + return nullptr; + } } }; -// using SupportedTypes = std::tuple, std::tuple>; -// using SupportedTypes = std::tuple, -// std::tuple, -// std::tuple, -// std::tuple, -// std::tuple, -// std::tuple, -// std::tuple>, -// std::tuple>, -// std::tuple>, -// std::tuple>, -// std::tuple, -// std::tuple, -// std::tuple, -// std::tuple, -// std::tuple>, -// std::tuple>, -// std::tuple, int64_t>, -// std::tuple, int64_t>, -// std::tuple, float>, -// std::tuple, double>>; - -using SupportedTypes = std::tuple< - // std::tuple, - // std::tuple, - std::tuple, - std::tuple>; - -} // namespace - -Interpolate::Interpolate() : dispatch_table("x", "idx") -{ - dispatch_table.populate_dispatch_table(); -} -std::tuple -Interpolate::call(const dpctl::tensor::usm_ndarray &x, - const dpctl::tensor::usm_ndarray &idx, - const dpctl::tensor::usm_ndarray &xp, - const dpctl::tensor::usm_ndarray &fp, - const size_t xp_size, - dpctl::tensor::usm_ndarray &output, - const std::vector &depends) +std::pair +py_interpolate(const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + dpctl::tensor::usm_ndarray &out, + sycl::queue &exec_q, + const std::vector &depends) { - // validate(x, xp, fp, output); + int typenum = x.get_typenum(); + auto array_types = td_ns::usm_ndarray_types(); + int type_id = array_types.typenum_to_lookup_id(typenum); - if (x.get_size() == 0) { - return {sycl::event(), sycl::event()}; + auto fn = interpolate_dispatch_table[type_id]; + if (!fn) { + throw py::type_error("Unsupported dtype."); } - const int x_typenum = x.get_typenum(); - const int idx_typenum = idx.get_typenum(); + std::size_t n = x.get_size(); + std::size_t xp_size = xp.get_size(); - auto interp_func = dispatch_table.get(x_typenum, idx_typenum); + sycl::event ev = fn(exec_q, x.get_data(), idx.get_data(), xp.get_data(), + fp.get_data(), out.get_data(), n, xp_size, depends); - auto exec_q = x.get_queue(); + sycl::event keep = dpctl::utils::keep_args_alive(exec_q, {x, idx, xp, fp, out}, {ev}); - // size_t n = x.get_size(); - // const size_t m = xp.get_size(); - - // std::vector x_h(n); - // std::vector xp_h(n); - // std::vector fp_h(n); - // std::vector output_h(n); - - // sycl::event copy_1 = exec_q.memcpy(x_h.data(), x.get_data(), n * sizeof(float)); - // sycl::event copy_2 = exec_q.memcpy(xp_h.data(), xp.get_data(), m * sizeof(float)); - // sycl::event copy_3 = exec_q.memcpy(fp_h.data(), fp.get_data(), m * sizeof(float)); - // sycl::event copy_4 = exec_q.memcpy(output_h.data(), output.get_data(), n * sizeof(float)); - - // copy_1.wait(); - // copy_2.wait(); - // copy_3.wait(); - // copy_4.wait(); - - // std::cout << "x: " << std::endl; - // for (size_t i = 0; i< n; i++){ - // std::cout << x_h[i] << " "; - // } - // std::cout << "\n"; - - // std::cout << "xp: " << std::endl; - // for (size_t i = 0; i< m; i++){ - // std::cout << xp_h[i] << " "; - // } - // std::cout << "\n"; - - - // std::cout << "fp: " << std::endl; - // for (size_t i = 0; i< m; i++){ - // std::cout << fp_h[i] << " "; - // } - // std::cout << "\n"; - - - // std::cout << "out: " << std::endl; - // for (size_t i = 0; i< n; i++){ - // std::cout << output_h[i] << " "; - // } - // std::cout << "\n"; - - auto ev = - interp_func(exec_q, x.get_data(), idx.get_data(), xp.get_data(), fp.get_data(), - output.get_data(), xp.get_size(), x.get_size(), depends); - - ev.wait(); - - auto args_ev = dpctl::utils::keep_args_alive( - exec_q, {x, idx, xp, fp, output}, {ev}); - - // size_t n = x.get_size(); - - // std::vector result(n); - // sycl::event copy_ev = exec_q.memcpy(result.data(), output.get_data(), n * sizeof(float), {ev}); - - // copy_ev.wait(); - - // std::cout << "out_host: " << std::endl; - // for (size_t i = 0; i< n; i++){ - // std::cout << result[i] << " "; - // } - - // std::cout << "\n"; - - return {args_ev, ev}; + return std::make_pair(keep, ev); } -std::unique_ptr interp; -void math::interpolate::populate_interpolate(py::module_ m) +void init_interpolate_dispatch_table() { - using namespace std::placeholders; + using namespace td_ns; + using kernels::interpolate_fn_ptr_t; - interp.reset(new Interpolate()); + DispatchVectorBuilder + dtb_interpolate; + dtb_interpolate.populate_dispatch_vector(interpolate_dispatch_table); +} - auto interp_func = - [interpp = interp.get()]( - const dpctl::tensor::usm_ndarray &x, - const dpctl::tensor::usm_ndarray &idx, - const dpctl::tensor::usm_ndarray &xp, - const dpctl::tensor::usm_ndarray &fp, - const size_t xp_size, - dpctl::tensor::usm_ndarray &output, - const std::vector &depends) { - return interpp->call(x, idx, xp, fp, xp_size, output, depends); - }; +void init_interpolate(py::module_ m) +{ + dpnp::extensions::math::init_interpolate_dispatch_table(); - m.def("interpolate", interp_func, "Perform linear interpolation.", + m.def("_interpolate", &py_interpolate, "", py::arg("x"), py::arg("idx"), py::arg("xp"), py::arg("fp"), - py::arg("xp_size"), py::arg("output"), py::arg("depends") = py::list()); - - auto interpolate_dtypes = [interpp = interp.get()]() { - return interpp->dispatch_table.get_all_supported_types(); - }; - - m.def("interpolate_dtypes", interpolate_dtypes, - "Get the supported data types for interpolation."); + py::arg("out"), py::arg("sycl_queue"), py::arg("depends") = py::list()); } + +} // namespace dpnp::extensions::math diff --git a/dpnp/backend/extensions/math/interpolate.hpp b/dpnp/backend/extensions/math/interpolate.hpp index 95c34d945dad..e5df239aabdf 100644 --- a/dpnp/backend/extensions/math/interpolate.hpp +++ b/dpnp/backend/extensions/math/interpolate.hpp @@ -25,42 +25,11 @@ #pragma once -#include +#include -#include "dispatch_table.hpp" +namespace py = pybind11; -namespace dpctl_td_ns = dpctl::tensor::type_dispatch; - -namespace math -{ -namespace interpolate +namespace dpnp::extensions::math { -struct Interpolate -{ - using FnT = sycl::event (*)(sycl::queue &, - const void *, - const void *, - const void *, - const void *, - void *, - const size_t, - const size_t, - const std::vector &); - - common::DispatchTable2 dispatch_table; - - Interpolate(); - - std::tuple - call(const dpctl::tensor::usm_ndarray &x, - const dpctl::tensor::usm_ndarray &idx, - const dpctl::tensor::usm_ndarray &xp, - const dpctl::tensor::usm_ndarray &fp, - const size_t xp_size, - dpctl::tensor::usm_ndarray &output, - const std::vector &depends); -}; - -void populate_interpolate(py::module_ m); -} // namespace interpolate -} // namespace math +void init_interpolate(py::module_ m); +} // namespace dpnp::extensions::math diff --git a/dpnp/backend/extensions/math/interpolate_kernel.hpp b/dpnp/backend/extensions/math/interpolate_kernel.hpp new file mode 100644 index 000000000000..bf63de08c280 --- /dev/null +++ b/dpnp/backend/extensions/math/interpolate_kernel.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include +#include +#include +#include + +#include + +#include "utils/type_utils.hpp" + +namespace dpnp::extensions::math::kernels +{ + +using interpolate_fn_ptr_t = sycl::event (*)(sycl::queue &, + const void *, // x + const void *, // idx + const void *, // xp + const void *, // fp + void *, // out + std::size_t, // n + std::size_t, // xp_size + const std::vector &); + +template +sycl::event interpolate_impl(sycl::queue &q, + const void *vx, + const void *vidx, + const void *vxp, + const void *vfp, + void *vout, + std::size_t n, + std::size_t xp_size, + const std::vector &depends) +{ + const T *x = static_cast(vx); + const std::size_t *idx = static_cast(vidx); + const T *xp = static_cast(vxp); + const T *fp = static_cast(vfp); + T *out = static_cast(vout); + + return q.submit([&](sycl::handler &h) { + h.depends_on(depends); + h.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { + T left = fp[0]; + T right = fp[xp_size - 1]; + std::size_t x_idx = idx[i] - 1; + + if (sycl::isnan(x[i])) { + out[i] = x[i]; + } + else if (x_idx < 0) { + out[i] = left; + } + else if (x[i] == xp[xp_size - 1]) { + out[i] = right; + } + else if (x_idx >= xp_size - 1) { + out[i] = right; + } + else if (x[i] == xp[x_idx]) { + out[i] = fp[x_idx]; + } + else { + T slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); + T res = slope * (x[i] - xp[x_idx]) + fp[x_idx]; + + if (sycl::isnan(res)) { + res = slope * (x[i] - xp[x_idx + 1]) + fp[x_idx + 1]; + if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { + res = fp[x_idx]; + } + } + out[i] = res; + } + }); + }); +} + +} // namespace dpnp::extensions::math::kernels diff --git a/dpnp/backend/extensions/math/math_py.cpp b/dpnp/backend/extensions/math/math_py.cpp index 6b8b5b9b5f28..29348d9e437c 100644 --- a/dpnp/backend/extensions/math/math_py.cpp +++ b/dpnp/backend/extensions/math/math_py.cpp @@ -33,5 +33,5 @@ PYBIND11_MODULE(_math_impl, m) { - math::interpolate::populate_interpolate(m); + dpnp::extensions::math::init_interpolate(m); } diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 08b830b1d360..72da99ca9e86 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -1200,7 +1200,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): # NumPy always returns float64 or complex128, so we upcast all values # on the fly in the kernel - out_dtype = 'f8' + out_dtype = x_dtype output = dpnp.empty(x.shape, dtype=out_dtype) idx = dpnp.searchsorted(xp, x, side='right') left = fp[0] if left is None else dpnp.array(left, fp.dtype) @@ -1210,15 +1210,15 @@ def interp(x, xp, fp, left=None, right=None, period=None): queue = x.sycl_queue _manager = dpu.SequentialOrderManager[queue] - mem_ev, ht_ev = math_ext.interpolate( + mem_ev, ht_ev = math_ext._interpolate( x.get_array(), idx.get_array(), xp.get_array(), fp.get_array(), - xp.size, # left, # right, output.get_array(), + queue, depends=_manager.submitted_events, ) _manager.add_event_pair(mem_ev, ht_ev) diff --git a/dpnp/tests/third_party/cupy/test_type_routines.py b/dpnp/tests/third_party/cupy/test_type_routines.py index e35b40d90841..bf5c7af9ded0 100644 --- a/dpnp/tests/third_party/cupy/test_type_routines.py +++ b/dpnp/tests/third_party/cupy/test_type_routines.py @@ -47,6 +47,10 @@ def test_can_cast(self, xp, from_dtype, to_dtype): return ret +<<<<<<< HEAD +======= +# @pytest.mark.skip("dpnp.common_type() is not implemented yet") +>>>>>>> e8871f0d797 (Second impl with dispatch_vector[only floating]) class TestCommonType(unittest.TestCase): @testing.numpy_cupy_equal() From f7d1da94dfbaf36ada8bf0d6ec73b75406970f0c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 2 Apr 2025 04:32:50 -0700 Subject: [PATCH 03/26] Implement interpolate_complex --- dpnp/backend/extensions/math/CMakeLists.txt | 1 - dpnp/backend/extensions/math/common.cpp | 75 ----- dpnp/backend/extensions/math/common.hpp | 93 ------ .../extensions/math/dispatch_table.hpp | 292 ------------------ dpnp/backend/extensions/math/interpolate.cpp | 79 ++--- .../extensions/math/interpolate_kernel.hpp | 132 ++++++-- dpnp/dpnp_iface_statistics.py | 30 +- 7 files changed, 166 insertions(+), 536 deletions(-) delete mode 100644 dpnp/backend/extensions/math/common.cpp delete mode 100644 dpnp/backend/extensions/math/common.hpp delete mode 100644 dpnp/backend/extensions/math/dispatch_table.hpp diff --git a/dpnp/backend/extensions/math/CMakeLists.txt b/dpnp/backend/extensions/math/CMakeLists.txt index fed91cfd3f9a..eed898b12496 100644 --- a/dpnp/backend/extensions/math/CMakeLists.txt +++ b/dpnp/backend/extensions/math/CMakeLists.txt @@ -26,7 +26,6 @@ set(python_module_name _math_impl) set(_module_src - ${CMAKE_CURRENT_SOURCE_DIR}/common.cpp ${CMAKE_CURRENT_SOURCE_DIR}/interpolate.cpp ${CMAKE_CURRENT_SOURCE_DIR}/math_py.cpp ) diff --git a/dpnp/backend/extensions/math/common.cpp b/dpnp/backend/extensions/math/common.cpp deleted file mode 100644 index 93723e838b21..000000000000 --- a/dpnp/backend/extensions/math/common.cpp +++ /dev/null @@ -1,75 +0,0 @@ -//***************************************************************************** -// 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 "common.hpp" -#include "utils/type_dispatch.hpp" -#include - -namespace dpctl_td_ns = dpctl::tensor::type_dispatch; - -namespace math -{ -namespace common -{ -pybind11::dtype dtype_from_typenum(int dst_typenum) -{ - dpctl_td_ns::typenum_t dst_typenum_t = - static_cast(dst_typenum); - switch (dst_typenum_t) { - case dpctl_td_ns::typenum_t::BOOL: - return py::dtype("?"); - case dpctl_td_ns::typenum_t::INT8: - return py::dtype("i1"); - case dpctl_td_ns::typenum_t::UINT8: - return py::dtype("u1"); - case dpctl_td_ns::typenum_t::INT16: - return py::dtype("i2"); - case dpctl_td_ns::typenum_t::UINT16: - return py::dtype("u2"); - case dpctl_td_ns::typenum_t::INT32: - return py::dtype("i4"); - case dpctl_td_ns::typenum_t::UINT32: - return py::dtype("u4"); - case dpctl_td_ns::typenum_t::INT64: - return py::dtype("i8"); - case dpctl_td_ns::typenum_t::UINT64: - return py::dtype("u8"); - case dpctl_td_ns::typenum_t::HALF: - return py::dtype("f2"); - case dpctl_td_ns::typenum_t::FLOAT: - return py::dtype("f4"); - case dpctl_td_ns::typenum_t::DOUBLE: - return py::dtype("f8"); - case dpctl_td_ns::typenum_t::CFLOAT: - return py::dtype("c8"); - case dpctl_td_ns::typenum_t::CDOUBLE: - return py::dtype("c16"); - default: - throw py::value_error("Unrecognized dst_typeid"); - } -} - -} // namespace common -} // namespace math diff --git a/dpnp/backend/extensions/math/common.hpp b/dpnp/backend/extensions/math/common.hpp deleted file mode 100644 index 1d436440c392..000000000000 --- a/dpnp/backend/extensions/math/common.hpp +++ /dev/null @@ -1,93 +0,0 @@ -//***************************************************************************** -// 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 -#include -#include - -// clang-format off -// math_utils.hpp doesn't include sycl header but uses sycl types -// so sycl.hpp must be included before math_utils.hpp -#include -#include "utils/math_utils.hpp" -// clang-format on - -namespace math -{ -namespace common -{ -template -struct Less -{ - bool operator()(const T &lhs, const T &rhs) const - { - return std::less{}(lhs, rhs); - } -}; - -template -struct Less> -{ - bool operator()(const std::complex &lhs, - const std::complex &rhs) const - { - return dpctl::tensor::math_utils::less_complex(lhs, rhs); - } -}; - -template -struct IsNan -{ - static bool isnan(const T &v) - { - if constexpr (std::is_floating_point_v || - std::is_same_v) { - return sycl::isnan(v); - } - - return false; - } -}; - -template -struct IsNan> -{ - static bool isnan(const std::complex &v) - { - T real1 = std::real(v); - T imag1 = std::imag(v); - return sycl::isnan(real1) || sycl::isnan(imag1); - } -}; - - -// This function is a copy from dpctl because it is not available in the public -// headers of dpctl. -pybind11::dtype dtype_from_typenum(int dst_typenum); - -} // namespace common -} // namespace math diff --git a/dpnp/backend/extensions/math/dispatch_table.hpp b/dpnp/backend/extensions/math/dispatch_table.hpp deleted file mode 100644 index 4cfd3d2a09a4..000000000000 --- a/dpnp/backend/extensions/math/dispatch_table.hpp +++ /dev/null @@ -1,292 +0,0 @@ -//***************************************************************************** -// 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 -#include - -#include "utils/type_dispatch.hpp" -#include -#include -#include -#include - -#include "common.hpp" - -namespace dpctl_td_ns = dpctl::tensor::type_dispatch; -namespace py = pybind11; - -namespace math -{ -namespace common -{ - -template -struct one_of -{ - static_assert(std::is_same_v>, - "one_of: second parameter cannot be empty std::tuple"); - static_assert(false, "one_of: second parameter must be std::tuple"); -}; - -template -struct one_of> -{ - static constexpr bool value = - std::is_same_v || one_of>::value; -}; - -template -struct one_of> -{ - static constexpr bool value = std::is_same_v; -}; - -template -constexpr bool one_of_v = one_of::value; - -template -using Table = FnT[dpctl_td_ns::num_types]; -template -using Table2 = Table[dpctl_td_ns::num_types]; - -using TypeId = int32_t; -using TypesPair = std::pair; - -struct int_pair_hash -{ - inline size_t operator()(const TypesPair &p) const - { - std::hash hasher; - return hasher(size_t(p.first) << (8 * sizeof(TypeId)) | - size_t(p.second)); - } -}; - -using SupportedTypesList = std::vector; -using SupportedTypesList2 = std::vector; -using SupportedTypesSet = std::unordered_set; -using SupportedTypesSet2 = std::unordered_set; - -using DType = py::dtype; -using DTypePair = std::pair; - -using SupportedDTypeList = std::vector; -using SupportedDTypeList2 = std::vector; - -template - typename Func> -struct TableBuilder2 -{ - template - struct impl - { - static constexpr bool is_defined = - one_of_v, SupportedTypes>; - - _FnT get() - { - if constexpr (is_defined) { - return Func::impl; - } - else { - return nullptr; - } - } - }; - - using type = - dpctl_td_ns::DispatchTableBuilder; -}; - -template -class DispatchTable2 -{ -public: - DispatchTable2(std::string first_name, std::string second_name) - : first_name(first_name), second_name(second_name) - { - } - - template - typename Func> - void populate_dispatch_table() - { - using TBulder = typename TableBuilder2::type; - TBulder builder; - - builder.populate_dispatch_table(table); - populate_supported_types(); - } - - FnT get_unsafe(int first_typenum, int second_typenum) const - { - auto array_types = dpctl_td_ns::usm_ndarray_types(); - const int first_type_id = - array_types.typenum_to_lookup_id(first_typenum); - const int second_type_id = - array_types.typenum_to_lookup_id(second_typenum); - - return table[first_type_id][second_type_id]; - } - - FnT get(int first_typenum, int second_typenum) const - { - auto fn = get_unsafe(first_typenum, second_typenum); - - if (fn == nullptr) { - auto array_types = dpctl_td_ns::usm_ndarray_types(); - const int first_type_id = - array_types.typenum_to_lookup_id(first_typenum); - const int second_type_id = - array_types.typenum_to_lookup_id(second_typenum); - - py::dtype first_dtype = dtype_from_typenum(first_type_id); - auto first_type_pos = - std::find(supported_first_type.begin(), - supported_first_type.end(), first_dtype); - if (first_type_pos == supported_first_type.end()) { - py::str types = py::str(py::cast(supported_first_type)); - py::str dtype = py::str(first_dtype); - - py::str err_msg = - py::str("'" + first_name + "' has unsupported type '") + - dtype + - py::str("'." - " Supported types are: ") + - types; - - throw py::value_error(static_cast(err_msg)); - } - - py::dtype second_dtype = dtype_from_typenum(second_type_id); - auto second_type_pos = - std::find(supported_second_type.begin(), - supported_second_type.end(), second_dtype); - if (second_type_pos == supported_second_type.end()) { - py::str types = py::str(py::cast(supported_second_type)); - py::str dtype = py::str(second_dtype); - - py::str err_msg = - py::str("'" + second_name + "' has unsupported type '") + - dtype + - py::str("'." - " Supported types are: ") + - types; - - throw py::value_error(static_cast(err_msg)); - } - - py::str first_dtype_str = py::str(first_dtype); - py::str second_dtype_str = py::str(second_dtype); - py::str types = py::str(py::cast(all_supported_types)); - - py::str err_msg = - py::str("'" + first_name + "' and '" + second_name + - "' has unsupported types combination: ('") + - first_dtype_str + py::str("', '") + second_dtype_str + - py::str("')." - " Supported types combinations are: ") + - types; - - throw py::value_error(static_cast(err_msg)); - } - - return fn; - } - - const SupportedDTypeList &get_supported_first_type() const - { - return supported_first_type; - } - - const SupportedDTypeList &get_supported_second_type() const - { - return supported_second_type; - } - - const SupportedDTypeList2 &get_all_supported_types() const - { - return all_supported_types; - } - -private: - void populate_supported_types() - { - SupportedTypesSet first_supported_types_set; - SupportedTypesSet second_supported_types_set; - SupportedTypesSet2 all_supported_types_set; - - for (int i = 0; i < dpctl_td_ns::num_types; ++i) { - for (int j = 0; j < dpctl_td_ns::num_types; ++j) { - if (table[i][j] != nullptr) { - all_supported_types_set.emplace(i, j); - first_supported_types_set.emplace(i); - second_supported_types_set.emplace(j); - } - } - } - - auto to_supported_dtype_list = [](const auto &supported_set, - auto &supported_list) { - SupportedTypesList lst(supported_set.begin(), supported_set.end()); - std::sort(lst.begin(), lst.end()); - supported_list.resize(supported_set.size()); - std::transform(lst.begin(), lst.end(), supported_list.begin(), - [](TypeId i) { return dtype_from_typenum(i); }); - }; - - to_supported_dtype_list(first_supported_types_set, - supported_first_type); - to_supported_dtype_list(second_supported_types_set, - supported_second_type); - - SupportedTypesList2 lst(all_supported_types_set.begin(), - all_supported_types_set.end()); - std::sort(lst.begin(), lst.end()); - all_supported_types.resize(all_supported_types_set.size()); - std::transform(lst.begin(), lst.end(), all_supported_types.begin(), - [](TypesPair p) { - return DTypePair(dtype_from_typenum(p.first), - dtype_from_typenum(p.second)); - }); - } - - std::string first_name; - std::string second_name; - - SupportedDTypeList supported_first_type; - SupportedDTypeList supported_second_type; - SupportedDTypeList2 all_supported_types; - - Table2 table; -}; - -} // namespace common -} // namespace math diff --git a/dpnp/backend/extensions/math/interpolate.cpp b/dpnp/backend/extensions/math/interpolate.cpp index a620299bfffa..c1549826e77c 100644 --- a/dpnp/backend/extensions/math/interpolate.cpp +++ b/dpnp/backend/extensions/math/interpolate.cpp @@ -23,10 +23,6 @@ // THE POSSIBILITY OF SUCH DAMAGE. //***************************************************************************** -#include -#include -#include -#include #include #include @@ -45,37 +41,26 @@ namespace dpnp::extensions::math namespace py = pybind11; namespace td_ns = dpctl::tensor::type_dispatch; -static kernels::interpolate_fn_ptr_t interpolate_dispatch_table[td_ns::num_types]; - -template -struct InterpolateFactory -{ - fnT get() - { - if constexpr (std::is_floating_point_v) { - return kernels::interpolate_impl; - } - else { - return nullptr; - } - } -}; - +static kernels::interpolate_fn_ptr_t + interpolate_dispatch_table[td_ns::num_types][td_ns::num_types]; std::pair -py_interpolate(const dpctl::tensor::usm_ndarray &x, - const dpctl::tensor::usm_ndarray &idx, - const dpctl::tensor::usm_ndarray &xp, - const dpctl::tensor::usm_ndarray &fp, - dpctl::tensor::usm_ndarray &out, - sycl::queue &exec_q, - const std::vector &depends) + py_interpolate(const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + dpctl::tensor::usm_ndarray &out, + sycl::queue &exec_q, + const std::vector &depends) { - int typenum = x.get_typenum(); + int xp_typenum = xp.get_typenum(); + int fp_typenum = fp.get_typenum(); + auto array_types = td_ns::usm_ndarray_types(); - int type_id = array_types.typenum_to_lookup_id(typenum); + int xp_type_id = array_types.typenum_to_lookup_id(xp_typenum); + int fp_type_id = array_types.typenum_to_lookup_id(fp_typenum); - auto fn = interpolate_dispatch_table[type_id]; + auto fn = interpolate_dispatch_table[xp_type_id][fp_type_id]; if (!fn) { throw py::type_error("Unsupported dtype."); } @@ -86,29 +71,51 @@ py_interpolate(const dpctl::tensor::usm_ndarray &x, sycl::event ev = fn(exec_q, x.get_data(), idx.get_data(), xp.get_data(), fp.get_data(), out.get_data(), n, xp_size, depends); - sycl::event keep = dpctl::utils::keep_args_alive(exec_q, {x, idx, xp, fp, out}, {ev}); + sycl::event keep = + dpctl::utils::keep_args_alive(exec_q, {x, idx, xp, fp, out}, {ev}); return std::make_pair(keep, ev); } +template +struct InterpolateFactory +{ + fnT get() + { + if constexpr (std::is_floating_point_v && + std::is_floating_point_v) + { + return kernels::interpolate_impl; + } + else if constexpr (std::is_floating_point_v && + (std::is_same_v> || + std::is_same_v>)) + { + return kernels::interpolate_complex_impl; + } + else { + return nullptr; + } + } +}; void init_interpolate_dispatch_table() { using namespace td_ns; using kernels::interpolate_fn_ptr_t; - DispatchVectorBuilder + DispatchTableBuilder dtb_interpolate; - dtb_interpolate.populate_dispatch_vector(interpolate_dispatch_table); + dtb_interpolate.populate_dispatch_table(interpolate_dispatch_table); } void init_interpolate(py::module_ m) { dpnp::extensions::math::init_interpolate_dispatch_table(); - m.def("_interpolate", &py_interpolate, "", - py::arg("x"), py::arg("idx"), py::arg("xp"), py::arg("fp"), - py::arg("out"), py::arg("sycl_queue"), py::arg("depends") = py::list()); + m.def("_interpolate", &py_interpolate, "", py::arg("x"), py::arg("idx"), + py::arg("xp"), py::arg("fp"), py::arg("out"), py::arg("sycl_queue"), + py::arg("depends") = py::list()); } } // namespace dpnp::extensions::math diff --git a/dpnp/backend/extensions/math/interpolate_kernel.hpp b/dpnp/backend/extensions/math/interpolate_kernel.hpp index bf63de08c280..82ce0c19a33d 100644 --- a/dpnp/backend/extensions/math/interpolate_kernel.hpp +++ b/dpnp/backend/extensions/math/interpolate_kernel.hpp @@ -2,10 +2,6 @@ #include #include -#include -#include - -#include #include "utils/type_utils.hpp" @@ -13,16 +9,16 @@ namespace dpnp::extensions::math::kernels { using interpolate_fn_ptr_t = sycl::event (*)(sycl::queue &, - const void *, // x - const void *, // idx - const void *, // xp - const void *, // fp - void *, // out - std::size_t, // n - std::size_t, // xp_size + const void *, // x + const void *, // idx + const void *, // xp + const void *, // fp + void *, // out + std::size_t, // n + std::size_t, // xp_size const std::vector &); -template +template sycl::event interpolate_impl(sycl::queue &q, const void *vx, const void *vidx, @@ -33,40 +29,43 @@ sycl::event interpolate_impl(sycl::queue &q, std::size_t xp_size, const std::vector &depends) { - const T *x = static_cast(vx); + const TCoord *x = static_cast(vx); const std::size_t *idx = static_cast(vidx); - const T *xp = static_cast(vxp); - const T *fp = static_cast(vfp); - T *out = static_cast(vout); + const TCoord *xp = static_cast(vxp); + const TValue *fp = static_cast(vfp); + TValue *out = static_cast(vout); return q.submit([&](sycl::handler &h) { h.depends_on(depends); h.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { - T left = fp[0]; - T right = fp[xp_size - 1]; + TValue left = fp[0]; + TValue right = fp[xp_size - 1]; + + TCoord x_val = x[i]; std::size_t x_idx = idx[i] - 1; - if (sycl::isnan(x[i])) { - out[i] = x[i]; + if (sycl::isnan(x_val)) { + out[i] = x_val; } else if (x_idx < 0) { out[i] = left; } - else if (x[i] == xp[xp_size - 1]) { + else if (x_val == xp[xp_size - 1]) { out[i] = right; } else if (x_idx >= xp_size - 1) { out[i] = right; } - else if (x[i] == xp[x_idx]) { + else if (x_val == xp[x_idx]) { out[i] = fp[x_idx]; } else { - T slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); - T res = slope * (x[i] - xp[x_idx]) + fp[x_idx]; + TValue slope = + (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); + TValue res = slope * (x_val - xp[x_idx]) + fp[x_idx]; if (sycl::isnan(res)) { - res = slope * (x[i] - xp[x_idx + 1]) + fp[x_idx + 1]; + res = slope * (x_val - xp[x_idx + 1]) + fp[x_idx + 1]; if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { res = fp[x_idx]; } @@ -77,4 +76,87 @@ sycl::event interpolate_impl(sycl::queue &q, }); } +template +sycl::event interpolate_complex_impl(sycl::queue &q, + const void *vx, + const void *vidx, + const void *vxp, + const void *vfp, + void *vout, + std::size_t n, + std::size_t xp_size, + const std::vector &depends) +{ + const TCoord *x = static_cast(vx); + const std::size_t *idx = static_cast(vidx); + const TCoord *xp = static_cast(vxp); + const TValue *fp = static_cast(vfp); + TValue *out = static_cast(vout); + + using realT = typename TValue::value_type; + + return q.submit([&](sycl::handler &h) { + h.depends_on(depends); + h.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { + realT left_r = fp[0].real(); + realT right_r = fp[xp_size - 1].real(); + realT left_i = fp[0].imag(); + realT right_i = fp[xp_size - 1].imag(); + + TCoord x_val = x[i]; + std::size_t x_idx = idx[i] - 1; + + realT res_r = 0.0; + realT res_i = 0.0; + + if (sycl::isnan(x_val)) { + res_r = x_val; + res_i = 0.0; + } + else if (x_idx < 0) { + res_r = left_r; + res_i = left_i; + } + else if (x_val == xp[xp_size - 1]) { + res_r = right_r; + res_i = right_i; + } + else if (x_idx >= xp_size - 1) { + res_r = right_r; + res_i = right_i; + } + else if (x_val == xp[x_idx]) { + res_r = fp[x_idx].real(); + res_i = fp[x_idx].imag(); + } + else { + realT dx = xp[x_idx + 1] - xp[x_idx]; + + realT slope_r = (fp[x_idx + 1].real() - fp[x_idx].real()) / dx; + res_r = slope_r * (x_val - xp[x_idx]) + fp[x_idx].real(); + if (sycl::isnan(res_r)) { + res_r = slope_r * (x_val - xp[x_idx + 1]) + + fp[x_idx + 1].real(); + if (sycl::isnan(res_r) && + fp[x_idx].real() == fp[x_idx + 1].real()) { + res_r = fp[x_idx].real(); + } + } + + realT slope_i = (fp[x_idx + 1].imag() - fp[x_idx].imag()) / dx; + res_i = slope_i * (x_val - xp[x_idx]) + fp[x_idx].imag(); + if (sycl::isnan(res_i)) { + res_i = slope_i * (x_val - xp[x_idx + 1]) + + fp[x_idx + 1].imag(); + if (sycl::isnan(res_i) && + fp[x_idx].imag() == fp[x_idx + 1].imag()) { + res_i = fp[x_idx].imag(); + } + } + } + out[i] = TValue(res_r, res_i); + }); + }); +} + } // namespace dpnp::extensions::math::kernels diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 72da99ca9e86..802243139de5 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -46,6 +46,7 @@ from dpctl.tensor._numpy_helper import normalize_axis_index import dpnp +import dpnp.backend.extensions.math._math_impl as math_ext # pylint: disable=no-name-in-module import dpnp.backend.extensions.statistics._statistics_impl as statistics_ext @@ -54,8 +55,6 @@ to_supported_dtypes, ) -import dpnp.backend.extensions.math._math_impl as math_ext - from .dpnp_utils import get_usm_allocations from .dpnp_utils.dpnp_utils_reduction import dpnp_wrap_reduction_call from .dpnp_utils.dpnp_utils_statistics import dpnp_cov, dpnp_median @@ -1161,19 +1160,22 @@ def interp(x, xp, fp, left=None, right=None, period=None): dpnp.check_supported_arrays_type(x, xp, fp) if xp.ndim != 1 or fp.ndim != 1: - raise ValueError('xp and fp must be 1D arrays') + raise ValueError("xp and fp must be 1D arrays") if xp.size != fp.size: - raise ValueError('fp and xp are not of the same length') + raise ValueError("fp and xp are not of the same length") if xp.size == 0: - raise ValueError('array of sample points is empty') + raise ValueError("array of sample points is empty") if not x.flags.c_contiguous: - raise NotImplementedError('Non-C-contiguous x is currently not ' - 'supported') + raise NotImplementedError( + "Non-C-contiguous x is currently not supported" + ) x_dtype = dpnp.common_type(x, xp) if not dpnp.can_cast(x_dtype, dpnp.default_float_type()): - raise TypeError('Cannot cast array data from' - ' {} to {} according to the rule \'safe\'' - .format(x_dtype, dpnp.default_float_type())) + raise TypeError( + "Cannot cast array data from" + f" {x_dtype} to {dpnp.default_float_type()} " + "according to the rule 'safe'" + ) if period is not None: # The handling of "period" below is modified from NumPy's @@ -1193,20 +1195,20 @@ def interp(x, xp, fp, left=None, right=None, period=None): asort_xp = dpnp.argsort(xp) xp = xp[asort_xp] fp = fp[asort_xp] - xp = dpnp.concatenate((xp[-1:]-period, xp, xp[0:1]+period)) + xp = dpnp.concatenate((xp[-1:] - period, xp, xp[0:1] + period)) fp = dpnp.concatenate((fp[-1:], fp, fp[0:1])) assert xp.flags.c_contiguous assert fp.flags.c_contiguous # NumPy always returns float64 or complex128, so we upcast all values # on the fly in the kernel - out_dtype = x_dtype + out_dtype = fp.dtype output = dpnp.empty(x.shape, dtype=out_dtype) - idx = dpnp.searchsorted(xp, x, side='right') + idx = dpnp.searchsorted(xp, x, side="right") left = fp[0] if left is None else dpnp.array(left, fp.dtype) right = fp[-1] if right is None else dpnp.array(right, fp.dtype) - idx = dpnp.array(idx, dtype='uint64') + idx = dpnp.array(idx, dtype="uint64") queue = x.sycl_queue _manager = dpu.SequentialOrderManager[queue] From e1b86982b19f1422c1d030eaaa01fcb35802935c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 2 Apr 2025 04:33:20 -0700 Subject: [PATCH 04/26] Move interpolate backend to ufunc --- dpnp/CMakeLists.txt | 1 - dpnp/backend/extensions/math/CMakeLists.txt | 86 ------------------- dpnp/backend/extensions/math/math_py.cpp | 37 -------- dpnp/backend/extensions/ufunc/CMakeLists.txt | 1 + .../ufunc/elementwise_functions/common.cpp | 2 + .../elementwise_functions}/interpolate.cpp | 39 ++++++--- .../elementwise_functions}/interpolate.hpp | 4 +- .../elementwise_functions/interpolate.hpp} | 14 +-- dpnp/dpnp_iface_statistics.py | 4 +- 9 files changed, 36 insertions(+), 152 deletions(-) delete mode 100644 dpnp/backend/extensions/math/CMakeLists.txt delete mode 100644 dpnp/backend/extensions/math/math_py.cpp rename dpnp/backend/extensions/{math => ufunc/elementwise_functions}/interpolate.cpp (76%) rename dpnp/backend/extensions/{math => ufunc/elementwise_functions}/interpolate.hpp (95%) rename dpnp/backend/{extensions/math/interpolate_kernel.hpp => kernels/elementwise_functions/interpolate.hpp} (88%) diff --git a/dpnp/CMakeLists.txt b/dpnp/CMakeLists.txt index 6c59141bd1d9..6be90d849dc4 100644 --- a/dpnp/CMakeLists.txt +++ b/dpnp/CMakeLists.txt @@ -60,7 +60,6 @@ add_subdirectory(backend/extensions/blas) add_subdirectory(backend/extensions/fft) add_subdirectory(backend/extensions/indexing) add_subdirectory(backend/extensions/lapack) -add_subdirectory(backend/extensions/math) add_subdirectory(backend/extensions/statistics) add_subdirectory(backend/extensions/ufunc) add_subdirectory(backend/extensions/vm) diff --git a/dpnp/backend/extensions/math/CMakeLists.txt b/dpnp/backend/extensions/math/CMakeLists.txt deleted file mode 100644 index eed898b12496..000000000000 --- a/dpnp/backend/extensions/math/CMakeLists.txt +++ /dev/null @@ -1,86 +0,0 @@ -# ***************************************************************************** -# Copyright (c) 2016-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. -# ***************************************************************************** - - -set(python_module_name _math_impl) -set(_module_src - ${CMAKE_CURRENT_SOURCE_DIR}/interpolate.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/math_py.cpp -) - -pybind11_add_module(${python_module_name} MODULE ${_module_src}) -add_sycl_to_target(TARGET ${python_module_name} SOURCES ${_module_src}) - -if(_dpnp_sycl_targets) - # make fat binary - target_compile_options( - ${python_module_name} - PRIVATE - -fsycl-targets=${_dpnp_sycl_targets} - ) - target_link_options( - ${python_module_name} - PRIVATE - -fsycl-targets=${_dpnp_sycl_targets} - ) -endif() - -if (WIN32) - if (${CMAKE_VERSION} VERSION_LESS "3.27") - # this is a work-around for target_link_options inserting option after -link option, cause - # linker to ignore it. - set(CMAKE_CXX_LINK_FLAGS "${CMAKE_CXX_LINK_FLAGS} -fsycl-device-code-split=per_kernel") - endif() -endif() - -set_target_properties(${python_module_name} PROPERTIES CMAKE_POSITION_INDEPENDENT_CODE ON) - -target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../include) -target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../src) - -target_include_directories(${python_module_name} PUBLIC ${Dpctl_INCLUDE_DIR}) -target_include_directories(${python_module_name} PUBLIC ${Dpctl_TENSOR_INCLUDE_DIR}) - -if (WIN32) - target_compile_options(${python_module_name} PRIVATE - /clang:-fno-approx-func - /clang:-fno-finite-math-only - ) -else() - target_compile_options(${python_module_name} PRIVATE - -fno-approx-func - -fno-finite-math-only - ) -endif() - -target_link_options(${python_module_name} PUBLIC -fsycl-device-code-split=per_kernel) - -if (DPNP_GENERATE_COVERAGE) - target_link_options(${python_module_name} PRIVATE -fprofile-instr-generate -fcoverage-mapping) -endif() - -install(TARGETS ${python_module_name} - DESTINATION "dpnp/backend/extensions/math" -) diff --git a/dpnp/backend/extensions/math/math_py.cpp b/dpnp/backend/extensions/math/math_py.cpp deleted file mode 100644 index 29348d9e437c..000000000000 --- a/dpnp/backend/extensions/math/math_py.cpp +++ /dev/null @@ -1,37 +0,0 @@ -//***************************************************************************** -// 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. -//***************************************************************************** -// -// This file defines functions of dpnp.backend._math_impl extensions -// -//***************************************************************************** - -#include - -#include "interpolate.hpp" - -PYBIND11_MODULE(_math_impl, m) -{ - dpnp::extensions::math::init_interpolate(m); -} diff --git a/dpnp/backend/extensions/ufunc/CMakeLists.txt b/dpnp/backend/extensions/ufunc/CMakeLists.txt index d363910f74df..99163539191c 100644 --- a/dpnp/backend/extensions/ufunc/CMakeLists.txt +++ b/dpnp/backend/extensions/ufunc/CMakeLists.txt @@ -36,6 +36,7 @@ set(_elementwise_sources ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/gcd.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/heaviside.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/i0.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/interpolate.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/lcm.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/ldexp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/logaddexp2.cpp diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp index 8ff89a1b03b6..c6dd3e038eb1 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp @@ -36,6 +36,7 @@ #include "gcd.hpp" #include "heaviside.hpp" #include "i0.hpp" +#include "interpolate.hpp" #include "lcm.hpp" #include "ldexp.hpp" #include "logaddexp2.hpp" @@ -64,6 +65,7 @@ void init_elementwise_functions(py::module_ m) init_gcd(m); init_heaviside(m); init_i0(m); + init_interpolate(m); init_lcm(m); init_ldexp(m); init_logaddexp2(m); diff --git a/dpnp/backend/extensions/math/interpolate.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp similarity index 76% rename from dpnp/backend/extensions/math/interpolate.cpp rename to dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp index c1549826e77c..4de077b87997 100644 --- a/dpnp/backend/extensions/math/interpolate.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp @@ -32,17 +32,29 @@ #include "dpctl4pybind11.hpp" #include "utils/type_dispatch.hpp" -#include "interpolate.hpp" -#include "interpolate_kernel.hpp" - -namespace dpnp::extensions::math -{ +#include "kernels/elementwise_functions/interpolate.hpp" namespace py = pybind11; namespace td_ns = dpctl::tensor::type_dispatch; -static kernels::interpolate_fn_ptr_t - interpolate_dispatch_table[td_ns::num_types][td_ns::num_types]; +namespace dpnp::extensions::ufunc +{ + +namespace impl +{ + +typedef sycl::event (*interpolate_fn_ptr_t)(sycl::queue &, + const void *, // x + const void *, // idx + const void *, // xp + const void *, // fp + void *, // out + std::size_t, // n + std::size_t, // xp_size + const std::vector &); + +interpolate_fn_ptr_t interpolate_dispatch_table[td_ns::num_types] + [td_ns::num_types]; std::pair py_interpolate(const dpctl::tensor::usm_ndarray &x, @@ -85,13 +97,14 @@ struct InterpolateFactory if constexpr (std::is_floating_point_v && std::is_floating_point_v) { - return kernels::interpolate_impl; + return dpnp::kernels::interpolate::interpolate_impl; } else if constexpr (std::is_floating_point_v && (std::is_same_v> || std::is_same_v>)) { - return kernels::interpolate_complex_impl; + return dpnp::kernels::interpolate::interpolate_complex_impl; } else { return nullptr; @@ -102,20 +115,22 @@ struct InterpolateFactory void init_interpolate_dispatch_table() { using namespace td_ns; - using kernels::interpolate_fn_ptr_t; DispatchTableBuilder dtb_interpolate; dtb_interpolate.populate_dispatch_table(interpolate_dispatch_table); } +} // namespace impl + void init_interpolate(py::module_ m) { - dpnp::extensions::math::init_interpolate_dispatch_table(); + impl::init_interpolate_dispatch_table(); + using impl::py_interpolate; m.def("_interpolate", &py_interpolate, "", py::arg("x"), py::arg("idx"), py::arg("xp"), py::arg("fp"), py::arg("out"), py::arg("sycl_queue"), py::arg("depends") = py::list()); } -} // namespace dpnp::extensions::math +} // namespace dpnp::extensions::ufunc diff --git a/dpnp/backend/extensions/math/interpolate.hpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.hpp similarity index 95% rename from dpnp/backend/extensions/math/interpolate.hpp rename to dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.hpp index e5df239aabdf..4ae1cb2c8958 100644 --- a/dpnp/backend/extensions/math/interpolate.hpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.hpp @@ -29,7 +29,7 @@ namespace py = pybind11; -namespace dpnp::extensions::math +namespace dpnp::extensions::ufunc { void init_interpolate(py::module_ m); -} // namespace dpnp::extensions::math +} // namespace dpnp::extensions::ufunc diff --git a/dpnp/backend/extensions/math/interpolate_kernel.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp similarity index 88% rename from dpnp/backend/extensions/math/interpolate_kernel.hpp rename to dpnp/backend/kernels/elementwise_functions/interpolate.hpp index 82ce0c19a33d..f562c355fc37 100644 --- a/dpnp/backend/extensions/math/interpolate_kernel.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -5,19 +5,9 @@ #include "utils/type_utils.hpp" -namespace dpnp::extensions::math::kernels +namespace dpnp::kernels::interpolate { -using interpolate_fn_ptr_t = sycl::event (*)(sycl::queue &, - const void *, // x - const void *, // idx - const void *, // xp - const void *, // fp - void *, // out - std::size_t, // n - std::size_t, // xp_size - const std::vector &); - template sycl::event interpolate_impl(sycl::queue &q, const void *vx, @@ -159,4 +149,4 @@ sycl::event interpolate_complex_impl(sycl::queue &q, }); } -} // namespace dpnp::extensions::math::kernels +} // namespace dpnp::kernels::interpolate diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 802243139de5..0f52a31b730f 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -46,10 +46,10 @@ from dpctl.tensor._numpy_helper import normalize_axis_index import dpnp -import dpnp.backend.extensions.math._math_impl as math_ext # pylint: disable=no-name-in-module import dpnp.backend.extensions.statistics._statistics_impl as statistics_ext +import dpnp.backend.extensions.ufunc._ufunc_impl as ufi from dpnp.dpnp_utils.dpnp_utils_common import ( result_type_for_device, to_supported_dtypes, @@ -1212,7 +1212,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): queue = x.sycl_queue _manager = dpu.SequentialOrderManager[queue] - mem_ev, ht_ev = math_ext._interpolate( + mem_ev, ht_ev = ufi._interpolate( x.get_array(), idx.get_array(), xp.get_array(), From 00374553f20fc445db8841483445ba44395fc4cb Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 2 Apr 2025 04:35:20 -0700 Subject: [PATCH 05/26] Move def interp()to dpnp_iface_mathematical --- dpnp/dpnp_iface_mathematical.py | 178 +++++++++++++++++++++++++++++++- dpnp/dpnp_iface_statistics.py | 177 ------------------------------- 2 files changed, 177 insertions(+), 178 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 9f9d5dcca082..551e28077f55 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -111,8 +111,9 @@ "gcd", "gradient", "heaviside", - "imag", "i0", + "imag", + "interp", "lcm", "ldexp", "maximum", @@ -2689,6 +2690,181 @@ def gradient(f, *varargs, axis=None, edge_order=1): ) +def interp(x, xp, fp, left=None, right=None, period=None): + """ + One-dimensional linear interpolation for monotonically increasing sample points. + + Returns the one-dimensional piecewise linear interpolant to a function + with given discrete data points (`xp`, `fp`), evaluated at `x`. + + Parameters + ---------- + x : array_like + The x-coordinates at which to evaluate the interpolated values. + + xp : 1-D sequence of floats + The x-coordinates of the data points, must be increasing if argument + `period` is not specified. Otherwise, `xp` is internally sorted after + normalizing the periodic boundaries with ``xp = xp % period``. + + fp : 1-D sequence of float or complex + The y-coordinates of the data points, same length as `xp`. + + left : optional float or complex corresponding to fp + Value to return for `x < xp[0]`, default is `fp[0]`. + + right : optional float or complex corresponding to fp + Value to return for `x > xp[-1]`, default is `fp[-1]`. + + period : None or float, optional + A period for the x-coordinates. This parameter allows the proper + interpolation of angular x-coordinates. Parameters `left` and `right` + are ignored if `period` is specified. + + Returns + ------- + y : float or complex (corresponding to fp) or ndarray + The interpolated values, same shape as `x`. + + Raises + ------ + ValueError + If `xp` and `fp` have different length + If `xp` or `fp` are not 1-D sequences + If `period == 0` + + See Also + -------- + scipy.interpolate + + Warnings + -------- + The x-coordinate sequence is expected to be increasing, but this is not + explicitly enforced. However, if the sequence `xp` is non-increasing, + interpolation results are meaningless. + + Note that, since NaN is unsortable, `xp` also cannot contain NaNs. + + A simple check for `xp` being strictly increasing is:: + + np.all(np.diff(xp) > 0) + + Examples + -------- + >>> import dpnp as np + >>> xp = np.array([1, 2, 3]) + >>> fp = np.array([3 ,2 ,0]) + >>> x = np.array([2.5]) + >>> np.interp(2.5, xp, fp) + 1.0 + >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp) + array([3. , 3. , 2.5 , 0.56, 0. ]) + >>> UNDEF = -99.0 + >>> np.interp(3.14, xp, fp, right=UNDEF) + -99.0 + + Plot an interpolant to the sine function: + + >>> x = np.linspace(0, 2*np.pi, 10) + >>> y = np.sin(x) + >>> xvals = np.linspace(0, 2*np.pi, 50) + >>> yinterp = np.interp(xvals, x, y) + >>> import matplotlib.pyplot as plt + >>> plt.plot(x, y, 'o') + [] + >>> plt.plot(xvals, yinterp, '-x') + [] + >>> plt.show() + + Interpolation with periodic x-coordinates: + + >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] + >>> xp = [190, -190, 350, -350] + >>> fp = [5, 10, 3, 4] + >>> np.interp(x, xp, fp, period=360) + array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75]) + + Complex interpolation: + + >>> x = [1.5, 4.0] + >>> xp = [2,3,5] + >>> fp = [1.0j, 0, 2+3j] + >>> np.interp(x, xp, fp) + array([0.+1.j , 1.+1.5j]) + + """ + + dpnp.check_supported_arrays_type(x, xp, fp) + + if xp.ndim != 1 or fp.ndim != 1: + raise ValueError("xp and fp must be 1D arrays") + if xp.size != fp.size: + raise ValueError("fp and xp are not of the same length") + if xp.size == 0: + raise ValueError("array of sample points is empty") + if not x.flags.c_contiguous: + raise NotImplementedError( + "Non-C-contiguous x is currently not supported" + ) + x_dtype = dpnp.common_type(x, xp) + if not dpnp.can_cast(x_dtype, dpnp.default_float_type()): + raise TypeError( + "Cannot cast array data from" + f" {x_dtype} to {dpnp.default_float_type()} " + "according to the rule 'safe'" + ) + + if period is not None: + # The handling of "period" below is modified from NumPy's + + if period == 0: + raise ValueError("period must be a non-zero value") + period = dpnp.abs(period) + left = None + right = None + + x = x.astype(dpnp.default_float_type()) + xp = xp.astype(dpnp.default_float_type()) + + # normalizing periodic boundaries + x %= period + xp %= period + asort_xp = dpnp.argsort(xp) + xp = xp[asort_xp] + fp = fp[asort_xp] + xp = dpnp.concatenate((xp[-1:] - period, xp, xp[0:1] + period)) + fp = dpnp.concatenate((fp[-1:], fp, fp[0:1])) + assert xp.flags.c_contiguous + assert fp.flags.c_contiguous + + # NumPy always returns float64 or complex128, so we upcast all values + # on the fly in the kernel + out_dtype = fp.dtype + output = dpnp.empty(x.shape, dtype=out_dtype) + idx = dpnp.searchsorted(xp, x, side="right") + left = fp[0] if left is None else dpnp.array(left, fp.dtype) + right = fp[-1] if right is None else dpnp.array(right, fp.dtype) + + idx = dpnp.array(idx, dtype="uint64") + + queue = x.sycl_queue + _manager = dpu.SequentialOrderManager[queue] + mem_ev, ht_ev = ufi._interpolate( + x.get_array(), + idx.get_array(), + xp.get_array(), + fp.get_array(), + # left, + # right, + output.get_array(), + queue, + depends=_manager.submitted_events, + ) + _manager.add_event_pair(mem_ev, ht_ev) + + return output + + _LCM_DOCSTRING = """ Returns the lowest common multiple of ``|x1|`` and ``|x2|``. diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 0f52a31b730f..3958127789ab 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -49,7 +49,6 @@ # pylint: disable=no-name-in-module import dpnp.backend.extensions.statistics._statistics_impl as statistics_ext -import dpnp.backend.extensions.ufunc._ufunc_impl as ufi from dpnp.dpnp_utils.dpnp_utils_common import ( result_type_for_device, to_supported_dtypes, @@ -67,7 +66,6 @@ "corrcoef", "correlate", "cov", - "interp", "max", "mean", "median", @@ -1053,181 +1051,6 @@ def cov( ) -def interp(x, xp, fp, left=None, right=None, period=None): - """ - One-dimensional linear interpolation for monotonically increasing sample points. - - Returns the one-dimensional piecewise linear interpolant to a function - with given discrete data points (`xp`, `fp`), evaluated at `x`. - - Parameters - ---------- - x : array_like - The x-coordinates at which to evaluate the interpolated values. - - xp : 1-D sequence of floats - The x-coordinates of the data points, must be increasing if argument - `period` is not specified. Otherwise, `xp` is internally sorted after - normalizing the periodic boundaries with ``xp = xp % period``. - - fp : 1-D sequence of float or complex - The y-coordinates of the data points, same length as `xp`. - - left : optional float or complex corresponding to fp - Value to return for `x < xp[0]`, default is `fp[0]`. - - right : optional float or complex corresponding to fp - Value to return for `x > xp[-1]`, default is `fp[-1]`. - - period : None or float, optional - A period for the x-coordinates. This parameter allows the proper - interpolation of angular x-coordinates. Parameters `left` and `right` - are ignored if `period` is specified. - - Returns - ------- - y : float or complex (corresponding to fp) or ndarray - The interpolated values, same shape as `x`. - - Raises - ------ - ValueError - If `xp` and `fp` have different length - If `xp` or `fp` are not 1-D sequences - If `period == 0` - - See Also - -------- - scipy.interpolate - - Warnings - -------- - The x-coordinate sequence is expected to be increasing, but this is not - explicitly enforced. However, if the sequence `xp` is non-increasing, - interpolation results are meaningless. - - Note that, since NaN is unsortable, `xp` also cannot contain NaNs. - - A simple check for `xp` being strictly increasing is:: - - np.all(np.diff(xp) > 0) - - Examples - -------- - >>> import dpnp as np - >>> xp = np.array([1, 2, 3]) - >>> fp = np.array([3 ,2 ,0]) - >>> x = np.array([2.5]) - >>> np.interp(2.5, xp, fp) - 1.0 - >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp) - array([3. , 3. , 2.5 , 0.56, 0. ]) - >>> UNDEF = -99.0 - >>> np.interp(3.14, xp, fp, right=UNDEF) - -99.0 - - Plot an interpolant to the sine function: - - >>> x = np.linspace(0, 2*np.pi, 10) - >>> y = np.sin(x) - >>> xvals = np.linspace(0, 2*np.pi, 50) - >>> yinterp = np.interp(xvals, x, y) - >>> import matplotlib.pyplot as plt - >>> plt.plot(x, y, 'o') - [] - >>> plt.plot(xvals, yinterp, '-x') - [] - >>> plt.show() - - Interpolation with periodic x-coordinates: - - >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] - >>> xp = [190, -190, 350, -350] - >>> fp = [5, 10, 3, 4] - >>> np.interp(x, xp, fp, period=360) - array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75]) - - Complex interpolation: - - >>> x = [1.5, 4.0] - >>> xp = [2,3,5] - >>> fp = [1.0j, 0, 2+3j] - >>> np.interp(x, xp, fp) - array([0.+1.j , 1.+1.5j]) - - """ - - dpnp.check_supported_arrays_type(x, xp, fp) - - if xp.ndim != 1 or fp.ndim != 1: - raise ValueError("xp and fp must be 1D arrays") - if xp.size != fp.size: - raise ValueError("fp and xp are not of the same length") - if xp.size == 0: - raise ValueError("array of sample points is empty") - if not x.flags.c_contiguous: - raise NotImplementedError( - "Non-C-contiguous x is currently not supported" - ) - x_dtype = dpnp.common_type(x, xp) - if not dpnp.can_cast(x_dtype, dpnp.default_float_type()): - raise TypeError( - "Cannot cast array data from" - f" {x_dtype} to {dpnp.default_float_type()} " - "according to the rule 'safe'" - ) - - if period is not None: - # The handling of "period" below is modified from NumPy's - - if period == 0: - raise ValueError("period must be a non-zero value") - period = dpnp.abs(period) - left = None - right = None - - x = x.astype(dpnp.default_float_type()) - xp = xp.astype(dpnp.default_float_type()) - - # normalizing periodic boundaries - x %= period - xp %= period - asort_xp = dpnp.argsort(xp) - xp = xp[asort_xp] - fp = fp[asort_xp] - xp = dpnp.concatenate((xp[-1:] - period, xp, xp[0:1] + period)) - fp = dpnp.concatenate((fp[-1:], fp, fp[0:1])) - assert xp.flags.c_contiguous - assert fp.flags.c_contiguous - - # NumPy always returns float64 or complex128, so we upcast all values - # on the fly in the kernel - out_dtype = fp.dtype - output = dpnp.empty(x.shape, dtype=out_dtype) - idx = dpnp.searchsorted(xp, x, side="right") - left = fp[0] if left is None else dpnp.array(left, fp.dtype) - right = fp[-1] if right is None else dpnp.array(right, fp.dtype) - - idx = dpnp.array(idx, dtype="uint64") - - queue = x.sycl_queue - _manager = dpu.SequentialOrderManager[queue] - mem_ev, ht_ev = ufi._interpolate( - x.get_array(), - idx.get_array(), - xp.get_array(), - fp.get_array(), - # left, - # right, - output.get_array(), - queue, - depends=_manager.submitted_events, - ) - _manager.add_event_pair(mem_ev, ht_ev) - - return output - - def max(a, axis=None, out=None, keepdims=False, initial=None, where=True): """ Return the maximum of an array or maximum along an axis. From 7866eb85bc342edaa2a7ca7168732e7f7073e1ca Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 2 Apr 2025 03:52:18 -0700 Subject: [PATCH 06/26] Use dispatch vector and remove interpolate_complex_impl --- .../elementwise_functions/interpolate.cpp | 91 +++++++++--- .../elementwise_functions/interpolate.hpp | 134 +++++------------- 2 files changed, 107 insertions(+), 118 deletions(-) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp index 4de077b87997..b2573d26f5b9 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp @@ -23,6 +23,7 @@ // THE POSSIBILITY OF SUCH DAMAGE. //***************************************************************************** +#include #include #include @@ -43,6 +44,21 @@ namespace dpnp::extensions::ufunc namespace impl { +template +struct value_type_of +{ + using type = T; +}; + +template +struct value_type_of> +{ + using type = T; +}; + +template +using value_type_of_t = typename value_type_of::type; + typedef sycl::event (*interpolate_fn_ptr_t)(sycl::queue &, const void *, // x const void *, // idx @@ -53,8 +69,34 @@ typedef sycl::event (*interpolate_fn_ptr_t)(sycl::queue &, std::size_t, // xp_size const std::vector &); -interpolate_fn_ptr_t interpolate_dispatch_table[td_ns::num_types] - [td_ns::num_types]; +template +sycl::event interpolate_call(sycl::queue &exec_q, + const void *vx, + const void *vidx, + const void *vxp, + const void *vfp, + void *vout, + std::size_t n, + std::size_t xp_size, + const std::vector &depends) +{ + using dpctl::tensor::type_utils::is_complex_v; + using TCoord = std::conditional_t, value_type_of_t, T>; + + const TCoord *x = static_cast(vx); + const std::size_t *idx = static_cast(vidx); + const TCoord *xp = static_cast(vxp); + const T *fp = static_cast(vfp); + T *out = static_cast(vout); + + using dpnp::kernels::interpolate::interpolate_impl; + sycl::event interpolate_ev = interpolate_impl( + exec_q, x, idx, xp, fp, out, n, xp_size, depends); + + return interpolate_ev; +} + +interpolate_fn_ptr_t interpolate_dispatch_vector[td_ns::num_types]; std::pair py_interpolate(const dpctl::tensor::usm_ndarray &x, @@ -72,7 +114,7 @@ std::pair int xp_type_id = array_types.typenum_to_lookup_id(xp_typenum); int fp_type_id = array_types.typenum_to_lookup_id(fp_typenum); - auto fn = interpolate_dispatch_table[xp_type_id][fp_type_id]; + auto fn = interpolate_dispatch_vector[fp_type_id]; if (!fn) { throw py::type_error("Unsupported dtype."); } @@ -89,43 +131,54 @@ std::pair return std::make_pair(keep, ev); } -template +/** + * @brief A factory to define pairs of supported types for which + * interpolate function is available. + * + * @tparam T Type of input vector `a` and of result vector `y`. + */ +template +struct InterpolateOutputType +{ + 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; +}; + +template struct InterpolateFactory { fnT get() { - if constexpr (std::is_floating_point_v && - std::is_floating_point_v) - { - return dpnp::kernels::interpolate::interpolate_impl; - } - else if constexpr (std::is_floating_point_v && - (std::is_same_v> || - std::is_same_v>)) + if constexpr (std::is_same_v< + typename InterpolateOutputType::value_type, void>) { - return dpnp::kernels::interpolate::interpolate_complex_impl; + return nullptr; } else { - return nullptr; + return interpolate_call; } } }; -void init_interpolate_dispatch_table() +void init_interpolate_dispatch_vectors() { using namespace td_ns; - DispatchTableBuilder + DispatchVectorBuilder dtb_interpolate; - dtb_interpolate.populate_dispatch_table(interpolate_dispatch_table); + dtb_interpolate.populate_dispatch_vector(interpolate_dispatch_vector); } } // namespace impl void init_interpolate(py::module_ m) { - impl::init_interpolate_dispatch_table(); + impl::init_interpolate_dispatch_vectors(); using impl::py_interpolate; m.def("_interpolate", &py_interpolate, "", py::arg("x"), py::arg("idx"), diff --git a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp index f562c355fc37..17497f98c4f4 100644 --- a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -5,26 +5,44 @@ #include "utils/type_utils.hpp" +namespace type_utils = dpctl::tensor::type_utils; + namespace dpnp::kernels::interpolate { +template +struct IsNan +{ + static bool isnan(const T &v) + { + if constexpr (type_utils::is_complex_v) { + using vT = typename T::value_type; + + const vT real1 = std::real(v); + const vT imag1 = std::imag(v); + + return IsNan::isnan(real1) || IsNan::isnan(imag1); + } + else if constexpr (std::is_floating_point_v || + std::is_same_v) { + return sycl::isnan(v); + } + + return false; + } +}; + template sycl::event interpolate_impl(sycl::queue &q, - const void *vx, - const void *vidx, - const void *vxp, - const void *vfp, - void *vout, - std::size_t n, - std::size_t xp_size, + const TCoord *x, + const std::size_t *idx, + const TCoord *xp, + const TValue *fp, + TValue *out, + const std::size_t n, + const std::size_t xp_size, const std::vector &depends) { - const TCoord *x = static_cast(vx); - const std::size_t *idx = static_cast(vidx); - const TCoord *xp = static_cast(vxp); - const TValue *fp = static_cast(vfp); - TValue *out = static_cast(vout); - return q.submit([&](sycl::handler &h) { h.depends_on(depends); h.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { @@ -34,7 +52,7 @@ sycl::event interpolate_impl(sycl::queue &q, TCoord x_val = x[i]; std::size_t x_idx = idx[i] - 1; - if (sycl::isnan(x_val)) { + if (IsNan::isnan(x_val)) { out[i] = x_val; } else if (x_idx < 0) { @@ -54,9 +72,10 @@ sycl::event interpolate_impl(sycl::queue &q, (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); TValue res = slope * (x_val - xp[x_idx]) + fp[x_idx]; - if (sycl::isnan(res)) { + if (IsNan::isnan(res)) { res = slope * (x_val - xp[x_idx + 1]) + fp[x_idx + 1]; - if (sycl::isnan(res) && (fp[x_idx] == fp[x_idx + 1])) { + if (IsNan::isnan(res) && + (fp[x_idx] == fp[x_idx + 1])) { res = fp[x_idx]; } } @@ -66,87 +85,4 @@ sycl::event interpolate_impl(sycl::queue &q, }); } -template -sycl::event interpolate_complex_impl(sycl::queue &q, - const void *vx, - const void *vidx, - const void *vxp, - const void *vfp, - void *vout, - std::size_t n, - std::size_t xp_size, - const std::vector &depends) -{ - const TCoord *x = static_cast(vx); - const std::size_t *idx = static_cast(vidx); - const TCoord *xp = static_cast(vxp); - const TValue *fp = static_cast(vfp); - TValue *out = static_cast(vout); - - using realT = typename TValue::value_type; - - return q.submit([&](sycl::handler &h) { - h.depends_on(depends); - h.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { - realT left_r = fp[0].real(); - realT right_r = fp[xp_size - 1].real(); - realT left_i = fp[0].imag(); - realT right_i = fp[xp_size - 1].imag(); - - TCoord x_val = x[i]; - std::size_t x_idx = idx[i] - 1; - - realT res_r = 0.0; - realT res_i = 0.0; - - if (sycl::isnan(x_val)) { - res_r = x_val; - res_i = 0.0; - } - else if (x_idx < 0) { - res_r = left_r; - res_i = left_i; - } - else if (x_val == xp[xp_size - 1]) { - res_r = right_r; - res_i = right_i; - } - else if (x_idx >= xp_size - 1) { - res_r = right_r; - res_i = right_i; - } - else if (x_val == xp[x_idx]) { - res_r = fp[x_idx].real(); - res_i = fp[x_idx].imag(); - } - else { - realT dx = xp[x_idx + 1] - xp[x_idx]; - - realT slope_r = (fp[x_idx + 1].real() - fp[x_idx].real()) / dx; - res_r = slope_r * (x_val - xp[x_idx]) + fp[x_idx].real(); - if (sycl::isnan(res_r)) { - res_r = slope_r * (x_val - xp[x_idx + 1]) + - fp[x_idx + 1].real(); - if (sycl::isnan(res_r) && - fp[x_idx].real() == fp[x_idx + 1].real()) { - res_r = fp[x_idx].real(); - } - } - - realT slope_i = (fp[x_idx + 1].imag() - fp[x_idx].imag()) / dx; - res_i = slope_i * (x_val - xp[x_idx]) + fp[x_idx].imag(); - if (sycl::isnan(res_i)) { - res_i = slope_i * (x_val - xp[x_idx + 1]) + - fp[x_idx + 1].imag(); - if (sycl::isnan(res_i) && - fp[x_idx].imag() == fp[x_idx + 1].imag()) { - res_i = fp[x_idx].imag(); - } - } - } - out[i] = TValue(res_r, res_i); - }); - }); -} - } // namespace dpnp::kernels::interpolate From 51b3bde065f47aab6ec47fb1ab1245348ceea1eb Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 2 Apr 2025 06:59:02 -0700 Subject: [PATCH 07/26] Add more backend checks --- .../elementwise_functions/interpolate.cpp | 37 ++++++++++++++++++- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp index b2573d26f5b9..cd6380babe2d 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp @@ -26,11 +26,12 @@ #include #include +#include "dpctl4pybind11.hpp" #include #include // dpctl tensor headers -#include "dpctl4pybind11.hpp" +#include "utils/output_validation.hpp" #include "utils/type_dispatch.hpp" #include "kernels/elementwise_functions/interpolate.hpp" @@ -107,16 +108,48 @@ std::pair sycl::queue &exec_q, const std::vector &depends) { + int x_typenum = x.get_typenum(); int xp_typenum = xp.get_typenum(); int fp_typenum = fp.get_typenum(); + int out_typenum = out.get_typenum(); auto array_types = td_ns::usm_ndarray_types(); + int x_type_id = array_types.typenum_to_lookup_id(x_typenum); int xp_type_id = array_types.typenum_to_lookup_id(xp_typenum); int fp_type_id = array_types.typenum_to_lookup_id(fp_typenum); + int out_type_id = array_types.typenum_to_lookup_id(out_typenum); + + if (x_type_id != xp_type_id) { + throw py::value_error("x and xp must have the same dtype"); + } + if (fp_type_id != out_type_id) { + throw py::value_error("fp and out must have the same dtype"); + } auto fn = interpolate_dispatch_vector[fp_type_id]; if (!fn) { - throw py::type_error("Unsupported dtype."); + throw py::type_error("Unsupported dtype"); + } + + if (!dpctl::utils::queues_are_compatible(exec_q, {x, idx, xp, fp, out})) { + throw py::value_error( + "Execution queue is not compatible with allocation queues"); + } + + dpctl::tensor::validation::CheckWritable::throw_if_not_writable(out); + + if (x.get_ndim() != 1 || xp.get_ndim() != 1 || fp.get_ndim() != 1 || + idx.get_ndim() != 1 || out.get_ndim() != 1) + { + throw py::value_error("All arrays must be one-dimensional"); + } + + if (xp.get_size() != fp.get_size()) { + throw py::value_error("xp and fp must have the same size"); + } + + if (x.get_size() != out.get_size() || x.get_size() != idx.get_size()) { + throw py::value_error("x, idx, and out must have the same size"); } std::size_t n = x.get_size(); From ecfa37dfafd493e2d92b6ee6f603241f7277df8d Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 10 Apr 2025 04:16:42 -0700 Subject: [PATCH 08/26] Add support left/right args --- .../elementwise_functions/interpolate.cpp | 44 ++++++++++++++++--- .../elementwise_functions/interpolate.hpp | 12 ++--- dpnp/dpnp_iface_mathematical.py | 12 +++-- 3 files changed, 52 insertions(+), 16 deletions(-) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp index cd6380babe2d..8190ffcd929c 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp @@ -65,6 +65,8 @@ typedef sycl::event (*interpolate_fn_ptr_t)(sycl::queue &, const void *, // idx const void *, // xp const void *, // fp + const void *, // left + const void *, // right void *, // out std::size_t, // n std::size_t, // xp_size @@ -76,6 +78,8 @@ sycl::event interpolate_call(sycl::queue &exec_q, const void *vidx, const void *vxp, const void *vfp, + const void *vleft, + const void *vright, void *vout, std::size_t n, std::size_t xp_size, @@ -88,11 +92,13 @@ sycl::event interpolate_call(sycl::queue &exec_q, const std::size_t *idx = static_cast(vidx); const TCoord *xp = static_cast(vxp); const T *fp = static_cast(vfp); + const T *left = static_cast(vleft); + const T *right = static_cast(vright); T *out = static_cast(vout); using dpnp::kernels::interpolate::interpolate_impl; sycl::event interpolate_ev = interpolate_impl( - exec_q, x, idx, xp, fp, out, n, xp_size, depends); + exec_q, x, idx, xp, fp, left, right, out, n, xp_size, depends); return interpolate_ev; } @@ -104,6 +110,8 @@ std::pair const dpctl::tensor::usm_ndarray &idx, const dpctl::tensor::usm_ndarray &xp, const dpctl::tensor::usm_ndarray &fp, + std::optional &left, + std::optional &right, dpctl::tensor::usm_ndarray &out, sycl::queue &exec_q, const std::vector &depends) @@ -155,13 +163,34 @@ std::pair std::size_t n = x.get_size(); std::size_t xp_size = xp.get_size(); - sycl::event ev = fn(exec_q, x.get_data(), idx.get_data(), xp.get_data(), - fp.get_data(), out.get_data(), n, xp_size, depends); + void *left_ptr = left.has_value() ? left.value().get_data() : nullptr; - sycl::event keep = - dpctl::utils::keep_args_alive(exec_q, {x, idx, xp, fp, out}, {ev}); + void *right_ptr = right.has_value() ? right.value().get_data() : nullptr; - return std::make_pair(keep, ev); + sycl::event ev = + fn(exec_q, x.get_data(), idx.get_data(), xp.get_data(), fp.get_data(), + left_ptr, right_ptr, out.get_data(), n, xp_size, depends); + + sycl::event args_ev; + + if (left.has_value() && right.has_value()) { + args_ev = dpctl::utils::keep_args_alive( + exec_q, {x, idx, xp, fp, out, left.value(), right.value()}, {ev}); + } + else if (left.has_value()) { + args_ev = dpctl::utils::keep_args_alive( + exec_q, {x, idx, xp, fp, out, left.value()}, {ev}); + } + else if (right.has_value()) { + args_ev = dpctl::utils::keep_args_alive( + exec_q, {x, idx, xp, fp, out, right.value()}, {ev}); + } + else { + args_ev = + dpctl::utils::keep_args_alive(exec_q, {x, idx, xp, fp, out}, {ev}); + } + + return std::make_pair(args_ev, ev); } /** @@ -215,7 +244,8 @@ void init_interpolate(py::module_ m) using impl::py_interpolate; m.def("_interpolate", &py_interpolate, "", py::arg("x"), py::arg("idx"), - py::arg("xp"), py::arg("fp"), py::arg("out"), py::arg("sycl_queue"), + py::arg("xp"), py::arg("fp"), py::arg("left"), py::arg("right"), + py::arg("out"), py::arg("sycl_queue"), py::arg("depends") = py::list()); } diff --git a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp index 17497f98c4f4..7eb974515f0e 100644 --- a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -38,6 +38,8 @@ sycl::event interpolate_impl(sycl::queue &q, const std::size_t *idx, const TCoord *xp, const TValue *fp, + const TValue *left, + const TValue *right, TValue *out, const std::size_t n, const std::size_t xp_size, @@ -46,8 +48,8 @@ sycl::event interpolate_impl(sycl::queue &q, return q.submit([&](sycl::handler &h) { h.depends_on(depends); h.parallel_for(sycl::range<1>(n), [=](sycl::id<1> i) { - TValue left = fp[0]; - TValue right = fp[xp_size - 1]; + TValue left_val = left ? *left : fp[0]; + TValue right_val = right ? *right : fp[xp_size - 1]; TCoord x_val = x[i]; std::size_t x_idx = idx[i] - 1; @@ -56,13 +58,13 @@ sycl::event interpolate_impl(sycl::queue &q, out[i] = x_val; } else if (x_idx < 0) { - out[i] = left; + out[i] = left_val; } else if (x_val == xp[xp_size - 1]) { - out[i] = right; + out[i] = right_val; } else if (x_idx >= xp_size - 1) { - out[i] = right; + out[i] = right_val; } else if (x_val == xp[x_idx]) { out[i] = fp[x_idx]; diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 551e28077f55..c72be3c1dd5b 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2842,8 +2842,12 @@ def interp(x, xp, fp, left=None, right=None, period=None): out_dtype = fp.dtype output = dpnp.empty(x.shape, dtype=out_dtype) idx = dpnp.searchsorted(xp, x, side="right") - left = fp[0] if left is None else dpnp.array(left, fp.dtype) - right = fp[-1] if right is None else dpnp.array(right, fp.dtype) + left_usm = ( + dpnp.array(left, fp.dtype).get_array() if left is not None else None + ) + right_usm = ( + dpnp.array(right, fp.dtype).get_array() if right is not None else None + ) idx = dpnp.array(idx, dtype="uint64") @@ -2854,8 +2858,8 @@ def interp(x, xp, fp, left=None, right=None, period=None): idx.get_array(), xp.get_array(), fp.get_array(), - # left, - # right, + left_usm, + right_usm, output.get_array(), queue, depends=_manager.submitted_events, From 5d53f9c1cdb8469961edd45e18db2199fc80c3a6 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 10 Apr 2025 04:25:27 -0700 Subject: [PATCH 09/26] Use get_usm_allocations in def interp --- dpnp/dpnp_iface_mathematical.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index c72be3c1dd5b..b5f0c16f3118 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2806,11 +2806,15 @@ def interp(x, xp, fp, left=None, right=None, period=None): raise NotImplementedError( "Non-C-contiguous x is currently not supported" ) + _, exec_q = get_usm_allocations([x, xp, fp]) + x_dtype = dpnp.common_type(x, xp) - if not dpnp.can_cast(x_dtype, dpnp.default_float_type()): + x_float_type = dpnp.default_float_type(exec_q) + + if not dpnp.can_cast(x_dtype, x_float_type): raise TypeError( "Cannot cast array data from" - f" {x_dtype} to {dpnp.default_float_type()} " + f" {x_dtype} to {x_float_type} " "according to the rule 'safe'" ) @@ -2823,8 +2827,8 @@ def interp(x, xp, fp, left=None, right=None, period=None): left = None right = None - x = x.astype(dpnp.default_float_type()) - xp = xp.astype(dpnp.default_float_type()) + x = x.astype(x_float_type) + xp = xp.astype(x_float_type) # normalizing periodic boundaries x %= period @@ -2849,10 +2853,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): dpnp.array(right, fp.dtype).get_array() if right is not None else None ) - idx = dpnp.array(idx, dtype="uint64") - - queue = x.sycl_queue - _manager = dpu.SequentialOrderManager[queue] + _manager = dpu.SequentialOrderManager[exec_q] mem_ev, ht_ev = ufi._interpolate( x.get_array(), idx.get_array(), @@ -2861,7 +2862,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): left_usm, right_usm, output.get_array(), - queue, + exec_q, depends=_manager.submitted_events, ) _manager.add_event_pair(mem_ev, ht_ev) From 9dbc2c5d5e3736f32463b857b91d53907e45330f Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 03:57:36 -0700 Subject: [PATCH 10/26] Pass idx as std::int64_t --- .../ufunc/elementwise_functions/interpolate.cpp | 2 +- .../kernels/elementwise_functions/interpolate.hpp | 9 +++------ 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp index 8190ffcd929c..7e6cf7b4613e 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp @@ -89,7 +89,7 @@ sycl::event interpolate_call(sycl::queue &exec_q, using TCoord = std::conditional_t, value_type_of_t, T>; const TCoord *x = static_cast(vx); - const std::size_t *idx = static_cast(vidx); + const std::int64_t *idx = static_cast(vidx); const TCoord *xp = static_cast(vxp); const T *fp = static_cast(vfp); const T *left = static_cast(vleft); diff --git a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp index 7eb974515f0e..1c12f051cad4 100644 --- a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -35,7 +35,7 @@ struct IsNan template sycl::event interpolate_impl(sycl::queue &q, const TCoord *x, - const std::size_t *idx, + const std::int64_t *idx, const TCoord *xp, const TValue *fp, const TValue *left, @@ -52,7 +52,7 @@ sycl::event interpolate_impl(sycl::queue &q, TValue right_val = right ? *right : fp[xp_size - 1]; TCoord x_val = x[i]; - std::size_t x_idx = idx[i] - 1; + std::int64_t x_idx = idx[i] - 1; if (IsNan::isnan(x_val)) { out[i] = x_val; @@ -63,12 +63,9 @@ sycl::event interpolate_impl(sycl::queue &q, else if (x_val == xp[xp_size - 1]) { out[i] = right_val; } - else if (x_idx >= xp_size - 1) { + else if (x_idx >= static_cast(xp_size - 1)) { out[i] = right_val; } - else if (x_val == xp[x_idx]) { - out[i] = fp[x_idx]; - } else { TValue slope = (fp[x_idx + 1] - fp[x_idx]) / (xp[x_idx + 1] - xp[x_idx]); From 1bafd7c614dc0290216ad29f4631c2f7b70565b9 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 04:01:09 -0700 Subject: [PATCH 11/26] Add proper casting input array --- dpnp/dpnp_iface_mathematical.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index b5f0c16f3118..ad17316c25a3 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2818,6 +2818,13 @@ def interp(x, xp, fp, left=None, right=None, period=None): "according to the rule 'safe'" ) + x = dpnp.asarray(x, dtype=x_float_type, order="C") + xp = dpnp.asarray(xp, dtype=x_float_type, order="C") + + out_dtype = dpnp.common_type(x, xp, fp) + + fp = dpnp.asarray(fp, dtype=out_dtype, order="C") + if period is not None: # The handling of "period" below is modified from NumPy's @@ -2841,9 +2848,6 @@ def interp(x, xp, fp, left=None, right=None, period=None): assert xp.flags.c_contiguous assert fp.flags.c_contiguous - # NumPy always returns float64 or complex128, so we upcast all values - # on the fly in the kernel - out_dtype = fp.dtype output = dpnp.empty(x.shape, dtype=out_dtype) idx = dpnp.searchsorted(xp, x, side="right") left_usm = ( From 2f43fd78c7c881fc361238a04dc5006100bb5f01 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 05:43:37 -0700 Subject: [PATCH 12/26] Update def interp to support period args --- dpnp/dpnp_iface_mathematical.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index ad17316c25a3..1fc2866c7f9f 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2806,7 +2806,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): raise NotImplementedError( "Non-C-contiguous x is currently not supported" ) - _, exec_q = get_usm_allocations([x, xp, fp]) + usm_type, exec_q = get_usm_allocations([x, xp, fp]) x_dtype = dpnp.common_type(x, xp) x_float_type = dpnp.default_float_type(exec_q) @@ -2828,6 +2828,15 @@ def interp(x, xp, fp, left=None, right=None, period=None): if period is not None: # The handling of "period" below is modified from NumPy's + if dpnp.is_supported_array_type(period): + if dpu.get_execution_queue([exec_q, period.sycl_queue]) is None: + raise ValueError( + "input arrays and period must be allocated " + "on the same SYCL queue" + ) + else: + period = dpnp.asarray(period, sycl_queue=exec_q, usm_type=usm_type) + if period == 0: raise ValueError("period must be a non-zero value") period = dpnp.abs(period) From ae65091bd5284e7b2b121a56557fcc141ef6163d Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 06:05:02 -0700 Subject: [PATCH 13/26] Return fp[-1] instead of right_val for x==xp[-1] --- dpnp/backend/kernels/elementwise_functions/interpolate.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp index 1c12f051cad4..fd9e64d5a75a 100644 --- a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -61,7 +61,7 @@ sycl::event interpolate_impl(sycl::queue &q, out[i] = left_val; } else if (x_val == xp[xp_size - 1]) { - out[i] = right_val; + out[i] = fp[xp_size - 1]; } else if (x_idx >= static_cast(xp_size - 1)) { out[i] = right_val; From 771d3ebb6e337bdc801a5ec900bf22fedb57e3bc Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 07:36:07 -0700 Subject: [PATCH 14/26] Unskip cupy tests for interp --- dpnp/tests/third_party/cupy/math_tests/test_misc.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/dpnp/tests/third_party/cupy/math_tests/test_misc.py b/dpnp/tests/third_party/cupy/math_tests/test_misc.py index 7746c56d3253..aeebd90142a0 100644 --- a/dpnp/tests/third_party/cupy/math_tests/test_misc.py +++ b/dpnp/tests/third_party/cupy/math_tests/test_misc.py @@ -367,7 +367,6 @@ def test_real_if_close_with_float_tol_false(self, xp, dtype): assert x.dtype == out.dtype return out - @pytest.mark.skip("interp() is not supported yet") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) @testing.numpy_cupy_allclose(atol=1e-5) @@ -378,7 +377,6 @@ def test_interp(self, xp, dtype_y, dtype_x): fy = xp.sin(fx).astype(dtype_y) return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) @testing.numpy_cupy_allclose(atol=1e-5) @@ -389,7 +387,6 @@ def test_interp_period(self, xp, dtype_y, dtype_x): fy = xp.sin(fx).astype(dtype_y) return xp.interp(x, fx, fy, period=5) - @pytest.mark.skip("interp() is not supported yet") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) @testing.numpy_cupy_allclose(atol=1e-5) @@ -402,7 +399,6 @@ def test_interp_left_right(self, xp, dtype_y, dtype_x): right = 20 return xp.interp(x, fx, fy, left, right) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_dtypes("efdFD", name="dtype_y") @@ -415,7 +411,6 @@ def test_interp_nan_fy(self, xp, dtype_y, dtype_x): fy[0] = fy[2] = fy[-1] = numpy.nan return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") @@ -428,7 +423,6 @@ def test_interp_nan_fx(self, xp, dtype_y, dtype_x): fx[-1] = numpy.nan # x and fx must remain sorted (NaNs are the last) return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") @@ -441,7 +435,6 @@ def test_interp_nan_x(self, xp, dtype_y, dtype_x): x[-1] = numpy.nan # x and fx must remain sorted (NaNs are the last) return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_dtypes("efdFD", name="dtype_y") @@ -454,7 +447,6 @@ def test_interp_inf_fy(self, xp, dtype_y, dtype_x): fy[0] = fy[2] = fy[-1] = numpy.inf return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") @@ -467,7 +459,6 @@ def test_interp_inf_fx(self, xp, dtype_y, dtype_x): fx[-1] = numpy.inf # x and fx must remain sorted return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") @@ -480,7 +471,6 @@ def test_interp_inf_x(self, xp, dtype_y, dtype_x): x[-1] = numpy.inf # x and fx must remain sorted return xp.interp(x, fx, fy) - @pytest.mark.skip("interp() is not supported yet") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) @testing.numpy_cupy_allclose(atol=1e-5) @@ -493,7 +483,6 @@ def test_interp_size1(self, xp, dtype_y, dtype_x): right = 20 return xp.interp(x, fx, fy, left, right) - @pytest.mark.skip("interp() is not supported yet") @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") From 5cda3d2cc07f0b01f0de6d64ac70b959a360b9f7 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 07:37:51 -0700 Subject: [PATCH 15/26] Add dpnp tests for interp --- dpnp/tests/test_mathematical.py | 138 ++++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) diff --git a/dpnp/tests/test_mathematical.py b/dpnp/tests/test_mathematical.py index c416a55e5796..91c13a579e0d 100644 --- a/dpnp/tests/test_mathematical.py +++ b/dpnp/tests/test_mathematical.py @@ -1143,6 +1143,144 @@ def test_complex(self, xp): assert_raises((ValueError, TypeError), xp.i0, a) +class TestInterp: + @pytest.mark.parametrize( + "dtype_x", get_all_dtypes(no_bool=True, no_complex=True) + ) + @pytest.mark.parametrize("dtype_y", get_all_dtypes(no_bool=True)) + def test_all_dtypes(self, dtype_x, dtype_y): + x = numpy.linspace(0.1, 9.9, 20).astype(dtype_x) + xp = numpy.linspace(0.0, 10.0, 5).astype(dtype_x) + fp = (xp * 1.5 + 1).astype(dtype_y) + + ix = dpnp.array(x) + ixp = dpnp.array(xp) + ifp = dpnp.array(fp) + + expected = numpy.interp(x, xp, fp) + result = dpnp.interp(ix, ixp, ifp) + assert_dtype_allclose(result, expected) + + @pytest.mark.parametrize( + "dtype_x", get_all_dtypes(no_bool=True, no_complex=True) + ) + @pytest.mark.parametrize("dtype_y", get_complex_dtypes()) + def test_complex_fp(self, dtype_x, dtype_y): + x = numpy.array([0.25, 0.75], dtype=dtype_x) + xp = numpy.array([0.0, 1.0], dtype=dtype_x) + fp = numpy.array([1 + 1j, 3 + 3j], dtype=dtype_y) + + ix = dpnp.array(x) + ixp = dpnp.array(xp) + ifp = dpnp.array(fp) + + expected = numpy.interp(x, xp, fp) + result = dpnp.interp(ix, ixp, ifp) + assert_dtype_allclose(result, expected) + + @pytest.mark.parametrize( + "dtype", get_all_dtypes(no_bool=True, no_complex=True) + ) + def test_left_right_args(self, dtype): + x = numpy.array([-1, 0, 1, 2, 3, 4, 5, 6], dtype=dtype) + xp = numpy.array([0, 3, 6], dtype=dtype) + fp = numpy.array([0, 9, 18], dtype=dtype) + + ix = dpnp.array(x) + ixp = dpnp.array(xp) + ifp = dpnp.array(fp) + + expected = numpy.interp(x, xp, fp, left=-40, right=40) + result = dpnp.interp(ix, ixp, ifp, left=-40, right=40) + assert_dtype_allclose(result, expected) + + @pytest.mark.parametrize("val", [numpy.nan, numpy.inf, -numpy.inf]) + def test_naninf(self, val): + x = numpy.array([0, 1, 2, val]) + xp = numpy.array([0, 1, 2]) + fp = numpy.array([10, 20, 30]) + + ix = dpnp.array(x) + ixp = dpnp.array(xp) + ifp = dpnp.array(fp) + + expected = numpy.interp(x, xp, fp) + result = dpnp.interp(ix, ixp, ifp) + assert_dtype_allclose(result, expected) + + def test_empty_x(self): + x = numpy.array([]) + xp = numpy.array([0, 1]) + fp = numpy.array([10, 20]) + + ix = dpnp.array(x) + ixp = dpnp.array(xp) + ifp = dpnp.array(fp) + + expected = numpy.interp(x, xp, fp) + result = dpnp.interp(ix, ixp, ifp) + assert_dtype_allclose(result, expected) + + @pytest.mark.parametrize("dtype", get_float_dtypes()) + def test_period(self, dtype): + x = numpy.array([-180, 0, 180], dtype=dtype) + xp = numpy.array([-90, 0, 90], dtype=dtype) + fp = numpy.array([0, 1, 0], dtype=dtype) + + ix = dpnp.array(x) + ixp = dpnp.array(xp) + ifp = dpnp.array(fp) + + expected = numpy.interp(x, xp, fp, period=180) + result = dpnp.interp(ix, ixp, ifp, period=180) + assert_dtype_allclose(result, expected) + + def test_errors(self): + x = dpnp.array([0.5]) + + # xp and fp have different lengths + xp = dpnp.array([0]) + fp = dpnp.array([1, 2]) + assert_raises(ValueError, dpnp.interp, x, xp, fp) + + # xp is not 1D + xp = dpnp.array([[0, 1]]) + fp = dpnp.array([1, 2]) + assert_raises(ValueError, dpnp.interp, x, xp, fp) + + # fp is not 1D + xp = dpnp.array([0, 1]) + fp = dpnp.array([[1, 2]]) + assert_raises(ValueError, dpnp.interp, x, xp, fp) + + # xp and fp are empty + xp = dpnp.array([]) + fp = dpnp.array([]) + assert_raises(ValueError, dpnp.interp, x, xp, fp) + + # x complex + x_complex = dpnp.array([1 + 2j]) + xp = dpnp.array([0.0, 2.0]) + fp = dpnp.array([0.0, 1.0]) + assert_raises(TypeError, dpnp.interp, x_complex, xp, fp) + + # period is zero + x = dpnp.array([1.0]) + xp = dpnp.array([0.0, 2.0]) + fp = dpnp.array([0.0, 1.0]) + assert_raises(ValueError, dpnp.interp, x, xp, fp, period=0) + + # period has a different SYCL queue + q1 = dpctl.SyclQueue() + q2 = dpctl.SyclQueue() + + x = dpnp.array([1.0], sycl_queue=q1) + xp = dpnp.array([0.0, 2.0], sycl_queue=q1) + fp = dpnp.array([0.0, 1.0], sycl_queue=q1) + period = dpnp.array([180], sycl_queue=q2) + assert_raises(ValueError, dpnp.interp, x, xp, fp, period=period) + + @pytest.mark.parametrize( "rhs", [[[1, 2, 3], [4, 5, 6]], [2.0, 1.5, 1.0], 3, 0.3] ) From a65a1dd92e40599f3db541f5bfbe76977cfea56e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 07:55:44 -0700 Subject: [PATCH 16/26] Update docstrings for def interp() --- dpnp/dpnp_iface_mathematical.py | 96 ++++++++++++++------------------- 1 file changed, 41 insertions(+), 55 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 1fc2866c7f9f..693b4583f33c 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2692,50 +2692,51 @@ def gradient(f, *varargs, axis=None, edge_order=1): def interp(x, xp, fp, left=None, right=None, period=None): """ - One-dimensional linear interpolation for monotonically increasing sample points. + One-dimensional linear interpolation. Returns the one-dimensional piecewise linear interpolant to a function with given discrete data points (`xp`, `fp`), evaluated at `x`. + For full documentation refer to :obj:`numpy.interp`. + Parameters ---------- - x : array_like - The x-coordinates at which to evaluate the interpolated values. + x : {dpnp.ndarray, usm_ndarray} + Input 1-D array. The x-coordinates at which to evaluate + the interpolated values. + + xp : {dpnp.ndarray, usm_ndarray} + Input 1-D array. The x-coordinates of the data points, + must be increasing if argument `period` is not specified. + Otherwise, `xp` is internally sorted after normalizing + the periodic boundaries with ``xp = xp % period``. + + fp : {dpnp.ndarray, usm_ndarray} + Input 1-D array. The y-coordinates of the data points, + same length as `xp`. - xp : 1-D sequence of floats - The x-coordinates of the data points, must be increasing if argument - `period` is not specified. Otherwise, `xp` is internally sorted after - normalizing the periodic boundaries with ``xp = xp % period``. + left : {None, scalar, dpnp.ndarray, usm_ndarray}, optional + Value to return for `x < xp[0]`. - fp : 1-D sequence of float or complex - The y-coordinates of the data points, same length as `xp`. + Default: ``fp[0]``. - left : optional float or complex corresponding to fp - Value to return for `x < xp[0]`, default is `fp[0]`. + right : {None, scalar, dpnp.ndarray, usm_ndarray}, optional + Value to return for `x > xp[-1]`. - right : optional float or complex corresponding to fp - Value to return for `x > xp[-1]`, default is `fp[-1]`. + Default: ``fp[-1]``. - period : None or float, optional + period : {None, scalar, dpnp.ndarray, usm_ndarray}, optional A period for the x-coordinates. This parameter allows the proper interpolation of angular x-coordinates. Parameters `left` and `right` are ignored if `period` is specified. + Default: ``None``. + Returns ------- - y : float or complex (corresponding to fp) or ndarray + y : {dpnp.ndarray, usm_ndarray} The interpolated values, same shape as `x`. - Raises - ------ - ValueError - If `xp` and `fp` have different length - If `xp` or `fp` are not 1-D sequences - If `period == 0` - - See Also - -------- - scipy.interpolate Warnings -------- @@ -2747,6 +2748,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): A simple check for `xp` being strictly increasing is:: + import dpnp as np np.all(np.diff(xp) > 0) Examples @@ -2755,40 +2757,29 @@ def interp(x, xp, fp, left=None, right=None, period=None): >>> xp = np.array([1, 2, 3]) >>> fp = np.array([3 ,2 ,0]) >>> x = np.array([2.5]) - >>> np.interp(2.5, xp, fp) - 1.0 - >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp) + >>> np.interp(x, xp, fp) + array([1.]) + >>> x = np.array([0, 1, 1.5, 2.72, 3.14]) + >>> np.interp(x, xp, fp) array([3. , 3. , 2.5 , 0.56, 0. ]) + >>> x = np.array([3.14]) >>> UNDEF = -99.0 - >>> np.interp(3.14, xp, fp, right=UNDEF) - -99.0 - - Plot an interpolant to the sine function: - - >>> x = np.linspace(0, 2*np.pi, 10) - >>> y = np.sin(x) - >>> xvals = np.linspace(0, 2*np.pi, 50) - >>> yinterp = np.interp(xvals, x, y) - >>> import matplotlib.pyplot as plt - >>> plt.plot(x, y, 'o') - [] - >>> plt.plot(xvals, yinterp, '-x') - [] - >>> plt.show() + >>> np.interp(x, xp, fp, right=UNDEF) + array([-99.]) Interpolation with periodic x-coordinates: - >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] - >>> xp = [190, -190, 350, -350] - >>> fp = [5, 10, 3, 4] + >>> x = np.array([-180, -170, -185, 185, -10, -5, 0, 365]) + >>> xp = np.array([190, -190, 350, -350]) + >>> fp = np.array([5, 10, 3, 4]) >>> np.interp(x, xp, fp, period=360) array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75]) Complex interpolation: - >>> x = [1.5, 4.0] - >>> xp = [2,3,5] - >>> fp = [1.0j, 0, 2+3j] + >>> x = np.array([1.5, 4.0]) + >>> xp = np.array([2,3,5]) + >>> fp = np.array([1.0j, 0, 2+3j]) >>> np.interp(x, xp, fp) array([0.+1.j , 1.+1.5j]) @@ -2802,10 +2793,7 @@ def interp(x, xp, fp, left=None, right=None, period=None): raise ValueError("fp and xp are not of the same length") if xp.size == 0: raise ValueError("array of sample points is empty") - if not x.flags.c_contiguous: - raise NotImplementedError( - "Non-C-contiguous x is currently not supported" - ) + usm_type, exec_q = get_usm_allocations([x, xp, fp]) x_dtype = dpnp.common_type(x, xp) @@ -2826,8 +2814,6 @@ def interp(x, xp, fp, left=None, right=None, period=None): fp = dpnp.asarray(fp, dtype=out_dtype, order="C") if period is not None: - # The handling of "period" below is modified from NumPy's - if dpnp.is_supported_array_type(period): if dpu.get_execution_queue([exec_q, period.sycl_queue]) is None: raise ValueError( From 99cc8b5147e5ddc9415f4aa94309f071c218d130 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 11 Apr 2025 08:53:41 -0700 Subject: [PATCH 17/26] Remove lines after merging --- dpnp/tests/third_party/cupy/test_type_routines.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/dpnp/tests/third_party/cupy/test_type_routines.py b/dpnp/tests/third_party/cupy/test_type_routines.py index bf5c7af9ded0..e35b40d90841 100644 --- a/dpnp/tests/third_party/cupy/test_type_routines.py +++ b/dpnp/tests/third_party/cupy/test_type_routines.py @@ -47,10 +47,6 @@ def test_can_cast(self, xp, from_dtype, to_dtype): return ret -<<<<<<< HEAD -======= -# @pytest.mark.skip("dpnp.common_type() is not implemented yet") ->>>>>>> e8871f0d797 (Second impl with dispatch_vector[only floating]) class TestCommonType(unittest.TestCase): @testing.numpy_cupy_equal() From 1263eb503d2e4a001ab9c5e3cf6ee38fffeb11f9 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 14 Apr 2025 02:45:18 -0700 Subject: [PATCH 18/26] Add type_check flag to cupy tests --- .../third_party/cupy/math_tests/test_misc.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/dpnp/tests/third_party/cupy/math_tests/test_misc.py b/dpnp/tests/third_party/cupy/math_tests/test_misc.py index aeebd90142a0..4542c51de33e 100644 --- a/dpnp/tests/third_party/cupy/math_tests/test_misc.py +++ b/dpnp/tests/third_party/cupy/math_tests/test_misc.py @@ -369,7 +369,7 @@ def test_real_if_close_with_float_tol_false(self, xp, dtype): @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -379,7 +379,7 @@ def test_interp(self, xp, dtype_y, dtype_x): @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_period(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -389,7 +389,7 @@ def test_interp_period(self, xp, dtype_y, dtype_x): @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_left_right(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -402,7 +402,7 @@ def test_interp_left_right(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_nan_fy(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -414,7 +414,7 @@ def test_interp_nan_fy(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_nan_fx(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -426,7 +426,7 @@ def test_interp_nan_fx(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_nan_x(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -438,7 +438,7 @@ def test_interp_nan_x(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_inf_fy(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -450,7 +450,7 @@ def test_interp_inf_fy(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_inf_fx(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -462,7 +462,7 @@ def test_interp_inf_fx(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_inf_x(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -473,7 +473,7 @@ def test_interp_inf_x(self, xp, dtype_y, dtype_x): @testing.for_all_dtypes(name="dtype_x", no_bool=True, no_complex=True) @testing.for_all_dtypes(name="dtype_y", no_bool=True) - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_size1(self, xp, dtype_y, dtype_x): # interpolate at points on and outside the boundaries x = xp.asarray([0, 1, 2, 4, 6, 8, 9, 10], dtype=dtype_x) @@ -486,7 +486,7 @@ def test_interp_size1(self, xp, dtype_y, dtype_x): @testing.with_requires("numpy>=1.17.0") @testing.for_float_dtypes(name="dtype_x") @testing.for_dtypes("efdFD", name="dtype_y") - @testing.numpy_cupy_allclose(atol=1e-5) + @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64()) def test_interp_inf_to_nan(self, xp, dtype_y, dtype_x): # from NumPy's test_non_finite_inf x = xp.asarray([0.5], dtype=dtype_x) From b84dd7e065e13688152d353507e2cdc82ba236d1 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 14 Apr 2025 07:14:42 -0700 Subject: [PATCH 19/26] Add common_interpolate_checks with common utils --- dpnp/backend/extensions/ufunc/CMakeLists.txt | 1 + .../elementwise_functions/interpolate.cpp | 110 +++++++++++++----- 2 files changed, 80 insertions(+), 31 deletions(-) diff --git a/dpnp/backend/extensions/ufunc/CMakeLists.txt b/dpnp/backend/extensions/ufunc/CMakeLists.txt index 99163539191c..bbc6881ffcd0 100644 --- a/dpnp/backend/extensions/ufunc/CMakeLists.txt +++ b/dpnp/backend/extensions/ufunc/CMakeLists.txt @@ -70,6 +70,7 @@ endif() set_target_properties(${python_module_name} PROPERTIES CMAKE_POSITION_INDEPENDENT_CODE ON) target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../) +target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../common) target_include_directories(${python_module_name} PUBLIC ${Dpctl_INCLUDE_DIR}) target_include_directories(${python_module_name} PUBLIC ${Dpctl_TENSOR_INCLUDE_DIR}) diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp index 7e6cf7b4613e..784cef224548 100644 --- a/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp +++ b/dpnp/backend/extensions/ufunc/elementwise_functions/interpolate.cpp @@ -36,9 +36,15 @@ #include "kernels/elementwise_functions/interpolate.hpp" +#include "ext/validation_utils.hpp" + namespace py = pybind11; namespace td_ns = dpctl::tensor::type_dispatch; +using ext::validation::array_names; +using ext::validation::array_ptr; +using ext::validation::common_checks; + namespace dpnp::extensions::ufunc { @@ -105,27 +111,22 @@ sycl::event interpolate_call(sycl::queue &exec_q, interpolate_fn_ptr_t interpolate_dispatch_vector[td_ns::num_types]; -std::pair - py_interpolate(const dpctl::tensor::usm_ndarray &x, - const dpctl::tensor::usm_ndarray &idx, - const dpctl::tensor::usm_ndarray &xp, - const dpctl::tensor::usm_ndarray &fp, - std::optional &left, - std::optional &right, - dpctl::tensor::usm_ndarray &out, - sycl::queue &exec_q, - const std::vector &depends) +void common_interpolate_checks( + const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + const dpctl::tensor::usm_ndarray &out, + const std::optional &left, + const std::optional &right) { - int x_typenum = x.get_typenum(); - int xp_typenum = xp.get_typenum(); - int fp_typenum = fp.get_typenum(); - int out_typenum = out.get_typenum(); + array_names names = {{&x, "x"}, {&xp, "xp"}, {&fp, "fp"}, {&out, "out"}}; auto array_types = td_ns::usm_ndarray_types(); - int x_type_id = array_types.typenum_to_lookup_id(x_typenum); - int xp_type_id = array_types.typenum_to_lookup_id(xp_typenum); - int fp_type_id = array_types.typenum_to_lookup_id(fp_typenum); - int out_type_id = array_types.typenum_to_lookup_id(out_typenum); + int x_type_id = array_types.typenum_to_lookup_id(x.get_typenum()); + int xp_type_id = array_types.typenum_to_lookup_id(xp.get_typenum()); + int fp_type_id = array_types.typenum_to_lookup_id(fp.get_typenum()); + int out_type_id = array_types.typenum_to_lookup_id(out.get_typenum()); if (x_type_id != xp_type_id) { throw py::value_error("x and xp must have the same dtype"); @@ -134,17 +135,37 @@ std::pair throw py::value_error("fp and out must have the same dtype"); } - auto fn = interpolate_dispatch_vector[fp_type_id]; - if (!fn) { - throw py::type_error("Unsupported dtype"); + if (left) { + const auto &l = left.value(); + names.insert({&l, "left"}); + if (l.get_ndim() != 0) { + throw py::value_error("left must be a zero-dimensional array"); + } + + int left_type_id = array_types.typenum_to_lookup_id(l.get_typenum()); + if (left_type_id != fp_type_id) { + throw py::value_error( + "left must have the same dtype as fp and out"); + } } - if (!dpctl::utils::queues_are_compatible(exec_q, {x, idx, xp, fp, out})) { - throw py::value_error( - "Execution queue is not compatible with allocation queues"); + if (right) { + const auto &r = right.value(); + names.insert({&r, "right"}); + if (r.get_ndim() != 0) { + throw py::value_error("right must be a zero-dimensional array"); + } + + int right_type_id = array_types.typenum_to_lookup_id(r.get_typenum()); + if (right_type_id != fp_type_id) { + throw py::value_error( + "right must have the same dtype as fp and out"); + } } - dpctl::tensor::validation::CheckWritable::throw_if_not_writable(out); + common_checks({&x, &xp, &fp, left ? &left.value() : nullptr, + right ? &right.value() : nullptr}, + {&out}, names); if (x.get_ndim() != 1 || xp.get_ndim() != 1 || fp.get_ndim() != 1 || idx.get_ndim() != 1 || out.get_ndim() != 1) @@ -159,13 +180,40 @@ std::pair if (x.get_size() != out.get_size() || x.get_size() != idx.get_size()) { throw py::value_error("x, idx, and out must have the same size"); } +} + +std::pair + py_interpolate(const dpctl::tensor::usm_ndarray &x, + const dpctl::tensor::usm_ndarray &idx, + const dpctl::tensor::usm_ndarray &xp, + const dpctl::tensor::usm_ndarray &fp, + std::optional &left, + std::optional &right, + dpctl::tensor::usm_ndarray &out, + sycl::queue &exec_q, + const std::vector &depends) +{ + if (x.get_size() == 0) { + return {sycl::event(), sycl::event()}; + } + + common_interpolate_checks(x, idx, xp, fp, out, left, right); + + int out_typenum = out.get_typenum(); + + auto array_types = td_ns::usm_ndarray_types(); + int out_type_id = array_types.typenum_to_lookup_id(out_typenum); + + auto fn = interpolate_dispatch_vector[out_type_id]; + if (!fn) { + throw py::type_error("Unsupported dtype"); + } std::size_t n = x.get_size(); std::size_t xp_size = xp.get_size(); - void *left_ptr = left.has_value() ? left.value().get_data() : nullptr; - - void *right_ptr = right.has_value() ? right.value().get_data() : nullptr; + void *left_ptr = left ? left.value().get_data() : nullptr; + void *right_ptr = right ? right.value().get_data() : nullptr; sycl::event ev = fn(exec_q, x.get_data(), idx.get_data(), xp.get_data(), fp.get_data(), @@ -173,15 +221,15 @@ std::pair sycl::event args_ev; - if (left.has_value() && right.has_value()) { + if (left && right) { args_ev = dpctl::utils::keep_args_alive( exec_q, {x, idx, xp, fp, out, left.value(), right.value()}, {ev}); } - else if (left.has_value()) { + else if (left) { args_ev = dpctl::utils::keep_args_alive( exec_q, {x, idx, xp, fp, out, left.value()}, {ev}); } - else if (right.has_value()) { + else if (right) { args_ev = dpctl::utils::keep_args_alive( exec_q, {x, idx, xp, fp, out, right.value()}, {ev}); } From e9e357c9d33601288c4c92d52eaade297d38d065 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 14 Apr 2025 07:30:11 -0700 Subject: [PATCH 20/26] Reuse IsNan from common utils --- .../elementwise_functions/interpolate.hpp | 29 ++++--------------- 1 file changed, 5 insertions(+), 24 deletions(-) diff --git a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp index fd9e64d5a75a..cde69266f1b5 100644 --- a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -3,35 +3,15 @@ #include #include +#include "ext/common.hpp" #include "utils/type_utils.hpp" namespace type_utils = dpctl::tensor::type_utils; -namespace dpnp::kernels::interpolate -{ +using ext::common::IsNan; -template -struct IsNan +namespace dpnp::kernels::interpolate { - static bool isnan(const T &v) - { - if constexpr (type_utils::is_complex_v) { - using vT = typename T::value_type; - - const vT real1 = std::real(v); - const vT imag1 = std::imag(v); - - return IsNan::isnan(real1) || IsNan::isnan(imag1); - } - else if constexpr (std::is_floating_point_v || - std::is_same_v) { - return sycl::isnan(v); - } - - return false; - } -}; - template sycl::event interpolate_impl(sycl::queue &q, const TCoord *x, @@ -74,7 +54,8 @@ sycl::event interpolate_impl(sycl::queue &q, if (IsNan::isnan(res)) { res = slope * (x_val - xp[x_idx + 1]) + fp[x_idx + 1]; if (IsNan::isnan(res) && - (fp[x_idx] == fp[x_idx + 1])) { + (fp[x_idx] == fp[x_idx + 1])) + { res = fp[x_idx]; } } From 50e451365dd66f0e5be6f96939734451d8a51627 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 14 Apr 2025 07:45:13 -0700 Subject: [PATCH 21/26] Remove dublicate copy --- dpnp/dpnp_iface_mathematical.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 7615bb2cf768..c9fec899f981 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2882,9 +2882,6 @@ def interp(x, xp, fp, left=None, right=None, period=None): left = None right = None - x = x.astype(x_float_type) - xp = xp.astype(x_float_type) - # normalizing periodic boundaries x %= period xp %= period From dbeb3131c52a92a293f1606d502ad3b90a6b446a Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 14 Apr 2025 10:58:15 -0700 Subject: [PATCH 22/26] Add _validate_interp_param() function --- dpnp/dpnp_iface_mathematical.py | 50 +++++++++++++++++++++++---------- 1 file changed, 35 insertions(+), 15 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index c9fec899f981..0e19cfcb6bda 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -349,6 +349,36 @@ def _process_ediff1d_args(arg, arg_name, ary_dtype, ary_sycl_queue, usm_type): return arg, usm_type +def _validate_interp_param(param, name, exec_q, usm_type): + """ + Validate and convert optional parameters for interpolation. + + Returns a USM array or None if the input is None. + """ + if param is None: + return None + + if dpnp.is_supported_array_type(param): + if param.ndim != 0: + raise ValueError( + f"a {name} value must be 0-dimensional, " + f"but got {param.ndim}-dim" + ) + if dpu.get_execution_queue([exec_q, param.sycl_queue]) is None: + raise ValueError( + "input arrays and {name} must be on the same SYCL queue" + ) + return param.get_array() + + if dpnp.isscalar(param): + return dpt.asarray(param, sycl_queue=exec_q, usm_type=usm_type) + + raise TypeError( + f"a {name} value must be a scalar or 0-d supported array, " + f"but got {type(param)}" + ) + + _ABS_DOCSTRING = """ Calculates the absolute value for each element :math:`x_i` of input array `x`. @@ -2867,18 +2897,12 @@ def interp(x, xp, fp, left=None, right=None, period=None): fp = dpnp.asarray(fp, dtype=out_dtype, order="C") if period is not None: - if dpnp.is_supported_array_type(period): - if dpu.get_execution_queue([exec_q, period.sycl_queue]) is None: - raise ValueError( - "input arrays and period must be allocated " - "on the same SYCL queue" - ) - else: - period = dpnp.asarray(period, sycl_queue=exec_q, usm_type=usm_type) - + period = _validate_interp_param(period, "period", exec_q, usm_type) if period == 0: raise ValueError("period must be a non-zero value") period = dpnp.abs(period) + + # left/right are ignored when period is specified left = None right = None @@ -2895,12 +2919,8 @@ def interp(x, xp, fp, left=None, right=None, period=None): output = dpnp.empty(x.shape, dtype=out_dtype) idx = dpnp.searchsorted(xp, x, side="right") - left_usm = ( - dpnp.array(left, fp.dtype).get_array() if left is not None else None - ) - right_usm = ( - dpnp.array(right, fp.dtype).get_array() if right is not None else None - ) + left_usm = _validate_interp_param(left, "left", exec_q, usm_type) + right_usm = _validate_interp_param(right, "right", exec_q, usm_type) _manager = dpu.SequentialOrderManager[exec_q] mem_ev, ht_ev = ufi._interpolate( From dbb1b55597c26b358372070133db7ad0021b61d7 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 15 Apr 2025 02:50:32 -0700 Subject: [PATCH 23/26] Impove code coverage --- dpnp/dpnp_iface_mathematical.py | 12 ++++++++---- dpnp/tests/test_mathematical.py | 16 ++++++++++++++++ 2 files changed, 24 insertions(+), 4 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index 0e19cfcb6bda..c15876684922 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -349,7 +349,7 @@ def _process_ediff1d_args(arg, arg_name, ary_dtype, ary_sycl_queue, usm_type): return arg, usm_type -def _validate_interp_param(param, name, exec_q, usm_type): +def _validate_interp_param(param, name, exec_q, usm_type, dtype=None): """ Validate and convert optional parameters for interpolation. @@ -371,7 +371,9 @@ def _validate_interp_param(param, name, exec_q, usm_type): return param.get_array() if dpnp.isscalar(param): - return dpt.asarray(param, sycl_queue=exec_q, usm_type=usm_type) + return dpt.asarray( + param, dtype=dtype, sycl_queue=exec_q, usm_type=usm_type + ) raise TypeError( f"a {name} value must be a scalar or 0-d supported array, " @@ -2919,8 +2921,10 @@ def interp(x, xp, fp, left=None, right=None, period=None): output = dpnp.empty(x.shape, dtype=out_dtype) idx = dpnp.searchsorted(xp, x, side="right") - left_usm = _validate_interp_param(left, "left", exec_q, usm_type) - right_usm = _validate_interp_param(right, "right", exec_q, usm_type) + left_usm = _validate_interp_param(left, "left", exec_q, usm_type, fp.dtype) + right_usm = _validate_interp_param( + right, "right", exec_q, usm_type, fp.dtype + ) _manager = dpu.SequentialOrderManager[exec_q] mem_ev, ht_ev = ufi._interpolate( diff --git a/dpnp/tests/test_mathematical.py b/dpnp/tests/test_mathematical.py index a88fcc85e42c..f125187d6bf4 100644 --- a/dpnp/tests/test_mathematical.py +++ b/dpnp/tests/test_mathematical.py @@ -1270,6 +1270,9 @@ def test_errors(self): fp = dpnp.array([0.0, 1.0]) assert_raises(ValueError, dpnp.interp, x, xp, fp, period=0) + # period is not scalar or 0-dim + assert_raises(TypeError, dpnp.interp, x, xp, fp, period=[180]) + # period has a different SYCL queue q1 = dpctl.SyclQueue() q2 = dpctl.SyclQueue() @@ -1280,6 +1283,19 @@ def test_errors(self): period = dpnp.array([180], sycl_queue=q2) assert_raises(ValueError, dpnp.interp, x, xp, fp, period=period) + # left is not scalar or 0-dim + left = dpnp.array([1.0]) + assert_raises(ValueError, dpnp.interp, x, xp, fp, left=left) + + # left is 1-d array + left = dpnp.array([1.0]) + assert_raises(ValueError, dpnp.interp, x, xp, fp, left=left) + + # left has a different SYCL queue + left = dpnp.array(1.0, sycl_queue=q2) + if q1 != q2: + assert_raises(ValueError, dpnp.interp, x, xp, fp, left=left) + @pytest.mark.parametrize( "rhs", [[[1, 2, 3], [4, 5, 6]], [2.0, 1.5, 1.0], 3, 0.3] From cbe7e7a4bf5202394f6ab22d0ee3707d0375d4cd Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 15 Apr 2025 03:29:36 -0700 Subject: [PATCH 24/26] Add sycl_queue tests for interp --- dpnp/dpnp_iface_mathematical.py | 6 +++++- dpnp/tests/test_sycl_queue.py | 18 ++++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index c15876684922..a3b6c8d4e6c8 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -368,6 +368,8 @@ def _validate_interp_param(param, name, exec_q, usm_type, dtype=None): raise ValueError( "input arrays and {name} must be on the same SYCL queue" ) + if dtype is not None: + param = param.astype(dtype) return param.get_array() if dpnp.isscalar(param): @@ -2919,7 +2921,9 @@ def interp(x, xp, fp, left=None, right=None, period=None): assert xp.flags.c_contiguous assert fp.flags.c_contiguous - output = dpnp.empty(x.shape, dtype=out_dtype) + output = dpnp.empty( + x.shape, dtype=out_dtype, sycl_queue=exec_q, usm_type=usm_type + ) idx = dpnp.searchsorted(xp, x, side="right") left_usm = _validate_interp_param(left, "left", exec_q, usm_type, fp.dtype) right_usm = _validate_interp_param( diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py index b0112702e308..1f015d6ab2dd 100644 --- a/dpnp/tests/test_sycl_queue.py +++ b/dpnp/tests/test_sycl_queue.py @@ -1453,6 +1453,24 @@ def test_choose(device): assert_sycl_queue_equal(result.sycl_queue, chc.sycl_queue) +@pytest.mark.parametrize("device", valid_dev, ids=dev_ids) +@pytest.mark.parametrize("left", [None, dpnp.array(-1.0)]) +@pytest.mark.parametrize("right", [None, dpnp.array(99.0)]) +@pytest.mark.parametrize("period", [None, dpnp.array(180.0)]) +def test_interp(device, left, right, period): + x = dpnp.linspace(0.1, 9.9, 20, device=device) + xp = dpnp.linspace(0.0, 10.0, 5, sycl_queue=x.sycl_queue) + fp = dpnp.array(xp * 2 + 1, sycl_queue=x.sycl_queue) + + l = None if left is None else dpnp.array(left, sycl_queue=x.sycl_queue) + r = None if right is None else dpnp.array(right, sycl_queue=x.sycl_queue) + p = None if period is None else dpnp.array(period, sycl_queue=x.sycl_queue) + + result = dpnp.interp(x, xp, fp, left=l, right=r, period=p) + + assert_sycl_queue_equal(result.sycl_queue, x.sycl_queue) + + @pytest.mark.parametrize("device", valid_dev, ids=dev_ids) class TestLinAlgebra: @pytest.mark.parametrize( From aa102bd5ffd7e40f7cdb802ee9d906d799613f8e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 15 Apr 2025 03:50:34 -0700 Subject: [PATCH 25/26] Add usm_type tests for interp() --- dpnp/dpnp_iface_mathematical.py | 10 ++++-- dpnp/tests/test_usm_type.py | 59 +++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+), 3 deletions(-) diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py index a3b6c8d4e6c8..5c5eb7b05244 100644 --- a/dpnp/dpnp_iface_mathematical.py +++ b/dpnp/dpnp_iface_mathematical.py @@ -2921,15 +2921,19 @@ def interp(x, xp, fp, left=None, right=None, period=None): assert xp.flags.c_contiguous assert fp.flags.c_contiguous - output = dpnp.empty( - x.shape, dtype=out_dtype, sycl_queue=exec_q, usm_type=usm_type - ) idx = dpnp.searchsorted(xp, x, side="right") left_usm = _validate_interp_param(left, "left", exec_q, usm_type, fp.dtype) right_usm = _validate_interp_param( right, "right", exec_q, usm_type, fp.dtype ) + usm_type, exec_q = get_usm_allocations( + [x, xp, fp, period, left_usm, right_usm] + ) + output = dpnp.empty( + x.shape, dtype=out_dtype, sycl_queue=exec_q, usm_type=usm_type + ) + _manager = dpu.SequentialOrderManager[exec_q] mem_ev, ht_ev = ufi._interpolate( x.get_array(), diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py index 1d512ce111a6..ad8b6ba9403f 100644 --- a/dpnp/tests/test_usm_type.py +++ b/dpnp/tests/test_usm_type.py @@ -1268,6 +1268,65 @@ def test_choose(usm_type_x, usm_type_ind): assert z.usm_type == du.get_coerced_usm_type([usm_type_x, usm_type_ind]) +class TestInterp: + @pytest.mark.parametrize("usm_type_x", list_of_usm_types) + @pytest.mark.parametrize("usm_type_xp", list_of_usm_types) + @pytest.mark.parametrize("usm_type_fp", list_of_usm_types) + def test_basic(self, usm_type_x, usm_type_xp, usm_type_fp): + x = dpnp.linspace(0.1, 9.9, 20, usm_type=usm_type_x) + xp = dpnp.linspace(0.0, 10.0, 5, usm_type=usm_type_xp) + fp = dpnp.array(xp * 2 + 1, usm_type=usm_type_fp) + + result = dpnp.interp(x, xp, fp) + + assert x.usm_type == usm_type_x + assert xp.usm_type == usm_type_xp + assert fp.usm_type == usm_type_fp + assert result.usm_type == du.get_coerced_usm_type( + [usm_type_x, usm_type_xp, usm_type_fp] + ) + + @pytest.mark.parametrize("usm_type_x", list_of_usm_types) + @pytest.mark.parametrize("usm_type_left", list_of_usm_types) + @pytest.mark.parametrize("usm_type_right", list_of_usm_types) + def test_left_right(self, usm_type_x, usm_type_left, usm_type_right): + x = dpnp.linspace(-1.0, 11.0, 5, usm_type=usm_type_x) + xp = dpnp.linspace(0.0, 10.0, 5, usm_type=usm_type_x) + fp = dpnp.array(xp * 2 + 1, usm_type=usm_type_x) + + left = dpnp.array(-100, usm_type=usm_type_left) + right = dpnp.array(100, usm_type=usm_type_right) + + result = dpnp.interp(x, xp, fp, left=left, right=right) + + assert left.usm_type == usm_type_left + assert right.usm_type == usm_type_right + assert result.usm_type == du.get_coerced_usm_type( + [ + x.usm_type, + xp.usm_type, + fp.usm_type, + left.usm_type, + right.usm_type, + ] + ) + + @pytest.mark.parametrize("usm_type_x", list_of_usm_types) + @pytest.mark.parametrize("usm_type_period", list_of_usm_types) + def test_period(self, usm_type_x, usm_type_period): + x = dpnp.linspace(0.1, 9.9, 20, usm_type=usm_type_x) + xp = dpnp.linspace(0.0, 10.0, 5, usm_type=usm_type_x) + fp = dpnp.array(xp * 2 + 1, usm_type=usm_type_x) + period = dpnp.array(10.0, usm_type=usm_type_period) + + result = dpnp.interp(x, xp, fp, period=period) + + assert period.usm_type == usm_type_period + assert result.usm_type == du.get_coerced_usm_type( + [x.usm_type, xp.usm_type, fp.usm_type, period.usm_type] + ) + + @pytest.mark.parametrize("usm_type", list_of_usm_types) class TestLinAlgebra: @pytest.mark.parametrize( From 82c657e1efa2c957a66257b1b56c57474d6735d6 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 15 Apr 2025 05:33:50 -0700 Subject: [PATCH 26/26] Fix pre-commit remark --- dpnp/backend/kernels/elementwise_functions/interpolate.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp index cde69266f1b5..76ad0946aa22 100644 --- a/dpnp/backend/kernels/elementwise_functions/interpolate.hpp +++ b/dpnp/backend/kernels/elementwise_functions/interpolate.hpp @@ -54,8 +54,7 @@ sycl::event interpolate_impl(sycl::queue &q, if (IsNan::isnan(res)) { res = slope * (x_val - xp[x_idx + 1]) + fp[x_idx + 1]; if (IsNan::isnan(res) && - (fp[x_idx] == fp[x_idx + 1])) - { + (fp[x_idx] == fp[x_idx + 1])) { res = fp[x_idx]; } }