Sign Up
Log In
Log In
or
Sign Up
Places
All Projects
Status Monitor
Collapse sidebar
devel:languages:python:backports
python-sherpa
numpy2.patch
Overview
Repositories
Revisions
Requests
Users
Attributes
Meta
File numpy2.patch of Package 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
Locations
Projects
Search
Status Monitor
Help
OpenBuildService.org
Documentation
API Documentation
Code of Conduct
Contact
Support
@OBShq
Terms
openSUSE Build Service is sponsored by
The Open Build Service is an
openSUSE project
.
Sign Up
Log In
Places
Places
All Projects
Status Monitor
