diff --git a/.github/workflows/linux_meson.yml b/.github/workflows/linux_meson.yml index 2ff6f8d36cae..141a1508e6ba 100644 --- a/.github/workflows/linux_meson.yml +++ b/.github/workflows/linux_meson.yml @@ -230,9 +230,9 @@ jobs: run: | pip install "numpy==1.21.6" && # build deps - pip install build meson-python ninja pythran "pybind11<2.11.0" "cython<3.0.0" "wheel<0.41.0" + pip install build "meson-python<0.14.0" ninja "pythran<0.14.0" "pybind11<2.11.0" "cython<3.0.0" "wheel<0.41.0" # test deps - pip install gmpy2 threadpoolctl mpmath pooch pythran "pybind11<2.11.0" pytest pytest-xdist==2.5.0 pytest-timeout + pip install gmpy2 threadpoolctl mpmath pooch "pythran<0.14.0" "pybind11<2.11.0" pytest pytest-xdist==2.5.0 pytest-timeout - name: Build wheel and install run: | diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index ec0f6fadef40..8af9a82e538d 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -131,14 +131,15 @@ jobs: CIBW_PRERELEASE_PYTHONS: True # required so that cp312 can grab a nightly wheel - # can remove when cp312 is release and we don't need to search + # can remove when cp312 is released and we don't need to search # other wheel locations. - CIBW_ENVIRONMENT: PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple PIP_PRE=1 + # PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple + CIBW_ENVIRONMENT: PIP_PRE=1 CIBW_ENVIRONMENT_WINDOWS: > PKG_CONFIG_PATH=c:/opt/64/lib/pkgconfig PIP_PRE=1 - PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple + # PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple # setting SDKROOT necessary when using the gfortran compiler # installed in cibw_before_build_macos.sh @@ -153,7 +154,7 @@ jobs: MACOS_DEPLOYMENT_TARGET=10.9 _PYTHON_HOST_PLATFORM=macosx-10.9-x86_64 PIP_PRE=1 - PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple + # PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple CIBW_REPAIR_WHEEL_COMMAND_MACOS: > DYLD_LIBRARY_PATH=/usr/local/lib delocate-listdeps {wheel} && diff --git a/.mailmap b/.mailmap index 56da1551d7bf..8fd76ba4da8a 100644 --- a/.mailmap +++ b/.mailmap @@ -167,6 +167,7 @@ Ed Schofield edschofield Egor Zemlyanoy egorz734 Egor Zemlyanoy Egorz734 Egor Zemlyanoy Egor +Ellie Litwack ellieLitwack Eric Larson Eric89GXL Eric Quintero e-q Eric Quintero Eric Quintero diff --git a/benchmarks/benchmarks/signal_filtering.py b/benchmarks/benchmarks/signal_filtering.py index b887bc90fd57..8c688b02caaa 100644 --- a/benchmarks/benchmarks/signal_filtering.py +++ b/benchmarks/benchmarks/signal_filtering.py @@ -6,7 +6,7 @@ with safe_import(): from scipy.signal import (lfilter, firwin, decimate, butter, sosfilt, - medfilt2d) + medfilt2d, freqz) class Decimate(Benchmark): @@ -102,3 +102,19 @@ def time_medfilt2d(self, threads): def peakmem_medfilt2d(self, threads): self._medfilt2d(threads) + + +class FreqzRfft(Benchmark): + param_names = ['whole', 'nyquist', 'worN'] + params = [ + [False, True], + [False, True], + [64, 65, 128, 129, 256, 257, 258, 512, 513, 65536, 65537, 65538], + ] + + def setup(self, whole, nyquist, worN): + self.y = np.zeros(worN) + self.y[worN//2] = 1.0 + + def time_freqz(self, whole, nyquist, worN): + freqz(self.y, whole=whole, include_nyquist=nyquist, worN=worN) diff --git a/ci/cirrus_wheels.yml b/ci/cirrus_wheels.yml index d3819fc94999..c53f2f983821 100644 --- a/ci/cirrus_wheels.yml +++ b/ci/cirrus_wheels.yml @@ -30,8 +30,8 @@ cirrus_wheels_linux_aarch64_task: env: CIBW_PRERELEASE_PYTHONS: True CIBW_ENVIRONMENT: > - PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple PIP_PRE=1 + # PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple build_script: | apt install -y python3-venv python-is-python3 @@ -59,8 +59,8 @@ cirrus_wheels_macos_arm64_task: CIBW_ENVIRONMENT: > MACOSX_DEPLOYMENT_TARGET=12.0 _PYTHON_HOST_PLATFORM="macosx-12.0-arm64" - PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple PIP_PRE=1 + # PIP_EXTRA_INDEX_URL=https://pypi.anaconda.org/scientific-python-nightly-wheels/simple PKG_CONFIG_PATH: /opt/arm64-builds/lib/pkgconfig # assumes that the cmake config is in /usr/local/lib/cmake CMAKE_PREFIX_PATH: /opt/arm64-builds/ diff --git a/dev.py b/dev.py index 099b9dc36d2d..e39c7b521c3b 100644 --- a/dev.py +++ b/dev.py @@ -111,14 +111,6 @@ def task_meta(cls, **kwargs): import traceback from concurrent.futures.process import _MAX_WINDOWS_WORKERS -# distutils is required to infer meson install path -# if this needs to be replaced for Python 3.12 support and there's no -# stdlib alternative, use CmdAction and the hack discussed in gh-16058 -with warnings.catch_warnings(): - warnings.filterwarnings("ignore", category=DeprecationWarning) - from distutils import dist - from distutils.command.install import INSTALL_SCHEMES - from pathlib import Path from collections import namedtuple from types import ModuleType as new_module @@ -337,14 +329,23 @@ def get_site_packages(self): Depending on whether we have debian python or not, return dist_packages path or site_packages path. """ - if 'deb_system' in INSTALL_SCHEMES: - # debian patched python in use - install_cmd = dist.Distribution().get_command_obj('install') - install_cmd.select_scheme('deb_system') - install_cmd.finalize_options() - plat_path = Path(install_cmd.install_platlib) - else: + if sys.version_info >= (3, 12): plat_path = Path(get_path('platlib')) + else: + # distutils is required to infer meson install path + # for python < 3.12 in debian patched python + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=DeprecationWarning) + from distutils import dist + from distutils.command.install import INSTALL_SCHEMES + if 'deb_system' in INSTALL_SCHEMES: + # debian patched python in use + install_cmd = dist.Distribution().get_command_obj('install') + install_cmd.select_scheme('deb_system') + install_cmd.finalize_options() + plat_path = Path(install_cmd.install_platlib) + else: + plat_path = Path(get_path('platlib')) return self.installed / plat_path.relative_to(sys.exec_prefix) diff --git a/doc/source/release/1.11.3-notes.rst b/doc/source/release/1.11.3-notes.rst index 65686bd93aa4..507ba41123d8 100644 --- a/doc/source/release/1.11.3-notes.rst +++ b/doc/source/release/1.11.3-notes.rst @@ -12,13 +12,65 @@ compared to 1.11.2. Authors ======= * Name (commits) +* Jake Bowhay (2) +* CJ Carey (1) +* Colin Carroll (1) + +* Anirudh Dagar (2) +* drestebon (1) + +* Ralf Gommers (5) +* Matt Haberland (2) +* Julien Jerphanion (1) +* Uwe L. Korn (1) + +* Ellie Litwack (2) +* Andrew Nelson (5) +* Bharat Raghunathan (1) +* Tyler Reddy (37) +* Søren Fuglede Jørgensen (2) +* Hielke Walinga (1) + +* Warren Weckesser (1) +* Bernhard M. Wiedemann (1) + +A total of 17 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.11.3 ------------------------ +* `#15093 `__: BUG: scipy.optimize's trust-constr algorithm hangs when keep-feasible... +* `#15273 `__: freqz: suboptimal performance for worN=2\*\*n+1, include_nyquist=True... +* `#17269 `__: Bug in scipy.sparse.csgraph.min_weight_full_bipartite_matching +* `#17289 `__: BUG: Different results between numpy.fft.rfft and scipy.signal.freqz +* `#18716 `__: Buffer dtype mismatch, expected 'ITYPE_t' but got 'long' +* `#18782 `__: BUG: johnsonsu distribution no longer accepts integer \`b\` parameter +* `#18922 `__: BUG: dev.py has \`distutils\` usage +* `#19101 `__: BUG: mesonpy embeds random path in .pyx files +* `#19103 `__: BUG: Regression in 1.11.2: optimize.least_squares with method='trf'... +* `#19132 `__: BUG: Build fails on latest commit +* `#19149 `__: BUG: scipy.sparse.csgraph.laplacian raises AttributeError on... +* `#19197 `__: BUG: Incorrect sampling from zero rank covariance Pull requests for 1.11.3 ------------------------ +* `#17633 `__: BUG: add infeasibility checks to min_weight_full_bipartite_matching +* `#18784 `__: BUG: Allow johnsonsu parameters to be floats +* `#18913 `__: BUG: sparse.csgraph: Support int64 indices in traversal.pyx +* `#18924 `__: BUG: Fix python3.12 distutils dev.py build +* `#18956 `__: BUG: trust-constr Bounds exclusive +* `#19076 `__: MAINT: should not be using np.float64() on arrays +* `#19084 `__: REL, MAINT: prep for 1.11.3 +* `#19111 `__: BUG: Fixes #19103 by adding back make_strictly_feasible to lsq... +* `#19123 `__: BLD: Avoid absolute pathnames in .pyx files +* `#19135 `__: MAINT: signal: Remove the cval parameter from the private function... +* `#19139 `__: BLD: revert to using published wheels [wheel build] +* `#19156 `__: BUG: Support sparse arrays in scipy.sparse.csgraph.laplacian +* `#19199 `__: MAINT: stats.CovViaEigendecomposition: fix \`_LA\` attribute... +* `#19200 `__: TST: fix \`TestODR.test_implicit\` test failure with tolerance... +* `#19208 `__: BUG: freqz rfft grid fix +* `#19280 `__: MAINT: newton, make sure x0 is an inexact type +* `#19286 `__: BUG: stats: fix build failure due to incorrect Boost policies... +* `#19290 `__: BLD: add float.h include to \`_fpumode.c\`, fixes Clang on Windows... +* `#19299 `__: MAINT: fix libquadmath licence diff --git a/pyproject.toml b/pyproject.toml index bc0519e40439..95d0ccbe91c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,13 +10,15 @@ [build-system] build-backend = 'mesonpy' requires = [ - "meson-python>=0.12.1,<0.14.0", # already working with 0.13.x series at branch time - "Cython>=0.29.35,<3.0", # when updating version, also update check in meson.build - "pybind11>=2.10.4,<2.11.0", - "pythran>=0.12.0,<0.14.0", # already working with 0.13.x series at branch time - # `wheel` is needed for non-isolated builds, given that `meson-python` - # doesn't list it as a runtime requirement (at least in 0.5.0) - "wheel<0.41.0", + # Reason for `<`: future-proofing (0.14.0 released and known to work) + "meson-python>=0.12.1,<0.15.0", + # Reason for `<`: we want to build the 1.11.x releases with Cython 0.29.x + # (too many changes in 3.0) + "Cython>=0.29.35,<3.0", # when updating version, also update check in meson.build + # Reason for `<`: future-proofing (2.11.1 is the most recent version) + "pybind11>=2.10.4,<2.11.1", + # Reason for `<`: future-proofing (0.14.0 released and known to work) + "pythran>=0.12.0,<0.15.0", # NumPy dependencies - to update these, sync from # https://github.com/scipy/oldest-supported-numpy/, and then @@ -41,13 +43,14 @@ requires = [ # (see oldest-supported-numpy issues gh-28 and gh-45) "numpy==1.21.6; python_version=='3.10' and (platform_system!='Windows' and platform_machine!='loongarch64') and platform_python_implementation != 'PyPy'", "numpy==1.23.2; python_version=='3.11' and platform_python_implementation != 'PyPy'", + "numpy>=1.26.0,<1.27.0; python_version=='3.12'", # no `==` because we expect relevant bug fixes in 1.26.1/2 # For Python versions which aren't yet officially supported, # we specify an unpinned NumPy which allows source distributions # to be used and allows wheels to be used as soon as they # become available. "numpy; python_version>='3.9' and platform_python_implementation=='PyPy'", - "numpy==1.26.0b1; python_version>='3.12.0.rc1'", + "numpy>=1.26.0; python_version>='3.13'", ] [project] @@ -64,8 +67,6 @@ maintainers = [ # https://scipy.github.io/devdocs/dev/core-dev/index.html#version-ranges-for-numpy-and-other-dependencies requires-python = ">=3.9,<3.13" dependencies = [ - # TODO: update to "pin-compatible" once possible, see - # https://github.com/FFY00/meson-python/issues/29 "numpy>=1.21.6,<1.28.0", ] readme = "README.rst" @@ -80,6 +81,7 @@ classifiers = [ "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", "Topic :: Software Development :: Libraries", "Topic :: Scientific/Engineering", "Operating System :: Microsoft :: Windows", diff --git a/scipy/_lib/_fpumode.c b/scipy/_lib/_fpumode.c index fb5cd7dc5b02..cb0e4440fcb3 100644 --- a/scipy/_lib/_fpumode.c +++ b/scipy/_lib/_fpumode.c @@ -4,6 +4,7 @@ #ifdef _MSC_VER +#include #pragma float_control(precise, on) #pragma fenv_access (on) #endif diff --git a/scipy/integrate/_quadrature.py b/scipy/integrate/_quadrature.py index 1048e352675e..1b9e2841b6b1 100644 --- a/scipy/integrate/_quadrature.py +++ b/scipy/integrate/_quadrature.py @@ -502,8 +502,8 @@ def _basic_simpson(y, start, stop, x, dx, axis): h = np.diff(x, axis=axis) sl0 = tupleset(slice_all, axis, slice(start, stop, step)) sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step)) - h0 = np.float64(h[sl0]) - h1 = np.float64(h[sl1]) + h0 = h[sl0].astype(float, copy=False) + h1 = h[sl1].astype(float, copy=False) hsum = h0 + h1 hprod = h0 * h1 h0divh1 = np.true_divide(h0, h1, out=np.zeros_like(h0), where=h1 != 0) diff --git a/scipy/odr/tests/test_odr.py b/scipy/odr/tests/test_odr.py index 6391eb3785d6..c05f764d4b02 100644 --- a/scipy/odr/tests/test_odr.py +++ b/scipy/odr/tests/test_odr.py @@ -5,7 +5,8 @@ import numpy as np from numpy import pi from numpy.testing import (assert_array_almost_equal, - assert_equal, assert_warns) + assert_equal, assert_warns, + assert_allclose) import pytest from pytest import raises as assert_raises @@ -127,7 +128,7 @@ def test_implicit(self): np.array([0.1113840353364371, 0.1097673310686467, 0.0041060738314314, 0.0027500347539902, 0.0034962501532468]), ) - assert_array_almost_equal( + assert_allclose( out.cov_beta, np.array([[2.1089274602333052e+00, -1.9437686411979040e+00, 7.0263550868344446e-02, -4.7175267373474862e-02, @@ -144,6 +145,7 @@ def test_implicit(self): [5.2515575927380355e-02, -5.8822307501391467e-02, 1.4528860663055824e-03, -1.2692942951415293e-03, 2.0778813389755596e-03]]), + rtol=1e-6, atol=2e-6, ) # Multi-variable Example diff --git a/scipy/optimize/_lsq/least_squares.py b/scipy/optimize/_lsq/least_squares.py index 40d03db77b08..2cb2475d5a34 100644 --- a/scipy/optimize/_lsq/least_squares.py +++ b/scipy/optimize/_lsq/least_squares.py @@ -12,7 +12,7 @@ from .trf import trf from .dogbox import dogbox -from .common import EPS, in_bounds +from .common import EPS, in_bounds, make_strictly_feasible TERMINATION_MESSAGES = { @@ -821,6 +821,9 @@ def least_squares( ftol, xtol, gtol = check_tolerance(ftol, xtol, gtol, method) + if method == 'trf': + x0 = make_strictly_feasible(x0, lb, ub, rstep=0) + def fun_wrapped(x): return np.atleast_1d(fun(x, *args, **kwargs)) diff --git a/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py b/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py index 1c86dd869b78..2835ea5445c0 100644 --- a/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py +++ b/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py @@ -3,7 +3,7 @@ from scipy.sparse.linalg import LinearOperator from .._differentiable_functions import VectorFunction from .._constraints import ( - NonlinearConstraint, LinearConstraint, PreparedConstraint, strict_bounds) + NonlinearConstraint, LinearConstraint, PreparedConstraint, Bounds, strict_bounds) from .._hessian_update_strategy import BFGS from .._optimize import OptimizeResult from .._differentiable_functions import ScalarFunction @@ -323,6 +323,11 @@ def _minimize_trustregion_constr(fun, x0, args, grad, verbose = 1 if bounds is not None: + modified_lb = np.nextafter(bounds.lb, -np.inf, where=bounds.lb > -np.inf) + modified_ub = np.nextafter(bounds.ub, np.inf, where=bounds.ub < np.inf) + modified_lb = np.where(np.isfinite(bounds.lb), modified_lb, bounds.lb) + modified_ub = np.where(np.isfinite(bounds.ub), modified_ub, bounds.ub) + bounds = Bounds(modified_lb, modified_ub, keep_feasible=bounds.keep_feasible) finite_diff_bounds = strict_bounds(bounds.lb, bounds.ub, bounds.keep_feasible, n_vars) else: diff --git a/scipy/optimize/_zeros_py.py b/scipy/optimize/_zeros_py.py index 61306f9ce889..af330880e798 100644 --- a/scipy/optimize/_zeros_py.py +++ b/scipy/optimize/_zeros_py.py @@ -290,7 +290,10 @@ class of similar problems can be solved together. full_output) # Convert to float (don't use float(x0); this works also for complex x0) - x0 = np.asarray(x0)[()] + # Use np.asarray because we want x0 to be a numpy object, not a Python + # object. e.g. np.complex(1+1j) > 0 is possible, but (1 + 1j) > 0 raises + # a TypeError + x0 = np.asarray(x0)[()] * 1.0 p0 = x0 funcalls = 0 if fprime is not None: diff --git a/scipy/optimize/tests/test_least_squares.py b/scipy/optimize/tests/test_least_squares.py index 667a2b765b87..438a2e12e4a3 100644 --- a/scipy/optimize/tests/test_least_squares.py +++ b/scipy/optimize/tests/test_least_squares.py @@ -822,3 +822,31 @@ def chi2(x): # if we choose an initial condition that is close to the solution # we shouldn't return an answer that is further away from the solution assert_allclose(res.x, answer, atol=initial_guess-answer) + + +def test_gh_19103(): + # Checks that least_squares trf method selects a strictly feasible point, + # and thus succeeds instead of failing, + # when the initial guess is reported exactly at a boundary point. + # This is a reduced example from gh191303 + + ydata = np.array([0.] * 66 + [ + 1., 0., 0., 0., 0., 0., 1., 1., 0., 0., 1., + 1., 1., 1., 0., 0., 0., 1., 0., 0., 2., 1., + 0., 3., 1., 6., 5., 0., 0., 2., 8., 4., 4., + 6., 9., 7., 2., 7., 8., 2., 13., 9., 8., 11., + 10., 13., 14., 19., 11., 15., 18., 26., 19., 32., 29., + 28., 36., 32., 35., 36., 43., 52., 32., 58., 56., 52., + 67., 53., 72., 88., 77., 95., 94., 84., 86., 101., 107., + 108., 118., 96., 115., 138., 137., + ]) + xdata = np.arange(0, ydata.size) * 0.1 + + def exponential_wrapped(params): + A, B, x0 = params + return A * np.exp(B * (xdata - x0)) - ydata + + x0 = [0.01, 1., 5.] + bounds = ((0.01, 0, 0), (np.inf, 10, 20.9)) + res = least_squares(exponential_wrapped, x0, method='trf', bounds=bounds) + assert res.success diff --git a/scipy/optimize/tests/test_minimize_constrained.py b/scipy/optimize/tests/test_minimize_constrained.py index 0c4a8d4676ad..7bb54a9e0ecb 100644 --- a/scipy/optimize/tests/test_minimize_constrained.py +++ b/scipy/optimize/tests/test_minimize_constrained.py @@ -615,6 +615,32 @@ def callback(x, info): # compatibility assert_(result.get('niter', -1) == 1) + def test_issue_15093(self): + # scipy docs define bounds as inclusive, so it shouldn't be + # an issue to set x0 on the bounds even if keep_feasible is + # True. Previously, trust-constr would treat bounds as + # exclusive. + + x0 = np.array([0., 0.5]) + + def obj(x): + x1 = x[0] + x2 = x[1] + return x1 ** 2 + x2 ** 2 + + bounds = Bounds(np.array([0., 0.]), np.array([1., 1.]), + keep_feasible=True) + + with suppress_warnings() as sup: + sup.filter(UserWarning, "delta_grad == 0.0") + result = minimize( + method='trust-constr', + fun=obj, + x0=x0, + bounds=bounds) + + assert result['success'] + class TestEmptyConstraint(TestCase): """ Here we minimize x^2+y^2 subject to x^2-y^2>1. diff --git a/scipy/optimize/tests/test_zeros.py b/scipy/optimize/tests/test_zeros.py index 8ce10049f235..e817d37a9d32 100644 --- a/scipy/optimize/tests/test_zeros.py +++ b/scipy/optimize/tests/test_zeros.py @@ -456,6 +456,20 @@ def test_gh17570_defaults(self): != res_newton_default.iterations == res_newton_default.function_calls/2) # newton 2-point diff + @pytest.mark.parametrize('method', ['secant', 'newton']) + def test_int_x0_gh19280(self, method): + # Originally, `newton` ensured that only floats were passed to the + # callable. This was indadvertently changed by gh-17669. Check that + # it has been changed back. + def f(x): + # an integer raised to a negative integer power would fail + return x**-2 - 2 + + res = root_scalar(f, x0=1, method=method) + assert res.converged + assert_allclose(abs(res.root), 2**-0.5) + assert res.root.dtype == np.dtype(np.float64) + def test_gh_5555(): root = 0.1 diff --git a/scipy/signal/_filter_design.py b/scipy/signal/_filter_design.py index 845b8dc9c67a..4f1fc2e5cf6f 100644 --- a/scipy/signal/_filter_design.py +++ b/scipy/signal/_filter_design.py @@ -447,12 +447,12 @@ def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi, lastpoint = 2 * pi if whole else pi # if include_nyquist is true and whole is false, w should # include end point - w = np.linspace(0, lastpoint, N, endpoint=include_nyquist and not whole) - if (a.size == 1 and N >= b.shape[0] and - sp_fft.next_fast_len(N) == N and - (b.ndim == 1 or (b.shape[-1] == 1))): - # if N is fast, 2 * N will be fast, too, so no need to check - n_fft = N if whole else N * 2 + w = np.linspace(0, lastpoint, N, + endpoint=include_nyquist and not whole) + n_fft = N if whole else 2 * (N - 1) if include_nyquist else 2 * N + if (a.size == 1 and (b.ndim == 1 or (b.shape[-1] == 1)) + and n_fft >= b.shape[0] + and n_fft > 0): # TODO: review threshold acc. to benchmark? if np.isrealobj(b) and np.isrealobj(a): fft_func = sp_fft.rfft else: diff --git a/scipy/signal/_upfirdn_apply.pyx b/scipy/signal/_upfirdn_apply.pyx index 4ba5e703b8f4..6419816b1984 100644 --- a/scipy/signal/_upfirdn_apply.pyx +++ b/scipy/signal/_upfirdn_apply.pyx @@ -234,7 +234,7 @@ cdef DTYPE_t _extend_right(DTYPE_t *x, np.intp_t idx, np.intp_t len_x, @cython.boundscheck(False) @cython.wraparound(False) cpdef _pad_test(np.ndarray[DTYPE_t] data, np.intp_t npre=0, np.intp_t npost=0, - object mode=0, DTYPE_t cval=0): + object mode=0): """1D test function for signal extension modes. Returns ``data extended by ``npre``, ``npost`` at the beginning, end. @@ -264,9 +264,9 @@ cpdef _pad_test(np.ndarray[DTYPE_t] data, np.intp_t npre=0, np.intp_t npost=0, with nogil: for idx in range(-npre, len_x + npost, 1): if idx < 0: - xval = _extend_left(data_ptr, idx, len_x, _mode, cval) + xval = _extend_left(data_ptr, idx, len_x, _mode, 0.0) elif idx >= len_x: - xval = _extend_right(data_ptr, idx, len_x, _mode, cval) + xval = _extend_right(data_ptr, idx, len_x, _mode, 0.0) else: xval = data_ptr[idx] out[cnt] = xval diff --git a/scipy/signal/tests/test_filter_design.py b/scipy/signal/tests/test_filter_design.py index 1a8d6ca553e1..28bb2bb5f561 100644 --- a/scipy/signal/tests/test_filter_design.py +++ b/scipy/signal/tests/test_filter_design.py @@ -893,6 +893,23 @@ def test_nyquist(self): assert_array_almost_equal(w1, w2) assert_array_almost_equal(h1, h2) + # https://github.com/scipy/scipy/issues/17289 + # https://github.com/scipy/scipy/issues/15273 + @pytest.mark.parametrize('whole,nyquist,worN', + [(False, False, 32), + (False, True, 32), + (True, False, 32), + (True, True, 32), + (False, False, 257), + (False, True, 257), + (True, False, 257), + (True, True, 257)]) + def test_17289(self, whole, nyquist, worN): + d = [0, 1] + w, Drfft = freqz(d, worN=32, whole=whole, include_nyquist=nyquist) + _, Dpoly = freqz(d, worN=w) + assert_allclose(Drfft, Dpoly) + class TestSOSFreqz: diff --git a/scipy/sparse/csgraph/_laplacian.py b/scipy/sparse/csgraph/_laplacian.py index 3ac929239a5c..36453f3243bd 100644 --- a/scipy/sparse/csgraph/_laplacian.py +++ b/scipy/sparse/csgraph/_laplacian.py @@ -403,11 +403,11 @@ def _laplacian_sparse_flo(graph, normed, axis, copy, form, dtype, symmetrized): if dtype is None: dtype = graph.dtype - graph_sum = graph.sum(axis=axis).getA1() + graph_sum = np.asarray(graph.sum(axis=axis)).ravel() graph_diagonal = graph.diagonal() diag = graph_sum - graph_diagonal if symmetrized: - graph_sum += graph.sum(axis=1 - axis).getA1() + graph_sum += np.asarray(graph.sum(axis=1 - axis)).ravel() diag = graph_sum - graph_diagonal - graph_diagonal if normed: @@ -456,7 +456,7 @@ def _laplacian_sparse(graph, normed, axis, copy, form, dtype, symmetrized): if symmetrized: m += m.T.conj() - w = m.sum(axis=axis).getA1() - m.diagonal() + w = np.asarray(m.sum(axis=axis)).ravel() - m.diagonal() if normed: m = m.tocoo(copy=needs_copy) isolated_node_mask = (w == 0) diff --git a/scipy/sparse/csgraph/_matching.pyx b/scipy/sparse/csgraph/_matching.pyx index 1a4f582a74e1..07477694b9de 100644 --- a/scipy/sparse/csgraph/_matching.pyx +++ b/scipy/sparse/csgraph/_matching.pyx @@ -471,11 +471,22 @@ def min_weight_full_bipartite_matching(biadjacency_matrix, maximize=False): a = np.arange(np.min(biadjacency_matrix.shape)) - # The algorithm expects more columns than rows in the graph. + # The algorithm expects more columns than rows in the graph, so + # we use the transpose if that is not already the case. We also + # ensure that we have a full matching. In principle, it should be + # possible to avoid this check for a performance improvement, by + # checking for infeasibility during the execution of _lapjvsp below + # instead, but some cases are not yet handled there. if j < i: biadjacency_matrix_t = biadjacency_matrix.T if biadjacency_matrix_t.format != "csr": biadjacency_matrix_t = biadjacency_matrix_t.tocsr() + matching, _ = _hopcroft_karp(biadjacency_matrix_t.indices, + biadjacency_matrix_t.indptr, + j, i) + matching = np.asarray(matching) + if np.sum(matching != -1) != min(i, j): + raise ValueError('no full matching exists') b = np.asarray(_lapjvsp(biadjacency_matrix_t.indptr, biadjacency_matrix_t.indices, biadjacency_matrix_t.data, @@ -485,6 +496,12 @@ def min_weight_full_bipartite_matching(biadjacency_matrix, maximize=False): else: if biadjacency_matrix.format != "csr": biadjacency_matrix = biadjacency_matrix.tocsr() + matching, _ = _hopcroft_karp(biadjacency_matrix.indices, + biadjacency_matrix.indptr, + i, j) + matching = np.asarray(matching) + if np.sum(matching != -1) != min(i, j): + raise ValueError('no full matching exists') b = np.asarray(_lapjvsp(biadjacency_matrix.indptr, biadjacency_matrix.indices, biadjacency_matrix.data, diff --git a/scipy/sparse/csgraph/_traversal.pyx b/scipy/sparse/csgraph/_traversal.pyx index ec5bd294cdc7..2a518da16a3b 100644 --- a/scipy/sparse/csgraph/_traversal.pyx +++ b/scipy/sparse/csgraph/_traversal.pyx @@ -184,7 +184,7 @@ def breadth_first_tree(csgraph, i_start, directed=True): Note that the resulting graph is a Directed Acyclic Graph which spans the graph. A breadth-first tree from a given node is unique. """ - node_list, predecessors = breadth_first_order(csgraph, i_start, + _node_list, predecessors = breadth_first_order(csgraph, i_start, directed, True) return reconstruct_path(csgraph, predecessors, directed) @@ -258,7 +258,7 @@ def depth_first_tree(csgraph, i_start, directed=True): had begun with the edge connecting nodes 0 and 3, the result would have been different. """ - node_list, predecessors = depth_first_order(csgraph, i_start, + _node_list, predecessors = depth_first_order(csgraph, i_start, directed, True) return reconstruct_path(csgraph, predecessors, directed) @@ -352,12 +352,22 @@ cpdef breadth_first_order(csgraph, i_start, return node_list[:length] -cdef unsigned int _breadth_first_directed( - unsigned int head_node, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr, - np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, - np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors) noexcept: +def _breadth_first_directed( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors): + return _breadth_first_directed2(head_node, indices, indptr, + node_list, predecessors) + + +cdef unsigned int _breadth_first_directed2( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors) noexcept: # Inputs: # head_node: (input) index of the node from which traversal starts # indices: (input) CSR indices of graph @@ -390,15 +400,25 @@ cdef unsigned int _breadth_first_directed( return i_nl - -cdef unsigned int _breadth_first_undirected( - unsigned int head_node, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices1, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr1, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices2, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr2, - np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, - np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors) noexcept: +def _breadth_first_undirected( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices1, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr1, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indices2, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indptr2, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors): + return _breadth_first_undirected2(head_node, indices1, indptr1, indices2, + indptr2, node_list, predecessors) + +cdef unsigned int _breadth_first_undirected2( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices1, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr1, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indices2, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indptr2, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors) noexcept: # Inputs: # head_node: (input) index of the node from which traversal starts # indices1: (input) CSR indices of graph @@ -538,14 +558,27 @@ cpdef depth_first_order(csgraph, i_start, return node_list[:length] -cdef unsigned int _depth_first_directed( - unsigned int head_node, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr, - np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, - np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors, - np.ndarray[ITYPE_t, ndim=1, mode='c'] root_list, - np.ndarray[ITYPE_t, ndim=1, mode='c'] flag) noexcept: +def _depth_first_directed( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors, + np.ndarray[ITYPE_t, ndim=1, mode='c'] root_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] flag): + return _depth_first_directed2(head_node, indices, indptr, + node_list, predecessors, + root_list, flag) + + +cdef unsigned int _depth_first_directed2( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors, + np.ndarray[ITYPE_t, ndim=1, mode='c'] root_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] flag) noexcept: cdef unsigned int i, i_nl_end, cnode, pnode cdef unsigned int N = node_list.shape[0] cdef int no_children, i_root @@ -582,16 +615,32 @@ cdef unsigned int _depth_first_directed( return i_nl_end -cdef unsigned int _depth_first_undirected( - unsigned int head_node, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices1, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr1, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices2, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr2, - np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, - np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors, - np.ndarray[ITYPE_t, ndim=1, mode='c'] root_list, - np.ndarray[ITYPE_t, ndim=1, mode='c'] flag) noexcept: +def _depth_first_undirected( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices1, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr1, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indices2, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indptr2, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors, + np.ndarray[ITYPE_t, ndim=1, mode='c'] root_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] flag): + return _depth_first_undirected2(head_node, indices1, indptr1, + indices2, indptr2, + node_list, predecessors, + root_list, flag) + + +cdef unsigned int _depth_first_undirected2( + unsigned int head_node, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices1, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr1, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indices2, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indptr2, + np.ndarray[ITYPE_t, ndim=1, mode='c'] node_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] predecessors, + np.ndarray[ITYPE_t, ndim=1, mode='c'] root_list, + np.ndarray[ITYPE_t, ndim=1, mode='c'] flag) noexcept: cdef unsigned int i, i_nl_end, cnode, pnode cdef unsigned int N = node_list.shape[0] cdef int no_children, i_root @@ -644,10 +693,17 @@ cdef unsigned int _depth_first_undirected( return i_nl_end -cdef int _connected_components_directed( - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr, - np.ndarray[ITYPE_t, ndim=1, mode='c'] labels) noexcept: +def _connected_components_directed( + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr, + np.ndarray[ITYPE_t, ndim=1, mode='c'] labels): + return _connected_components_directed2(indices, indptr, labels) + + +cdef int _connected_components_directed2( + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr, + np.ndarray[ITYPE_t, ndim=1, mode='c'] labels) noexcept: """ Uses an iterative version of Tarjan's algorithm to find the strongly connected components of a directed graph represented as a @@ -766,12 +822,24 @@ cdef int _connected_components_directed( labels += (N - 1) return (N - 1) - label -cdef int _connected_components_undirected( - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices1, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr1, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indices2, - np.ndarray[ITYPE_t, ndim=1, mode='c'] indptr2, - np.ndarray[ITYPE_t, ndim=1, mode='c'] labels) noexcept: + +def _connected_components_undirected( + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices1, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr1, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indices2, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indptr2, + np.ndarray[ITYPE_t, ndim=1, mode='c'] labels): + return _connected_components_undirected2(indices1, indptr1, + indices2, indptr2, + labels) + + +cdef int _connected_components_undirected2( + np.ndarray[int32_or_int64, ndim=1, mode='c'] indices1, + np.ndarray[int32_or_int64, ndim=1, mode='c'] indptr1, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indices2, + np.ndarray[int32_or_int64_b, ndim=1, mode='c'] indptr2, + np.ndarray[ITYPE_t, ndim=1, mode='c'] labels) noexcept: cdef int v, w, j, label, SS_head cdef int N = labels.shape[0] diff --git a/scipy/sparse/csgraph/parameters.pxi b/scipy/sparse/csgraph/parameters.pxi index 91cbdd640441..5391bedf4d4f 100644 --- a/scipy/sparse/csgraph/parameters.pxi +++ b/scipy/sparse/csgraph/parameters.pxi @@ -10,6 +10,11 @@ ctypedef fused int32_or_int64: np.int32_t np.int64_t +# Another copy of the same fused type, for working with mixed-type functions. +ctypedef fused int32_or_int64_b: + np.int32_t + np.int64_t + # EPS is the precision of DTYPE DEF DTYPE_EPS = 1E-15 diff --git a/scipy/sparse/csgraph/tests/test_connected_components.py b/scipy/sparse/csgraph/tests/test_connected_components.py index 681f29fd48a3..0b190a24deb9 100644 --- a/scipy/sparse/csgraph/tests/test_connected_components.py +++ b/scipy/sparse/csgraph/tests/test_connected_components.py @@ -1,6 +1,6 @@ import numpy as np from numpy.testing import assert_equal, assert_array_almost_equal -from scipy.sparse import csgraph +from scipy.sparse import csgraph, csr_array def test_weak_connections(): @@ -97,3 +97,23 @@ def test_fully_connected_graph(): g = np.ones((4, 4)) n_components, labels = csgraph.connected_components(g) assert_equal(n_components, 1) + + +def test_int64_indices_undirected(): + # See https://github.com/scipy/scipy/issues/18716 + g = csr_array(([1], np.array([[0], [1]], dtype=np.int64)), shape=(2, 2)) + assert g.indices.dtype == np.int64 + n, labels = csgraph.connected_components(g, directed=False) + assert n == 1 + assert_array_almost_equal(labels, [0, 0]) + + +def test_int64_indices_directed(): + # See https://github.com/scipy/scipy/issues/18716 + g = csr_array(([1], np.array([[0], [1]], dtype=np.int64)), shape=(2, 2)) + assert g.indices.dtype == np.int64 + n, labels = csgraph.connected_components(g, directed=True, + connection='strong') + assert n == 2 + assert_array_almost_equal(labels, [1, 0]) + diff --git a/scipy/sparse/csgraph/tests/test_graph_laplacian.py b/scipy/sparse/csgraph/tests/test_graph_laplacian.py index e5648ba77718..e2ec68cd1d8b 100644 --- a/scipy/sparse/csgraph/tests/test_graph_laplacian.py +++ b/scipy/sparse/csgraph/tests/test_graph_laplacian.py @@ -170,7 +170,9 @@ def _check_laplacian_dtype( @pytest.mark.parametrize("dtype", DTYPES) @pytest.mark.parametrize("arr_type", [np.array, sparse.csr_matrix, - sparse.coo_matrix]) + sparse.coo_matrix, + sparse.csr_array, + sparse.coo_array]) @pytest.mark.parametrize("copy", [True, False]) @pytest.mark.parametrize("normed", [True, False]) @pytest.mark.parametrize("use_out_degree", [True, False]) @@ -244,7 +246,11 @@ def test_sparse_formats(fmt, normed, copy): @pytest.mark.parametrize( - "arr_type", [np.asarray, sparse.csr_matrix, sparse.coo_matrix] + "arr_type", [np.asarray, + sparse.csr_matrix, + sparse.coo_matrix, + sparse.csr_array, + sparse.coo_array] ) @pytest.mark.parametrize("form", ["array", "function", "lo"]) def test_laplacian_symmetrized(arr_type, form): @@ -301,7 +307,11 @@ def test_laplacian_symmetrized(arr_type, form): @pytest.mark.parametrize( - "arr_type", [np.asarray, sparse.csr_matrix, sparse.coo_matrix] + "arr_type", [np.asarray, + sparse.csr_matrix, + sparse.coo_matrix, + sparse.csr_array, + sparse.coo_array] ) @pytest.mark.parametrize("dtype", DTYPES) @pytest.mark.parametrize("normed", [True, False]) diff --git a/scipy/sparse/csgraph/tests/test_matching.py b/scipy/sparse/csgraph/tests/test_matching.py index 387aa5e07159..87e2920fe971 100644 --- a/scipy/sparse/csgraph/tests/test_matching.py +++ b/scipy/sparse/csgraph/tests/test_matching.py @@ -161,6 +161,7 @@ def test_min_weight_full_matching_trivial_graph(num_rows, num_cols): [ [[1, 1, 1], [1, 0, 0], [1, 0, 0]], [[1, 1, 1], [0, 0, 1], [0, 0, 1]], + [[1, 0, 0, 1], [1, 1, 0, 1], [0, 0, 0, 0]], [[1, 0, 0], [2, 0, 0]], [[0, 1, 0], [0, 2, 0]], [[1, 0], [2, 0], [5, 0]] @@ -170,6 +171,60 @@ def test_min_weight_full_matching_infeasible_problems(biadjacency_matrix): min_weight_full_bipartite_matching(csr_matrix(biadjacency_matrix)) +def test_min_weight_full_matching_large_infeasible(): + # Regression test for GitHub issue #17269 + a = np.asarray([ + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001], + [0.0, 0.11687445, 0.0, 0.0, 0.01319788, 0.07509257, 0.0, + 0.0, 0.0, 0.74228317, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.81087935, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.8408466, 0.0, 0.0, 0.0, 0.0, 0.01194389, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.82994211, 0.0, 0.0, 0.0, 0.11468516, 0.0, 0.0, 0.0, + 0.11173505, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0], + [0.18796507, 0.0, 0.04002318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.75883335, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.71545464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.02748488, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.78470564, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.14829198, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.10870609, 0.0, 0.0, 0.0, 0.8918677, 0.0, 0.0, 0.0, 0.06306644, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.63844085, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7442354, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.09850549, 0.0, 0.0, 0.18638258, + 0.2769244, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.73182464, 0.0, 0.0, 0.46443561, + 0.38589284, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.29510278, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.09666032, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + ]) + with pytest.raises(ValueError, match='no full matching exists'): + min_weight_full_bipartite_matching(csr_matrix(a)) + + def test_explicit_zero_causes_warning(): with pytest.warns(UserWarning): biadjacency_matrix = csr_matrix(((2, 0, 3), (0, 1, 1), (0, 2, 3))) diff --git a/scipy/sparse/csgraph/tests/test_traversal.py b/scipy/sparse/csgraph/tests/test_traversal.py index 026fbe275484..414e2d14864d 100644 --- a/scipy/sparse/csgraph/tests/test_traversal.py +++ b/scipy/sparse/csgraph/tests/test_traversal.py @@ -1,5 +1,7 @@ import numpy as np +import pytest from numpy.testing import assert_array_almost_equal +from scipy.sparse import csr_array from scipy.sparse.csgraph import (breadth_first_tree, depth_first_tree, csgraph_to_dense, csgraph_from_dense) @@ -66,3 +68,14 @@ def test_graph_depth_first_trivial_graph(): bfirst_test = depth_first_tree(csgraph, 0, directed) assert_array_almost_equal(csgraph_to_dense(bfirst_test), bfirst) + + +@pytest.mark.parametrize('directed', [True, False]) +@pytest.mark.parametrize('tree_func', [breadth_first_tree, depth_first_tree]) +def test_int64_indices(tree_func, directed): + # See https://github.com/scipy/scipy/issues/18716 + g = csr_array(([1], np.array([[0], [1]], dtype=np.int64)), shape=(2, 2)) + assert g.indices.dtype == np.int64 + tree = tree_func(g, 0, directed=directed) + assert_array_almost_equal(csgraph_to_dense(tree), [[0, 1], [0, 0]]) + diff --git a/scipy/special/_generate_pyx.py b/scipy/special/_generate_pyx.py index 7c8f487df8f6..2486a7ae7960 100644 --- a/scipy/special/_generate_pyx.py +++ b/scipy/special/_generate_pyx.py @@ -1232,6 +1232,7 @@ def get_declaration(ufunc, c_name, c_proto, cy_proto, header, # redeclare the function, so that the assumed # signature is checked at compile time new_name = f"{ufunc.cython_func_name(c_name)} \"{c_name}\"" + proto_h_filename = os.path.basename(proto_h_filename) defs.append(f'cdef extern from r"{proto_h_filename}":') defs.append(" cdef %s" % (cy_proto.replace('(*)', new_name))) defs_h.append(f'#include "{header}"') diff --git a/scipy/stats/_boost/include/func_defs.hpp b/scipy/stats/_boost/include/func_defs.hpp index 2ec998c09f8d..ff373d39ed3a 100644 --- a/scipy/stats/_boost/include/func_defs.hpp +++ b/scipy/stats/_boost/include/func_defs.hpp @@ -11,12 +11,6 @@ typedef boost::math::policies::policy< boost::math::policies::discrete_quantile< boost::math::policies::integer_round_up > > Policy; -// Run user_error function when evaluation_errors and overflow_errors are encountered -typedef boost::math::policies::policy< - boost::math::policies::evaluation_error, - boost::math::policies::overflow_error > user_error_policy; -BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS(user_error_policy) - // Raise a RuntimeWarning making users aware that something went wrong during // evaluation of the function, but return the best guess diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py index c21b0a710e2c..3a02b49a1f67 100644 --- a/scipy/stats/_continuous_distns.py +++ b/scipy/stats/_continuous_distns.py @@ -5275,7 +5275,7 @@ def _stats(self, a, b, moments='mv'): # Numerical improvements left to future enhancements. mu, mu2, g1, g2 = None, None, None, None - bn2 = b**-2 + bn2 = b**-2. expbn2 = np.exp(bn2) a_b = a / b diff --git a/scipy/stats/_covariance.py b/scipy/stats/_covariance.py index 79d5009f18f0..ec012c64da4d 100644 --- a/scipy/stats/_covariance.py +++ b/scipy/stats/_covariance.py @@ -579,7 +579,7 @@ def __init__(self, eigendecomposition): psuedo_reciprocals[i_zero] = 0 self._LP = eigenvectors * psuedo_reciprocals - self._LA = eigenvectors * np.sqrt(positive_eigenvalues) + self._LA = eigenvectors * np.sqrt(eigenvalues) self._rank = positive_eigenvalues.shape[-1] - i_zero.sum(axis=-1) self._w = eigenvalues self._v = eigenvectors diff --git a/scipy/stats/_mstats_basic.py b/scipy/stats/_mstats_basic.py index 2d10561d27ea..812da31cdc74 100644 --- a/scipy/stats/_mstats_basic.py +++ b/scipy/stats/_mstats_basic.py @@ -108,7 +108,7 @@ def argstoarray(*args): Notes ----- - `numpy.ma.row_stack` has identical behavior, but is called with a sequence + `numpy.ma.vstack` has identical behavior, but is called with a sequence of sequences. Examples diff --git a/scipy/stats/meson.build b/scipy/stats/meson.build index 6df5577a86aa..3591714fea5e 100644 --- a/scipy/stats/meson.build +++ b/scipy/stats/meson.build @@ -109,7 +109,16 @@ _stats_gen_pyx = custom_target('_stats_gen_pyx', ], input: '_generate_pyx.py', command: [py3, '@INPUT@', '-o', '@OUTDIR@'], - depends: _stats_pxd + depends: _stats_pxd, + depend_files: [ + '_boost/include/code_gen.py', + '_boost/include/gen_func_defs_pxd.py', + '_boost/include/_info.py', + # TODO: once setup.py is gone, remove this .pxd file from here and from the + # output of this custom_target, and stop copying it manually in + # code_gen.py. Instead, use fs.copyfile at the top of this file. + '_boost/include/templated_pyufunc.pxd', + ] ) # Only needed to establish a dependency; see comments in linalg/meson.build diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py index c9cb800a74c5..b94f37252ce1 100644 --- a/scipy/stats/tests/test_multivariate.py +++ b/scipy/stats/tests/test_multivariate.py @@ -150,6 +150,15 @@ def test_covariance(self, matrix_type, cov_type_name): if hasattr(cov_object, "_colorize") and "singular" not in matrix_type: assert_close(cov_object.colorize(res), x) + # gh-19197 reported that multivariate normal `rvs` produced incorrect + # results when a singular Covariance object was produce using + # `from_eigenvalues`. This was due to an issue in `colorize` with + # singular covariance matrices. Check this edge case, which is skipped + # in the previous tests. + if hasattr(cov_object, "_colorize"): + res = cov_object.colorize(np.eye(len(A))) + assert_close(res.T @ res, A) + @pytest.mark.parametrize("size", [None, tuple(), 1, (2, 4, 3)]) @pytest.mark.parametrize("matrix_type", list(_matrices)) @pytest.mark.parametrize("cov_type_name", _all_covariance_types) @@ -244,6 +253,24 @@ def test_gh9942(self): ref = multivariate_normal.rvs(mean, cov, random_state=rng2) assert_equal(res, ref) + def test_gh19197(self): + # gh-19197 reported that multivariate normal `rvs` produced incorrect + # results when a singular Covariance object was produce using + # `from_eigenvalues`. Check that this specific issue is resolved; + # a more general test is included in `test_covariance`. + mean = np.ones(2) + cov = Covariance.from_eigendecomposition((np.zeros(2), np.eye(2))) + dist = scipy.stats.multivariate_normal(mean=mean, cov=cov) + rvs = dist.rvs(size=None) + assert_equal(rvs, mean) + + cov = scipy.stats.Covariance.from_eigendecomposition( + (np.array([1., 0.]), np.array([[1., 0.], [0., 400.]]))) + dist = scipy.stats.multivariate_normal(mean=mean, cov=cov) + rvs = dist.rvs(size=None) + assert rvs[0] != mean[0] + assert rvs[1] == mean[1] + def _random_covariance(dim, evals, rng, singular=False): # Generates random covariance matrix with dimensionality `dim` and diff --git a/tools/wheels/LICENSE_linux.txt b/tools/wheels/LICENSE_linux.txt index e4f810cd9cc1..eb1c901d66b9 100644 --- a/tools/wheels/LICENSE_linux.txt +++ b/tools/wheels/LICENSE_linux.txt @@ -5,10 +5,10 @@ This binary distribution of SciPy also bundles the following software: Name: OpenBLAS -Files: .libs/libopenb*.so +Files: scipy.libs/libopenblas*.so Description: bundled as a dynamically linked library -Availability: https://github.com/xianyi/OpenBLAS/ -License: 3-clause BSD +Availability: https://github.com/OpenMathLib/OpenBLAS/ +License: BSD-3-Clause-Attribution Copyright (c) 2011-2014, The OpenBLAS Project All rights reserved. @@ -41,10 +41,10 @@ License: 3-clause BSD Name: LAPACK -Files: .libs/libopenb*.so +Files: scipy.libs/libopenblas*.so Description: bundled in OpenBLAS -Availability: https://github.com/xianyi/OpenBLAS/ -License 3-clause BSD +Availability: https://github.com/OpenMathLib/OpenBLAS/ +License: BSD-3-Clause-Attribution Copyright (c) 1992-2013 The University of Tennessee and The University of Tennessee Research Foundation. All rights reserved. @@ -96,10 +96,10 @@ License 3-clause BSD Name: GCC runtime library -Files: .libs/libgfortran*.so +Files: scipy.libs/libgfortran*.so Description: dynamically linked to files compiled with gcc -Availability: https://gcc.gnu.org/viewcvs/gcc/ -License: GPLv3 + runtime exception +Availability: https://gcc.gnu.org/git/?p=gcc.git;a=tree;f=libgfortran +License: GPL-3.0-with-GCC-exception Copyright (C) 2002-2017 Free Software Foundation, Inc. Libgfortran is free software; you can redistribute it and/or modify @@ -878,3 +878,26 @@ may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read . + + +Name: libquadmath +Files: scipy.libs/libquadmath*.so +Description: dynamically linked to files compiled with gcc +Availability: https://gcc.gnu.org/git/?p=gcc.git;a=tree;f=libquadmath +License: LGPL-2.1-or-later + + GCC Quad-Precision Math Library + Copyright (C) 2010-2019 Free Software Foundation, Inc. + Written by Francois-Xavier Coudert + + This file is part of the libquadmath library. + Libquadmath is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + Libquadmath is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html diff --git a/tools/wheels/LICENSE_osx.txt b/tools/wheels/LICENSE_osx.txt index 63a5497ed7c6..447dbdb1d7e0 100644 --- a/tools/wheels/LICENSE_osx.txt +++ b/tools/wheels/LICENSE_osx.txt @@ -4,11 +4,102 @@ This binary distribution of SciPy also bundles the following software: +Name: OpenBLAS +Files: scipy/.dylibs/libopenblas*.so +Description: bundled as a dynamically linked library +Availability: https://github.com/OpenMathLib/OpenBLAS/ +License: BSD-3-Clause-Attribution + Copyright (c) 2011-2014, The OpenBLAS Project + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. 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. + 3. Neither the name of the OpenBLAS project nor the names of + its contributors 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 OWNER 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. + + +Name: LAPACK +Files: scipy/.dylibs/libopenblas*.so +Description: bundled in OpenBLAS +Availability: https://github.com/OpenMathLib/OpenBLAS/ +License: BSD-3-Clause-Attribution + Copyright (c) 1992-2013 The University of Tennessee and The University + of Tennessee Research Foundation. All rights + reserved. + Copyright (c) 2000-2013 The University of California Berkeley. All + rights reserved. + Copyright (c) 2006-2013 The University of Colorado Denver. All rights + reserved. + + $COPYRIGHT$ + + Additional copyrights may follow + + $HEADER$ + + 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 listed + in this license in the documentation and/or other materials + provided with the distribution. + + - Neither the name of the copyright holders nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + + The copyright holders provide no reassurances that the source code + provided does not infringe any patent, copyright, or any other + intellectual property rights of third parties. The copyright holders + disclaim any liability to any recipient for claims brought against + recipient by any third party for infringement of that parties + intellectual property rights. + + 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 + OWNER 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. + + Name: GCC runtime library -Files: .dylibs/* +Files: scipy/.dylibs/libgfortran*, scipy/.dylibs/libgcc* Description: dynamically linked to files compiled with gcc -Availability: https://gcc.gnu.org/viewcvs/gcc/ -License: GPLv3 + runtime exception +Availability: https://gcc.gnu.org/git/?p=gcc.git;a=tree;f=libgfortran +License: GPL-3.0-with-GCC-exception Copyright (C) 2002-2017 Free Software Foundation, Inc. Libgfortran is free software; you can redistribute it and/or modify @@ -787,3 +878,26 @@ may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read . + + +Name: libquadmath +Files: scipy/.dylibs/libquadmath*.so +Description: dynamically linked to files compiled with gcc +Availability: https://gcc.gnu.org/git/?p=gcc.git;a=tree;f=libquadmath +License: LGPL-2.1-or-later + + GCC Quad-Precision Math Library + Copyright (C) 2010-2019 Free Software Foundation, Inc. + Written by Francois-Xavier Coudert + + This file is part of the libquadmath library. + Libquadmath is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + Libquadmath is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html diff --git a/tools/wheels/LICENSE_win32.txt b/tools/wheels/LICENSE_win32.txt index c93727261e45..3d7e17e516b1 100644 --- a/tools/wheels/LICENSE_win32.txt +++ b/tools/wheels/LICENSE_win32.txt @@ -5,10 +5,10 @@ This binary distribution of SciPy also bundles the following software: Name: OpenBLAS -Files: extra-dll\libopenb*.dll +Files: scipy.libs\libopenblas*.dll Description: bundled as a dynamically linked library -Availability: https://github.com/xianyi/OpenBLAS/ -License: 3-clause BSD +Availability: https://github.com/OpenMathLib/OpenBLAS/ +License: BSD-3-Clause-Attribution Copyright (c) 2011-2014, The OpenBLAS Project All rights reserved. @@ -41,10 +41,10 @@ License: 3-clause BSD Name: LAPACK -Files: extra-dll\libopenb*.dll +Files: scipy.libs\libopenblas*.dll Description: bundled in OpenBLAS -Availability: https://github.com/xianyi/OpenBLAS/ -License 3-clause BSD +Availability: https://github.com/OpenMathLib/OpenBLAS/ +License: BSD-3-Clause-Attribution Copyright (c) 1992-2013 The University of Tennessee and The University of Tennessee Research Foundation. All rights reserved. @@ -96,10 +96,10 @@ License 3-clause BSD Name: GCC runtime library -Files: extra-dll\*.dll -Description: statically linked, in DLL files compiled with gfortran only -Availability: https://gcc.gnu.org/viewcvs/gcc/ -License: GPLv3 + runtime exception +Files: scipy.libs\libopenblas*.dll +Description: statically linked to files compiled with gcc +Availability: https://gcc.gnu.org/git/?p=gcc.git;a=tree;f=libgfortran +License: GPL-3.0-with-GCC-exception Copyright (C) 2002-2017 Free Software Foundation, Inc. Libgfortran is free software; you can redistribute it and/or modify @@ -879,3 +879,26 @@ may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read . + + +Name: libquadmath +Files: scipy.libs\libopenblas*.dll +Description: statically linked to files compiled with gcc +Availability: https://gcc.gnu.org/git/?p=gcc.git;a=tree;f=libquadmath +License: LGPL-2.1-or-later + + GCC Quad-Precision Math Library + Copyright (C) 2010-2019 Free Software Foundation, Inc. + Written by Francois-Xavier Coudert + + This file is part of the libquadmath library. + Libquadmath is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + Libquadmath is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html