From e4c90de548903fc2484a6b462aa2416ca0ed9041 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Wed, 14 Apr 2021 17:39:21 +0200 Subject: [PATCH 01/48] prevent implicit np.array -> scalar conversion --- scipy/io/_fortran.py | 2 +- scipy/linalg/tests/test_basic.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scipy/io/_fortran.py b/scipy/io/_fortran.py index de139b9a8f4c..193f6776b1f3 100644 --- a/scipy/io/_fortran.py +++ b/scipy/io/_fortran.py @@ -131,7 +131,7 @@ def _read_size(self, eof_ok=False): elif len(b) < n: raise FortranFormattingError( "End of file in the middle of the record size") - return int(np.frombuffer(b, dtype=self._header_dtype, count=1)) + return int(np.frombuffer(b, dtype=self._header_dtype, count=1)[0]) def write_record(self, *items): """ diff --git a/scipy/linalg/tests/test_basic.py b/scipy/linalg/tests/test_basic.py index 8ff7e7b05c66..831ceab8eac1 100644 --- a/scipy/linalg/tests/test_basic.py +++ b/scipy/linalg/tests/test_basic.py @@ -1566,7 +1566,7 @@ def test_scaling(self): # Pre/post LAPACK 3.5.0 gives the same result up to an offset # since in each case col norm is x1000 greater and # 1000 / 32 ~= 1 * 32 hence balanced with 2 ** 5. - assert_allclose(int(np.diff(np.log2(np.diag(y)))), 5) + assert_allclose(int(np.diff(np.log2(np.diag(y)))[0]), 5) def test_scaling_order(self): A = np.array([[1, 0, 1e-4], [1, 1, 1e-2], [1e4, 1e2, 1]]) @@ -1576,7 +1576,7 @@ def test_scaling_order(self): def test_separate(self): _, (y, z) = matrix_balance(np.array([[1000, 1], [1000, 0]]), separate=1) - assert_equal(int(np.diff(np.log2(y))), 5) + assert_equal(int(np.diff(np.log2(y))[0]), 5) assert_allclose(z, np.arange(2)) def test_permutation(self): From 5446bfd13a50a78d94f6dee964022c09242d8d05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Wed, 14 Apr 2021 17:54:02 +0200 Subject: [PATCH 02/48] better f(x) handling in slsqp() --- scipy/optimize/slsqp.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/scipy/optimize/slsqp.py b/scipy/optimize/slsqp.py index ecd1a0c8b344..ec7988ecc398 100644 --- a/scipy/optimize/slsqp.py +++ b/scipy/optimize/slsqp.py @@ -413,11 +413,13 @@ def cjac(x, *args): # mode is zero on entry, so call objective, constraints and gradients # there should be no func evaluations here because it's cached from # ScalarFunction - fx = wrapped_fun(x) - try: - fx = float(np.asarray(fx)) - except (TypeError, ValueError) as e: - raise ValueError("Objective function must return a scalar") from e + fx = np.asarray(wrapped_fun(x)) + # Ideally, we'd like to a have a true scalar returned from f(x). For + # backwards-compatility, also allow np.array([1.3]), np.array([[1.3]]) etc. + if fx.size != 1: + raise ValueError("Objective function must return a scalar") + fx = fx.flat[0] + g = append(wrapped_grad(x), 0.0) c = _eval_constraint(x, cons) a = _eval_con_normals(x, cons, la, n, m, meq, mieq) From 66d130eb183e89081fa1e3eb7b418c95cb1e1003 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Wed, 14 Apr 2021 18:15:14 +0200 Subject: [PATCH 03/48] nelder-mead: some return value handling --- scipy/optimize/optimize.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index 0c09027d6aa3..e2b8c7efa740 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -315,11 +315,11 @@ def rosen(x): >>> X = 0.1 * np.arange(10) >>> rosen(X) 76.56 - + For higher-dimensional input ``rosen`` broadcasts. In the following example, we use this to plot a 2D landscape. Note that ``rosen_hess`` does not broadcast in this manner. - + >>> import matplotlib.pyplot as plt >>> from mpl_toolkits.mplot3d import Axes3D >>> x = np.linspace(-1, 1, 50) @@ -744,10 +744,13 @@ def _minimize_neldermead(func, x0, args=(), callback=None, sim = np.clip(sim, lower_bound, upper_bound) one2np1 = list(range(1, N + 1)) - fsim = np.empty((N + 1,), float) - for k in range(N + 1): - fsim[k] = func(sim[k]) + fsim = np.array([func(sim[k]) for k in range(N + 1)], dtype=float) + if fsim.size != N + 1: + raise ValueError("Objective function must return a scalar") + # Ideally, we'd like to a have a true scalar returned from f(x). For + # backwards-compatility, also allow np.array([1.3]), np.array([[1.3]]) etc. + fsim = fsim.reshape(N + 1) ind = np.argsort(fsim) fsim = np.take(fsim, ind, 0) @@ -772,7 +775,11 @@ def _minimize_neldermead(func, x0, args=(), callback=None, xe = (1 + rho * chi) * xbar - rho * chi * sim[-1] if bounds is not None: xe = np.clip(xe, lower_bound, upper_bound) - fxe = func(xe) + + fxe = np.asarray(func(xe)) + if fxe.size != 1: + raise ValueError("Objective function must return a scalar") + fxe = fxe.flat[0] if fxe < fxr: sim[-1] = xe From 18ce26732a4303594b573c35253920dd30ae7baf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Wed, 14 Apr 2021 18:17:58 +0200 Subject: [PATCH 04/48] avoid implicit array conversion in _trustregion_constr/report.py --- scipy/optimize/_trustregion_constr/report.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/optimize/_trustregion_constr/report.py b/scipy/optimize/_trustregion_constr/report.py index 8d35dab2792d..78cc5ccae5d6 100644 --- a/scipy/optimize/_trustregion_constr/report.py +++ b/scipy/optimize/_trustregion_constr/report.py @@ -23,7 +23,7 @@ def print_iteration(cls, *args): # trust-constr typically provides a length 1 array. We have to coerce # it to a float, otherwise the string format doesn't work. args = list(args) - args[3] = float(args[3]) + args[3] = args[3][0] iteration_format = ["{{:{}}}".format(x) for x in cls.ITERATION_FORMATS] fmt = "|" + "|".join(iteration_format) + "|" From d6967d069b23ab07c876d557a8108675f1a59e68 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Wed, 14 Apr 2021 20:29:10 +0200 Subject: [PATCH 05/48] flat[0] -> item() --- scipy/optimize/optimize.py | 2 +- scipy/optimize/slsqp.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index e2b8c7efa740..7d6b0aab0a82 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -779,7 +779,7 @@ def _minimize_neldermead(func, x0, args=(), callback=None, fxe = np.asarray(func(xe)) if fxe.size != 1: raise ValueError("Objective function must return a scalar") - fxe = fxe.flat[0] + fxe = fxe.item() if fxe < fxr: sim[-1] = xe diff --git a/scipy/optimize/slsqp.py b/scipy/optimize/slsqp.py index ec7988ecc398..00f029fbc5ee 100644 --- a/scipy/optimize/slsqp.py +++ b/scipy/optimize/slsqp.py @@ -418,7 +418,7 @@ def cjac(x, *args): # backwards-compatility, also allow np.array([1.3]), np.array([[1.3]]) etc. if fx.size != 1: raise ValueError("Objective function must return a scalar") - fx = fx.flat[0] + fx = fx.item() g = append(wrapped_grad(x), 0.0) c = _eval_constraint(x, cons) From e665bb85b6b0eb2514e3e4e4c61ffe7e3c2aa76e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Wed, 14 Apr 2021 20:30:27 +0200 Subject: [PATCH 06/48] comment typos --- scipy/optimize/optimize.py | 2 +- scipy/optimize/slsqp.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index 7d6b0aab0a82..e8df0b31ae18 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -749,7 +749,7 @@ def _minimize_neldermead(func, x0, args=(), callback=None, if fsim.size != N + 1: raise ValueError("Objective function must return a scalar") # Ideally, we'd like to a have a true scalar returned from f(x). For - # backwards-compatility, also allow np.array([1.3]), np.array([[1.3]]) etc. + # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. fsim = fsim.reshape(N + 1) ind = np.argsort(fsim) diff --git a/scipy/optimize/slsqp.py b/scipy/optimize/slsqp.py index 00f029fbc5ee..c8a32fb8d41b 100644 --- a/scipy/optimize/slsqp.py +++ b/scipy/optimize/slsqp.py @@ -415,7 +415,7 @@ def cjac(x, *args): # ScalarFunction fx = np.asarray(wrapped_fun(x)) # Ideally, we'd like to a have a true scalar returned from f(x). For - # backwards-compatility, also allow np.array([1.3]), np.array([[1.3]]) etc. + # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. if fx.size != 1: raise ValueError("Objective function must return a scalar") fx = fx.item() From 275d33bd7645a4c7f21e22fb7357a158c0513705 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 08:51:01 +0200 Subject: [PATCH 07/48] Update scipy/linalg/tests/test_basic.py Co-authored-by: peterbell10 --- scipy/linalg/tests/test_basic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/linalg/tests/test_basic.py b/scipy/linalg/tests/test_basic.py index 831ceab8eac1..24c761679ee6 100644 --- a/scipy/linalg/tests/test_basic.py +++ b/scipy/linalg/tests/test_basic.py @@ -1576,7 +1576,7 @@ def test_scaling_order(self): def test_separate(self): _, (y, z) = matrix_balance(np.array([[1000, 1], [1000, 0]]), separate=1) - assert_equal(int(np.diff(np.log2(y))[0]), 5) + assert_equal(np.diff(np.log2(y)), [5]) assert_allclose(z, np.arange(2)) def test_permutation(self): From e59c062949ae448016a41298574b59fdac59a51f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 08:52:21 +0200 Subject: [PATCH 08/48] test_basic: allclose(int(x[0]), y) -> allclose(x, [y]) --- scipy/linalg/tests/test_basic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/linalg/tests/test_basic.py b/scipy/linalg/tests/test_basic.py index 24c761679ee6..1da183f65b15 100644 --- a/scipy/linalg/tests/test_basic.py +++ b/scipy/linalg/tests/test_basic.py @@ -1566,7 +1566,7 @@ def test_scaling(self): # Pre/post LAPACK 3.5.0 gives the same result up to an offset # since in each case col norm is x1000 greater and # 1000 / 32 ~= 1 * 32 hence balanced with 2 ** 5. - assert_allclose(int(np.diff(np.log2(np.diag(y)))[0]), 5) + assert_allclose(np.diff(np.log2(np.diag(y))), [5]) def test_scaling_order(self): A = np.array([[1, 0, 1e-4], [1, 1, 1e-2], [1e4, 1e2, 1]]) From 14e4391060f5a83f2a977bf657fa10c26a331453 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 09:04:54 +0200 Subject: [PATCH 09/48] report: more robust float conversion --- scipy/optimize/_trustregion_constr/report.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/scipy/optimize/_trustregion_constr/report.py b/scipy/optimize/_trustregion_constr/report.py index 78cc5ccae5d6..e2b335ffab8f 100644 --- a/scipy/optimize/_trustregion_constr/report.py +++ b/scipy/optimize/_trustregion_constr/report.py @@ -1,6 +1,7 @@ """Progress report printers.""" from __future__ import annotations +import numpy as np from typing import List class ReportBase: @@ -19,11 +20,11 @@ def print_header(cls): @classmethod def print_iteration(cls, *args): - # args[3] is obj func. It should really be a float. However, - # trust-constr typically provides a length 1 array. We have to coerce - # it to a float, otherwise the string format doesn't work. + # args[3] is the value of the user-defined objective function. It should really + # be a float, but often a length-1 array is provided. We have to coerce it to a + # float, otherwise the string format doesn't work. args = list(args) - args[3] = args[3][0] + args[3] = np.asarray(args[3]).item() iteration_format = ["{{:{}}}".format(x) for x in cls.ITERATION_FORMATS] fmt = "|" + "|".join(iteration_format) + "|" From ef8bae522b1c32945317c3aad390620a8d50e98f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 10:33:39 +0200 Subject: [PATCH 10/48] move scalar check to function wrapper --- scipy/optimize/optimize.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index e8df0b31ae18..92a8b01a9f81 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -454,14 +454,24 @@ def rosen_hess_prod(x, p): def _wrap_function(function, args): # wraps a minimizer function to count number of evaluations # and to easily provide an args kwd. - # A copy of x is sent to the user function (gh13740) ncalls = [0] if function is None: return ncalls, None def function_wrapper(x, *wrapper_args): ncalls[0] += 1 - return function(np.copy(x), *(wrapper_args + args)) + # A copy of x is sent to the user function (gh13740) + fx = function(np.copy(x), *(wrapper_args + args)) + # Ideally, we'd like to a have a true scalar returned from f(x). For + # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. + print() + print(x, fx) + try: + fx = np.asarray(fx).item() + except (TypeError, ValueError) as e: + raise ValueError("Objective function must return a scalar") from e + print(fx) + return fx return ncalls, function_wrapper @@ -746,11 +756,6 @@ def _minimize_neldermead(func, x0, args=(), callback=None, one2np1 = list(range(1, N + 1)) fsim = np.array([func(sim[k]) for k in range(N + 1)], dtype=float) - if fsim.size != N + 1: - raise ValueError("Objective function must return a scalar") - # Ideally, we'd like to a have a true scalar returned from f(x). For - # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. - fsim = fsim.reshape(N + 1) ind = np.argsort(fsim) fsim = np.take(fsim, ind, 0) @@ -776,10 +781,7 @@ def _minimize_neldermead(func, x0, args=(), callback=None, if bounds is not None: xe = np.clip(xe, lower_bound, upper_bound) - fxe = np.asarray(func(xe)) - if fxe.size != 1: - raise ValueError("Objective function must return a scalar") - fxe = fxe.item() + fxe = func(xe) if fxe < fxr: sim[-1] = xe From 2fcb4a0632f5f45f78c7503ea9bbb592e18deebf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 10:59:57 +0200 Subject: [PATCH 11/48] move scalar check into ScalarFunction --- scipy/optimize/_differentiable_functions.py | 12 +++++++++++- scipy/optimize/optimize.py | 16 ---------------- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/scipy/optimize/_differentiable_functions.py b/scipy/optimize/_differentiable_functions.py index 9550b7eefda8..57e8a3127ddc 100644 --- a/scipy/optimize/_differentiable_functions.py +++ b/scipy/optimize/_differentiable_functions.py @@ -131,7 +131,17 @@ def fun_wrapped(x): # Send a copy because the user may overwrite it. # Overwriting results in undefined behaviour because # fun(self.x) will change self.x, with the two no longer linked. - return fun(np.copy(x), *args) + fx = fun(np.copy(x), *args) + # Make sure the function returns a true scalar + if not np.isscalar(fx): + try: + fx = np.asarray(fx).item() + except (TypeError, ValueError) as e: + raise ValueError( + "The user-provided objective function " + "must return a scalar value." + ) from e + return fx def update_fun(): self.f = fun_wrapped(self.x) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index 92a8b01a9f81..a0085f0d6330 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -1187,14 +1187,6 @@ def _minimize_bfgs(fun, x0, args=(), jac=None, callback=None, old_fval = f(x0) gfk = myfprime(x0) - if not np.isscalar(old_fval): - try: - old_fval = old_fval.item() - except (ValueError, AttributeError) as e: - raise ValueError("The user-provided " - "objective function must " - "return a scalar value.") from e - k = 0 N = len(x0) I = np.eye(N, dtype=int) @@ -1508,14 +1500,6 @@ def _minimize_cg(fun, x0, args=(), jac=None, callback=None, old_fval = f(x0) gfk = myfprime(x0) - if not np.isscalar(old_fval): - try: - old_fval = old_fval.item() - except (ValueError, AttributeError) as e: - raise ValueError("The user-provided " - "objective function must " - "return a scalar value.") from e - k = 0 xk = x0 # Sets the initial step guess to dx ~ 1 From 899e50a4c59e81e2cc984bd9fddc0bbd6564caee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 11:31:25 +0200 Subject: [PATCH 12/48] separate function wrapper for hessian --- scipy/optimize/_trustregion.py | 17 ++++++++++++++++- scipy/optimize/optimize.py | 9 +++------ 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/scipy/optimize/_trustregion.py b/scipy/optimize/_trustregion.py index 79a5fc31203a..b390419cd0de 100644 --- a/scipy/optimize/_trustregion.py +++ b/scipy/optimize/_trustregion.py @@ -3,13 +3,28 @@ import numpy as np import scipy.linalg -from .optimize import (_check_unknown_options, _wrap_function, _status_message, +from .optimize import (_check_unknown_options, _status_message, OptimizeResult, _prepare_scalar_function) from scipy.optimize._hessian_update_strategy import HessianUpdateStrategy from scipy.optimize._differentiable_functions import FD_METHODS __all__ = [] +def _wrap_function(function, args): + # wraps a minimizer function to count number of evaluations + # and to easily provide an args kwd. + ncalls = [0] + if function is None: + return ncalls, None + + def function_wrapper(x, *wrapper_args): + ncalls[0] += 1 + # A copy of x is sent to the user function (gh13740) + return function(np.copy(x), *(wrapper_args + args)) + + return ncalls, function_wrapper + + class BaseQuadraticSubproblem: """ Base/abstract class defining the quadratic model for trust-region diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index a0085f0d6330..13551363ea19 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -451,7 +451,7 @@ def rosen_hess_prod(x, p): return Hp -def _wrap_function(function, args): +def _wrap_scalar_function(function, args): # wraps a minimizer function to count number of evaluations # and to easily provide an args kwd. ncalls = [0] @@ -464,13 +464,10 @@ def function_wrapper(x, *wrapper_args): fx = function(np.copy(x), *(wrapper_args + args)) # Ideally, we'd like to a have a true scalar returned from f(x). For # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. - print() - print(x, fx) try: fx = np.asarray(fx).item() except (TypeError, ValueError) as e: raise ValueError("Objective function must return a scalar") from e - print(fx) return fx return ncalls, function_wrapper @@ -679,7 +676,7 @@ def _minimize_neldermead(func, x0, args=(), callback=None, maxfun = maxfev retall = return_all - fcalls, func = _wrap_function(func, args) + fcalls, func = _wrap_scalar_function(func, args) if adaptive: dim = float(len(x0)) @@ -2937,7 +2934,7 @@ def _minimize_powell(func, x0, args=(), callback=None, bounds=None, retall = return_all # we need to use a mutable object here that we can update in the # wrapper function - fcalls, func = _wrap_function(func, args) + fcalls, func = _wrap_scalar_function(func, args) x = asarray(x0).flatten() if retall: allvecs = [x] From 242361e5c4c1f7c129882ff5015309a3f17bb85e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 11:35:34 +0200 Subject: [PATCH 13/48] remove workaround in print_iteration --- scipy/optimize/_trustregion_constr/report.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/scipy/optimize/_trustregion_constr/report.py b/scipy/optimize/_trustregion_constr/report.py index e2b335ffab8f..1631bf21e504 100644 --- a/scipy/optimize/_trustregion_constr/report.py +++ b/scipy/optimize/_trustregion_constr/report.py @@ -1,7 +1,6 @@ """Progress report printers.""" from __future__ import annotations -import numpy as np from typing import List class ReportBase: @@ -20,12 +19,6 @@ def print_header(cls): @classmethod def print_iteration(cls, *args): - # args[3] is the value of the user-defined objective function. It should really - # be a float, but often a length-1 array is provided. We have to coerce it to a - # float, otherwise the string format doesn't work. - args = list(args) - args[3] = np.asarray(args[3]).item() - iteration_format = ["{{:{}}}".format(x) for x in cls.ITERATION_FORMATS] fmt = "|" + "|".join(iteration_format) + "|" print(fmt.format(*args)) From 567613c8ad6d898c2068a742fbe43076a6ae3015 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 11:40:05 +0200 Subject: [PATCH 14/48] less code modification --- scipy/optimize/optimize.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index 13551363ea19..d63313104702 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -752,7 +752,9 @@ def _minimize_neldermead(func, x0, args=(), callback=None, one2np1 = list(range(1, N + 1)) - fsim = np.array([func(sim[k]) for k in range(N + 1)], dtype=float) + fsim = np.empty((N + 1,), float) + for k in range(N + 1): + fsim[k] = func(sim[k]) ind = np.argsort(fsim) fsim = np.take(fsim, ind, 0) @@ -777,7 +779,6 @@ def _minimize_neldermead(func, x0, args=(), callback=None, xe = (1 + rho * chi) * xbar - rho * chi * sim[-1] if bounds is not None: xe = np.clip(xe, lower_bound, upper_bound) - fxe = func(xe) if fxe < fxr: From 47f46e1ad51222f7312e059a1f2ef7675564fcb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 11:41:48 +0200 Subject: [PATCH 15/48] move newline --- scipy/optimize/optimize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index d63313104702..926fb10369bf 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -751,8 +751,8 @@ def _minimize_neldermead(func, x0, args=(), callback=None, sim = np.clip(sim, lower_bound, upper_bound) one2np1 = list(range(1, N + 1)) - fsim = np.empty((N + 1,), float) + for k in range(N + 1): fsim[k] = func(sim[k]) From 31eba20cca12f101a3e26f24d526473f2ceec4be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 11:45:51 +0200 Subject: [PATCH 16/48] remove redundant workaround for slsqp() --- scipy/optimize/slsqp.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/scipy/optimize/slsqp.py b/scipy/optimize/slsqp.py index c8a32fb8d41b..669775698576 100644 --- a/scipy/optimize/slsqp.py +++ b/scipy/optimize/slsqp.py @@ -413,12 +413,7 @@ def cjac(x, *args): # mode is zero on entry, so call objective, constraints and gradients # there should be no func evaluations here because it's cached from # ScalarFunction - fx = np.asarray(wrapped_fun(x)) - # Ideally, we'd like to a have a true scalar returned from f(x). For - # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. - if fx.size != 1: - raise ValueError("Objective function must return a scalar") - fx = fx.item() + fx = wrapped_fun(x) g = append(wrapped_grad(x), 0.0) c = _eval_constraint(x, cons) From b4ef943ece9196f6ac25b3b45f7a3bd404c87170 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 11:46:35 +0200 Subject: [PATCH 17/48] remove newline --- scipy/optimize/slsqp.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scipy/optimize/slsqp.py b/scipy/optimize/slsqp.py index 669775698576..b068c1e77fdc 100644 --- a/scipy/optimize/slsqp.py +++ b/scipy/optimize/slsqp.py @@ -414,7 +414,6 @@ def cjac(x, *args): # there should be no func evaluations here because it's cached from # ScalarFunction fx = wrapped_fun(x) - g = append(wrapped_grad(x), 0.0) c = _eval_constraint(x, cons) a = _eval_con_normals(x, cons, la, n, m, meq, mieq) From 08c7434d1dc395182364a1a2dea49dedb0cd127c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 15 Apr 2021 12:14:04 +0200 Subject: [PATCH 18/48] dual_annealing: get scalar from visit_fn() --- scipy/optimize/_dual_annealing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/optimize/_dual_annealing.py b/scipy/optimize/_dual_annealing.py index bc3db4eba5ed..89c4a37c0ef6 100644 --- a/scipy/optimize/_dual_annealing.py +++ b/scipy/optimize/_dual_annealing.py @@ -92,7 +92,7 @@ def visiting(self, x, step, temperature): # Changing only one coordinate at a time based on strategy # chain step x_visit = np.copy(x) - visit = self.visit_fn(temperature, 1) + visit = self.visit_fn(temperature, 1)[0] if visit > self.TAIL_LIMIT: visit = self.TAIL_LIMIT * self.rand_gen.uniform() elif visit < -self.TAIL_LIMIT: From bedbc5037d20e2fd9e14d9879e87d376b32fa390 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 27 May 2021 17:35:45 +0200 Subject: [PATCH 19/48] optimize: change error text Co-authored-by: peterbell10 --- scipy/optimize/optimize.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index 926fb10369bf..eeb294af41c3 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -467,7 +467,8 @@ def function_wrapper(x, *wrapper_args): try: fx = np.asarray(fx).item() except (TypeError, ValueError) as e: - raise ValueError("Objective function must return a scalar") from e + raise ValueError("The user-provided objective function " + "must return a scalar value.") from e return fx return ncalls, function_wrapper From d04da6c2ce604d7f9db9ae1f9b8eec65474415a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 27 May 2021 17:39:10 +0200 Subject: [PATCH 20/48] optimiztion: short-circuit value check --- scipy/optimize/optimize.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index eeb294af41c3..0617a9427223 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -464,11 +464,12 @@ def function_wrapper(x, *wrapper_args): fx = function(np.copy(x), *(wrapper_args + args)) # Ideally, we'd like to a have a true scalar returned from f(x). For # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc. - try: - fx = np.asarray(fx).item() - except (TypeError, ValueError) as e: - raise ValueError("The user-provided objective function " - "must return a scalar value.") from e + if not np.isscalar(fx): + try: + fx = np.asarray(fx).item() + except (TypeError, ValueError) as e: + raise ValueError("The user-provided objective function " + "must return a scalar value.") from e return fx return ncalls, function_wrapper From b868eb7c117ebcc9d3f4bedfe92c6c27516f7d89 Mon Sep 17 00:00:00 2001 From: Itrimel Date: Wed, 9 Jun 2021 09:32:26 +0200 Subject: [PATCH 21/48] BUG: Fix Nelder-Mead logic when using a non-1D x0 and adapative --- scipy/optimize/optimize.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scipy/optimize/optimize.py b/scipy/optimize/optimize.py index 7cfd8ab07645..aedad176a877 100644 --- a/scipy/optimize/optimize.py +++ b/scipy/optimize/optimize.py @@ -671,6 +671,8 @@ def _minimize_neldermead(func, x0, args=(), callback=None, fcalls, func = _wrap_function(func, args) + x0 = asfarray(x0).flatten() + if adaptive: dim = float(len(x0)) rho = 1 @@ -686,8 +688,6 @@ def _minimize_neldermead(func, x0, args=(), callback=None, nonzdelt = 0.05 zdelt = 0.00025 - x0 = asfarray(x0).flatten() - if bounds is not None: lower_bound, upper_bound = bounds.lb, bounds.ub # check bounds From aefa2e02113d8af023a346a977e9aa21e1683375 Mon Sep 17 00:00:00 2001 From: Bas van Beek <43369155+BvB93@users.noreply.github.com> Date: Fri, 18 Jun 2021 20:07:54 +0200 Subject: [PATCH 22/48] BUG: Fixed an issue wherein `SphericalVoronoi` could raise for certain array-likes --- scipy/spatial/_spherical_voronoi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/spatial/_spherical_voronoi.py b/scipy/spatial/_spherical_voronoi.py index 4d7af19a01a3..d88f476233d4 100644 --- a/scipy/spatial/_spherical_voronoi.py +++ b/scipy/spatial/_spherical_voronoi.py @@ -176,7 +176,7 @@ def __init__(self, points, radius=1, center=None, threshold=1e-06): self.radius = float(radius) self.points = np.array(points).astype(np.double) - self._dim = len(points[0]) + self._dim = self.points.shape[1] if center is None: self.center = np.zeros(self._dim) else: From 919f0de2803477d07335c17ee134c526d0fd39ce Mon Sep 17 00:00:00 2001 From: Bas van Beek <43369155+BvB93@users.noreply.github.com> Date: Fri, 18 Jun 2021 20:12:05 +0200 Subject: [PATCH 23/48] TST: Add a `SphericalVoronoi` test for a non-sequence-based array-like --- scipy/spatial/tests/test_spherical_voronoi.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/scipy/spatial/tests/test_spherical_voronoi.py b/scipy/spatial/tests/test_spherical_voronoi.py index 356d03834cae..10729ebc4e9c 100644 --- a/scipy/spatial/tests/test_spherical_voronoi.py +++ b/scipy/spatial/tests/test_spherical_voronoi.py @@ -128,6 +128,11 @@ def test_constructor(self): assert_array_equal(s4.center, center) assert_equal(s4.radius, radius) + # Test a non-sequence/-ndarray based array-like + s5 = SphericalVoronoi(memoryview(self.points)) + assert_array_equal(s5.center, np.array([0, 0, 0])) + assert_equal(s5.radius, 1) + def test_vertices_regions_translation_invariance(self): sv_origin = SphericalVoronoi(self.points) center = np.array([1, 1, 1]) From b3dc4b071cd38034372ae3784e0816b01003a0c8 Mon Sep 17 00:00:00 2001 From: Bas van Beek <43369155+BvB93@users.noreply.github.com> Date: Fri, 18 Jun 2021 23:45:41 +0200 Subject: [PATCH 24/48] TYP: Silence a mypy error (false positive) Related to the fact that `ndarray` is not considered a "valid" implementation of the buffer protocol due to typeshed's use of an inexhaustive union --- scipy/spatial/tests/test_spherical_voronoi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/spatial/tests/test_spherical_voronoi.py b/scipy/spatial/tests/test_spherical_voronoi.py index 10729ebc4e9c..d0277d6d4589 100644 --- a/scipy/spatial/tests/test_spherical_voronoi.py +++ b/scipy/spatial/tests/test_spherical_voronoi.py @@ -129,7 +129,7 @@ def test_constructor(self): assert_equal(s4.radius, radius) # Test a non-sequence/-ndarray based array-like - s5 = SphericalVoronoi(memoryview(self.points)) + s5 = SphericalVoronoi(memoryview(self.points)) # type: ignore[arg-type] assert_array_equal(s5.center, np.array([0, 0, 0])) assert_equal(s5.radius, 1) From 4f485c6a19506f91c5890bda6ca174772d73a589 Mon Sep 17 00:00:00 2001 From: czgdp1807 Date: Sat, 19 Jun 2021 13:29:20 +0530 Subject: [PATCH 25/48] fixed deprecated API calls --- scipy/optimize/minpack.h | 2 +- scipy/optimize/tnc/moduleTNC.c | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scipy/optimize/minpack.h b/scipy/optimize/minpack.h index 9872e4dc20b7..d53e5f49e240 100644 --- a/scipy/optimize/minpack.h +++ b/scipy/optimize/minpack.h @@ -146,7 +146,7 @@ static PyObject *call_python_function(PyObject *func, npy_intp n, double *x, PyO /* Call function object --- variable passed to routine. Extra arguments are in another passed variable. */ - if ((result = PyEval_CallObject(func, arglist))==NULL) { + if ((result = PyObject_CallObject(func, arglist))==NULL) { goto fail; } diff --git a/scipy/optimize/tnc/moduleTNC.c b/scipy/optimize/tnc/moduleTNC.c index 2c03a856a1bc..e8bee90d761b 100644 --- a/scipy/optimize/tnc/moduleTNC.c +++ b/scipy/optimize/tnc/moduleTNC.c @@ -67,7 +67,7 @@ static int function(double x[], double *f, double g[], void *state) memcpy(PyArray_DATA(py_x), x, (py_state->n)*sizeof(double)); arglist = Py_BuildValue("(N)", py_x); - result = PyEval_CallObject(py_state->py_function, arglist); + result = PyObject_CallObject(py_state->py_function, arglist); Py_DECREF(arglist); if (result == NULL) @@ -128,7 +128,7 @@ static void callback(double x[], void *state) memcpy(PyArray_DATA(py_x), x, (py_state->n)*sizeof(double)); arglist = Py_BuildValue("(N)", py_x); - result = PyEval_CallObject(py_state->py_callback, arglist); + result = PyObject_CallObject(py_state->py_callback, arglist); Py_DECREF(arglist); Py_XDECREF(result); } From a362cbf99cbeaa7c7603eedbcf0f2d7d1114a43e Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Sat, 19 Jun 2021 11:19:21 +0200 Subject: [PATCH 26/48] DOC: fix stats.pearsonr example that was failing in CI --- scipy/stats/mstats_basic.py | 2 +- scipy/stats/stats.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scipy/stats/mstats_basic.py b/scipy/stats/mstats_basic.py index 508bbddf56a8..8325bd920bc8 100644 --- a/scipy/stats/mstats_basic.py +++ b/scipy/stats/mstats_basic.py @@ -473,6 +473,7 @@ def pearsonr(x, y): e follow a normal distribution with mean zero and standard deviation s>0. >>> s = 0.5 + >>> x = stats.norm.rvs(size=500) >>> e = stats.norm.rvs(scale=s, size=500) >>> y = x + e >>> mstats.pearsonr(x, y) @@ -495,7 +496,6 @@ def pearsonr(x, y): cov(x, y) = E[x*y]. By definition, this equals E[x*abs(x)] which is zero by symmetry. The following lines of code illustrate this observation: - >>> x = stats.norm.rvs(size=500) >>> y = np.abs(x) >>> mstats.pearsonr(x, y) (-0.016172891856853524, 0.7182823678751942) # may vary diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py index b36ab053f30a..a4142b064c61 100644 --- a/scipy/stats/stats.py +++ b/scipy/stats/stats.py @@ -3997,6 +3997,7 @@ def pearsonr(x, y): e follow a normal distribution with mean zero and standard deviation s>0. >>> s = 0.5 + >>> x = stats.norm.rvs(size=500) >>> e = stats.norm.rvs(scale=s, size=500) >>> y = x + e >>> stats.pearsonr(x, y) @@ -4019,7 +4020,6 @@ def pearsonr(x, y): cov(x, y) = E[x*y]. By definition, this equals E[x*abs(x)] which is zero by symmetry. The following lines of code illustrate this observation: - >>> x = stats.norm.rvs(size=500) >>> y = np.abs(x) >>> stats.pearsonr(x, y) (-0.016172891856853524, 0.7182823678751942) # may vary From f860573eae038b0ad80cc27611a702136e55043d Mon Sep 17 00:00:00 2001 From: Tirth Patel Date: Sat, 19 Jun 2021 15:28:38 +0530 Subject: [PATCH 27/48] BUG: stats: fix discrete `.isf` to work at boundaries when `loc != 0` (#14242) --- scipy/stats/_distn_infrastructure.py | 7 +++++-- scipy/stats/tests/test_discrete_basic.py | 25 ++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/scipy/stats/_distn_infrastructure.py b/scipy/stats/_distn_infrastructure.py index ee37edb86b91..ae6037ba866b 100644 --- a/scipy/stats/_distn_infrastructure.py +++ b/scipy/stats/_distn_infrastructure.py @@ -3414,13 +3414,16 @@ def isf(self, q, *args, **kwds): cond0 = self._argcheck(*args) & (loc == loc) cond1 = (q > 0) & (q < 1) cond2 = (q == 1) & cond0 + cond3 = (q == 0) & cond0 cond = cond0 & cond1 # same problem as with ppf; copied from ppf and changed output = np.full(shape(cond), fill_value=self.badvalue, dtype='d') # output type 'd' to handle nin and inf - place(output, (q == 0)*(cond == cond), _b) - place(output, cond2, _a-1) + lower_bound = _a - 1 + loc + upper_bound = _b + loc + place(output, cond2*(cond == cond), lower_bound) + place(output, cond3*(cond == cond), upper_bound) # call place only if at least 1 valid argument if np.any(cond): diff --git a/scipy/stats/tests/test_discrete_basic.py b/scipy/stats/tests/test_discrete_basic.py index ce2ed7635cd4..c18b8408b893 100644 --- a/scipy/stats/tests/test_discrete_basic.py +++ b/scipy/stats/tests/test_discrete_basic.py @@ -148,6 +148,31 @@ def test_ppf_with_loc(dist, args): ) +@pytest.mark.parametrize('dist, args', distdiscrete) +def test_isf_with_loc(dist, args): + try: + distfn = getattr(stats, dist) + except TypeError: + distfn = dist + # check with a negative, no and positive relocation. + np.random.seed(1942349) + re_locs = [np.random.randint(-10, -1), 0, np.random.randint(1, 10)] + _a, _b = distfn.support(*args) + for loc in re_locs: + expected = _b + loc, _a - 1 + loc + res = distfn.isf(0., *args, loc=loc), distfn.isf(1., *args, loc=loc) + npt.assert_array_equal(expected, res) + # test broadcasting behaviour + re_locs = [np.random.randint(-10, -1, size=(5, 3)), + np.zeros((5, 3)), + np.random.randint(1, 10, size=(5, 3))] + _a, _b = distfn.support(*args) + for loc in re_locs: + expected = _b + loc, _a - 1 + loc + res = distfn.isf(0., *args, loc=loc), distfn.isf(1., *args, loc=loc) + npt.assert_array_equal(expected, res) + + def check_cdf_ppf(distfn, arg, supp, msg): # cdf is a step function, and ppf(q) = min{k : cdf(k) >= q, k integer} npt.assert_array_equal(distfn.ppf(distfn.cdf(supp, *arg), *arg), From 740c2d32a98e9c636c632f50c1e78a876c2c1471 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Sat, 19 Jun 2021 13:00:41 +0200 Subject: [PATCH 28/48] BLD: optimize: fix some warnings in moduleTNC and minpack.h The warnings were: ``` [2/15] Compiling C object scipy/optimize/moduleTNC.cpython-39-x86_64-linux-gnu.so.p/tnc_moduleTNC.c.o ../scipy/optimize/tnc/moduleTNC.c:26:19: warning: 'rcsid' defined but not used [-Wunused-const-variable=] 26 | static char const rcsid[] = | ^~~~~ [6/15] Compiling C object scipy/optimize/_minpack.cpython-39-x86_64-linux-gnu.so.p/_minpackmodule.c.o In file included from ../scipy/optimize/_minpackmodule.c:5: ../scipy/optimize/__minpack.h: In function 'minpack_hybrd': ../scipy/optimize/minpack.h:38:89: warning: unused variable 'jac_callback_info' [-Wunused-variable] 38 | #define STORE_VARS() ccallback_t callback; int callback_inited = 0; jac_callback_info_t jac_callback_info; | ^~~~~~~~~~~~~~~~~ ../scipy/optimize/__minpack.h:258:3: note: in expansion of macro 'STORE_VARS' 258 | STORE_VARS(); /* Define storage variables for global variables. */ | ^~~~~~~~~~ ../scipy/optimize/__minpack.h: In function 'minpack_lmdif': ../scipy/optimize/minpack.h:38:89: warning: unused variable 'jac_callback_info' [-Wunused-variable] 38 | #define STORE_VARS() ccallback_t callback; int callback_inited = 0; jac_callback_info_t jac_callback_info; | ^~~~~~~~~~~~~~~~~ ../scipy/optimize/__minpack.h:467:3: note: in expansion of macro 'STORE_VARS' 467 | STORE_VARS(); | ^~~~~~~~~~ ``` --- scipy/optimize/__minpack.h | 4 ++-- scipy/optimize/minpack.h | 1 + scipy/optimize/tnc/moduleTNC.c | 6 ++++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/scipy/optimize/__minpack.h b/scipy/optimize/__minpack.h index 627699245fef..8197ccc842da 100644 --- a/scipy/optimize/__minpack.h +++ b/scipy/optimize/__minpack.h @@ -255,7 +255,7 @@ static PyObject *minpack_hybrd(PyObject *dummy, PyObject *args) { int allocated = 0; double *wa = NULL; - STORE_VARS(); /* Define storage variables for global variables. */ + STORE_VARS_NO_INFO(); /* Define storage variables for global variables. */ if (!PyArg_ParseTuple(args, "OO|OidiiiddO", &fcn, &x0, &extra_args, &full_output, &xtol, &maxfev, &ml, &mu, &epsfcn, &factor, &o_diag)) return NULL; @@ -464,7 +464,7 @@ static PyObject *minpack_lmdif(PyObject *dummy, PyObject *args) { int allocated = 0; double *wa = NULL; - STORE_VARS(); + STORE_VARS_NO_INFO(); if (!PyArg_ParseTuple(args, "OO|OidddiddO", &fcn, &x0, &extra_args, &full_output, &ftol, &xtol, >ol, &maxfev, &epsfcn, &factor, &o_diag)) return NULL; diff --git a/scipy/optimize/minpack.h b/scipy/optimize/minpack.h index d53e5f49e240..f30035f87eff 100644 --- a/scipy/optimize/minpack.h +++ b/scipy/optimize/minpack.h @@ -36,6 +36,7 @@ the result tuple when the full_output argument is non-zero. #define PYERR2(errobj,message) {PyErr_Print(); PyErr_SetString(errobj, message); goto fail;} #define STORE_VARS() ccallback_t callback; int callback_inited = 0; jac_callback_info_t jac_callback_info; +#define STORE_VARS_NO_INFO() ccallback_t callback; int callback_inited = 0; #define INIT_FUNC(fun,arg,errobj) do { /* Get extra arguments or set to zero length tuple */ \ if (arg == NULL) { \ diff --git a/scipy/optimize/tnc/moduleTNC.c b/scipy/optimize/tnc/moduleTNC.c index e8bee90d761b..1e0e652be5e5 100644 --- a/scipy/optimize/tnc/moduleTNC.c +++ b/scipy/optimize/tnc/moduleTNC.c @@ -23,8 +23,10 @@ * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ -static char const rcsid[] = - "@(#) $Jeannot: moduleTNC.c,v 1.12 2005/01/28 18:27:31 js Exp $"; +/* + * static char const rcsid[] = + * "@(#) $Jeannot: moduleTNC.c,v 1.12 2005/01/28 18:27:31 js Exp $"; + */ #include "Python.h" #include "numpy/arrayobject.h" From b8d6a0a0d6a275555a58136b0bc35c0b6c872aa9 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Sat, 19 Jun 2021 14:15:52 +0200 Subject: [PATCH 29/48] BLD: fix include order and build warnings for `optimize/_trlib` `Python.h` must always be included before standard system headers, and numpy includes pull in `Python.h` so must come first too. This was the cause of build warnings like (repeated quite a few times): ``` [10/17] Compiling C object scipy/optimize/_trlib/_trlib.cpython-39-x86_64-linux-gnu.so.p/trlib_krylov.c.o In file included from /home/rgommers/anaconda3/envs/scipy-meson/include/python3.9/Python.h:85, from /home/rgommers/anaconda3/envs/scipy-meson/lib/python3.9/site-packages/numpy/core/include/numpy/ndarrayobject.h:11, from /home/rgommers/anaconda3/envs/scipy-meson/lib/python3.9/site-packages/numpy/core/include/numpy/arrayobject.h:4, from ../scipy/optimize/_trlib/trlib_private.h:38, from ../scipy/optimize/_trlib/trlib_krylov.c:26: /home/rgommers/anaconda3/envs/scipy-meson/include/python3.9/pytime.h:153:60: warning: 'struct timespec' declared inside parameter list will not be visible outside of this definition or declaration 153 | PyAPI_FUNC(int) _PyTime_FromTimespec(_PyTime_t *tp, struct timespec *ts); | ^~~~~~~~ /home/rgommers/anaconda3/envs/scipy-meson/include/python3.9/pytime.h:158:56: warning: 'struct timespec' declared inside parameter list will not be visible outside of this definition or declaration 158 | PyAPI_FUNC(int) _PyTime_AsTimespec(_PyTime_t t, struct timespec *ts); | ^~~~~~~~ ``` The unused variable warnings were: ``` ../scipy/optimize/_trlib/trlib_krylov.c: In function 'trlib_krylov_min': ../scipy/optimize/_trlib/trlib_krylov.c:725:25: warning: unused variable 'obj' [-Wunused-variable] 725 | trlib_flt_t obj = fwork[8]; | ^~~ ../scipy/optimize/_trlib/trlib_krylov.c:734:25: warning: unused variable 'lam' [-Wunused-variable] 734 | trlib_flt_t lam = fwork[7]; | ^~~ ``` --- scipy/optimize/_trlib/trlib_krylov.c | 2 -- scipy/optimize/_trlib/trlib_private.h | 6 +++--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/scipy/optimize/_trlib/trlib_krylov.c b/scipy/optimize/_trlib/trlib_krylov.c index d38ac442f933..95832e208690 100644 --- a/scipy/optimize/_trlib/trlib_krylov.c +++ b/scipy/optimize/_trlib/trlib_krylov.c @@ -722,7 +722,6 @@ trlib_int_t trlib_krylov_min( if( ret < 10 && *outerstatus < 100 && convexify ) { // exit, check if we should convexify trlib_flt_t lam = fwork[7]; - trlib_flt_t obj = fwork[8]; if( lam > 1e-2*fmax(1.0, fwork[13]) && fwork[14] < 0.0 && fabs(fwork[14]) < 1e-8 * fwork[13]) { // do only if it seems to make sense based on eigenvalue estimation // ask caller to compute objective value *outerstatus = 200 + ret; @@ -731,7 +730,6 @@ trlib_int_t trlib_krylov_min( } } if( *outerstatus > 200 && *outerstatus < 300 ) { - trlib_flt_t lam = fwork[7]; trlib_flt_t obj = fwork[8]; trlib_flt_t realobj = g_dot_g; if( fabs(obj - realobj) > fmax(1e-6, 1e-1*fabs(realobj)) || realobj > 0.0) { diff --git a/scipy/optimize/_trlib/trlib_private.h b/scipy/optimize/_trlib/trlib_private.h index c77a81dd3176..402db57fdfb3 100644 --- a/scipy/optimize/_trlib/trlib_private.h +++ b/scipy/optimize/_trlib/trlib_private.h @@ -28,6 +28,9 @@ /* #undef TRLIB_MEASURE_TIME */ /* #undef TRLIB_MEASURE_SUBTIME */ +#include "numpy/arrayobject.h" +#include "npy_cblas.h" + #include "trlib.h" #include #include @@ -35,9 +38,6 @@ #include #include -#include "numpy/arrayobject.h" -#include "npy_cblas.h" - // blas void BLAS_FUNC(daxpy)(CBLAS_INT *n, double *alpha, double *x, CBLAS_INT *incx, double *y, CBLAS_INT *incy); void BLAS_FUNC(dscal)(CBLAS_INT *n, double *alpha, double *x, CBLAS_INT *incx); From f6d7d396641665fa6ee4dac85727d7aca5aa3bd9 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Sat, 19 Jun 2021 12:32:43 +0200 Subject: [PATCH 30/48] CI: pin mypy to 0.902 and fix one CI failure Newer mypy releases often come with new failures (because Mypy got smarter in this case, so the `#type: ignore` is no longer needed), so we want to pin the mypy version in CI. --- mypy_requirements.txt | 2 +- scipy/stats/_qmc.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/mypy_requirements.txt b/mypy_requirements.txt index 570dd3bddba9..5a307d119dce 100644 --- a/mypy_requirements.txt +++ b/mypy_requirements.txt @@ -1,4 +1,4 @@ # Note: this should disappear at some point. For now, please keep it # in sync with the dev dependencies in pyproject.toml -mypy +mypy==0.902 typing_extensions diff --git a/scipy/stats/_qmc.py b/scipy/stats/_qmc.py index 61290c62dbb0..c81c2b58e299 100644 --- a/scipy/stats/_qmc.py +++ b/scipy/stats/_qmc.py @@ -884,9 +884,9 @@ def random(self, n: IntNumber = 1) -> np.ndarray: """ if self.centered: - samples = 0.5 + samples: np.ndarray | float = 0.5 else: - samples = self.rng.uniform(size=(n, self.d)) # type: ignore[assignment] + samples = self.rng.uniform(size=(n, self.d)) perms = np.tile(np.arange(1, n + 1), (self.d, 1)) for i in range(self.d): From d0d7eeadc592e5dac787d582fb9b08750e043825 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Sun, 20 Jun 2021 13:51:27 -0600 Subject: [PATCH 31/48] DOC: forward port 1.7.0 relnotes * forward port the SciPy `1.7.0` release notes following the final release --- doc/release/1.7.0-notes.rst | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/doc/release/1.7.0-notes.rst b/doc/release/1.7.0-notes.rst index f527ba3c311d..6fef55c2a5d5 100644 --- a/doc/release/1.7.0-notes.rst +++ b/doc/release/1.7.0-notes.rst @@ -28,7 +28,7 @@ Highlights of this release - A new submodule for quasi-Monte Carlo, `scipy.stats.qmc`, was added - The documentation design was updated to use the same PyData-Sphinx theme as - other NumFOCUS packages like NumPy. + NumPy and other ecosystem libraries. - We now vendor and leverage the Boost C++ library to enable numerous improvements for long-standing weaknesses in `scipy.stats` - `scipy.stats` has six new distributions, eight new (or overhauled) @@ -612,6 +612,7 @@ Issues closed for 1.7.0 * `#13793 `__: CI: CircleCI doc build failure * `#13840 `__: manylinux1 builds are failing because of C99 usage in \`special/_cosine.c\` * `#13850 `__: CI: Homebrew is failing due to bintray +* `#13875 `__: BUG: chi2_contingency with Yates correction * `#13878 `__: BUG: \`signal.get_window\` argument handling issue * `#13880 `__: Remove all usages of numpy.compat * `#13896 `__: Boschloo's Test for More Powerful Hypothesis Testing of 2x2 Contingency... @@ -628,6 +629,10 @@ Issues closed for 1.7.0 * `#14048 `__: DOC: missing git submodule information * `#14055 `__: linalg.solve: Unclear error when using assume_a='her' with real... * `#14093 `__: DOC: Inconsistency in the definition of default values in the... +* `#14158 `__: TST, BUG: test_rbfinterp.py -- test_interpolation_misfit_1d fails... +* `#14170 `__: TST: signal submodule test_filtfilt_gust failing on 32-bit amd64... +* `#14194 `__: MAINT: download-wheels.py missing import +* `#14199 `__: Generated sources for biasedurn extension are broken in 1.7.0rc1 *********************** @@ -931,6 +936,7 @@ Pull requests for 1.7.0 * `#13910 `__: DOC: update Readme * `#13911 `__: MAINT: use dict built-in rather than OrderedDict * `#13920 `__: BUG: Reactivate conda environment in init +* `#13925 `__: BUG: stats: magnitude of Yates' correction <= abs(observed-expected)... * `#13926 `__: DOC: correct return type in disjoint_set.subsets docstring * `#13927 `__: DOC/MAINT: Add copyright notice to qmc.primes_from_2_to * `#13928 `__: BUG: DOC: signal: fix need argument config and add missing doc... @@ -1011,9 +1017,19 @@ Pull requests for 1.7.0 * `#14110 `__: DOC: mailmap update * `#14113 `__: ENH: stats: bootstrap: add \`paired\` parameter * `#14116 `__: MAINT: fix deprecated Python C API usage in odr +* `#14118 `__: DOC: 1.7.0 release notes * `#14125 `__: DOC: fix typo * `#14126 `__: ENH: stats: bootstrap: add \`batch\` parameter to control batch... * `#14127 `__: CI: upgrade pip in benchmarks CI run * `#14130 `__: BUG: Fix trust-constr report TypeError if verbose is set to 2... * `#14133 `__: MAINT: interpolate: raise NotImplementedError not ValueError * `#14139 `__: FIX/DOC: lsqr doctests print failure +* `#14145 `__: MAINT: 1.7.x version pins ("backport") +* `#14146 `__: MAINT: commit count if no tag +* `#14164 `__: TST, BUG: fix rbf matrix value +* `#14166 `__: CI, MAINT: restrictions on pre-release CI +* `#14171 `__: TST: signal: Bump tolerances for a test of Gustafsson's... +* `#14175 `__: TST: stats: Loosen tolerance in some binomtest tests. +* `#14182 `__: MAINT: stats: Update ppcc_plot and ppcc_max docstring. +* `#14195 `__: MAINT: download-wheels missing import +* `#14230 `__: REL: stop shipping generated Cython sources in sdist From edb50b09e92cc1b493242076b3f63e89397032e4 Mon Sep 17 00:00:00 2001 From: Gagandeep Singh Date: Tue, 22 Jun 2021 12:22:58 +0530 Subject: [PATCH 32/48] MAINT: Replaced direct field access in `PyArrayObject*` with wrapper functions (#14268) Fixed code to access data and dimensions from PyArrayObject --- scipy/integrate/__quadpack.h | 96 ++++++++++++++++++------------------ 1 file changed, 48 insertions(+), 48 deletions(-) diff --git a/scipy/integrate/__quadpack.h b/scipy/integrate/__quadpack.h index b2d88ac2ffab..9d92a21da115 100644 --- a/scipy/integrate/__quadpack.h +++ b/scipy/integrate/__quadpack.h @@ -400,11 +400,11 @@ static PyObject *quadpack_qagse(PyObject *dummy, PyObject *args) { ap_rlist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); ap_elist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); if (setjmp(callback.error_buf) != 0) { goto fail; @@ -483,11 +483,11 @@ static PyObject *quadpack_qagie(PyObject *dummy, PyObject *args) { ap_elist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); if (setjmp(callback.error_buf) != 0) { goto fail; @@ -562,9 +562,9 @@ static PyObject *quadpack_qagpe(PyObject *dummy, PyObject *args) { ap_points = (PyArrayObject *)PyArray_ContiguousFromObject(o_points, NPY_DOUBLE, 1, 1); if (ap_points == NULL) goto fail; - npts2 = ap_points->dimensions[0]; + npts2 = PyArray_DIMS(ap_points)[0]; npts2_shape[0] = npts2; - points = (double *)ap_points->data; + points = (double *)PyArray_DATA(ap_points); /* Set up iwork and work arrays */ ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,F_INT_NPY); @@ -576,14 +576,14 @@ static PyObject *quadpack_qagpe(PyObject *dummy, PyObject *args) { ap_level = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,F_INT_NPY); ap_ndin = (PyArrayObject *)PyArray_SimpleNew(1,npts2_shape,F_INT_NPY); if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL || ap_pts == NULL || ap_level == NULL || ap_ndin == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; - pts = (double *)ap_pts->data; - level = (F_INT *)ap_level->data; - ndin = (F_INT *)ap_ndin->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); + pts = (double *)PyArray_DATA(ap_pts); + level = (F_INT *)PyArray_DATA(ap_level); + ndin = (F_INT *)PyArray_DATA(ap_ndin); if (setjmp(callback.error_buf) != 0) { goto fail; @@ -671,7 +671,7 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) { if (o_chebmo != NULL) { ap_chebmo = (PyArrayObject *)PyArray_ContiguousFromObject(o_chebmo, NPY_DOUBLE, 2, 2); if (ap_chebmo == NULL) goto fail; - if (ap_chebmo->dimensions[1] != maxp1 || ap_chebmo->dimensions[0] != 25) + if (PyArray_DIMS(ap_chebmo)[1] != maxp1 || PyArray_DIMS(ap_chebmo)[0] != 25) PYERR(quadpack_error,"Chebyshev moment array has the wrong size."); } else { @@ -680,7 +680,7 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) { ap_chebmo = (PyArrayObject *)PyArray_SimpleNew(2,sz,NPY_DOUBLE); if (ap_chebmo == NULL) goto fail; } - chebmo = (double *) ap_chebmo->data; + chebmo = (double *) PyArray_DATA(ap_chebmo); /* Set up iwork and work arrays */ ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,F_INT_NPY); @@ -690,12 +690,12 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) { ap_rlist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); ap_elist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); if (ap_iord == NULL || ap_nnlog == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - nnlog = (F_INT *)ap_nnlog->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + nnlog = (F_INT *)PyArray_DATA(ap_nnlog); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); if (setjmp(callback.error_buf) != 0) { goto fail; @@ -781,7 +781,7 @@ static PyObject *quadpack_qawfe(PyObject *dummy, PyObject *args) { sz[1] = maxp1; ap_chebmo = (PyArrayObject *)PyArray_SimpleNew(2,sz,NPY_DOUBLE); if (ap_chebmo == NULL) goto fail; - chebmo = (double *) ap_chebmo->data; + chebmo = (double *) PyArray_DATA(ap_chebmo); /* Set up iwork and work arrays */ ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,F_INT_NPY); @@ -794,15 +794,15 @@ static PyObject *quadpack_qawfe(PyObject *dummy, PyObject *args) { ap_erlst = (PyArrayObject *)PyArray_SimpleNew(1,limlst_shape,NPY_DOUBLE); ap_ierlst = (PyArrayObject *)PyArray_SimpleNew(1,limlst_shape,F_INT_NPY); if (ap_iord == NULL || ap_nnlog == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL || ap_rslst == NULL || ap_erlst == NULL || ap_ierlst == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - nnlog = (F_INT *)ap_nnlog->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; - rslst = (double *)ap_rslst->data; - erlst = (double *)ap_erlst->data; - ierlst = (F_INT *)ap_ierlst->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + nnlog = (F_INT *)PyArray_DATA(ap_nnlog); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); + rslst = (double *)PyArray_DATA(ap_rslst); + erlst = (double *)PyArray_DATA(ap_erlst); + ierlst = (F_INT *)PyArray_DATA(ap_ierlst); if (setjmp(callback.error_buf) != 0) { goto fail; @@ -889,11 +889,11 @@ static PyObject *quadpack_qawce(PyObject *dummy, PyObject *args) { ap_rlist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); ap_elist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); if (setjmp(callback.error_buf) != 0) { goto fail; @@ -972,11 +972,11 @@ static PyObject *quadpack_qawse(PyObject *dummy, PyObject *args) { ap_rlist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); ap_elist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE); if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL) goto fail; - iord = (F_INT *)ap_iord->data; - alist = (double *)ap_alist->data; - blist = (double *)ap_blist->data; - rlist = (double *)ap_rlist->data; - elist = (double *)ap_elist->data; + iord = (F_INT *)PyArray_DATA(ap_iord); + alist = (double *)PyArray_DATA(ap_alist); + blist = (double *)PyArray_DATA(ap_blist); + rlist = (double *)PyArray_DATA(ap_rlist); + elist = (double *)PyArray_DATA(ap_elist); if (setjmp(callback.error_buf) != 0) { goto fail; From bf73d92e9ff02b25543d45c2bf30bd04fe3c47fd Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 13:50:18 +0530 Subject: [PATCH 33/48] MAINT: Update uarray vendor lib --- scipy/_lib/_uarray/__init__.py | 12 +- scipy/_lib/_uarray/_backend.py | 322 ++++++- scipy/_lib/_uarray/_uarray_dispatch.cxx | 1009 ++++++++++++++++++---- scipy/_lib/_uarray/small_dynamic_array.h | 228 +++++ scipy/_lib/_uarray/vectorcall.cxx | 235 +++++ scipy/_lib/_uarray/vectorcall.h | 48 + 6 files changed, 1670 insertions(+), 184 deletions(-) create mode 100644 scipy/_lib/_uarray/small_dynamic_array.h create mode 100644 scipy/_lib/_uarray/vectorcall.cxx create mode 100644 scipy/_lib/_uarray/vectorcall.h diff --git a/scipy/_lib/_uarray/__init__.py b/scipy/_lib/_uarray/__init__.py index 9bcaf4029907..b3c395d0aa62 100644 --- a/scipy/_lib/_uarray/__init__.py +++ b/scipy/_lib/_uarray/__init__.py @@ -1,10 +1,10 @@ """ -.. note:: +.. note: If you are looking for overrides for NumPy-specific methods, see the documentation for :obj:`unumpy`. This page explains how to write back-ends and multimethods. -``uarray`` is built around a back-end protocol and overridable multimethods. +``uarray`` is built around a back-end protocol, and overridable multimethods. It is necessary to define multimethods for back-ends to be able to override them. See the documentation of :obj:`generate_multimethod` on how to write multimethods. @@ -97,7 +97,7 @@ ('override_me', (1, '2'), {}) You also have the option to return ``NotImplemented``, in which case processing moves on -to the next back-end, which, in this case, doesn't exist. The same applies to +to the next back-end, which in this case, doesn't exist. The same applies to ``__ua_convert__``. >>> be.__ua_function__ = lambda *a, **kw: NotImplemented @@ -105,7 +105,7 @@ ... overridden_me(1, "2") Traceback (most recent call last): ... -uarray.backend.BackendNotImplementedError: ... +uarray.BackendNotImplementedError: ... The last possibility is if we don't have ``__ua_convert__``, in which case the job is left up to ``__ua_function__``, but putting things back into arrays after conversion will not be @@ -113,5 +113,7 @@ """ from ._backend import * +from ._version import get_versions -__version__ = '0.5.1+49.g4c3f1d7.scipy' +__version__ = get_versions()["version"] +del get_versions diff --git a/scipy/_lib/_uarray/_backend.py b/scipy/_lib/_uarray/_backend.py index 930f1c243c7e..3a38d74c82a7 100644 --- a/scipy/_lib/_uarray/_backend.py +++ b/scipy/_lib/_uarray/_backend.py @@ -1,10 +1,11 @@ import typing +import types import inspect import functools -from . import _uarray # type: ignore -import copyreg # type: ignore -import atexit +from . import _uarray +import copyreg import pickle +import contextlib ArgumentExtractorType = typing.Callable[..., typing.Tuple["Dispatchable", ...]] ArgumentReplacerType = typing.Callable[ @@ -16,6 +17,7 @@ _Function, _SkipBackendContext, _SetBackendContext, + _BackendState, ) __all__ = [ @@ -23,6 +25,8 @@ "set_global_backend", "skip_backend", "register_backend", + "determine_backend", + "determine_backend_multi", "clear_backends", "create_multimethod", "generate_multimethod", @@ -30,17 +34,31 @@ "BackendNotImplementedError", "Dispatchable", "wrap_single_convertor", + "wrap_single_convertor_instance", "all_of_type", "mark_as", + "set_state", + "get_state", + "reset_state", + "_BackendState", + "_SkipBackendContext", + "_SetBackendContext", ] -def unpickle_function(mod_name, qname): +def unpickle_function(mod_name, qname, self_): import importlib try: module = importlib.import_module(mod_name) - func = getattr(module, qname) + qname = qname.split(".") + func = module + for q in qname: + func = getattr(func, q) + + if self_ is not None: + func = types.MethodType(func, self_) + return func except (ImportError, AttributeError) as e: from pickle import UnpicklingError @@ -51,9 +69,10 @@ def unpickle_function(mod_name, qname): def pickle_function(func): mod_name = getattr(func, "__module__", None) qname = getattr(func, "__qualname__", None) + self_ = getattr(func, "__self__", None) try: - test = unpickle_function(mod_name, qname) + test = unpickle_function(mod_name, qname, self_) except pickle.UnpicklingError: test = None @@ -62,11 +81,73 @@ def pickle_function(func): "Can't pickle {}: it's not the same object as {}".format(func, test) ) - return unpickle_function, (mod_name, qname) + return unpickle_function, (mod_name, qname, self_) + + +def pickle_state(state): + return _uarray._BackendState._unpickle, state._pickle() + + +def pickle_set_backend_context(ctx): + return _SetBackendContext, ctx._pickle() + + +def pickle_skip_backend_context(ctx): + return _SkipBackendContext, ctx._pickle() copyreg.pickle(_Function, pickle_function) -atexit.register(_uarray.clear_all_globals) +copyreg.pickle(_uarray._BackendState, pickle_state) +copyreg.pickle(_SetBackendContext, pickle_set_backend_context) +copyreg.pickle(_SkipBackendContext, pickle_skip_backend_context) + + +def get_state(): + """ + Returns an opaque object containing the current state of all the backends. + + Can be used for synchronization between threads/processes. + + See Also + -------- + set_state + Sets the state returned by this function. + """ + return _uarray.get_state() + + +@contextlib.contextmanager +def reset_state(): + """ + Returns a context manager that resets all state once exited. + + See Also + -------- + set_state + Context manager that sets the backend state. + get_state + Gets a state to be set by this context manager. + """ + with set_state(get_state()): + yield + + +@contextlib.contextmanager +def set_state(state): + """ + A context manager that sets the state of the backends to one returned by :obj:`get_state`. + + See Also + -------- + get_state + Gets a state to be set by this context manager. + """ + old_state = get_state() + _uarray.set_state(state) + try: + yield + finally: + _uarray.set_state(old_state, True) def create_multimethod(*args, **kwargs): @@ -80,7 +161,8 @@ def create_multimethod(*args, **kwargs): See Also -------- - generate_multimethod : Generates a multimethod. + generate_multimethod + Generates a multimethod. """ def wrapper(a): @@ -93,7 +175,7 @@ def generate_multimethod( argument_extractor: ArgumentExtractorType, argument_replacer: ArgumentReplacerType, domain: str, - default: typing.Optional[typing.Callable] = None + default: typing.Optional[typing.Callable] = None, ): """ Generates a multimethod. @@ -109,7 +191,7 @@ def generate_multimethod( return an (args, kwargs) pair with the dispatchables replaced inside the args/kwargs. domain : str A string value indicating the domain of this multimethod. - default : Optional[Callable], optional + default: Optional[Callable], optional The default implementation of this multimethod, where ``None`` (the default) specifies there is no default implementation. @@ -138,7 +220,8 @@ def generate_multimethod( >>> overridden_me(1, "a") Traceback (most recent call last): ... - uarray.backend.BackendNotImplementedError: ... + uarray.BackendNotImplementedError: ... + >>> overridden_me2 = generate_multimethod( ... override_me, override_replacer, "ua_examples", default=lambda x, y: (x, y) ... ) @@ -147,7 +230,7 @@ def generate_multimethod( See Also -------- - uarray : + uarray See the module documentation for how to override the method by creating backends. """ kw_defaults, arg_defaults, opts = get_defaults(argument_extractor) @@ -178,8 +261,8 @@ def set_backend(backend, coerce=False, only=False): See Also -------- - skip_backend : A context manager that allows skipping of backends. - set_global_backend : Set a single, global backend for a domain. + skip_backend: A context manager that allows skipping of backends. + set_global_backend: Set a single, global backend for a domain. """ try: return backend.__ua_cache__["set", coerce, only] @@ -206,8 +289,8 @@ def skip_backend(backend): See Also -------- - set_backend : A context manager that allows setting of backends. - set_global_backend : Set a single, global backend for a domain. + set_backend: A context manager that allows setting of backends. + set_global_backend: Set a single, global backend for a domain. """ try: return backend.__ua_cache__["skip"] @@ -239,7 +322,7 @@ def get_defaults(f): return kw_defaults, tuple(arg_defaults), opts -def set_global_backend(backend, coerce=False, only=False): +def set_global_backend(backend, coerce=False, only=False, *, try_last=False): """ This utility method replaces the default backend for permanent use. It will be tried in the list of backends automatically, unless the @@ -258,13 +341,20 @@ def set_global_backend(backend, coerce=False, only=False): ---------- backend The backend to register. + coerce : bool + Whether to coerce input types when trying this backend. + only : bool + If ``True``, no more backends will be tried if this fails. + Implied by ``coerce=True``. + try_last : bool + If ``True``, the global backend is tried after registered backends. See Also -------- - set_backend : A context manager that allows setting of backends. - skip_backend : A context manager that allows skipping of backends. + set_backend: A context manager that allows setting of backends. + skip_backend: A context manager that allows skipping of backends. """ - _uarray.set_global_backend(backend, coerce, only) + _uarray.set_global_backend(backend, coerce, only, try_last) def register_backend(backend): @@ -290,7 +380,7 @@ def clear_backends(domain, registered=True, globals=False): .. warning:: We caution library authors against using this function in their code. We do *not* support this use-case. This function - is meant to be used only by the users themselves. + is meant to be used only by users themselves. .. warning:: Do NOT use this method inside a multimethod call, or the @@ -423,3 +513,191 @@ def __ua_convert__(dispatchables, coerce): return converted return __ua_convert__ + + +def wrap_single_convertor_instance(convert_single): + """ + Wraps a ``__ua_convert__`` defined for a single element to all elements. + If any of them return ``NotImplemented``, the operation is assumed to be + undefined. + + Accepts a signature of (value, type, coerce). + """ + + @functools.wraps(convert_single) + def __ua_convert__(self, dispatchables, coerce): + converted = [] + for d in dispatchables: + c = convert_single(self, d.value, d.type, coerce and d.coercible) + + if c is NotImplemented: + return NotImplemented + + converted.append(c) + + return converted + + return __ua_convert__ + + +def determine_backend(value, dispatch_type, *, domain, only=True, coerce=False): + """Set the backend to the first active backend that supports ``value`` + + This is useful for functions that call multimethods without any dispatchable + arguments. You can use :func:`determine_backend` to ensure the same backend + is used everywhere in a block of multimethod calls. + + Parameters + ---------- + value + The value being tested + dispatch_type + The dispatch type associated with ``value``, aka + ":ref:`marking `". + domain: string + The domain to query for backends and set. + coerce: bool + Whether or not to allow coercion to the backend's types. Implies ``only``. + only: bool + Whether or not this should be the last backend to try. + + See Also + -------- + set_backend: For when you know which backend to set + + Notes + ----- + + Support is determined by the ``__ua_convert__`` protocol. Backends not + supporting the type must return ``NotImplemented`` from their + ``__ua_convert__`` if they don't support input of that type. + + Examples + -------- + + Suppose we have two backends ``BackendA`` and ``BackendB`` each supporting + different types, ``TypeA`` and ``TypeB``. Neither supporting the other type: + + >>> with ua.set_backend(ex.BackendA): + ... ex.call_multimethod(ex.TypeB(), ex.TypeB()) + Traceback (most recent call last): + ... + uarray.BackendNotImplementedError: ... + + Now consider a multimethod that creates a new object of ``TypeA``, or + ``TypeB`` depending on the active backend. + + >>> with ua.set_backend(ex.BackendA), ua.set_backend(ex.BackendB): + ... res = ex.creation_multimethod() + ... ex.call_multimethod(res, ex.TypeA()) + Traceback (most recent call last): + ... + uarray.BackendNotImplementedError: ... + + ``res`` is an object of ``TypeB`` because ``BackendB`` is set in the + innermost with statement. So, ``call_multimethod`` fails since the types + don't match. + + Instead, we need to first find a backend suitable for all of our objects. + + >>> with ua.set_backend(ex.BackendA), ua.set_backend(ex.BackendB): + ... x = ex.TypeA() + ... with ua.determine_backend(x, "mark", domain="ua_examples"): + ... res = ex.creation_multimethod() + ... ex.call_multimethod(res, x) + TypeA + + """ + dispatchables = (Dispatchable(value, dispatch_type, coerce),) + backend = _uarray.determine_backend(domain, dispatchables, coerce) + + return set_backend(backend, coerce=coerce, only=only) + + +def determine_backend_multi( + dispatchables, *, domain, only=True, coerce=False, **kwargs +): + """Set a backend supporting all ``dispatchables`` + + This is useful for functions that call multimethods without any dispatchable + arguments. You can use :func:`determine_backend_multi` to ensure the same + backend is used everywhere in a block of multimethod calls involving + multiple arrays. + + Parameters + ---------- + dispatchables: Sequence[Union[uarray.Dispatchable, Any]] + The dispatchables that must be supported + domain: string + The domain to query for backends and set. + coerce: bool + Whether or not to allow coercion to the backend's types. Implies ``only``. + only: bool + Whether or not this should be the last backend to try. + dispatch_type: Optional[Any] + The default dispatch type associated with ``dispatchables``, aka + ":ref:`marking `". + + See Also + -------- + determine_backend: For a single dispatch value + set_backend: For when you know which backend to set + + Notes + ----- + + Support is determined by the ``__ua_convert__`` protocol. Backends not + supporting the type must return ``NotImplemented`` from their + ``__ua_convert__`` if they don't support input of that type. + + Examples + -------- + + :func:`determine_backend` allows the backend to be set from a single + object. :func:`determine_backend_multi` allows multiple objects to be + checked simultaneously for support in the backend. Suppose we have a + ``BackendAB`` which supports ``TypeA`` and ``TypeB`` in the same call, + and a ``BackendBC`` that doesn't support ``TypeA``. + + >>> with ua.set_backend(ex.BackendAB), ua.set_backend(ex.BackendBC): + ... a, b = ex.TypeA(), ex.TypeB() + ... with ua.determine_backend_multi( + ... [ua.Dispatchable(a, "mark"), ua.Dispatchable(b, "mark")], + ... domain="ua_examples" + ... ): + ... res = ex.creation_multimethod() + ... ex.call_multimethod(res, a, b) + TypeA + + This won't call ``BackendBC`` because it doesn't support ``TypeA``. + + We can also use leave out the ``ua.Dispatchable`` if we specify the + default ``dispatch_type`` for the ``dispatchables`` argument. + + >>> with ua.set_backend(ex.BackendAB), ua.set_backend(ex.BackendBC): + ... a, b = ex.TypeA(), ex.TypeB() + ... with ua.determine_backend_multi( + ... [a, b], dispatch_type="mark", domain="ua_examples" + ... ): + ... res = ex.creation_multimethod() + ... ex.call_multimethod(res, a, b) + TypeA + + """ + if "dispatch_type" in kwargs: + disp_type = kwargs.pop("dispatch_type") + dispatchables = tuple( + d if isinstance(d, Dispatchable) else Dispatchable(d, disp_type) + for d in dispatchables + ) + else: + dispatchables = tuple(dispatchables) + if not all(isinstance(d, Dispatchable) for d in dispatchables): + raise TypeError("dispatchables must be instances of uarray.Dispatchable") + + if len(kwargs) != 0: + raise TypeError("Received unexpected keyword arguments: {}".format(kwargs)) + + backend = _uarray.determine_backend(domain, dispatchables, coerce) + + return set_backend(backend, coerce=coerce, only=only) diff --git a/scipy/_lib/_uarray/_uarray_dispatch.cxx b/scipy/_lib/_uarray/_uarray_dispatch.cxx index 0ea4272b3689..da02a06de2b6 100644 --- a/scipy/_lib/_uarray/_uarray_dispatch.cxx +++ b/scipy/_lib/_uarray/_uarray_dispatch.cxx @@ -1,14 +1,17 @@ +#include "small_dynamic_array.h" +#include "vectorcall.h" + +#include #include #include #include +#include #include #include #include #include -#include - namespace { @@ -88,9 +91,17 @@ py_ref py_make_tuple(const Ts &... args) { return py_ref::steal(PyTuple_Pack(sizeof...(args), py_get(args)...)); } +py_ref py_bool(bool input) { return py_ref::ref(input ? Py_True : Py_False); } + +template +constexpr size_t array_size(const T (&array)[N]) { + return N; +} + struct backend_options { py_ref backend; - bool coerce, only; + bool coerce = false; + bool only = false; bool operator==(const backend_options & other) const { return ( @@ -106,6 +117,7 @@ struct backend_options { struct global_backends { backend_options global; std::vector registered; + bool try_global_backend_last = false; }; struct local_backends { @@ -113,15 +125,19 @@ struct local_backends { std::vector preferred; }; +using global_state_t = std::unordered_map; +using local_state_t = std::unordered_map; static py_ref BackendNotImplementedError; -static std::unordered_map global_domain_map; -thread_local std::unordered_map local_domain_map; +static global_state_t global_domain_map; +thread_local global_state_t * current_global_state = &global_domain_map; +thread_local global_state_t thread_local_domain_map; +thread_local local_state_t local_domain_map; /** Constant Python string identifiers Using these with PyObject_GetAttr is faster than PyObject_GetAttrString which -has to create a new Python string internally. +has to create a new python string internally. */ struct { py_ref ua_convert; @@ -151,9 +167,23 @@ struct { } } identifiers; -std::string domain_to_string(PyObject * domain) { +bool domain_validate(PyObject * domain) { if (!PyUnicode_Check(domain)) { PyErr_SetString(PyExc_TypeError, "__ua_domain__ must be a string"); + return false; + } + + auto size = PyUnicode_GetLength(domain); + if (size == 0) { + PyErr_SetString(PyExc_ValueError, "__ua_domain__ must be non-empty"); + return false; + } + + return true; +} + +std::string domain_to_string(PyObject * domain) { + if (!domain_validate(domain)) { return {}; } @@ -170,44 +200,399 @@ std::string domain_to_string(PyObject * domain) { return std::string(str, size); } -std::string backend_to_domain_string(PyObject * backend) { +Py_ssize_t backend_get_num_domains(PyObject * backend) { auto domain = py_ref::steal(PyObject_GetAttr(backend, identifiers.ua_domain.get())); if (!domain) - return {}; + return -1; + + if (PyUnicode_Check(domain.get())) { + return 1; + } + + if (!PySequence_Check(domain.get())) { + PyErr_SetString( + PyExc_TypeError, + "__ua_domain__ must be a string or sequence of strings"); + return -1; + } - return domain_to_string(domain.get()); + return PySequence_Size(domain.get()); } +enum class LoopReturn { Continue, Break, Error }; -/** Use to clean up Python references before the interpreter is finalized. - * - * This must be installed in a Python atexit handler. This prevents Py_DECREF - * being called after the interpreter has already shut down. - */ -PyObject * clear_all_globals(PyObject * /* self */, PyObject * /* args */) { +template +LoopReturn backend_for_each_domain(PyObject * backend, Func f) { + auto domain = + py_ref::steal(PyObject_GetAttr(backend, identifiers.ua_domain.get())); + if (!domain) + return LoopReturn::Error; + + if (PyUnicode_Check(domain.get())) { + return f(domain.get()); + } + + if (!PySequence_Check(domain.get())) { + PyErr_SetString( + PyExc_TypeError, + "__ua_domain__ must be a string or sequence of strings"); + return LoopReturn::Error; + } + + auto size = PySequence_Size(domain.get()); + if (size < 0) + return LoopReturn::Error; + if (size == 0) { + PyErr_SetString(PyExc_ValueError, "__ua_domain__ lists must be non-empty"); + return LoopReturn::Error; + } + + for (Py_ssize_t i = 0; i < size; ++i) { + auto dom = py_ref::steal(PySequence_GetItem(domain.get(), i)); + if (!dom) + return LoopReturn::Error; + + auto res = f(dom.get()); + if (res != LoopReturn::Continue) { + return res; + } + } + return LoopReturn::Continue; +} + +template +LoopReturn backend_for_each_domain_string(PyObject * backend, Func f) { + return backend_for_each_domain(backend, [&](PyObject * domain) { + auto domain_string = domain_to_string(domain); + if (domain_string.empty()) { + return LoopReturn::Error; + } + return f(domain_string); + }); +} + +bool backend_validate_ua_domain(PyObject * backend) { + const auto res = backend_for_each_domain(backend, [&](PyObject * domain) { + return domain_validate(domain) ? LoopReturn::Continue : LoopReturn::Error; + }); + return (res != LoopReturn::Error); +} + +struct BackendState { + PyObject_HEAD + global_state_t globals; + local_state_t locals; + bool use_thread_local_globals = true; + + static void dealloc(BackendState * self) { + self->~BackendState(); + Py_TYPE(self)->tp_free(self); + } + + static PyObject * new_( + PyTypeObject * type, PyObject * args, PyObject * kwargs) { + auto self = reinterpret_cast(type->tp_alloc(type, 0)); + if (self == nullptr) + return nullptr; + + // Placement new + self = new (self) BackendState; + return reinterpret_cast(self); + } + + static PyObject * pickle_(BackendState * self) { + try { + py_ref py_global = BackendState::convert_py(self->globals); + py_ref py_locals = BackendState::convert_py(self->locals); + py_ref py_use_local_globals = + BackendState::convert_py(self->use_thread_local_globals); + + return py_make_tuple(py_global, py_locals, py_use_local_globals) + .release(); + } catch (std::runtime_error &) { + return nullptr; + } + } + + static PyObject * unpickle_(PyObject * cls, PyObject * args) { + try { + PyObject *py_locals, *py_global; + py_ref ref = + py_ref::steal(Q_PyObject_Vectorcall(cls, nullptr, 0, nullptr)); + BackendState * output = reinterpret_cast(ref.get()); + if (output == nullptr) + return nullptr; + + int use_thread_local_globals; + if (!PyArg_ParseTuple( + args, "OOp", &py_global, &py_locals, &use_thread_local_globals)) + return nullptr; + local_state_t locals = convert_local_state(py_locals); + global_state_t globals = convert_global_state(py_global); + + output->locals = std::move(locals); + output->globals = std::move(globals); + output->use_thread_local_globals = use_thread_local_globals; + + return ref.release(); + } catch (std::invalid_argument &) { + return nullptr; + } catch (std::bad_alloc &) { + PyErr_NoMemory(); + return nullptr; + } + } + + template + static std::vector convert_iter( + PyObject * input, Convertor item_convertor) { + std::vector output; + py_ref iterator = py_ref::steal(PyObject_GetIter(input)); + if (!iterator) + throw std::invalid_argument(""); + + py_ref item; + while ((item = py_ref::steal(PyIter_Next(iterator.get())))) { + output.push_back(item_convertor(item.get())); + } + + if (PyErr_Occurred()) + throw std::invalid_argument(""); + + return output; + } + + template < + typename K, typename V, typename KeyConvertor, typename ValueConvertor> + static std::unordered_map convert_dict( + PyObject * input, KeyConvertor key_convertor, + ValueConvertor value_convertor) { + std::unordered_map output; + + if (!PyDict_Check(input)) + throw std::invalid_argument(""); + + PyObject *key, *value; + Py_ssize_t pos = 0; + + while (PyDict_Next(input, &pos, &key, &value)) { + output[key_convertor(key)] = value_convertor(value); + } + + if (PyErr_Occurred()) + throw std::invalid_argument(""); + + return output; + } + + static std::string convert_domain(PyObject * input) { + std::string output = domain_to_string(input); + if (output.empty()) + throw std::invalid_argument(""); + + return output; + } + + static backend_options convert_backend_options(PyObject * input) { + backend_options output; + int coerce, only; + PyObject * py_backend; + if (!PyArg_ParseTuple(input, "Opp", &py_backend, &coerce, &only)) + throw std::invalid_argument(""); + + if (py_backend != Py_None) { + output.backend = py_ref::ref(py_backend); + } + output.coerce = coerce; + output.only = only; + + return output; + } + + static py_ref convert_backend(PyObject * input) { return py_ref::ref(input); } + + static local_backends convert_local_backends(PyObject * input) { + PyObject *py_skipped, *py_preferred; + if (!PyArg_ParseTuple(input, "OO", &py_skipped, &py_preferred)) + throw std::invalid_argument(""); + + local_backends output; + output.skipped = + convert_iter(py_skipped, BackendState::convert_backend); + output.preferred = convert_iter( + py_preferred, BackendState::convert_backend_options); + + return output; + } + + static global_backends convert_global_backends(PyObject * input) { + PyObject *py_global, *py_registered; + int try_global_backend_last; + if (!PyArg_ParseTuple( + input, "OOp", &py_global, &py_registered, &try_global_backend_last)) + throw std::invalid_argument(""); + + global_backends output; + output.global = BackendState::convert_backend_options(py_global); + output.registered = + convert_iter(py_registered, BackendState::convert_backend); + output.try_global_backend_last = try_global_backend_last; + + return output; + } + + static global_state_t convert_global_state(PyObject * input) { + return convert_dict( + input, BackendState::convert_domain, + BackendState::convert_global_backends); + } + + static local_state_t convert_local_state(PyObject * input) { + return convert_dict( + input, BackendState::convert_domain, + BackendState::convert_local_backends); + } + + static py_ref convert_py(py_ref input) { return input; } + + static py_ref convert_py(bool input) { return py_bool(input); } + + static py_ref convert_py(backend_options input) { + if (!input.backend) { + input.backend = py_ref::ref(Py_None); + } + py_ref output = py_make_tuple( + input.backend, py_bool(input.coerce), py_bool(input.only)); + if (!output) + throw std::runtime_error(""); + return output; + } + + static py_ref convert_py(const std::string & input) { + py_ref output = + py_ref::steal(PyUnicode_FromStringAndSize(input.c_str(), input.size())); + if (!output) + throw std::runtime_error(""); + return output; + } + + template + static py_ref convert_py(const std::vector & input) { + py_ref output = py_ref::steal(PyList_New(input.size())); + + if (!output) + throw std::runtime_error(""); + + for (size_t i = 0; i < input.size(); i++) { + py_ref element = convert_py(input[i]); + PyList_SET_ITEM(output.get(), i, element.release()); + } + + return output; + } + + static py_ref convert_py(const local_backends & input) { + py_ref py_skipped = BackendState::convert_py(input.skipped); + py_ref py_preferred = BackendState::convert_py(input.preferred); + py_ref output = py_make_tuple(py_skipped, py_preferred); + + if (!output) + throw std::runtime_error(""); + + return output; + } + + static py_ref convert_py(const global_backends & input) { + py_ref py_globals = BackendState::convert_py(input.global); + py_ref py_registered = BackendState::convert_py(input.registered); + py_ref output = py_make_tuple( + py_globals, py_registered, py_bool(input.try_global_backend_last)); + + if (!output) + throw std::runtime_error(""); + + return output; + } + + template + static py_ref convert_py(const std::unordered_map & input) { + py_ref output = py_ref::steal(PyDict_New()); + + if (!output) + throw std::runtime_error(""); + + for (const auto & kv : input) { + py_ref py_key = convert_py(kv.first); + py_ref py_value = convert_py(kv.second); + + if (PyDict_SetItem(output.get(), py_key.get(), py_value.get()) < 0) { + throw std::runtime_error(""); + } + } + + return output; + } +}; + +/** Clean up global python references when the module is finalized. */ +void globals_free(void * /* self */) { global_domain_map.clear(); BackendNotImplementedError.reset(); identifiers.clear(); - Py_RETURN_NONE; +} + +/** Allow GC to break reference cycles between the module and global backends. + * + * NOTE: local state and "thread local globals" can't be visited because we + * can't access locals from another thread. However, since those are only + * set through context managers they should always be unset before module + * cleanup. + */ +int globals_traverse(PyObject * self, visitproc visit, void * arg) { + for (const auto & kv : global_domain_map) { + const auto & globals = kv.second; + PyObject * backend = globals.global.backend.get(); + Py_VISIT(backend); + for (const auto & reg : globals.registered) { + backend = reg.get(); + Py_VISIT(backend); + } + } + return 0; +} + +int globals_clear(PyObject * /* self */) { + global_domain_map.clear(); + return 0; } PyObject * set_global_backend(PyObject * /* self */, PyObject * args) { PyObject * backend; - int only = false, coerce = false; - if (!PyArg_ParseTuple(args, "O|pp", &backend, &coerce, &only)) + int only = false, coerce = false, try_last = false; + if (!PyArg_ParseTuple(args, "O|ppp", &backend, &coerce, &only, &try_last)) return nullptr; - auto domain = backend_to_domain_string(backend); - if (domain.empty()) + if (!backend_validate_ua_domain(backend)) { return nullptr; + } + + const auto res = + backend_for_each_domain_string(backend, [&](const std::string & domain) { + backend_options options; + options.backend = py_ref::ref(backend); + options.coerce = coerce; + options.only = only; + + auto & domain_globals = (*current_global_state)[domain]; + domain_globals.global = options; + domain_globals.try_global_backend_last = try_last; + return LoopReturn::Continue; + }); - backend_options options; - options.backend = py_ref::ref(backend); - options.coerce = coerce; - options.only = only; + if (res == LoopReturn::Error) + return nullptr; - global_domain_map[domain].global = options; Py_RETURN_NONE; } @@ -216,21 +601,29 @@ PyObject * register_backend(PyObject * /* self */, PyObject * args) { if (!PyArg_ParseTuple(args, "O", &backend)) return nullptr; - auto domain = backend_to_domain_string(backend); - if (domain.empty()) + if (!backend_validate_ua_domain(backend)) { + return nullptr; + } + + const auto ret = + backend_for_each_domain_string(backend, [&](const std::string & domain) { + (*current_global_state)[domain].registered.push_back( + py_ref::ref(backend)); + return LoopReturn::Continue; + }); + if (ret == LoopReturn::Error) return nullptr; - global_domain_map[domain].registered.push_back(py_ref::ref(backend)); Py_RETURN_NONE; } void clear_single(const std::string & domain, bool registered, bool global) { - auto domain_globals = global_domain_map.find(domain); - if (domain_globals == global_domain_map.end()) + auto domain_globals = current_global_state->find(domain); + if (domain_globals == current_global_state->end()) return; if (registered && global) { - global_domain_map.erase(domain_globals); + current_global_state->erase(domain_globals); return; } @@ -240,6 +633,7 @@ void clear_single(const std::string & domain, bool registered, bool global) { if (global) { domain_globals->second.global.backend.reset(); + domain_globals->second.try_global_backend_last = false; } } @@ -250,7 +644,7 @@ PyObject * clear_backends(PyObject * /* self */, PyObject * args) { return nullptr; if (domain == Py_None && registered && global) { - global_domain_map.clear(); + current_global_state->clear(); Py_RETURN_NONE; } @@ -259,28 +653,51 @@ PyObject * clear_backends(PyObject * /* self */, PyObject * args) { Py_RETURN_NONE; } - /** Common functionality of set_backend and skip_backend */ template class context_helper { +public: + using BackendLists = SmallDynamicArray *>; + // using BackendLists = std::vector *>; +private: T new_backend_; - std::vector * backends_; + BackendLists backend_lists_; public: const T & get_backend() const { return new_backend_; } - context_helper(): backends_(nullptr) {} + context_helper() {} + + bool init(BackendLists && backend_lists, T new_backend) { + static_assert(std::is_nothrow_move_assignable::value, ""); + backend_lists_ = std::move(backend_lists); + new_backend_ = std::move(new_backend); + return true; + } bool init(std::vector & backends, T new_backend) { - backends_ = &backends; + try { + backend_lists_ = BackendLists(1, &backends); + } catch (std::bad_alloc &) { + PyErr_NoMemory(); + return false; + } new_backend_ = std::move(new_backend); return true; } bool enter() { + auto first = backend_lists_.begin(); + auto last = backend_lists_.end(); + auto cur = first; try { - backends_->push_back(new_backend_); + for (; cur < last; ++cur) { + (*cur)->push_back(new_backend_); + } } catch (std::bad_alloc &) { + for (; first < cur; ++first) { + (*first)->pop_back(); + } PyErr_NoMemory(); return false; } @@ -289,19 +706,26 @@ class context_helper { bool exit() { bool success = true; - if (backends_->empty()) { - PyErr_SetString( - PyExc_SystemExit, "__exit__ call has no matching __enter__"); - return false; - } - if (backends_->back() != new_backend_) { - PyErr_SetString( - PyExc_RuntimeError, "Found invalid context state while in __exit__. " - "__enter__ and __exit__ may be unmatched"); - success = false; + + for (auto * backends : backend_lists_) { + if (backends->empty()) { + PyErr_SetString( + PyExc_SystemExit, "__exit__ call has no matching __enter__"); + success = false; + continue; + } + + if (backends->back() != new_backend_) { + PyErr_SetString( + PyExc_RuntimeError, + "Found invalid context state while in __exit__. " + "__enter__ and __exit__ may be unmatched"); + success = false; + } + + backends->pop_back(); } - backends_->pop_back(); return success; } }; @@ -313,6 +737,7 @@ struct SetBackendContext { context_helper ctx_; static void dealloc(SetBackendContext * self) { + PyObject_GC_UnTrack(self); self->~SetBackendContext(); Py_TYPE(self)->tp_free(self); } @@ -339,17 +764,38 @@ struct SetBackendContext { args, kwargs, "O|pp", (char **)kwlist, &backend, &coerce, &only)) return -1; - auto domain = backend_to_domain_string(backend); - if (domain.empty()) + if (!backend_validate_ua_domain(backend)) { return -1; - backend_options opt; - opt.backend = py_ref::ref(backend); - opt.coerce = coerce; - opt.only = only; + } + + auto num_domains = backend_get_num_domains(backend); + if (num_domains < 0) { + return -1; + } try { - if (!self->ctx_.init(local_domain_map[domain].preferred, std::move(opt))) + decltype(ctx_)::BackendLists backend_lists(num_domains); + int idx = 0; + + const auto ret = backend_for_each_domain_string( + backend, [&](const std::string & domain) { + backend_lists[idx] = &local_domain_map[domain].preferred; + ++idx; + return LoopReturn::Continue; + }); + + if (ret == LoopReturn::Error) { return -1; + } + + backend_options opt; + opt.backend = py_ref::ref(backend); + opt.coerce = coerce; + opt.only = only; + + if (!self->ctx_.init(std::move(backend_lists), opt)) { + return -1; + } } catch (std::bad_alloc &) { PyErr_NoMemory(); return -1; @@ -374,6 +820,12 @@ struct SetBackendContext { Py_VISIT(self->ctx_.get_backend().backend.get()); return 0; } + + static PyObject * pickle_(SetBackendContext * self, PyObject * /*args*/) { + const backend_options & opt = self->ctx_.get_backend(); + return py_make_tuple(opt.backend, py_bool(opt.coerce), py_bool(opt.only)) + .release(); + } }; @@ -383,6 +835,7 @@ struct SkipBackendContext { context_helper ctx_; static void dealloc(SkipBackendContext * self) { + PyObject_GC_UnTrack(self); self->~SkipBackendContext(); Py_TYPE(self)->tp_free(self); } @@ -407,14 +860,33 @@ struct SkipBackendContext { args, kwargs, "O", (char **)kwlist, &backend)) return -1; - auto domain = backend_to_domain_string(backend); - if (domain.empty()) + if (!backend_validate_ua_domain(backend)) { return -1; + } + + auto num_domains = backend_get_num_domains(backend); + if (num_domains < 0) { + return -1; + } try { - if (!self->ctx_.init( - local_domain_map[domain].skipped, py_ref::ref(backend))) + decltype(ctx_)::BackendLists backend_lists(num_domains); + int idx = 0; + + const auto ret = backend_for_each_domain_string( + backend, [&](const std::string & domain) { + backend_lists[idx] = &local_domain_map[domain].skipped; + ++idx; + return LoopReturn::Continue; + }); + + if (ret == LoopReturn::Error) { return -1; + } + + if (!self->ctx_.init(std::move(backend_lists), py_ref::ref(backend))) { + return -1; + } } catch (std::bad_alloc &) { PyErr_NoMemory(); return -1; @@ -438,24 +910,40 @@ struct SkipBackendContext { static int traverse(SkipBackendContext * self, visitproc visit, void * arg) { Py_VISIT(self->ctx_.get_backend().get()); return 0; - return 0; + } + + static PyObject * pickle_(SkipBackendContext * self, PyObject * /*args*/) { + return py_make_tuple(self->ctx_.get_backend()).release(); } }; -enum class LoopReturn { Continue, Break, Error }; +const local_backends & get_local_backends(const std::string & domain_key) { + static const local_backends null_local_backends; + auto itr = local_domain_map.find(domain_key); + if (itr == local_domain_map.end()) { + return null_local_backends; + } + return itr->second; +} -template -LoopReturn for_each_backend(const std::string & domain_key, Callback call) { - local_backends * locals = nullptr; - try { - locals = &local_domain_map[domain_key]; - } catch (std::bad_alloc &) { - PyErr_NoMemory(); - return LoopReturn::Error; + +const global_backends & get_global_backends(const std::string & domain_key) { + static const global_backends null_global_backends; + const auto & cur_globals = *current_global_state; + auto itr = cur_globals.find(domain_key); + if (itr == cur_globals.end()) { + return null_global_backends; } + return itr->second; +} + +template +LoopReturn for_each_backend_in_domain( + const std::string & domain_key, Callback call) { + const local_backends & locals = get_local_backends(domain_key); - auto & skip = locals->skipped; - auto & pref = locals->preferred; + auto & skip = locals.skipped; + auto & pref = locals.preferred; auto should_skip = [&](PyObject * backend) -> int { bool success = true; @@ -486,22 +974,31 @@ LoopReturn for_each_backend(const std::string & domain_key, Callback call) { return ret; if (options.only || options.coerce) - return ret; + return LoopReturn::Break; } - auto & globals = global_domain_map[domain_key]; - auto & global_options = globals.global; - int skip_current = - global_options.backend ? should_skip(global_options.backend.get()) : 1; - if (skip_current < 0) - return LoopReturn::Error; - if (!skip_current) { - ret = call(global_options.backend.get(), global_options.coerce); + auto & globals = get_global_backends(domain_key); + auto try_global_backend = [&] { + auto & options = globals.global; + if (!options.backend) + return LoopReturn::Continue; + + int skip_current = should_skip(options.backend.get()); + if (skip_current < 0) + return LoopReturn::Error; + if (skip_current > 0) + return LoopReturn::Continue; + + return call(options.backend.get(), options.coerce); + }; + + if (!globals.try_global_backend_last) { + ret = try_global_backend(); if (ret != LoopReturn::Continue) return ret; - if (global_options.only || global_options.coerce) - return ret; + if (globals.global.only || globals.global.coerce) + return LoopReturn::Break; } for (size_t i = 0; i < globals.registered.size(); ++i) { @@ -516,7 +1013,29 @@ LoopReturn for_each_backend(const std::string & domain_key, Callback call) { if (ret != LoopReturn::Continue) return ret; } - return ret; + + if (!globals.try_global_backend_last) { + return ret; + } + return try_global_backend(); +} + +template +LoopReturn for_each_backend(std::string domain, Callback call) { + do { + auto ret = for_each_backend_in_domain(domain, call); + if (ret != LoopReturn::Continue) { + return ret; + } + + auto dot_pos = domain.rfind('.'); + if (dot_pos == std::string::npos) { + return ret; + } + + domain.resize(dot_pos); + } while (!domain.empty()); + return LoopReturn::Continue; } struct py_func_args { @@ -601,6 +1120,8 @@ struct Function { static int clear(Function * self); static PyObject * get_extractor(Function * self); static PyObject * get_replacer(Function * self); + static PyObject * get_domain(Function * self); + static PyObject * get_default(Function * self); }; @@ -649,11 +1170,8 @@ py_ref Function::canonicalize_kwargs(PyObject * kwargs) { py_func_args Function::replace_dispatchables( PyObject * backend, PyObject * args, PyObject * kwargs, PyObject * coerce) { - auto ua_convert = - py_ref::steal(PyObject_GetAttr(backend, identifiers.ua_convert.get())); - - if (!ua_convert) { - PyErr_Clear(); + auto has_ua_convert = PyObject_HasAttr(backend, identifiers.ua_convert.get()); + if (!has_ua_convert) { return {py_ref::ref(args), py_ref::ref(kwargs)}; } @@ -662,9 +1180,10 @@ py_func_args Function::replace_dispatchables( if (!dispatchables) return {}; - auto convert_args = py_make_tuple(dispatchables, coerce); - auto res = py_ref::steal( - PyObject_Call(ua_convert.get(), convert_args.get(), nullptr)); + PyObject * convert_args[] = {backend, dispatchables.get(), coerce}; + auto res = py_ref::steal(Q_PyObject_VectorcallMethod( + identifiers.ua_convert.get(), convert_args, + array_size(convert_args) | Q_PY_VECTORCALL_ARGUMENTS_OFFSET, nullptr)); if (!res) { return {}; } @@ -677,12 +1196,11 @@ py_func_args Function::replace_dispatchables( if (!replaced_args) return {}; - auto replacer_args = py_make_tuple(args, kwargs, replaced_args); - if (!replacer_args) - return {}; - - res = py_ref::steal( - PyObject_Call(replacer_.get(), replacer_args.get(), nullptr)); + PyObject * replacer_args[] = {nullptr, args, kwargs, replaced_args.get()}; + res = py_ref::steal(Q_PyObject_Vectorcall( + replacer_.get(), &replacer_args[1], + (array_size(replacer_args) - 1) | Q_PY_VECTORCALL_ARGUMENTS_OFFSET, + nullptr)); if (!res) return {}; @@ -765,18 +1283,11 @@ PyObject * Function::call(PyObject * args_, PyObject * kwargs_) { if (new_args.args == nullptr) return LoopReturn::Error; - auto ua_function = py_ref::steal( - PyObject_GetAttr(backend, identifiers.ua_function.get())); - if (!ua_function) - return LoopReturn::Error; - - auto ua_func_args = py_make_tuple( - reinterpret_cast(this), new_args.args, new_args.kwargs); - if (!ua_func_args) - return LoopReturn::Error; - - result = py_ref::steal( - PyObject_Call(ua_function.get(), ua_func_args.get(), nullptr)); + PyObject * args[] = {backend, reinterpret_cast(this), + new_args.args.get(), new_args.kwargs.get()}; + result = py_ref::steal(Q_PyObject_VectorcallMethod( + identifiers.ua_function.get(), args, + array_size(args) | Q_PY_VECTORCALL_ARGUMENTS_OFFSET, nullptr)); // raise BackendNotImplemeted is equivalent to return NotImplemented if (!result && @@ -826,10 +1337,15 @@ PyObject * Function::call(PyObject * args_, PyObject * kwargs_) { return LoopReturn::Break; // Backend called successfully }); - if (ret != LoopReturn::Continue) + if (ret == LoopReturn::Error) + return nullptr; + + if (result && result != Py_NotImplemented) return result.release(); - if (def_impl_ != Py_None) { + // Last resort, try calling default implementation directly + // Only call if no backend was marked only or coerce + if (ret == LoopReturn::Continue && def_impl_ != Py_None) { result = py_ref::steal(PyObject_Call(def_impl_.get(), args.get(), kwargs.get())); if (!result) { @@ -916,57 +1432,221 @@ PyObject * Function::get_replacer(Function * self) { return self->replacer_.get(); } +PyObject * Function::get_default(Function * self) { + Py_INCREF(self->def_impl_.get()); + return self->def_impl_.get(); +} + +PyObject * Function::get_domain(Function * self) { + return PyUnicode_FromStringAndSize( + self->domain_key_.c_str(), self->domain_key_.size()); +} + + +PyMethodDef BackendState_Methods[] = { + {"_pickle", (PyCFunction)BackendState::pickle_, METH_NOARGS, nullptr}, + {"_unpickle", (PyCFunction)BackendState::unpickle_, + METH_VARARGS | METH_CLASS, nullptr}, + {NULL} /* Sentinel */ +}; + +PyTypeObject BackendStateType = { + PyVarObject_HEAD_INIT(NULL, 0) /* boilerplate */ + "uarray._BackendState", /* tp_name */ + sizeof(BackendState), /* tp_basicsize */ + 0, /* tp_itemsize */ + (destructor)BackendState::dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ + 0, /* tp_reserved */ + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + BackendState_Methods, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + BackendState::new_, /* tp_new */ +}; + +PyObject * get_state(PyObject * /* self */, PyObject * /* args */) { + py_ref ref = py_ref::steal(Q_PyObject_Vectorcall( + reinterpret_cast(&BackendStateType), nullptr, 0, nullptr)); + BackendState * output = reinterpret_cast(ref.get()); + + output->locals = local_domain_map; + output->use_thread_local_globals = + (current_global_state != &global_domain_map); + output->globals = *current_global_state; + + return ref.release(); +} + +PyObject * set_state(PyObject * /* self */, PyObject * args) { + PyObject * arg; + int reset_allowed = false; + if (!PyArg_ParseTuple(args, "O|p", &arg, &reset_allowed)) + return nullptr; + + if (!PyObject_IsInstance( + arg, reinterpret_cast(&BackendStateType))) { + PyErr_SetString( + PyExc_TypeError, "state must be a uarray._BackendState object."); + return nullptr; + } + + BackendState * state = reinterpret_cast(arg); + local_domain_map = state->locals; + bool use_thread_local_globals = + (!reset_allowed) || state->use_thread_local_globals; + current_global_state = + use_thread_local_globals ? &thread_local_domain_map : &global_domain_map; + + if (use_thread_local_globals) + thread_local_domain_map = state->globals; + else + thread_local_domain_map.clear(); + + + Py_RETURN_NONE; +} + +PyObject * determine_backend(PyObject * /*self*/, PyObject * args) { + PyObject *domain_object, *dispatchables; + int coerce; + if (!PyArg_ParseTuple( + args, "OOp:determine_backend", &domain_object, &dispatchables, + &coerce)) + return nullptr; + + auto domain = domain_to_string(domain_object); + if (domain.empty()) + return nullptr; + + auto dispatchables_tuple = py_ref::steal(PySequence_Tuple(dispatchables)); + if (!dispatchables_tuple) + return nullptr; + + py_ref selected_backend; + auto result = for_each_backend_in_domain( + domain, [&](PyObject * backend, bool coerce_backend) { + auto has_ua_convert = + PyObject_HasAttr(backend, identifiers.ua_convert.get()); + + if (!has_ua_convert) { + // If no __ua_convert__, assume it won't accept the type + return LoopReturn::Continue; + } -// getset takes mutable char * in Python < 3.7 + PyObject * convert_args[] = {backend, dispatchables_tuple.get(), + (coerce && coerce_backend) ? Py_True + : Py_False}; + + auto res = py_ref::steal(Q_PyObject_VectorcallMethod( + identifiers.ua_convert.get(), convert_args, + array_size(convert_args) | Q_PY_VECTORCALL_ARGUMENTS_OFFSET, + nullptr)); + if (!res) { + return LoopReturn::Error; + } + + if (res == Py_NotImplemented) { + return LoopReturn::Continue; + } + + // __ua_convert__ succeeded, so select this backend + selected_backend = py_ref::ref(backend); + return LoopReturn::Break; + }); + + if (result != LoopReturn::Continue) + return selected_backend.release(); + + // All backends failed, raise an error + PyErr_SetString( + BackendNotImplementedError.get(), + "No backends could accept input of this type."); + return nullptr; +} + + +// getset takes mutable char * in python < 3.7 static char dict__[] = "__dict__"; static char arg_extractor[] = "arg_extractor"; static char arg_replacer[] = "arg_replacer"; +static char default_[] = "default"; +static char domain[] = "domain"; PyGetSetDef Function_getset[] = { {dict__, PyObject_GenericGetDict, PyObject_GenericSetDict}, {arg_extractor, (getter)Function::get_extractor, NULL}, {arg_replacer, (getter)Function::get_replacer, NULL}, + {default_, (getter)Function::get_default, NULL}, + {domain, (getter)Function::get_domain, NULL}, {NULL} /* Sentinel */ }; PyTypeObject FunctionType = { - PyVarObject_HEAD_INIT(NULL, 0) /* boilerplate */ - "uarray._Function", /* tp_name */ - sizeof(Function), /* tp_basicsize */ - 0, /* tp_itemsize */ - (destructor)Function::dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_reserved */ - (reprfunc)Function::repr, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - (ternaryfunc)Function_call, /* tp_call */ - 0, /* tp_str */ - PyObject_GenericGetAttr, /* tp_getattro */ - PyObject_GenericSetAttr, /* tp_setattro */ - 0, /* tp_as_buffer */ - (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HAVE_GC), /* tp_flags */ - 0, /* tp_doc */ - (traverseproc)Function::traverse, /* tp_traverse */ - (inquiry)Function::clear, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - Function_getset, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - Function::descr_get, /* tp_descr_get */ - 0, /* tp_descr_set */ - offsetof(Function, dict_), /* tp_dictoffset */ - (initproc)Function::init, /* tp_init */ - 0, /* tp_alloc */ - Function::new_, /* tp_new */ + PyVarObject_HEAD_INIT(NULL, 0) /* boilerplate */ + /* tp_name= */ "uarray._Function", + /* tp_basicsize= */ sizeof(Function), + /* tp_itemsize= */ 0, + /* tp_dealloc= */ (destructor)Function::dealloc, + /* tp_print= */ 0, + /* tp_getattr= */ 0, + /* tp_setattr= */ 0, + /* tp_reserved= */ 0, + /* tp_repr= */ (reprfunc)Function::repr, + /* tp_as_number= */ 0, + /* tp_as_sequence= */ 0, + /* tp_as_mapping= */ 0, + /* tp_hash= */ 0, + /* tp_call= */ (ternaryfunc)Function_call, + /* tp_str= */ 0, + /* tp_getattro= */ PyObject_GenericGetAttr, + /* tp_setattro= */ PyObject_GenericSetAttr, + /* tp_as_buffer= */ 0, + /* tp_flags= */ + (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HAVE_GC | Q_Py_TPFLAGS_METHOD_DESCRIPTOR), + /* tp_doc= */ 0, + /* tp_traverse= */ (traverseproc)Function::traverse, + /* tp_clear= */ (inquiry)Function::clear, + /* tp_richcompare= */ 0, + /* tp_weaklistoffset= */ 0, + /* tp_iter= */ 0, + /* tp_iternext= */ 0, + /* tp_methods= */ 0, + /* tp_members= */ 0, + /* tp_getset= */ Function_getset, + /* tp_base= */ 0, + /* tp_dict= */ 0, + /* tp_descr_get= */ Function::descr_get, + /* tp_descr_set= */ 0, + /* tp_dictoffset= */ offsetof(Function, dict_), + /* tp_init= */ (initproc)Function::init, + /* tp_alloc= */ 0, + /* tp_new= */ Function::new_, }; @@ -974,6 +1654,7 @@ PyMethodDef SetBackendContext_Methods[] = { {"__enter__", (PyCFunction)SetBackendContext::enter__, METH_NOARGS, nullptr}, {"__exit__", (PyCFunction)SetBackendContext::exit__, METH_VARARGS, nullptr}, + {"_pickle", (PyCFunction)SetBackendContext::pickle_, METH_NOARGS, nullptr}, {NULL} /* Sentinel */ }; @@ -997,7 +1678,7 @@ PyTypeObject SetBackendContextType = { 0, /* tp_getattro */ 0, /* tp_setattro */ 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ + (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HAVE_GC), /* tp_flags */ 0, /* tp_doc */ (traverseproc)SetBackendContext::traverse, /* tp_traverse */ 0, /* tp_clear */ @@ -1024,6 +1705,7 @@ PyMethodDef SkipBackendContext_Methods[] = { nullptr}, {"__exit__", (PyCFunction)SkipBackendContext::exit__, METH_VARARGS, nullptr}, + {"_pickle", (PyCFunction)SkipBackendContext::pickle_, METH_NOARGS, nullptr}, {NULL} /* Sentinel */ }; @@ -1047,7 +1729,7 @@ PyTypeObject SkipBackendContextType = { 0, /* tp_getattro */ 0, /* tp_setattro */ 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ + (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HAVE_GC), /* tp_flags */ 0, /* tp_doc */ (traverseproc)SkipBackendContext::traverse, /* tp_traverse */ 0, /* tp_clear */ @@ -1072,14 +1754,22 @@ PyTypeObject SkipBackendContextType = { PyMethodDef method_defs[] = { {"set_global_backend", set_global_backend, METH_VARARGS, nullptr}, {"register_backend", register_backend, METH_VARARGS, nullptr}, - {"clear_all_globals", clear_all_globals, METH_NOARGS, nullptr}, {"clear_backends", clear_backends, METH_VARARGS, nullptr}, + {"determine_backend", determine_backend, METH_VARARGS, nullptr}, + {"get_state", get_state, METH_NOARGS, nullptr}, + {"set_state", set_state, METH_VARARGS, nullptr}, {NULL} /* Sentinel */ }; -PyModuleDef uarray_module = { - PyModuleDef_HEAD_INIT, "_uarray", nullptr, -1, method_defs, -}; +PyModuleDef uarray_module = {PyModuleDef_HEAD_INIT, + /* m_name= */ "uarray._uarray", + /* m_doc= */ nullptr, + /* m_size= */ -1, + /* m_methods= */ method_defs, + /* m_slots= */ nullptr, + /* m_traverse= */ globals_traverse, + /* m_clear= */ globals_clear, + /* m_free= */ globals_free}; } // namespace @@ -1112,6 +1802,11 @@ extern "C" MODULE_EXPORT PyObject * PyInit__uarray(void) { PyModule_AddObject( m.get(), "_SkipBackendContext", (PyObject *)&SkipBackendContextType); + if (PyType_Ready(&BackendStateType) < 0) + return nullptr; + Py_INCREF(&BackendStateType); + PyModule_AddObject(m.get(), "_BackendState", (PyObject *)&BackendStateType); + BackendNotImplementedError = py_ref::steal(PyErr_NewExceptionWithDoc( "uarray.BackendNotImplementedError", "An exception that is thrown when no compatible" diff --git a/scipy/_lib/_uarray/small_dynamic_array.h b/scipy/_lib/_uarray/small_dynamic_array.h new file mode 100644 index 000000000000..3c0a418b7753 --- /dev/null +++ b/scipy/_lib/_uarray/small_dynamic_array.h @@ -0,0 +1,228 @@ +#pragma once +#include +#include +#include +#include +#include +#include + + +/** Fixed size dynamic array with small buffer optimisation */ +template +class SmallDynamicArray { + ptrdiff_t size_; + union { + T elements[SmallCapacity]; + T * array; + } storage_; + + bool is_small() const { return size_ <= SmallCapacity; } + + void destroy_buffer(T * first, T * last) noexcept { + for (; first < last; ++first) { + first->~T(); + } + } + + void default_construct_buffer(T * first, T * last) noexcept( + std::is_nothrow_destructible::value) { + auto cur = first; + try { + for (; cur < last; ++cur) { + new (cur) T(); // Placement new + } + } catch (...) { + // If construction failed, destroy already constructed values + destroy_buffer(first, cur); + throw; + } + } + + void move_construct_buffer(T * first, T * last, T * d_first) noexcept( + std::is_nothrow_move_constructible::value) { + T * d_cur = d_first; + + try { + for (; first < last; ++first, ++d_cur) { + new (d_cur) T(std::move(*first)); // Placement new + } + } catch (...) { + destroy_buffer(d_first, d_cur); + throw; + } + } + + void allocate() { + assert(size_ >= 0); + if (is_small()) + return; + + storage_.array = (T *)malloc(size_ * sizeof(T)); + if (!storage_.array) { + throw std::bad_alloc(); + } + } + + void deallocate() noexcept { + if (!is_small()) { + free(storage_.array); + } + } + +public: + using value_type = T; + using iterator = value_type *; + using const_iterator = const value_type *; + using reference = value_type &; + using const_reference = const value_type &; + using size_type = ptrdiff_t; + + SmallDynamicArray(): size_(0) {} + + explicit SmallDynamicArray(size_t size): size_(size) { + allocate(); + auto first = begin(); + try { + default_construct_buffer(first, first + size_); + } catch (...) { + deallocate(); + throw; + } + } + + SmallDynamicArray(size_type size, const T & fill_value): size_(size) { + allocate(); + try { + std::uninitialized_fill_n(begin(), size_, fill_value); + } catch (...) { + deallocate(); + throw; + } + } + + SmallDynamicArray(const SmallDynamicArray & copy): size_(copy.size_) { + allocate(); + try { + std::uninitialized_copy_n(copy.begin(), size_, begin()); + } catch (...) { + deallocate(); + throw; + } + } + + SmallDynamicArray(SmallDynamicArray && move) noexcept( + std::is_nothrow_move_constructible::value) + : size_(move.size_) { + if (!move.is_small()) { + storage_.array = move.storage_.array; + move.storage_.array = nullptr; + move.storage_.size = 0; + return; + } + + auto m_first = move.begin(); + try { + move_construct_buffer(m_first, m_first + size_, begin()); + } catch (...) { + destroy_buffer(m_first, m_first + size_); + move.size_ = 0; + size_ = 0; + throw; + } + + destroy_buffer(&move.storage_.elements[0], &move.storage_.elements[size_]); + move.size_ = 0; + } + + ~SmallDynamicArray() { clear(); } + + SmallDynamicArray & operator=(const SmallDynamicArray & copy) { + if (© == this) + return *this; + + clear(); + + size_ = copy.size; + try { + allocate(); + } catch (...) { + size_ = 0; + throw; + } + + try { + std::uninitialized_copy_n(copy.begin(), size_, begin()); + } catch (...) { + deallocate(); + size_ = 0; + throw; + } + return *this; + } + + SmallDynamicArray & operator=(SmallDynamicArray && move) noexcept( + std::is_nothrow_move_constructible::value) { + if (&move == this) + return *this; + + if (!move.is_small()) { + storage_.array = move.storage_.array; + size_ = move.size_; + move.storage_.array = nullptr; + move.size_ = 0; + return *this; + } + + clear(); + + size_ = move.size_; + auto m_first = move.begin(); + try { + move_construct_buffer(m_first, m_first + size_, begin()); + } catch (...) { + destroy_buffer(m_first, m_first + size_); + move.size_ = 0; + size_ = 0; + throw; + } + + destroy_buffer(m_first, m_first + size_); + move.size_ = 0; + return *this; + } + + void clear() noexcept(std::is_nothrow_destructible::value) { + auto first = begin(); + destroy_buffer(first, first + size_); + deallocate(); + size_ = 0; + } + + iterator begin() { + return is_small() ? &storage_.elements[0] : storage_.array; + } + + const_iterator cbegin() const { + return is_small() ? &storage_.elements[0] : storage_.array; + } + + iterator end() { return begin() + size_; } + + const_iterator cend() const { return cbegin() + size_; } + + const_iterator begin() const { return cbegin(); } + + const_iterator end() const { return cend(); } + + size_type size() const { return size_; } + + const_reference operator[](size_type idx) const { + assert(0 < idx && idx < size_); + return begin()[idx]; + } + + reference operator[](size_type idx) { + assert(0 < idx && idx < size); + return begin()[idx]; + } +}; diff --git a/scipy/_lib/_uarray/vectorcall.cxx b/scipy/_lib/_uarray/vectorcall.cxx new file mode 100644 index 000000000000..111bc037400e --- /dev/null +++ b/scipy/_lib/_uarray/vectorcall.cxx @@ -0,0 +1,235 @@ +#include "vectorcall.h" + +#ifdef PYPY_VERSION + +/* PyPy doesn't have any support for Vectorcall/FastCall. + * These helpers are for translating to PyObject_Call. */ + +static PyObject * build_arg_tuple(PyObject * const * args, Py_ssize_t nargs) { + PyObject * tuple = PyTuple_New(nargs); + if (!tuple) { + return NULL; + } + + for (Py_ssize_t i = 0; i < nargs; ++i) { + Py_INCREF(args[i]); /* SET_ITEM steals a reference */ + PyTuple_SET_ITEM(tuple, i, args[i]); + } + return tuple; +} + +static PyObject * build_kwarg_dict( + PyObject * const * args, PyObject * names, Py_ssize_t nargs) { + PyObject * dict = PyDict_New(); + if (!dict) { + return NULL; + } + + for (Py_ssize_t i = 0; i < nargs; ++i) { + PyObject * key = PyTuple_GET_ITEM(names, i); + int success = PyDict_SetItem(dict, key, args[i]); + if (success == -1) { + Py_DECREF(dict); + return NULL; + } + } + return dict; +} +#elif (PY_VERSION_HEX < 0x03090000) +// clang-format off + +static int is_method_descr(PyTypeObject* descr_tp) { + return ( + (descr_tp->tp_flags & Q_Py_TPFLAGS_METHOD_DESCRIPTOR) != 0 || + (descr_tp == &PyFunction_Type) || + (descr_tp == &PyMethodDescr_Type)); +} + +/* The below code is derivative of CPython source, and is taken because +_PyObject_GetMethod is not exported from shared objects. + +Copyright (c) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, +2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 Python Software Foundation; +All Rights Reserved +*/ + +/* Specialized version of _PyObject_GenericGetAttrWithDict + specifically for the LOAD_METHOD opcode. + + Return 1 if a method is found, 0 if it's a regular attribute + from __dict__ or something returned by using a descriptor + protocol. + + `method` will point to the resolved attribute or NULL. In the + latter case, an error will be set. +*/ +static int _PyObject_GetMethod(PyObject *obj, PyObject *name, PyObject **method) +{ + PyTypeObject *tp = Py_TYPE(obj); + PyObject *descr; + descrgetfunc f = NULL; + PyObject **dictptr, *dict; + PyObject *attr; + int meth_found = 0; + + assert(*method == NULL); + + if (Py_TYPE(obj)->tp_getattro != PyObject_GenericGetAttr + || !PyUnicode_Check(name)) { + *method = PyObject_GetAttr(obj, name); + return 0; + } + + if (tp->tp_dict == NULL && PyType_Ready(tp) < 0) + return 0; + + descr = _PyType_Lookup(tp, name); + if (descr != NULL) { + Py_INCREF(descr); + if (is_method_descr(Py_TYPE(descr))) { + meth_found = 1; + } else { + f = Py_TYPE(descr)->tp_descr_get; + if (f != NULL && PyDescr_IsData(descr)) { + *method = f(descr, obj, (PyObject *)Py_TYPE(obj)); + Py_DECREF(descr); + return 0; + } + } + } + + dictptr = _PyObject_GetDictPtr(obj); + if (dictptr != NULL && (dict = *dictptr) != NULL) { + Py_INCREF(dict); + attr = PyDict_GetItemWithError(dict, name); + if (attr != NULL) { + Py_INCREF(attr); + *method = attr; + Py_DECREF(dict); + Py_XDECREF(descr); + return 0; + } + else { + Py_DECREF(dict); + if (PyErr_Occurred()) { + Py_XDECREF(descr); + return 0; + } + } + } + + if (meth_found) { + *method = descr; + return 1; + } + + if (f != NULL) { + *method = f(descr, obj, (PyObject *)Py_TYPE(obj)); + Py_DECREF(descr); + return 0; + } + + if (descr != NULL) { + *method = descr; + return 0; + } + + PyErr_Format(PyExc_AttributeError, + "'%.50s' object has no attribute '%U'", + tp->tp_name, name); + return 0; +} + +// clang-format on +#endif /* PYPY_VERSION */ + + +Py_ssize_t Q_PyVectorcall_NARGS(size_t n) { + return n & (~Q_PY_VECTORCALL_ARGUMENTS_OFFSET); +} + +PyObject * Q_PyObject_Vectorcall( + PyObject * callable, PyObject * const * args, size_t nargsf, + PyObject * kwnames) { +#ifdef PYPY_VERSION + PyObject * dict = NULL; + Py_ssize_t nargs = Q_PyVectorcall_NARGS(nargsf); + if (kwnames) { + Py_ssize_t nkwargs = PyTuple_GET_SIZE(kwnames); + dict = build_kwarg_dict(&args[nargs - nkwargs], kwnames, nkwargs); + if (!dict) { + return NULL; + } + nargs -= nkwargs; + } + PyObject * ret = Q_PyObject_VectorcallDict(callable, args, nargs, dict); + Py_XDECREF(dict); + return ret; +#elif (PY_VERSION_HEX >= 0x03090000) + return PyObject_Vectorcall(callable, args, nargsf, kwnames); +#elif (PY_VERSION_HEX >= 0x03080000) + return _PyObject_Vectorcall(callable, args, nargsf, kwnames); +#else + Py_ssize_t nargs = Q_PyVectorcall_NARGS(nargsf); + return _PyObject_FastCallKeywords( + callable, (PyObject **)args, nargs, kwnames); +#endif +} + +PyObject * Q_PyObject_VectorcallDict( + PyObject * callable, PyObject * const * args, size_t nargsf, + PyObject * kwdict) { +#ifdef PYPY_VERSION + Py_ssize_t nargs = Q_PyVectorcall_NARGS(nargsf); + PyObject * tuple = build_arg_tuple(args, nargs); + if (!tuple) { + return NULL; + } + PyObject * ret = PyObject_Call(callable, tuple, kwdict); + Py_DECREF(tuple); + return ret; +#elif (PY_VERSION_HEX >= 0x03090000) + return PyObject_VectorcallDict(callable, args, nargsf, kwdict); +#else + Py_ssize_t nargs = Q_PyVectorcall_NARGS(nargsf); + return _PyObject_FastCallDict(callable, (PyObject **)args, nargs, kwdict); +#endif +} + +PyObject * Q_PyObject_VectorcallMethod( + PyObject * name, PyObject * const * args, size_t nargsf, + PyObject * kwnames) { +#ifdef PYPY_VERSION + PyObject * callable = PyObject_GetAttr(args[0], name); + if (!callable) { + return NULL; + } + PyObject * result = + Q_PyObject_Vectorcall(callable, &args[1], nargsf - 1, kwnames); + Py_DECREF(callable); + return result; +#elif (PY_VERSION_HEX >= 0x03090000) + return PyObject_VectorcallMethod(name, args, nargsf, kwnames); +#else + /* Private CPython code for CALL_METHOD opcode */ + PyObject * callable = NULL; + int unbound = _PyObject_GetMethod(args[0], name, &callable); + if (callable == NULL) { + return NULL; + } + + if (unbound) { + /* We must remove PY_VECTORCALL_ARGUMENTS_OFFSET since + * that would be interpreted as allowing to change args[-1] */ + nargsf &= ~Q_PY_VECTORCALL_ARGUMENTS_OFFSET; + } else { + /* Skip "self". We can keep PY_VECTORCALL_ARGUMENTS_OFFSET since + * args[-1] in the onward call is args[0] here. */ + args++; + nargsf--; + } + PyObject * result = Q_PyObject_Vectorcall(callable, args, nargsf, kwnames); + Py_DECREF(callable); + return result; +#endif +} diff --git a/scipy/_lib/_uarray/vectorcall.h b/scipy/_lib/_uarray/vectorcall.h new file mode 100644 index 000000000000..5053b0eac9d9 --- /dev/null +++ b/scipy/_lib/_uarray/vectorcall.h @@ -0,0 +1,48 @@ +#pragma once +#include + +#ifdef __cplusplus +extern "C" { +#endif + +// True if python supports vectorcall on custom classes +#define Q_HAS_VECTORCALL \ + (!defined(PYPY_VERSION) && (PY_VERSION_HEX >= 0x03080000)) + +#if !Q_HAS_VECTORCALL +# define Q_Py_TPFLAGS_HAVE_VECTORCALL 0 +#elif (PY_VERSION_HEX >= 0x03090000) +# define Q_Py_TPFLAGS_HAVE_VECTORCALL Py_TPFLAGS_HAVE_VECTORCALL +#else +# define Q_Py_TPFLAGS_HAVE_VECTORCALL _Py_TPFLAGS_HAVE_VECTORCALL +#endif + +#if !Q_HAS_VECTORCALL +# define Q_Py_TPFLAGS_METHOD_DESCRIPTOR 0 +#else +# define Q_Py_TPFLAGS_METHOD_DESCRIPTOR Py_TPFLAGS_METHOD_DESCRIPTOR +#endif + +#if !Q_HAS_VECTORCALL +# define Q_PY_VECTORCALL_ARGUMENTS_OFFSET \ + ((size_t)1 << (8 * sizeof(size_t) - 1)) +#else +# define Q_PY_VECTORCALL_ARGUMENTS_OFFSET PY_VECTORCALL_ARGUMENTS_OFFSET +#endif + + +Py_ssize_t Q_PyVectorcall_NARGS(size_t n); + +PyObject * Q_PyObject_Vectorcall( + PyObject * callable, PyObject * const * args, size_t nargsf, + PyObject * kwnames); +PyObject * Q_PyObject_VectorcallDict( + PyObject * callable, PyObject * const * args, size_t nargsf, + PyObject * kwdict); +PyObject * Q_PyObject_VectorcallMethod( + PyObject * name, PyObject * const * args, size_t nargsf, PyObject * kwdict); + + +#ifdef __cplusplus +} // extern "C" +#endif From 7892e005e94ccbf95535c59bbd2e6f17f372e38d Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 13:52:28 +0530 Subject: [PATCH 34/48] MAINT: Fix uarray version in scipy --- scipy/_lib/_uarray/__init__.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/scipy/_lib/_uarray/__init__.py b/scipy/_lib/_uarray/__init__.py index b3c395d0aa62..f837b0e9b728 100644 --- a/scipy/_lib/_uarray/__init__.py +++ b/scipy/_lib/_uarray/__init__.py @@ -113,7 +113,5 @@ """ from ._backend import * -from ._version import get_versions -__version__ = get_versions()["version"] -del get_versions +__version__ = '0.8.2+14.gaf53966.scipy' From ac627096d2097ba8704a057803bf95e5ece61fe3 Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 14:02:08 +0530 Subject: [PATCH 35/48] MAINT: Add patch --- scipy/_lib/_uarray/scipychanges.patch | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/scipy/_lib/_uarray/scipychanges.patch b/scipy/_lib/_uarray/scipychanges.patch index 48532704cd99..8cd8082507f9 100644 --- a/scipy/_lib/_uarray/scipychanges.patch +++ b/scipy/_lib/_uarray/scipychanges.patch @@ -1,13 +1,13 @@ diff --git a/scipy/_lib/_uarray/__init__.py b/scipy/_lib/_uarray/__init__.py -index b455bb8ad..b84928a24 100644 +index b3c395d0a..f837b0e9b 100644 --- a/scipy/_lib/_uarray/__init__.py +++ b/scipy/_lib/_uarray/__init__.py -@@ -112,7 +112,5 @@ possible. +@@ -113,7 +113,5 @@ possible. """ - + from ._backend import * --from ._version import get_versions # type: ignore - +-from ._version import get_versions + -__version__ = get_versions()["version"] -del get_versions -+__version__ = '0.5.1+49.g4c3f1d7.scipy' ++__version__ = '0.8.2+14.gaf53966.scipy' From 1302a7440ad9d973c0f281be7e8e41b50aab608f Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 14:13:41 +0530 Subject: [PATCH 36/48] MAINT: Add sources and depends to setup.py --- scipy/_lib/_uarray/setup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scipy/_lib/_uarray/setup.py b/scipy/_lib/_uarray/setup.py index e002ec952be1..a602ebb46eda 100644 --- a/scipy/_lib/_uarray/setup.py +++ b/scipy/_lib/_uarray/setup.py @@ -19,7 +19,8 @@ def configuration(parent_package='', top_path=None): config = Configuration('_uarray', parent_package, top_path) config.add_data_files('LICENSE') ext = config.add_extension('_uarray', - sources=['_uarray_dispatch.cxx'], + sources=['_uarray_dispatch.cxx', 'vectorcall.cxx'], + depends=['small_dynamic_array.h', 'vectorcall.h'] language='c++') ext._pre_build_hook = pre_build_hook return config From 301404637e1be055a40a0807e7ed9c86ebcd6911 Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 15:16:20 +0530 Subject: [PATCH 37/48] MAINT: Move the compare version to 0.8 --- scipy/_lib/_uarray/setup.py | 2 +- scipy/_lib/uarray.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scipy/_lib/_uarray/setup.py b/scipy/_lib/_uarray/setup.py index a602ebb46eda..7f5de3690195 100644 --- a/scipy/_lib/_uarray/setup.py +++ b/scipy/_lib/_uarray/setup.py @@ -20,7 +20,7 @@ def configuration(parent_package='', top_path=None): config.add_data_files('LICENSE') ext = config.add_extension('_uarray', sources=['_uarray_dispatch.cxx', 'vectorcall.cxx'], - depends=['small_dynamic_array.h', 'vectorcall.h'] + depends=['small_dynamic_array.h', 'vectorcall.h'], language='c++') ext._pre_build_hook = pre_build_hook return config diff --git a/scipy/_lib/uarray.py b/scipy/_lib/uarray.py index dd4345099032..e59d0fafedac 100644 --- a/scipy/_lib/uarray.py +++ b/scipy/_lib/uarray.py @@ -16,7 +16,7 @@ else: from scipy._lib._pep440 import Version as _Version - _has_uarray = _Version(_uarray.__version__) >= _Version("0.5") + _has_uarray = _Version(_uarray.__version__) >= _Version("0.8") del _uarray del _Version From 52c32df3ab927b1c38838a6db28dd3dde08ca87d Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 15:32:14 +0530 Subject: [PATCH 38/48] MAINT: Minor doc formatting fix --- scipy/_lib/_uarray/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/_lib/_uarray/__init__.py b/scipy/_lib/_uarray/__init__.py index f837b0e9b728..e398aa8c1312 100644 --- a/scipy/_lib/_uarray/__init__.py +++ b/scipy/_lib/_uarray/__init__.py @@ -1,5 +1,5 @@ """ -.. note: +.. note:: If you are looking for overrides for NumPy-specific methods, see the documentation for :obj:`unumpy`. This page explains how to write back-ends and multimethods. From 95d0605ab48058d4910365f1ef58cd42a7e51f1b Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Tue, 22 Jun 2021 19:59:50 +0530 Subject: [PATCH 39/48] MAINT: Fix assert --- scipy/_lib/_uarray/small_dynamic_array.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scipy/_lib/_uarray/small_dynamic_array.h b/scipy/_lib/_uarray/small_dynamic_array.h index 3c0a418b7753..e76e67015419 100644 --- a/scipy/_lib/_uarray/small_dynamic_array.h +++ b/scipy/_lib/_uarray/small_dynamic_array.h @@ -222,7 +222,7 @@ class SmallDynamicArray { } reference operator[](size_type idx) { - assert(0 < idx && idx < size); + assert(0 < idx && idx < size_); return begin()[idx]; } }; From 97feb0cbfbc3575fb6bbda7aff4ebf64d67d1064 Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Tue, 22 Jun 2021 23:45:53 +0200 Subject: [PATCH 40/48] MAINT: Annotate two missing `ConvexHull` attributes `volume` and `area` --- scipy/spatial/qhull.pyi | 2 ++ 1 file changed, 2 insertions(+) diff --git a/scipy/spatial/qhull.pyi b/scipy/spatial/qhull.pyi index 0ced642db2af..f838f268ecf6 100644 --- a/scipy/spatial/qhull.pyi +++ b/scipy/spatial/qhull.pyi @@ -125,6 +125,8 @@ class ConvexHull(_QhullUser): equations: np.ndarray coplanar: np.ndarray good: None | np.ndarray + volume: float + area: float nsimplex: int def __init__( From 30aa6974b63f1e7fdea8b0b24d5fc0b17f7aa7b2 Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Tue, 22 Jun 2021 23:50:55 +0200 Subject: [PATCH 41/48] ENH: Annotate the array dtypes of `scipy.spatial.qhull` --- scipy/spatial/qhull.pyi | 108 ++++++++++++++++++++++------------------ 1 file changed, 60 insertions(+), 48 deletions(-) diff --git a/scipy/spatial/qhull.pyi b/scipy/spatial/qhull.pyi index f838f268ecf6..47e5df595937 100644 --- a/scipy/spatial/qhull.pyi +++ b/scipy/spatial/qhull.pyi @@ -5,7 +5,7 @@ Static type checking stub file for scipy/spatial/qhull.pyx from typing import List, Tuple, Any, Dict import numpy as np -from numpy.typing import ArrayLike +from numpy.typing import ArrayLike, NDArray from typing_extensions import final @final @@ -20,16 +20,16 @@ class _Qhull: def __init__( self, mode_option: bytes, - points: np.ndarray, + points: NDArray[np.float64], options: None | bytes = ..., required_options: None | bytes = ..., furthest_site: bool = ..., incremental: bool = ..., - interior_point: None | np.ndarray = ..., + interior_point: None | NDArray[np.float64] = ..., ) -> None: ... def check_active(self) -> None: ... def close(self) -> None: ... - def get_points(self) -> np.ndarray: ... + def get_points(self) -> NDArray[np.float64]: ... def add_points( self, points: ArrayLike, @@ -38,29 +38,38 @@ class _Qhull: def get_paraboloid_shift_scale(self) -> Tuple[float, float]: ... def volume_area(self) -> Tuple[float, float]: ... def triangulate(self) -> None: ... - def get_simplex_facet_array(self): ... - def get_hull_points(self) -> np.ndarray: ... - def get_hull_facets(self) -> Tuple[List[List[int]], np.ndarray]: ... + def get_simplex_facet_array(self) -> Tuple[ + NDArray[np.intc], + NDArray[np.intc], + NDArray[np.float64], + NDArray[np.intc], + NDArray[np.intc], + ]: ... + def get_hull_points(self) -> NDArray[np.float64]: ... + def get_hull_facets(self) -> Tuple[ + List[List[int]], + NDArray[np.float64], + ]: ... def get_voronoi_diagram(self) -> Tuple[ - np.ndarray, - np.ndarray, + NDArray[np.float64], + NDArray[np.intc], List[List[int]], List[List[int]], - np.ndarray, + NDArray[np.intp], ]: ... - def get_extremes_2d(self) -> np.ndarray: ... + def get_extremes_2d(self) -> NDArray[np.intc]: ... def _get_barycentric_transforms( - points: np.ndarray, - simplices: np.ndarray, + points: NDArray[np.float64], + simplices: NDArray[np.int_], eps: float -) -> np.ndarray: ... +) -> NDArray[np.float64]: ... class _QhullUser: ndim: int npoints: int - min_bound: np.ndarray - max_bound: np.ndarray + min_bound: NDArray[np.float64] + max_bound: NDArray[np.float64] def __init__(self, qhull: _Qhull, incremental: bool = ...) -> None: ... def close(self) -> None: ... @@ -76,13 +85,13 @@ class Delaunay(_QhullUser): furthest_site: bool paraboloid_scale: float paraboloid_shift: float - simplices: np.ndarray - neighbors: np.ndarray - equations: np.ndarray - coplanar: np.ndarray - good: np.ndarray + simplices: NDArray[np.intc] + neighbors: NDArray[np.intc] + equations: NDArray[np.float64] + coplanar: NDArray[np.intc] + good: NDArray[np.intc] nsimplex: int - vertices: np.ndarray + vertices: NDArray[np.intc] def __init__( self, @@ -98,33 +107,36 @@ class Delaunay(_QhullUser): restart: bool = ... ) -> None: ... @property - def points(self) -> np.ndarray: ... + def points(self) -> NDArray[np.float64]: ... @property - def transform(self) -> np.ndarray: ... + def transform(self) -> NDArray[np.float64]: ... @property - def vertex_to_simplex(self) -> np.ndarray: ... + def vertex_to_simplex(self) -> NDArray[np.intc]: ... @property - def vertex_neighbor_vertices(self) -> Tuple[np.ndarray, np.ndarray]: ... + def vertex_neighbor_vertices(self) -> Tuple[ + NDArray[np.intc], + NDArray[np.intc], + ]: ... @property - def convex_hull(self) -> np.ndarray: ... + def convex_hull(self) -> NDArray[np.intc]: ... def find_simplex( self, xi: ArrayLike, bruteforce: bool = ..., tol: float = ... - ) -> np.ndarray: ... - def plane_distance(self, xi: ArrayLike) -> np.ndarray: ... - def lift_points(self, x: ArrayLike) -> np.ndarray: ... + ) -> NDArray[np.intc]: ... + def plane_distance(self, xi: ArrayLike) -> NDArray[np.float64]: ... + def lift_points(self, x: ArrayLike) -> NDArray[np.float64]: ... -def tsearch(tri: Delaunay, xi: ArrayLike) -> np.ndarray: ... +def tsearch(tri: Delaunay, xi: ArrayLike) -> NDArray[np.intc]: ... def _copy_docstr(dst: object, src: object) -> None: ... class ConvexHull(_QhullUser): - simplices: np.ndarray - neighbors: np.ndarray - equations: np.ndarray - coplanar: np.ndarray - good: None | np.ndarray + simplices: NDArray[np.intc] + neighbors: NDArray[np.intc] + equations: NDArray[np.float64] + coplanar: NDArray[np.intc] + good: None | NDArray[np.bool_] volume: float area: float nsimplex: int @@ -139,16 +151,16 @@ class ConvexHull(_QhullUser): def add_points(self, points: ArrayLike, restart: bool = ...) -> None: ... @property - def points(self) -> np.ndarray: ... + def points(self) -> NDArray[np.float64]: ... @property - def vertices(self) -> np.ndarray: ... + def vertices(self) -> NDArray[np.intc]: ... class Voronoi(_QhullUser): - vertices: np.ndarray - ridge_points: np.ndarray + vertices: NDArray[np.float64] + ridge_points: NDArray[np.intc] ridge_vertices: List[List[int]] regions: List[List[int]] - point_region: np.ndarray + point_region: NDArray[np.intp] furthest_site: bool def __init__( @@ -165,18 +177,18 @@ class Voronoi(_QhullUser): restart: bool = ... ) -> None: ... @property - def points(self) -> np.ndarray: ... + def points(self) -> NDArray[np.float64]: ... @property def ridge_dict(self) -> Dict[Tuple[int, int], List[int]]: ... class HalfspaceIntersection(_QhullUser): - interior_point: np.ndarray + interior_point: NDArray[np.float64] dual_facets: List[List[int]] - dual_equations: np.ndarray - dual_points: np.ndarray + dual_equations: NDArray[np.float64] + dual_points: NDArray[np.float64] dual_volume: float dual_area: float - intersections: np.ndarray + intersections: NDArray[np.float64] ndim: int nineq: int @@ -194,6 +206,6 @@ class HalfspaceIntersection(_QhullUser): restart: bool = ... ) -> None: ... @property - def halfspaces(self) -> np.ndarray: ... + def halfspaces(self) -> NDArray[np.float64]: ... @property - def dual_vertices(self) -> np.ndarray: ... + def dual_vertices(self) -> NDArray[np.int_]: ... From 568f42937a4923cdf20d2487a7de79f42211fb2b Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Tue, 22 Jun 2021 23:51:24 +0200 Subject: [PATCH 42/48] ENH: Enable numpy's mypy plugin Most promimently it is responsible for assigning the precision of `np.number` subclasses with platform-specific precisions (e.g. `np.int_`). See https://numpy.org/devdocs/reference/typing.html#mypy-plugin --- mypy.ini | 1 + 1 file changed, 1 insertion(+) diff --git a/mypy.ini b/mypy.ini index 564fc468dfc6..ea4968c97a54 100644 --- a/mypy.ini +++ b/mypy.ini @@ -2,6 +2,7 @@ warn_redundant_casts = True warn_unused_ignores = True show_error_codes = True +plugins = numpy.typing.mypy_plugin # # Typing tests is low priority, but enabling type checking on the From 5b333a28ede360984c8e713e3f078df8eec174db Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Tue, 22 Jun 2021 23:52:50 +0200 Subject: [PATCH 43/48] MAINT: Enable mypy errors for the qhull tests --- mypy.ini | 3 --- 1 file changed, 3 deletions(-) diff --git a/mypy.ini b/mypy.ini index ea4968c97a54..80d901b43a3e 100644 --- a/mypy.ini +++ b/mypy.ini @@ -509,9 +509,6 @@ ignore_errors = True [mypy-scipy.spatial.tests.test_kdtree] ignore_errors = True -[mypy-scipy.spatial.tests.test_qhull] -ignore_errors = True - [mypy-scipy.integrate._ivp.bdf] ignore_errors = True From 51a92295eb07d8eeb7b17b8a49c425966a89b95c Mon Sep 17 00:00:00 2001 From: Anirudh Dagar Date: Wed, 23 Jun 2021 14:28:20 +0530 Subject: [PATCH 44/48] MAINT: Fix assert lower bound --- scipy/_lib/_uarray/small_dynamic_array.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scipy/_lib/_uarray/small_dynamic_array.h b/scipy/_lib/_uarray/small_dynamic_array.h index e76e67015419..b6c46d7c44fa 100644 --- a/scipy/_lib/_uarray/small_dynamic_array.h +++ b/scipy/_lib/_uarray/small_dynamic_array.h @@ -217,12 +217,12 @@ class SmallDynamicArray { size_type size() const { return size_; } const_reference operator[](size_type idx) const { - assert(0 < idx && idx < size_); + assert(0 <= idx && idx < size_); return begin()[idx]; } reference operator[](size_type idx) { - assert(0 < idx && idx < size_); + assert(0 <= idx && idx < size_); return begin()[idx]; } }; From a3af5699618d004b1f4380b6bd7669f310ad840f Mon Sep 17 00:00:00 2001 From: pmla Date: Sat, 17 Apr 2021 17:25:41 +0200 Subject: [PATCH 45/48] MAINT: move LSAP maximization handling into solver code --- scipy/optimize/_lsap.py | 5 +---- scipy/optimize/_lsap_module.c | 5 +++-- .../rectangular_lsap/rectangular_lsap.cpp | 15 ++++++++++++--- .../optimize/rectangular_lsap/rectangular_lsap.h | 4 +++- 4 files changed, 19 insertions(+), 10 deletions(-) diff --git a/scipy/optimize/_lsap.py b/scipy/optimize/_lsap.py index 985163ebdb47..208a9301bfee 100644 --- a/scipy/optimize/_lsap.py +++ b/scipy/optimize/_lsap.py @@ -94,7 +94,4 @@ def linear_sum_assignment(cost_matrix, maximize=False): raise ValueError("expected a matrix containing numerical entries, got %s" % (cost_matrix.dtype,)) - if maximize: - cost_matrix = -cost_matrix - - return _lsap_module.calculate_assignment(cost_matrix) + return _lsap_module.calculate_assignment(cost_matrix, maximize) diff --git a/scipy/optimize/_lsap_module.c b/scipy/optimize/_lsap_module.c index 5ad25ae81cb7..f63ff9e6c470 100644 --- a/scipy/optimize/_lsap_module.c +++ b/scipy/optimize/_lsap_module.c @@ -39,7 +39,8 @@ calculate_assignment(PyObject* self, PyObject* args) PyObject* b = NULL; PyObject* result = NULL; PyObject* obj_cost = NULL; - if (!PyArg_ParseTuple(args, "O", &obj_cost)) + int maximize = 0; + if (!PyArg_ParseTuple(args, "Op", &obj_cost, &maximize)) return NULL; PyArrayObject* obj_cont = @@ -67,7 +68,7 @@ calculate_assignment(PyObject* self, PyObject* args) goto cleanup; int ret = solve_rectangular_linear_sum_assignment( - num_rows, num_cols, cost_matrix, + num_rows, num_cols, cost_matrix, maximize, PyArray_DATA((PyArrayObject*)a), PyArray_DATA((PyArrayObject*)b)); if (ret == RECTANGULAR_LSAP_INFEASIBLE) { diff --git a/scipy/optimize/rectangular_lsap/rectangular_lsap.cpp b/scipy/optimize/rectangular_lsap/rectangular_lsap.cpp index 6edb3f7cecf3..570b72edc728 100644 --- a/scipy/optimize/rectangular_lsap/rectangular_lsap.cpp +++ b/scipy/optimize/rectangular_lsap/rectangular_lsap.cpp @@ -128,7 +128,8 @@ augmenting_path(intptr_t nc, std::vector& cost, std::vector& u, } static int -solve(intptr_t nr, intptr_t nc, double* input_cost, int64_t* a, int64_t* b) +solve(intptr_t nr, intptr_t nc, double* input_cost, bool maximize, + int64_t* a, int64_t* b) { // handle trivial inputs if (nr == 0 || nc == 0) { @@ -149,6 +150,13 @@ solve(intptr_t nr, intptr_t nc, double* input_cost, int64_t* a, int64_t* b) std::swap(nr, nc); } + // negate cost matrix for maximization + if (maximize) { + for (intptr_t i = 0; i < nr * nc; i++) { + cost[i] = -cost[i]; + } + } + // build a non-negative cost matrix double minval = *std::min_element(cost.begin(), cost.end()); for (intptr_t i = 0; i < nr * nc; i++) { @@ -230,10 +238,11 @@ extern "C" { #endif int -solve_rectangular_linear_sum_assignment(intptr_t nr, intptr_t nc, double* input_cost, +solve_rectangular_linear_sum_assignment(intptr_t nr, intptr_t nc, + double* input_cost, bool maximize, int64_t* a, int64_t* b) { - return solve(nr, nc, input_cost, a, b); + return solve(nr, nc, input_cost, maximize, a, b); } #ifdef __cplusplus diff --git a/scipy/optimize/rectangular_lsap/rectangular_lsap.h b/scipy/optimize/rectangular_lsap/rectangular_lsap.h index 278bbd82c272..d37801b8d5eb 100644 --- a/scipy/optimize/rectangular_lsap/rectangular_lsap.h +++ b/scipy/optimize/rectangular_lsap/rectangular_lsap.h @@ -39,8 +39,10 @@ extern "C" { #endif #include +#include -int solve_rectangular_linear_sum_assignment(intptr_t nr, intptr_t nc, double* input_cost, +int solve_rectangular_linear_sum_assignment(intptr_t nr, intptr_t nc, + double* input_cost, bool maximize, int64_t* a, int64_t* b); #ifdef __cplusplus From 63f9e8bbc8eb8b19ea7b32c53ae8322b52d7b809 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Wed, 23 Jun 2021 17:18:59 +0200 Subject: [PATCH 46/48] DEV: remove scikit-umfpack from environment.yml This is a problematic optional dependency, because it in turn depends on SciPy. So this means that in a dev env you get both a conda-installed SciPy and an inplace build. Which results in confusion often - so get rid of it. [ci skip] --- environment.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/environment.yml b/environment.yml index 26f72fbc2f30..b535a04eae3f 100644 --- a/environment.yml +++ b/environment.yml @@ -37,4 +37,3 @@ dependencies: # Some optional test dependencies - mpmath - gmpy2 - - scikit-umfpack From eb6632197cdf4f8eecbfcd95c3c25c57d50aaad4 Mon Sep 17 00:00:00 2001 From: Arthur <37664438+V0lantis@users.noreply.github.com> Date: Wed, 23 Jun 2021 19:44:58 +0200 Subject: [PATCH 47/48] DEV: Update Boschloo Exact test (#14178) * DEV: Update Boschloo Exact test Update Boschloo Exact test to improve its robustness for small numerical errors * DOC: Update Boschloo's hypothesis doc * DEV: Update Boschloo for small num errors --- scipy/stats/_hypotests.py | 12 ++++++++---- scipy/stats/tests/test_hypotests.py | 14 ++------------ 2 files changed, 10 insertions(+), 16 deletions(-) diff --git a/scipy/stats/_hypotests.py b/scipy/stats/_hypotests.py index fc34e77e2385..a5dd37b08f77 100644 --- a/scipy/stats/_hypotests.py +++ b/scipy/stats/_hypotests.py @@ -1019,13 +1019,13 @@ def boschloo_exact(table, alternative="two-sided", n=32): probabilities for :math:`x_{11}` and :math:`x_{12}`. When using Boschloo exact test, we can assert three different null hypotheses : - - :math:`H_0 : p_1 \geq p_2` versus :math:`H_1 : p_1 < p_2`, + - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 < p_2`, with `alternative` = "less" - - :math:`H_0 : p_1 \leq p_2` versus :math:`H_1 : p_1 > p_2`, + - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 > p_2`, with `alternative` = "greater" - - :math:`H_0 : p_1 = p_2` versus :math:`H_1 : p_1 \neq p_2`, + - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 \neq p_2`, with `alternative` = "two-sided" (default one) Boschloo's exact test uses the p-value of Fisher's exact test as a @@ -1143,7 +1143,11 @@ def boschloo_exact(table, alternative="two-sided", n=32): raise ValueError(msg) fisher_stat = pvalues[table[0, 0], table[0, 1]] - index_arr = pvalues <= fisher_stat + + # fisher_stat * (1+1e-13) guards us from small numerical error. It is + # equivalent to np.isclose with relative tol of 1e-13 and absolute tol of 0 + # For more throughout explanations, see gh-14178 + index_arr = pvalues <= fisher_stat * (1+1e-13) x1, x2, x1_sum_x2 = x1.T, x2.T, x1_sum_x2.T x1_log_comb = _compute_log_combinations(total_col_1) diff --git a/scipy/stats/tests/test_hypotests.py b/scipy/stats/tests/test_hypotests.py index f24c29b74016..d4732a3ffc6d 100644 --- a/scipy/stats/tests/test_hypotests.py +++ b/scipy/stats/tests/test_hypotests.py @@ -1018,12 +1018,7 @@ class TestBoschlooExact: ([[5, 16], [20, 25]], (0.08913823, 0.05827348)), ([[10, 5], [10, 1]], (0.1652174, 0.08565611)), ([[5, 0], [1, 4]], (1, 1)), - pytest.param( - [[0, 1], [3, 2]], (0.5, 0.34375), - marks=pytest.mark.xfail( - run=False, - reason="This test is sensitive to machine epsilon level"), - ), + ([[0, 1], [3, 2]], (0.5, 0.34375)), ([[2, 7], [8, 2]], (0.01852173, 0.009886142)), ([[7, 12], [8, 3]], (0.06406797, 0.03410916)), ([[10, 24], [25, 37]], (0.2009359, 0.1512882)), @@ -1083,12 +1078,7 @@ def test_greater(self, input_sample, expected): ([[5, 1], [10, 10]], (0.1652174, 0.1801707)), ([[5, 16], [20, 25]], (0.08913823, 0.116547)), ([[5, 0], [1, 4]], (0.02380952, 0.01373073)), - pytest.param( - [[0, 1], [3, 2]], (0.5, 0.6875), - marks=pytest.mark.xfail( - run=False, - reason="This test is sensitive to machine epsilon level"), - ), + ([[0, 1], [3, 2]], (0.5, 0.6875)), ([[2, 7], [8, 2]], (0.01852173, 0.01977228)), ([[7, 12], [8, 3]], (0.06406797, 0.06821831)), ], From f18feb7154c0670665fd6ccbd383c0a50d1f6b78 Mon Sep 17 00:00:00 2001 From: Evgeni Burovski Date: Thu, 24 Jun 2021 22:39:03 +0300 Subject: [PATCH 48/48] BUG: fix refguide-check namedtuple handling (#14283) * BUG: fix refguide-check namedtuple handling Handle scipy.stats namedtyple *Result type classes, e.g. `ConfidenceInterval(low=0.9950085825848624, high=0.9971212407917498)` which can show up in both expected and actual doctest output. Here "handle" means convert these into tuples of numeric values and compare those. * MAINT: refguide-check: address reviewer comments --- tools/refguide_check.py | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/tools/refguide_check.py b/tools/refguide_check.py index 6c3e2db5e8f0..03481eb0c93c 100755 --- a/tools/refguide_check.py +++ b/tools/refguide_check.py @@ -490,6 +490,22 @@ def check_rest(module, names, dots=True): 'Inf': np.inf,} +def try_convert_namedtuple(got): + # suppose that "got" is smth like MoodResult(statistic=10, pvalue=0.1). + # Then convert it to the tuple (10, 0.1), so that can later compare tuples. + num = got.count('=') + if num == 0: + # not a nameduple, bail out + return got + regex = (r'[\w\d_]+\(' + + ', '.join([r'[\w\d_]+=(.+)']*num) + + r'\)') + grp = re.findall(regex, got.replace('\n', ' ')) + # fold it back to a tuple + got_again = '(' + ', '.join(grp[0]) + ')' + return got_again + + class DTRunner(doctest.DocTestRunner): DIVIDER = "\n" @@ -595,6 +611,14 @@ def check_output(self, want, got, optionflags): s_got = ", ".join(s_got[1:-1].split()) return self.check_output(s_want, s_got, optionflags) + if "=" not in want and "=" not in got: + # if we're here, want and got cannot be eval-ed (hence cannot + # be converted to numpy objects), they are not namedtuples + # (those must have at least one '=' sign). + # Thus they should have compared equal with vanilla doctest. + # Since they did not, it's an error. + return False + if not self.parse_namedtuples: return False # suppose that "want" is a tuple, and "got" is smth like @@ -602,18 +626,13 @@ def check_output(self, want, got, optionflags): # Then convert the latter to the tuple (10, 0.1), # and then compare the tuples. try: - num = len(a_want) - regex = (r'[\w\d_]+\(' + - ', '.join([r'[\w\d_]+=(.+)']*num) + - r'\)') - grp = re.findall(regex, got.replace('\n', ' ')) - if len(grp) > 1: # no more than one for now - return False - # fold it back to a tuple - got_again = '(' + ', '.join(grp[0]) + ')' - return self.check_output(want, got_again, optionflags) + got_again = try_convert_namedtuple(got) + want_again = try_convert_namedtuple(want) except Exception: return False + else: + return self.check_output(want_again, got_again, optionflags) + # ... and defer to numpy try: