MAINT: utilise array_api_compat v1.8

[skip cirrus] [skip circle]

Author: lucascolley

scipy/_lib/_array_api.py

@@ -126,7 +126,11 @@ def array_namespace(*arrays: Array) -> ModuleType:
 
     _arrays = compliance_scipy(_arrays)
 
-    return array_api_compat.array_namespace(*_arrays)
+    # data-apis/array-api-compat#168
+    try: # return the wrapped namespace for NumPy arrays
+        return array_api_compat.array_namespace(*_arrays, use_compat=True)
+    except ValueError: # if the library is not wrapped, like array-api-strict
+        return array_api_compat.array_namespace(*_arrays, use_compat=None)
 
 
 def _asarray(
@@ -493,41 +497,6 @@ def scipy_namespace_for(xp: ModuleType) -> ModuleType | None:
     return None
 
 
-# temporary substitute for xp.minimum, which is not yet in all backends
-# or covered by array_api_compat.
-def xp_minimum(x1: Array, x2: Array, /) -> Array:
-    # xp won't be passed in because it doesn't need to be passed in to xp.minimum
-    xp = array_namespace(x1, x2)
-    if hasattr(xp, 'minimum'):
-        return xp.minimum(x1, x2)
-    x1, x2 = xp.broadcast_arrays(x1, x2)
-    i = (x2 < x1) | xp.isnan(x2)
-    res = xp.where(i, x2, x1)
-    return res[()] if res.ndim == 0 else res
-
-
-# temporary substitute for xp.clip, which is not yet in all backends
-# or covered by array_api_compat.
-def xp_clip(
-        x: Array,
-        /,
-        min: int | float | Array | None = None,
-        max: int | float | Array | None = None,
-        *,
-        xp: ModuleType | None = None) -> Array:
-    xp = array_namespace(x) if xp is None else xp
-    a, b = xp.asarray(min, dtype=x.dtype), xp.asarray(max, dtype=x.dtype)
-    if hasattr(xp, 'clip'):
-        return xp.clip(x, a, b)
-    x, a, b = xp.broadcast_arrays(x, a, b)
-    y = xp.asarray(x, copy=True)
-    ia = y < a
-    y[ia] = a[ia]
-    ib = y > b
-    y[ib] = b[ib]
-    return y[()] if y.ndim == 0 else y
-
-
 # temporary substitute for xp.moveaxis, which is not yet in all backends
 # or covered by array_api_compat.
 def xp_moveaxis_to_end(

scipy/_lib/array_api_compat

@@ -2,8 +2,7 @@
 import numpy as np
 import scipy._lib._elementwise_iterative_method as eim
 from scipy._lib._util import _RichResult
-from scipy._lib._array_api import (xp_clip, xp_minimum, xp_sign, copy as xp_copy,
-                                   xp_take_along_axis)
+from scipy._lib._array_api import xp_sign, copy as xp_copy, xp_take_along_axis
 
 # TODO:
 # - (maybe?) don't use fancy indexing assignment
@@ -142,7 +141,7 @@ def _chandrupatla(func, a, b, *, args=(), xatol=None, xrtol=None,
     xatol = 4*finfo.smallest_normal if xatol is None else xatol
     xrtol = 4*finfo.eps if xrtol is None else xrtol
     fatol = finfo.smallest_normal if fatol is None else fatol
-    frtol = frtol * xp_minimum(xp.abs(f1), xp.abs(f2))
+    frtol = frtol * xp.minimum(xp.abs(f1), xp.abs(f2))
     maxiter = (math.log2(finfo.max) - math.log2(finfo.smallest_normal)
                if maxiter is None else maxiter)
     work = _RichResult(x1=x1, f1=f1, x2=x2, f2=f2, x3=None, f3=None, t=0.5,
@@ -226,7 +225,7 @@ def post_termination_check(work):
         # [1] Figure 1 (last box; see also BASIC in appendix with comment
         # "Adjust T Away from the Interval Boundary")
         tl = 0.5 * work.tol / work.dx
-        work.t = xp_clip(t, tl, 1 - tl)
+        work.t = xp.clip(t, tl, 1 - tl)
 
     def customize_result(res, shape):
         xl, xr, fl, fr = res['xl'], res['xr'], res['fl'], res['fr']

scipy/optimize/_chandrupatla.py

@@ -6,7 +6,7 @@
 import scipy._lib._elementwise_iterative_method as eim
 from scipy.conftest import array_api_compatible
 from scipy._lib._array_api import (array_namespace, xp_assert_close, xp_assert_equal,
-                                   xp_assert_less, xp_minimum, is_numpy, is_cupy,
+                                   xp_assert_less, is_numpy, is_cupy,
                                    xp_ravel, size as xp_size)
 
 from scipy.optimize.elementwise import find_minimum, find_root
@@ -655,7 +655,7 @@ def f(*args, **kwargs):
         xp_assert_equal(res.fr, self.f(res.xr, *args_xp))
 
         assert xp.all(xp.abs(res.fun[finite]) ==
-                      xp_minimum(xp.abs(res.fl[finite]),
+                      xp.minimum(xp.abs(res.fl[finite]),
                                  xp.abs(res.fr[finite])))
 
     def test_flags(self, xp):
@@ -727,7 +727,7 @@ def test_convergence(self, xp):
         kwargs = kwargs0.copy()
         kwargs['frtol'] = 1e-3
         x1, x2 = bracket
-        f0 = xp_minimum(xp.abs(self.f(x1, *args)), xp.abs(self.f(x2, *args)))
+        f0 = xp.minimum(xp.abs(self.f(x1, *args)), xp.abs(self.f(x2, *args)))
         res1 = _chandrupatla_root(self.f, *bracket, **kwargs)
         xp_assert_less(xp.abs(res1.fun), 1e-3*f0)
         kwargs['frtol'] = 1e-6

scipy/optimize/tests/test_chandrupatla.py

@@ -15,10 +15,8 @@
 from scipy._lib._array_api import (
     array_namespace,
     size as xp_size,
-    xp_minimum,
     xp_moveaxis_to_end,
     xp_vector_norm,
-    xp_clip
 )
 
 from ._ansari_swilk_statistics import gscale, swilk
@@ -2899,7 +2897,7 @@ def bartlett(*samples, axis=0):
     chi2 = _SimpleChi2(xp.asarray(k-1))
     pvalue = _get_pvalue(T, chi2, alternative='greater', symmetric=False, xp=xp)
 
-    T = xp_clip(T, min=0., max=xp.inf)
+    T = xp.clip(T, min=0., max=xp.inf)
     T = T[()] if T.ndim == 0 else T
     pvalue = pvalue[()] if pvalue.ndim == 0 else pvalue
 
@@ -4103,7 +4101,7 @@ def circvar(samples, high=2*pi, low=0, axis=None, nan_policy='propagate'):
     cos_mean = xp.mean(cos_samp, axis=axis)
     hypotenuse = (sin_mean**2. + cos_mean**2.)**0.5
     # hypotenuse can go slightly above 1 due to rounding errors
-    R = xp_minimum(xp.asarray(1.), hypotenuse)
+    R = xp.clip(hypotenuse, max=1.)
 
     res = 1. - R
     return res
@@ -4207,7 +4205,7 @@ def circstd(samples, high=2*pi, low=0, axis=None, nan_policy='propagate', *,
     cos_mean = xp.mean(cos_samp, axis=axis)  # [1] (2.2.3)
     hypotenuse = (sin_mean**2. + cos_mean**2.)**0.5
     # hypotenuse can go slightly above 1 due to rounding errors
-    R = xp_minimum(xp.asarray(1.), hypotenuse)  # [1] (2.2.4)
+    R = xp.clip(hypotenuse, max=1.)  # [1] (2.2.4)
 
     res = (-2*xp.log(R))**0.5+0.0  # torch.pow returns -0.0 if R==1
     if not normalize:

scipy/stats/_morestats.py

@@ -8,8 +8,7 @@
 import inspect
 
 from scipy._lib._util import check_random_state, _rename_parameter, rng_integers
-from scipy._lib._array_api import (array_namespace, is_numpy, xp_minimum,
-                                   xp_clip, xp_moveaxis_to_end)
+from scipy._lib._array_api import array_namespace, is_numpy, xp_moveaxis_to_end
 from scipy.special import ndtr, ndtri, comb, factorial
 
 from ._common import ConfidenceInterval
@@ -996,15 +995,15 @@ def greater(null_distribution, observed):
     def two_sided(null_distribution, observed):
         pvalues_less = less(null_distribution, observed)
         pvalues_greater = greater(null_distribution, observed)
-        pvalues = xp_minimum(pvalues_less, pvalues_greater) * 2
+        pvalues = xp.minimum(pvalues_less, pvalues_greater) * 2
         return pvalues
 
     compare = {"less": less,
                "greater": greater,
                "two-sided": two_sided}
 
     pvalues = compare[alternative](null_distribution, observed)
-    pvalues = xp_clip(pvalues, 0., 1., xp=xp)
+    pvalues = xp.clip(pvalues, 0., 1.)
 
     return MonteCarloTestResult(observed, pvalues, null_distribution)
 

scipy/stats/_resampling.py

@@ -71,8 +71,7 @@
 from scipy.optimize import root_scalar
 from scipy._lib._util import normalize_axis_index
 from scipy._lib._array_api import (array_namespace, is_numpy, atleast_nd,
-                                   xp_clip, xp_moveaxis_to_end, xp_sign,
-                                   xp_minimum, xp_vector_norm)
+                                   xp_moveaxis_to_end, xp_sign, xp_vector_norm)
 from scipy._lib.array_api_compat import size as xp_size
 from scipy._lib.deprecation import _deprecated
 
@@ -1518,9 +1517,8 @@ def _get_pvalue(statistic, distribution, alternative, symmetric=True, xp=None):
         pvalue = distribution.sf(statistic)
     elif alternative == 'two-sided':
         pvalue = 2 * (distribution.sf(xp.abs(statistic)) if symmetric
-                      else xp_minimum(distribution.cdf(statistic),
-                                      distribution.sf(statistic),
-                                      xp=xp))
+                      else xp.minimum(distribution.cdf(statistic),
+                                      distribution.sf(statistic)))
     else:
         message = "`alternative` must be 'less', 'greater', or 'two-sided'."
         raise ValueError(message)
@@ -4625,9 +4623,7 @@ def statistic(x, y, axis):
     # Presumably, if abs(r) > 1, then it is only some small artifact of
     # floating point arithmetic.
     one = xp.asarray(1, dtype=dtype)
-    # `clip` only recently added to array API, so it's not yet available in
-    # array_api_strict. Replace with e.g. `xp.clip(r, -one, one)` when available.
-    r = xp.asarray(xp_clip(r, -one, one, xp=xp))
+    r = xp.asarray(xp.clip(r, -one, one))
     r[const_xy] = xp.nan
 
     # As explained in the docstring, the distribution of `r` under the null
@@ -7073,7 +7069,7 @@ def warn_masked(arg):
             f_obs_sum = _m_sum(f_obs_float, axis=axis, preserve_mask=False, xp=xp)
             f_exp_sum = _m_sum(f_exp, axis=axis, preserve_mask=False, xp=xp)
             relative_diff = (xp.abs(f_obs_sum - f_exp_sum) /
-                             xp_minimum(f_obs_sum, f_exp_sum))
+                             xp.minimum(f_obs_sum, f_exp_sum))
             diff_gt_tol = xp.any(relative_diff > rtol, axis=None)
         if diff_gt_tol:
             msg = (f"For each axis slice, the sum of the observed "