Resort to median fitting of the sky to avoid an infinite loop in subt…

…ractSky. Triggered by EXTHD-193.

geminidr/core/parameters_preprocess.py

@@ -128,6 +128,7 @@ class subtractSkyConfig(config.Config):
     offset_sky = config.Field("Apply offset to sky frame to match science frame?", bool, False)
     sky = config.ListField("Sky frame to subtract", (str, AstroData), None, optional=True, single=True)
     save_sky = config.Field("Save sky frame to disk?", bool, False)
+    debug_threshold = config.Field("Convergence threshold when scaling", float, 0.001)
 
 
 class skyCorrectConfig(parameters_stack.stackSkyFramesConfig, subtractSkyConfig):

geminidr/core/primitives_preprocess.py

@@ -30,6 +30,8 @@
 from scipy.interpolate import interp1d
 from scipy.ndimage import binary_dilation
 
+from gempy.utils.errors import ConvergenceError
+
 from . import parameters_preprocess
 
 
@@ -1626,6 +1628,7 @@ def subtractSky(self, adinputs=None, **params):
         reset_sky = params["reset_sky"]
         scale = params["scale_sky"]
         zero = params["offset_sky"]
+        debug_threshold = params["debug_threshold"]
         if scale and zero:
             log.warning("Both the scale_sky and offset_sky parameters are set. "
                         "Setting offset_sky=False.")
@@ -1635,7 +1638,7 @@ def subtractSky(self, adinputs=None, **params):
         # in gt.measure_bg_from_image()
         sampling = 1 if adinputs[0].instrument() == 'GNIRS' else 10
         skyfunc = partial(gt.measure_bg_from_image, value_only=True,
-                          sampling=sampling)
+                          sampling=sampling, gaussfit=True)
 
         for ad, ad_sky in zip(*gt.make_lists(adinputs, params["sky"],
                                              force_ad=True)):
@@ -1652,8 +1655,26 @@ def subtractSky(self, adinputs=None, **params):
                             f"the science frame {ad.filename}")
                 if scale or zero:
                     # This actually does the sky subtraction as well
-                    factors = [gt.sky_factor(ext, ext_sky, skyfunc, multiplicative=scale)
-                               for ext, ext_sky in zip(ad, ad_sky)]
+                    try:
+                        factors = [gt.sky_factor(ext, ext_sky, skyfunc,
+                                                 multiplicative=scale,
+                                                 threshold=debug_threshold)
+                                   for ext, ext_sky in zip(ad, ad_sky)]
+                    except ConvergenceError as error:
+                        log.warning(f"The scaling of sky using a gaussian fit "
+                                    f"did not converge.  \n"
+                                    f"Using the median method instead.")
+                        skyfunc = partial(gt.measure_bg_from_image,
+                                          value_only=True,
+                                          sampling=sampling, gaussfit=False)
+                        try:
+                            factors = [gt.sky_factor(ext, ext_sky, skyfunc,
+                                                     multiplicative=scale)
+                                       for ext, ext_sky in zip(ad, ad_sky)]
+                        except ConvergenceError as error:
+                            log.error(f"Failed to scaled sky.")
+                            raise(error)
+
                     for ext_sky, factor in zip(ad_sky, factors):
                         log.fullinfo("Applying {} of {} to extension {}".
                                 format("scaling" if scale else "offset",

gempy/gemini/gemini_tools.py

@@ -30,6 +30,8 @@
 import astrodata
 from astrodata import Section
 
+from gempy.utils.errors import ConvergenceError
+
 ArrayInfo = namedtuple("ArrayInfo", "detector_shape origins array_shapes "
                                     "extensions")
 
@@ -1993,21 +1995,34 @@ def sky_factor(nd1, nd2, skyfunc, multiplicative=False, threshold=0.001):
     -------
     float : factor to apply to "sky" to match "science"
     """
+    log = logutils.get_logger(__name__)
+
     factor = 0
     if multiplicative:
-        current_sky = 1
+        current_sky = 1.
         # A subtlety here: deepcopy-ing an AD slice will create a full AD
         # object, and so skyfunc() will return a list instead of the single
         # float value we want. So make sure the copy is a single slice too
         if isinstance(nd1, astrodata.AstroData) and nd1.is_single:
             ndcopy = deepcopy(nd1)[0]
         else:
             ndcopy = deepcopy(nd1)
-        while abs(current_sky) > threshold:
+        iter = 1
+        max_iter = 100   # normally converges in < 10 iterations
+        while abs(current_sky) > threshold and iter <= max_iter:
             f = skyfunc(ndcopy) / skyfunc(nd2)
             ndcopy.subtract(nd2.multiply(f))
             current_sky *= f
             factor += current_sky
+            iter += 1
+        #print('iter upon exit: ', iter)
+        if iter > max_iter:
+            log.warning(f"Failed to converge.\n"
+                        f"   Reached: {abs(current_sky)} while threshold = {threshold}\n"
+                        f"   Final factor = {factor}")
+            nd2.divide(current_sky)  # reset to original value
+            raise(ConvergenceError)
+
         nd1.subtract(nd2.multiply(factor / current_sky))
         nd2.divide(factor)  # reset to original value
     else:

gempy/utils/errors.py

@@ -0,0 +1,2 @@
+class ConvergenceError(RuntimeError):
+    pass