Repository: transientskp/pyse

Issue:
Ground truth data from testSingleSourceExtraction unit test was based on older algorithms
The [TestSimpleImageSourceFind.testSingleSourceExtraction](https://github.com/transientskp/pyse/blob/8bd77c9ae181ade28a0f78ceeed02fe25e90ccd6/test/test_image.py#L293C9-L293C35) unit test compares source measurements with [ground truth data](https://github.com/transientskp/pyse/blob/8bd77c9ae181ade28a0f78ceeed02fe25e90ccd6/test/test_image.py#L301):

```
        known_result = (
            136.89603241069054, 14.022184792492785,  # RA, DEC
            5.341819139061954e-4, 1.3428186757078464e-3,  # Err, Err
            7.226590529214518e-4, 1.0918184742211533e-4,  # Peak flux, err
            6.067963179204716e-4, 1.7037685531724465e-4,
            # Integrated flux, err
            6.192259965962862, 25.516190123153514,
            # Significance level, Beam semimajor-axis width (arcsec)
            10.718798843620489, 178.62899212789304,
            # Beam semiminor-axis width (arcsec), Beam parallactic angle
            ew_sys_err, ns_sys_err,
            5.181697175052841,  # error_radius
            1,  # fit_type
            0.59184643302,  # chisq
            0.67199741142,  # reduced chisq
        )
```

However, it is unclear whether this _ground truth_ corresponds with the most accurate measurement possible.
This would be different if an artificial source were inserted in the image, i.e. what is measured should be close to what is inserted, but that is not the case here, since this is not an artificial image.

How to proceed?
Currently, the test is passed only when a [minimal resemblance is imposed (self.assertAlmostEqual(r[i], known_result[i], places=0))](https://github.com/transientskp/pyse/blob/8bd77c9ae181ade28a0f78ceeed02fe25e90ccd6/test/test_image.py#L329C13-L329C68).

A long time ago, with the older algorithms, a [much higher resemblance could be imposed (self.assertAlmostEqual(r[i], known_result[i], places=5))](https://github.com/transientskp/pyse/blob/a8310c2967bf597d265376e7e5bd8ca598bb34be/test/test_image.py#L308C13-L308C68), but that does not mean that the older algorithms were more accurate, on the contrary.

Not sure what to do here. Perhaps replacing with an artificial image would be best.
But that would boil down to a complete replacement of this unit test.

Or investigate how the old algorithms differ from the current ones, i.e. explain how the changes in the source measurement emerged. Perhaps the explanation from [commit 2ee8447](https://github.com/transientskp/pyse/commit/2ee84472343eed2d1b367ca2a4bd85d32fa8543c) is sufficient:
```
I made some of the unit tests less tight. 
They were based on noise levels from the original version of stats.sigma_clip, 
with very little margin, as if the original stats.sigma_clip was perfect.
This original version, however, failed for images with a very high source density. 
The present version gives slightly lower noise levels in these circumstances 
and closer to the ground truth.
```



Discussion:

Patch:
diff --git a/sourcefinder/extract.py b/sourcefinder/extract.py
index ae45c69..af4f924 100644
--- a/sourcefinder/extract.py
+++ b/sourcefinder/extract.py
@@ -210,7 +210,7 @@ class Island(object):
         try:
             measurement, gauss_island, gauss_residual = \
                 source_profile_and_errors(self.data, self.threshold(),
-                                          self.noise(), self.beam,
+                                          self.rms, self.noise(), self.beam,
                                           fudge_max_pix_factor,
                                           max_pix_variance_factor, beamsize,
                                           correlation_lengths, fixed=fixed)
@@ -681,7 +681,7 @@ class ParamSet(MutableMapping):
         return self
 
 
-def source_profile_and_errors(data, threshold, noise,
+def source_profile_and_errors(data, threshold, rms, noise,
                               beam, fudge_max_pix_factor, max_pix_variance_factor,
                               beamsize, correlation_lengths, fixed=None):
     """Return a number of measurable properties with errorbars
@@ -706,7 +706,11 @@ def source_profile_and_errors(data, threshold, noise,
         threshold (float): Threshold used for selecting pixels for the
             source (ie, building an island)
 
-        noise (float): Noise level in data
+        rms (np.ndarray or np.ma.MaskedArray): noise levels at pixel positions
+            corresponding to the data array, determined from the interpolated
+            grid of standard deviations of the background pixels.
+
+        noise (float): rms (noise level) at the maximum pixel position
 
         beam (tuple): beam parameters (semimaj,semimin,theta)
 
@@ -847,7 +851,7 @@ def source_profile_and_errors(data, threshold, noise,
         gauss_resid_filled = gauss_resid_masked
 
     param.chisq, param.reduced_chisq = fitting.goodness_of_fit(
-        gauss_resid_masked, noise, correlation_lengths)
+        gauss_resid_masked, rms, correlation_lengths)
 
     return param, gauss_island_filled, gauss_resid_filled
 
diff --git a/sourcefinder/fitting.py b/sourcefinder/fitting.py
index c63081c..75ea2b3 100644
--- a/sourcefinder/fitting.py
+++ b/sourcefinder/fitting.py
@@ -143,12 +143,12 @@ def moments(data, fudge_max_pix_factor, beamsize, threshold=0):
 
                 # Set limits for the root finder similar to the bounds for
                 # Gaussian fits.
-                if peak > 0:
-                    low_bound = 0.5 * peak
-                    upp_bound = 1.5 * peak
+                if basevalue > 0:
+                    low_bound = 0.5 * basevalue
+                    upp_bound = 1.5 * basevalue
                 else:
-                    low_bound = 1.5 * peak
-                    upp_bound = 0.5 * peak
+                    low_bound = 1.5 * basevalue
+                    upp_bound = 0.5 * basevalue
 
                 # The number of iterations used for the root finder is also
                 # returned, but not used here.
@@ -443,13 +443,12 @@ def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
                 # position.
                 for i in (posx==rounded_barycenter[0]).nonzero()[0]:
                     if posy[i] == rounded_barycenter[1]:
-                        baseindex = i
-                        basevalue = source_island[baseindex]
+                        basevalue = source_island[i]
                         break
             else:
                 # The rounded barycenter position is not in source_island, so we
                 # revert to the maximum pixel position, which is always included.
-                basepos = maxpos[0], maxpos[1]
+                basepos = maxpos[0] - chunkpos[0], maxpos[1] - chunkpos[1]
                 # In this case we do not need to figure out the source_island
                 # index corresponding to the maximum pixel position, because
                 # maxi, the maximum pixel value has already been provided, as an
@@ -472,12 +471,12 @@ def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
 
                 # Set limits for the root finder similar to the bounds for
                 # Gaussian fits.
-                if peak > 0:
-                    low_bound = 0.5 * peak
-                    upp_bound = 1.5 * peak
+                if basevalue > 0:
+                    low_bound = 0.5 * basevalue
+                    upp_bound = 1.5 * basevalue
                 else:
-                    low_bound = 1.5 * peak
-                    upp_bound = 0.5 * peak
+                    low_bound = 1.5 * basevalue
+                    upp_bound = 0.5 * basevalue
 
                 # The number of iterations used for the root finder is also
                 # returned, but not used here.
@@ -1017,10 +1016,6 @@ def goodness_of_fit(masked_residuals, noise, correlation_lengths):
     of datapoints less than 1! As a result, it doesn't really make sense to try
     and apply the 'degrees-of-freedom' correction, as this would likely result
     in a negative ``reduced_chisq`` value.
-    
-    Finally, note that when called from extract.source_profile_and_errors, the
-    noise-estimate at the peak-pixel is supplied, so will typically over-estimate
-    the noise and hence under-estimate the chi-squared values.
     """
     gauss_resid_normed = (masked_residuals / noise).compressed()
     chisq = np.sum(gauss_resid_normed * gauss_resid_normed)
diff --git a/sourcefinder/image.py b/sourcefinder/image.py
index 31c19fb..e9faf18 100644
--- a/sourcefinder/image.py
+++ b/sourcefinder/image.py
@@ -745,7 +745,8 @@ class ImageData(object):
 
         try:
             measurement, _, _ = extract.source_profile_and_errors(
-                fitme, threshold_at_pixel, self.rmsmap[int(x), int(y)],
+                fitme, threshold_at_pixel, self.rmsmap[chunk],
+                self.rmsmap[int(x), int(y)],
                 self.beam, self.fudge_max_pix_factor,
                 self.max_pix_variance_factor,
                 self.beamsize, self.correlation_lengths, fixed=fixed)
@@ -1265,7 +1266,9 @@ class ImageData(object):
                 chunk = slices[label - 1]
 
                 param = extract.ParamSet()
-                param.sig = maxis[count] / self.rmsmap.data[tuple(maxposs[count])]
+                param.sig = sig[count]
+                param.chisq = chisq[count]
+                param.reduced_chisq = reduced_chisq[count]
 
                 param["peak"] = Uncertain(moments_of_sources[count, 0, 0],
                                           moments_of_sources[count, 1, 0])


Test patch:
diff --git a/test/test_L15_12h_const.py b/test/test_L15_12h_const.py
index 5439495..2b3664f 100644
--- a/test/test_L15_12h_const.py
+++ b/test/test_L15_12h_const.py
@@ -79,10 +79,10 @@ class L15_12hConstCor(unittest.TestCase):
         self.assertEqual(len(self.results), 5)
 
     @requires_data(corrected_fits)
-    def testFluxes(self):
-        # All sources in this image are supposed to have the same flux.
-        # But they don't, because the simulation is broken, so this test
-        # checks they fall in a vaguely plausible range.
+    def testBrightnesses(self):
+        # All sources in this image are supposed to have the same peak spectral
+        # brightness. But they don't, because the simulation is broken, so this
+        # test checks they fall in a vaguely plausible range.
         for mysource in self.results:
             self.assertTrue(mysource.peak.value > 0.35)
             self.assertTrue(mysource.peak.value < 0.60)
diff --git a/test/test_gaussian.py b/test/test_gaussian.py
index 3b111b8..a8c8d85 100644
--- a/test/test_gaussian.py
+++ b/test/test_gaussian.py
@@ -3,7 +3,7 @@ Tests for elliptical Gaussian fitting code in the TKP pipeline.
 """
 import unittest
 
-import numpy
+import numpy as np
 
 from sourcefinder.extract import source_profile_and_errors
 from sourcefinder.fitting import moments, fitgaussian, FIT_PARAMS
@@ -21,14 +21,14 @@ class SimpleGaussTest(unittest.TestCase):
     """Generic, easy-to-fit elliptical Gaussian"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = 10
         self.x = 250
         self.y = 250
         self.maj = 40
         self.min = 20
         self.theta = 0
-        self.mygauss = numpy.ma.array(
+        self.mygauss = np.ma.array(
             gaussian(self.height, self.x, self.y, self.maj, self.min,
                      self.theta)(Xin, Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -69,14 +69,14 @@ class NegativeGaussTest(SimpleGaussTest):
     """Negative Gaussian"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = -10
         self.x = 250
         self.y = 250
         self.maj = 40
         self.min = 20
         self.theta = 0
-        self.mygauss = numpy.ma.array(
+        self.mygauss = np.ma.array(
             gaussian(self.height, self.x, self.y, self.maj, self.min,
                      self.theta)(Xin, Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -92,14 +92,14 @@ class CircularGaussTest(SimpleGaussTest):
     """Circular Gaussian: it makes no sense to measure a rotation angle"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = 10
         self.x = 250
         self.y = 250
         self.maj = 40
         self.min = 40
         self.theta = 0
-        self.mygauss = numpy.ma.array(
+        self.mygauss = np.ma.array(
             gaussian(self.height, self.x, self.y, self.maj, self.min,
                      self.theta)(Xin, Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -118,14 +118,14 @@ class NarrowGaussTest(SimpleGaussTest):
     """Only 1 pixel wide"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = 10
         self.x = 250
         self.y = 250
         self.maj = 40
         self.min = 1
         self.theta = 0
-        self.mygauss = numpy.ma.array(gaussian(
+        self.mygauss = np.ma.array(gaussian(
             self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
                                                                          Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -138,14 +138,14 @@ class RotatedGaussTest(SimpleGaussTest):
     """Rotated by an angle < pi/2"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = 10
         self.x = 250
         self.y = 250
         self.maj = 40
         self.min = 20
-        self.theta = numpy.pi / 4
-        self.mygauss = numpy.ma.array(gaussian(
+        self.theta = np.pi / 4
+        self.mygauss = np.ma.array(gaussian(
             self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
                                                                          Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -161,14 +161,14 @@ class RotatedGaussTest2(SimpleGaussTest):
     """Rotated by an angle > pi/2; theta becomes negative"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = 10
         self.x = 250
         self.y = 250
         self.maj = 40
         self.min = 20
-        self.theta = 3 * numpy.pi / 4
-        self.mygauss = numpy.ma.array(gaussian(
+        self.theta = 3 * np.pi / 4
+        self.mygauss = np.ma.array(gaussian(
             self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
                                                                          Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -177,10 +177,10 @@ class RotatedGaussTest2(SimpleGaussTest):
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testMomentAngle(self):
-        self.assertAlmostEqual(self.moments["theta"], self.theta - numpy.pi)
+        self.assertAlmostEqual(self.moments["theta"], self.theta - np.pi)
 
     def testFitAngle(self):
-        self.assertAlmostEqual(self.fit["theta"], self.theta - numpy.pi)
+        self.assertAlmostEqual(self.fit["theta"], self.theta - np.pi)
 
 
 class AxesSwapGaussTest(SimpleGaussTest):
@@ -188,14 +188,14 @@ class AxesSwapGaussTest(SimpleGaussTest):
     change the angle"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.height = 10
         self.x = 250
         self.y = 250
         self.maj = 20
         self.min = 40
         self.theta = 0
-        self.mygauss = numpy.ma.array(gaussian(
+        self.mygauss = np.ma.array(gaussian(
             self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
                                                                          Yin))
         self.beamsize = calculate_beamsize(self.maj, self.min)
@@ -209,8 +209,8 @@ class AxesSwapGaussTest(SimpleGaussTest):
         # Presumably there's some numerical quirk causing different,
         # but equivalent, convergence in the optimization.
         if theta < 0:
-            theta = theta + numpy.pi
-        self.assertAlmostEqual(theta, numpy.pi / 2)
+            theta = theta + np.pi
+        self.assertAlmostEqual(theta, np.pi / 2)
 
     def testFitAngle(self):
         theta = self.fit["theta"]
@@ -218,8 +218,8 @@ class AxesSwapGaussTest(SimpleGaussTest):
         # Presumably there's some numerical quirk causing different,
         # but equivalent, convergence in the optimization.
         if theta < 0:
-            theta = theta + numpy.pi
-        self.assertAlmostEqual(theta, numpy.pi / 2)
+            theta = theta + np.pi
+        self.assertAlmostEqual(theta, np.pi / 2)
 
     def testMomentSize(self):
         self.assertAlmostEqual(self.moments["semiminor"], self.maj, 5)
@@ -235,8 +235,8 @@ class RandomGaussTest(unittest.TestCase):
     measure moments, though -- things should be fairly evenly distributed."""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
-        self.mygauss = numpy.random.random(Xin.shape)
+        Xin, Yin = np.indices((500, 500))
+        self.mygauss = np.random.random(Xin.shape)
         self.beamsize = calculate_beamsize(beam[0], beam[1])
         self.fudge_max_pix_factor = fudge_max_pix(beam[0], beam[1], beam[2])
 
@@ -251,14 +251,14 @@ class NoisyGaussTest(unittest.TestCase):
     """Test calculation of chi-sq and fitting in presence of (artificial) noise"""
 
     def setUp(self):
-        Xin, Yin = numpy.indices((500, 500))
+        Xin, Yin = np.indices((500, 500))
         self.peak = 10
         self.xbar = 250
         self.ybar = 250
         self.semimajor = 40
         self.semiminor = 20
         self.theta = 0
-        self.mygauss = numpy.ma.array(
+        self.mygauss = np.ma.array(
             gaussian(self.peak, self.xbar, self.ybar,
                      self.semimajor, self.semiminor, self.theta)(Xin, Yin))
         self.beamsize = calculate_beamsize(self.semimajor, self.semiminor)
@@ -287,7 +287,7 @@ class NoisyGaussTest(unittest.TestCase):
 
     def test_noisy_background(self):
         # Use a fixed random state seed, so unit-test is reproducible:
-        rstate = numpy.random.RandomState(42)
+        rstate = np.random.RandomState(42)
         pixel_noise = 0.5
         self.mygauss += rstate.normal(scale=pixel_noise,
                                       size=len(self.mygauss.ravel())).reshape(
@@ -310,11 +310,15 @@ class NoisyGaussTest(unittest.TestCase):
         self.fit_w_errs, _, _ = source_profile_and_errors(
             data=self.mygauss,
             threshold=0.,
+            rms=np.ones(self.mygauss.shape) * pixel_noise,
             noise=pixel_noise,
             beam=(self.semimajor, self.semiminor, self.theta),
-            fudge_max_pix_factor=fudge_max_pix(self.semimajor, self.semiminor, self.theta),
-            max_pix_variance_factor=maximum_pixel_method_variance(self.semimajor, self.semiminor, self.theta),
-            correlation_lengths=calculate_correlation_lengths(self.semimajor, self.semiminor),
+            fudge_max_pix_factor=fudge_max_pix(self.semimajor, self.semiminor,
+                                               self.theta),
+            max_pix_variance_factor=maximum_pixel_method_variance(
+                self.semimajor, self.semiminor, self.theta),
+            correlation_lengths=calculate_correlation_lengths(
+                self.semimajor, self.semiminor),
             beamsize=calculate_beamsize(self.semimajor, self.semiminor)
         )
 
diff --git a/test/test_image.py b/test/test_image.py
index c42315a..4853a8b 100644
--- a/test/test_image.py
+++ b/test/test_image.py
@@ -298,35 +298,49 @@ class TestSimpleImageSourceFind(unittest.TestCase):
         at around 5 or 6 sigma detection level."""
 
         ew_sys_err, ns_sys_err = 0.0, 0.0
-        known_result = (
-            136.89603241069054, 14.022184792492785,  # RA, DEC
-            5.341819139061954e-4, 1.3428186757078464e-3,  # Err, Err
-            7.226590529214518e-4, 1.0918184742211533e-4,  # Peak flux, err
-            6.067963179204716e-4, 1.7037685531724465e-4,
-            # Integrated flux, err
-            6.192259965962862, 25.516190123153514,
-            # Significance level, Beam semimajor-axis width (arcsec)
-            10.718798843620489, 178.62899212789304,
-            # Beam semiminor-axis width (arcsec), Beam parallactic angle
-            ew_sys_err, ns_sys_err,
-            5.181697175052841,  # error_radius
-            1,  # fit_type
-            0.59184643302,  # chisq
-            0.67199741142,  # reduced chisq
-        )
+
+        known_result_fit = \
+            [1.36896042e+02, 1.40221872e+01,   # RA (deg), DEC (deg)
+             5.06084005e-04, 1.29061600e-03,  # Err, err
+             7.24671176e-04, 1.04806706e-04,  # Peak spectral brightness, err
+             6.03179622e-04, 1.62549622e-04,  # Flux density, err
+             6.44646215e+00, 2.55194168e+01,
+             # Significance level, beam semimajor-axis width (arcsec)
+             1.06461773e+01, 1.78499710e+02,
+             # Beam semiminor-axis width (arcsec), beam position angle (deg)
+             ew_sys_err, ns_sys_err,
+             4.97109604e+00, 1.00000000e+00,  # error_radius (arcsec), fit_type
+             6.03417635e-01, 6.67105734e-01]  # chisq, reduced chisq
+
+        known_result_moments = \
+            [1.3689603e+02, 1.4022377e+01,  # RA (deg), DEC (deg)
+             5.5378844e-04, 1.1825778e-03,  # Err, err
+             7.3612988e-04, 1.1431403e-04,  # Peak spectral brightness, err
+             6.0276804e-04, 1.6508212e-04,  # Flux density, err
+             6.4464622e+00, 2.4559519e+01,
+             # Significance level, beam semimajor-axis width (arcsec)
+             1.1146187e+01, 1.7876042e+02,  # Beam semiminor-axis width (arcsec),
+             # Beam position angle (deg).
+             ew_sys_err, ns_sys_err,
+             4.6760769e+00, 0.0000000e+00,  # error_radius (arcsec), fit_type
+             8.3038670e-01, 9.1803038e-01]  # chisq, reduced chisq
+
         self.image = accessors.sourcefinder_image_from_accessor(
             FitsImage(GRB120422A))
 
         results = self.image.extract(det=5, anl=3)
         results = [result.serialize(ew_sys_err, ns_sys_err) for result in
                    results]
-        # Our modified kappa,sigma clipper gives a slightly lower noise
-        # which catches an extra noise peak at the 5 sigma level.
         self.assertEqual(len(results), 2)
-        r = results[1]
-        self.assertEqual(len(r), len(known_result))
-        for i in range(len(r)):
-            self.assertAlmostEqual(r[i], known_result[i], places=0)
+        r = np.array(results[1], dtype=np.float32)
+        # Check if we derived source parameters from a fit or from moments.
+        if r[-3] == 1:
+            known_result = np.array(known_result_fit, dtype=np.float32)
+        else:
+            known_result = np.array(known_result_moments, dtype=np.float32)
+        self.assertEqual(r.size, known_result.size)
+        self.assertTrue(np.allclose(r, known_result, rtol=1e-4, atol=0,
+                                    equal_nan=True))
 
     @requires_data(GRB120422A)
     def testForceSourceShape(self):