Add implementation of dpnp.float_power (#1957)

* Add implementation of dpnp.float_power()

* Updated third party tests

* Added more tests to cover different scenarios

* Mute back test for dpnp.floor_divide

* Update docstring with another way to state NaN value

* Address comment to update tests

Author: antonwolfy

doc/reference/ufunc.rst

@@ -20,21 +20,27 @@ Math operations
    dpnp.add
    dpnp.subtract
    dpnp.multiply
+   dpnp.matmul
    dpnp.divide
    dpnp.logaddexp
    dpnp.logaddexp2
    dpnp.true_divide
    dpnp.floor_divide
    dpnp.negative
+   dpnp.positive
    dpnp.power
+   dpnp.float_power
    dpnp.remainder
    dpnp.mod
    dpnp.fmod
-   dpnp.abs
+   dpnp.divmod
    dpnp.absolute
    dpnp.fabs
    dpnp.rint
    dpnp.sign
+   dpnp.heaviside
+   dpnp.conj
+   dpnp.conjugate
    dpnp.exp
    dpnp.exp2
    dpnp.log
@@ -44,13 +50,24 @@ Math operations
    dpnp.log1p
    dpnp.proj
    dpnp.sqrt
-   dpnp.cbrt
    dpnp.square
+   dpnp.cbrt
    dpnp.reciprocal
    dpnp.rsqrt
    dpnp.gcd
    dpnp.lcm
 
+.. tip::
+
+   The optional output arguments can be used to help you save memory
+   for large calculations. If your arrays are large, complicated
+   expressions can take longer than absolutely necessary due to the
+   creation and (later) destruction of temporary calculation
+   spaces. For example, the expression ``G = A * B + C`` is equivalent to
+   ``T1 = A * B; G = T1 + C; del T1``. It will be more quickly executed
+   as ``G = A * B; add(G, C, G)`` which is the same as
+   ``G = A * B; G += C``.
+
 
 Trigonometric functions
 ~~~~~~~~~~~~~~~~~~~~~~~

dpnp/backend/extensions/ufunc/CMakeLists.txt

@@ -27,6 +27,7 @@ set(_elementwise_sources
     ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/common.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/degrees.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/fabs.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/float_power.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/fmax.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/fmin.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/elementwise_functions/fmod.cpp

dpnp/backend/extensions/ufunc/elementwise_functions/common.cpp

@@ -27,6 +27,7 @@
 
 #include "degrees.hpp"
 #include "fabs.hpp"
+#include "float_power.hpp"
 #include "fmax.hpp"
 #include "fmin.hpp"
 #include "fmod.hpp"
@@ -44,6 +45,7 @@ void init_elementwise_functions(py::module_ m)
 {
     init_degrees(m);
     init_fabs(m);
+    init_float_power(m);
     init_fmax(m);
     init_fmin(m);
     init_fmod(m);

dpnp/backend/extensions/ufunc/elementwise_functions/float_power.cpp

@@ -0,0 +1,100 @@
+//*****************************************************************************
+// Copyright (c) 2024, Intel Corporation
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// maxification, 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 <sycl/sycl.hpp>
+
+#include "dpctl4pybind11.hpp"
+
+#include "float_power.hpp"
+
+// include a local copy of elementwise common header from dpctl tensor:
+// dpctl/tensor/libtensor/source/elementwise_functions/elementwise_functions.hpp
+// TODO: replace by including dpctl header once available
+#include "../../elementwise_functions/elementwise_functions.hpp"
+
+// dpctl tensor headers
+#include "utils/type_dispatch.hpp"
+
+namespace dpnp::extensions::ufunc
+{
+namespace py = pybind11;
+namespace py_int = dpnp::extensions::py_internal;
+
+namespace impl
+{
+namespace td_ns = dpctl::tensor::type_dispatch;
+
+// Supports only float and complex types
+template <typename T1, typename T2>
+struct OutputType
+{
+    using value_type = typename std::disjunction<
+        td_ns::BinaryTypeMapResultEntry<T1, float, T2, float, float>,
+        td_ns::BinaryTypeMapResultEntry<T1, double, T2, double, double>,
+        td_ns::BinaryTypeMapResultEntry<T1,
+                                        std::complex<float>,
+                                        T2,
+                                        std::complex<float>,
+                                        std::complex<float>>,
+        td_ns::BinaryTypeMapResultEntry<T1,
+                                        std::complex<double>,
+                                        T2,
+                                        std::complex<double>,
+                                        std::complex<double>>,
+        td_ns::DefaultResultEntry<void>>::result_type;
+};
+
+static int float_power_output_typeid_table[td_ns::num_types][td_ns::num_types];
+
+template <typename fnT, typename T1, typename T2>
+struct TypeMapFactory
+{
+    std::enable_if_t<std::is_same<fnT, int>::value, int> get()
+    {
+        using rT = typename OutputType<T1, T2>::value_type;
+        return td_ns::GetTypeid<rT>{}.get();
+    }
+};
+
+void populate_float_power_dispatch_tables(void)
+{
+    td_ns::DispatchTableBuilder<int, TypeMapFactory, td_ns::num_types> dvb;
+    dvb.populate_dispatch_table(float_power_output_typeid_table);
+}
+} // namespace impl
+
+void init_float_power(py::module_ m)
+{
+    impl::populate_float_power_dispatch_tables();
+    using impl::float_power_output_typeid_table;
+
+    auto float_power_result_type_pyapi = [&](const py::dtype &dtype1,
+                                             const py::dtype &dtype2) {
+        return py_int::py_binary_ufunc_result_type(
+            dtype1, dtype2, float_power_output_typeid_table);
+    };
+    m.def("_float_power_result_type", float_power_result_type_pyapi);
+}
+} // namespace dpnp::extensions::ufunc

dpnp/backend/extensions/ufunc/elementwise_functions/float_power.hpp

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

dpnp/dpnp_iface_mathematical.py

@@ -98,6 +98,7 @@
     "divide",
     "ediff1d",
     "fabs",
+    "float_power",
     "floor",
     "floor_divide",
     "fmax",
@@ -1364,6 +1365,105 @@ def ediff1d(x1, to_end=None, to_begin=None):
 )
 
 
+_FLOAT_POWER_DOCSTRING = """
+Calculates `x1_i` raised to `x2_i` for each element `x1_i` of the input array
+`x1` with the respective element `x2_i` of the input array `x2`.
+
+This differs from the power function in that boolean, integers, and float16 are
+promoted to floats with a minimum precision of float32 so that the result is
+always inexact. The intent is that the function will return a usable result for
+negative powers and seldom overflow for positive powers.
+
+Negative values raised to a non-integral value will return ``NaN``. To get
+complex results, cast the input to complex, or specify the ``dtype`` to be one
+of complex dtype.
+
+For full documentation refer to :obj:`numpy.float_power`.
+
+Parameters
+----------
+x1 : {dpnp.ndarray, usm_ndarray, scalar}
+    First input array, expected to have floating-point data types.
+    Both inputs `x1` and `x2` can not be scalars at the same time.
+x2 : {dpnp.ndarray, usm_ndarray, scalar}
+    Second input array, also expected to floating-point data types.
+    Both inputs `x1` and `x2` can not be scalars at the same time.
+out : {None, dpnp.ndarray, usm_ndarray}, optional
+    Output array to populate. Array must have the correct shape and
+    the expected data type.
+    Default: ``None``.
+order : {"C", "F", "A", "K"}, optional
+    Memory layout of the newly output array, if parameter `out` is ``None``.
+    Default: ``"K"``.
+
+Returns
+-------
+out : dpnp.ndarray
+    An array containing the bases in `x1` raised to the exponents in `x2`
+    element-wise.
+
+Limitations
+-----------
+Parameters `where` and `subok` are supported with their default values.
+Keyword argument `kwargs` is currently unsupported.
+Otherwise ``NotImplementedError`` exception will be raised.
+
+See Also
+--------
+:obj:`dpnp.power` : Power function that preserves type.
+
+Examples
+--------
+>>> import dpnp as np
+
+Cube each element in an array:
+
+>>> x1 = np.arange(6)
+>>> x1
+array([0, 1, 2, 3, 4, 5])
+>>> np.float_power(x1, 3)
+array([  0.,   1.,   8.,  27.,  64., 125.])
+
+Raise the bases to different exponents:
+
+>>> x2 = np.array([1.0, 2.0, 3.0, 3.0, 2.0, 1.0])
+>>> np.float_power(x1, x2)
+array([ 0.,  1.,  8., 27., 16.,  5.])
+
+The effect of broadcasting:
+
+>>> x2 = np.array([[1, 2, 3, 3, 2, 1], [1, 2, 3, 3, 2, 1]])
+>>> x2
+array([[1, 2, 3, 3, 2, 1],
+       [1, 2, 3, 3, 2, 1]])
+>>> np.float_power(x1, x2)
+array([[ 0.,  1.,  8., 27., 16.,  5.],
+       [ 0.,  1.,  8., 27., 16.,  5.]])
+
+Negative values raised to a non-integral value will result in ``NaN``:
+
+>>> x3 = np.array([-1, -4])
+>>> np.float_power(x3, 1.5)
+array([nan, nan])
+
+To get complex results, give the argument one of complex dtype, i.e.
+``dtype=np.complex64``:
+
+>>> np.float_power(x3, 1.5, dtype=np.complex64)
+array([1.1924881e-08-1.j, 9.5399045e-08-8.j], dtype=complex64)
+"""
+
+float_power = DPNPBinaryFunc(
+    "float_power",
+    ufi._float_power_result_type,
+    ti._pow,
+    _FLOAT_POWER_DOCSTRING,
+    mkl_fn_to_call=vmi._mkl_pow_to_call,
+    mkl_impl_fn=vmi._pow,
+    binary_inplace_fn=ti._pow_inplace,
+)
+
+
 _FLOOR_DOCSTRING = """
 Returns the floor for each element `x_i` for input array `x`.
 
@@ -2468,6 +2568,7 @@ def modf(x1, **kwargs):
 :obj:`dpnp.fmax` : Element-wise maximum of array elements.
 :obj:`dpnp.fmin` : Element-wise minimum of array elements.
 :obj:`dpnp.fmod` : Calculate the element-wise remainder of division.
+:obj:`dpnp.float_power` : Power function that promotes integers to floats.
 
 Examples
 --------

tests/skipped_tests.tbl

@@ -39,7 +39,6 @@ tests/test_umath.py::test_umaths[('divmod', 'ii')]
 tests/test_umath.py::test_umaths[('divmod', 'll')]
 tests/test_umath.py::test_umaths[('divmod', 'ff')]
 tests/test_umath.py::test_umaths[('divmod', 'dd')]
-tests/test_umath.py::test_umaths[('float_power', 'dd')]
 tests/test_umath.py::test_umaths[('frexp', 'f')]
 tests/test_umath.py::test_umaths[('frexp', 'd')]
 tests/test_umath.py::test_umaths[('gcd', 'ii')]

tests/skipped_tests_gpu.tbl

@@ -21,7 +21,6 @@ tests/test_umath.py::test_umaths[('divmod', 'ii')]
 tests/test_umath.py::test_umaths[('divmod', 'll')]
 tests/test_umath.py::test_umaths[('divmod', 'ff')]
 tests/test_umath.py::test_umaths[('divmod', 'dd')]
-tests/test_umath.py::test_umaths[('float_power', 'dd')]
 tests/test_umath.py::test_umaths[('floor_divide', 'ff')]
 tests/test_umath.py::test_umaths[('frexp', 'f')]
 tests/test_umath.py::test_umaths[('frexp', 'd')]

tests/test_sycl_queue.py

@@ -655,6 +655,11 @@ def test_reduce_hypot(device):
         pytest.param("dot", [3, 4, 5], [1, 2, 3]),
         pytest.param("dot", [3 + 2j, 4 + 1j, 5], [1, 2 + 3j, 3]),
         pytest.param("extract", [False, True, True, False], [0, 1, 2, 3]),
+        pytest.param(
+            "float_power",
+            [0, 1, 2, 3, 4, 5],
+            [1.0, 2.0, 3.0, 3.0, 2.0, 1.0],
+        ),
         pytest.param(
             "floor_divide", [1.0, 2.0, 3.0, 4.0], [2.5, 2.5, 2.5, 2.5]
         ),