Merge pull request #83 from lsst/tickets/DM-14155

DM-14155 make jointcal compatible with other source selectors

Readme.md

@@ -4,9 +4,13 @@ This package performs joint astrometry and photometry across multiple visits of
 
 It was derived from lsst_france/meas_simastrom
 
-## A note on tests
+## A note on numbers in tests
 
 All of the `tests/test_jointcal_<obs>.py` tests run jointcal on their respective `testdata_jointcal/<obs>` dataset.
-The numeric values being tested (the metrics, pa1, and relative/absolute rms) were not computed in advance, but are empirical, from the first time that test was successfully run.
-They are updated when changes to the algorithm warrants, e.g. a change to outlier rejection might affect the final chi2/ndof for some of the testdata.
+The numeric values being tested (the `metrics` dictionary, and the `pa1`/`relative_error`/`dist_absolute_rms` values) were not computed in advance, but are empirical, from the first time that test was successfully run.
+The one exception is the `cfht_minimal` dataset, which was designed to be small enough (3 `MeasuredStar`, 2 `FittedStar`, 1 `RefStar`) that the "correct" answers could be computed in advance with e.g. Mathematica.
+They are updated when changes to the algorithm warrants, e.g. a change to outlier rejection that affects the final chi2/ndof for some of the testdata.
+The `pa1` (see LPM-17 table 14) and relative/absolute rms (related to, but not the same as LPM-17's AM1 in table 18) values should be considered regression tests: algorithmic changes shouldn't result in larger values, but could produce smaller ones.
+They are "measured value must be less than" tests, and the calculations are in `lsst.jointcal.utils.JointcalStatistics.py`.
+
 More sophisticated testing of jointcal's validity is handled by `validate_drp` on much larger datasets that have multiple full-focal plane visits.

python/lsst/jointcal/jointcal.py

@@ -208,14 +208,19 @@ class JointcalConfig(pexConfig.Config):
         doc="Write initial/final fit files containing the contributions to chi2.",
         default=False
     )
+    sourceFluxType = pexConfig.Field(
+        dtype=str,
+        doc="Source flux field to use in source selection and to get fluxes from the catalog.",
+        default='Calib'
+    )
 
     def setDefaults(self):
         sourceSelector = self.sourceSelector["astrometry"]
         sourceSelector.setDefaults()
         # don't want to lose existing flags, just add to them.
         sourceSelector.badFlags.extend(["slot_Shape_flag"])
         # This should be used to set the FluxField value in jointcal::JointcalControl
-        sourceSelector.sourceFluxType = 'Calib'
+        sourceSelector.sourceFluxType = self.sourceFluxType
 
 
 class JointcalTask(pipeBase.CmdLineTask):
@@ -308,7 +313,7 @@ def _build_ccdImage(self, dataRef, associations, jointcalControl):
         fluxMag0 = calib.getFluxMag0()
         photoCalib = afwImage.PhotoCalib(1.0/fluxMag0[0], fluxMag0[1]/fluxMag0[0]**2, bbox)
 
-        goodSrc = self.sourceSelector.selectSources(src)
+        goodSrc = self.sourceSelector.run(src)
 
         if len(goodSrc.sourceCat) == 0:
             self.log.warn("No sources selected in visit %s ccd %s", visit, ccdId)
@@ -354,7 +359,7 @@ def run(self, dataRefs, profile_jointcal=False):
 
         exitStatus = 0  # exit status for shell
 
-        sourceFluxField = "slot_%sFlux" % (self.sourceSelector.config.sourceFluxType,)
+        sourceFluxField = "slot_%sFlux" % (self.config.sourceFluxType,)
         jointcalControl = lsst.jointcal.JointcalControl(sourceFluxField)
         associations = lsst.jointcal.Associations()
 

tests/test_astrometryModel.py

@@ -86,7 +86,7 @@ def setUp(self):
             dataRef = self.butler.dataRef('calexp', visit=visit, ccd=ccd)
 
             src = dataRef.get("src", flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS, immediate=True)
-            goodSrc = sourceSelector.selectSources(src)
+            goodSrc = sourceSelector.run(src)
             # Need memory contiguity to do vector-like things on the sourceCat.
             goodSrc = goodSrc.sourceCat.copy(deep=True)
 

tests/test_jointcal_cfht.py

@@ -30,8 +30,7 @@ def setUpClass(cls):
     def setUp(self):
         # We don't want the absolute astrometry to become significantly worse
         # than the single-epoch astrometry (about 0.040").
-        # This value was empirically determined from the first run of jointcal on
-        # this data, and will likely vary from survey to survey.
+        # See Readme for an explanation of this empirical value.
         self.dist_rms_absolute = 48.6e-3*u.arcsecond
 
         do_plot = False
@@ -58,9 +57,7 @@ def test_jointcalTask_2_visits(self):
         # to test whether we got the expected chi2 contribution files.
         self.other_args.extend(['--config', 'writeChi2ContributionFiles=True'])
 
-        # NOTE: The relative RMS limit was empirically determined from the
-        # first run of jointcal on this data. We should always do better than
-        # this in the future!
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 11e-3*u.arcsecond
         pa1 = 0.014
         metrics = {'collected_astrometry_refStars': 825,
@@ -100,9 +97,7 @@ def test_jointcalTask_2_visits_constrainedAstrometry_no_photometry(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
 
-        # NOTE: The relative RMS limit was empirically determined from the
-        # first run of jointcal on this data. We should always do better than
-        # this in the future!
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 12e-3*u.arcsecond
         dist_rms_absolute = 48e-3*u.arcsecond
         pa1 = None
@@ -125,6 +120,7 @@ def test_jointcalTask_2_visits_constrainedPhotometry_no_astrometry(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_astrometry = False
 
+        # See Readme for an explanation of these empirical values.
         pa1 = 0.017
         metrics = {'collected_photometry_refStars': 825,
                    'selected_photometry_refStars': 350,
@@ -137,6 +133,32 @@ def test_jointcalTask_2_visits_constrainedPhotometry_no_astrometry(self):
 
         self._testJointcalTask(2, None, None, pa1, metrics=metrics)
 
+    def test_jointcalTask_2_visits_constrainedPhotometry_flagged(self):
+        """Test the use of the FlaggedSourceSelector."""
+        self.config = lsst.jointcal.jointcal.JointcalConfig()
+        self.config.photometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
+        self.config.photometryModel = "constrained"
+        self.config.sourceSelector.name = "flagged"
+        # Reduce warnings due to flaggedSourceSelector having fewer sources than astrometrySourceSelector.
+        self.config.minMeasuredStarsPerCcd = 30
+        self.config.minRefStarsPerCcd = 20
+        self.config.doAstrometry = False
+        self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
+        self.jointcalStatistics.do_astrometry = False
+
+        # See Readme for an explanation of these empirical values.
+        pa1 = 0.026
+        metrics = {'collected_photometry_refStars': 825,
+                   'selected_photometry_refStars': 212,
+                   'associated_photometry_fittedStars': 270,
+                   'selected_photometry_fittedStars': 244,
+                   'selected_photometry_ccdImages': 12,
+                   'photometry_final_chi2': 369.96,
+                   'photometry_final_ndof': 252
+                   }
+
+        self._testJointcalTask(2, None, None, pa1, metrics=metrics)
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass

tests/test_jointcal_decam.py

@@ -28,8 +28,7 @@ def setUpClass(cls):
             raise unittest.SkipTest("testdata_jointcal not setup")
 
     def setUp(self):
-        # This value was empirically determined from the first run of jointcal on
-        # this data, and will likely vary from survey to survey.
+        # See Readme for an explanation of this empirical value.
         self.dist_rms_absolute = 62.5e-3*u.arcsecond
 
         do_plot = False
@@ -56,12 +55,9 @@ def test_jointcalTask_2_visits(self):
         self.config.astrometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
 
-        # NOTE: The relative RMS limit was empirically determined from the
-        # first run of jointcal on this data. We should always do better than
-        # this in the future!
+        # See Readme for an explanation of these empirical values.
         relative_error = 19e-3*u.arcsecond
         pa1 = 0.14
-        # NOTE: decam fits are currently not converging; the chi2 jumps around, so skip that Metric.
         metrics = {'collected_astrometry_refStars': 4866,
                    'collected_photometry_refStars': 4865,
                    'selected_astrometry_refStars': 661,
@@ -89,9 +85,7 @@ def test_jointcalTask_2_visits_constrainedAstrometry_no_photometry(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
 
-        # NOTE: The relative RMS limit was empirically determined from the
-        # first run of jointcal on this data. We should always do better than
-        # this in the future!
+        # See Readme for an explanation of these empirical values.
         relative_error = 17e-3*u.arcsecond
         pa1 = None
         metrics = {'collected_astrometry_refStars': 4866,
@@ -105,6 +99,55 @@ def test_jointcalTask_2_visits_constrainedAstrometry_no_photometry(self):
 
         self._testJointcalTask(2, relative_error, self.dist_rms_absolute, pa1, metrics=metrics)
 
+    def setup_jointcalTask_2_visits_constrainedPhotometry_no_astrometry(self):
+        """Help keep the constrainedPhotometry tests consistent and make
+        the differences between them more obvious.
+        """
+        self.config = lsst.jointcal.jointcal.JointcalConfig()
+        self.config.astrometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
+        self.config.photometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
+        self.config.photometryModel = "constrained"
+        self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
+        self.config.doAstrometry = False
+        self.jointcalStatistics.do_astrometry = False
+
+        # See Readme for an explanation of these empirical values.
+        pa1 = 0.11
+        metrics = {'collected_photometry_refStars': 4865,
+                   'selected_photometry_refStars': 661,
+                   'associated_photometry_fittedStars': 6749,
+                   'selected_photometry_fittedStars': 2044,
+                   'selected_photometry_ccdImages': 14,
+                   'photometry_final_chi2': 3066.22,
+                   'photometry_final_ndof': 1998,
+                   }
+
+        return pa1, metrics
+
+    @unittest.skip("DM-14439 : This test produces different chi2/ndof on Linux and macOS.")
+    def test_jointcalTask_2_visits_constrainedPhotometry_no_astrometry(self):
+        pa1, metrics = self.setup_jointcalTask_2_visits_constrainedPhotometry_no_astrometry()
+        self._testJointcalTask(2, None, None, pa1, metrics=metrics)
+
+    @unittest.skip("DM-14439 : This test produces different chi2/ndof on Linux and macOS.")
+    def test_jointcalTask_2_visits_constrainedPhotometry_flagged_selector(self):
+        pa1, metrics = self.setup_jointcalTask_2_visits_constrainedPhotometry_no_astrometry()
+        self.config.sourceSelector.name = 'flagged'
+        # Reduce warnings due to flaggedSourceSelector having fewer sources than astrometrySourceSelector.
+        self.config.minMeasuredStarsPerCcd = 40
+
+        # See Readme for an explanation of these empirical values.
+        metrics = {'collected_photometry_refStars': 4865,
+                   'selected_photometry_refStars': 551,
+                   'associated_photometry_fittedStars': 860,
+                   'selected_photometry_fittedStars': 593,
+                   'selected_photometry_ccdImages': 14,
+                   'photometry_final_chi2': 817.124,
+                   'photometry_final_ndof': 607,
+                   }
+
+        self._testJointcalTask(2, None, None, pa1, metrics=metrics)
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass

tests/test_jointcal_hsc.py

@@ -30,8 +30,7 @@ def setUpClass(cls):
             raise unittest.SkipTest("testdata_jointcal not setup")
 
     def setUp(self):
-        # This value was empirically determined from the first run of jointcal on
-        # this data, and will likely vary from survey to survey.
+        # See Readme for an explanation of this empirical value.
         self.dist_rms_absolute = 53e-3*u.arcsecond
 
         do_plot = False
@@ -54,9 +53,7 @@ def test_jointcalTask_2_visits(self):
         self.config.astrometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
 
-        # NOTE: The relative RMS limit was empirically determined from the
-        # first run of jointcal on this data. We should always do better than
-        # this in the future!
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 17e-3*u.arcsecond
         pa1 = 0.024
         metrics = {'collected_astrometry_refStars': 2187,
@@ -87,9 +84,7 @@ def test_jointcalTask_11_visits_no_photometry(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
 
-        # NOTE: The relative RMS limit was empirically determined from the
-        # first run of jointcal on this data. We should always do better than
-        # this in the future!
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 17e-3*u.arcsecond
         pa1 = None  # pa1 = 0.134
         metrics = {'collected_astrometry_refStars': 3649,
@@ -104,6 +99,7 @@ def test_jointcalTask_11_visits_no_photometry(self):
 
     def testJointcalTask_2_visits_no_astrometry(self):
         """Test turning off fitting astrometry."""
+        # See Readme for an explanation of these empirical values.
         pa1 = 0.024
         metrics = {'collected_photometry_refStars': 2187,
                    'selected_photometry_refStars': 515,
@@ -139,6 +135,7 @@ def testJointcalTask_2_visits_no_astrometry(self):
 
     def testJointcalTask_2_visits_no_photometry(self):
         """Test turning off fitting photometry."""
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 17e-3*u.arcsecond
         metrics = {'collected_astrometry_refStars': 2187,
                    'selected_astrometry_refStars': 515,
@@ -182,6 +179,8 @@ def test_jointcalTask_2_visits_gaia_refcat(self):
         test_config = os.path.join(lsst.utils.getPackageDir('jointcal'), 'tests/config/hsc-config.py')
         self.other_args.extend(['--configfile', test_config])
         dist_rms_relative = 17e-3*u.arcsecond
+
+        # See Readme for an explanation of these empirical values.
         # NOTE: PA1 is slightly different here, because the number of SDSS
         # cross-matches within 0.1" goes down after we apply the GAIA-fit WCS.
         pa1 = 0.02405
@@ -214,6 +213,7 @@ def test_jointcalTask_2_visits_no_photometry_match_cut_10(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
 
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 17e-3*u.arcsecond
         metrics = {'collected_astrometry_refStars': 2187,
                    'selected_astrometry_refStars': 546,
@@ -235,6 +235,7 @@ def test_jointcalTask_3_visits_no_photometry(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
 
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 17e-3*u.arcsecond
         metrics = {'collected_astrometry_refStars': 2187,
                    'selected_astrometry_refStars': 541,
@@ -258,6 +259,7 @@ def test_jointcalTask_3_visits_no_photometry_min_measurements_3(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
 
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 17e-3*u.arcsecond
         metrics = {'collected_astrometry_refStars': 2187,
                    'selected_astrometry_refStars': 541,

tests/test_jointcal_lsstSim.py

@@ -33,8 +33,7 @@ def setUpClass(cls):
     def setUp(self):
         # We don't want the absolute astrometry to become significantly worse
         # than the single-epoch astrometry (about 0.040").
-        # This value was empirically determined from the first run of jointcal on
-        # this data, and will likely vary from survey to survey.
+        # See Readme for an explanation of this empirical value.
         self.dist_rms_absolute = 42e-3*u.arcsecond
 
         do_plot = False
@@ -78,6 +77,8 @@ def testJointcalTask_2_visits(self):
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.config.doPhotometry = False
         self.jointcalStatistics.do_photometry = False
+
+        # See Readme for an explanation of these empirical values.
         # pa1 = 2.64e-3
         # 'collected_photometry_refStars': 1686,
         # 'selected_photometry_refStars': 1686,
@@ -107,6 +108,8 @@ def testJointcalTask_10_visits(self):
         self.config.doPhotometry = False
         self.config.sourceSelector['astrometry'].badFlags.append("base_PixelFlags_flag_interpolated")
         self.jointcalStatistics.do_photometry = False
+
+        # See Readme for an explanation of these empirical values.
         # pa1 = 2.64e-3
         # 'collected_photometry_refStars': 1686,
         # 'selected_photometry_refStars': 1686,
@@ -128,6 +131,7 @@ def testJointcalTask_10_visits(self):
     @unittest.skip("A.net reference catalog missing flux errors. Unskip once DM-11397 is fixed.")
     def testJointcalTask_2_visits_no_astrometry(self):
         """Test turning off fitting astrometry."""
+        # See Readme for an explanation of these empirical values.
         pa1 = 2.64e-3
         metrics = {'collected_photometry_refStars': 1686,
                    'selected_photometry_refStars': 176,
@@ -163,6 +167,7 @@ def testJointcalTask_2_visits_no_astrometry(self):
 
     def testJointcalTask_2_visits_no_photometry(self):
         """Test turning off fitting photometry."""
+        # See Readme for an explanation of these empirical values.
         dist_rms_relative = 9.7e-3*u.arcsecond
         metrics = {'collected_astrometry_refStars': 1686,
                    'selected_astrometry_refStars': 176,