Merge pull request scipy#22233 from tylerjereddy/treddy_backports_1.15.0_final

* DEP: linalg.solve_toeplitz/matmul_toeplitz: warn on n-D `c`, `r` (scipy#22193)

* DEP: linalg.solve_toeplitz: warn on n-D c, r

* MAINT: linalg.matmul_toeplitz: add admonition to documentation; test warning

* DOC: update release notes

[docs only]

* DOC: differentiate.jacobian: correct/improve documentation about callable interface

[docs only]

* Update scipy/differentiate/_differentiate.py

[docs only]

* DOC: Update 1.15.0 relnotes

* Update the SciPy `1.15.0` release notes following
additional backport activity.

* MAINT: PR 22233 wheel build [wheel build]

* Test wheel builds in above PR before release
process starts.

[wheel build]

---------

Co-authored-by: Matt Haberland <[email protected]>

Author: tylerjereddy

doc/source/release/1.15.0-notes.rst

@@ -2,8 +2,6 @@
 SciPy 1.15.0 Release Notes
 ==========================
 
-.. note:: SciPy 1.15.0 is not released yet!
-
 .. contents::
 
 SciPy 1.15.0 is the culmination of 6 months of hard work. It contains
@@ -376,9 +374,9 @@ with support added for SciPy ``1.15.0`` include:
   and for other backends will transit via NumPy arrays on the host.
 
 
-*******************
-Deprecated features
-*******************
+**************************************
+Deprecated features and future changes
+**************************************
 - Functions `scipy.linalg.interpolative.rand` and
   `scipy.linalg.interpolative.seed` have been deprecated and will be removed
   in SciPy ``1.17.0``.
@@ -388,16 +386,15 @@ Deprecated features
 - `scipy.spatial.distance.kulczynski1` and
   `scipy.spatial.distance.sokalmichener` were deprecated and will be removed
   in SciPy ``1.17.0``.
-- `scipy.stats.find_repeats` is deprecated as of SciPy ``1.15.0`` and will be
+- `scipy.stats.find_repeats` is deprecated and will be
   removed in SciPy ``1.17.0``. Please use
   ``numpy.unique``/``numpy.unique_counts`` instead.
 - `scipy.linalg.kron` is deprecated in favour of ``numpy.kron``.
 - Using object arrays and longdouble arrays in `scipy.signal`
   convolution/correlation functions (`scipy.signal.correlate`,
   `scipy.signal.convolve` and `scipy.signal.choose_conv_method`) and
   filtering functions (`scipy.signal.lfilter`, `scipy.signal.sosfilt`) has
-  been deprecated as of SciPy ``1.15.0`` and will be removed in SciPy
-  ``1.17.0``.
+  been deprecated and will be removed in SciPy ``1.17.0``.
 - `scipy.stats.linregress` has deprecated one-argument use; the two
   variables must be specified as separate arguments.
 - ``scipy.stats.trapz`` is deprecated in favor of `scipy.stats.trapezoid`.
@@ -406,8 +403,9 @@ Deprecated features
   `scipy.special.assoc_legendre_p_all`.
 - `scipy.special.sph_harm` has been deprecated in favor of
   `scipy.special.sph_harm_y`.
-- The raveling of multi-dimensional input by `scipy.linalg.toeplitz` has
-  been deprecated. It will support batching in SciPy ``1.17.0``.
+- Multi-dimensional ``r`` and ``c`` arrays passed to `scipy.linalg.toeplitz`,
+  `scipy.linalg.matmul_toeplitz`, or `scipy.linalg.solve_toeplitz` will be
+  treated as batches of 1-D coefficients beginning in SciPy ``1.17.0``.
 - The ``random_state`` and ``permutations`` arguments of
   `scipy.stats.ttest_ind` are deprecated. Use ``method`` to perform a
   permutation test, instead.
@@ -557,11 +555,11 @@ Authors (commits)
 * Sahil Garje (1) +
 * Gabriel Gerlero (2)
 * Yotam Gingold (1) +
-* Ralf Gommers (110)
+* Ralf Gommers (111)
 * Rohit Goswami (62)
 * Anil Gurses (1) +
 * Oscar Gustafsson (1) +
-* Matt Haberland (389)
+* Matt Haberland (392)
 * Matt Hall (1) +
 * Joren Hammudoglu (6) +
 * CY Han (1) +
@@ -616,7 +614,7 @@ Authors (commits)
 * Tom M. Ragonneau (2)
 * Peter Ralph (1) +
 * Stephan Rave (1) +
-* Tyler Reddy (182)
+* Tyler Reddy (192)
 * redha2404 (2) +
 * Ritvik1sharma (1) +
 * Érico Nogueira Rolim (1) +
@@ -939,9 +937,11 @@ Issues closed for 1.15.0
 * `#22097 <https://github.com/scipy/scipy/issues/22097>`__: DEP: ``interpolate.interpnd.GradientEstimationWarning`` still...
 * `#22112 <https://github.com/scipy/scipy/issues/22112>`__: BUG/DOC: sparse: ND COO unexpected behaviour 1.15.0rc1
 * `#22123 <https://github.com/scipy/scipy/issues/22123>`__: DOC: stats: random variable transition guide launches wrong notebook
+* `#22128 <https://github.com/scipy/scipy/issues/22128>`__: BUG/DOC: it's not clear how to use ``differentiate.jacobian``...
 * `#22137 <https://github.com/scipy/scipy/issues/22137>`__: BUG: ``stats._distribution_infrastructure._Domain.symbols`` class...
 * `#22143 <https://github.com/scipy/scipy/issues/22143>`__: BUG: Fail to call ``BSpline`` after unpickling with ``mmap_mode="r"``
 * `#22146 <https://github.com/scipy/scipy/issues/22146>`__: BUG:``stats.ContinuousDistribution.llf``\ : should not be public
+* `#22204 <https://github.com/scipy/scipy/issues/22204>`__: BUG: signal.ShortTimeFFT: ``istft`` with ``mfft > len(win)``...
 
 ************************
 Pull requests for 1.15.0
@@ -1368,6 +1368,7 @@ Pull requests for 1.15.0
 * `#21668 <https://github.com/scipy/scipy/pull/21668>`__: BUG: fft.fht: set ``u.imag[-1] = 0`` only when ``n`` is even
 * `#21672 <https://github.com/scipy/scipy/pull/21672>`__: BUG: ndimage: fix 0d arrays in ``_normalize_sequence``
 * `#21673 <https://github.com/scipy/scipy/pull/21673>`__: BUG: signal.ShortTimeFFT: fix multichannel roundtrip with ``mfft``...
+* `#21678 <https://github.com/scipy/scipy/pull/21678>`__: BUG: fix ``nan`` output of ``special.betaincinv``
 * `#21680 <https://github.com/scipy/scipy/pull/21680>`__: MAINT:integrate: Silence a few QUADPACK compiler warnings
 * `#21682 <https://github.com/scipy/scipy/pull/21682>`__: DOC: Reduce duplication in user guide
 * `#21686 <https://github.com/scipy/scipy/pull/21686>`__: BUG: signal: int handling for ``resample_poly``
@@ -1573,11 +1574,15 @@ Pull requests for 1.15.0
 * `#22135 <https://github.com/scipy/scipy/pull/22135>`__: MAINT: _lib: add missing f string to _deprecate_positional_args
 * `#22139 <https://github.com/scipy/scipy/pull/22139>`__: MAINT: stats._SimpleDomain: ensure that instances do not share...
 * `#22149 <https://github.com/scipy/scipy/pull/22149>`__: MAINT: stats.ContinuousDistribution.llf: remove method
+* `#22150 <https://github.com/scipy/scipy/pull/22150>`__: MAINT: SciPy 1.15.0rc2 backports
 * `#22156 <https://github.com/scipy/scipy/pull/22156>`__: DEP: deprecation warnings for ``special.lpn`` and ``[c]lpmn``
 * `#22158 <https://github.com/scipy/scipy/pull/22158>`__: MAINT: accept ndarray subclasses in interpolate._dierckx
 * `#22162 <https://github.com/scipy/scipy/pull/22162>`__: TYP: temporarily ignore the ``numpy==2.2.1`` mypy errors
 * `#22167 <https://github.com/scipy/scipy/pull/22167>`__: DEP: special: deprecation warning for ``sph_harm`` + comments
 * `#22168 <https://github.com/scipy/scipy/pull/22168>`__: BUG: fix incorrect values in factorial for 0 with uint dtypes
 * `#22175 <https://github.com/scipy/scipy/pull/22175>`__: MAINT: stats: fix thread-safety issues under free-threaded CPython
 * `#22177 <https://github.com/scipy/scipy/pull/22177>`__: MAINT: fix extension module not declaring free-threading support,...
+* `#22181 <https://github.com/scipy/scipy/pull/22181>`__: REL: set 1.15.0rc3 unreleased
+* `#22193 <https://github.com/scipy/scipy/pull/22193>`__: DEP: linalg.solve_toeplitz/matmul_toeplitz: warn on n-D ``c``\...
+* `#22225 <https://github.com/scipy/scipy/pull/22225>`__: DOC: differentiate.jacobian: correct/improve documentation about...
 

scipy/differentiate/_differentiate.py

@@ -802,27 +802,47 @@ def jacobian(f, x, *, tolerances=None, maxiter=10, order=8, initial_step=0.5,
     Notes
     -----
     Suppose we wish to evaluate the Jacobian of a function
-    :math:`f: \mathbf{R}^m \rightarrow \mathbf{R}^n`, and assign to variables
+    :math:`f: \mathbf{R}^m \rightarrow \mathbf{R}^n`. Assign to variables
     ``m`` and ``n`` the positive integer values of :math:`m` and :math:`n`,
-    respectively. If we wish to evaluate the Jacobian at a single point,
-    then:
+    respectively, and let ``...`` represent an arbitrary tuple of integers.
+    If we wish to evaluate the Jacobian at a single point, then:
 
     - argument `x` must be an array of shape ``(m,)``
-    - argument `f` must be vectorized to accept an array of shape ``(m, p)``.
-      The first axis represents the :math:`m` inputs of :math:`f`; the second
-      is for evaluating the function at multiple points in a single call.
-    - argument `f` must return an array of shape ``(n, p)``. The first
-      axis represents the :math:`n` outputs of :math:`f`; the second
-      is for the result of evaluating the function at multiple points.
+    - argument `f` must be vectorized to accept an array of shape ``(m, ...)``.
+      The first axis represents the :math:`m` inputs of :math:`f`; the remainder
+      are for evaluating the function at multiple points in a single call.
+    - argument `f` must return an array of shape ``(n, ...)``. The first
+      axis represents the :math:`n` outputs of :math:`f`; the remainder
+      are for the result of evaluating the function at multiple points.
     - attribute ``df`` of the result object will be an array of shape ``(n, m)``,
       the Jacobian.
 
     This function is also vectorized in the sense that the Jacobian can be
     evaluated at ``k`` points in a single call. In this case, `x` would be an
     array of shape ``(m, k)``, `f` would accept an array of shape
-    ``(m, k, p)`` and return an array of shape ``(n, k, p)``, and the ``df``
+    ``(m, k, ...)`` and return an array of shape ``(n, k, ...)``, and the ``df``
     attribute of the result would have shape ``(n, m, k)``.
 
+    Suppose the desired callable ``f_not_vectorized`` is not vectorized; it can
+    only accept an array of shape ``(m,)``. A simple solution to satisfy the required
+    interface is to wrap ``f_not_vectorized`` as follows::
+
+        def f(x):
+            return np.apply_along_axis(f_not_vectorized, axis=0, arr=x)
+
+    Alternatively, suppose the desired callable ``f_vec_q`` is vectorized, but
+    only for 2-D arrays of shape ``(m, q)``. To satisfy the required interface,
+    consider::
+
+        def f(x):
+            m, batch = x.shape[0], x.shape[1:]  # x.shape is (m, ...)
+            x = np.reshape(x, (m, -1))  # `-1` is short for q = prod(batch)
+            res = f_vec_q(x)  # pass shape (m, q) to function
+            n = res.shape[0]
+            return np.reshape(res, (n,) + batch)  # return shape (n, ...)
+
+    Then pass the wrapped callable ``f`` as the first argument of `jacobian`.
+
     References
     ----------
     .. [1] Jacobian matrix and determinant, *Wikipedia*,

scipy/linalg/_basic.py

@@ -4,6 +4,7 @@
 # w/ additions by Travis Oliphant, March 2002
 #              and Jake Vanderplas, August 2012
 
+import warnings
 from warnings import warn
 from itertools import product
 import numpy as np
@@ -782,21 +783,25 @@ def solveh_banded(ab, b, overwrite_ab=False, overwrite_b=False, lower=False,
 
 
 def solve_toeplitz(c_or_cr, b, check_finite=True):
-    """Solve a Toeplitz system using Levinson Recursion
+    r"""Solve a Toeplitz system using Levinson Recursion
 
     The Toeplitz matrix has constant diagonals, with c as its first column
     and r as its first row. If r is not given, ``r == conjugate(c)`` is
     assumed.
 
+    .. warning::
+
+        Beginning in SciPy 1.17, multidimensional input will be treated as a batch,
+        not ``ravel``\ ed. To preserve the existing behavior, ``ravel`` arguments
+        before passing them to `solve_toeplitz`.
+
     Parameters
     ----------
     c_or_cr : array_like or tuple of (array_like, array_like)
-        The vector ``c``, or a tuple of arrays (``c``, ``r``). Whatever the
-        actual shape of ``c``, it will be converted to a 1-D array. If not
+        The vector ``c``, or a tuple of arrays (``c``, ``r``). If not
         supplied, ``r = conjugate(c)`` is assumed; in this case, if c[0] is
         real, the Toeplitz matrix is Hermitian. r[0] is ignored; the first row
-        of the Toeplitz matrix is ``[c[0], r[1:]]``. Whatever the actual shape
-        of ``r``, it will be converted to a 1-D array.
+        of the Toeplitz matrix is ``[c[0], r[1:]]``.
     b : (M,) or (M, K) array_like
         Right-hand side in ``T x = b``.
     check_finite : bool, optional
@@ -1915,12 +1920,21 @@ def _validate_args_for_toeplitz_ops(c_or_cr, b, check_finite, keep_b_shape,
 
     if isinstance(c_or_cr, tuple):
         c, r = c_or_cr
-        c = _asarray_validated(c, check_finite=check_finite).ravel()
-        r = _asarray_validated(r, check_finite=check_finite).ravel()
+        c = _asarray_validated(c, check_finite=check_finite)
+        r = _asarray_validated(r, check_finite=check_finite)
     else:
-        c = _asarray_validated(c_or_cr, check_finite=check_finite).ravel()
+        c = _asarray_validated(c_or_cr, check_finite=check_finite)
         r = c.conjugate()
 
+    if c.ndim > 1 or r.ndim > 1:
+        msg = ("Beginning in SciPy 1.17, multidimensional input will be treated as a "
+               "batch, not `ravel`ed. To preserve the existing behavior and silence "
+               "this warning, `ravel` arguments before passing them to "
+               "`toeplitz`, `matmul_toeplitz`, and `solve_toeplitz`.")
+        warnings.warn(msg, FutureWarning, stacklevel=2)
+        c = c.ravel()
+        r = r.ravel()
+
     if b is None:
         raise ValueError('`b` must be an array, not None.')
 
@@ -1944,7 +1958,7 @@ def _validate_args_for_toeplitz_ops(c_or_cr, b, check_finite, keep_b_shape,
 
 
 def matmul_toeplitz(c_or_cr, x, check_finite=False, workers=None):
-    """Efficient Toeplitz Matrix-Matrix Multiplication using FFT
+    r"""Efficient Toeplitz Matrix-Matrix Multiplication using FFT
 
     This function returns the matrix multiplication between a Toeplitz
     matrix and a dense matrix.
@@ -1953,15 +1967,19 @@ def matmul_toeplitz(c_or_cr, x, check_finite=False, workers=None):
     and r as its first row. If r is not given, ``r == conjugate(c)`` is
     assumed.
 
+    .. warning::
+
+        Beginning in SciPy 1.17, multidimensional input will be treated as a batch,
+        not ``ravel``\ ed. To preserve the existing behavior, ``ravel`` arguments
+        before passing them to `matmul_toeplitz`.
+
     Parameters
     ----------
     c_or_cr : array_like or tuple of (array_like, array_like)
-        The vector ``c``, or a tuple of arrays (``c``, ``r``). Whatever the
-        actual shape of ``c``, it will be converted to a 1-D array. If not
+        The vector ``c``, or a tuple of arrays (``c``, ``r``). If not
         supplied, ``r = conjugate(c)`` is assumed; in this case, if c[0] is
         real, the Toeplitz matrix is Hermitian. r[0] is ignored; the first row
-        of the Toeplitz matrix is ``[c[0], r[1:]]``. Whatever the actual shape
-        of ``r``, it will be converted to a 1-D array.
+        of the Toeplitz matrix is ``[c[0], r[1:]]``.
     x : (M,) or (M, K) array_like
         Matrix with which to multiply.
     check_finite : bool, optional

scipy/linalg/_special_matrices.py

@@ -27,19 +27,17 @@ def toeplitz(c, r=None):
     Parameters
     ----------
     c : array_like
-        First column of the matrix.  Whatever the actual shape of `c`, it
-        will be converted to a 1-D array.
+        First column of the matrix.
     r : array_like, optional
         First row of the matrix. If None, ``r = conjugate(c)`` is assumed;
         in this case, if c[0] is real, the result is a Hermitian matrix.
         r[0] is ignored; the first row of the returned matrix is
-        ``[c[0], r[1:]]``.  Whatever the actual shape of `r`, it will be
-        converted to a 1-D array.
+        ``[c[0], r[1:]]``.
 
         .. warning::
 
             Beginning in SciPy 1.17, multidimensional input will be treated as a batch,
-            not ``ravel``\ ed. To preserve the existing behavior, ``ravel`` aruments
+            not ``ravel``\ ed. To preserve the existing behavior, ``ravel`` arguments
             before passing them to `toeplitz`.
 
     Returns
@@ -81,7 +79,7 @@ def toeplitz(c, r=None):
     if c.ndim > 1 or r.ndim > 1:
         msg = ("Beginning in SciPy 1.17, multidimensional input will be treated as a "
                "batch, not `ravel`ed. To preserve the existing behavior and silence "
-               "this warning, `ravel` aruments before passing them to `toeplitz`.")
+               "this warning, `ravel` arguments before passing them to `toeplitz`.")
         warnings.warn(msg, FutureWarning, stacklevel=2)
 
     c, r = c.ravel(), r.ravel()

scipy/linalg/tests/test_matmul_toeplitz.py

@@ -125,11 +125,12 @@ def test_exceptions(self):
 
     # For toeplitz matrices, matmul_toeplitz() should be equivalent to @.
     def do(self, x, c, r=None, check_finite=False, workers=None):
+        c = np.ravel(c)
         if r is None:
             actual = matmul_toeplitz(c, x, check_finite, workers)
         else:
             r = np.ravel(r)
             actual = matmul_toeplitz((c, r), x, check_finite)
-        desired = toeplitz(np.ravel(c), r) @ x
+        desired = toeplitz(c, r) @ x
         assert_allclose(actual, desired,
             rtol=self.tolerance, atol=self.tolerance)

scipy/linalg/tests/test_solve_toeplitz.py

@@ -2,7 +2,7 @@
 """
 import numpy as np
 from scipy.linalg._solve_toeplitz import levinson
-from scipy.linalg import solve, toeplitz, solve_toeplitz
+from scipy.linalg import solve, toeplitz, solve_toeplitz, matmul_toeplitz
 from numpy.testing import assert_equal, assert_allclose
 
 import pytest
@@ -134,3 +134,17 @@ def test_empty(dt_c, dt_b):
     x1 = solve_toeplitz(c, b)
     assert x1.shape == (0, 0)
     assert x1.dtype == x.dtype
+
+
+@pytest.mark.parametrize('fun', [solve_toeplitz, matmul_toeplitz])
+def test_nd_FutureWarning(fun):
+    # Test future warnings with n-D `c`/`r`
+    rng = np.random.default_rng(283592436523456)
+    c = rng.random((2, 3, 4))
+    r = rng.random((2, 3, 4))
+    b_or_x = rng.random(24)
+    message = "Beginning in SciPy 1.17, multidimensional input will be..."
+    with pytest.warns(FutureWarning, match=message):
+         fun(c, b_or_x)
+    with pytest.warns(FutureWarning, match=message):
+         fun((c, r), b_or_x)