forked from pool/python-sherpa
Matej Cepl
7ae14f2442
OBS-URL: https://build.opensuse.org/package/show/devel:languages:python:numeric/python-sherpa?expand=0&rev=44
298 lines
12 KiB
Diff
298 lines
12 KiB
Diff
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
|