diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 9dc3d0cc5cd3..a47715f4eb9b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -113,7 +113,7 @@ repos: "--disable=redefined-builtin", "--disable=unused-wildcard-import" ] - files: '^dpnp/(dpnp_iface.*|fft|linalg|special|dpnp_array)' + files: '^dpnp/(dpnp_iface.*|fft|linalg|scipy|dpnp_array)' - repo: https://github.com/macisamuele/language-formatters-pre-commit-hooks rev: v2.15.0 hooks: diff --git a/CHANGELOG.md b/CHANGELOG.md index 4fe8515d1346..e69c47fd5ed1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ This release is compatible with NumPy 2.3.3. * Added implementation of `dpnp.piecewise` [#2550](https://github.com/IntelPython/dpnp/pull/2550) * Added implementation of `dpnp.linalg.lu_solve` for 2D inputs (SciPy-compatible) [#2575](https://github.com/IntelPython/dpnp/pull/2575) * Added implementation of `dpnp.special.erfc` [#2588](https://github.com/IntelPython/dpnp/pull/2588) +* Added `dpnp.scipy` submodule to aggregate new SciPy-compatible functions from `linalg` and `special` namespaces [#2603](https://github.com/IntelPython/dpnp/pull/2603) ### Changed diff --git a/doc/reference/routines.rst b/doc/reference/routines.rst index 5cc64b74246a..1cfca82a15de 100644 --- a/doc/reference/routines.rst +++ b/doc/reference/routines.rst @@ -1,6 +1,5 @@ --------- -Routines --------- +Routines (NumPy) +================ The following pages describe NumPy-compatible routines. These functions cover a subset of diff --git a/doc/reference/scipy.rst b/doc/reference/scipy.rst new file mode 100644 index 000000000000..e95b199710a0 --- /dev/null +++ b/doc/reference/scipy.rst @@ -0,0 +1,12 @@ +Routines (SciPy) (:mod:`dpnp.scipy`) +==================================== + +The following pages describe SciPy-compatible routines. +These functions cover a subset of +`SciPy routines `_. + +.. toctree:: + :maxdepth: 2 + + scipy_linalg + scipy_special diff --git a/doc/reference/scipy_linalg.rst b/doc/reference/scipy_linalg.rst new file mode 100644 index 000000000000..0ab8d5a248ff --- /dev/null +++ b/doc/reference/scipy_linalg.rst @@ -0,0 +1,17 @@ +.. currentmodule:: dpnp.scipy.linalg + +Linear algebra (:mod:`dpnp.scipy.linalg`) +========================================= + +.. Hint:: `SciPy API Reference: Linear algebra (scipy.linalg) `_ + +Decompositions +-------------- + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + lu + lu_factor + lu_solve diff --git a/doc/reference/special.rst b/doc/reference/scipy_special.rst similarity index 80% rename from doc/reference/special.rst rename to doc/reference/scipy_special.rst index 5e6735a6afbe..0df23dbd6ce2 100644 --- a/doc/reference/special.rst +++ b/doc/reference/scipy_special.rst @@ -1,6 +1,6 @@ -.. currentmodule:: dpnp.special +.. currentmodule:: dpnp.scipy.special -Special functions (:mod:`dpnp.special`) +Special functions (:mod:`dpnp.scipy.special`) ======================================= .. Hint:: `SciPy API Reference: Special functions (scipy.special) `_ diff --git a/dpnp/__init__.py b/dpnp/__init__.py index 4fbac978ec99..3d8145b5362c 100644 --- a/dpnp/__init__.py +++ b/dpnp/__init__.py @@ -71,12 +71,13 @@ from .dpnp_iface_utils import __all__ as _ifaceutils__all__ from ._version import get_versions from . import linalg as linalg +from . import scipy as scipy __all__ = _iface__all__ __all__ += _ifaceutils__all__ # add submodules -__all__ += ["linalg"] +__all__ += ["linalg", "scipy"] __version__ = get_versions()["version"] diff --git a/dpnp/dpnp_iface.py b/dpnp/dpnp_iface.py index f7f4f470be48..d0c7ad9ad1ee 100644 --- a/dpnp/dpnp_iface.py +++ b/dpnp/dpnp_iface.py @@ -54,7 +54,6 @@ from dpnp.fft import * from dpnp.memory import * from dpnp.random import * -from dpnp.special import * __all__ = [ "are_same_logical_tensors", diff --git a/dpnp/dpnp_iface_arraycreation.py b/dpnp/dpnp_iface_arraycreation.py index 0c256d8d58de..55eae24c3b90 100644 --- a/dpnp/dpnp_iface_arraycreation.py +++ b/dpnp/dpnp_iface_arraycreation.py @@ -38,6 +38,8 @@ """ +# pylint: disable=duplicate-code + import operator import dpctl.tensor as dpt diff --git a/dpnp/dpnp_iface_manipulation.py b/dpnp/dpnp_iface_manipulation.py index b994819fa390..2175b920233c 100644 --- a/dpnp/dpnp_iface_manipulation.py +++ b/dpnp/dpnp_iface_manipulation.py @@ -38,6 +38,8 @@ """ +# pylint: disable=duplicate-code + import math import operator import warnings diff --git a/dpnp/linalg/__init__.py b/dpnp/linalg/__init__.py index 3139f5dfb84e..e3f0bdc6fe25 100644 --- a/dpnp/linalg/__init__.py +++ b/dpnp/linalg/__init__.py @@ -50,8 +50,6 @@ eigvalsh, inv, lstsq, - lu_factor, - lu_solve, matmul, matrix_norm, matrix_power, diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index c31a56c24957..ec48ccf381f6 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -58,8 +58,6 @@ dpnp_eigh, dpnp_inv, dpnp_lstsq, - dpnp_lu_factor, - dpnp_lu_solve, dpnp_matrix_power, dpnp_matrix_rank, dpnp_multi_dot, @@ -84,8 +82,6 @@ "eigvalsh", "inv", "lstsq", - "lu_factor", - "lu_solve", "matmul", "matrix_norm", "matrix_power", @@ -911,149 +907,6 @@ def lstsq(a, b, rcond=None): return dpnp_lstsq(a, b, rcond=rcond) -def lu_factor(a, overwrite_a=False, check_finite=True): - """ - Compute the pivoted LU decomposition of `a` matrix. - - The decomposition is:: - - A = P @ L @ U - - where `P` is a permutation matrix, `L` is lower triangular with unit - diagonal elements, and `U` is upper triangular. - - For full documentation refer to :obj:`scipy.linalg.lu_factor`. - - Parameters - ---------- - a : (..., M, N) {dpnp.ndarray, usm_ndarray} - Input array to decompose. - overwrite_a : {None, bool}, optional - Whether to overwrite data in `a` (may increase performance). - - Default: ``False``. - check_finite : {None, bool}, optional - Whether to check that the input matrix contains only finite numbers. - Disabling may give a performance gain, but may result in problems - (crashes, non-termination) if the inputs do contain infinities or NaNs. - - Default: ``True``. - - Returns - ------- - lu : (..., M, N) dpnp.ndarray - Matrix containing `U` in its upper triangle, - and `L` in its lower triangle. - The unit diagonal elements of `L` are not stored. - piv : (..., K) dpnp.ndarray - Pivot indices representing the permutation matrix `P`: - row i of matrix was interchanged with row piv[i]. - Where ``K = min(M, N)``. - - Warning - ------- - This function synchronizes in order to validate array elements - when ``check_finite=True``. - - See Also - -------- - :obj:`dpnp.linalg.lu_solve` : Solve an equation system using - the LU factorization of `a` matrix. - - Examples - -------- - >>> import dpnp as np - >>> a = np.array([[4., 3.], [6., 3.]]) - >>> lu, piv = np.linalg.lu_factor(a) - >>> lu - array([[6. , 3. ], - [0.66666667, 1. ]]) - >>> piv - array([1, 1]) - - """ - - dpnp.check_supported_arrays_type(a) - assert_stacked_2d(a) - - return dpnp_lu_factor(a, overwrite_a=overwrite_a, check_finite=check_finite) - - -def lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True): - """ - Solve a linear system, :math:`a x = b`, given the LU factorization of `a`. - - For full documentation refer to :obj:`scipy.linalg.lu_solve`. - - Parameters - ---------- - lu, piv : {tuple of dpnp.ndarrays or usm_ndarrays} - LU factorization of matrix `a` (M, M) together with pivot indices. - b : {(M,), (..., M, K)} {dpnp.ndarray, usm_ndarray} - Right-hand side - trans : {0, 1, 2} , optional - Type of system to solve: - - ===== ================= - trans system - ===== ================= - 0 :math:`a x = b` - 1 :math:`a^T x = b` - 2 :math:`a^H x = b` - ===== ================= - - Default: ``0``. - overwrite_b : {None, bool}, optional - Whether to overwrite data in `b` (may increase performance). - - Default: ``False``. - check_finite : {None, bool}, optional - Whether to check that the input matrix contains only finite numbers. - Disabling may give a performance gain, but may result in problems - (crashes, non-termination) if the inputs do contain infinities or NaNs. - - Default: ``True``. - - Returns - ------- - x : {(M,), (M, K)} dpnp.ndarray - Solution to the system - - Warning - ------- - This function synchronizes in order to validate array elements - when ``check_finite=True``. - - See Also - -------- - :obj:`dpnp.linalg.lu_factor` : LU factorize a matrix. - - Examples - -------- - >>> import dpnp as np - >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]]) - >>> b = np.array([1, 1, 1, 1]) - >>> lu, piv = np.linalg.lu_factor(A) - >>> x = np.linalg.lu_solve((lu, piv), b) - >>> np.allclose(A @ x - b, np.zeros((4,))) - array(True) - - """ - - (lu, piv) = lu_and_piv - dpnp.check_supported_arrays_type(lu, piv, b) - assert_stacked_2d(lu) - - return dpnp_lu_solve( - lu, - piv, - b, - trans=trans, - overwrite_b=overwrite_b, - check_finite=check_finite, - ) - - def matmul(x1, x2, /): """ Computes the matrix product. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index ce36b39f86a7..ee7650d4600d 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -33,13 +33,13 @@ """ +# pylint: disable=duplicate-code # pylint: disable=invalid-name # pylint: disable=no-name-in-module # pylint: disable=protected-access # pylint: disable=useless-import-alias from typing import NamedTuple -from warnings import warn import dpctl.tensor._tensor_impl as ti import dpctl.utils as dpu @@ -360,131 +360,6 @@ def _batched_lu_factor(a, res_type): return (out_a, out_ipiv, out_dev_info) -def _batched_lu_factor_scipy(a, res_type): # pylint: disable=too-many-locals - """SciPy-compatible LU factorization for batched inputs.""" - - # TODO: Find out at which array sizes the best performance is obtained - # getrf_batch can be slow on large GPU arrays. - # Use getrf_batch only on CPU. - # On GPU fall back to calling getrf per 2D slice. - use_batch = a.sycl_device.has_aspect_cpu - - a_sycl_queue = a.sycl_queue - a_usm_type = a.usm_type - _manager = dpu.SequentialOrderManager[a_sycl_queue] - - m, n = a.shape[-2:] - k = min(m, n) - orig_shape = a.shape - batch_shape = orig_shape[:-2] - - # handle empty input - if a.size == 0: - lu = dpnp.empty_like(a) - piv = dpnp.empty( - (*batch_shape, k), - dtype=dpnp.int64, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, - ) - return lu, piv - - # get 3d input arrays by reshape - a = dpnp.reshape(a, (-1, m, n)) - batch_size = a.shape[0] - - # Move batch axis to the end (m, n, batch) in Fortran order: - # required by getrf_batch - # and ensures each a[..., i] is F-contiguous for getrf - a = dpnp.moveaxis(a, 0, -1) - - a_usm_arr = dpnp.get_usm_ndarray(a) - - # `a` must be copied because getrf/getrf_batch destroys the input matrix - a_h = dpnp.empty_like(a, order="F", dtype=res_type) - ht_ev, copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=a_usm_arr, - dst=a_h.get_array(), - sycl_queue=a_sycl_queue, - depends=_manager.submitted_events, - ) - _manager.add_event_pair(ht_ev, copy_ev) - - ipiv_h = dpnp.empty( - (batch_size, k), - dtype=dpnp.int64, - order="C", - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, - ) - - if use_batch: - dev_info_h = [0] * batch_size - - ipiv_stride = k - a_stride = a_h.strides[-1] - - # Call the LAPACK extension function _getrf_batch - # to perform LU decomposition of a batch of general matrices - ht_ev, getrf_ev = li._getrf_batch( - a_sycl_queue, - a_h.get_array(), - ipiv_h.get_array(), - dev_info_h, - m, - n, - a_stride, - ipiv_stride, - batch_size, - depends=[copy_ev], - ) - _manager.add_event_pair(ht_ev, getrf_ev) - - if any(dev_info_h): - diag_nums = ", ".join(str(v) for v in dev_info_h if v > 0) - warn( - f"Diagonal numbers {diag_nums} are exactly zero. " - "Singular matrix.", - RuntimeWarning, - stacklevel=2, - ) - else: - dev_info_vecs = [[0] for _ in range(batch_size)] - - # Sequential LU factorization using getrf per slice - for i in range(batch_size): - ht_ev, getrf_ev = li._getrf( - a_sycl_queue, - a_h[..., i].get_array(), - ipiv_h[i].get_array(), - dev_info_vecs[i], - depends=[copy_ev], - ) - _manager.add_event_pair(ht_ev, getrf_ev) - - diag_nums = ", ".join( - str(v) for info in dev_info_vecs for v in info if v > 0 - ) - if diag_nums: - warn( - f"Diagonal number {diag_nums} are exactly zero. " - "Singular matrix.", - RuntimeWarning, - stacklevel=2, - ) - - # Restore original shape: move batch axis back and reshape - a_h = dpnp.moveaxis(a_h, -1, 0).reshape(orig_shape) - ipiv_h = ipiv_h.reshape((*batch_shape, k)) - - # oneMKL LAPACK uses 1-origin while SciPy uses 0-origin - ipiv_h -= 1 - - # Return a tuple containing the factorized matrix 'a_h', - # pivot indices 'ipiv_h' - return (a_h, ipiv_h) - - def _batched_solve(a, b, exec_q, res_usm_type, res_type): """ _batched_solve(a, b, exec_q, res_usm_type, res_type) @@ -1076,23 +951,6 @@ def _hermitian_svd(a, compute_uv): return dpnp.sort(s)[..., ::-1] -def _is_copy_required(a, res_type): - """ - Determine if `a` needs to be copied before LU decomposition. - This matches SciPy behavior: copy is needed unless input is suitable - for in-place modification. - """ - - if a.dtype != res_type: - return True - if not a.flags["F_CONTIGUOUS"]: - return True - if not a.flags["WRITABLE"]: - return True - - return False - - def _is_empty_2d(arr): # check size first for efficiency return arr.size == 0 and numpy.prod(arr.shape[-2:]) == 0 @@ -2382,215 +2240,6 @@ def dpnp_lstsq(a, b, rcond=None): return x, resids, rank, s -def dpnp_lu_factor(a, overwrite_a=False, check_finite=True): - """ - dpnp_lu_factor(a, overwrite_a=False, check_finite=True) - - Compute pivoted LU decomposition (SciPy-compatible behavior). - - This function mimics the behavior of `scipy.linalg.lu_factor` including - support for `overwrite_a`, `check_finite` and 0-based pivot indexing. - - """ - - res_type = _common_type(a) - a_sycl_queue = a.sycl_queue - a_usm_type = a.usm_type - - if check_finite: - if not dpnp.isfinite(a).all(): - raise ValueError("array must not contain infs or NaNs") - - if a.ndim > 2: - # SciPy always copies each 2D slice, - # so `overwrite_a` is ignored here - return _batched_lu_factor_scipy(a, res_type) - - # accommodate empty arrays - if a.size == 0: - lu = dpnp.empty_like(a) - piv = dpnp.arange( - 0, dtype=dpnp.int64, usm_type=a_usm_type, sycl_queue=a_sycl_queue - ) - return lu, piv - - _manager = dpu.SequentialOrderManager[a_sycl_queue] - a_usm_arr = dpnp.get_usm_ndarray(a) - - # SciPy-compatible behavior - # Copy is required if: - # - overwrite_a is False (always copy), - # - dtype mismatch, - # - not F-contiguous, - # - not writeable - if not overwrite_a or _is_copy_required(a, res_type): - a_h = dpnp.empty_like(a, order="F", dtype=res_type) - ht_ev, dep_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=a_usm_arr, - dst=a_h.get_array(), - sycl_queue=a_sycl_queue, - depends=_manager.submitted_events, - ) - _manager.add_event_pair(ht_ev, dep_ev) - dep_ev = [dep_ev] - else: - # input is suitable for in-place modification - a_h = a - dep_ev = _manager.submitted_events - - m, n = a.shape - - ipiv_h = dpnp.empty( - min(m, n), - dtype=dpnp.int64, - order="C", - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, - ) - dev_info_h = [0] - - # Call the LAPACK extension function _getrf - # to perform LU decomposition on the input matrix - ht_ev, getrf_ev = li._getrf( - a_sycl_queue, - a_h.get_array(), - ipiv_h.get_array(), - dev_info_h, - depends=dep_ev, - ) - _manager.add_event_pair(ht_ev, getrf_ev) - - if any(dev_info_h): - diag_nums = ", ".join(str(v) for v in dev_info_h if v > 0) - warn( - f"Diagonal number {diag_nums} is exactly zero. Singular matrix.", - RuntimeWarning, - stacklevel=2, - ) - - # MKL lapack uses 1-origin while SciPy uses 0-origin - ipiv_h -= 1 - - # Return a tuple containing the factorized matrix 'a_h', - # pivot indices 'ipiv_h' - return (a_h, ipiv_h) - - -def dpnp_lu_solve(lu, piv, b, trans=0, overwrite_b=False, check_finite=True): - """ - dpnp_lu_solve(lu, piv, b, trans=0, overwrite_b=False, check_finite=True) - - Solve an equation system (SciPy-compatible behavior). - - This function mimics the behavior of `scipy.linalg.lu_solve` including - support for `trans`, `overwrite_b`, `check_finite`, - and 0-based pivot indexing. - - """ - - res_usm_type, exec_q = get_usm_allocations([lu, piv, b]) - - res_type = _common_type(lu, b) - - # TODO: add broadcasting - if lu.shape[0] != b.shape[0]: - raise ValueError( - f"Shapes of lu {lu.shape} and b {b.shape} are incompatible" - ) - - if b.size == 0: - return dpnp.empty_like(b, dtype=res_type, usm_type=res_usm_type) - - if lu.ndim > 2: - raise NotImplementedError("Batched matrices are not supported") - - if check_finite: - if not dpnp.isfinite(lu).all(): - raise ValueError( - "LU factorization array must not contain infs or NaNs.\n" - "Note that when a singular matrix is given, unlike " - "dpnp.linalg.lu_factor returns an array containing NaN." - ) - if not dpnp.isfinite(b).all(): - raise ValueError( - "Right-hand side array must not contain infs or NaNs" - ) - - lu_usm_arr = dpnp.get_usm_ndarray(lu) - b_usm_arr = dpnp.get_usm_ndarray(b) - - # dpnp.linalg.lu_factor() returns 0-based pivots to match SciPy, - # convert to 1-based for oneMKL getrs - piv_h = piv + 1 - - _manager = dpu.SequentialOrderManager[exec_q] - dep_evs = _manager.submitted_events - - # oneMKL LAPACK getrs overwrites `lu`. - lu_h = dpnp.empty_like(lu, order="F", dtype=res_type, usm_type=res_usm_type) - - # use DPCTL tensor function to fill the сopy of the input array - # from the input array - ht_ev, lu_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=lu_usm_arr, - dst=lu_h.get_array(), - sycl_queue=lu.sycl_queue, - depends=dep_evs, - ) - _manager.add_event_pair(ht_ev, lu_copy_ev) - - # SciPy-compatible behavior - # Copy is required if: - # - overwrite_b is False (always copy), - # - dtype mismatch, - # - not F-contiguous, - # - not writeable - if not overwrite_b or _is_copy_required(b, res_type): - b_h = dpnp.empty_like( - b, order="F", dtype=res_type, usm_type=res_usm_type - ) - ht_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=b_usm_arr, - dst=b_h.get_array(), - sycl_queue=b.sycl_queue, - depends=dep_evs, - ) - _manager.add_event_pair(ht_ev, b_copy_ev) - dep_evs = [lu_copy_ev, b_copy_ev] - else: - # input is suitable for in-place modification - b_h = b - dep_evs = [lu_copy_ev] - - if not isinstance(trans, int): - raise TypeError("`trans` must be an integer") - - # Map SciPy-style trans codes (0, 1, 2) to MKL transpose enums - if trans == 0: - trans_mkl = li.Transpose.N - elif trans == 1: - trans_mkl = li.Transpose.T - elif trans == 2: - trans_mkl = li.Transpose.C - else: - raise ValueError("`trans` must be 0 (N), 1 (T), or 2 (C)") - - # Call the LAPACK extension function _getrs - # to solve the system of linear equations with an LU-factored - # coefficient square matrix, with multiple right-hand sides. - ht_ev, getrs_ev = li._getrs( - exec_q, - lu_h.get_array(), - piv_h.get_array(), - b_h.get_array(), - trans_mkl, - depends=dep_evs, - ) - _manager.add_event_pair(ht_ev, getrs_ev) - - return b_h - - def dpnp_matrix_power(a, n): """ dpnp_matrix_power(a, n) diff --git a/dpnp/scipy/__init__.py b/dpnp/scipy/__init__.py new file mode 100644 index 000000000000..989dc6eccd0e --- /dev/null +++ b/dpnp/scipy/__init__.py @@ -0,0 +1,39 @@ +# -*- coding: utf-8 -*- +# ***************************************************************************** +# 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. +# ***************************************************************************** + +""" +``dpnp.scipy`` +============== + +The SciPy-compatible interface of DPNP. + +This namespace provides submodules that mimic parts of ``SciPy`` on top of +DPNP functionality, reusing DPNP and oneMKL implementations underneath. +""" + +from . import linalg, special + +__all__ = ["linalg", "special"] diff --git a/dpnp/scipy/linalg/__init__.py b/dpnp/scipy/linalg/__init__.py new file mode 100644 index 000000000000..7109ee06bed4 --- /dev/null +++ b/dpnp/scipy/linalg/__init__.py @@ -0,0 +1,42 @@ +# -*- coding: utf-8 -*- +# ***************************************************************************** +# 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. +# ***************************************************************************** + +""" +``dpnp.scipy.linalg`` +===================== + +The SciPy-compatible linear algebra functions in DPNP rely on LAPACK +to provide efficient low-level implementations of standard algorithms. + +""" + + +from ._decomp_lu import lu_factor, lu_solve + +__all__ = [ + "lu_factor", + "lu_solve", +] diff --git a/dpnp/scipy/linalg/_decomp_lu.py b/dpnp/scipy/linalg/_decomp_lu.py new file mode 100644 index 000000000000..50a824a822cf --- /dev/null +++ b/dpnp/scipy/linalg/_decomp_lu.py @@ -0,0 +1,195 @@ +# -*- coding: utf-8 -*- +# ***************************************************************************** +# 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. +# ***************************************************************************** + + +""" +Interface of the SciPy-compatible Linear Algebra subset for DPNP. + +Notes +----- +This module exposes the public API for ``dpnp.scipy.linalg``. +It contains: + - SciPy-like interface functions + - documentation for the functions + +""" + + +import dpnp +from dpnp.linalg.dpnp_utils_linalg import assert_stacked_2d + +from ._utils import ( + dpnp_lu_factor, + dpnp_lu_solve, +) + +__all__ = [ + "lu_factor", + "lu_solve", +] + + +def lu_factor(a, overwrite_a=False, check_finite=True): + """ + Compute the pivoted LU decomposition of `a` matrix. + + The decomposition is:: + + A = P @ L @ U + + where `P` is a permutation matrix, `L` is lower triangular with unit + diagonal elements, and `U` is upper triangular. + + For full documentation refer to :obj:`scipy.linalg.lu_factor`. + + Parameters + ---------- + a : (..., M, N) {dpnp.ndarray, usm_ndarray} + Input array to decompose. + overwrite_a : {None, bool}, optional + Whether to overwrite data in `a` (may increase performance). + + Default: ``False``. + check_finite : {None, bool}, optional + Whether to check that the input matrix contains only finite numbers. + Disabling may give a performance gain, but may result in problems + (crashes, non-termination) if the inputs do contain infinities or NaNs. + + Default: ``True``. + + Returns + ------- + lu : (..., M, N) dpnp.ndarray + Matrix containing `U` in its upper triangle, + and `L` in its lower triangle. + The unit diagonal elements of `L` are not stored. + piv : (..., K) dpnp.ndarray + Pivot indices representing the permutation matrix `P`: + row i of matrix was interchanged with row piv[i]. + Where ``K = min(M, N)``. + + Warning + ------- + This function synchronizes in order to validate array elements + when ``check_finite=True``. + + See Also + -------- + :obj:`dpnp.scipy.linalg.lu_solve` : Solve an equation system using + the LU factorization of `a` matrix. + + Examples + -------- + >>> import dpnp as np + >>> a = np.array([[4., 3.], [6., 3.]]) + >>> lu, piv = np.scipy.linalg.lu_factor(a) + >>> lu + array([[6. , 3. ], + [0.66666667, 1. ]]) + >>> piv + array([1, 1]) + + """ + + dpnp.check_supported_arrays_type(a) + assert_stacked_2d(a) + + return dpnp_lu_factor(a, overwrite_a=overwrite_a, check_finite=check_finite) + + +def lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True): + """ + Solve a linear system, :math:`a x = b`, given the LU factorization of `a`. + + For full documentation refer to :obj:`scipy.linalg.lu_solve`. + + Parameters + ---------- + lu, piv : {tuple of dpnp.ndarrays or usm_ndarrays} + LU factorization of matrix `a` (M, M) together with pivot indices. + b : {(M,), (..., M, K)} {dpnp.ndarray, usm_ndarray} + Right-hand side + trans : {0, 1, 2} , optional + Type of system to solve: + + ===== ================= + trans system + ===== ================= + 0 :math:`a x = b` + 1 :math:`a^T x = b` + 2 :math:`a^H x = b` + ===== ================= + + Default: ``0``. + overwrite_b : {None, bool}, optional + Whether to overwrite data in `b` (may increase performance). + + Default: ``False``. + check_finite : {None, bool}, optional + Whether to check that the input matrix contains only finite numbers. + Disabling may give a performance gain, but may result in problems + (crashes, non-termination) if the inputs do contain infinities or NaNs. + + Default: ``True``. + + Returns + ------- + x : {(M,), (M, K)} dpnp.ndarray + Solution to the system + + Warning + ------- + This function synchronizes in order to validate array elements + when ``check_finite=True``. + + See Also + -------- + :obj:`dpnp.scipy.linalg.lu_factor` : LU factorize a matrix. + + Examples + -------- + >>> import dpnp as np + >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]]) + >>> b = np.array([1, 1, 1, 1]) + >>> lu, piv = np.scipy.linalg.lu_factor(A) + >>> x = np.scipy.linalg.lu_solve((lu, piv), b) + >>> np.allclose(A @ x - b, np.zeros((4,))) + array(True) + + """ + + (lu, piv) = lu_and_piv + dpnp.check_supported_arrays_type(lu, piv, b) + assert_stacked_2d(lu) + + return dpnp_lu_solve( + lu, + piv, + b, + trans=trans, + overwrite_b=overwrite_b, + check_finite=check_finite, + ) diff --git a/dpnp/scipy/linalg/_utils.py b/dpnp/scipy/linalg/_utils.py new file mode 100644 index 000000000000..377c41ac35c4 --- /dev/null +++ b/dpnp/scipy/linalg/_utils.py @@ -0,0 +1,406 @@ +# -*- coding: utf-8 -*- +# ***************************************************************************** +# 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. +# ***************************************************************************** + + +""" +Utility functions for the SciPy-compatible linear algebra interface. + +These include helper functions to check array properties and +functions with the main implementation part to fulfill the interface. +The main computational work is performed by enabling LAPACK functions +available as a pybind11 extension. + +""" + + +# pylint: disable=no-name-in-module +# pylint: disable=protected-access + +from warnings import warn + +import dpctl.tensor._tensor_impl as ti +import dpctl.utils as dpu + +import dpnp +import dpnp.backend.extensions.lapack._lapack_impl as li +from dpnp.dpnp_utils import get_usm_allocations +from dpnp.linalg.dpnp_utils_linalg import _common_type + +__all__ = [ + "dpnp_lu_factor", + "dpnp_lu_solve", +] + + +def _batched_lu_factor_scipy(a, res_type): # pylint: disable=too-many-locals + """SciPy-compatible LU factorization for batched inputs.""" + + # TODO: Find out at which array sizes the best performance is obtained + # getrf_batch can be slow on large GPU arrays. + # Use getrf_batch only on CPU. + # On GPU fall back to calling getrf per 2D slice. + use_batch = a.sycl_device.has_aspect_cpu + + a_sycl_queue = a.sycl_queue + a_usm_type = a.usm_type + _manager = dpu.SequentialOrderManager[a_sycl_queue] + + m, n = a.shape[-2:] + k = min(m, n) + orig_shape = a.shape + batch_shape = orig_shape[:-2] + + # handle empty input + if a.size == 0: + lu = dpnp.empty_like(a) + piv = dpnp.empty( + (*batch_shape, k), + dtype=dpnp.int64, + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + return lu, piv + + # get 3d input arrays by reshape + a = dpnp.reshape(a, (-1, m, n)) + batch_size = a.shape[0] + + # Move batch axis to the end (m, n, batch) in Fortran order: + # required by getrf_batch + # and ensures each a[..., i] is F-contiguous for getrf + a = dpnp.moveaxis(a, 0, -1) + + a_usm_arr = dpnp.get_usm_ndarray(a) + + # `a` must be copied because getrf/getrf_batch destroys the input matrix + a_h = dpnp.empty_like(a, order="F", dtype=res_type) + ht_ev, copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, + dst=a_h.get_array(), + sycl_queue=a_sycl_queue, + depends=_manager.submitted_events, + ) + _manager.add_event_pair(ht_ev, copy_ev) + + ipiv_h = dpnp.empty( + (batch_size, k), + dtype=dpnp.int64, + order="C", + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + + if use_batch: + dev_info_h = [0] * batch_size + + ipiv_stride = k + a_stride = a_h.strides[-1] + + # Call the LAPACK extension function _getrf_batch + # to perform LU decomposition of a batch of general matrices + ht_ev, getrf_ev = li._getrf_batch( + a_sycl_queue, + a_h.get_array(), + ipiv_h.get_array(), + dev_info_h, + m, + n, + a_stride, + ipiv_stride, + batch_size, + depends=[copy_ev], + ) + _manager.add_event_pair(ht_ev, getrf_ev) + + if any(dev_info_h): + diag_nums = ", ".join(str(v) for v in dev_info_h if v > 0) + warn( + f"Diagonal numbers {diag_nums} are exactly zero. " + "Singular matrix.", + RuntimeWarning, + stacklevel=2, + ) + else: + dev_info_vecs = [[0] for _ in range(batch_size)] + + # Sequential LU factorization using getrf per slice + for i in range(batch_size): + ht_ev, getrf_ev = li._getrf( + a_sycl_queue, + a_h[..., i].get_array(), + ipiv_h[i].get_array(), + dev_info_vecs[i], + depends=[copy_ev], + ) + _manager.add_event_pair(ht_ev, getrf_ev) + + diag_nums = ", ".join( + str(v) for info in dev_info_vecs for v in info if v > 0 + ) + if diag_nums: + warn( + f"Diagonal number {diag_nums} are exactly zero. " + "Singular matrix.", + RuntimeWarning, + stacklevel=2, + ) + + # Restore original shape: move batch axis back and reshape + a_h = dpnp.moveaxis(a_h, -1, 0).reshape(orig_shape) + ipiv_h = ipiv_h.reshape((*batch_shape, k)) + + # oneMKL LAPACK uses 1-origin while SciPy uses 0-origin + ipiv_h -= 1 + + # Return a tuple containing the factorized matrix 'a_h', + # pivot indices 'ipiv_h' + return (a_h, ipiv_h) + + +def _is_copy_required(a, res_type): + """ + Determine if `a` needs to be copied before LU decomposition. + This matches SciPy behavior: copy is needed unless input is suitable + for in-place modification. + """ + + if a.dtype != res_type: + return True + if not a.flags["F_CONTIGUOUS"]: + return True + if not a.flags["WRITABLE"]: + return True + + return False + + +def dpnp_lu_factor(a, overwrite_a=False, check_finite=True): + """ + dpnp_lu_factor(a, overwrite_a=False, check_finite=True) + + Compute pivoted LU decomposition (SciPy-compatible behavior). + + This function mimics the behavior of `scipy.linalg.lu_factor` including + support for `overwrite_a`, `check_finite` and 0-based pivot indexing. + + """ + + res_type = _common_type(a) + a_sycl_queue = a.sycl_queue + a_usm_type = a.usm_type + + if check_finite: + if not dpnp.isfinite(a).all(): + raise ValueError("array must not contain infs or NaNs") + + if a.ndim > 2: + # SciPy always copies each 2D slice, + # so `overwrite_a` is ignored here + return _batched_lu_factor_scipy(a, res_type) + + # accommodate empty arrays + if a.size == 0: + lu = dpnp.empty_like(a) + piv = dpnp.arange( + 0, dtype=dpnp.int64, usm_type=a_usm_type, sycl_queue=a_sycl_queue + ) + return lu, piv + + _manager = dpu.SequentialOrderManager[a_sycl_queue] + a_usm_arr = dpnp.get_usm_ndarray(a) + + # SciPy-compatible behavior + # Copy is required if: + # - overwrite_a is False (always copy), + # - dtype mismatch, + # - not F-contiguous, + # - not writeable + if not overwrite_a or _is_copy_required(a, res_type): + a_h = dpnp.empty_like(a, order="F", dtype=res_type) + ht_ev, dep_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, + dst=a_h.get_array(), + sycl_queue=a_sycl_queue, + depends=_manager.submitted_events, + ) + _manager.add_event_pair(ht_ev, dep_ev) + dep_ev = [dep_ev] + else: + # input is suitable for in-place modification + a_h = a + dep_ev = _manager.submitted_events + + m, n = a.shape + + ipiv_h = dpnp.empty( + min(m, n), + dtype=dpnp.int64, + order="C", + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + dev_info_h = [0] + + # Call the LAPACK extension function _getrf + # to perform LU decomposition on the input matrix + ht_ev, getrf_ev = li._getrf( + a_sycl_queue, + a_h.get_array(), + ipiv_h.get_array(), + dev_info_h, + depends=dep_ev, + ) + _manager.add_event_pair(ht_ev, getrf_ev) + + if any(dev_info_h): + diag_nums = ", ".join(str(v) for v in dev_info_h if v > 0) + warn( + f"Diagonal number {diag_nums} is exactly zero. Singular matrix.", + RuntimeWarning, + stacklevel=2, + ) + + # MKL lapack uses 1-origin while SciPy uses 0-origin + ipiv_h -= 1 + + # Return a tuple containing the factorized matrix 'a_h', + # pivot indices 'ipiv_h' + return (a_h, ipiv_h) + + +def dpnp_lu_solve(lu, piv, b, trans=0, overwrite_b=False, check_finite=True): + """ + dpnp_lu_solve(lu, piv, b, trans=0, overwrite_b=False, check_finite=True) + + Solve an equation system (SciPy-compatible behavior). + + This function mimics the behavior of `scipy.linalg.lu_solve` including + support for `trans`, `overwrite_b`, `check_finite`, + and 0-based pivot indexing. + + """ + + res_usm_type, exec_q = get_usm_allocations([lu, piv, b]) + + res_type = _common_type(lu, b) + + # TODO: add broadcasting + if lu.shape[0] != b.shape[0]: + raise ValueError( + f"Shapes of lu {lu.shape} and b {b.shape} are incompatible" + ) + + if b.size == 0: + return dpnp.empty_like(b, dtype=res_type, usm_type=res_usm_type) + + if lu.ndim > 2: + raise NotImplementedError("Batched matrices are not supported") + + if check_finite: + if not dpnp.isfinite(lu).all(): + raise ValueError( + "LU factorization array must not contain infs or NaNs.\n" + "Note that when a singular matrix is given, unlike " + "dpnp.scipy.linalg.lu_factor returns an array containing NaN." + ) + if not dpnp.isfinite(b).all(): + raise ValueError( + "Right-hand side array must not contain infs or NaNs" + ) + + lu_usm_arr = dpnp.get_usm_ndarray(lu) + b_usm_arr = dpnp.get_usm_ndarray(b) + + # dpnp.scipy.linalg.lu_factor() returns 0-based pivots to match SciPy, + # convert to 1-based for oneMKL getrs + piv_h = piv + 1 + + _manager = dpu.SequentialOrderManager[exec_q] + dep_evs = _manager.submitted_events + + # oneMKL LAPACK getrs overwrites `lu`. + lu_h = dpnp.empty_like(lu, order="F", dtype=res_type, usm_type=res_usm_type) + + # use DPCTL tensor function to fill the сopy of the input array + # from the input array + ht_ev, lu_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=lu_usm_arr, + dst=lu_h.get_array(), + sycl_queue=lu.sycl_queue, + depends=dep_evs, + ) + _manager.add_event_pair(ht_ev, lu_copy_ev) + + # SciPy-compatible behavior + # Copy is required if: + # - overwrite_b is False (always copy), + # - dtype mismatch, + # - not F-contiguous, + # - not writeable + if not overwrite_b or _is_copy_required(b, res_type): + b_h = dpnp.empty_like( + b, order="F", dtype=res_type, usm_type=res_usm_type + ) + ht_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=b_usm_arr, + dst=b_h.get_array(), + sycl_queue=b.sycl_queue, + depends=dep_evs, + ) + _manager.add_event_pair(ht_ev, b_copy_ev) + dep_evs = [lu_copy_ev, b_copy_ev] + else: + # input is suitable for in-place modification + b_h = b + dep_evs = [lu_copy_ev] + + if not isinstance(trans, int): + raise TypeError("`trans` must be an integer") + + # Map SciPy-style trans codes (0, 1, 2) to MKL transpose enums + if trans == 0: + trans_mkl = li.Transpose.N + elif trans == 1: + trans_mkl = li.Transpose.T + elif trans == 2: + trans_mkl = li.Transpose.C + else: + raise ValueError("`trans` must be 0 (N), 1 (T), or 2 (C)") + + # Call the LAPACK extension function _getrs + # to solve the system of linear equations with an LU-factored + # coefficient square matrix, with multiple right-hand sides. + ht_ev, getrs_ev = li._getrs( + exec_q, + lu_h.get_array(), + piv_h.get_array(), + b_h.get_array(), + trans_mkl, + depends=dep_evs, + ) + _manager.add_event_pair(ht_ev, getrs_ev) + + return b_h diff --git a/dpnp/special/__init__.py b/dpnp/scipy/special/__init__.py similarity index 97% rename from dpnp/special/__init__.py rename to dpnp/scipy/special/__init__.py index 8e1a3fcf24e0..50c8592a37d8 100644 --- a/dpnp/special/__init__.py +++ b/dpnp/scipy/special/__init__.py @@ -25,8 +25,8 @@ # ***************************************************************************** """ -``dpnp.special`` -================ +``dpnp.scipy.special`` +====================== The submodule provides a large collection of mathematical functions that are widely used in science and engineering. It includes special functions of diff --git a/dpnp/special/_erf.py b/dpnp/scipy/special/_erf.py similarity index 84% rename from dpnp/special/_erf.py rename to dpnp/scipy/special/_erf.py index 6060f09e943b..2b443d2628c2 100644 --- a/dpnp/special/_erf.py +++ b/dpnp/scipy/special/_erf.py @@ -29,7 +29,7 @@ Notes ----- -This module is a face or public interface file for the library +This module exposes the public interface for ``dpnp.scipy.special``. it contains: - Interface functions - documentation for the functions @@ -92,11 +92,11 @@ def __call__(self, x, out=None): # pylint: disable=signature-differs See Also -------- -:obj:`dpnp.special.erfc` : Complementary error function. -:obj:`dpnp.special.erfinv` : Inverse of the error function. -:obj:`dpnp.special.erfcinv` : Inverse of the complementary error function. -:obj:`dpnp.special.erfcx` : Scaled complementary error function. -:obj:`dpnp.special.erfi` : Imaginary error function. +:obj:`dpnp.scipy.special.erfc` : Complementary error function. +:obj:`dpnp.scipy.special.erfinv` : Inverse of the error function. +:obj:`dpnp.scipy.special.erfcinv` : Inverse of the complementary error function. +:obj:`dpnp.scipy.special.erfcx` : Scaled complementary error function. +:obj:`dpnp.scipy.special.erfi` : Imaginary error function. Notes ----- @@ -113,7 +113,7 @@ def __call__(self, x, out=None): # pylint: disable=signature-differs -------- >>> import dpnp as np >>> x = np.linspace(-3, 3, num=5) ->>> np.special.erf(x) +>>> np.scipy.special.erf(x) array([[-0.99997791, -0.96610515, 0. , 0.96610515, 0.99997791]) """ @@ -148,17 +148,17 @@ def __call__(self, x, out=None): # pylint: disable=signature-differs See Also -------- -:obj:`dpnp.special.erf` : Gauss error function. -:obj:`dpnp.special.erfinv` : Inverse of the error function. -:obj:`dpnp.special.erfcinv` : Inverse of the complementary error function. -:obj:`dpnp.special.erfcx` : Scaled complementary error function. -:obj:`dpnp.special.erfi` : Imaginary error function. +:obj:`dpnp.scipy.special.erf` : Gauss error function. +:obj:`dpnp.scipy.special.erfinv` : Inverse of the error function. +:obj:`dpnp.scipy.special.erfcinv` : Inverse of the complementary error function. +:obj:`dpnp.scipy.special.erfcx` : Scaled complementary error function. +:obj:`dpnp.scipy.special.erfi` : Imaginary error function. Examples -------- >>> import dpnp as np >>> x = np.linspace(-3, 3, num=5) ->>> np.special.erfc(x) +>>> np.scipy.special.erfc(x) array([[-0.99997791, -0.96610515, 0. , 0.96610515, 0.99997791]) """ diff --git a/dpnp/tests/test_linalg.py b/dpnp/tests/test_linalg.py index 0180ecc7d6de..7c0753b96fe6 100644 --- a/dpnp/tests/test_linalg.py +++ b/dpnp/tests/test_linalg.py @@ -1902,7 +1902,7 @@ def test_lu_factor(self, shape, order, dtype): a_np = self._make_nonsingular_np(shape, dtype, order) a_dp = dpnp.array(a_np, order=order) - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, check_finite=False, overwrite_a=False ) @@ -1926,7 +1926,7 @@ def test_lu_factor(self, shape, order, dtype): def test_overwrite_inplace(self, dtype): a_dp = dpnp.array([[4, 3], [6, 3]], dtype=dtype, order="F") a_dp_orig = a_dp.copy() - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, overwrite_a=True, check_finite=False ) @@ -1944,7 +1944,7 @@ def test_overwrite_inplace(self, dtype): def test_overwrite_copy(self, dtype): a_dp = dpnp.array([[4, 3], [6, 3]], dtype=dtype, order="C") a_dp_orig = a_dp.copy() - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, overwrite_a=True, check_finite=False ) @@ -1971,7 +1971,7 @@ def test_overwrite_copy_special(self): a2.flags["WRITABLE"] = False for a_dp, a_orig in zip((a1, a2), (a1_orig, a2_orig)): - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, overwrite_a=True, check_finite=False ) @@ -1989,7 +1989,7 @@ def test_overwrite_copy_special(self): @pytest.mark.parametrize("shape", [(0, 0), (0, 2), (2, 0)]) def test_empty_inputs(self, shape): a_dp = dpnp.empty(shape, dtype=dpnp.default_float_type(), order="F") - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) assert lu.shape == shape assert piv.shape == (min(shape),) @@ -2005,7 +2005,7 @@ def test_strided(self, sl): a_np = base[sl] a_dp = dpnp.array(a_np) - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) L, U = self._split_lu(lu, *a_dp.shape) PA = self._apply_pivots_rows(a_dp, piv) LU = L @ U @@ -2015,13 +2015,13 @@ def test_strided(self, sl): def test_singular_matrix(self): a_dp = dpnp.array([[1.0, 2.0], [2.0, 4.0]]) with pytest.warns(RuntimeWarning, match="Singular matrix"): - dpnp.linalg.lu_factor(a_dp, check_finite=False) + dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) @pytest.mark.parametrize("bad", [numpy.inf, -numpy.inf, numpy.nan]) def test_check_finite_raises(self, bad): a_dp = dpnp.array([[1.0, 2.0], [3.0, bad]], order="F") assert_raises( - ValueError, dpnp.linalg.lu_factor, a_dp, check_finite=True + ValueError, dpnp.scipy.linalg.lu_factor, a_dp, check_finite=True ) @@ -2077,7 +2077,7 @@ def test_lu_factor_batched(self, shape, order, dtype): a_np = self._make_nonsingular_nd_np(shape, dtype, order) a_dp = dpnp.array(a_np, order=order) - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, check_finite=False, overwrite_a=False ) @@ -2101,7 +2101,7 @@ def test_overwrite(self, dtype, order): a_np = self._make_nonsingular_nd_np((3, 2, 2), dtype, order) a_dp = dpnp.array(a_np, order=order) a_dp_orig = a_dp.copy() - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, overwrite_a=True, check_finite=False ) @@ -2124,7 +2124,7 @@ def test_overwrite(self, dtype, order): def test_empty_inputs(self, shape): a = dpnp.empty(shape, dtype=dpnp.default_float_type(), order="F") - lu, piv = dpnp.linalg.lu_factor(a, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a, check_finite=False) assert lu.shape == shape m, n = shape[-2:] assert piv.shape == (*shape[:-2], min(m, n)) @@ -2135,7 +2135,7 @@ def test_strided(self): ) a_dp = dpnp.array(a_np, order="F") a_stride = a_dp[::2] - lu, piv = dpnp.linalg.lu_factor(a_stride, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_stride, check_finite=False) for i in range(a_stride.shape[0]): L, U = self._split_lu(lu[i], 3, 3) PA = self._apply_pivots_rows( @@ -2149,12 +2149,14 @@ def test_singular_matrix(self): a[1] = dpnp.eye(2) a[2] = dpnp.array([[1.0, 1.0], [1.0, 1.0]]) with pytest.warns(RuntimeWarning, match="Singular matrix"): - dpnp.linalg.lu_factor(a, check_finite=False) + dpnp.scipy.linalg.lu_factor(a, check_finite=False) def test_check_finite_raises(self): a = dpnp.ones((2, 3, 3), dtype=dpnp.default_float_type(), order="F") a[1, 0, 0] = dpnp.nan - assert_raises(ValueError, dpnp.linalg.lu_factor, a, check_finite=True) + assert_raises( + ValueError, dpnp.scipy.linalg.lu_factor, a, check_finite=True + ) class TestLuSolve: @@ -2185,8 +2187,8 @@ def test_lu_solve(self, shape, rhs_cols, order, dtype): b_np = generate_random_numpy_array((n, rhs_cols), dtype, order) b_dp = dpnp.array(b_np, order=order) - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) - x = dpnp.linalg.lu_solve( + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) + x = dpnp.scipy.linalg.lu_solve( (lu, piv), b_dp, trans=0, overwrite_b=False, check_finite=False ) @@ -2202,8 +2204,8 @@ def test_trans(self, trans, dtype): a_dp = dpnp.array(a_np, order="F") b_dp = dpnp.array(generate_random_numpy_array((n, 2), dtype, "F")) - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) - x = dpnp.linalg.lu_solve( + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) + x = dpnp.scipy.linalg.lu_solve( (lu, piv), b_dp, trans=trans, overwrite_b=False, check_finite=False ) @@ -2222,10 +2224,10 @@ def test_overwrite_inplace(self, dtype): b_dp = dpnp.array([1, 0], dtype=dtype, order="F") b_orig = b_dp.copy() - lu, piv = dpnp.linalg.lu_factor( + lu, piv = dpnp.scipy.linalg.lu_factor( a_dp, overwrite_a=False, check_finite=False ) - x = dpnp.linalg.lu_solve( + x = dpnp.scipy.linalg.lu_solve( (lu, piv), b_dp, trans=0, overwrite_b=True, check_finite=False ) @@ -2236,11 +2238,11 @@ def test_overwrite_inplace(self, dtype): @pytest.mark.parametrize("dtype", get_float_complex_dtypes()) def test_overwrite_copy_special(self, dtype): a_dp = dpnp.array([[4, 3], [6, 3]], dtype=dtype, order="F") - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) # F-contig but dtype != res_type b1 = dpnp.array([1, 0], dtype=dpnp.int32, order="F") - x1 = dpnp.linalg.lu_solve( + x1 = dpnp.scipy.linalg.lu_solve( (lu, piv), b1, overwrite_b=True, check_finite=False ) assert x1 is not b1 @@ -2249,7 +2251,7 @@ def test_overwrite_copy_special(self, dtype): # F-contig, match dtype but read-only input b2 = dpnp.array([1, 0], dtype=dtype, order="F") b2.flags["WRITABLE"] = False - x2 = dpnp.linalg.lu_solve( + x2 = dpnp.scipy.linalg.lu_solve( (lu, piv), b2, overwrite_b=True, check_finite=False ) assert x2 is not b2 @@ -2276,8 +2278,8 @@ def test_diff_type(self, dtype_a, dtype_b): b_np = generate_random_numpy_array((3,), dtype_b, order="F") b_dp = dpnp.array(b_np, order="F") - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) - x = dpnp.linalg.lu_solve((lu, piv), b_dp, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) + x = dpnp.scipy.linalg.lu_solve((lu, piv), b_dp, check_finite=False) assert dpnp.allclose( a_dp @ x, b_dp.astype(x.dtype, copy=False), rtol=1e-5, atol=1e-5 ) @@ -2297,8 +2299,8 @@ def test_strided_rhs(self): ) b_dp = rhs_full[:, ::2][:, :3] - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) - x = dpnp.linalg.lu_solve( + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) + x = dpnp.scipy.linalg.lu_solve( (lu, piv), b_dp, overwrite_b=False, check_finite=False ) @@ -2324,8 +2326,8 @@ def test_broadcast_rhs(self, b_shape): b_np = generate_random_numpy_array(b_shape, dtype, order="F") b_dp = dpnp.array(b_np, order="F") - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) - x = dpnp.linalg.lu_solve( + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) + x = dpnp.scipy.linalg.lu_solve( (lu, piv), b_dp, overwrite_b=False, check_finite=False ) @@ -2348,20 +2350,20 @@ def test_empty_shapes(self, shape, rhs_cols): b_shape = (n, rhs_cols) b_dp = dpnp.empty(b_shape, dtype=dpnp.default_float_type(), order="F") - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) - x = dpnp.linalg.lu_solve((lu, piv), b_dp, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) + x = dpnp.scipy.linalg.lu_solve((lu, piv), b_dp, check_finite=False) assert x.shape == b_shape @pytest.mark.parametrize("bad", [numpy.inf, -numpy.inf, numpy.nan]) def test_check_finite_raises(self, bad): a_dp = dpnp.array([[1.0, 0.0], [0.0, 1.0]], order="F") - lu, piv = dpnp.linalg.lu_factor(a_dp, check_finite=False) + lu, piv = dpnp.scipy.linalg.lu_factor(a_dp, check_finite=False) b_bad = dpnp.array([1.0, bad], order="F") assert_raises( ValueError, - dpnp.linalg.lu_solve, + dpnp.scipy.linalg.lu_solve, (lu, piv), b_bad, check_finite=True, diff --git a/dpnp/tests/test_special.py b/dpnp/tests/test_special.py index e894c3cc46b5..e55fa001786d 100644 --- a/dpnp/tests/test_special.py +++ b/dpnp/tests/test_special.py @@ -25,7 +25,7 @@ def test_basic(self, func, dt): a = generate_random_numpy_array((2, 5), dtype=dt) ia = dpnp.array(a) - result = getattr(dpnp.special, func)(ia) + result = getattr(dpnp.scipy.special, func)(ia) expected = getattr(scipy.special, func)(a) # scipy >= 0.16.0 returns float64, but dpnp returns float32 @@ -41,7 +41,7 @@ def test_nan_inf(self, func): a = numpy.array([numpy.nan, -numpy.inf, numpy.inf]) ia = dpnp.array(a) - result = getattr(dpnp.special, func)(ia) + result = getattr(dpnp.scipy.special, func)(ia) expected = getattr(scipy.special, func)(a) assert_allclose(result, expected) @@ -51,7 +51,7 @@ def test_zeros(self, func): a = numpy.array([0.0, -0.0]) ia = dpnp.array(a) - result = getattr(dpnp.special, func)(ia) + result = getattr(dpnp.scipy.special, func)(ia) expected = getattr(scipy.special, func)(a) assert_allclose(result, expected) assert_equal(dpnp.signbit(result), numpy.signbit(expected)) @@ -60,7 +60,7 @@ def test_zeros(self, func): def test_complex(self, func, dt): x = dpnp.empty(5, dtype=dt) with pytest.raises(ValueError): - getattr(dpnp.special, func)(x) + getattr(dpnp.scipy.special, func)(x) class TestConsistency: @@ -73,11 +73,11 @@ def test_erfc(self): a = rng.pareto(0.02, n) * (2 * rng.randint(0, 2, n) - 1) a = dpnp.array(a) - res = 1 - dpnp.special.erf(a) + res = 1 - dpnp.scipy.special.erf(a) mask = dpnp.isfinite(res) a = a[mask] tol = 8 * dpnp.finfo(a).resolution assert dpnp.allclose( - dpnp.special.erfc(a), res[mask], rtol=tol, atol=tol + dpnp.scipy.special.erfc(a), res[mask], rtol=tol, atol=tol ) diff --git a/dpnp/tests/test_strides.py b/dpnp/tests/test_strides.py index 607934c472ef..3d91231bcb7b 100644 --- a/dpnp/tests/test_strides.py +++ b/dpnp/tests/test_strides.py @@ -175,7 +175,7 @@ def test_erf_funcs(func, stride): x = generate_random_numpy_array(10) a, ia = x[::stride], dpnp.array(x)[::stride] - result = getattr(dpnp.special, func)(ia) + result = getattr(dpnp.scipy.special, func)(ia) expected = getattr(scipy.special, func)(a) assert_dtype_allclose(result, expected) diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py index 9b2634fad7e1..bce599c22790 100644 --- a/dpnp/tests/test_sycl_queue.py +++ b/dpnp/tests/test_sycl_queue.py @@ -1493,7 +1493,7 @@ def test_interp(device, left, right, period): def test_erf_funcs(func, device): x = dpnp.linspace(-3, 3, num=5, device=device) - result = getattr(dpnp.special, func)(x) + result = getattr(dpnp.scipy.special, func)(x) assert_sycl_queue_equal(result.sycl_queue, x.sycl_queue) @@ -1603,7 +1603,7 @@ def test_lstsq(self, m, n, nrhs, device): ) def test_lu_factor(self, data, device): a = dpnp.array(data, device=device) - result = dpnp.linalg.lu_factor(a) + result = dpnp.scipy.linalg.lu_factor(a) for param in result: param_queue = param.sycl_queue @@ -1615,10 +1615,10 @@ def test_lu_factor(self, data, device): ) def test_lu_solve(self, b_data, device): a = dpnp.array([[1.0, 2.0], [3.0, 5.0]], device=device) - lu, piv = dpnp.linalg.lu_factor(a) + lu, piv = dpnp.scipy.linalg.lu_factor(a) b = dpnp.array(b_data, device=device) - result = dpnp.linalg.lu_solve((lu, piv), b) + result = dpnp.scipy.linalg.lu_solve((lu, piv), b) assert_sycl_queue_equal(result.sycl_queue, a.sycl_queue) assert_sycl_queue_equal(result.sycl_queue, b.sycl_queue) diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py index 90d2dcf38f41..eb059e335a3e 100644 --- a/dpnp/tests/test_usm_type.py +++ b/dpnp/tests/test_usm_type.py @@ -1300,7 +1300,7 @@ def test_choose(usm_type_x, usm_type_ind): @pytest.mark.parametrize("usm_type", list_of_usm_types) def test_erf_funcs(func, usm_type): x = dpnp.linspace(-3, 3, num=5, usm_type=usm_type) - y = getattr(dpnp.special, func)(x) + y = getattr(dpnp.scipy.special, func)(x) assert x.usm_type == y.usm_type == usm_type @@ -1480,7 +1480,7 @@ def test_lstsq(self, m, n, nrhs, usm_type, usm_type_other): ) def test_lu_factor(self, data, usm_type): a = dpnp.array(data, usm_type=usm_type) - result = dpnp.linalg.lu_factor(a) + result = dpnp.scipy.linalg.lu_factor(a) assert a.usm_type == usm_type for param in result: @@ -1493,10 +1493,10 @@ def test_lu_factor(self, data, usm_type): ) def test_lu_solve(self, b_data, usm_type, usm_type_rhs): a = dpnp.array([[1.0, 2.0], [3.0, 5.0]], usm_type=usm_type) - lu, piv = dpnp.linalg.lu_factor(a) + lu, piv = dpnp.scipy.linalg.lu_factor(a) b = dpnp.array(b_data, usm_type=usm_type_rhs) - result = dpnp.linalg.lu_solve((lu, piv), b) + result = dpnp.scipy.linalg.lu_solve((lu, piv), b) assert lu.usm_type == usm_type assert b.usm_type == usm_type_rhs diff --git a/dpnp/tests/third_party/cupy/testing/_loops.py b/dpnp/tests/third_party/cupy/testing/_loops.py index 946455857d7b..4e459ac4a02f 100644 --- a/dpnp/tests/third_party/cupy/testing/_loops.py +++ b/dpnp/tests/third_party/cupy/testing/_loops.py @@ -71,7 +71,7 @@ def _call_func_cupy(impl, args, kw, name, sp_name, scipy_name): if sp_name: kw[sp_name] = cupyx.scipy.sparse if scipy_name: - kw[scipy_name] = cupy + kw[scipy_name] = cupy.scipy kw[name] = cupy result, error = _call_func(impl, args, kw) return result, error diff --git a/dpnp/tests/third_party/cupyx/scipy_tests/linalg_tests/test_decomp_lu.py b/dpnp/tests/third_party/cupyx/scipy_tests/linalg_tests/test_decomp_lu.py index 2e0da0044138..fb51c3e39244 100644 --- a/dpnp/tests/third_party/cupyx/scipy_tests/linalg_tests/test_decomp_lu.py +++ b/dpnp/tests/third_party/cupyx/scipy_tests/linalg_tests/test_decomp_lu.py @@ -49,8 +49,7 @@ def test_lu_factor(self, dtype): a_cpu = testing.shaped_random(self.shape, numpy, dtype=dtype) a_gpu = cupy.asarray(a_cpu) result_cpu = scipy.linalg.lu_factor(a_cpu) - # Originally used cupyx.scipy.linalg.lu_factor - result_gpu = cupy.linalg.lu_factor(a_gpu) + result_gpu = cupy.scipy.linalg.lu_factor(a_gpu) assert len(result_cpu) == len(result_gpu) assert result_cpu[0].dtype == result_gpu[0].dtype # DPNP returns pivot indices as int64, while SciPy returns int32. @@ -63,7 +62,7 @@ def test_lu_factor(self, dtype): def check_lu_factor_reconstruction(self, A): m, n = self.shape - lu, piv = cupy.linalg.lu_factor(A) + lu, piv = cupy.scipy.linalg.lu_factor(A) # extract ``L`` and ``U`` from ``lu`` L = cupy.tril(lu, k=-1) cupy.fill_diagonal(L, 1.0) @@ -191,7 +190,7 @@ def test_lu_solve_backend(self, xp, dtype): lu = scipy.linalg.lu_factor(A) backend = "scipy" else: - lu = cupy.linalg.lu_factor(A) + lu = cupy.scipy.linalg.lu_factor(A) backend = cupy.linalg with scipy.linalg.set_backend(backend): out = scipy.linalg.lu_solve(lu, b, trans=self.trans) diff --git a/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py b/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py index c4a640a97d9e..02e505a51ed4 100644 --- a/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py +++ b/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py @@ -6,7 +6,7 @@ import pytest import dpnp as cupy -import dpnp.special +import dpnp.scipy.special from dpnp.tests.third_party.cupy import testing diff --git a/setup.py b/setup.py index 9c798f14b37b..5b8671ebcaa0 100644 --- a/setup.py +++ b/setup.py @@ -37,7 +37,9 @@ "dpnp.linalg", "dpnp.memory", "dpnp.random", - "dpnp.special", + "dpnp.scipy", + "dpnp.scipy.linalg", + "dpnp.scipy.special", ], package_data={ "dpnp": [