implement dpnp.logaddexp

Author: vtavana

dpnp/dpnp_algo/dpnp_elementwise_common.py

@@ -71,6 +71,7 @@
     "dpnp_less",
     "dpnp_less_equal",
     "dpnp_log",
+    "dpnp_logaddexp",
     "dpnp_logical_and",
     "dpnp_logical_not",
     "dpnp_logical_or",
@@ -1599,6 +1600,58 @@ def dpnp_log(x, out=None, order="K"):
     return dpnp_array._create_from_usm_ndarray(res_usm)
 
 
+_logaddexp_docstring_ = """
+logaddexp(x1, x2, out=None, order="K")
+
+Calculates the natural logarithm of the sum of exponentiations for each element
+`x1_i` of the input array `x1` with the respective element `x2_i` of the input
+array `x2`.
+
+This function calculates `log(exp(x1) + exp(x2))` more accurately for small
+values of `x`.
+
+Args:
+    x1 (dpnp.ndarray):
+        First input array, expected to have a real-valued floating-point
+        data type.
+    x2 (dpnp.ndarray):
+        Second input array, also expected to have a real-valued
+        floating-point data type.
+    out ({None, dpnp.ndarray}, optional):
+        Output array to populate.
+        Array have the correct shape and the expected data type.
+    order ("C","F","A","K", None, optional):
+        Memory layout of the newly output array, if parameter `out` is `None`.
+        Default: "K".
+Returns:
+    dpnp.ndarray:
+        An array containing the result of element-wise result. The data type
+        of the returned array is determined by the Type Promotion Rules.
+"""
+
+
+logaddexp_func = BinaryElementwiseFunc(
+    "logaddexp",
+    ti._logaddexp_result_type,
+    ti._logaddexp,
+    _logaddexp_docstring_,
+)
+
+
+def dpnp_logaddexp(x1, x2, out=None, order="K"):
+    """Invokes logaddexp() from dpctl.tensor implementation for logaddexp() function."""
+
+    # dpctl.tensor only works with usm_ndarray or scalar
+    x1_usm_or_scalar = dpnp.get_usm_ndarray_or_scalar(x1)
+    x2_usm_or_scalar = dpnp.get_usm_ndarray_or_scalar(x2)
+    out_usm = None if out is None else dpnp.get_usm_ndarray(out)
+
+    res_usm = logaddexp_func(
+        x1_usm_or_scalar, x2_usm_or_scalar, out=out_usm, order=order
+    )
+    return dpnp_array._create_from_usm_ndarray(res_usm)
+
+
 _logical_and_docstring_ = """
 logical_and(x1, x2, out=None, order='K')
 

dpnp/dpnp_iface_trigonometric.py

@@ -58,6 +58,7 @@
     dpnp_cos,
     dpnp_cosh,
     dpnp_log,
+    dpnp_logaddexp,
     dpnp_sin,
     dpnp_sinh,
     dpnp_sqrt,
@@ -87,6 +88,7 @@
     "log10",
     "log1p",
     "log2",
+    "logaddexp",
     "rad2deg",
     "radians",
     "reciprocal",
@@ -1089,6 +1091,69 @@ def log2(x1):
     return call_origin(numpy.log2, x1)
 
 
+def logaddexp(
+    x1,
+    x2,
+    /,
+    out=None,
+    *,
+    where=True,
+    order="K",
+    dtype=None,
+    subok=True,
+    **kwargs,
+):
+    """
+    Calculates ``log(exp(x1) + exp(x2))``, element-wise.
+
+    For full documentation refer to :obj:`numpy.logaddexp`.
+
+    Returns
+    -------
+    out : dpnp.ndarray
+        Logarithm of ``exp(x1) + exp(x2)``, element-wise.
+
+    Limitations
+    -----------
+    Parameters `x1` and `x2` are supported as either scalar, :class:`dpnp.ndarray`
+    or :class:`dpctl.tensor.usm_ndarray`, but both `x1` and `x2` can not be scalars at the same time.
+    Parameters `where`, `dtype` and `subok` are supported with their default values.
+    Keyword arguments `kwargs` are currently unsupported.
+    Otherwise the function will be executed sequentially on CPU.
+    Input array data types are limited by supported DPNP :ref:`Data types`.
+
+    See Also
+    --------
+    :obj:`dpnp.log` : Natural logarithm, element-wise.
+    :obj:`dpnp.exp` : Exponential, element-wise.
+
+    Examples
+    --------
+    >>> import dpnp as np
+    >>> prob1 = np.log(np.array(1e-50))
+    >>> prob2 = np.log(np.array(2.5e-50))
+    >>> prob12 = np.logaddexp(prob1, prob2)
+    >>> prob12
+    -113.87649168120691
+    >>> np.exp(prob12)
+    3.5000000000000057e-50
+
+    """
+
+    return check_nd_call_func(
+        numpy.logaddexp,
+        dpnp_logaddexp,
+        x1,
+        x2,
+        out=out,
+        where=where,
+        order=order,
+        dtype=dtype,
+        subok=subok,
+        **kwargs,
+    )
+
+
 def reciprocal(x1, **kwargs):
     """
     Return the reciprocal of the argument, element-wise.

tests/skipped_tests.tbl

@@ -113,8 +113,6 @@ tests/test_umath.py::test_umaths[('ldexp', 'fi')]
 tests/test_umath.py::test_umaths[('ldexp', 'fl')]
 tests/test_umath.py::test_umaths[('ldexp', 'di')]
 tests/test_umath.py::test_umaths[('ldexp', 'dl')]
-tests/test_umath.py::test_umaths[('logaddexp', 'ff')]
-tests/test_umath.py::test_umaths[('logaddexp', 'dd')]
 tests/test_umath.py::test_umaths[('logaddexp2', 'ff')]
 tests/test_umath.py::test_umaths[('logaddexp2', 'dd')]
 tests/test_umath.py::test_umaths[('nextafter', 'ff')]
@@ -513,7 +511,6 @@ tests/third_party/cupy/math_tests/test_arithmetic.py::TestArithmeticModf::test_m
 
 tests/third_party/cupy/math_tests/test_arithmetic.py::TestArithmeticRaisesWithNumpyInput_param_3_{name='angle', nargs=1}::test_raises_with_numpy_input
 
-tests/third_party/cupy/math_tests/test_explog.py::TestExplog::test_logaddexp
 tests/third_party/cupy/math_tests/test_explog.py::TestExplog::test_logaddexp2
 tests/third_party/cupy/math_tests/test_explog.py::TestExplog::test_logaddexp2_infinities
 tests/third_party/cupy/math_tests/test_floating.py::TestFloating::test_copysign_float

tests/skipped_tests_gpu.tbl

@@ -53,8 +53,6 @@ tests/test_umath.py::test_umaths[('ldexp', 'fi')]
 tests/test_umath.py::test_umaths[('ldexp', 'fl')]
 tests/test_umath.py::test_umaths[('ldexp', 'di')]
 tests/test_umath.py::test_umaths[('ldexp', 'dl')]
-tests/test_umath.py::test_umaths[('logaddexp', 'ff')]
-tests/test_umath.py::test_umaths[('logaddexp', 'dd')]
 tests/test_umath.py::test_umaths[('logaddexp2', 'ff')]
 tests/test_umath.py::test_umaths[('logaddexp2', 'dd')]
 tests/test_umath.py::test_umaths[('nextafter', 'ff')]
@@ -642,7 +640,6 @@ tests/third_party/cupy/manipulation_tests/test_tiling.py::TestTile_param_5_{reps
 
 tests/third_party/cupy/math_tests/test_arithmetic.py::TestArithmeticRaisesWithNumpyInput_param_3_{name='angle', nargs=1}::test_raises_with_numpy_input
 
-tests/third_party/cupy/math_tests/test_explog.py::TestExplog::test_logaddexp
 tests/third_party/cupy/math_tests/test_explog.py::TestExplog::test_logaddexp2
 tests/third_party/cupy/math_tests/test_explog.py::TestExplog::test_logaddexp2_infinities
 tests/third_party/cupy/math_tests/test_floating.py::TestFloating::test_copysign_float

tests/test_strides.py

@@ -138,6 +138,7 @@ def test_strides_reciprocal(dtype, shape):
         "arctan2",
         "divide",
         "hypot",
+        "logaddexp",
         "maximum",
         "minimum",
         "multiply",

tests/test_sycl_queue.py

@@ -353,6 +353,11 @@ def test_proj(device):
             [-3.0, -2.0, -1.0, 1.0, 2.0, 3.0],
             [2.0, 2.0, 2.0, 2.0, 2.0, 2.0],
         ),
+        pytest.param(
+            "logaddexp",
+            [[-1, 2, 5, 9]],
+            [[4, -3, 2, -8]],
+        ),
         pytest.param(
             "matmul", [[1.0, 0.0], [0.0, 1.0]], [[4.0, 1.0], [1.0, 2.0]]
         ),

tests/test_usm_type.py

@@ -337,6 +337,11 @@ def test_1in_1out(func, data, usm_type):
             [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0]],
             [[4.0, 4.0], [4.0, 4.0], [4.0, 4.0]],
         ),
+        pytest.param(
+            "logaddexp",
+            [[-1, 2, 5, 9]],
+            [[4, -3, 2, -8]],
+        ),
     ],
 )
 @pytest.mark.parametrize("usm_type_x", list_of_usm_types, ids=list_of_usm_types)

tests/third_party/cupy/math_tests/test_explog.py

@@ -3,6 +3,7 @@
 import numpy
 import pytest
 
+from tests.helper import has_support_aspect64
 from tests.third_party.cupy import testing
 
 
@@ -17,7 +18,7 @@ def check_unary(self, name, xp, dtype, no_complex=False):
         return getattr(xp, name)(a)
 
     @testing.for_all_dtypes()
-    @testing.numpy_cupy_allclose(atol=1e-5)
+    @testing.numpy_cupy_allclose(atol=1e-5, type_check=has_support_aspect64())
     def check_binary(self, name, xp, dtype, no_complex=False):
         if no_complex:
             if numpy.dtype(dtype).kind == "c":
@@ -62,3 +63,11 @@ def test_logaddexp2(self):
     def test_logaddexp2_infinities(self, xp, dtype, val):
         a = xp.full((2, 3), val, dtype=dtype)
         return xp.logaddexp2(a, a)
+
+
+@pytest.mark.parametrize("val", [numpy.inf, -numpy.inf])
+@testing.for_float_dtypes()
+@testing.numpy_cupy_allclose()
+def test_logaddexp_infinities(xp, dtype, val):
+    a = xp.full((2, 3), val, dtype=dtype)
+    return xp.logaddexp(a, a)