File numpy2.patch of Package failed_python-sherpa
From 72028ffe7ce2566a8f1e88c2c06d79cf5f0be9c1 Mon Sep 17 00:00:00 2001
From: Douglas Burke <dburke.gw@gmail.com>
Date: Thu, 27 Jun 2024 12:42:52 -0400
Subject: [PATCH 1/7] root: internal code cleanup
The root-finding code is not documented well. This adds a small
wrapper routine to avoid some replicated code, but could we
just add this to transformed_quad_coef() instead - which is
not explicitly marked as an external routine?
Several comments have been added for potential future work.
---
sherpa/utils/__init__.py | 38 ++++++++++++++++++++++-----------
sherpa/utils/tests/test_root.py | 5 +++++
2 files changed, 30 insertions(+), 13 deletions(-)
Index: sherpa-4.16.1/sherpa/utils/__init__.py
===================================================================
--- sherpa-4.16.1.orig/sherpa/utils/__init__.py
+++ sherpa-4.16.1/sherpa/utils/__init__.py
@@ -1480,7 +1480,7 @@ def create_expr_integrated(lovals, hival
delim : str, optional
The separator for a range.
eps : number, optional
- The tolerance for comparing two numbers with sao_fcmp.
+ This value is unused.
Raises
------
@@ -3389,6 +3389,7 @@ def bisection(fcn, xa, xb, fa=None, fb=N
return [[None, None], [[xa, fa], [xb, fb]], nfev[0]]
+# Is this used at all?
def quad_coef(x, f):
"""
p( x ) = f( xc ) + A ( x - xc ) + B ( x - xc ) ( x - xb )
@@ -3461,6 +3462,11 @@ def transformed_quad_coef(x, f):
xa, xb, xc = x[0], x[1], x[2]
fa, fb, fc = f[0], f[1], f[2]
+ # What happens if xb_xa or xc_xa are 0? That is, either
+ # xa == xb
+ # xc == xa
+ # Is the assumption that this just never happen?
+ #
xc_xb = xc - xb
fc_fb = fc - fb
A = fc_fb / xc_xb
@@ -3472,6 +3478,21 @@ def transformed_quad_coef(x, f):
return [B, C]
+def _get_discriminant(xa, xb, xc, fa, fb, fc):
+ """Wrap up code to transformed_quad_coef.
+
+ This is common code that could be added to transformed_quad_coef
+ but is left out at the moment, to make it easier to look back
+ at code changes. There is no description of the parameters as
+ the existing code has none.
+
+ """
+
+ [B, C] = transformed_quad_coef([xa, xb, xc], [fa, fb, fc])
+ discriminant = max(C * C - 4.0 * fc * B, 0.0)
+ return B, C, discriminant
+
+
def demuller(fcn, xa, xb, xc, fa=None, fb=None, fc=None, args=(),
maxfev=32, tol=1.0e-6):
"""A root-finding algorithm using Muller's method.
@@ -3578,10 +3599,7 @@ def demuller(fcn, xa, xb, xc, fa=None, f
while nfev[0] < maxfev:
- [B, C] = transformed_quad_coef([xa, xb, xc], [fa, fb, fc])
-
- discriminant = max(C * C - 4.0 * fc * B, 0.0)
-
+ B, C, discriminant = _get_discriminant(xa, xb, xc, fa, fb, fc)
if is_nan(B) or is_nan(C) or \
0.0 == C + mysgn(C) * np.sqrt(discriminant):
return [[None, None], [[None, None], [None, None]], nfev[0]]
@@ -3685,11 +3703,7 @@ def new_muller(fcn, xa, xb, fa=None, fb=
if abs(fc) <= tol:
return [[xc, fc], [[xa, fa], [xb, fb]], nfev[0]]
- tran = transformed_quad_coef([xa, xb, xc], [fa, fb, fc])
- B = tran[0]
- C = tran[1]
-
- discriminant = max(C * C - 4.0 * fc * B, 0.0)
+ B, C, discriminant = _get_discriminant(xa, xb, xc, fa, fb, fc)
xd = xc - 2.0 * fc / (C + mysgn(C) * np.sqrt(discriminant))
@@ -3827,11 +3841,9 @@ def apache_muller(fcn, xa, xb, fa=None,
oldx = 1.0e128
while nfev[0] < maxfev:
- tran = transformed_quad_coef([xa, xb, xc], [fa, fb, fc])
- B = tran[0]
- C = tran[1]
- discriminant = max(C * C - 4.0 * fc * B, 0.0)
- den = mysgn(C) * np.sqrt(discriminant)
+
+ B, C, discriminant = _get_discriminant(xa, xb, xc, fa, fb, fc)
+ den = np.sign(C) * np.sqrt(discriminant)
xplus = xc - 2.0 * fc / (C + den)
if C != den:
xminus = xc - 2.0 * fc / (C - den)
@@ -4008,9 +4020,13 @@ def zeroin(fcn, xa, xb, fa=None, fb=None
warning('%s: %s fa * fb < 0 is not met', __name__, fcn.__name__)
return [[None, None], [[None, None], [None, None]], nfev[0]]
+ # With NumPy 2.0 the casting rules changed, leading to some
+ # behavioural changes in this code. The simplest fix was to
+ # make sure DBL_EPSILON did not remain a np.float32 value.
+ #
xc = xa
fc = fa
- DBL_EPSILON = np.finfo(np.float32).eps
+ DBL_EPSILON = float(np.finfo(np.float32).eps)
while nfev[0] < maxfev:
prev_step = xb - xa
Index: sherpa-4.16.1/sherpa/utils/tests/test_root.py
===================================================================
--- sherpa-4.16.1.orig/sherpa/utils/tests/test_root.py
+++ sherpa-4.16.1/sherpa/utils/tests/test_root.py
@@ -1,5 +1,6 @@
#
-# Copyright (C) 2007, 2016, 2018, 2020, 2021 Smithsonian Astrophysical Observatory
+# Copyright (C) 2007, 2016, 2018, 2020, 2021, 2024
+# Smithsonian Astrophysical Observatory
#
#
# This program is free software; you can redistribute it and/or modify
@@ -27,7 +28,7 @@ from sherpa.utils import demuller, bisec
zeroin
-def sqr(x, *args):
+def sqr(x):
return x * x
@@ -177,9 +178,7 @@ def prob34(x, *args):
return 1.0 / x - numpy.sin(x) + 1.0
-def prob35(x, *args):
- return (x*x - 2.0) * x - 5.0
-
+# prob35 was the same as prob16
def prob36(x, *args):
return 1.0 / x - 1.0
@@ -288,7 +287,6 @@ def demuller2(fcn, xa, xb, fa=None, fb=N
(prob32, 0.1, 0.9),
(prob33, 2.8, 3.1),
(prob34, -1.3, -0.5),
- (prob35, 2.0, 3.0),
(prob36, 0.5, 1.5),
(prob37, 0.5, 5.0),
(prob38, 1.0, 4.0),
Index: sherpa-4.16.1/sherpa/estmethods/__init__.py
===================================================================
--- sherpa-4.16.1.orig/sherpa/estmethods/__init__.py
+++ sherpa-4.16.1/sherpa/estmethods/__init__.py
@@ -380,6 +380,11 @@ def covariance(pars, parmins, parmaxes,
eflag = est_success
ubound = diag[num]
lbound = -diag[num]
+
+ # What happens when lbound or ubound is NaN? This is
+ # presumably why the code is written as it is below (e.g. a
+ # pass if the values can be added to pars[num]).
+ #
if pars[num] + ubound < parhardmaxes[num]:
pass
else:
@@ -1093,6 +1098,7 @@ def confidence(pars, parmins, parmaxes,
print_status(myblog.blogger.info, verbose, status_prefix[dirn],
delta_zero, lock)
+ # This should really set the error flag appropriately.
error_flags.append(est_success)
#
Index: sherpa-4.16.1/sherpa/fit.py
===================================================================
--- sherpa-4.16.1.orig/sherpa/fit.py
+++ sherpa-4.16.1/sherpa/fit.py
@@ -277,7 +277,7 @@ class FitResults(NoNewAttributesAfterIni
self.succeeded = results[0]
self.parnames = tuple(p.fullname for p in fit.model.get_thawed_pars())
- self.parvals = tuple(results[1])
+ self.parvals = tuple(float(r) for r in results[1])
self.istatval = init_stat
self.statval = results[2]
self.dstatval = np.abs(results[2] - init_stat)
@@ -439,25 +439,28 @@ class ErrorEstResults(NoNewAttributesAft
self.sigma = fit.estmethod.sigma
self.percent = erf(self.sigma / sqrt(2.0)) * 100.0
self.parnames = tuple(p.fullname for p in parlist if not p.frozen)
- self.parvals = tuple(p.val for p in parlist if not p.frozen)
+ self.parvals = tuple(float(p.val) for p in parlist if not p.frozen)
self.parmins = ()
self.parmaxes = ()
- self.nfits = 0
for i in range(len(parlist)):
if (results[2][i] == est_hardmin or
- results[2][i] == est_hardminmax):
+ results[2][i] == est_hardminmax or
+ results[0][i] is None # It looks like confidence does not set the flag
+ ):
self.parmins = self.parmins + (None,)
warning("hard minimum hit for parameter %s", self.parnames[i])
else:
- self.parmins = self.parmins + (results[0][i],)
+ self.parmins = self.parmins + (float(results[0][i]),)
if (results[2][i] == est_hardmax or
- results[2][i] == est_hardminmax):
+ results[2][i] == est_hardminmax or
+ results[1][i] is None # It looks like confidence does not set the flag
+ ):
self.parmaxes = self.parmaxes + (None,)
warning("hard maximum hit for parameter %s", self.parnames[i])
else:
- self.parmaxes = self.parmaxes + (results[1][i],)
+ self.parmaxes = self.parmaxes + (float(results[1][i]),)
self.nfits = results[3]
self.extra_output = results[4]
Index: sherpa-4.16.1/sherpa/astro/tests/test_astro.py
===================================================================
--- sherpa-4.16.1.orig/sherpa/astro/tests/test_astro.py
+++ sherpa-4.16.1/sherpa/astro/tests/test_astro.py
@@ -206,7 +206,7 @@ def test_sourceandbg(parallel, run_threa
assert fit_results.numpoints == 1330
assert fit_results.dof == 1325
- assert covarerr[0] == approx(0.012097, rel=1e-3)
+ assert covarerr[0] == approx(0.012097, rel=1.05e-3)
assert covarerr[1] == approx(0, rel=1e-3)
assert covarerr[2] == approx(0.000280678, rel=1e-3)
assert covarerr[3] == approx(0.00990783, rel=1e-3)
Index: sherpa-4.16.1/docs/developer/index.rst
===================================================================
--- sherpa-4.16.1.orig/docs/developer/index.rst
+++ sherpa-4.16.1/docs/developer/index.rst
@@ -100,6 +100,17 @@ files and shows exactly which lines were
Run doctests locally
--------------------
+
+.. note::
+ The documentation tests are known to fail if NumPy 2.0 is installed
+ because the representation of NumPy types such as ``np.float64``
+ have changed, leading to errors like::
+
+ Expected:
+ 2.5264364698914e-06
+ Got:
+ np.float64(2.5264364698914e-06)
+
If `doctestplus <https://pypi.org/project/pytest-doctestplus/>` is installed
(and it probably is because it's part of
`sphinx-astropy <https://pypi.org/project/sphinx-astropy/>`,
Index: sherpa-4.16.1/docs/install.rst
===================================================================
--- sherpa-4.16.1.orig/docs/install.rst
+++ sherpa-4.16.1/docs/install.rst
@@ -34,17 +34,14 @@ Requirements
Sherpa has the following requirements:
* Python 3.9 to 3.11
-* NumPy (the exact lower limit has not been determined,
- 1.21.0 or later will work, earlier version may work)
+* NumPy (version 2.0 should work but there has been limited testing)
* Linux or OS-X (patches to add Windows support are welcome)
Sherpa can take advantage of the following Python packages
if installed:
* :term:`Astropy`: for reading and writing files in
- :term:`FITS` format. The minimum required version of astropy
- is version 1.3, although only versions 2 and higher are used in testing
- (version 3.2 is known to cause problems, but version 3.2.1 is okay).
+ :term:`FITS` format.
* :term:`matplotlib`: for visualisation of
one-dimensional data or models, one- or two- dimensional
error analysis, and the results of Monte-Carlo Markov Chain
