Merge pull request #20375 from tylerjereddy/treddy_prep_1_13_0_final

MAINT, REL: Prepare for SciPy 1.13.0 "final" (proposing to skip RC2 for Numpy 2 series support)

.github/workflows/wheels.yml

@@ -177,7 +177,7 @@ jobs:
           fi
 
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.16.5
+        uses: pypa/cibuildwheel@v2.17.0
         env:
           CIBW_BUILD: ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }}*
           CIBW_ARCHS: ${{ matrix.buildplat[2] }}

ci/cirrus_wheels.yml

@@ -1,6 +1,6 @@
 build_and_store_wheels: &BUILD_AND_STORE_WHEELS
   install_cibuildwheel_script:
-    - python -m pip install cibuildwheel==2.16.5
+    - python -m pip install cibuildwheel==2.17.0
   cibuildwheel_script:
     - cibuildwheel
   wheels_artifacts:

doc/source/reference/index.rst

@@ -84,12 +84,6 @@ section that the submodule in question is public. Of course you can still use::
 .. note:: SciPy is using a lazy loading mechanism which means that modules
           are only loaded in memory when you first try to access them.
 
-.. note::
-
-    The ``scipy`` namespace itself also contains functions imported from ``numpy``.
-    These functions still exist for backwards compatibility, but should be
-    imported from ``numpy`` directly.
-
 API definition
 --------------
 

meson.build

@@ -4,7 +4,7 @@ project(
   # Note that the git commit hash cannot be added dynamically here (it is added
   # in the dynamically generated and installed `scipy/version.py` though - see
   # tools/version_utils.py
-  version: '1.13.0rc2',
+  version: '1.13.0',
   license: 'BSD-3',
   meson_version: '>= 1.1.0',
   default_options: [

pyproject.toml

@@ -40,7 +40,7 @@ requires = [
 
 [project]
 name = "scipy"
-version = "1.13.0rc2"
+version = "1.13.0"
 # TODO: add `license-files` once PEP 639 is accepted (see meson-python#88)
 #       at that point, no longer include them in `py3.install_sources()`
 license = {file = "LICENSE.txt"}

requirements/build.txt

@@ -2,7 +2,7 @@
 # Do not edit this file; modify `pyproject.toml` instead and run `python tools/generate_requirements.py`.
 meson-python>=0.15.0
 Cython>=3.0.8
-pybind11>=2.10.4
+pybind11>=2.12.0
 pythran>=0.14.0
 ninja
 numpy>=1.22.4

scipy/_lib/tests/test_warnings.py

@@ -26,6 +26,7 @@ def visit_Attribute(self, node):
     def visit_Name(self, node):
         self.ls.append(node.id)
 
+
 class FindFuncs(ast.NodeVisitor):
     def __init__(self, filename):
         super().__init__()
@@ -96,6 +97,7 @@ def test_warning_calls_filters(warning_calls):
         os.path.join('datasets', '__init__.py'),
         os.path.join('optimize', '_optimize.py'),
         os.path.join('optimize', '_constraints.py'),
+        os.path.join('optimize', '_nnls.py'),
         os.path.join('signal', '_ltisys.py'),
         os.path.join('sparse', '__init__.py'),  # np.matrix pending-deprecation
         os.path.join('stats', '_discrete_distns.py'),  # gh-14901

scipy/interpolate/_rbfinterp.py

@@ -342,10 +342,11 @@ def __init__(self, y, d,
             degree = int(degree)
             if degree < -1:
                 raise ValueError("`degree` must be at least -1.")
-            elif degree < min_degree:
+            elif -1 < degree < min_degree:
                 warnings.warn(
-                    f"`degree` should not be below {min_degree} when `kernel` "
-                    f"is '{kernel}'. The interpolant may not be uniquely "
+                    f"`degree` should not be below {min_degree} except -1 "
+                    f"when `kernel` is '{kernel}'."
+                    f"The interpolant may not be uniquely "
                     f"solvable, and the smoothing parameter may have an "
                     f"unintuitive effect.",
                     UserWarning, stacklevel=2

scipy/interpolate/tests/test_rbfinterp.py

@@ -369,9 +369,18 @@ def test_degree_warning(self):
         y = np.linspace(0, 1, 5)[:, None]
         d = np.zeros(5)
         for kernel, deg in _NAME_TO_MIN_DEGREE.items():
-            match = f'`degree` should not be below {deg}'
-            with pytest.warns(Warning, match=match):
-                self.build(y, d, epsilon=1.0, kernel=kernel, degree=deg-1)
+            # Only test for kernels that its minimum degree is not 0.
+            if deg >= 1:
+                match = f'`degree` should not be below {deg}'
+                with pytest.warns(Warning, match=match):
+                    self.build(y, d, epsilon=1.0, kernel=kernel, degree=deg-1)
+
+    def test_minus_one_degree(self):
+        # Make sure a degree of -1 is accepted without any warning.
+        y = np.linspace(0, 1, 5)[:, None]
+        d = np.zeros(5)
+        for kernel, _ in _NAME_TO_MIN_DEGREE.items():
+            self.build(y, d, epsilon=1.0, kernel=kernel, degree=-1)
 
     def test_rank_error(self):
         # An error should be raised when `kernel` is "thin_plate_spline" and

scipy/linalg/_decomp_svd.py

@@ -116,6 +116,20 @@ def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False,
     if lapack_driver not in ('gesdd', 'gesvd'):
         message = f'lapack_driver must be "gesdd" or "gesvd", not "{lapack_driver}"'
         raise ValueError(message)
+
+    if lapack_driver == 'gesdd' and compute_uv:
+        # XXX: revisit int32 when ILP64 lapack becomes a thing
+        max_mn, min_mn = (m, n) if m > n else (n, m)
+        if full_matrices:
+            if max_mn*max_mn > numpy.iinfo(numpy.int32).max:
+                raise ValueError(f"Indexing a matrix size {max_mn} x {max_mn} "
+                                  " would incur integer overflow in LAPACK.")
+        else:
+            sz = max(m * min_mn, n * min_mn)
+            if max(m * min_mn, n * min_mn) > numpy.iinfo(numpy.int32).max:
+                raise ValueError(f"Indexing a matrix of {sz} elements would "
+                                  "incur an in integer overflow in LAPACK.")
+
     funcs = (lapack_driver, lapack_driver + '_lwork')
     gesXd, gesXd_lwork = get_lapack_funcs(funcs, (a1,), ilp64='preferred')
 

scipy/linalg/fblas_l1.pyf.src

@@ -388,7 +388,7 @@ function <prefix3>nrm2(n,x,offx,incx) result(n2)
 
   <ftype3> dimension(*),intent(in) :: x
 
-  integer optional, intent(in),check(incx>0||incx<0) :: incx = 1
+  integer optional, intent(in),check(incx>0) :: incx = 1
 
   integer optional,intent(in),depend(x) :: offx=0
   check(offx>=0 && offx<len(x)) :: offx
@@ -410,7 +410,7 @@ function <prefix4>nrm2(n,x,offx,incx) result(n2)
 
   <ftype4> dimension(*),intent(in) :: x
 
-  integer optional, intent(in),check(incx>0||incx<0) :: incx = 1
+  integer optional, intent(in),check(incx>0) :: incx = 1
 
   integer optional,intent(in),depend(x) :: offx=0
   check(offx>=0 && offx<len(x)) :: offx

scipy/linalg/tests/test_blas.py

@@ -1102,3 +1102,13 @@ def test_gh_169309():
     actual = scipy.linalg.blas.dnrm2(x, 5, 3, -1)
     expected = math.sqrt(500)
     assert_allclose(actual, expected)
+
+
+def test_dnrm2_neg_incx():
+    # check that dnrm2(..., incx < 0) raises
+    # XXX: remove the test after the lowest supported BLAS implements
+    # negative incx (new in LAPACK 3.10)
+    x = np.repeat(10, 9)
+    incx = -1
+    with assert_raises(fblas.__fblas_error):
+        scipy.linalg.blas.dnrm2(x, 5, 3, incx)

scipy/linalg/tests/test_decomp.py

@@ -59,22 +59,6 @@ def _random_hermitian_matrix(n, posdef=False, dtype=float):
 DTYPES = REAL_DTYPES + COMPLEX_DTYPES
 
 
-def clear_fuss(ar, fuss_binary_bits=7):
-    """Clears trailing `fuss_binary_bits` of mantissa of a floating number"""
-    x = np.asanyarray(ar)
-    if np.iscomplexobj(x):
-        return clear_fuss(x.real) + 1j * clear_fuss(x.imag)
-
-    significant_binary_bits = np.finfo(x.dtype).nmant
-    x_mant, x_exp = np.frexp(x)
-    f = 2.0**(significant_binary_bits - fuss_binary_bits)
-    x_mant *= f
-    np.rint(x_mant, out=x_mant)
-    x_mant /= f
-
-    return np.ldexp(x_mant, x_exp)
-
-
 # XXX: This function should not be defined here, but somewhere in
 #      scipy.linalg namespace
 def symrand(dim_or_eigv, rng):
@@ -238,11 +222,18 @@ def _check_gen_eig(self, A, B, atol_homog=1e-13, rtol_homog=1e-13):
                 assert_allclose(res[:, i], 0,
                                 rtol=1e-13, atol=1e-13, err_msg=msg)
 
+        # try to consistently order eigenvalues, including complex conjugate pairs
         w_fin = w[isfinite(w)]
         wt_fin = wt[isfinite(wt)]
-        perm = argsort(clear_fuss(w_fin))
-        permt = argsort(clear_fuss(wt_fin))
-        assert_allclose(w[perm], wt[permt],
+
+        # prune noise in the real parts
+        w_fin = -1j * np.real_if_close(1j*w_fin, tol=1e-10)
+        wt_fin = -1j * np.real_if_close(1j*wt_fin, tol=1e-10)
+
+        perm = argsort(w_fin)
+        permt = argsort(wt_fin)
+
+        assert_allclose(w_fin[perm], wt_fin[permt],
                         atol=1e-7, rtol=1e-7, err_msg=msg)
 
         length = np.empty(len(vr))
@@ -1099,6 +1090,14 @@ class TestSVD_GESVD(TestSVD_GESDD):
     lapack_driver = 'gesvd'
 
 
+def test_svd_gesdd_nofegfault():
+    # svd(a) with {U,VT}.size > INT_MAX does not segfault
+    # cf https://github.com/scipy/scipy/issues/14001
+    df=np.ones((4799, 53130), dtype=np.float64)
+    with assert_raises(ValueError):
+        svd(df)
+
+
 class TestSVDVals:
 
     def test_empty(self):

scipy/meson.build

@@ -85,7 +85,7 @@ npyrandom_path = _incdir_numpy_abs / '..' / '..' / 'random' / 'lib'
 npymath_lib = cc.find_library('npymath', dirs: npymath_path)
 npyrandom_lib = cc.find_library('npyrandom', dirs: npyrandom_path)
 
-pybind11_dep = dependency('pybind11', version: '>=2.10.4')
+pybind11_dep = dependency('pybind11', version: '>=2.12.0')
 
 # Pythran include directory and build flags
 if use_pythran

scipy/optimize/_linesearch.py

@@ -414,7 +414,7 @@ def extra_condition(alpha, phi):
             return True
 
     for i in range(maxiter):
-        if alpha1 == 0 or (amax is not None and alpha0 == amax):
+        if alpha1 == 0 or (amax is not None and alpha0 > amax):
             # alpha1 == 0: This shouldn't happen. Perhaps the increment has
             # slipped below machine precision?
             alpha_star = None

scipy/optimize/_nnls.py

@@ -1,6 +1,6 @@
 import numpy as np
-from scipy.linalg import solve
-
+from scipy.linalg import solve, LinAlgWarning
+import warnings
 
 __all__ = ['nnls']
 
@@ -112,43 +112,47 @@ def _nnls(A, b, maxiter=None, tol=None):
 
     # Initialize vars
     x = np.zeros(n, dtype=np.float64)
-
+    s = np.zeros(n, dtype=np.float64)
     # Inactive constraint switches
     P = np.zeros(n, dtype=bool)
 
     # Projected residual
-    resid = Atb.copy().astype(np.float64)  # x=0. Skip (-AtA @ x) term
+    w = Atb.copy().astype(np.float64)  # x=0. Skip (-AtA @ x) term
 
     # Overall iteration counter
     # Outer loop is not counted, inner iter is counted across outer spins
     iter = 0
 
-    while (not P.all()) and (resid[~P] > tol).any():  # B
+    while (not P.all()) and (w[~P] > tol).any():  # B
         # Get the "most" active coeff index and move to inactive set
-        resid[P] = -np.inf
-        k = np.argmax(resid)  # B.2
+        k = np.argmax(w * (~P))  # B.2
         P[k] = True  # B.3
 
         # Iteration solution
-        s = np.zeros(n, dtype=np.float64)
-        P_ind = P.nonzero()[0]
-        s[P] = solve(AtA[P_ind[:, None], P_ind[None, :]], Atb[P],
-                     assume_a='sym', check_finite=False)  # B.4
+        s[:] = 0.
+        # B.4
+        with warnings.catch_warnings():
+            warnings.filterwarnings('ignore', message='Ill-conditioned matrix',
+                                    category=LinAlgWarning)
+            s[P] = solve(AtA[np.ix_(P, P)], Atb[P], assume_a='sym', check_finite=False)
 
         # Inner loop
-        while (iter < maxiter) and (s[P].min() <= tol):  # C.1
-            alpha_ind = ((s < tol) & P).nonzero()
-            alpha = (x[alpha_ind] / (x[alpha_ind] - s[alpha_ind])).min()  # C.2
+        while (iter < maxiter) and (s[P].min() < 0):  # C.1
+            iter += 1
+            inds = P * (s < 0)
+            alpha = (x[inds] / (x[inds] - s[inds])).min()  # C.2
             x *= (1 - alpha)
             x += alpha*s
-            P[x < tol] = False
-            s[P] = solve(AtA[np.ix_(P, P)], Atb[P], assume_a='sym',
-                         check_finite=False)
+            P[x <= tol] = False
+            with warnings.catch_warnings():
+                warnings.filterwarnings('ignore', message='Ill-conditioned matrix',
+                                        category=LinAlgWarning)
+                s[P] = solve(AtA[np.ix_(P, P)], Atb[P], assume_a='sym',
+                             check_finite=False)
             s[~P] = 0  # C.6
-            iter += 1
 
         x[:] = s[:]
-        resid = Atb - AtA @ x
+        w[:] = Atb - AtA @ x
 
         if iter == maxiter:
             # Typically following line should return

scipy/optimize/_optimize.py

@@ -797,6 +797,15 @@ def _minimize_neldermead(func, x0, args=(), callback=None,
             maxfun = np.inf
 
     if bounds is not None:
+        # The default simplex construction may make all entries (for a given
+        # parameter) greater than an upper bound if x0 is very close to the
+        # upper bound. If one simply clips the simplex to the bounds this could
+        # make the simplex entries degenerate. If that occurs reflect into the
+        # interior.
+        msk = sim > upper_bound
+        # reflect into the interior
+        sim = np.where(msk, 2*upper_bound - sim, sim)
+        # but make sure the reflection is no less than the lower_bound
         sim = np.clip(sim, lower_bound, upper_bound)
 
     one2np1 = list(range(1, N + 1))
@@ -1153,7 +1162,7 @@ def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval,
 
 def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=np.inf,
               epsilon=_epsilon, maxiter=None, full_output=0, disp=1,
-              retall=0, callback=None, xrtol=0, c1=1e-4, c2=0.9, 
+              retall=0, callback=None, xrtol=0, c1=1e-4, c2=0.9,
               hess_inv0=None):
     """
     Minimize a function using the BFGS algorithm.
@@ -1226,7 +1235,7 @@ def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=np.inf,
     Optimize the function, `f`, whose gradient is given by `fprime`
     using the quasi-Newton method of Broyden, Fletcher, Goldfarb,
     and Shanno (BFGS).
-    
+
     Parameters `c1` and `c2` must satisfy ``0 < c1 < c2 < 1``.
 
     See Also
@@ -1298,7 +1307,7 @@ def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=np.inf,
 def _minimize_bfgs(fun, x0, args=(), jac=None, callback=None,
                    gtol=1e-5, norm=np.inf, eps=_epsilon, maxiter=None,
                    disp=False, return_all=False, finite_diff_rel_step=None,
-                   xrtol=0, c1=1e-4, c2=0.9, 
+                   xrtol=0, c1=1e-4, c2=0.9,
                    hess_inv0=None, **unknown_options):
     """
     Minimization of scalar function of one or more variables using the
@@ -2025,7 +2034,8 @@ def terminate(warnflag, msg):
     cg_maxiter = 20*len(x0)
 
     xtol = len(x0) * avextol
-    update_l1norm = 2 * xtol
+    # Make sure we enter the while loop.
+    update_l1norm = np.finfo(float).max
     xk = np.copy(x0)
     if retall:
         allvecs = [xk]
@@ -2116,7 +2126,7 @@ def terminate(warnflag, msg):
         update_l1norm = np.linalg.norm(update, ord=1)
 
     else:
-        if np.isnan(old_fval) or np.isnan(update).any():
+        if np.isnan(old_fval) or np.isnan(update_l1norm):
             return terminate(3, _status_message['nan'])
 
         msg = _status_message['success']