Merge sherpa#2108 (DougBurke) - Asymmetric settings

Asymmetric settings

sherpa/astro/ui/utils.py

@@ -1790,9 +1790,13 @@ def unpack_data(self, filename, *args, **kwargs):
                     # If this errors out then so be it
                     return self.unpack_ascii(filename, *args, **kwargs)
 
-    def load_ascii_with_errors(self, id, filename=None, colkeys=None, sep=' ',
-                               comment='#', func=np.average,
-                               delta=False) -> None:
+    def load_ascii_with_errors(self, id, filename=None,
+                               colkeys: Optional[Sequence[str]] = None,
+                               sep: str = ' ',
+                               comment: str = '#',
+                               func: Callable = np.average,
+                               delta: bool = False
+                               ) -> None:
         """Load an ASCII file with asymmetric errors as a data set.
 
         Create a dataset with asymmetric error bars which can be used
@@ -14406,7 +14410,9 @@ def plot_bkg_fit_delchi(self,
 
     def resample_data(self,
                       id: Optional[IdType] = None,
-                      niter=1000, seed=None):
+                      niter: int = 1000,
+                      seed: Optional[int] = None
+                      ) -> dict[str, np.ndarray]:
         """Resample data with asymmetric error bars.
 
         The function performs a parametric bootstrap assuming a skewed
@@ -14418,6 +14424,12 @@ def resample_data(self,
         each realization, and displays the average and standard
         deviation for each parameter.
 
+        .. versionchanged:: 4.17.0
+           The resampling now uses the chosen statistic and optimizer
+           (set with set_stat and set_method). Previously the
+           least-squares statistic and Levenberg-Marquardt method were
+           always used.
+
         .. versionchanged:: 4.16.0
            The random number generation is now controlled by the
            `set_rng` routine. The seed argument is now deprecated.
@@ -14457,6 +14469,8 @@ def resample_data(self,
         Account for of asymmetric errors when calculating parameter
         uncertainties:
 
+        >>> set_stat("leastsq")
+        >>> set_method("levmar")
         >>> load_ascii_with_errors(1, 'test.dat')
         >>> set_model(polynom1d.p0)
         >>> thaw(p0.c1)
@@ -14500,8 +14514,12 @@ def resample_data(self,
         """
         data = self.get_data(id)
         model = self.get_model(id)
+        stat = self.get_stat()
+        method = self.get_method()
+
         resampledata = sherpa.sim.ReSampleData(data, model)
         return resampledata(niter=niter, seed=seed,
+                            stat=stat, method=method,
                             rng=self.get_rng())
 
     def sample_photon_flux(self, lo=None, hi=None,

sherpa/sim/__init__.py

@@ -1,5 +1,5 @@
 #
-#  Copyright (C) 2011, 2015, 2016, 2018, 2019, 2020, 2021, 2023
+#  Copyright (C) 2011, 2015, 2016, 2018 - 2021, 2023, 2024
 #  Smithsonian Astrophysical Observatory
 #
 #
@@ -193,8 +193,13 @@
 """
 
 import logging
+from typing import Optional
 
-import numpy
+import numpy as np
+
+from sherpa.data import Data1D, Data1DAsymmetricErrs
+from sherpa.fit import Fit
+from sherpa.optmethods import LevMar, OptMethod
 
 # Although all this module needs is the following import
 #   from sherpa.sim.mh import LimitError, MetropolisMH, MH, Sampler, Walk
@@ -204,19 +209,16 @@
 from sherpa.sim.simulate import *
 from sherpa.sim.sample import *
 from sherpa.sim.mh import *
-from sherpa.stats import Cash, CStat, WStat, LeastSq
+
+from sherpa.stats import Cash, CStat, WStat, LeastSq, Stat
 from sherpa.utils import NoNewAttributesAfterInit, get_keyword_defaults, \
     sao_fcmp
+from sherpa.utils.logging import SherpaVerbosity
 from sherpa.utils import random
 
-from sherpa.fit import Fit
-from sherpa.data import Data1D, Data1DAsymmetricErrs
-from sherpa.optmethods import LevMar
-
 info = logging.getLogger("sherpa").info
-_log = logging.getLogger("sherpa")
 
-_tol = numpy.finfo(float).eps
+_tol = np.finfo(float).eps
 
 string_types = (str, )
 
@@ -692,7 +694,6 @@ def get_draws(self, fit, sigma, niter=1000, cache=True, rng=None):
             raise ValueError("Fit statistic must be cash, cstat or "
                              f"wstat, not {fit.stat.name}")
 
-        _level = _log.getEffectiveLevel()
         mu = fit.model.thawedpars
         dof = len(mu)
 
@@ -714,7 +715,7 @@ def get_draws(self, fit, sigma, niter=1000, cache=True, rng=None):
         sampler_kwargs = self._sampler_opt.copy()
         sampler_kwargs['priors'] = priors
 
-        oldthawedpars = numpy.array(mu)
+        oldthawedpars = np.array(mu)
         thawedparmins = fit.model.thawedparhardmins
         thawedparmaxes = fit.model.thawedparhardmaxes
 
@@ -726,28 +727,19 @@ def calc_stat(proposed_params):
             if -1 in mins or -1 in maxes:
                 raise LimitError('Sherpa parameter hard limit exception')
 
-            level = _log.getEffectiveLevel()
-
-            try:
-                # ignore warning from Sherpa about hard limits
-                _log.setLevel(logging.CRITICAL)
-
-                # soft limits are ignored, hard limits rejected.
-                # proposed values beyond hard limit default to limit.
-                fit.model.thawedpars = proposed_params
-
-                # Calculate statistic on proposal, use likelihood
-                proposed_stat = -0.5 * fit.calc_stat()
+            with SherpaVerbosity(logging.CRITICAL):
+                try:
+                    # soft limits are ignored, hard limits rejected.
+                    # proposed values beyond hard limit default to limit.
+                    fit.model.thawedpars = proposed_params
 
-                # _log.setLevel(level)
+                    # Calculate statistic on proposal, use likelihood
+                    proposed_stat = -0.5 * fit.calc_stat()
 
-            except:
-                # set the model back to original state on exception
-                fit.model.thawedpars = oldthawedpars
-                raise
-            finally:
-                # set the logger back to previous level
-                _log.setLevel(level)
+                except:
+                    # set the model back to original state on exception
+                    fit.model.thawedpars = oldthawedpars
+                    raise
 
             return proposed_stat
 
@@ -762,9 +754,6 @@ def calc_stat(proposed_params):
             # set the model back to original state
             fit.model.thawedpars = oldthawedpars
 
-            # set the logger back to previous level
-            _log.setLevel(_level)
-
         # Change to Sherpa statistic convention
         stats = -2.0 * stats
 
@@ -845,14 +834,31 @@ def __init__(self, data, model):
 
         self.data = data
         self.model = model
-        NoNewAttributesAfterInit.__init__(self)
-
-    def __call__(self, niter=1000, seed=None, rng=None):
-        return self.call(niter, seed, rng=rng)
-
-    def call(self, niter, seed=None, rng=None):
+        super().__init__()
+
+    def __call__(self, niter=1000, seed=None, rng=None,
+                 *,
+                 stat: Optional[Stat] = None,
+                 method: Optional[OptMethod] = None
+                 ) -> dict[str, np.ndarray]:
+        return self.call(niter, seed=seed, rng=rng, stat=stat,
+                         method=method)
+
+    def call(self, niter, seed=None, rng=None,
+             *,
+             # Mark these as keyword-only as they are additions to the
+             # interface and the interface is complicated enough it
+             # is worth marking them keyword only.
+             #
+             stat: Optional[Stat] = None,
+             method: Optional[OptMethod] = None
+             ) -> dict[str, np.ndarray]:
         """Resample the data and fit the model to each iteration.
 
+        .. versionadded:: 4.17.0
+           The stat and method parameter were added. These are
+           keyword-only arguments.
+
         .. versionadded:: 4.16.0
            The rng parameter was added.
 
@@ -872,6 +878,10 @@ def call(self, niter, seed=None, rng=None):
            Determines how random numbers are created. If set to None then
            the routines from `numpy.random` are used, and so can be
            controlled by calling `numpy.random.seed`.
+        stat : Stat or None, optional
+           If None then the `LeastSq` statistic is used.
+        method: OptMethod or None, optional
+           If None then the `LevMar` method is used.
 
         Returns
         -------
@@ -889,6 +899,9 @@ def call(self, niter, seed=None, rng=None):
 
         """
 
+        chosen_stat = LeastSq() if stat is None else stat
+        chosen_meth = LevMar() if method is None else method
+
         # Each fit is reset to this set of values as the starting point
         orig_pars = self.model.thawedpars
 
@@ -900,7 +913,7 @@ def call(self, niter, seed=None, rng=None):
 
             name = par.fullname
             pars_index.append(name)
-            pars[name] = numpy.zeros(niter)
+            pars[name] = np.zeros(niter)
 
         data = self.data
         x = data.x
@@ -915,13 +928,13 @@ def call(self, niter, seed=None, rng=None):
         ny = len(y)
 
         # TODO: we do not properly handle Data1DInt here
-        fake_data = Data1D('tmp', x, numpy.zeros(ny))
+        fake_data = Data1D('tmp', x, np.zeros(ny))
 
         if rng is None:
-            numpy.random.seed(seed)
+            np.random.seed(seed)
 
-        ry_all = numpy.zeros((niter, ny), dtype=y_l.dtype)
-        stats = numpy.zeros(niter)
+        ry_all = np.zeros((niter, ny), dtype=y_l.dtype)
+        stats = np.zeros(niter)
         for j in range(niter):
             ry = ry_all[j]
             for i in range(ny):
@@ -972,7 +985,8 @@ def call(self, niter, seed=None, rng=None):
             # start the fit (by making sure we always reset after a fit).
             #
             fake_data.y = ry
-            fit = Fit(fake_data, self.model, LeastSq(), LevMar())
+            fit = Fit(data=fake_data, model=self.model,
+                      stat=chosen_stat, method=chosen_meth)
             try:
                 fit_result = fit.fit()
             finally:
@@ -984,8 +998,8 @@ def call(self, niter, seed=None, rng=None):
 
         result = {'samples': ry_all, 'statistic': stats}
         for name in pars_index:
-            avg = numpy.average(pars[name])
-            std = numpy.std(pars[name])
+            avg = np.average(pars[name])
+            std = np.std(pars[name])
             info('%s : avg = %s , std = %s', name, avg, std)
             result[name] = pars[name]
 

sherpa/ui/utils.py

@@ -97,6 +97,8 @@ def _check_str_type(arg: str, argname: str) -> None:
     raise ArgumentTypeErr('badarg', argname, "a string")
 
 
+# With Python 3.10 these can return TypeGuard annotations.
+#
 def _is_integer(val):
     return isinstance(val, (int, np.integer))
 
@@ -2344,7 +2346,10 @@ def get_method_name(self) -> str:
 
     # DOC-TODO: remove the list of supported methods once the
     # relevant documentation has been updated.
-    def set_method(self, meth: str) -> None:
+    #
+    def set_method(self,
+                   meth: Union[OptMethod, str]
+                   ) -> None:
         """Set the optimization method.
 
         The primary task of Sherpa is to fit a model M(p) to a set of
@@ -2417,7 +2422,7 @@ def set_method(self, meth: str) -> None:
         if _is_str(meth):
             method = self._get_method_by_name(meth)
         else:
-            _check_type(meth, sherpa.optmethods.OptMethod, 'meth',
+            _check_type(meth, OptMethod, 'meth',
                         'a method name or object')
             method = meth
 
@@ -2442,7 +2447,9 @@ def _check_method_opt(self, optname: str) -> None:
         if optname not in self._current_method.config:
             raise ArgumentErr('badopt', optname, self.get_method_name())
 
-    def get_method_opt(self, optname=None):
+    def get_method_opt(self,
+                       optname: Optional[str] = None
+                       ) -> Any:
         """Return one or all of the options for the current optimization
         method.
 
@@ -2486,7 +2493,7 @@ def get_method_opt(self, optname=None):
         self._check_method_opt(optname)
         return self._current_method.config[optname]
 
-    def set_method_opt(self, optname, val) -> None:
+    def set_method_opt(self, optname: str, val: Any) -> None:
         """Set an option for the current optimization method.
 
         This is a helper function since the optimization options can also
@@ -2546,7 +2553,9 @@ def get_iter_method_name(self) -> str:
         """
         return self._current_itermethod['name']
 
-    def get_iter_method_opt(self, optname=None):
+    def get_iter_method_opt(self,
+                            optname: Optional[str] = None
+                            ) -> Any:
         """Return one or all options for the iterative-fitting scheme.
 
         The options available for the iterative-fitting methods are
@@ -2725,7 +2734,7 @@ def set_iter_method(self, meth: str) -> None:
 
         self._current_itermethod = self._itermethods[meth]
 
-    def set_iter_method_opt(self, optname, val) -> None:
+    def set_iter_method_opt(self, optname: str, val: Any) -> None:
         """Set an option for the iterative-fitting scheme.
 
         Parameters