File numpy2.patch of Package python-sherpa

Overview Repositories Revisions Requests Users Attributes Meta

File numpy2.patch of Package python-sherpa (Revision 1bb503837b5659ab07ce1bc11f4437d5)

Currently displaying revision 1bb503837b5659ab07ce1bc11f4437d5 , Show latest

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

Places

File numpy2.patch of Package python-sherpa (Revision 1bb503837b5659ab07ce1bc11f4437d5)

Places