Add a couple fixes related to low snr objects

piff/gsobject_model.py

@@ -187,9 +187,9 @@ def _resid(self, params, star, convert_func, draw_method):
         image, weight, image_pos = star.data.getImage()
         flux, du, dv, scale, g1, g2 = params
 
-        # Make sure the shear is sane.
+        # Make sure the size and shear are sane.
         g = g1 + 1j * g2
-        if np.abs(g) >= 1.:
+        if scale <= 0. or np.abs(g) >= 1.:
             # Return "infinity"
             return np.ones_like(image.array.ravel(), dtype=np.float64) * 1.e300
 

piff/input.py

@@ -887,6 +887,9 @@ def _makeStarsFromImage(image_kwargs, cat_kwargs, wcs, chipnum,
             pos = galsim.PositionD(x,y)
             data = StarData(stamp, pos, weight=wt_stamp, pointing=pointing,
                             properties=props, property_types=prop_types)
+            # Attach an attribute so we can tell later whether the weight is real or
+            # if it's just the default w=1.0.
+            data._piff_default_weight = getattr(wt, '_piff_default_weight', False)
             star = Star(data, None)
 
             objects.append(star)
@@ -955,13 +958,12 @@ def readImage(image_file_name, image_hdu, weight_file_name, weight_hdu,
         :param sky_file_name:   A file to use for a sky background to subtract from the image
                                 (if any).
         :param sky_hdu:         The hdu to use in the sky_file_name (if any).
-        :param noise:           Either a float constant noise value to use in lieu of a weight
-                                map or a str keyword to use to read a value from FITS header.
         :param sky:             If this is 'median', then treat the median as the sky level.
                                 Otherwise, the sky level is set to the value passed here,
                                 or to the value in the fits header associated to the
                                 keyword that's passed here.
-        :param noise:           A constant noise value to use in lieu of a weight map.
+        :param noise:           Either a float constant noise value to use in lieu of a weight
+                                map or a str keyword to use to read a value from FITS header.
         :param logger:          A logger object for logging debug info.
 
         :returns: image, weight
@@ -989,6 +991,7 @@ def readImage(image_file_name, image_hdu, weight_file_name, weight_hdu,
                 logger.warning("Warning: weight map has invalid negative-valued pixels. "+
                                "Taking them to be 0.0")
                 weight.array[weight.array < 0] = 0.
+            weight._piff_default_weight = False
         elif noise is not None:
             try:
                 noise = float(noise)
@@ -998,9 +1001,11 @@ def readImage(image_file_name, image_hdu, weight_file_name, weight_hdu,
                 noise = float(image.header[noise])
             logger.debug("Making uniform weight image based on noise variance = %f", noise)
             weight = galsim.ImageF(image.bounds, init_value=1./noise)
+            weight._piff_default_weight = False
         else:
             logger.debug("Making trivial (wt==1) weight image")
             weight = galsim.ImageF(image.bounds, init_value=1)
+            weight._piff_default_weight = True
 
         # Make sure to do this before anything that could modify the input image.
         if badpix_zeros:

piff/select.py

@@ -331,6 +331,18 @@ def rejectStars(self, stars, logger=None):
 
             good_stars.append(star)
 
+        zero_snr_count = sum(s['raw_snr'] == 0 for s in good_stars)
+        if zero_snr_count > 0 and self.min_snr is None:
+            logger.warning("Found %d candidate stars with S/N <= 0. Consider setting "
+                           "select.min_snr to remove noise-dominated objects.",
+                           zero_snr_count)
+            default_weight = any(getattr(s.data, '_piff_default_weight', False)
+                                 for s in good_stars if s['raw_snr']== 0)
+            if default_weight:
+                logger.warning("Note -- no variance information (either noise value or weight map) "
+                               "was provided, so these estimates are not likely to be accurate. "
+                               "Consider setting input.noise or providing a weight image.")
+
         # Calculate the hsm size for each star and throw out extreme outliers.
         if self.hsm_size_reject != 0:
             sigma = [star.hsm[3] for star in good_stars]

tests/test_gsobject_model.py

@@ -1127,6 +1127,16 @@ def test_fail():
         model3.fit(star3)
     print('7')
 
+    # Non-positive scale values are invalid optimizer probes and should be rejected before
+    # trying to draw them with GalSim.
+    model = piff.Gaussian()
+    star = model.initialize(star1)
+    params = model._get_params(star)
+    params[3] = -0.1
+    resid = model._resid(params, star, convert_func=None, draw_method=None)
+    np.testing.assert_array_equal(resid, np.ones_like(resid) * 1.e300)
+    print('7b')
+
     # reflux is harder to make fail.  Rather than try something even more contrived,
     # just mock np.linalg.solve to simulate the case that AtA ends up singular.
     with mock.patch('numpy.linalg.solve', side_effect=np.linalg.LinAlgError) as raising_solve:

tests/test_input.py

@@ -1805,6 +1805,26 @@ def test_stars():
     snr0 = piff.util.calculateSNR(star0.data.image, star0.data.orig_weight)
     assert snr0 == 0.
 
+    # raw_snr=0 stars are allowed by default, but warn that they are noise-dominated.
+    image0 = galsim.Image(8, 8, init_value=0., scale=1.)
+    image1 = galsim.Image(8, 8, init_value=10., scale=1.)
+    weight0 = galsim.Image(8, 8, init_value=1., scale=1.)
+    weight1 = galsim.Image(8, 8, init_value=1., scale=1.)
+    zero_snr_star = piff.Star(piff.StarData(image0, image0.true_center, weight0), None)
+    positive_snr_star = piff.Star(piff.StarData(image1, image1.true_center, weight1), None)
+    select = piff.FlagSelect({}, logger=logger)
+    with CaptureLog(1) as cl:
+        stars = select.rejectStars([zero_snr_star, positive_snr_star], logger=cl.logger)
+    assert len(stars) == 2
+    assert stars[0]['raw_snr'] == 0.
+    assert "Found 1 candidate stars with S/N <= 0" in cl.output
+    assert "no variance information" not in cl.output
+
+    zero_snr_star.data._piff_default_weight = True
+    with CaptureLog(1) as cl:
+        select.rejectStars([zero_snr_star, positive_snr_star], logger=cl.logger)
+    assert "no variance information" in cl.output
+
 
 @timer
 def test_select_min_sep():