DM-11620: Numerical problem with Gaussian test in meas_base with pytest

python/lsst/meas/base/SConscript

@@ -1,4 +1,4 @@
-## -*- python -*-
+# -*- python -*-
 from lsst.sconsUtils import scripts
 scripts.BasicSConscript.pybind11(['algorithm',
                                   'apertureFlux',
@@ -20,4 +20,4 @@ scripts.BasicSConscript.pybind11(['algorithm',
                                   'sdssShape',
                                   'sincCoeffs',
                                   'shapeUtilities',
-                                  'transform',], addUnderscore=False)
+                                  'transform', ], addUnderscore=False)

python/lsst/meas/base/forcedMeasurement.py

@@ -303,7 +303,8 @@ def run(self, measCat, exposure, refCat, refWcs, exposureId=None, beginOrder=Non
         footprints = {ref.getId(): (ref.getParent(), measRecord.getFootprint())
                       for (ref, measRecord) in zip(refCat, measCat)}
 
-        self.log.info("Performing forced measurement on %d sources" % (len(refCat),))
+        self.log.info("Performing forced measurement on %d source%s", len(refCat),
+                      "" if len(refCat) == 1 else "s")
 
         if self.config.doReplaceWithNoise:
             noiseReplacer = NoiseReplacer(self.config.noiseReplacer, exposure,

python/lsst/meas/base/sfm.py

@@ -340,8 +340,12 @@ def runPlugins(self, noiseReplacer, measCat, exposure, beginOrder=None, endOrder
         # Loop through all the parent sources, first processing the children, then the parent
         measParentCat = measCat.getChildren(0)
 
-        self.log.info("Measuring %d sources (%d parents, %d children) "
-                      % (len(measCat), len(measParentCat), len(measCat) - len(measParentCat)))
+        nMeasCat = len(measCat)
+        nMeasParentCat = len(measParentCat)
+        self.log.info("Measuring %d source%s (%d parent%s, %d child%s) ",
+                      nMeasCat, ("" if nMeasCat == 1 else "s"),
+                      nMeasParentCat, ("" if nMeasParentCat == 1 else "s"),
+                      nMeasCat - nMeasParentCat, ("" if nMeasCat - nMeasParentCat == 1 else "ren"))
 
         for parentIdx, measParentRecord in enumerate(measParentCat):
             # first get all the children of this parent, insert footprint in turn, and measure

python/lsst/meas/base/tests.py

@@ -181,7 +181,7 @@ def makeMinimalSchema(cls):
     def makePerturbedWcs(oldWcs, minScaleFactor=1.2, maxScaleFactor=1.5,
                          minRotation=None, maxRotation=None,
                          minRefShift=None, maxRefShift=None,
-                         minPixShift=2.0, maxPixShift=4.0):
+                         minPixShift=2.0, maxPixShift=4.0, randomSeed=None):
         """!
         Create a new undistorted TanWcs that is similar but not identical to another, with random
         scaling, rotation, and offset (in both pixel position and reference position).
@@ -195,7 +195,11 @@ def makePerturbedWcs(oldWcs, minScaleFactor=1.2, maxScaleFactor=1.5,
         The default range for rotation is 30-60 degrees, and the default range for reference shift
         is 0.5-1.0 arcseconds (these cannot be safely included directly as default values because Angle
         objects are mutable).
+
+        The random number generator is primed with the seed given. If ``None``, a seed is
+        automatically chosen.
         """
+        random_state = np.random.RandomState(randomSeed)
         if minRotation is None:
             minRotation = 30.0*lsst.afw.geom.degrees
         if maxRotation is None:
@@ -210,10 +214,10 @@ def splitRandom(min1, max1, min2=None, max2=None):
                 min2 = -max1
             if max2 is None:
                 max2 = -min1
-            if np.random.uniform() > 0.5:
-                return float(np.random.uniform(min1, max1))
+            if random_state.uniform() > 0.5:
+                return float(random_state.uniform(min1, max1))
             else:
-                return float(np.random.uniform(min2, max2))
+                return float(random_state.uniform(min2, max2))
         # Generate random perturbations
         scaleFactor = splitRandom(minScaleFactor, maxScaleFactor, 1.0/maxScaleFactor, 1.0/minScaleFactor)
         rotation = splitRandom(minRotation.asRadians(), maxRotation.asRadians())*lsst.afw.geom.radians
@@ -415,7 +419,7 @@ def transform(self, wcs, **kwds):
             result.addSource(newFlux, newCentroid, newDeconvolvedShape)
         return result
 
-    def realize(self, noise, schema):
+    def realize(self, noise, schema, randomSeed=None):
         """!
         Create a simulated with noise and a simulated post-detection catalog with (Heavy)Footprints.
 
@@ -424,16 +428,18 @@ def realize(self, noise, schema):
         @param[in]   schema     Schema of the new catalog to be created.  Must start with self.schema (i.e.
                                 schema.contains(self.schema) must be True), but typically contains fields for
                                 already-configured measurement algorithms as well.
+        @param[in]   randomSeed Seed for the random number generator. If None, a seed is chosen automatically.
 
         @return a tuple of (exposure, catalog)
         """
+        random_state = np.random.RandomState(randomSeed)
         assert schema.contains(self.schema)
         mapper = lsst.afw.table.SchemaMapper(self.schema)
         mapper.addMinimalSchema(self.schema, True)
         exposure = self.exposure.clone()
         exposure.getMaskedImage().getVariance().getArray()[:, :] = noise**2
         exposure.getMaskedImage().getImage().getArray()[:, :] \
-            += np.random.randn(exposure.getHeight(), exposure.getWidth())*noise
+            += random_state.randn(exposure.getHeight(), exposure.getWidth())*noise
         catalog = lsst.afw.table.SourceCatalog(schema)
         catalog.extend(self.catalog, mapper=mapper)
         # Loop over sources and generate new HeavyFootprints that divide up the noisy pixels, not the
@@ -466,20 +472,6 @@ def realize(self, noise, schema):
 
 
 class AlgorithmTestCase(object):
-    # Some tests depend on the noise realization in the test data or from the
-    # np.random number generator. In most cases, they are testing that the
-    # measured flux lies within 2 sigma of the correct value, which we should
-    # expect to fail sometimes. Some -- but sadly not all -- of these cases
-    # have been marked with an "rng dependent" comment.
-    #
-    # We ensure these tests are provided with data which causes them to pass
-    # by seeding the np.RNG with this value. It can be over-ridden as
-    # necessary in subclasses.
-    randomSeed = 1234
-
-    @classmethod
-    def setUpClass(cls):
-        np.random.seed(cls.randomSeed)
 
     def makeSingleFrameMeasurementConfig(self, plugin=None, dependencies=()):
         """Convenience function to create a Config instance for SingleFrameMeasurementTask

tests/SConscript

@@ -1,11 +1,11 @@
 # -*- python -*-
 from lsst.sconsUtils import scripts
 
-pybind11_test_modules = ['sillyCentroid',]
+pybind11_test_modules = ['sillyCentroid']
 noBuildList = [name + '.cc' for name in pybind11_test_modules]
 ignoreList = [name + '.py' for name in pybind11_test_modules]
 ignoreList.append('testLib.py')
 
 scripts.BasicSConscript.pybind11(pybind11_test_modules)
 scripts.BasicSConscript.tests(noBuildList=noBuildList,
-                              ignoreList=ignoreList)
+                              ignoreList=ignoreList, pyList=[])

tests/test_ApertureFlux.py

@@ -171,7 +171,7 @@ def testSingleFramePlugin(self):
         ctrl = config.plugins[baseName].makeControl()
         algMetadata = lsst.daf.base.PropertyList()
         task = self.makeSingleFrameMeasurementTask(config=config, algMetadata=algMetadata)
-        exposure, catalog = self.dataset.realize(10.0, task.schema)
+        exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=0)
         task.run(catalog, exposure)
         radii = algMetadata.get("%s_radii" % (baseName,))
         self.assertEqual(list(radii), list(ctrl.radii))
@@ -222,9 +222,9 @@ def testForcedPlugin(self):
         algMetadata = lsst.daf.base.PropertyList()
         task = self.makeForcedMeasurementTask(baseName, algMetadata=algMetadata)
         radii = algMetadata.get("%s_radii" % (baseName,))
-        measWcs = self.dataset.makePerturbedWcs(self.dataset.exposure.getWcs())
+        measWcs = self.dataset.makePerturbedWcs(self.dataset.exposure.getWcs(), randomSeed=1)
         measDataset = self.dataset.transform(measWcs)
-        exposure, truthCatalog = measDataset.realize(10.0, measDataset.makeMinimalSchema())
+        exposure, truthCatalog = measDataset.realize(10.0, measDataset.makeMinimalSchema(), randomSeed=1)
         refCat = self.dataset.catalog
         refWcs = self.dataset.exposure.getWcs()
         measCat = task.generateMeasCat(exposure, refCat, refWcs)

tests/test_Blendedness.py

@@ -89,7 +89,7 @@ def testBlendedness(self):
         Check that we measure a positive blendedness for two overlapping sources
         """
         task = self.makeSingleFrameMeasurementTask("base_Blendedness")
-        exposure, catalog = self.dataset.realize(10.0, task.schema)
+        exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=0)
         task.run(catalog, exposure)
         self.assertGreater(catalog[1].get('base_Blendedness_abs_flux'), 0)
         self.assertGreater(catalog[2].get('base_Blendedness_abs_flux'), 0)

tests/test_CentroidChecker.py

@@ -128,7 +128,7 @@ def testNoError(self):
         config = self.makeConfig()
         config.slots.centroid = "truth"
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=0)
         task.run(cat, exposure)
         source = cat[0]
         self.assertFalse(source.get("test_Centroider_flag"))
@@ -145,7 +145,7 @@ def testCentroidDistance(self):
         config.plugins[self.algName].moveX = -2
         config.plugins[self.algName].dist = 1
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=1)
         source = cat[0]
         task.run(cat, exposure)
         self.assertTrue(source.get("test_Centroider_flag"))
@@ -162,7 +162,7 @@ def testCentroidOutsideFootprint(self):
         config.slots.centroid = "truth"
         config.plugins[self.algName].moveX = -30
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=2)
         source = cat[0]
         task.run(cat, exposure)
         self.assertTrue(source.get("test_Centroider_flag"))
@@ -177,7 +177,7 @@ def testNaiveCentroid(self):
         config = self.makeConfig("base_NaiveCentroid")
         config.plugins["base_NaiveCentroid"].maxDistToPeak = .0001
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=3)
         source = cat[0]
         task.run(cat, exposure)
         self.assertTrue(source.get("base_NaiveCentroid_flag"))
@@ -191,7 +191,7 @@ def testSdssCentroid(self):
         config = self.makeConfig("base_SdssCentroid")
         config.plugins["base_SdssCentroid"].maxDistToPeak = .0001
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=4)
         source = cat[0]
         task.run(cat, exposure)
         self.assertTrue(source.get("base_SdssCentroid_flag"))
@@ -205,7 +205,7 @@ def testGaussianCentroid(self):
         config = self.makeConfig("base_GaussianCentroid")
         config.plugins["base_GaussianCentroid"].maxDistToPeak = .0001
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=5)
         source = cat[0]
         task.run(cat, exposure)
         self.assertTrue(source.get("base_GaussianCentroid_flag"))

tests/test_Classification.py

@@ -52,7 +52,7 @@ def testSingleFramePlugin(self):
         config.slots.modelFlux = "truth"
         task = self.makeSingleFrameMeasurementTask(config=config)
         abTask = catCalc.CatalogCalculationTask(schema=task.schema)
-        exposure, catalog = self.dataset.realize(10.0, task.schema)
+        exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=0)
         task.run(catalog, exposure)
         abTask.run(catalog)
         self.assertLess(catalog[0].get("base_ClassificationExtendedness_value"), 0.5)
@@ -79,7 +79,7 @@ def runFlagTest(psfFlux=100.0, modelFlux=200.0,
                         psfFluxFlag=False, modelFluxFlag=False):
             task = self.makeSingleFrameMeasurementTask(config=config)
             abTask = catCalc.CatalogCalculationTask(schema=task.schema, config=abConfig)
-            exposure, catalog = self.dataset.realize(10.0, task.schema)
+            exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=1)
             source = catalog[0]
             source.set("base_PsfFlux_flux", psfFlux)
             source.set("base_PsfFlux_fluxSigma", psfFluxSigma)

tests/test_FPPosition.py

@@ -53,7 +53,7 @@ def tearDown(self):
 
     def testFPPosition(self):
         task = self.makeSingleFrameMeasurementTask("base_FPPosition")
-        exposure, catalog = self.dataset.realize(10.0, task.schema)
+        exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=0)
         exposure.setDetector(self.dw.detector)
         task.run(catalog, exposure)
         pointKey = lsst.afw.table.Point2DKey(catalog.schema["base_FPPosition"])

tests/test_FlagHandler.py

@@ -253,7 +253,7 @@ def testPluginNoError(self):
         """
         schema = self.dataset.makeMinimalSchema()
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=self.config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=0)
         task.run(cat, exposure)
         source = cat[0]
         self.assertFalse(source.get(self.algName + "_flag"))
@@ -269,7 +269,7 @@ def testPluginUnexpectedError(self):
         self.config.plugins[self.algName].flux0 = 0.0     # this causes a divide by zero
         schema = self.dataset.makeMinimalSchema()
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=self.config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=1)
         task.log.setLevel(task.log.FATAL)
         task.run(cat, exposure)
         source = cat[0]
@@ -283,7 +283,7 @@ def testPluginContainsNan(self):
         """
         schema = self.dataset.makeMinimalSchema()
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=self.config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=2)
         source = cat[0]
         exposure.getMaskedImage().getImage().getArray()[int(source.getY()), int(source.getX())] = np.nan
         task.run(cat, exposure)
@@ -297,7 +297,7 @@ def testPluginEdgeError(self):
         """
         schema = self.dataset.makeMinimalSchema()
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=self.config)
-        exposure, cat = self.dataset.realize(noise=100.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=100.0, schema=schema, randomSeed=3)
         # Set the size large enough to trigger the edge error
         self.config.plugins[self.algName].size = exposure.getDimensions()[1]//2
         task.log.setLevel(task.log.FATAL)
@@ -315,7 +315,7 @@ def testSafeCentroider(self):
         schema = self.dataset.makeMinimalSchema()
         task = lsst.meas.base.SingleFrameMeasurementTask(schema=schema, config=self.config)
         task.log.setLevel(task.log.FATAL)
-        exposure, cat = self.dataset.realize(noise=0.0, schema=schema)
+        exposure, cat = self.dataset.realize(noise=0.0, schema=schema, randomSeed=4)
         source = cat[0]
         task.run(cat, exposure)
         self.assertFalse(source.get(self.algName + "_flag"))

tests/test_FootprintArea.py

@@ -48,7 +48,7 @@ def testSingleFramePlugin(self):
         config.plugins.names = set(["base_FootprintArea"])
         task = lsst.meas.base.CatalogCalculationTask(config=config,
                                                      schema=schema)
-        exposure, catalog = self.dataset.realize(10.0, task.schema)
+        exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=0)
         task.run(catalog)
         record = catalog[0]
         self.assertEqual(record.getFootprint().getArea(),

tests/test_GaussianCentroid.py

@@ -45,7 +45,7 @@ def tearDown(self):
 
     def testSingleFramePlugin(self):
         task = self.makeSingleFrameMeasurementTask("base_GaussianCentroid")
-        exposure, catalog = self.dataset.realize(10.0, task.schema)
+        exposure, catalog = self.dataset.realize(10.0, task.schema, randomSeed=0)
         task.run(catalog, exposure)
         record = catalog[0]
         x = record.get("base_GaussianCentroid_x")