Bug#1105940: unblock: scipy/1.15.3-1
Package: release.debian.org
Followup-For: Bug #1105940
In the big scheme of things the debdiff is not so large,
though larger than usual for a typical unblock request.
Attached here.
diff -Nru scipy-1.15.2/benchmarks/benchmarks/io_mm.py scipy-1.15.3/benchmarks/benchmarks/io_mm.py
--- scipy-1.15.2/benchmarks/benchmarks/io_mm.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/benchmarks/benchmarks/io_mm.py 2025-05-08 04:14:18.000000000 +0200
@@ -46,7 +46,7 @@
def params(self):
return [
list(self._get_size().keys()),
- ['scipy.io', 'scipy.io._mmio', 'scipy.io._fast_matrix_market'],
+ ['scipy.io'],
['dense', 'coo'] # + ['csr']
]
@@ -54,8 +54,12 @@
size = {
'1M': int(1e6),
'10M': int(10e6),
- '100M': int(100e6),
- '300M': int(300e6),
+ '25M': int(10e6),
+ # Note: the below sizes should work locally but cause issues on CircleCI
+ # it fails to allocate memory even though there is easily enough
+ # available (see gh-22574).
+ #'100M': int(100e6),
+ #'300M': int(300e6),
# '500M': int(500e6),
# '1000M': int(1000e6),
}
diff -Nru scipy-1.15.2/debian/changelog scipy-1.15.3/debian/changelog
--- scipy-1.15.2/debian/changelog 2025-04-14 23:05:55.000000000 +0200
+++ scipy-1.15.3/debian/changelog 2025-05-11 11:41:19.000000000 +0200
@@ -1,3 +1,11 @@
+scipy (1.15.3-1) unstable; urgency=medium
+
+ * New upstream bug-fix release
+ - applies debian patches ndimage_pytest_safety_PR22828.patch
+ and ndimage_segfault_filter_mirror_PR22608.patch
+
+ -- Drew Parsons <dparsons@debian.org> Sun, 11 May 2025 11:41:19 +0200
+
scipy (1.15.2-8) unstable; urgency=medium
* debian patch ndimage_pytest_safety_PR22828.patch applies upstream
diff -Nru scipy-1.15.2/debian/patches/ndimage_pytest_safety_PR22828.patch scipy-1.15.3/debian/patches/ndimage_pytest_safety_PR22828.patch
--- scipy-1.15.2/debian/patches/ndimage_pytest_safety_PR22828.patch 2025-04-14 23:05:55.000000000 +0200
+++ scipy-1.15.3/debian/patches/ndimage_pytest_safety_PR22828.patch 1970-01-01 01:00:00.000000000 +0100
@@ -1,22 +0,0 @@
-From dd595d7d9fdce495e0210ba79de8508704d5ea96 Mon Sep 17 00:00:00 2001
-From: Lucas Colley <lucas.colley8@gmail.com>
-Date: Sat, 12 Apr 2025 17:56:03 +0100
-Subject: [PATCH] BUG: add workaround for pytest assertion rewriting overreach
-
-closes gh-22236
-see discussion of upstream bug from https://github.com/pytest-dev/pytest/issues/10845#issuecomment-2575032929
----
- scipy/ndimage/_ndimage_api.py | 3 ++-
- 1 file changed, 2 insertions(+), 1 deletion(-)
-
-diff --git a/scipy/ndimage/_ndimage_api.py b/scipy/ndimage/_ndimage_api.py
-index d19e4594ba34..1673391726a4 100644
---- a/scipy/ndimage/_ndimage_api.py
-+++ b/scipy/ndimage/_ndimage_api.py
-@@ -12,4 +12,5 @@
- from ._measurements import * # noqa: F403
- from ._morphology import * # noqa: F403
-
--__all__ = [s for s in dir() if not s.startswith('_')]
-+# '@' due to pytest bug, scipy/scipy#22236
-+__all__ = [s for s in dir() if not s.startswith(('_', '@'))]
diff -Nru scipy-1.15.2/debian/patches/ndimage_segfault_filter_mirror_PR22608.patch scipy-1.15.3/debian/patches/ndimage_segfault_filter_mirror_PR22608.patch
--- scipy-1.15.2/debian/patches/ndimage_segfault_filter_mirror_PR22608.patch 2025-04-14 23:05:55.000000000 +0200
+++ scipy-1.15.3/debian/patches/ndimage_segfault_filter_mirror_PR22608.patch 1970-01-01 01:00:00.000000000 +0100
@@ -1,106 +0,0 @@
-From 45c3825b83192f6701bee88cdf50bbe62103af29 Mon Sep 17 00:00:00 2001
-From: Roger Aiudi <aiudirog@gmail.com>
-Date: Fri, 28 Feb 2025 14:53:40 -0500
-Subject: [PATCH 1/3] BUG:ndimage.median_filter: fix segfault when using the
- mode='mirror'
-
-It appears there was an off-by-one error in the loop that finalizes the mirror mode. Additionally, added a test to directly compare the new 1D optimized implementation to the general ND one. Closes gh-22586
----
- scipy/ndimage/src/_rank_filter_1d.cpp | 9 +++----
- scipy/ndimage/tests/test_filters.py | 37 +++++++++++++++++++++++++++
- 2 files changed, 40 insertions(+), 6 deletions(-)
-
-Index: scipy/scipy/ndimage/src/_rank_filter_1d.cpp
-===================================================================
---- scipy.orig/scipy/ndimage/src/_rank_filter_1d.cpp 2025-03-06 17:48:47.919991800 +0100
-+++ scipy/scipy/ndimage/src/_rank_filter_1d.cpp 2025-03-06 17:48:47.907991689 +0100
-@@ -157,11 +157,7 @@
- int lim2 = arr_len - lim;
- int offset;
- Mediator *m = MediatorNew(win_len, rank);
-- T *data = new T[win_len];
-- for (int i = 0; i < win_len; ++i) {
-- data[i] = 0;
-- }
--
-+ T *data = new T[win_len]();
-
- switch (mode) {
- case REFLECT:
-@@ -204,6 +200,7 @@
- MediatorInsert(data, m, in_arr[i]);
- out_arr[i - lim] = data[m->heap[0]];
- }
-+
- switch (mode) {
- case REFLECT:
- arr_len_thresh = arr_len - 1;
-@@ -227,7 +224,7 @@
- break;
- case MIRROR:
- arr_len_thresh = arr_len - 2;
-- for (i = 0; i < lim + 1; i++) {
-+ for (i = 0; i < lim; i++) {
- MediatorInsert(data, m, in_arr[arr_len_thresh - i]);
- out_arr[lim2 + i] = data[m->heap[0]];
- }
-Index: scipy/scipy/ndimage/tests/test_filters.py
-===================================================================
---- scipy.orig/scipy/ndimage/tests/test_filters.py 2025-03-06 17:48:47.919991800 +0100
-+++ scipy/scipy/ndimage/tests/test_filters.py 2025-03-06 17:48:47.911991726 +0100
-@@ -6,6 +6,10 @@
- import numpy as np
- import pytest
- from numpy.testing import suppress_warnings, assert_allclose, assert_array_equal
-+from numpy.typing import NDArray
-+from hypothesis import strategies as st
-+from hypothesis import given
-+import hypothesis.extra.numpy as npst
- from pytest import raises as assert_raises
- from scipy import ndimage
- from scipy._lib._array_api import (
-@@ -2904,3 +2908,44 @@
- expected = [58, 67, 87, 108, 163, 108, 108, 108, 87]
- actual = ndimage.median_filter(x, size=9, mode='constant')
- assert_array_equal(actual, expected)
-+
-+
-+@pytest.mark.slow
-+@given(
-+ arr=npst.arrays(
-+ dtype=np.float64,
-+ shape=st.integers(min_value=1, max_value=20).map(lambda x: x * 50),
-+ ),
-+ size=st.integers(min_value=1, max_value=25),
-+ mode=st.sampled_from(['reflect', 'constant', 'nearest', 'mirror', 'wrap']),
-+)
-+def test_gh_22586_median_filter_no_segfault_and_1d_matches_nd(
-+ arr: NDArray[np.float64],
-+ size: int,
-+ mode: str,
-+ monkeypatch: pytest.MonkeyPatch,
-+):
-+ size = min(size, len(arr))
-+ mf = functools.partial(
-+ ndimage.median_filter,
-+ size=size,
-+ mode=mode,
-+ )
-+ actual = mf(arr)
-+
-+ # Monkeypatch _rank_filter_1d.rank_filter() to actually call
-+ # _nd_image.rank_filter() so that 1D results from the optimized
-+ # implementation can be compared to the general ND implementation
-+
-+ def _1d_to_nd_rank_filter(i, r, s, *args):
-+ from scipy.ndimage import _nd_image
-+ return _nd_image.rank_filter(i, r, np.ones(s, dtype=bool), *args)
-+
-+ with monkeypatch.context():
-+ monkeypatch.setattr(
-+ 'scipy.ndimage._rank_filter_1d.rank_filter',
-+ _1d_to_nd_rank_filter,
-+ )
-+ desired = mf(arr)
-+
-+ np.testing.assert_allclose(actual, desired)
diff -Nru scipy-1.15.2/debian/patches/series scipy-1.15.3/debian/patches/series
--- scipy-1.15.2/debian/patches/series 2025-04-14 23:05:55.000000000 +0200
+++ scipy-1.15.3/debian/patches/series 2025-05-11 11:41:19.000000000 +0200
@@ -10,5 +10,3 @@
special_libsf_error_SOABI.patch
scipy_config_SOABI.patch
test_nan_policy_array_22360.patch
-ndimage_segfault_filter_mirror_PR22608.patch
-ndimage_pytest_safety_PR22828.patch
diff -Nru scipy-1.15.2/doc/source/release/1.15.3-notes.rst scipy-1.15.3/doc/source/release/1.15.3-notes.rst
--- scipy-1.15.2/doc/source/release/1.15.3-notes.rst 1970-01-01 01:00:00.000000000 +0100
+++ scipy-1.15.3/doc/source/release/1.15.3-notes.rst 2025-05-08 04:14:18.000000000 +0200
@@ -0,0 +1,108 @@
+==========================
+SciPy 1.15.3 Release Notes
+==========================
+
+.. contents::
+
+SciPy 1.15.3 is a bug-fix release with no new features
+compared to 1.15.2.
+
+
+
+Authors
+=======
+* Name (commits)
+* aiudirog (1) +
+* Nickolai Belakovski (1)
+* Florian Bourgey (1) +
+* Richard Strong Bowen (2) +
+* Jake Bowhay (1)
+* Dietrich Brunn (2)
+* Evgeni Burovski (1)
+* Lucas Colley (1)
+* Ralf Gommers (1)
+* Saarthak Gupta (1) +
+* Matt Haberland (4)
+* Chengyu Han (1) +
+* Lukas Huber (1) +
+* Nick ODell (2)
+* Ilhan Polat (4)
+* Tyler Reddy (52)
+* Neil Schemenauer (1) +
+* Dan Schult (1)
+* sildater (1) +
+* Gagandeep Singh (4)
+* Albert Steppi (2)
+* Matthias Urlichs (1) +
+* David Varela (1) +
+* ਗਗਨਦੀਪ ਸਿੰਘ (Gagandeep Singh) (3)
+
+A total of 24 people contributed to this release.
+People with a "+" by their names contributed a patch for the first time.
+This list of names is automatically generated, and may not be fully complete.
+
+
+Issues closed for 1.15.3
+------------------------
+
+* `#10634 <https://github.com/scipy/scipy/issues/10634>`__: BUG: optimize: ``least_squares`` with ``'trf'`` and ``'trf_sover=lsmr'``...
+* `#18146 <https://github.com/scipy/scipy/issues/18146>`__: BUG: scipy.sparse.linalg.expm_multiply fails with sparse matrices
+* `#19418 <https://github.com/scipy/scipy/issues/19418>`__: BUG: integrate.solve_ivp fails for some step sizes if dense_output=True...
+* `#19865 <https://github.com/scipy/scipy/issues/19865>`__: BUG: HalfspaceIntersection.add_halfspaces() does not seem to...
+* `#20988 <https://github.com/scipy/scipy/issues/20988>`__: BUG: special.hyp2f1: wrong result for extreme inputs
+* `#22236 <https://github.com/scipy/scipy/issues/22236>`__: BUG: scipy v1.15 breaking for pytest when assert-rewrite is on
+* `#22400 <https://github.com/scipy/scipy/issues/22400>`__: BUG: stats.genextreme.stats: Spurious warning from ``genextreme.stats(0.0,``...
+* `#22451 <https://github.com/scipy/scipy/issues/22451>`__: BUG: interpolative svd broken for non-square linear operators
+* `#22515 <https://github.com/scipy/scipy/issues/22515>`__: CI: Some GitHub workflows failing due to check on ``actions/cache``...
+* `#22547 <https://github.com/scipy/scipy/issues/22547>`__: BUG: _lib: Data race reported by TSAN in ``ccallback`` mechanism
+* `#22558 <https://github.com/scipy/scipy/issues/22558>`__: BUG: linalg.expm: bug on Windows / conda
+* `#22574 <https://github.com/scipy/scipy/issues/22574>`__: CI: benchmark job on CircleCI is failing on ``io.mmread`` memory...
+* `#22586 <https://github.com/scipy/scipy/issues/22586>`__: BUG: ndimage.median_filter: additional hard crashes
+* `#22589 <https://github.com/scipy/scipy/issues/22589>`__: BUG: spatial: ``Rotation`` no longer supports zero-length collections
+* `#22599 <https://github.com/scipy/scipy/issues/22599>`__: DOC: sparse.linalg.ArpackError: entire default ``infodict`` displayed
+* `#22615 <https://github.com/scipy/scipy/issues/22615>`__: CI: oneAPI job: ``Not enough disk space.``
+* `#22637 <https://github.com/scipy/scipy/issues/22637>`__: BUG: Transposed LinearOperator fails on vector multiplication
+* `#22655 <https://github.com/scipy/scipy/issues/22655>`__: BUG: optimize.linprog: 40x slower in v1.15 compared to v1.14
+* `#22681 <https://github.com/scipy/scipy/issues/22681>`__: DOC: integrate.tanhsinh: documentation refers to non-existent...
+* `#22684 <https://github.com/scipy/scipy/issues/22684>`__: BUG: signal.resample_poly: dtype not preserved
+* `#22720 <https://github.com/scipy/scipy/issues/22720>`__: MAINT, CI: floating point exceptions activated in NumPy
+* `#22868 <https://github.com/scipy/scipy/issues/22868>`__: BUG: re-importing ``scipy`` fails
+* `#22903 <https://github.com/scipy/scipy/issues/22903>`__: BUG: special.logsumexp: nan in 1.15
+
+
+Pull requests for 1.15.3
+------------------------
+
+* `#20035 <https://github.com/scipy/scipy/pull/20035>`__: BUG: spatial.HalfspaceIntersection: raise on non-feasible half...
+* `#22502 <https://github.com/scipy/scipy/pull/22502>`__: BUG: special: Fix typo in specfun::chgu
+* `#22517 <https://github.com/scipy/scipy/pull/22517>`__: CI: Use actions/cache 4.2.0
+* `#22532 <https://github.com/scipy/scipy/pull/22532>`__: BUG: Remove warning for genextreme.stats(0.0, moments='mvsk')
+* `#22543 <https://github.com/scipy/scipy/pull/22543>`__: REL, MAINT: prep for 1.15.3
+* `#22555 <https://github.com/scipy/scipy/pull/22555>`__: BUG: ``scipy.sparse.linalg``\ : Fix ``expm_multiply`` if both...
+* `#22561 <https://github.com/scipy/scipy/pull/22561>`__: BUG: _lib: Fix data race found by TSAN, use SCIPY_TLS.
+* `#22567 <https://github.com/scipy/scipy/pull/22567>`__: BUG: optimize: Fix ``bracket_root`` termination check and default...
+* `#22582 <https://github.com/scipy/scipy/pull/22582>`__: BUG: ``integrate.solve_ivp``\ : Avoid duplicate time stamps in...
+* `#22587 <https://github.com/scipy/scipy/pull/22587>`__: BUG: Pin jupyterlite-sphinx to >= 0.19.1
+* `#22588 <https://github.com/scipy/scipy/pull/22588>`__: BUG/BLD: xsf: force defining the mdspan parenthesis operator...
+* `#22590 <https://github.com/scipy/scipy/pull/22590>`__: BENCH: remove triple run of mmread/mmwrite benchmark, limit sizes
+* `#22600 <https://github.com/scipy/scipy/pull/22600>`__: BUG: Fix ArpackError default argument
+* `#22608 <https://github.com/scipy/scipy/pull/22608>`__: BUG: ndimage.median_filter: fix segfault when using ``mode='mirror'``
+* `#22617 <https://github.com/scipy/scipy/pull/22617>`__: CI: minimise disk space usage for oneAPI jobs
+* `#22642 <https://github.com/scipy/scipy/pull/22642>`__: BUG: sparse: sparse sum/mean out parameter shape not enforced...
+* `#22643 <https://github.com/scipy/scipy/pull/22643>`__: BUG: spatial.transform.Rotation: support 0-length rotations
+* `#22660 <https://github.com/scipy/scipy/pull/22660>`__: BUG: optimize: avoid expensive access of ``basis.col_status``...
+* `#22689 <https://github.com/scipy/scipy/pull/22689>`__: BUG: signal.resample_poly: fix dtype preservation
+* `#22690 <https://github.com/scipy/scipy/pull/22690>`__: MAINT/DOC: integrate.tanhsinh: lightly refactor error estimate...
+* `#22693 <https://github.com/scipy/scipy/pull/22693>`__: BUG: spatial.HalfspaceIntersection: fix ``add_halfspaces`` batch...
+* `#22726 <https://github.com/scipy/scipy/pull/22726>`__: MAINT: compensate for dot exceptions
+* `#22763 <https://github.com/scipy/scipy/pull/22763>`__: BUG: sparse: Remove reference cycle to improve memory use
+* `#22772 <https://github.com/scipy/scipy/pull/22772>`__: BUG: sparse.linalg: Transposed ``LinearOperator`` multiplication...
+* `#22784 <https://github.com/scipy/scipy/pull/22784>`__: BUG: signal._short_time_fft: incorrect index computation in ``upper_border_begin``...
+* `#22792 <https://github.com/scipy/scipy/pull/22792>`__: BUG: signal.ShortTimeFFT.upper_border_begin: Document parameter...
+* `#22801 <https://github.com/scipy/scipy/pull/22801>`__: BUG: ``signal.windows._windows.kaiser_bessel_derived``\ : use...
+* `#22810 <https://github.com/scipy/scipy/pull/22810>`__: BUG: special.hyp2f1: fix for extreme inputs
+* `#22822 <https://github.com/scipy/scipy/pull/22822>`__: BUG: linalg.expm: Fix noncompliant compiler branch typos in C...
+* `#22828 <https://github.com/scipy/scipy/pull/22828>`__: BUG: add workaround for pytest assertion rewriting overreach
+* `#22834 <https://github.com/scipy/scipy/pull/22834>`__: BUG: linalg: Fix shape mismatch in interpolative.svd
+* `#22869 <https://github.com/scipy/scipy/pull/22869>`__: BUG: optimize._highspy: don't import from inside a C module
+* `#22910 <https://github.com/scipy/scipy/pull/22910>`__: MAINT: special.logsumexp: improvement when weight of largest...
diff -Nru scipy-1.15.2/doc/source/release.rst scipy-1.15.3/doc/source/release.rst
--- scipy-1.15.2/doc/source/release.rst 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/doc/source/release.rst 2025-05-08 04:14:18.000000000 +0200
@@ -8,6 +8,7 @@
.. toctree::
:maxdepth: 1
+ release/1.15.3-notes
release/1.15.2-notes
release/1.15.1-notes
release/1.15.0-notes
diff -Nru scipy-1.15.2/pyproject.toml scipy-1.15.3/pyproject.toml
--- scipy-1.15.2/pyproject.toml 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/pyproject.toml 2025-05-08 04:14:18.000000000 +0200
@@ -43,7 +43,7 @@
[project]
name = "scipy"
-version = "1.15.2"
+version = "1.15.3"
# 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" }
@@ -111,7 +111,7 @@
"jupytext",
"myst-nb",
"pooch",
- "jupyterlite-sphinx>=0.16.5",
+ "jupyterlite-sphinx>=0.19.1",
"jupyterlite-pyodide-kernel",
]
dev = [
diff -Nru scipy-1.15.2/requirements/doc.txt scipy-1.15.3/requirements/doc.txt
--- scipy-1.15.2/requirements/doc.txt 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/requirements/doc.txt 2025-05-08 04:14:18.000000000 +0200
@@ -10,5 +10,5 @@
jupytext
myst-nb
pooch
-jupyterlite-sphinx>=0.16.5
+jupyterlite-sphinx>=0.19.1
jupyterlite-pyodide-kernel
diff -Nru scipy-1.15.2/scipy/integrate/_ivp/common.py scipy-1.15.3/scipy/integrate/_ivp/common.py
--- scipy-1.15.2/scipy/integrate/_ivp/common.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/integrate/_ivp/common.py 2025-05-08 04:14:18.000000000 +0200
@@ -65,7 +65,7 @@
return np.linalg.norm(x) / x.size ** 0.5
-def select_initial_step(fun, t0, y0, t_bound,
+def select_initial_step(fun, t0, y0, t_bound,
max_step, f0, direction, order, rtol, atol):
"""Empirically select a good initial step.
@@ -80,7 +80,7 @@
y0 : ndarray, shape (n,)
Initial value of the dependent variable.
t_bound : float
- End-point of integration interval; used to ensure that t0+step<=tbound
+ End-point of integration interval; used to ensure that t0+step<=tbound
and that fun is only evaluated in the interval [t0,tbound]
max_step : float
Maximum allowable step size.
@@ -112,7 +112,7 @@
interval_length = abs(t_bound - t0)
if interval_length == 0.0:
return 0.0
-
+
scale = atol + np.abs(y0) * rtol
d0 = norm(y0 / scale)
d1 = norm(f0 / scale)
diff -Nru scipy-1.15.2/scipy/integrate/_ivp/ivp.py scipy-1.15.3/scipy/integrate/_ivp/ivp.py
--- scipy-1.15.2/scipy/integrate/_ivp/ivp.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/integrate/_ivp/ivp.py 2025-05-08 04:14:18.000000000 +0200
@@ -694,8 +694,15 @@
g = g_new
if t_eval is None:
- ts.append(t)
- ys.append(y)
+ donot_append = (len(ts) > 1 and
+ ts[-1] == t and
+ dense_output)
+ if not donot_append:
+ ts.append(t)
+ ys.append(y)
+ else:
+ if len(interpolants) > 0:
+ interpolants.pop()
else:
# The value in t_eval equal to t will be included.
if solver.direction > 0:
diff -Nru scipy-1.15.2/scipy/integrate/_ivp/tests/test_ivp.py scipy-1.15.3/scipy/integrate/_ivp/tests/test_ivp.py
--- scipy-1.15.2/scipy/integrate/_ivp/tests/test_ivp.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/integrate/_ivp/tests/test_ivp.py 2025-05-08 04:14:18.000000000 +0200
@@ -151,6 +151,25 @@
e = (y - y_true) / (atol + rtol * np.abs(y_true))
return np.linalg.norm(e, axis=0) / np.sqrt(e.shape[0])
+def test_duplicate_timestamps():
+ def upward_cannon(t, y):
+ return [y[1], -9.80665]
+
+ def hit_ground(t, y):
+ return y[0]
+
+ hit_ground.terminal = True
+ hit_ground.direction = -1
+
+ sol = solve_ivp(upward_cannon, [0, np.inf], [0, 0.01],
+ max_step=0.05 * 0.001 / 9.80665,
+ events=hit_ground, dense_output=True)
+ assert_allclose(sol.sol(0.01), np.asarray([-0.00039033, -0.08806632]),
+ rtol=1e-5, atol=1e-8)
+ assert_allclose(sol.t_events, np.asarray([[0.00203943]]), rtol=1e-5, atol=1e-8)
+ assert_allclose(sol.y_events, [np.asarray([[ 0.0, -0.01 ]])], atol=1e-9)
+ assert sol.success
+ assert_equal(sol.status, 1)
@pytest.mark.thread_unsafe
def test_integration():
diff -Nru scipy-1.15.2/scipy/integrate/_tanhsinh.py scipy-1.15.3/scipy/integrate/_tanhsinh.py
--- scipy-1.15.2/scipy/integrate/_tanhsinh.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/integrate/_tanhsinh.py 2025-05-08 04:14:18.000000000 +0200
@@ -98,8 +98,9 @@
Absolute termination tolerance (default: 0) and relative termination
tolerance (default: ``eps**0.75``, where ``eps`` is the precision of
the result dtype), respectively. Iteration will stop when
- ``res.error < atol + rtol * abs(res.df)``. The error estimate is as
- described in [1]_ Section 5. While not theoretically rigorous or
+ ``res.error < atol`` or ``res.error < res.integral * rtol``. The error
+ estimate is as described in [1]_ Section 5 but with a lower bound of
+ ``eps * res.integral``. While not theoretically rigorous or
conservative, it is said to work well in practice. Must be non-negative
and finite if `log` is False, and must be expressed as the log of a
non-negative and finite number if `log` is True.
@@ -445,9 +446,9 @@
stop[i] = True
else:
# Terminate if convergence criterion is met
- work.rerr, work.aerr = _estimate_error(work, xp)
- i = ((work.rerr < rtol) | (work.rerr + xp_real(work.Sn) < atol) if log
- else (work.rerr < rtol) | (work.rerr * xp.abs(work.Sn) < atol))
+ rerr, aerr = _estimate_error(work, xp)
+ i = (rerr < rtol) | (aerr < atol)
+ work.aerr = xp.reshape(xp.astype(aerr, work.dtype), work.Sn.shape)
work.status[i] = eim._ECONVERGED
stop[i] = True
@@ -772,22 +773,23 @@
d2 = xp_real(special.logsumexp(xp.stack([work.Sn, Snm2 + work.pi*1j]), axis=0))
d3 = log_e1 + xp.max(xp_real(work.fjwj), axis=-1)
d4 = work.d4
- ds = xp.stack([d1 ** 2 / d2, 2 * d1, d3, d4])
+ d5 = log_e1 + xp.real(work.Sn)
+ temp = xp.where(d1 > -xp.inf, d1 ** 2 / d2, -xp.inf)
+ ds = xp.stack([temp, 2 * d1, d3, d4, d5])
aerr = xp.max(ds, axis=0)
- rerr = xp.maximum(log_e1, aerr - xp_real(work.Sn))
+ rerr = aerr - xp.real(work.Sn)
else:
# Note: explicit computation of log10 of each of these is unnecessary.
d1 = xp.abs(work.Sn - Snm1)
d2 = xp.abs(work.Sn - Snm2)
d3 = e1 * xp.max(xp.abs(work.fjwj), axis=-1)
d4 = work.d4
- # If `d1` is 0, no need to warn. This does the right thing.
- # with np.errstate(divide='ignore'):
- ds = xp.stack([d1**(xp.log(d1)/xp.log(d2)), d1**2, d3, d4])
+ d5 = e1 * xp.abs(work.Sn)
+ temp = xp.where(d1 > 0, d1**(xp.log(d1)/xp.log(d2)), 0)
+ ds = xp.stack([temp, d1**2, d3, d4, d5])
aerr = xp.max(ds, axis=0)
- rerr = xp.maximum(e1, aerr/xp.abs(work.Sn))
+ rerr = aerr/xp.abs(work.Sn)
- aerr = xp.reshape(xp.astype(aerr, work.dtype), work.Sn.shape)
return rerr, aerr
diff -Nru scipy-1.15.2/scipy/integrate/tests/test_tanhsinh.py scipy-1.15.3/scipy/integrate/tests/test_tanhsinh.py
--- scipy-1.15.2/scipy/integrate/tests/test_tanhsinh.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/integrate/tests/test_tanhsinh.py 2025-05-08 04:14:18.000000000 +0200
@@ -753,6 +753,16 @@
x[-1] = 1000
_tanhsinh(np.sin, 1, x)
+ def test_gh_22681_finite_error(self, xp):
+ # gh-22681 noted a case in which the error was NaN on some platforms;
+ # check that this does in fact fail in CI.
+ a = complex(12, -10)
+ b = complex(12, 39)
+ def f(t):
+ return xp.sin(a * (1 - t) + b * t)
+ res = _tanhsinh(f, xp.asarray(0.), xp.asarray(1.), atol=0, rtol=0, maxlevel=10)
+ assert xp.isfinite(res.error)
+
@array_api_compatible
@pytest.mark.usefixtures("skip_xp_backends")
diff -Nru scipy-1.15.2/scipy/_lib/_array_api.py scipy-1.15.3/scipy/_lib/_array_api.py
--- scipy-1.15.2/scipy/_lib/_array_api.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/_lib/_array_api.py 2025-05-08 04:14:18.000000000 +0200
@@ -593,3 +593,14 @@
arr = xp.astype(arr, xp.complex128)
return arr
+
+
+def xp_default_dtype(xp):
+ """Query the namespace-dependent default floating-point dtype.
+ """
+ if is_torch(xp):
+ # historically, we allow pytorch to keep its default of float32
+ return xp.get_default_dtype()
+ else:
+ # we default to float64
+ return xp.float64
diff -Nru scipy-1.15.2/scipy/_lib/src/ccallback.h scipy-1.15.3/scipy/_lib/src/ccallback.h
--- scipy-1.15.2/scipy/_lib/src/ccallback.h 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/_lib/src/ccallback.h 2025-05-08 04:14:18.000000000 +0200
@@ -244,7 +244,7 @@
static int ccallback_prepare(ccallback_t *callback, ccallback_signature_t *signatures,
PyObject *callback_obj, int flags)
{
- static PyTypeObject *lowlevelcallable_type = NULL;
+ static SCIPY_TLS PyTypeObject *lowlevelcallable_type = NULL;
PyObject *callback_obj2 = NULL;
PyObject *capsule = NULL;
diff -Nru scipy-1.15.2/scipy/_lib/tests/test_array_api.py scipy-1.15.3/scipy/_lib/tests/test_array_api.py
--- scipy-1.15.2/scipy/_lib/tests/test_array_api.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/_lib/tests/test_array_api.py 2025-05-08 04:14:18.000000000 +0200
@@ -4,7 +4,7 @@
from scipy.conftest import array_api_compatible
from scipy._lib._array_api import (
_GLOBAL_CONFIG, array_namespace, _asarray, xp_copy, xp_assert_equal, is_numpy,
- np_compat,
+ np_compat, xp_default_dtype
)
from scipy._lib._array_api_no_0d import xp_assert_equal as xp_assert_equal_no_0d
@@ -185,3 +185,7 @@
# scalars-vs-0d passes (if values match) also with regular python objects
xp_assert_equal_no_0d(0., xp.asarray(0.))
xp_assert_equal_no_0d(42, xp.asarray(42))
+
+ @array_api_compatible
+ def test_default_dtype(self, xp):
+ assert xp_default_dtype(xp) == xp.asarray(1.).dtype
diff -Nru scipy-1.15.2/scipy/linalg/_decomp_interpolative.pyx scipy-1.15.3/scipy/linalg/_decomp_interpolative.pyx
--- scipy-1.15.2/scipy/linalg/_decomp_interpolative.pyx 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/linalg/_decomp_interpolative.pyx 2025-05-08 04:14:18.000000000 +0200
@@ -1057,14 +1057,14 @@
def iddr_rsvd(A: LinearOperator, int krank, *, rng):
- cdef int n = A.shape[1], j
+ cdef int m = A.shape[0], n = A.shape[1], j
cdef cnp.ndarray[cnp.int64_t, mode='c', ndim=1] perms
cdef cnp.ndarray[cnp.float64_t, ndim=2] proj
cdef cnp.ndarray[cnp.float64_t, mode='c', ndim=2] col
perms, proj = iddr_rid(A, krank, rng=rng)
# idd_getcols
- col = cnp.PyArray_EMPTY(2, [n, krank], cnp.NPY_FLOAT64, 0)
+ col = cnp.PyArray_EMPTY(2, [m, krank], cnp.NPY_FLOAT64, 0)
x = cnp.PyArray_ZEROS(1, [n], cnp.NPY_FLOAT64, 0)
for j in range(krank):
x[perms[j]] = 1.
diff -Nru scipy-1.15.2/scipy/linalg/interpolative.pyf scipy-1.15.3/scipy/linalg/interpolative.pyf
--- scipy-1.15.2/scipy/linalg/interpolative.pyf 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/linalg/interpolative.pyf 1970-01-01 01:00:00.000000000 +0100
@@ -1,693 +0,0 @@
-!*******************************************************************************
-! Copyright (C) 2013 Kenneth L. Ho
-!
-! 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.
-!
-! None of the names of the copyright holders may be used to endorse or promote
-! products derived from this software without specific prior written
-! permission.
-!
-! 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.
-!*******************************************************************************
-
-!===============================================================================
-! main module
-!===============================================================================
-
- python module _interpolative
- interface
-
-! ----------------------------------------------------------------------------
-! id_rand.f
-! ----------------------------------------------------------------------------
-
- subroutine id_srand(n,r)
- integer, intent(in) :: n
- real*8, intent(out), depend(n) :: r(n)
- real*8, intent(in) :: t(55)
- entry id_srandi(t)
- entry id_srando()
- end subroutine id_srand
-
-! ----------------------------------------------------------------------------
-! idd_frm.f
-! ----------------------------------------------------------------------------
-
- subroutine idd_frm(m,n,w,x,y)
- integer, integer(in), optional, depend(x) :: m = len(x)
- integer, intent(in) :: n
- real*8, intent(in), depend(m) :: w(17*m+70)
- real*8, intent(in) :: x(m)
- real*8, intent(out), depend(n) :: y(n)
- end subroutine idd_frm
-
- subroutine idd_sfrm(l,m,n,w,x,y)
- integer, intent(in), check(l<=n), depend(n) :: l
- integer, intent(in), optional, depend(x) :: m = len(x)
- integer, intent(in) :: n
- real*8, intent(in), depend(m), :: w(27*m+90)
- real*8, intent(in) :: x(m)
- real*8, intent(out), depend(l) :: y(l)
- end subroutine idd_sfrm
-
- subroutine idd_frmi(m,n,w)
- integer, intent(in) :: m
- integer, intent(out) :: n
- real*8, intent(out), depend(m) :: w(17*m+70)
- end subroutine idd_frmi
-
- subroutine idd_sfrmi(l,m,n,w)
- integer, intent(in) :: l, m
- integer, intent(out) :: n
- real*8, intent(out), depend(m) :: w(27*m+90)
- end subroutine idd_sfrmi
-
-! ----------------------------------------------------------------------------
-! idd_id.f
-! ----------------------------------------------------------------------------
-
- subroutine iddp_id(eps,m,n,a,krank,list,rnorms)
- real*8, intent(in) :: eps, a(m,n)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- integer, intent(out) :: krank
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(out), depend(n) :: rnorms(n)
- end subroutine iddp_id
-
- subroutine iddr_id(m,n,a,krank,list,rnorms)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(out), depend(n) :: rnorms(n)
- end subroutine iddr_id
-
- subroutine idd_reconid(m,krank,col,n,list,proj,approx)
- integer, intent(in), optional, depend(col) :: m = shape(col,0), krank = shape(col,1)
- real*8, intent(in) :: col(m,krank)
- integer, intent(in), depend(list) :: n = len(list)
- integer, intent(in) :: list(n)
- real*8, intent(in), depend(krank,n) :: proj(krank,n-krank)
- real*8, intent(out), depend(m,n) :: approx(m,n)
- end subroutine idd_reconid
-
- subroutine idd_reconint(n,list,krank,proj,p)
- integer, intent(in), optional, depend(list) :: n = len(list)
- integer, intent(in) :: list(n)
- integer, intent(in), optional, depend(proj) :: krank = shape(proj,0)
- real*8, intent(in), depend(n) :: proj(krank,n-krank)
- real*8, intent(out), depend(krank,n) :: p(krank,n)
- end subroutine idd_reconint
-
- subroutine idd_copycols(m,n,a,krank,list,col)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in) :: a(m,n)
- integer, intent(in) :: krank, list(*)
- real*8, intent(out), depend(m,krank) :: col(m,krank)
- end subroutine idd_copycols
-
-! ----------------------------------------------------------------------------
-! idd_id2svd.f
-! ----------------------------------------------------------------------------
-
- subroutine idd_id2svd(m,krank,b,n,list,proj,u,v,s,ier,w)
- integer, intent(in), optional, depend(b) :: m = shape(b,0), krank = shape(b,1)
- real*8, intent(in) :: b(m,krank)
- integer, intent(in), optional, depend(list) :: n = len(list)
- integer, intent(in) :: list(n)
- real*8, intent(in), depend(krank,n) :: proj(krank,n-krank)
- real*8, intent(out), depend(m,krank) :: u(m,krank)
- real*8, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- real*8, intent(in), optional, depend(m,n,krank) :: w((krank+1)*(m+3*n)+26*pow(krank,2))
- end subroutine idd_id2svd
-
-! ----------------------------------------------------------------------------
-! idd_snorm.f
-! ----------------------------------------------------------------------------
-
- subroutine idd_snorm(m,n,matvect,p1t,p2t,p3t,p4t,matvec,p1,p2,p3,p4,its,snorm,v,u)
- use idd__user__routines
- integer, intent(in) :: m, n, its
- external matvect, matvec
- real*8, intent(in), optional :: p1t, p2t, p3t, p4t, p1, p2, p3, p4
- real*8, intent(out) :: snorm
- real*8, intent(out), depend(n) :: v(n)
- real*8, intent(in), optional, depend(m) :: u(m)
- end subroutine idd_snorm
-
- subroutine idd_diffsnorm(m,n,matvect,p1t,p2t,p3t,p4t,matvect2,p1t2,p2t2,p3t2,p4t2,matvec,p1,p2,p3,p4,matvec2,p12,p22,p32,p42,its,snorm,w)
- use idd__user__routines
- integer, intent(in) :: m, n, its
- external matvect, matvect2, matvec, matvec2
- real*8, intent(in), optional :: p1t, p2t, p3t, p4t, p1t2, p2t2, p3t2, p4t2, p1, p2, p3, p4, p12, p22, p32, p42
- real*8, intent(out) :: snorm
- real*8, intent(in), optional, depend(m,n) :: w(3*(m+n))
- end subroutine idd_diffsnorm
-
-! ----------------------------------------------------------------------------
-! idd_svd.f
-! ----------------------------------------------------------------------------
-
- subroutine iddr_svd(m,n,a,krank,u,v,s,ier,r)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- real*8, intent(out), depend(m,krank) :: u(m,krank)
- real*8, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- real*8, intent(in), optional, depend(m,n,krank) :: r((krank+2)*n+8*min(m,n)+15*pow(krank,2)+8*krank)
- end subroutine iddr_svd
-
- subroutine iddp_svd(lw,eps,m,n,a,krank,iu,iv,is,w,ier)
- integer, intent(hide), optional, depend(m,n) :: lw = (min(m,n)+1)*(m+2*n+9)+8*min(m,n)+15*pow(min(m,n),2)
- real*8, intent(in) :: eps, a(m,n)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- integer, intent(out) :: krank, iu, iv, is, ier
- real*8, intent(out), depend(m,n) :: w((min(m,n)+1)*(m+2*n+9)+8*min(m,n)+15*pow(min(m,n),2))
- end subroutine iddp_svd
-
-! ----------------------------------------------------------------------------
-! iddp_aid.f
-! ----------------------------------------------------------------------------
-
- subroutine iddp_aid(eps,m,n,a,work,krank,list,proj)
- real*8, intent(in) :: eps, a(m,n)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in), depend(m) :: work(17*m+70)
- integer, intent(out) :: krank
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(in,out) :: proj(*)
- end subroutine iddp_aid
-
- subroutine idd_estrank(eps,m,n,a,w,krank,ra)
- real*8, intent(in) :: eps, a(m,n)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in), depend(m) :: w(17*m+70)
- integer, intent(out) :: krank
- real*8, intent(in,out) :: ra(*)
- end subroutine idd_estrank
-
-! ----------------------------------------------------------------------------
-! iddp_asvd.f
-! ----------------------------------------------------------------------------
-
- subroutine iddp_asvd(lw,eps,m,n,a,winit,krank,iu,iv,is,w,ier)
- integer, intent(hide), optional, depend(m,n) :: lw = max((min(m,n)+1)*(3*m+5*n+1)+25*pow(min(m,n),2),(2*n+1)*(m+1))
- real*8, intent(in) :: eps, a(m,n)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in), depend(m) :: winit(17*m+70)
- integer, intent(out) :: krank, iu, iv, is, ier
- real*8, intent(in,out) :: w(*)
- end subroutine iddp_asvd
-
-! ----------------------------------------------------------------------------
-! iddp_rid.f
-! ----------------------------------------------------------------------------
-
- subroutine iddp_rid(lproj,eps,m,n,matvect,p1,p2,p3,p4,krank,list,proj,ier)
- use idd__user__routines
- integer, intent(hide), optional, depend(m,n) :: lproj = m+1+2*n*(min(m,n)+1)
- real*8, intent(in) :: eps
- integer, intent(in) :: m, n
- external matvect
- real*8, intent(in), optional :: p1, p2, p3, p4
- integer, intent(out) :: krank, ier
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(in,out) :: proj(*)
- end subroutine iddp_rid
-
- subroutine idd_findrank(lra,eps,m,n,matvect,p1,p2,p3,p4,krank,ra,ier,w)
- use idd__user__routines
- integer, intent(hide), optional, depend(m,n) :: lra = 2*n*min(m,n)
- real*8, intent(in) :: eps
- integer, intent(in) :: m, n
- external matvect
- real*8, intent(in), optional :: p1, p2, p3, p4
- integer, intent(out) :: krank, ier
- real*8, intent(out), depend(m,n) :: ra(2*n*min(m,n))
- real*8, intent(in), optional, depend(m,n) :: w(m+2*n+1)
- end subroutine idd_findrank
-
-! ----------------------------------------------------------------------------
-! iddp_rsvd.f
-! ----------------------------------------------------------------------------
-
- subroutine iddp_rsvd(lw,eps,m,n,matvect,p1t,p2t,p3t,p4t,matvec,p1,p2,p3,p4,krank,iu,iv,is,w,ier)
- use idd__user__routines
- integer, intent(hide), optional, depend(m,n) :: lw = (min(m,n)+1)*(3*m+5*n+1)+25*pow(min(m,n),2)
- real*8, intent(in) :: eps
- integer, intent(in) :: m ,n
- external matvect, matvec
- real*8, intent(in), optional :: p1t, p2t, p3t, p4t, p1, p2, p3, p4
- integer, intent(out) :: krank, iu, iv, is, ier
- real*8, intent(out), depend(m,n) :: w((min(m,n)+1)*(3*m+5*n+1)+25*pow(min(m,n),2))
- end subroutine iddp_rsvd
-
-! ----------------------------------------------------------------------------
-! iddr_aid.f
-! ----------------------------------------------------------------------------
-
- subroutine iddr_aid(m,n,a,krank,w,list,proj)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- real*8, intent(in), depend(m,n,krank) :: w((2*krank+17)*n+27*m+100)
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(out), depend(n,krank) :: proj(max(krank*(n-krank),1))
- end subroutine iddr_aid
-
- subroutine iddr_aidi(m,n,krank,w)
- integer, intent(in) :: m, n, krank
- real*8, intent(out), depend(m,n,krank) :: w((2*krank+17)*n+27*m+100)
- end subroutine iddr_aidi
-
-! ----------------------------------------------------------------------------
-! iddr_asvd.f
-! ----------------------------------------------------------------------------
-
- subroutine iddr_asvd(m,n,a,krank,w,u,v,s,ier)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- real*8, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- real*8, intent(in), depend(m,n,krank) :: w((2*krank+28)*m+(6*krank+21)*n+25*pow(krank,2)+100)
- real*8, intent(out), depend(m,krank) :: u(m,krank)
- real*8, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- end subroutine iddr_asvd
-
-! ----------------------------------------------------------------------------
-! iddr_rid.f
-! ----------------------------------------------------------------------------
-
- subroutine iddr_rid(m,n,matvect,p1,p2,p3,p4,krank,list,proj)
- use idd__user__routines
- integer, intent(in) :: m, n, krank
- external matvect
- real*8, intent(in), optional :: p1, p2, p3, p4
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(out), depend(m,n,krank) :: proj(m+(krank+3)*n)
- end subroutine iddr_rid
-
-! ----------------------------------------------------------------------------
-! iddr_rsvd.f
-! ----------------------------------------------------------------------------
-
- subroutine iddr_rsvd(m,n,matvect,p1t,p2t,p3t,p4t,matvec,p1,p2,p3,p4,krank,u,v,s,ier,w)
- use idd__user__routines
- integer, intent(in) :: m, n, krank
- external matvect, matvec
- real*8, intent(in), optional :: p1t, p2t, p3t, p4t, p1, p2, p3, p4
- real*8, intent(out), depend(m,krank) :: u(m,krank)
- real*8, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- real*8, intent(in), optional, depend(m,n,krank) :: w((krank+1)*(2*m+4*n)+25*pow(krank,2))
- end subroutine iddr_rsvd
-
-! ----------------------------------------------------------------------------
-! idz_frm.f
-! ----------------------------------------------------------------------------
-
- subroutine idz_frm(m,n,w,x,y)
- integer, integer(in), optional, depend(x) :: m = len(x)
- integer, intent(in) :: n
- complex*16, intent(in), depend(m) :: w(17*m+70)
- complex*16, intent(in) :: x(m)
- complex*16, intent(out), depend(n) :: y(n)
- end subroutine idz_frm
-
- subroutine idz_sfrm(l,m,n,w,x,y)
- integer, intent(in), check(l<=n), depend(n) :: l
- integer, intent(in), optional, depend(x) :: m = len(x)
- integer, intent(in) :: n
- complex*16, intent(in), depend(m), :: w(27*m+90)
- complex*16, intent(in) :: x(m)
- complex*16, intent(out), depend(l) :: y(l)
- end subroutine idz_sfrm
-
- subroutine idz_frmi(m,n,w)
- integer, intent(in) :: m
- integer, intent(out) :: n
- complex*16, intent(out), depend(m) :: w(17*m+70)
- end subroutine idz_frmi
-
- subroutine idz_sfrmi(l,m,n,w)
- integer, intent(in) :: l, m
- integer, intent(out) :: n
- complex*16, intent(out), depend(m) :: w(27*m+90)
- end subroutine idz_sfrmi
-
-! ----------------------------------------------------------------------------
-! idz_id.f
-! ----------------------------------------------------------------------------
-
- subroutine idzp_id(eps,m,n,a,krank,list,rnorms)
- real*8, intent(in) :: eps
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(out) :: krank
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(out), depend(n) :: rnorms(n)
- end subroutine idzp_id
-
- subroutine idzr_id(m,n,a,krank,list,rnorms)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- integer, intent(out), depend(n) :: list(n)
- real*8, intent(out), depend(n) :: rnorms(n)
- end subroutine idzr_id
-
- subroutine idz_reconid(m,krank,col,n,list,proj,approx)
- integer, intent(in), optional, depend(col) :: m = shape(col,0), krank = shape(col,1)
- complex*16, intent(in) :: col(m,krank)
- integer, intent(in), depend(list) :: n = len(list)
- integer, intent(in) :: list(n)
- complex*16, intent(in), depend(krank,n) :: proj(krank,n-krank)
- complex*16, intent(out), depend(m,n) :: approx(m,n)
- end subroutine idz_reconid
-
- subroutine idz_reconint(n,list,krank,proj,p)
- integer, intent(in), optional, depend(list) :: n = len(list)
- integer, intent(in) :: list(n)
- integer, intent(in), optional, depend(proj) :: krank = shape(proj,0)
- complex*16, intent(in), depend(n) :: proj(krank,n-krank)
- complex*16, intent(out), depend(krank,n) :: p(krank,n)
- end subroutine idz_reconint
-
- subroutine idz_copycols(m,n,a,krank,list,col)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(in) :: krank, list(*)
- complex*16, intent(out), depend(m,krank) :: col(m,krank)
- end subroutine idz_copycols
-
-! ----------------------------------------------------------------------------
-! idz_id2svd.f
-! ----------------------------------------------------------------------------
-
- subroutine idz_id2svd(m,krank,b,n,list,proj,u,v,s,ier,w)
- integer, intent(in), optional, depend(b) :: m = shape(b,0), krank = shape(b,1)
- complex*16, intent(in) :: b(m,krank)
- integer, intent(in), optional, depend(list) :: n = len(list)
- integer, intent(in) :: list(n)
- complex*16, intent(in), depend(krank,n) :: proj(krank,n-krank)
- complex*16, intent(out), depend(m,krank) :: u(m,krank)
- complex*16, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- complex*16, intent(in), optional, depend(m,n,krank) :: w((krank+1)*(m+3*n+10)+9*pow(krank,2))
- end subroutine idz_id2svd
-
-! ----------------------------------------------------------------------------
-! idz_snorm.f
-! ----------------------------------------------------------------------------
-
- subroutine idz_snorm(m,n,matveca,p1a,p2a,p3a,p4a,matvec,p1,p2,p3,p4,its,snorm,v,u)
- use idz__user__routines
- integer, intent(in) :: m, n, its
- external matveca, matvec
- complex*16, intent(in), optional :: p1a, p2a, p3a, p4a, p1, p2, p3, p4
- real*8, intent(out) :: snorm
- complex*16, intent(out), depend(n) :: v(n)
- complex*16, intent(in), optional, depend(m) :: u(m)
- end subroutine idz_snorm
-
- subroutine idz_diffsnorm(m,n,matveca,p1a,p2a,p3a,p4a,matveca2,p1a2,p2a2,p3a2,p4a2,matvec,p1,p2,p3,p4,matvec2,p12,p22,p32,p42,its,snorm,w)
- use idz__user__routines
- integer, intent(in) :: m, n, its
- external matveca, matveca2, matvec, matvec2
- complex*16, intent(in), optional :: p1a, p2a, p3a, p4a, p1a2, p2a2, p3a2, p4a2, p1, p2, p3, p4, p12, p22, p32, p42
- real*8, intent(out) :: snorm
- complex*16, intent(in), optional, depend(m,n) :: w(3*(m+n))
- end subroutine idz_diffsnorm
-
-! ----------------------------------------------------------------------------
-! idz_svd.f
-! ----------------------------------------------------------------------------
-
- subroutine idzr_svd(m,n,a,krank,u,v,s,ier,r)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- complex*16, intent(out), depend(m,krank) :: u(m,krank)
- complex*16, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- complex*16, intent(in), optional, depend(m,n,krank) :: r((krank+2)*n+8*min(m,n)+6*pow(krank,2)+8*krank)
- end subroutine idzr_svd
-
- subroutine idzp_svd(lw,eps,m,n,a,krank,iu,iv,is,w,ier)
- integer, intent(hide), optional, depend(m,n) :: lw = (min(m,n)+1)*(m+2*n+9)+8*min(m,n)+6*pow(min(m,n),2)
- real*8, intent(in) :: eps
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(out) :: krank, iu, iv, is, ier
- complex*16, intent(out), depend(m,n) :: w((min(m,n)+1)*(m+2*n+9)+8*min(m,n)+6*pow(min(m,n),2))
- end subroutine idzp_svd
-
-! ----------------------------------------------------------------------------
-! idzp_aid.f
-! ----------------------------------------------------------------------------
-
- subroutine idzp_aid(eps,m,n,a,work,krank,list,proj)
- real*8, intent(in) :: eps
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- complex*16, intent(in), depend(m) :: work(17*m+70)
- integer, intent(out) :: krank
- integer, intent(out), depend(n) :: list(n)
- complex*16, intent(in,out) :: proj(*)
- end subroutine idzp_aid
-
- subroutine idz_estrank(eps,m,n,a,w,krank,ra)
- real*8, intent(in) :: eps
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- complex*16, intent(in), depend(m) :: w(17*m+70)
- integer, intent(out) :: krank
- complex*16, intent(in,out) :: ra(*)
- end subroutine idz_estrank
-
-! ----------------------------------------------------------------------------
-! idzp_asvd.f
-! ----------------------------------------------------------------------------
-
- subroutine idzp_asvd(lw,eps,m,n,a,winit,krank,iu,iv,is,w,ier)
- integer, intent(hide), optional, depend(m,n) :: lw = max((min(m,n)+1)*(3*m+5*n+11)+8*pow(min(m,n),2),(2*n+1)*(m+1))
- real*8, intent(in) :: eps
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- complex*16, intent(in), depend(m) :: winit(17*m+70)
- integer, intent(out) :: krank, iu, iv, is, ier
- complex*16, intent(in,out) :: w(*)
- end subroutine idzp_asvd
-
-! ----------------------------------------------------------------------------
-! idzp_rid.f
-! ----------------------------------------------------------------------------
-
- subroutine idzp_rid(lproj,eps,m,n,matveca,p1,p2,p3,p4,krank,list,proj,ier)
- use idz__user__routines
- integer, intent(hide), optional, depend(m,n) :: lproj = m+1+2*n*(min(m,n)+1)
- real*8, intent(in) :: eps
- integer, intent(in) :: m, n
- external matveca
- complex*16, intent(in), optional :: p1, p2, p3, p4
- integer, intent(out) :: krank, ier
- integer, intent(out), depend(n) :: list(n)
- complex*16, intent(in,out) :: proj(*)
- end subroutine idzp_rid
-
- subroutine idz_findrank(lra,eps,m,n,matveca,p1,p2,p3,p4,krank,ra,ier,w)
- use idz__user__routines
- integer, intent(hide), optional, depend(m,n) :: lra = 2*n*min(m,n)
- real*8, intent(in) :: eps
- integer, intent(in) :: m, n
- external matveca
- complex*16, intent(in), optional :: p1, p2, p3, p4
- integer, intent(out) :: krank, ier
- complex*16, intent(out), depend(m,n) :: ra(2*n*min(m,n))
- complex*16, intent(in), optional, depend(m,n) :: w(m+2*n+1)
- end subroutine idz_findrank
-
-! ----------------------------------------------------------------------------
-! idzp_rsvd.f
-! ----------------------------------------------------------------------------
-
- subroutine idzp_rsvd(lw,eps,m,n,matveca,p1a,p2a,p3a,p4a,matvec,p1,p2,p3,p4,krank,iu,iv,is,w,ier)
- use idz__user__routines
- integer, intent(hide), optional, depend(m,n) :: lw = (min(m,n)+1)*(3*m+5*n+11)+8*pow(min(m,n),2)
- real*8, intent(in) :: eps
- integer, intent(in) :: m ,n
- external matveca, matvec
- complex*16, intent(in), optional :: p1a, p2a, p3a, p4a, p1, p2, p3, p4
- integer, intent(out) :: krank, iu, iv, is, ier
- complex*16, intent(out), depend(m,n) :: w((min(m,n)+1)*(3*m+5*n+11)+8*pow(min(m,n),2))
- end subroutine idzp_rsvd
-
-! ----------------------------------------------------------------------------
-! idzr_aid.f
-! ----------------------------------------------------------------------------
-
- subroutine idzr_aid(m,n,a,krank,w,list,proj)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- complex*16, intent(in), depend(m,n,krank) :: w((2*krank+17)*n+21*m+80)
- integer, intent(out), depend(n) :: list(n)
- complex*16, intent(out), depend(n,krank) :: proj(max(krank*(n-krank),1))
- end subroutine idzr_aid
-
- subroutine idzr_aidi(m,n,krank,w)
- integer, intent(in) :: m, n, krank
- complex*16, intent(out), depend(m,n,krank) :: w((2*krank+17)*n+21*m+80)
- end subroutine idzr_aidi
-
-! ----------------------------------------------------------------------------
-! idzr_asvd.f
-! ----------------------------------------------------------------------------
-
- subroutine idzr_asvd(m,n,a,krank,w,u,v,s,ier)
- integer, intent(in), optional, depend(a) :: m = shape(a,0), n = shape(a,1)
- complex*16, intent(in) :: a(m,n)
- integer, intent(in) :: krank
- complex*16, intent(in), depend(m,n,krank) :: w((2*krank+22)*m+(6*krank+21)*n+8*pow(krank,2)+10*krank+90)
- complex*16, intent(out), depend(m,krank) :: u(m,krank)
- complex*16, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- end subroutine idzr_asvd
-
-! ----------------------------------------------------------------------------
-! idzr_rid.f
-! ----------------------------------------------------------------------------
-
- subroutine idzr_rid(m,n,matveca,p1,p2,p3,p4,krank,list,proj)
- use idz__user__routines
- integer, intent(in) :: m, n, krank
- external matveca
- complex*16, intent(in), optional :: p1, p2, p3, p4
- integer, intent(out), depend(n) :: list(n)
- complex*16, intent(out), depend(m,n,krank) :: proj(m+(krank+3)*n)
- end subroutine idzr_rid
-
-! ----------------------------------------------------------------------------
-! idzr_rsvd.f
-! ----------------------------------------------------------------------------
-
- subroutine idzr_rsvd(m,n,matveca,p1a,p2a,p3a,p4a,matvec,p1,p2,p3,p4,krank,u,v,s,ier,w)
- use idz__user__routines
- integer, intent(in) :: m, n, krank
- external matveca, matvec
- complex*16, intent(in), optional :: p1a, p2a, p3a, p4a, p1, p2, p3, p4
- complex*16, intent(out), depend(m,krank) :: u(m,krank)
- complex*16, intent(out), depend(n,krank) :: v(n,krank)
- real*8, intent(out), depend(krank) :: s(krank)
- integer, intent(out) :: ier
- complex*16, intent(in), optional, depend(m,n,krank) :: w((krank+1)*(2*m+4*n+10)+8*pow(krank,2))
- end subroutine idzr_rsvd
-
- end interface
- end python module _interpolative
-
-!===============================================================================
-! auxiliary modules
-!===============================================================================
-
- python module idd__user__routines
- interface idd_user_interface
-
- subroutine matvect(m,x,n,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(u) :: m = len(x)
- real*8, intent(in) :: x(m)
- integer, intent(in), optional :: n
- real*8, intent(out), depend(n) :: y(n)
- real*8, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvect
-
- subroutine matvec(n,x,m,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(v) :: n = len(x)
- real*8, intent(in) :: x(n)
- integer, intent(in), optional :: m
- real*8, intent(out), depend(m) :: y(m)
- real*8, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvec
-
- subroutine matvect2(m,x,n,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(u) :: m = len(x)
- real*8, intent(in) :: x(m)
- integer, intent(in), optional :: n
- real*8, intent(out), depend(n) :: y(n)
- real*8, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvect2
-
- subroutine matvec2(n,x,m,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(v) :: n = len(x)
- real*8, intent(in) :: x(n)
- integer, intent(in), optional :: m
- real*8, intent(out), depend(m) :: y(m)
- real*8, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvec2
-
- end interface idd_user_interface
- end python module idd__user__routines
-
- python module idz__user__routines
- interface idz_user_interface
-
- subroutine matveca(m,x,n,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(u) :: m = len(x)
- complex*16, intent(in) :: x(m)
- integer, intent(in), optional :: n
- complex*16, intent(out), depend(n) :: y(n)
- complex*16, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvect
-
- subroutine matvec(n,x,m,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(v) :: n = len(x)
- complex*16, intent(in) :: x(n)
- integer, intent(in), optional :: m
- complex*16, intent(out), depend(m) :: y(m)
- complex*16, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvec
-
- subroutine matveca2(m,x,n,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(u) :: m = len(x)
- complex*16, intent(in) :: x(m)
- integer, intent(in), optional :: n
- complex*16, intent(out), depend(n) :: y(n)
- complex*16, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvect2
-
- subroutine matvec2(n,x,m,y,p1,p2,p3,p4)
- integer, intent(in), optional, depend(v) :: n = len(x)
- complex*16, intent(in) :: x(n)
- integer, intent(in), optional :: m
- complex*16, intent(out), depend(m) :: y(m)
- complex*16, intent(in), optional :: p1, p2, p3, p4
- end subroutine matvec2
-
- end interface idz_user_interface
- end python module idz__user__routines
diff -Nru scipy-1.15.2/scipy/linalg/_matfuncs_expm.c scipy-1.15.3/scipy/linalg/_matfuncs_expm.c
--- scipy-1.15.2/scipy/linalg/_matfuncs_expm.c 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/linalg/_matfuncs_expm.c 2025-05-08 04:14:18.000000000 +0200
@@ -1843,8 +1843,8 @@
inter1 = _FCmulcr(Am[n2 + n*i + j], b[2]);
inter2 = _FCmulcr(Am[2*n2 + n*i + j], b[4]);
inter3 = _FCmulcr(Am[3*n2 + n*i + j], b[6]);
- Am[4*n2 + n*i + j] = CPLX_C(crealf(inter1) + crealf(inter2) + crealf(inter3),
- cimagf(inter1) + cimagf(inter2) + cimagf(inter3));
+ Am[n2 + n*i + j] = CPLX_C(crealf(inter1) + crealf(inter2) + crealf(inter3),
+ cimagf(inter1) + cimagf(inter2) + cimagf(inter3));
if (i == j) { Am[n2 + n*i + j] = CPLX_C(crealf(Am[n2 + n*i + j]) + b[0], cimagf(Am[n2 + n*i + j])); }
#else
Am[n2 + n*i + j] = b[6]*Am[3*n2 + n*i + j] + b[4]*Am[2*n2 + n*i + j] + b[2]*Am[n2 + n*i + j];
@@ -2032,7 +2032,7 @@
csscal_(&n2, &two, &Am[2*n2], &int1);
#if defined(_MSC_VER)
- for (i = 0; i < n; i++) { Am[2*n2 + n*i + i] = CPLX_C(crealf(Am[n2 + n*i + j]) + b[0], cimagf(Am[n2 + n*i + j])); }
+ for (i = 0; i < n; i++) { Am[2*n2 + n*i + i] = CPLX_C(crealf(Am[2*n2 + n*i + i]) + 1.0f, cimagf(Am[2*n2 + n*i + i])); }
#else
for (i = 0; i < n; i++) { Am[2*n2 + n*i + i] += cdbl1; }
#endif
@@ -2179,8 +2179,8 @@
inter1 = _Cmulcr(Am[n2 + n*i + j], b[2]);
inter2 = _Cmulcr(Am[2*n2 + n*i + j], b[4]);
inter3 = _Cmulcr(Am[3*n2 + n*i + j], b[6]);
- Am[4*n2 + n*i + j] = CPLX_Z(creal(inter1) + creal(inter2) + creal(inter3),
- cimag(inter1) + cimag(inter2) + cimag(inter3));
+ Am[n2 + n*i + j] = CPLX_Z(creal(inter1) + creal(inter2) + creal(inter3),
+ cimag(inter1) + cimag(inter2) + cimag(inter3));
if (i == j) { Am[n2 + n*i + j] = CPLX_Z(creal(Am[n2 + n*i + j]) + b[0], cimag(Am[n2 + n*i + j])); }
#else
Am[n2 + n*i + j] = b[6]*Am[3*n2 + n*i + j] + b[4]*Am[2*n2 + n*i + j] + b[2]*Am[n2 + n*i + j];
diff -Nru scipy-1.15.2/scipy/linalg/tests/test_interpolative.py scipy-1.15.3/scipy/linalg/tests/test_interpolative.py
--- scipy-1.15.2/scipy/linalg/tests/test_interpolative.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/linalg/tests/test_interpolative.py 2025-05-08 04:14:18.000000000 +0200
@@ -213,3 +213,20 @@
B = A.copy()
interp_decomp(A.T, eps, rand=rand)
assert_array_equal(A, B)
+
+ def test_svd_aslinearoperator_shape_check(self):
+ # See gh-issue #22451
+ rng = np.random.default_rng(1744580941832515)
+ x = rng.uniform(size=[7, 5])
+ xl = aslinearoperator(x)
+ u, s, v = pymatrixid.svd(xl, 3)
+ assert_equal(u.shape, (7, 3))
+ assert_equal(s.shape, (3,))
+ assert_equal(v.shape, (5, 3))
+
+ x = rng.uniform(size=[4, 9])
+ xl = aslinearoperator(x)
+ u, s, v = pymatrixid.svd(xl, 2)
+ assert_equal(u.shape, (4, 2))
+ assert_equal(s.shape, (2,))
+ assert_equal(v.shape, (9, 2))
diff -Nru scipy-1.15.2/scipy/ndimage/_ndimage_api.py scipy-1.15.3/scipy/ndimage/_ndimage_api.py
--- scipy-1.15.2/scipy/ndimage/_ndimage_api.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/ndimage/_ndimage_api.py 2025-05-08 04:14:18.000000000 +0200
@@ -12,4 +12,5 @@
from ._measurements import * # noqa: F403
from ._morphology import * # noqa: F403
-__all__ = [s for s in dir() if not s.startswith('_')]
+# '@' due to pytest bug, scipy/scipy#22236
+__all__ = [s for s in dir() if not s.startswith(('_', '@'))]
diff -Nru scipy-1.15.2/scipy/ndimage/src/_rank_filter_1d.cpp scipy-1.15.3/scipy/ndimage/src/_rank_filter_1d.cpp
--- scipy-1.15.2/scipy/ndimage/src/_rank_filter_1d.cpp 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/ndimage/src/_rank_filter_1d.cpp 2025-05-08 04:14:18.000000000 +0200
@@ -155,13 +155,10 @@
int mode, T cval, int origin) {
int i, arr_len_thresh, lim = (win_len - 1) / 2 - origin;
int lim2 = arr_len - lim;
+ if (lim2 < 0) return;
int offset;
Mediator *m = MediatorNew(win_len, rank);
- T *data = new T[win_len];
- for (int i = 0; i < win_len; ++i) {
- data[i] = 0;
- }
-
+ T *data = new T[win_len]();
switch (mode) {
case REFLECT:
@@ -227,7 +224,7 @@
break;
case MIRROR:
arr_len_thresh = arr_len - 2;
- for (i = 0; i < lim + 1; i++) {
+ for (i = 0; i < lim; i++) {
MediatorInsert(data, m, in_arr[arr_len_thresh - i]);
out_arr[lim2 + i] = data[m->heap[0]];
}
diff -Nru scipy-1.15.2/scipy/ndimage/tests/test_filters.py scipy-1.15.3/scipy/ndimage/tests/test_filters.py
--- scipy-1.15.2/scipy/ndimage/tests/test_filters.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/ndimage/tests/test_filters.py 2025-05-08 04:14:18.000000000 +0200
@@ -6,6 +6,9 @@
import numpy as np
import pytest
from numpy.testing import suppress_warnings, assert_allclose, assert_array_equal
+from hypothesis import strategies as st
+from hypothesis import given
+import hypothesis.extra.numpy as npst
from pytest import raises as assert_raises
from scipy import ndimage
from scipy._lib._array_api import (
@@ -2904,3 +2907,14 @@
expected = [58, 67, 87, 108, 163, 108, 108, 108, 87]
actual = ndimage.median_filter(x, size=9, mode='constant')
assert_array_equal(actual, expected)
+
+
+@given(x=npst.arrays(dtype=np.float64,
+ shape=st.integers(min_value=1, max_value=1000)),
+ size=st.integers(min_value=1, max_value=50),
+ mode=st.sampled_from(["constant", "mirror", "wrap", "reflect",
+ "nearest"]),
+ )
+def test_gh_22586_crash_property(x, size, mode):
+ # property-based test for median_filter resilience to hard crashing
+ ndimage.median_filter(x, size=size, mode=mode)
diff -Nru scipy-1.15.2/scipy/optimize/_bracket.py scipy-1.15.3/scipy/optimize/_bracket.py
--- scipy-1.15.2/scipy/optimize/_bracket.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/optimize/_bracket.py 2025-05-08 04:14:18.000000000 +0200
@@ -1,7 +1,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 array_namespace, xp_ravel
+from scipy._lib._array_api import array_namespace, xp_ravel, xp_default_dtype
_ELIMITS = -1 # used in _bracket_root
_ESTOPONESIDE = 2 # used in _bracket_root
@@ -19,8 +19,17 @@
if (not xp.isdtype(xl0.dtype, "numeric")
or xp.isdtype(xl0.dtype, "complex floating")):
raise ValueError('`xl0` must be numeric and real.')
+ if not xp.isdtype(xl0.dtype, "real floating"):
+ xl0 = xp.asarray(xl0, dtype=xp_default_dtype(xp))
+
+ # If xr0 is not supplied, fill with a dummy value for the sake of
+ # broadcasting. We need to wait until xmax has been validated to
+ # compute the default value.
+ xr0_not_supplied = False
+ if xr0 is None:
+ xr0 = xp.nan
+ xr0_not_supplied = True
- xr0 = xl0 + 1 if xr0 is None else xr0
xmin = -xp.inf if xmin is None else xmin
xmax = xp.inf if xmax is None else xmax
factor = 2. if factor is None else factor
@@ -45,6 +54,12 @@
if not xp.all(factor > 1):
raise ValueError('All elements of `factor` must be greater than 1.')
+ # Calculate the default value of xr0 if a value has not been supplied.
+ # Be careful to ensure xr0 is not larger than xmax.
+ if xr0_not_supplied:
+ xr0 = xl0 + xp.minimum((xmax - xl0)/ 8, xp.asarray(1.0))
+ xr0 = xp.astype(xr0, xl0.dtype, copy=False)
+
maxiter = xp.asarray(maxiter)
message = '`maxiter` must be a non-negative integer.'
if (not xp.isdtype(maxiter.dtype, "numeric") or maxiter.shape != tuple()
@@ -281,14 +296,22 @@
# `work.status`.
# Get the integer indices of the elements that can also stop
also_stop = (work.active[i] + work.n) % (2*work.n)
- # Check whether they are still active.
- # To start, we need to find out where in `work.active` they would
- # appear if they are indeed there.
+ # Check whether they are still active. We want to find the indices
+ # in work.active where the associated values in work.active are
+ # contained in also_stop. xp.searchsorted let's us take advantage
+ # of work.active being sorted, but requires some hackery because
+ # searchsorted solves the separate but related problem of finding
+ # the indices where the values in also_stop should be added to
+ # maintain sorted order.
j = xp.searchsorted(work.active, also_stop)
# If the location exceeds the length of the `work.active`, they are
- # not there.
- j = j[j < work.active.shape[0]]
- # Check whether they are still there.
+ # not there. This happens when a value in also_stop is larger than
+ # the greatest value in work.active. This case needs special handling
+ # because we cannot simply check that also_stop == work.active[j].
+ mask = j < work.active.shape[0]
+ # Note that we also have to use the mask to filter also_stop to ensure
+ # that also_stop and j will still have the same shape.
+ j, also_stop = j[mask], also_stop[mask]
j = j[also_stop == work.active[j]]
# Now convert these to boolean indices to use with `work.status`.
i = xp.zeros_like(stop)
@@ -407,6 +430,8 @@
if (not xp.isdtype(xm0.dtype, "numeric")
or xp.isdtype(xm0.dtype, "complex floating")):
raise ValueError('`xm0` must be numeric and real.')
+ if not xp.isdtype(xm0.dtype, "real floating"):
+ xm0 = xp.asarray(xm0, dtype=xp_default_dtype(xp))
xmin = -xp.inf if xmin is None else xmin
xmax = xp.inf if xmax is None else xmax
@@ -457,8 +482,10 @@
# of (xmin, xmax).
if xl0_not_supplied:
xl0 = xm0 - xp.minimum((xm0 - xmin)/16, xp.asarray(0.5))
+ xl0 = xp.astype(xl0, xm0.dtype, copy=False)
if xr0_not_supplied:
xr0 = xm0 + xp.minimum((xmax - xm0)/16, xp.asarray(0.5))
+ xr0 = xp.astype(xr0, xm0.dtype, copy=False)
maxiter = xp.asarray(maxiter)
message = '`maxiter` must be a non-negative integer.'
diff -Nru scipy-1.15.2/scipy/optimize/_highspy/_highs_wrapper.py scipy-1.15.3/scipy/optimize/_highspy/_highs_wrapper.py
--- scipy-1.15.2/scipy/optimize/_highspy/_highs_wrapper.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/optimize/_highspy/_highs_wrapper.py 2025-05-08 04:14:18.000000000 +0200
@@ -264,11 +264,13 @@
# Lagrangians for bounds based on column statuses
marg_bnds = np.zeros((2, numcol))
+ basis_col_status = basis.col_status
+ solution_col_dual = solution.col_dual
for ii in range(numcol):
- if basis.col_status[ii] == _h.HighsBasisStatus.kLower:
- marg_bnds[0, ii] = solution.col_dual[ii]
- elif basis.col_status[ii] == _h.HighsBasisStatus.kUpper:
- marg_bnds[1, ii] = solution.col_dual[ii]
+ if basis_col_status[ii] == _h.HighsBasisStatus.kLower:
+ marg_bnds[0, ii] = solution_col_dual[ii]
+ elif basis_col_status[ii] == _h.HighsBasisStatus.kUpper:
+ marg_bnds[1, ii] = solution_col_dual[ii]
res.update(
{
diff -Nru scipy-1.15.2/scipy/optimize/_linprog_highs.py scipy-1.15.3/scipy/optimize/_linprog_highs.py
--- scipy-1.15.2/scipy/optimize/_linprog_highs.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/optimize/_linprog_highs.py 2025-05-08 04:14:18.000000000 +0200
@@ -23,13 +23,13 @@
HighsDebugLevel,
ObjSense,
HighsModelStatus,
-)
-from ._highspy._core.simplex_constants import (
- SimplexStrategy,
- SimplexEdgeWeightStrategy,
+ simplex_constants as s_c, # [1]
)
from scipy.sparse import csc_matrix, vstack, issparse
+# [1]: Directly importing from "._highspy._core.simplex_constants"
+# causes problems when reloading.
+# See https://github.com/scipy/scipy/pull/22869 for details.
def _highs_to_scipy_status_message(highs_status, highs_message):
"""Converts HiGHS status number/message to SciPy status number/message"""
@@ -293,13 +293,13 @@
simplex_dual_edge_weight_strategy,
'simplex_dual_edge_weight_strategy',
choices={'dantzig': \
- SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategyDantzig,
+ s_c.SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategyDantzig,
'devex': \
- SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategyDevex,
+ s_c.SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategyDevex,
'steepest-devex': \
- SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategyChoose,
+ s_c.SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategyChoose,
'steepest': \
- SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategySteepestEdge,
+ s_c.SimplexEdgeWeightStrategy.kSimplexEdgeWeightStrategySteepestEdge,
None: None})
c, A_ub, b_ub, A_eq, b_eq, bounds, x0, integrality = lp
@@ -334,7 +334,7 @@
'primal_feasibility_tolerance': primal_feasibility_tolerance,
'simplex_dual_edge_weight_strategy':
simplex_dual_edge_weight_strategy_enum,
- 'simplex_strategy': SimplexStrategy.kSimplexStrategyDual,
+ 'simplex_strategy': s_c.SimplexStrategy.kSimplexStrategyDual,
'ipm_iteration_limit': maxiter,
'simplex_iteration_limit': maxiter,
'mip_rel_gap': mip_rel_gap,
diff -Nru scipy-1.15.2/scipy/optimize/_linprog_util.py scipy-1.15.3/scipy/optimize/_linprog_util.py
--- scipy-1.15.2/scipy/optimize/_linprog_util.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/optimize/_linprog_util.py 2025-05-08 04:14:18.000000000 +0200
@@ -1409,9 +1409,10 @@
x = rev(x)
fun = x.dot(c)
- slack = b_ub - A_ub.dot(x) # report slack for ORIGINAL UB constraints
- # report residuals of ORIGINAL EQ constraints
- con = b_eq - A_eq.dot(x)
+ with np.errstate(invalid="ignore"):
+ slack = b_ub - A_ub.dot(x) # report slack for ORIGINAL UB constraints
+ # report residuals of ORIGINAL EQ constraints
+ con = b_eq - A_eq.dot(x)
return x, fun, slack, con
diff -Nru scipy-1.15.2/scipy/optimize/tests/test_bracket.py scipy-1.15.3/scipy/optimize/tests/test_bracket.py
--- scipy-1.15.2/scipy/optimize/tests/test_bracket.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/optimize/tests/test_bracket.py 2025-05-08 04:14:18.000000000 +0200
@@ -352,6 +352,41 @@
xmin=1)
assert not res.success
+ def test_bug_fixes(self):
+ # 1. Bug in double sided bracket search.
+ # Happened in some cases where there are terminations on one side
+ # after corresponding searches on other side failed due to reaching the
+ # boundary.
+
+ # https://github.com/scipy/scipy/pull/22560#discussion_r1962853839
+ def f(x, p):
+ return np.exp(x) - p
+
+ p = np.asarray([0.29, 0.35])
+ res = _bracket_root(f, xl0=-1, xmin=-np.inf, xmax=0, args=(p, ))
+
+ # https://github.com/scipy/scipy/pull/22560/files#r1962952517
+ def f(x, p, c):
+ return np.exp(x*c) - p
+
+ p = [0.32061201, 0.39175242, 0.40047535, 0.50527218, 0.55654373,
+ 0.11911647, 0.37507896, 0.66554191]
+ c = [1., -1., 1., 1., -1., 1., 1., 1.]
+ xl0 = [-7.63108551, 3.27840947, -8.36968526, -1.78124372,
+ 0.92201295, -2.48930123, -0.66733533, -0.44606749]
+ xr0 = [-6.63108551, 4.27840947, -7.36968526, -0.78124372,
+ 1.92201295, -1.48930123, 0., 0.]
+ xmin = [-np.inf, 0., -np.inf, -np.inf, 0., -np.inf, -np.inf,
+ -np.inf]
+ xmax = [0., np.inf, 0., 0., np.inf, 0., 0., 0.]
+
+ res = _bracket_root(f, xl0=xl0, xr0=xr0, xmin=xmin, xmax=xmax, args=(p, c))
+
+ # 2. Default xl0 + 1 for xr0 exceeds xmax.
+ # https://github.com/scipy/scipy/pull/22560#discussion_r1962947434
+ res = _bracket_root(lambda x: x + 0.25, xl0=-0.5, xmin=-np.inf, xmax=0)
+ assert res.success
+
@pytest.mark.skip_xp_backends('array_api_strict', reason=array_api_strict_skip_reason)
@pytest.mark.skip_xp_backends('jax.numpy', reason=jax_skip_reason)
diff -Nru scipy-1.15.2/scipy/signal/_short_time_fft.py scipy-1.15.3/scipy/signal/_short_time_fft.py
--- scipy-1.15.2/scipy/signal/_short_time_fft.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/signal/_short_time_fft.py 2025-05-08 04:14:18.000000000 +0200
@@ -206,14 +206,14 @@
It is possible to calculate the SFT of signal parts:
- >>> p_q = SFT.nearest_k_p(N // 2)
- >>> Sx0 = SFT.stft(x[:p_q])
- >>> Sx1 = SFT.stft(x[p_q:])
+ >>> N2 = SFT.nearest_k_p(N // 2)
+ >>> Sx0 = SFT.stft(x[:N2])
+ >>> Sx1 = SFT.stft(x[N2:])
When assembling sequential STFT parts together, the overlap needs to be
considered:
- >>> p0_ub = SFT.upper_border_begin(p_q)[1] - SFT.p_min
+ >>> p0_ub = SFT.upper_border_begin(N2)[1] - SFT.p_min
>>> p1_le = SFT.lower_border_end[1] - SFT.p_min
>>> Sx01 = np.hstack((Sx0[:, :p0_ub],
... Sx0[:, p0_ub:] + Sx1[:, :p1_le],
@@ -1226,7 +1226,15 @@
@lru_cache(maxsize=256)
def _post_padding(self, n: int) -> tuple[int, int]:
- """Largest signal index and slice index due to padding."""
+ """Largest signal index and slice index due to padding.
+
+ Parameters
+ ----------
+ n : int
+ Number of samples of input signal (must be ≥ half of the window length).
+ """
+ if not (n >= (m2p := self.m_num - self.m_num_mid)):
+ raise ValueError(f"Parameter n must be >= ceil(m_num/2) = {m2p}!")
w2 = self.win.real**2 + self.win.imag**2
# move window to the right until the overlap for t < t[n] vanishes:
q1 = n // self.hop # last slice index with t[p1] <= t[n]
@@ -1247,6 +1255,11 @@
A detailed example is provided in the :ref:`tutorial_stft_sliding_win`
section of the :ref:`user_guide`.
+ Parameters
+ ----------
+ n : int
+ Number of samples of input signal (must be ≥ half of the window length).
+
See Also
--------
k_min: The smallest possible signal index.
@@ -1348,6 +1361,19 @@
A detailed example is given :ref:`tutorial_stft_sliding_win` section
of the :ref:`user_guide`.
+ Parameters
+ ----------
+ n : int
+ Number of samples of input signal (must be ≥ half of the window length).
+
+ Returns
+ -------
+ k_ub : int
+ Lowest signal index, where a touching time slice sticks out past the
+ signal end.
+ p_ub : int
+ Lowest index of time slice of which the end sticks out past the signal end.
+
See Also
--------
k_min: The smallest possible signal index.
@@ -1359,13 +1385,15 @@
p_range: Determine and validate slice index range.
ShortTimeFFT: Class this method belongs to.
"""
+ if not (n >= (m2p := self.m_num - self.m_num_mid)):
+ raise ValueError(f"Parameter n must be >= ceil(m_num/2) = {m2p}!")
w2 = self.win.real**2 + self.win.imag**2
q2 = n // self.hop + 1 # first t[q] >= t[n]
q1 = max((n-self.m_num) // self.hop - 1, -1)
# move window left until does not stick out to the right:
for q_ in range(q2, q1, -1):
k_ = q_ * self.hop + (self.m_num - self.m_num_mid)
- if k_ < n or all(w2[n-k_:] == 0):
+ if k_ <= n or all(w2[n-k_:] == 0):
return (q_ + 1) * self.hop - self.m_num_mid, q_ + 1
return 0, 0 # border starts at first slice
diff -Nru scipy-1.15.2/scipy/signal/_signaltools.py scipy-1.15.3/scipy/signal/_signaltools.py
--- scipy-1.15.2/scipy/signal/_signaltools.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/signal/_signaltools.py 2025-05-08 04:14:18.000000000 +0200
@@ -3703,7 +3703,10 @@
max_rate = max(up, down)
f_c = 1. / max_rate # cutoff of FIR filter (rel. to Nyquist)
half_len = 10 * max_rate # reasonable cutoff for sinc-like function
- if np.issubdtype(x.dtype, np.floating):
+ if np.issubdtype(x.dtype, np.complexfloating):
+ h = firwin(2 * half_len + 1, f_c,
+ window=window).astype(x.dtype) # match dtype of x
+ elif np.issubdtype(x.dtype, np.floating):
h = firwin(2 * half_len + 1, f_c,
window=window).astype(x.dtype) # match dtype of x
else:
diff -Nru scipy-1.15.2/scipy/signal/tests/_scipy_spectral_test_shim.py scipy-1.15.3/scipy/signal/tests/_scipy_spectral_test_shim.py
--- scipy-1.15.2/scipy/signal/tests/_scipy_spectral_test_shim.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/signal/tests/_scipy_spectral_test_shim.py 2025-05-08 04:14:18.000000000 +0200
@@ -103,7 +103,7 @@
# This is an edge case where shortTimeFFT returns one more time slice
# than the Scipy stft() shorten to remove last time slice:
if n % 2 == 1 and nperseg % 2 == 1 and noverlap % 2 == 1:
- x = x[..., :axis - 1]
+ x = x[..., : -1]
nadd = (-(x.shape[-1]-nperseg) % nstep) % nperseg
zeros_shape = list(x.shape[:-1]) + [nadd]
@@ -124,11 +124,8 @@
k_off = nperseg // 2
p0 = 0 # ST.lower_border_end[1] + 1
nn = x.shape[axis] if padded else n+k_off+1
- p1 = ST.upper_border_begin(nn)[1] # ST.p_max(n) + 1
-
- # This is bad hack to pass the test test_roundtrip_boundary_extension():
- if padded is True and nperseg - noverlap == 1:
- p1 -= nperseg // 2 - 1 # the reasoning behind this is not clear to me
+ # number of frames akin to legacy stft computation
+ p1 = (x.shape[axis] - nperseg) // nstep + 1
detr = None if detrend is False else detrend
Sxx = ST.stft_detrend(x, detr, p0, p1, k_offset=k_off, axis=axis)
@@ -136,11 +133,6 @@
if x.dtype in (np.float32, np.complex64):
Sxx = Sxx.astype(np.complex64)
- # workaround for test_average_all_segments() - seems to be buggy behavior:
- if boundary is None and padded is False:
- t, Sxx = t[1:-1], Sxx[..., :-2]
- t -= k_off / fs
-
return ST.f, t, Sxx
diff -Nru scipy-1.15.2/scipy/signal/tests/test_short_time_fft.py scipy-1.15.3/scipy/signal/tests/test_short_time_fft.py
--- scipy-1.15.2/scipy/signal/tests/test_short_time_fft.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/signal/tests/test_short_time_fft.py 2025-05-08 04:14:18.000000000 +0200
@@ -332,7 +332,11 @@
assert SFT.p_max(10) == 4
assert SFT.k_max(10) == 16
assert SFT.upper_border_begin(10) == (4, 2)
-
+ # Raise exceptions:
+ with pytest.raises(ValueError, match="^Parameter n must be"):
+ SFT.upper_border_begin(3)
+ with pytest.raises(ValueError, match="^Parameter n must be"):
+ SFT._post_padding(3)
def test_border_values_exotic():
"""Ensure that the border calculations are correct for windows with
@@ -342,7 +346,11 @@
assert SFT.lower_border_end == (0, 0)
SFT = ShortTimeFFT(np.flip(w), hop=20, fs=1)
- assert SFT.upper_border_begin(4) == (0, 0)
+ assert SFT.upper_border_begin(4) == (16, 1)
+ assert SFT.upper_border_begin(5) == (16, 1)
+ assert SFT.upper_border_begin(23) == (36, 2)
+ assert SFT.upper_border_begin(24) == (36, 2)
+ assert SFT.upper_border_begin(25) == (36, 2)
SFT._hop = -1 # provoke unreachable line
with pytest.raises(RuntimeError):
diff -Nru scipy-1.15.2/scipy/signal/tests/test_signaltools.py scipy-1.15.3/scipy/signal/tests/test_signaltools.py
--- scipy-1.15.2/scipy/signal/tests/test_signaltools.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/signal/tests/test_signaltools.py 2025-05-08 04:14:18.000000000 +0200
@@ -3926,3 +3926,8 @@
unique, multiplicity = unique_roots(p, 2)
assert_almost_equal(unique, [np.min(p)], decimal=15)
assert_equal(multiplicity, [100])
+
+
+def test_gh_22684():
+ actual = signal.resample_poly(np.arange(2000, dtype=np.complex64), 6, 4)
+ assert actual.dtype == np.complex64
diff -Nru scipy-1.15.2/scipy/sparse/_base.py scipy-1.15.3/scipy/sparse/_base.py
--- scipy-1.15.2/scipy/sparse/_base.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/_base.py 2025-05-08 04:14:18.000000000 +0200
@@ -1141,13 +1141,8 @@
if self.ndim == 1:
if axis not in (None, -1, 0):
raise ValueError("axis must be None, -1 or 0")
- ret = (self @ np.ones(self.shape, dtype=res_dtype)).astype(dtype)
-
- if out is not None:
- if any(dim != 1 for dim in out.shape):
- raise ValueError("dimensions do not match")
- out[...] = ret
- return ret
+ res = self @ np.ones(self.shape, dtype=res_dtype)
+ return res.sum(dtype=dtype, out=out)
# We use multiplication by a matrix of ones to achieve this.
# For some sparse array formats more efficient methods are
@@ -1175,14 +1170,6 @@
np.ones((N, 1), dtype=res_dtype)
)
- if out is not None:
- if isinstance(self, sparray):
- ret_shape = ret.shape[:axis] + ret.shape[axis + 1:]
- else:
- ret_shape = ret.shape
- if out.shape != ret_shape:
- raise ValueError("dimensions do not match")
-
return ret.sum(axis=axis, dtype=dtype, out=out)
def mean(self, axis=None, dtype=None, out=None):
diff -Nru scipy-1.15.2/scipy/sparse/_compressed.py scipy-1.15.3/scipy/sparse/_compressed.py
--- scipy-1.15.2/scipy/sparse/_compressed.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/_compressed.py 2025-05-08 04:14:18.000000000 +0200
@@ -677,9 +677,6 @@
if axis % 2 == 1:
ret = ret.T
- if out is not None and out.shape != ret.shape:
- raise ValueError('dimensions do not match')
-
return ret.sum(axis=(), dtype=dtype, out=out)
else:
# _spbase handles the situations when axis is in {None, -2, -1, 0, 1}
diff -Nru scipy-1.15.2/scipy/sparse/_dia.py scipy-1.15.3/scipy/sparse/_dia.py
--- scipy-1.15.2/scipy/sparse/_dia.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/_dia.py 2025-05-08 04:14:18.000000000 +0200
@@ -174,9 +174,6 @@
ret = self._ascontainer(row_sums.sum(axis=axis))
- if out is not None and out.shape != ret.shape:
- raise ValueError("dimensions do not match")
-
return ret.sum(axis=(), dtype=dtype, out=out)
sum.__doc__ = _spbase.sum.__doc__
diff -Nru scipy-1.15.2/scipy/sparse/linalg/_eigen/arpack/arpack.py scipy-1.15.3/scipy/sparse/linalg/_eigen/arpack/arpack.py
--- scipy-1.15.2/scipy/sparse/linalg/_eigen/arpack/arpack.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/_eigen/arpack/arpack.py 2025-05-08 04:14:18.000000000 +0200
@@ -46,7 +46,6 @@
from scipy.sparse.linalg import gmres, splu
from scipy._lib._util import _aligned_zeros
from scipy._lib._threadsafety import ReentrancyLock
-
from . import _arpack
arpack_int = _arpack.timing.nbx.dtype
@@ -277,9 +276,12 @@
ARPACK error
"""
- def __init__(self, info, infodict=_NAUPD_ERRORS):
+ def __init__(self, info, infodict=None):
+ if infodict is None:
+ infodict = _NAUPD_ERRORS
+
msg = infodict.get(info, "Unknown error")
- RuntimeError.__init__(self, "ARPACK error %d: %s" % (info, msg))
+ super().__init__(f"ARPACK error {info}: {msg}")
class ArpackNoConvergence(ArpackError):
diff -Nru scipy-1.15.2/scipy/sparse/linalg/_expm_multiply.py scipy-1.15.3/scipy/sparse/linalg/_expm_multiply.py
--- scipy-1.15.2/scipy/sparse/linalg/_expm_multiply.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/_expm_multiply.py 2025-05-08 04:14:18.000000000 +0200
@@ -115,7 +115,7 @@
----------
A : transposable linear operator
The operator whose exponential is of interest.
- B : ndarray
+ B : ndarray, sparse array
The matrix or vector to be multiplied by the matrix exponential of A.
start : scalar, optional
The starting time point of the sequence.
@@ -443,7 +443,7 @@
def d(self, p):
"""
- Lazily estimate :math:`d_p(A) ~= || A^p ||^(1/p)`
+ Lazily estimate :math:`d_p(A) ~= || A^p ||^(1/p)`
where :math:`||.||` is the 1-norm.
"""
if p not in self._d:
@@ -702,7 +702,12 @@
m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell)
# Compute the expm action up to the initial time point.
- X[0] = _expm_multiply_simple_core(A, B, t_0, mu, m_star, s)
+ action_t0 = _expm_multiply_simple_core(A, B, t_0, mu, m_star, s)
+ if scipy.sparse.issparse(action_t0):
+ action_t0 = action_t0.toarray()
+ elif is_pydata_spmatrix(action_t0):
+ action_t0 = action_t0.todense()
+ X[0] = action_t0
# Compute the expm action at the rest of the time points.
if q <= s:
diff -Nru scipy-1.15.2/scipy/sparse/linalg/_interface.py scipy-1.15.3/scipy/sparse/linalg/_interface.py
--- scipy-1.15.2/scipy/sparse/linalg/_interface.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/_interface.py 2025-05-08 04:14:18.000000000 +0200
@@ -322,6 +322,10 @@
"""Default implementation of _rmatvec; defers to adjoint."""
if type(self)._adjoint == LinearOperator._adjoint:
# _adjoint not overridden, prevent infinite recursion
+ if (hasattr(self, "_rmatmat")
+ and type(self)._rmatmat != LinearOperator._rmatmat):
+ # Try to use _rmatmat as a fallback
+ return self._rmatmat(x.reshape(-1, 1)).reshape(-1)
raise NotImplementedError
else:
return self.H.matvec(x)
@@ -822,22 +826,22 @@
def _adjoint(self):
if self.__adj is None:
- self.__adj = _AdjointMatrixOperator(self)
+ self.__adj = _AdjointMatrixOperator(self.A)
return self.__adj
+
class _AdjointMatrixOperator(MatrixLinearOperator):
- def __init__(self, adjoint):
- self.A = adjoint.A.T.conj()
- self.__adjoint = adjoint
- self.args = (adjoint,)
- self.shape = adjoint.shape[1], adjoint.shape[0]
+ def __init__(self, adjoint_array):
+ self.A = adjoint_array.T.conj()
+ self.args = (adjoint_array,)
+ self.shape = adjoint_array.shape[1], adjoint_array.shape[0]
@property
def dtype(self):
- return self.__adjoint.dtype
+ return self.args[0].dtype
def _adjoint(self):
- return self.__adjoint
+ return MatrixLinearOperator(self.args[0])
class IdentityOperator(LinearOperator):
diff -Nru scipy-1.15.2/scipy/sparse/linalg/_isolve/_gcrotmk.py scipy-1.15.3/scipy/sparse/linalg/_isolve/_gcrotmk.py
--- scipy-1.15.2/scipy/sparse/linalg/_isolve/_gcrotmk.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/_isolve/_gcrotmk.py 2025-05-08 04:14:18.000000000 +0200
@@ -433,7 +433,8 @@
ux = axpy(u, ux, ux.shape[0], -byc) # ux -= u*byc
# cx := V H y
- hy = Q.dot(R.dot(y))
+ with np.errstate(invalid="ignore"):
+ hy = Q.dot(R.dot(y))
cx = vs[0] * hy[0]
for v, hyc in zip(vs[1:], hy[1:]):
cx = axpy(v, cx, cx.shape[0], hyc) # cx += v*hyc
diff -Nru scipy-1.15.2/scipy/sparse/linalg/tests/test_expm_multiply.py scipy-1.15.3/scipy/sparse/linalg/tests/test_expm_multiply.py
--- scipy-1.15.2/scipy/sparse/linalg/tests/test_expm_multiply.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/tests/test_expm_multiply.py 2025-05-08 04:14:18.000000000 +0200
@@ -7,6 +7,7 @@
from numpy.testing import (assert_allclose, assert_, assert_equal,
suppress_warnings)
from scipy.sparse import SparseEfficiencyWarning
+import scipy.sparse
from scipy.sparse.linalg import aslinearoperator
import scipy.linalg
from scipy.sparse.linalg import expm as sp_expm
@@ -260,19 +261,28 @@
A = scipy.sparse.diags_array(np.arange(5),format='csr', dtype=int)
B = np.ones(5, dtype=int)
Aexpm = scipy.sparse.diags_array(np.exp(np.arange(5)),format='csr')
+ BI = np.identity(5, dtype=int)
+ BI_sparse = scipy.sparse.csr_array(BI)
assert_allclose(expm_multiply(A,B,0,1)[-1], Aexpm.dot(B))
+ assert_allclose(np.diag(expm_multiply(A, BI_sparse, 0, 1)[-1]), Aexpm.dot(B))
# Test A complex, B int
A = scipy.sparse.diags_array(-1j*np.arange(5),format='csr', dtype=complex)
B = np.ones(5, dtype=int)
Aexpm = scipy.sparse.diags_array(np.exp(-1j*np.arange(5)),format='csr')
assert_allclose(expm_multiply(A,B,0,1)[-1], Aexpm.dot(B))
+ assert_allclose(np.diag(expm_multiply(A, BI_sparse, 0, 1)[-1]), Aexpm.dot(B))
# Test A int, B complex
A = scipy.sparse.diags_array(np.arange(5),format='csr', dtype=int)
B = np.full(5, 1j, dtype=complex)
Aexpm = scipy.sparse.diags_array(np.exp(np.arange(5)),format='csr')
assert_allclose(expm_multiply(A,B,0,1)[-1], Aexpm.dot(B))
+ BI = np.identity(5, dtype=complex)*1j
+ assert_allclose(
+ np.diag(expm_multiply(A, scipy.sparse.csr_array(BI), 0, 1)[-1]),
+ Aexpm.dot(B)
+ )
def test_expm_multiply_interval_status_0(self):
self._help_test_specific_expm_interval_status(0)
diff -Nru scipy-1.15.2/scipy/sparse/linalg/tests/test_interface.py scipy-1.15.3/scipy/sparse/linalg/tests/test_interface.py
--- scipy-1.15.2/scipy/sparse/linalg/tests/test_interface.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/tests/test_interface.py 2025-05-08 04:14:18.000000000 +0200
@@ -13,6 +13,7 @@
import scipy.sparse.linalg._interface as interface
from scipy.sparse._sputils import matrix
+from scipy._lib._gcutils import assert_deallocated, IS_PYPY
class TestLinearOperator:
@@ -512,6 +513,30 @@
assert_equal(B.dot(v), Y.dot(v))
assert_equal(B.T.dot(v), Y.T.dot(v))
+def test_transpose_multiplication():
+ class MyMatrix(interface.LinearOperator):
+ def __init__(self, A):
+ super().__init__(A.dtype, A.shape)
+ self.A = A
+ def _matmat(self, other): return self.A @ other
+ def _rmatmat(self, other): return self.A.T @ other
+
+ A = MyMatrix(np.array([[1, 2], [3, 4]]))
+ X = np.array([1, 2])
+ B = np.array([[10, 20], [30, 40]])
+ X2 = X.reshape(-1, 1)
+ Y = np.array([[1, 2], [3, 4]])
+
+ assert_equal(A @ B, Y @ B)
+ assert_equal(B.T @ A, B.T @ Y)
+ assert_equal(A.T @ B, Y.T @ B)
+ assert_equal(A @ X, Y @ X)
+ assert_equal(X.T @ A, X.T @ Y)
+ assert_equal(A.T @ X, Y.T @ X)
+ assert_equal(A @ X2, Y @ X2)
+ assert_equal(X2.T @ A, X2.T @ Y)
+ assert_equal(A.T @ X2, Y.T @ X2)
+
def test_sparse_matmat_exception():
A = interface.LinearOperator((2, 2), matvec=lambda x: x)
B = sparse.eye_array(2)
@@ -524,3 +549,13 @@
A @ np.identity(4)
with assert_raises(ValueError):
np.identity(4) @ A
+
+
+@pytest.mark.skipif(IS_PYPY, reason="Test not meaningful on PyPy")
+def test_MatrixLinearOperator_refcycle():
+ # gh-10634
+ # Test that MatrixLinearOperator can be automatically garbage collected
+ A = np.eye(2)
+ with assert_deallocated(interface.MatrixLinearOperator, A) as op:
+ op.adjoint()
+ del op
diff -Nru scipy-1.15.2/scipy/sparse/linalg/tests/test_pydata_sparse.py scipy-1.15.3/scipy/sparse/linalg/tests/test_pydata_sparse.py
--- scipy-1.15.2/scipy/sparse/linalg/tests/test_pydata_sparse.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/linalg/tests/test_pydata_sparse.py 2025-05-08 04:14:18.000000000 +0200
@@ -239,6 +239,10 @@
x = splin.expm_multiply(A_sparse, b)
assert_allclose(x, x0)
+ x0 = splin.expm_multiply(A_dense, A_dense)
+ x = splin.expm_multiply(A_sparse, A_sparse)
+ assert_allclose(x.todense(), x0)
+
def test_eq(same_matrix):
sp_sparse, pd_sparse = same_matrix
diff -Nru scipy-1.15.2/scipy/sparse/tests/test_base.py scipy-1.15.3/scipy/sparse/tests/test_base.py
--- scipy-1.15.2/scipy/sparse/tests/test_base.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/tests/test_base.py 2025-05-08 04:14:18.000000000 +0200
@@ -1087,6 +1087,12 @@
datsp.sum(axis=1, out=datsp_out)
assert_array_almost_equal(dat_out, datsp_out)
+ # check that wrong shape out parameter raises
+ with assert_raises(ValueError, match="output parameter.*wrong.*dimension"):
+ datsp.sum(out=array([0]))
+ with assert_raises(ValueError, match="output parameter.*wrong.*dimension"):
+ datsp.sum(out=array([[0]] if self.is_array_test else 0))
+
def test_numpy_sum(self):
# See gh-5987
dat = array([[0, 1, 2],
@@ -1184,6 +1190,12 @@
datsp.mean(axis=1, out=datsp_out)
assert_array_almost_equal(dat_out, datsp_out)
+ # check that wrong shape out parameter raises
+ with assert_raises(ValueError, match="output parameter.*wrong.*dimension"):
+ datsp.mean(out=array([0]))
+ with assert_raises(ValueError, match="output parameter.*wrong.*dimension"):
+ datsp.mean(out=array([[0]] if self.is_array_test else 0))
+
def test_numpy_mean(self):
# See gh-5987
dat = array([[0, 1, 2],
diff -Nru scipy-1.15.2/scipy/sparse/tests/test_common1d.py scipy-1.15.3/scipy/sparse/tests/test_common1d.py
--- scipy-1.15.2/scipy/sparse/tests/test_common1d.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/sparse/tests/test_common1d.py 2025-05-08 04:14:18.000000000 +0200
@@ -142,7 +142,7 @@
datsp.sum(axis=(0, 1))
with pytest.raises(TypeError, match='axis must be an integer'):
datsp.sum(axis=1.5)
- with pytest.raises(ValueError, match='dimensions do not match'):
+ with pytest.raises(ValueError, match='output parameter.*wrong.*dimension'):
datsp.sum(axis=0, out=out)
def test_numpy_sum(self, spcreator):
@@ -180,7 +180,7 @@
datsp.mean(axis=(0, 1))
with pytest.raises(TypeError, match='axis must be an integer'):
datsp.mean(axis=1.5)
- with pytest.raises(ValueError, match='dimensions do not match'):
+ with pytest.raises(ValueError, match='output parameter.*wrong.*dimension'):
datsp.mean(axis=1, out=out)
def test_sum_dtype(self, spcreator):
@@ -209,17 +209,22 @@
dat = np.array([0, 1, 2])
datsp = spcreator(dat)
- dat_out = np.array([0])
- datsp_out = np.array([0])
+ dat_out = np.array(0)
+ datsp_out = np.array(0)
- dat.mean(out=dat_out, keepdims=True)
+ dat.mean(out=dat_out)
datsp.mean(out=datsp_out)
assert_allclose(dat_out, datsp_out)
- dat.mean(axis=0, out=dat_out, keepdims=True)
+ dat.mean(axis=0, out=dat_out)
datsp.mean(axis=0, out=datsp_out)
assert_allclose(dat_out, datsp_out)
+ with pytest.raises(ValueError, match="output parameter.*dimension"):
+ datsp.mean(out=np.array([0]))
+ with pytest.raises(ValueError, match="output parameter.*dimension"):
+ datsp.mean(out=np.array([[0]]))
+
def test_numpy_mean(self, spcreator):
dat = np.array([0, 1, 2])
datsp = spcreator(dat)
diff -Nru scipy-1.15.2/scipy/spatial/_qhull.pyx scipy-1.15.3/scipy/spatial/_qhull.pyx
--- scipy-1.15.2/scipy/spatial/_qhull.pyx 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/spatial/_qhull.pyx 2025-05-08 04:14:18.000000000 +0200
@@ -453,7 +453,7 @@
#Store the halfspaces in _points and the dual points in _dual_points later
self._point_arrays.append(np.array(points, copy=True))
dists = points[:, :-1].dot(interior_point)+points[:, -1]
- arr = np.array(-points[:, :-1]/dists, dtype=np.double, order="C", copy=True)
+ arr = np.array(-points[:, :-1]/dists[:, np.newaxis], dtype=np.double, order="C", copy=True)
else:
arr = np.array(points, dtype=np.double, order="C", copy=True)
@@ -2905,9 +2905,11 @@
Parameters
----------
- halfspaces : ndarray
- New halfspaces to add. The dimensionality should match that of the
- initial halfspaces.
+ halfspaces : ndarray of double, shape (n_new_ineq, ndim+1)
+ New halfspaces to add. The dimensionality (ndim) should match that of the
+ initial halfspaces. Like in the constructor, these are stacked
+ inequalites of the form Ax + b <= 0 in format [A; b]. The original
+ feasible point must also be feasible for these new inequalities.
restart : bool, optional
Whether to restart processing from scratch, rather than
adding halfspaces incrementally.
@@ -2929,6 +2931,22 @@
of halfspaces is also not possible after `close` has been called.
"""
+ if halfspaces.ndim > 2:
+ raise ValueError("`halfspaces` should be provided as a 2D array")
+ # We check for non-feasibility of incremental additions
+ # in a manner similar to `qh_sethalfspace`
+ halfspaces = np.atleast_2d(halfspaces)
+ dists = np.dot(halfspaces[:, :self.ndim], self.interior_point) + halfspaces[:, -1]
+ # HalfspaceIntersection uses closed half spaces so
+ # the feasible point also cannot be directly on the boundary
+ viols = dists >= 0
+ if viols.any():
+ # error out with an indication of the first violating
+ # half space discovered
+ first_viol = np.nonzero(viols)[0].min()
+ bad_hs = halfspaces[first_viol, :]
+ msg = f"feasible point is not clearly inside halfspace: {bad_hs}"
+ raise QhullError(msg)
self._add_points(halfspaces, restart, self.interior_point)
@property
diff -Nru scipy-1.15.2/scipy/spatial/tests/test_qhull.py scipy-1.15.3/scipy/spatial/tests/test_qhull.py
--- scipy-1.15.2/scipy/spatial/tests/test_qhull.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/spatial/tests/test_qhull.py 2025-05-08 04:14:18.000000000 +0200
@@ -1181,6 +1181,105 @@
assert_allclose(hs.dual_points, qhalf_points)
+ @pytest.mark.parametrize("k", range(1,4))
+ def test_halfspace_batch(self, k):
+ # Test that we can add halfspaces a few at a time
+ big_square = np.array([[ 1., 0., -2.],
+ [-1., 0., -2.],
+ [ 0., 1., -2.],
+ [ 0., -1., -2.]])
+
+ small_square = np.array([[ 1., 0., -1.],
+ [-1., 0., -1.],
+ [ 0., 1., -1.],
+ [ 0., -1., -1.]])
+
+ hs = qhull.HalfspaceIntersection(big_square,
+ np.array([0.3141, 0.2718]),
+ incremental=True)
+
+ hs.add_halfspaces(small_square[0:k,:])
+ hs.add_halfspaces(small_square[k:4,:])
+ hs.close()
+
+ # Check the intersections are correct (they are the corners of the small square)
+ expected_intersections = np.array([[1., 1.],
+ [1., -1.],
+ [-1., 1.],
+ [-1., -1.]])
+ actual_intersections = hs.intersections
+ # They may be in any order, so just check that under some permutation
+ # expected=actual.
+
+ ind1 = np.lexsort((actual_intersections[:, 1], actual_intersections[:, 0]))
+ ind2 = np.lexsort((expected_intersections[:, 1], expected_intersections[:, 0]))
+ assert_allclose(actual_intersections[ind1], expected_intersections[ind2])
+
+
+ @pytest.mark.parametrize("halfspaces", [
+ (np.array([-0.70613882, -0.45589431, 0.04178256])),
+ (np.array([[-0.70613882, -0.45589431, 0.04178256],
+ [0.70807342, -0.45464871, -0.45969769],
+ [0., 0.76515026, -0.35614825]])),
+ ])
+ def test_gh_19865(self, halfspaces):
+ # starting off with a feasible interior point and
+ # adding halfspaces for which it is no longer feasible
+ # should result in an error rather than a problematic
+ # intersection polytope
+ initial_square = np.array(
+ [[1, 0, -1], [0, 1, -1], [-1, 0, -1], [0, -1, -1]]
+ )
+ incremental_intersector = qhull.HalfspaceIntersection(initial_square,
+ np.zeros(2),
+ incremental=True)
+ with pytest.raises(qhull.QhullError, match="feasible.*-0.706.*"):
+ incremental_intersector.add_halfspaces(halfspaces)
+
+
+ def test_gh_19865_3d(self):
+ # 3d case where closed half space is enforced for
+ # feasibility
+ halfspaces = np.array([[1, 1, 1, -1], # doesn't exclude origin
+ [-1, -1, -1, -1], # doesn't exclude origin
+ [1, 0, 0, 0]]) # the origin is on the line
+ initial_cube = np.array([[1, 0, 0, -1],
+ [-1, 0, 0, -1],
+ [0, 1, 0, -1],
+ [0, -1, 0, -1],
+ [0, 0, 1, -1],
+ [0, 0, -1, -1]])
+ incremental_intersector = qhull.HalfspaceIntersection(initial_cube,
+ np.zeros(3),
+ incremental=True)
+ with pytest.raises(qhull.QhullError, match="feasible.*[1 0 0 0]"):
+ incremental_intersector.add_halfspaces(halfspaces)
+
+
+ def test_2d_add_halfspace_input(self):
+ # incrementally added halfspaces should respect the 2D
+ # array shape requirement
+ initial_square = np.array(
+ [[1, 0, -1], [0, 1, -1], [-1, 0, -1], [0, -1, -1]]
+ )
+ incremental_intersector = qhull.HalfspaceIntersection(initial_square,
+ np.zeros(2),
+ incremental=True)
+ with pytest.raises(ValueError, match="2D array"):
+ incremental_intersector.add_halfspaces(np.ones((4, 4, 4)))
+
+ def test_1d_add_halfspace_input(self):
+ # we do allow 1D `halfspaces` input to add_halfspaces()
+ initial_square = np.array(
+ [[1, 0, -1], [0, 1, -1], [-1, 0, -1], [0, -1, -1]]
+ )
+ incremental_intersector = qhull.HalfspaceIntersection(initial_square,
+ np.zeros(2),
+ incremental=True)
+ assert_allclose(incremental_intersector.dual_vertices, np.arange(4))
+ incremental_intersector.add_halfspaces(np.array([2, 2, -1]))
+ assert_allclose(incremental_intersector.dual_vertices, np.arange(5))
+
@pytest.mark.parametrize("diagram_type", [Voronoi, qhull.Delaunay])
def test_gh_20623(diagram_type):
diff -Nru scipy-1.15.2/scipy/spatial/transform/_rotation.pyx scipy-1.15.3/scipy/spatial/transform/_rotation.pyx
--- scipy-1.15.2/scipy/spatial/transform/_rotation.pyx 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/spatial/transform/_rotation.pyx 2025-05-08 04:14:18.000000000 +0200
@@ -15,12 +15,20 @@
# utilities for empty array initialization
cdef inline double[:] _empty1(int n) noexcept:
+ if n == 0:
+ return array(shape=(1,), itemsize=sizeof(double), format=b"d")[:0]
return array(shape=(n,), itemsize=sizeof(double), format=b"d")
cdef inline double[:, :] _empty2(int n1, int n2) noexcept :
+ if n1 == 0:
+ return array(shape=(1, n2), itemsize=sizeof(double), format=b"d")[:0]
return array(shape=(n1, n2), itemsize=sizeof(double), format=b"d")
cdef inline double[:, :, :] _empty3(int n1, int n2, int n3) noexcept:
+ if n1 == 0:
+ return array(shape=(1, n2, n3), itemsize=sizeof(double), format=b"d")[:0]
return array(shape=(n1, n2, n3), itemsize=sizeof(double), format=b"d")
cdef inline double[:, :] _zeros2(int n1, int n2) noexcept:
+ if n1 == 0:
+ return array(shape=(1, n2), itemsize=sizeof(double), format=b"d")[:0]
cdef double[:, :] arr = array(shape=(n1, n2),
itemsize=sizeof(double), format=b"d")
arr[:, :] = 0
@@ -448,7 +456,8 @@
cdef inline double[:, :] _compose_quat(
const double[:, :] p, const double[:, :] q
) noexcept:
- cdef Py_ssize_t n = max(p.shape[0], q.shape[0])
+ cdef Py_ssize_t n = q.shape[0] if p.shape[0] == 1 else p.shape[0]
+
cdef double[:, :] product = _empty2(n, 4)
# dealing with broadcasting
@@ -832,8 +841,7 @@
quat = np.asarray(quat, dtype=float)
if (quat.ndim not in [1, 2]
- or quat.shape[len(quat.shape) - 1] != 4
- or quat.shape[0] == 0):
+ or quat.shape[len(quat.shape) - 1] != 4):
raise ValueError("Expected `quat` to have shape (4,) or (N, 4), "
f"got {quat.shape}.")
@@ -2938,6 +2946,9 @@
>>> r.mean().as_euler('zyx', degrees=True)
array([0.24945696, 0.25054542, 0.24945696])
"""
+ if self._quat.shape[0] == 0:
+ raise ValueError("Mean of an empty rotation set is undefined.")
+
if weights is None:
weights = np.ones(len(self))
else:
@@ -3734,7 +3745,7 @@
if not isinstance(rotations, Rotation):
raise TypeError("`rotations` must be a `Rotation` instance.")
- if rotations.single or len(rotations) == 1:
+ if rotations.single or len(rotations) <= 1:
raise ValueError("`rotations` must be a sequence of at least 2 rotations.")
times = np.asarray(times)
diff -Nru scipy-1.15.2/scipy/spatial/transform/tests/test_rotation.py scipy-1.15.3/scipy/spatial/transform/tests/test_rotation.py
--- scipy-1.15.2/scipy/spatial/transform/tests/test_rotation.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/spatial/transform/tests/test_rotation.py 2025-05-08 04:14:18.000000000 +0200
@@ -168,10 +168,6 @@
[[4, 5, 6, 7]]
]))
- # 0-length 2d array
- with pytest.raises(ValueError, match='Expected `quat` to have shape'):
- Rotation.from_quat(np.array([]).reshape((0, 4)))
-
def test_zero_norms_from_quat():
x = np.array([
@@ -1611,18 +1607,23 @@
t = np.array([0, 1])
Slerp(t, r)
+SLERP_EXCEPTION_MESSAGE = "must be a sequence of at least 2 rotations"
def test_slerp_single_rot():
- msg = "must be a sequence of at least 2 rotations"
- with pytest.raises(ValueError, match=msg):
- r = Rotation.from_quat([1, 2, 3, 4])
+ r = Rotation.from_quat([1, 2, 3, 4])
+ with pytest.raises(ValueError, match=SLERP_EXCEPTION_MESSAGE):
Slerp([1], r)
+def test_slerp_rot_len0():
+ r = Rotation.random()
+ with pytest.raises(ValueError, match=SLERP_EXCEPTION_MESSAGE):
+ Slerp([], r)
+
+
def test_slerp_rot_len1():
- msg = "must be a sequence of at least 2 rotations"
- with pytest.raises(ValueError, match=msg):
- r = Rotation.from_quat([[1, 2, 3, 4]])
+ r = Rotation.random(1)
+ with pytest.raises(ValueError, match=SLERP_EXCEPTION_MESSAGE):
Slerp([1], r)
@@ -2015,3 +2016,173 @@
eul = rot.as_euler(seq)
dav = rot.as_davenport(ax, order)
assert_allclose(eul, dav, rtol=1e-12)
+
+
+def test_zero_rotation_construction():
+ r = Rotation.random(num=0)
+ assert len(r) == 0
+
+ r_ide = Rotation.identity(num=0)
+ assert len(r_ide) == 0
+
+ r_get = Rotation.random(num=3)[[]]
+ assert len(r_get) == 0
+
+ r_quat = Rotation.from_quat(np.zeros((0, 4)))
+ assert len(r_quat) == 0
+
+ r_matrix = Rotation.from_matrix(np.zeros((0, 3, 3)))
+ assert len(r_matrix) == 0
+
+ r_euler = Rotation.from_euler("xyz", np.zeros((0, 3)))
+ assert len(r_euler) == 0
+
+ r_vec = Rotation.from_rotvec(np.zeros((0, 3)))
+ assert len(r_vec) == 0
+
+ r_dav = Rotation.from_davenport(np.eye(3), "extrinsic", np.zeros((0, 3)))
+ assert len(r_dav) == 0
+
+ r_mrp = Rotation.from_mrp(np.zeros((0, 3)))
+ assert len(r_mrp) == 0
+
+
+def test_zero_rotation_representation():
+ r = Rotation.random(num=0)
+ assert r.as_quat().shape == (0, 4)
+ assert r.as_matrix().shape == (0, 3, 3)
+ assert r.as_euler("xyz").shape == (0, 3)
+ assert r.as_rotvec().shape == (0, 3)
+ assert r.as_mrp().shape == (0, 3)
+ assert r.as_davenport(np.eye(3), "extrinsic").shape == (0, 3)
+
+
+def test_zero_rotation_array_rotation():
+ r = Rotation.random(num=0)
+
+ v = np.array([1, 2, 3])
+ v_rotated = r.apply(v)
+ assert v_rotated.shape == (0, 3)
+
+ v0 = np.zeros((0, 3))
+ v0_rot = r.apply(v0)
+ assert v0_rot.shape == (0, 3)
+
+ v2 = np.ones((2, 3))
+ with pytest.raises(
+ ValueError, match="Expected equal numbers of rotations and vectors"):
+ r.apply(v2)
+
+
+def test_zero_rotation_multiplication():
+ r = Rotation.random(num=0)
+
+ r_single = Rotation.random()
+ r_mult_left = r * r_single
+ assert len(r_mult_left) == 0
+
+ r_mult_right = r_single * r
+ assert len(r_mult_right) == 0
+
+ r0 = Rotation.random(0)
+ r_mult = r * r0
+ assert len(r_mult) == 0
+
+ msg_rotation_error = "Expected equal number of rotations"
+ r2 = Rotation.random(2)
+ with pytest.raises(ValueError, match=msg_rotation_error):
+ r0 * r2
+
+ with pytest.raises(ValueError, match=msg_rotation_error):
+ r2 * r0
+
+
+def test_zero_rotation_concatentation():
+ r = Rotation.random(num=0)
+
+ r0 = Rotation.concatenate([r, r])
+ assert len(r0) == 0
+
+ r1 = r.concatenate([Rotation.random(), r])
+ assert len(r1) == 1
+
+ r3 = r.concatenate([Rotation.random(3), r])
+ assert len(r3) == 3
+
+ r4 = r.concatenate([r, Rotation.random(4)])
+ assert len(r4) == 4
+
+
+def test_zero_rotation_power():
+ r = Rotation.random(num=0)
+ for pp in [-1.5, -1, 0, 1, 1.5]:
+ pow0 = r**pp
+ assert len(pow0) == 0
+
+
+def test_zero_rotation_inverse():
+ r = Rotation.random(num=0)
+ r_inv = r.inv()
+ assert len(r_inv) == 0
+
+
+def test_zero_rotation_magnitude():
+ r = Rotation.random(num=0)
+ magnitude = r.magnitude()
+ assert magnitude.shape == (0,)
+
+
+def test_zero_rotation_mean():
+ r = Rotation.random(num=0)
+ with pytest.raises(ValueError, match="Mean of an empty rotation set is undefined."):
+ r.mean()
+
+
+def test_zero_rotation_approx_equal():
+ r = Rotation.random(0)
+ assert r.approx_equal(Rotation.random(0)).shape == (0,)
+ assert r.approx_equal(Rotation.random()).shape == (0,)
+ assert Rotation.random().approx_equal(r).shape == (0,)
+
+ approx_msg = "Expected equal number of rotations"
+ r3 = Rotation.random(2)
+ with pytest.raises(ValueError, match=approx_msg):
+ r.approx_equal(r3)
+
+ with pytest.raises(ValueError, match=approx_msg):
+ r3.approx_equal(r)
+
+
+def test_zero_rotation_get_set():
+ r = Rotation.random(0)
+
+ r_get = r[[]]
+ assert len(r_get) == 0
+
+ r_slice = r[:0]
+ assert len(r_slice) == 0
+
+ with pytest.raises(IndexError):
+ r[[0]]
+
+ with pytest.raises(IndexError):
+ r[[True]]
+
+ with pytest.raises(IndexError):
+ r[0] = Rotation.random()
+
+
+def test_boolean_indexes():
+ r = Rotation.random(3)
+
+ r0 = r[[False, False, False]]
+ assert len(r0) == 0
+
+ r1 = r[[False, True, False]]
+ assert len(r1) == 1
+
+ r3 = r[[True, True, True]]
+ assert len(r3) == 3
+
+ with pytest.raises(IndexError):
+ r[[True, True]]
diff -Nru scipy-1.15.2/scipy/special/_logsumexp.py scipy-1.15.3/scipy/special/_logsumexp.py
--- scipy-1.15.2/scipy/special/_logsumexp.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/special/_logsumexp.py 2025-05-08 04:14:18.000000000 +0200
@@ -5,9 +5,9 @@
array_namespace,
xp_size,
xp_broadcast_promote,
- xp_real,
xp_copy,
xp_float_to_complex,
+ is_complex,
)
from scipy._lib import array_api_extra as xpx
@@ -219,25 +219,30 @@
s = xp.where(s == 0, s, s/m)
# Separate sign/magnitude information
- sgn = None
- if return_sign:
- # Use the numpy>=2.0 convention for sign.
- # When all array libraries agree, this can become sng = xp.sign(s).
- sgn = _sign(s + 1, xp=xp) * _sign(m, xp=xp)
-
- if xp.isdtype(s.dtype, "real floating"):
- # The log functions need positive arguments
- s = xp.where(s < -1, -s - 2, s)
- m = xp.abs(m)
- else:
- # `a_max` can have a sign component for complex input
- j = xp.asarray(1j, dtype=a_max.dtype)
- sgn = sgn * xp.exp(xp.imag(a_max) * j)
+ # Originally, this was only performed if `return_sign=True`.
+ # However, this is also needed if any elements of `m < 0` or `s < -1`.
+ # An improvement would be to perform the calculations only on these entries.
+
+ # Use the numpy>=2.0 convention for sign.
+ # When all array libraries agree, this can become sng = xp.sign(s).
+ sgn = _sign(s + 1, xp=xp) * _sign(m, xp=xp)
+
+ if xp.isdtype(s.dtype, "real floating"):
+ # The log functions need positive arguments
+ s = xp.where(s < -1, -s - 2, s)
+ m = xp.abs(m)
+ else:
+ # `a_max` can have a sign component for complex input
+ sgn = sgn * xp.exp(xp.imag(a_max) * xp.asarray(1.0j, dtype=a_max.dtype))
# Take log and undo shift
out = xp.log1p(s) + xp.log(m) + a_max
- out = xp_real(out) if return_sign else out
+ if return_sign:
+ if is_complex(out, xp):
+ out = xp.real(out)
+ elif xp.isdtype(out.dtype, 'real floating'):
+ out[sgn < 0] = xp.nan
return out, sgn
diff -Nru scipy-1.15.2/scipy/special/tests/test_hyp2f1.py scipy-1.15.3/scipy/special/tests/test_hyp2f1.py
--- scipy-1.15.2/scipy/special/tests/test_hyp2f1.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/special/tests/test_hyp2f1.py 2025-05-08 04:14:18.000000000 +0200
@@ -2543,3 +2543,24 @@
if mark.name == 'parametrize'
for case in mark.args[1]
]
+
+class TestHyp2f1ExtremeInputs:
+
+ @pytest.mark.parametrize("a", [1.0, 2.0, 3.0, -np.inf, np.inf])
+ @pytest.mark.parametrize("b", [3.0, 4.0, 5.0, -np.inf, np.inf])
+ @pytest.mark.parametrize("c", [3.0, 5.0, 6.0, 7.0])
+ @pytest.mark.parametrize("z", [4.0 + 1.0j])
+ def test_inf_a_b(self, a, b, c, z):
+ if np.any(np.isinf(np.asarray([a, b]))):
+ assert(np.isnan(hyp2f1(a, b, c, z)))
+
+ def test_large_a_b(self):
+ assert(np.isnan(hyp2f1(10**7, 1.0, 3.0, 4.0 + 1.0j)))
+ assert(np.isnan(hyp2f1(-10**7, 1.0, 3.0, 4.0 + 1.0j)))
+
+ assert(np.isnan(hyp2f1(1.0, 10**7, 3.0, 4.0 + 1.0j)))
+ assert(np.isnan(hyp2f1(1.0, -10**7, 3.0, 4.0 + 1.0j)))
+
+ # Already correct in main but testing for surety
+ assert(np.isnan(hyp2f1(np.inf, 1.0, 3.0, 4.0)))
+ assert(np.isnan(hyp2f1(1.0, np.inf, 3.0, 4.0)))
diff -Nru scipy-1.15.2/scipy/special/tests/test_logsumexp.py scipy-1.15.3/scipy/special/tests/test_logsumexp.py
--- scipy-1.15.2/scipy/special/tests/test_logsumexp.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/special/tests/test_logsumexp.py 2025-05-08 04:14:18.000000000 +0200
@@ -257,6 +257,20 @@
xp_assert_close(xp.real(res), xp.real(ref))
xp_assert_close(xp.imag(res), xp.imag(ref), atol=0, rtol=1e-15)
+ def test_gh22903(self, xp):
+ # gh-22903 reported that `logsumexp` produced NaN where the weight associated
+ # with the max magnitude element was negative and `return_sign=False`, even if
+ # the net result should be the log of a positive number.
+
+ # result is log of positive number
+ a = xp.asarray([3.06409428, 0.37251854, 3.87471931])
+ b = xp.asarray([1.88190708, 2.84174795, -0.85016884])
+ xp_assert_close(logsumexp(a, b=b), logsumexp(a, b=b, return_sign=True)[0])
+
+ # result is log of negative number
+ b = xp.asarray([1.88190708, 2.84174795, -3.85016884])
+ xp_assert_close(logsumexp(a, b=b), xp.asarray(xp.nan))
+
class TestSoftmax:
def test_softmax_fixtures(self):
diff -Nru scipy-1.15.2/scipy/special/xsf/hyp2f1.h scipy-1.15.3/scipy/special/xsf/hyp2f1.h
--- scipy-1.15.2/scipy/special/xsf/hyp2f1.h 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/special/xsf/hyp2f1.h 2025-05-08 04:14:18.000000000 +0200
@@ -553,7 +553,7 @@
* the series at a or b of smaller magnitude. This is to ensure proper
* handling of situations like a < c < b <= 0, a, b, c all non-positive
* integers, where terminating at a would lead to a term of the form 0 / 0. */
- std::uint64_t max_degree;
+ double max_degree;
if (a_neg_int || b_neg_int) {
if (a_neg_int && b_neg_int) {
max_degree = a > b ? std::abs(a) : std::abs(b);
@@ -562,7 +562,7 @@
} else {
max_degree = std::abs(b);
}
- if (max_degree <= UINT64_MAX) {
+ if (max_degree <= (double) UINT64_MAX) {
auto series_generator = detail::HypergeometricSeriesGenerator(a, b, c, z);
return detail::series_eval_fixed_length(series_generator, std::complex<double>{0.0, 0.0}, max_degree + 1);
} else {
@@ -583,7 +583,7 @@
* (DLMF 15.8.1) */
if (c_minus_a_neg_int || c_minus_b_neg_int) {
max_degree = c_minus_b_neg_int ? std::abs(c - b) : std::abs(c - a);
- if (max_degree <= UINT64_MAX) {
+ if (max_degree <= (double) UINT64_MAX) {
result = std::pow(1.0 - z, c - a - b);
auto series_generator = detail::HypergeometricSeriesGenerator(c - a, c - b, c, z);
result *=
diff -Nru scipy-1.15.2/scipy/special/xsf/numpy.h scipy-1.15.3/scipy/special/xsf/numpy.h
--- scipy-1.15.2/scipy/special/xsf/numpy.h 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/special/xsf/numpy.h 2025-05-08 04:14:18.000000000 +0200
@@ -19,6 +19,9 @@
#include "dual.h"
#include "error.h"
+// Force defining the parenthesis operator even when compiling with a compiler
+// defaulting to C++ >= 23.
+#define MDSPAN_USE_PAREN_OPERATOR 1
#include "third_party/kokkos/mdspan.hpp"
/* PyUFunc_getfperr gets bits for current floating point error (fpe) status codes so we
diff -Nru scipy-1.15.2/scipy/special/xsf/specfun/specfun.h scipy-1.15.3/scipy/special/xsf/specfun/specfun.h
--- scipy-1.15.2/scipy/special/xsf/specfun/specfun.h 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/special/xsf/specfun/specfun.h 2025-05-08 04:14:18.000000000 +0200
@@ -887,7 +887,7 @@
// =======================================================
int il1, il2, il3, bl1, bl2, bl3, bn, id1 = 0, id;
- double aa, hu = 0.0, hu1;
+ double aa, hu = 0.0, hu1, b00;
aa = a - b + 1.0;
*isfer = 0;
@@ -927,10 +927,11 @@
}
} else {
if (b <= a) {
+ b00 = b;
a -= b - 1.0;
b = 2.0 - b;
hu = chguit(x, a, b, &id);
- hu *= pow(x, 1.0 - b);
+ hu *= pow(x, 1.0 - b00);
*md = 4;
} else if (bn && (~il1)) {
hu = chgubi(x, a, b, &id);
diff -Nru scipy-1.15.2/scipy/stats/_continuous_distns.py scipy-1.15.3/scipy/stats/_continuous_distns.py
--- scipy-1.15.2/scipy/stats/_continuous_distns.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/stats/_continuous_distns.py 2025-05-08 04:14:18.000000000 +0200
@@ -3432,29 +3432,39 @@
g3 = g(3)
g4 = g(4)
g2mg12 = np.where(abs(c) < 1e-7, (c*np.pi)**2.0/6.0, g2-g1**2.0)
- gam2k = np.where(abs(c) < 1e-7, np.pi**2.0/6.0,
- sc.expm1(sc.gammaln(2.0*c+1.0)-2*sc.gammaln(c + 1.0))/c**2.0)
+ def gam2k_f(c):
+ return sc.expm1(sc.gammaln(2.0*c+1.0)-2*sc.gammaln(c + 1.0))/c**2.0
+ gam2k = _lazywhere(abs(c) >= 1e-7, (c,), f=gam2k_f, fillvalue=np.pi**2.0/6.0)
eps = 1e-14
- gamk = np.where(abs(c) < eps, -_EULER, sc.expm1(sc.gammaln(c + 1))/c)
+ def gamk_f(c):
+ return sc.expm1(sc.gammaln(c + 1))/c
+ gamk = _lazywhere(abs(c) >= eps, (c,), f=gamk_f, fillvalue=-_EULER)
+ # mean
m = np.where(c < -1.0, np.nan, -gamk)
+
+ # variance
v = np.where(c < -0.5, np.nan, g1**2.0*gam2k)
# skewness
- sk1 = _lazywhere(c >= -1./3,
- (c, g1, g2, g3, g2mg12),
- lambda c, g1, g2, g3, g2mg12:
- np.sign(c)*(-g3 + (g2 + 2*g2mg12)*g1)/g2mg12**1.5,
- fillvalue=np.nan)
- sk = np.where(abs(c) <= eps**0.29, 12*np.sqrt(6)*_ZETA3/np.pi**3, sk1)
+ def sk1_eval(c, *args):
+ def sk1_eval_f(c, g1, g2, g3, g2mg12):
+ return np.sign(c)*(-g3 + (g2 + 2*g2mg12)*g1)/g2mg12**1.5
+ return _lazywhere(c >= -1./3, (c,)+args, f=sk1_eval_f, fillvalue=np.nan)
+
+ sk_fill = 12*np.sqrt(6)*_ZETA3/np.pi**3
+ args = (g1, g2, g3, g2mg12)
+ sk = _lazywhere(abs(c) > eps**0.29, (c,)+args, f=sk1_eval, fillvalue=sk_fill)
# kurtosis
- ku1 = _lazywhere(c >= -1./4,
- (g1, g2, g3, g4, g2mg12),
- lambda g1, g2, g3, g4, g2mg12:
- (g4 + (-4*g3 + 3*(g2 + g2mg12)*g1)*g1)/g2mg12**2,
- fillvalue=np.nan)
- ku = np.where(abs(c) <= (eps)**0.23, 12.0/5.0, ku1-3.0)
+ def ku1_eval(c, *args):
+ def ku1_eval_f(g1, g2, g3, g4, g2mg12):
+ return (g4 + (-4*g3 + 3*(g2 + g2mg12)*g1)*g1)/g2mg12**2 - 3
+ return _lazywhere(c >= -1./4, args, ku1_eval_f, fillvalue=np.nan)
+
+ args = (g1, g2, g3, g4, g2mg12)
+ ku = _lazywhere(abs(c) > eps**0.23, (c,)+args, f=ku1_eval, fillvalue=12.0/5.0)
+
return m, v, sk, ku
def _fitstart(self, data):
diff -Nru scipy-1.15.2/scipy/stats/tests/test_distributions.py scipy-1.15.3/scipy/stats/tests/test_distributions.py
--- scipy-1.15.2/scipy/stats/tests/test_distributions.py 2025-02-16 23:01:46.000000000 +0100
+++ scipy-1.15.3/scipy/stats/tests/test_distributions.py 2025-05-08 04:14:18.000000000 +0200
@@ -9251,6 +9251,28 @@
assert_equal(number_of_warnings_thrown, 0)
+def test_moments_gh22400():
+ # Regression test for gh-22400
+ # Check for correct results at c=0 with no warnings. While we're at it,
+ # check that NaN and sufficiently negative input produce NaNs, and output
+ # with `c=1` also agrees with reference values.
+ res = np.asarray(stats.genextreme.stats([0.0, np.nan, 1, -1.5], moments='mvsk'))
+
+ # Reference values for c=0 (Wikipedia)
+ mean = np.euler_gamma
+ var = np.pi**2 / 6
+ skew = 12 * np.sqrt(6) * special.zeta(3) / np.pi**3
+ kurt = 12 / 5
+ ref_0 = [mean, var, skew, kurt]
+ ref_1 = ref_3 = [np.nan]*4
+ ref_2 = [0, 1, -2, 6] # Wolfram Alpha, MaxStableDistribution[0, 1, -1]
+
+ assert_allclose(res[:, 0], ref_0, rtol=1e-14)
+ assert_equal(res[:, 1], ref_1)
+ assert_allclose(res[:, 2], ref_2, rtol=1e-14)
+ assert_equal(res[:, 3], ref_3)
+
+
def test_genextreme_entropy():
# regression test for gh-5181
euler_gamma = 0.5772156649015329
Reply to: