DM-12602 fitting bugfixes and other cleanups

include/lsst/jointcal/ConstrainedPolyModel.h

@@ -26,12 +26,12 @@ namespace jointcal {
  * transformation depending on the chip number (instrument model) and a
  * transformation per visit (anamorphism). The two-transformation Mapping
  * required for this model is TwoTransfoMapping. This modeling of distortions
- * is meant for set of images from a single mosaic imager.
+ * is meant for a set of images from a single mosaic imager.
  */
 class ConstrainedPolyModel : public AstrometryModel {
 public:
     ConstrainedPolyModel(CcdImageList const &ccdImageList, ProjectionHandler const *projectionHandler,
-                         bool initFromWCS, unsigned nNotFit = 0);
+                         bool initFromWCS, unsigned nNotFit = 0, int chipDegree = 3, int visitDegree = 2);
 
     /// No copy or move: there is only ever one instance of a given model (i.e. per ccd+visit)
     ConstrainedPolyModel(ConstrainedPolyModel const &) = delete;

python/lsst/jointcal/astrometryModels.cc

@@ -55,8 +55,9 @@ void declareConstrainedPolyModel(py::module &mod) {
     py::class_<ConstrainedPolyModel, std::shared_ptr<ConstrainedPolyModel>, AstrometryModel> cls(
             mod, "ConstrainedPolyModel");
 
-    cls.def(py::init<CcdImageList const &, const ProjectionHandler *, bool, unsigned>(), "ccdImageList"_a,
-            "projectionHandler"_a, "initFromWcs"_a, "nNotFit"_a = 0);
+    cls.def(py::init<CcdImageList const &, const ProjectionHandler *, bool, unsigned, int, int>(),
+            "ccdImageList"_a, "projectionHandler"_a, "initFromWcs"_a, "nNotFit"_a = 0, "chipDegree"_a = 3,
+            "visitDegree"_a = 2);
 }
 
 PYBIND11_PLUGIN(astrometryModels) {

python/lsst/jointcal/dataIds.py

@@ -98,7 +98,8 @@ def makeDataRefList(self, namespace):
                 tract = dataId.pop("tract")
                 # making a DataRef for src fills out any missing keys and allows us to iterate
                 for ref in namespace.butler.subset("src", dataId=dataId):
-                    self._addDataRef(namespace, ref.dataId, tract)
+                    if ref.datasetExists():
+                        self._addDataRef(namespace, ref.dataId, tract)
 
         # Ensure all components of a visit are kept together by putting them all in the same set of tracts
         # NOTE: sorted() here is to keep py2 and py3 dataRefs in the same order.

python/lsst/jointcal/jointcal.py

@@ -130,8 +130,18 @@ class JointcalConfig(pexConfig.Config):
         dtype=int,
         default=2,
     )
-    polyOrder = pexConfig.Field(
-        doc="Polynomial order for fitting distorsion",
+    astrometrySimpleDegree = pexConfig.Field(
+        doc="Polynomial degree for fitting the simple astrometry model.",
+        dtype=int,
+        default=3,
+    )
+    astrometryChipDegree = pexConfig.Field(
+        doc="Degree of the per-chip transform for the constrained astrometry model.",
+        dtype=int,
+        default=2,
+    )
+    astrometryVisitDegree = pexConfig.Field(
+        doc="Degree of the per-visit transform for the constrained astrometry model.",
         dtype=int,
         default=3,
     )
@@ -262,7 +272,7 @@ def _build_ccdImage(self, dataRef, associations, jointcalControl):
 
         visitInfo = dataRef.get('calexp_visitInfo')
         detector = dataRef.get('calexp_detector')
-        ccdname = detector.getId()
+        ccdId = detector.getId()
         calib = dataRef.get('calexp_calib')
         tanWcs = dataRef.get('calexp_wcs')
         bbox = dataRef.get('calexp_bbox')
@@ -274,15 +284,15 @@ def _build_ccdImage(self, dataRef, associations, jointcalControl):
         goodSrc = self.sourceSelector.selectSources(src)
 
         if len(goodSrc.sourceCat) == 0:
-            self.log.warn("No sources selected in visit %s ccd %s", visit, ccdname)
+            self.log.warn("No sources selected in visit %s ccd %s", visit, ccdId)
         else:
-            self.log.info("%d sources selected in visit %d ccd %d", len(goodSrc.sourceCat), visit, ccdname)
+            self.log.info("%d sources selected in visit %d ccd %d", len(goodSrc.sourceCat), visit, ccdId)
         associations.addImage(goodSrc.sourceCat, tanWcs, visitInfo, bbox, filterName, photoCalib, detector,
-                              visit, ccdname, jointcalControl)
+                              visit, ccdId, jointcalControl)
 
         Result = collections.namedtuple('Result_from_build_CcdImage', ('wcs', 'key', 'filter'))
         Key = collections.namedtuple('Key', ('visit', 'ccd'))
-        return Result(tanWcs, Key(visit, ccdname), filterName)
+        return Result(tanWcs, Key(visit, ccdId), filterName)
 
     @pipeBase.timeMethod
     def run(self, dataRefs, profile_jointcal=False):
@@ -572,11 +582,13 @@ def _fit_astrometry(self, associations, dataName=None):
 
         if self.config.astrometryModel == "constrainedPoly":
             model = lsst.jointcal.ConstrainedPolyModel(associations.getCcdImageList(),
-                                                       sky_to_tan_projection, True, 0)
+                                                       sky_to_tan_projection, True, 0,
+                                                       chipDegree=self.config.astrometryChipDegree,
+                                                       visitDegree=self.config.astrometryVisitDegree)
         elif self.config.astrometryModel == "simplePoly":
             model = lsst.jointcal.SimplePolyModel(associations.getCcdImageList(),
                                                   sky_to_tan_projection,
-                                                  True, 0, self.config.polyOrder)
+                                                  True, 0, self.config.astrometrySimpleDegree)
 
         fit = lsst.jointcal.AstrometryFit(associations, model, self.config.posError)
         chi2 = fit.computeChi2()

src/ConstrainedPolyModel.cc

@@ -25,20 +25,14 @@ routines AstrometryFit needs to what is needed for this two-transfo model.
 The two-transfo mappings are implemented using two one-transfo
 mappings.*/
 
-// TODO : separate the polynomial degrees for chip and visit transfos.
-// TODO propagate those into python:
-static int DistortionDegree = 3;
-
 using namespace std;
 
 ConstrainedPolyModel::ConstrainedPolyModel(CcdImageList const &ccdImageList,
                                            ProjectionHandler const *projectionHandler, bool initFromWCS,
-                                           unsigned nNotFit)
+                                           unsigned nNotFit, int chipDegree, int visitDegree)
         : _sky2TP(projectionHandler)
 
 {
-    // from datacards (or default)
-    unsigned degree = DistortionDegree;
     // first loop to initialize all visit  and chip transfos.
     for (auto &ccdImage : ccdImageList) {
         const CcdImage &im = *ccdImage;
@@ -51,15 +45,15 @@ ConstrainedPolyModel::ConstrainedPolyModel(CcdImageList const &ccdImageList,
                         std::unique_ptr<SimpleGtransfoMapping>(new SimpleGtransfoMapping(GtransfoIdentity()));
             } else {
                 _visitMap[visit] = std::unique_ptr<SimpleGtransfoMapping>(
-                        new SimplePolyMapping(GtransfoLin(), GtransfoPoly(degree)));
+                        new SimplePolyMapping(GtransfoLin(), GtransfoPoly(visitDegree)));
             }
         }
         auto chipp = _chipMap.find(chip);
         if (chipp == _chipMap.end()) {
             const Frame &frame = im.getImageFrame();
 
             _tpFrame += applyTransfo(frame, *im.getPix2CommonTangentPlane(), LargeFrame);
-            GtransfoPoly pol(im.getPix2TangentPlane(), frame, degree);
+            GtransfoPoly pol(im.getPix2TangentPlane(), frame, chipDegree);
             GtransfoLin shiftAndNormalize = normalizeCoordinatesTransfo(frame);
 
             _chipMap[chip] = std::unique_ptr<SimplePolyMapping>(
@@ -77,7 +71,7 @@ ConstrainedPolyModel::ConstrainedPolyModel(CcdImageList const &ccdImageList,
             LOGLS_WARN(_log, "Chip " << chip << " is missing in the reference exposure, expect troubles.");
             GtransfoLin norm = normalizeCoordinatesTransfo(im.getImageFrame());
             _chipMap[chip] =
-                    std::unique_ptr<SimplePolyMapping>(new SimplePolyMapping(norm, GtransfoPoly(degree)));
+                    std::unique_ptr<SimplePolyMapping>(new SimplePolyMapping(norm, GtransfoPoly(chipDegree)));
         }
         _mappings[&im] = std::unique_ptr<TwoTransfoMapping>(
                 new TwoTransfoMapping(_chipMap[chip].get(), _visitMap[visit].get()));
@@ -206,7 +200,7 @@ std::shared_ptr<TanSipPix2RaDec> ConstrainedPolyModel::produceSipWcs(CcdImage co
     } catch (std::bad_cast &) {
         try {
             const GtransfoPoly &t2_poly = dynamic_cast<const GtransfoPoly &>(mapping->getTransfo2());
-            pix2Tp = t1 * t2_poly;
+            pix2Tp = t2_poly * t1;
         } catch (std::bad_cast &) {
             LOGLS_ERROR(_log, "Problem with transform 2 of ccd/visit " << ccdImage.getCcdId() << "/"
                                                                        << ccdImage.getVisit() << ": T2 "

src/FitterBase.cc

@@ -153,7 +153,7 @@ MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaC
         offsetParams(delta);
         Chi2Statistic currentChi2(computeChi2());
         LOGLS_DEBUG(_log, currentChi2);
-        if (currentChi2.chi2 > oldChi2) {
+        if (currentChi2.chi2 > oldChi2 && totalOutliers != 0) {
             LOGL_WARN(_log, "chi2 went up, skipping outlier rejection loop");
             returnCode = MinimizeResult::Chi2Increased;
             break;

tests/test_jointcal_cfht.py

@@ -105,18 +105,19 @@ def test_jointcalTask_2_visits_constrainedPoly(self):
         # 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!
-        dist_rms_relative = 12e-3*u.arcsecond
+        dist_rms_relative = 16e-3*u.arcsecond
+        dist_rms_absolute = 54e-3*u.arcsecond
         pa1 = None
         metrics = {'collectedAstrometryRefStars': 825,
                    'selectedAstrometryRefStars': 825,
                    'associatedAstrometryFittedStars': 2269,
                    'selectedAstrometryFittedStars': 1239,
                    'selectedAstrometryCcdImageList': 12,
-                   'astrometryFinalChi2': 1241.6,
-                   'astrometryFinalNdof': 2640,
+                   'astrometryFinalChi2': 1376.24,
+                   'astrometryFinalNdof': 2630,
                    }
 
-        self._testJointcalTask(2, dist_rms_relative, self.dist_rms_absolute, pa1, metrics=metrics)
+        self._testJointcalTask(2, dist_rms_relative, dist_rms_absolute, pa1, metrics=metrics)
 
     def test_jointcalTask_2_visits_constrainedPhotometry_no_astrometry(self):
         self.config = lsst.jointcal.jointcal.JointcalConfig()

tests/test_jointcal_decam.py

@@ -79,7 +79,7 @@ def test_jointcalTask_2_visits(self):
                    'astrometryFinalChi2': None,
                    'astrometryFinalNdof': 4306,
                    'photometryFinalChi2': None,
-                   'photometryFinalNdof': 2391,
+                   'photometryFinalNdof': 2333,
                    }
 
         self._testJointcalTask(2, relative_error, self.dist_rms_absolute, pa1, metrics=metrics)