DM-16235: Jointcal PhotoCalib returns negative mean calibrations

include/lsst/jointcal/AstrometryFit.h

@@ -101,6 +101,9 @@ class AstrometryFit : public FitterBase {
     /// @copydoc FitterBase::saveChi2RefContributions
     void saveChi2RefContributions(std::string const &baseName) const override;
 
+    /// Return the model being fit.
+    std::shared_ptr<AstrometryModel> getModel() const { return _astrometryModel; }
+
     /**
      * DEBUGGING routine
      */

include/lsst/jointcal/AstrometryModel.h

@@ -4,6 +4,7 @@
 
 #include "memory"
 
+#include "lsst/jointcal/CcdImage.h"
 #include "lsst/jointcal/Gtransfo.h"
 #include "lsst/jointcal/Eigenstuff.h"
 #include "lsst/jointcal/AstrometryMapping.h"
@@ -61,6 +62,17 @@ class AstrometryModel {
 
     virtual ~AstrometryModel(){};
 
+    /**
+     * Return true if this is a "reasonable" model.
+     *
+     * Not yet implemented: DM-16324
+     * An example might be that the model produces finite RA/Dec on each sensor's bounding box.
+     *
+     * @param ccdImageList The ccdImages to test the model validity on.
+     * @return True if the model is valid on all ccdImages.
+     */
+    bool validate(CcdImageList const &ccdImageList) const { return true; }
+
 protected:
     /// Return a pointer to the mapping associated with this ccdImage.
     virtual AstrometryMapping *findMapping(CcdImage const &ccdImage) const = 0;

include/lsst/jointcal/ConstrainedPhotometryModel.h

@@ -30,12 +30,14 @@ class ConstrainedPhotometryModel : public PhotometryModel {
      * @param      ccdImageList   The list of CCDImages to construct the model for.
      * @param      focalPlaneBBox The bounding box of the camera's focal plane, defining the domain of the
      *                            visit polynomial.
+     * @param[in]  log An lsst::log::Log instance to log messages to.
      * @param[in]  visitOrder    The order of the visit polynomial.
+     * @param[in]  errorPedestal_ A pedestal in flux or magnitude to apply to all MeasuredStar flux errors.
      */
     explicit ConstrainedPhotometryModel(CcdImageList const &ccdImageList,
-                                        afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7,
-                                        double errorPedestal_ = 0)
-            : PhotometryModel(errorPedestal_), _fittingChips(false), _fittingVisits(false) {
+                                        afw::geom::Box2D const &focalPlaneBBox, LOG_LOGGER log,
+                                        int visitOrder = 7, double errorPedestal_ = 0)
+            : PhotometryModel(log, errorPedestal_), _fittingChips(false), _fittingVisits(false) {
         _chipVisitMap.reserve(ccdImageList.size());
     }
 
@@ -103,7 +105,9 @@ class ConstrainedFluxModel : public ConstrainedPhotometryModel {
 public:
     explicit ConstrainedFluxModel(CcdImageList const &ccdImageList, afw::geom::Box2D const &focalPlaneBBox,
                                   int visitOrder = 7, double errorPedestal_ = 0)
-            : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder, errorPedestal_) {
+            : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox,
+                                         LOG_GET("jointcal.ConstrainedFluxModel"), visitOrder,
+                                         errorPedestal_) {
         initialize<FluxTransfoSpatiallyInvariant, FluxTransfoChebyshev, ChipVisitFluxMapping>(
                 ccdImageList, focalPlaneBBox, visitOrder);
     }
@@ -145,7 +149,9 @@ class ConstrainedMagnitudeModel : public ConstrainedPhotometryModel {
     explicit ConstrainedMagnitudeModel(CcdImageList const &ccdImageList,
                                        afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7,
                                        double errorPedestal_ = 0)
-            : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder, errorPedestal_) {
+            : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox,
+                                         LOG_GET("jointcal.ConstrainedMagnitudeModel"), visitOrder,
+                                         errorPedestal_) {
         initialize<MagnitudeTransfoSpatiallyInvariant, MagnitudeTransfoChebyshev, ChipVisitMagnitudeMapping>(
                 ccdImageList, focalPlaneBBox, visitOrder);
     }

include/lsst/jointcal/MeasuredStar.h

@@ -41,8 +41,8 @@ class MeasuredStar : public BaseStar {
               _instMag(0.),
               _instMagErr(0.) {}
 
-    /// No move, allow copy constructor: we may copy the fitted StarLists when associating and matching
-    /// catalogs, otherwise Stars should be managed by shared_ptr only.
+    /// No move, allow copy constructor: we may copy the fitted StarLists when associating and
+    /// matching catalogs, otherwise Stars should be managed by shared_ptr only.
     MeasuredStar(MeasuredStar const &) = default;
     MeasuredStar(MeasuredStar &&) = delete;
     MeasuredStar &operator=(MeasuredStar const &) = delete;

include/lsst/jointcal/PhotometryFit.h

@@ -25,8 +25,6 @@ class PhotometryFit : public FitterBase {
      *
      * @param associations The associations catalog to use in the fitter.
      * @param photometryModel The model to build the fitter for.
-     * @param fluxError The systematic error pedestal to apply to measured instFlux errors
-     *                  (as a percent of instFlux).
      */
     PhotometryFit(std::shared_ptr<Associations> associations,
                   std::shared_ptr<PhotometryModel> photometryModel)
@@ -68,6 +66,9 @@ class PhotometryFit : public FitterBase {
     /// @copydoc FitterBase::saveChi2RefContributions
     void saveChi2RefContributions(std::string const &baseName) const override;
 
+    /// Return the model being fit.
+    std::shared_ptr<PhotometryModel> getModel() const { return _photometryModel; }
+
 private:
     bool _fittingModel, _fittingFluxes;
     std::shared_ptr<PhotometryModel> _photometryModel;

include/lsst/jointcal/PhotometryModel.h

@@ -2,23 +2,25 @@
 #ifndef LSST_JOINTCAL_PHOTOMETRY_MODEL_H
 #define LSST_JOINTCAL_PHOTOMETRY_MODEL_H
 
+#include <string>
+#include <vector>
+
 #include "lsst/afw/image/PhotoCalib.h"
+#include "lsst/log/Log.h"
 
 #include "lsst/jointcal/CcdImage.h"
 #include "lsst/jointcal/Eigenstuff.h"
 #include "lsst/jointcal/PhotometryMapping.h"
 #include "lsst/jointcal/FittedStar.h"
 #include "lsst/jointcal/MeasuredStar.h"
 #include "lsst/jointcal/RefStar.h"
-#include <string>
-#include <vector>
 
 namespace lsst {
 namespace jointcal {
 
 class PhotometryModel {
 public:
-    PhotometryModel(double errorPedestal_ = 0) : errorPedestal(errorPedestal_) {}
+    PhotometryModel(LOG_LOGGER log, double errorPedestal_ = 0) : _log(log), errorPedestal(errorPedestal_) {}
 
     /**
      * Assign indices in the full matrix to the parameters being fit in the mappings, starting at firstIndex.
@@ -134,6 +136,24 @@ class PhotometryModel {
     /// Dump the contents of the transfos, for debugging.
     virtual void dump(std::ostream &stream = std::cout) const = 0;
 
+    /**
+     * Return true if this is a "reasonable" model.
+     *
+     * A valid photometry model is positive within each sensor's bounding box.
+     *
+     * @param ccdImageList The ccdImages to test the model validity on.
+     * @return True if the model is valid on all ccdImages.
+     */
+    bool validate(CcdImageList const &ccdImageList) const;
+
+    /**
+     * Check that the model is positive on the ccdImage bbox.
+     *
+     * @param ccdImage The ccdImage to test.
+     * @return True if the image is positive on a sampling of points of the ccdImage bbox.
+     */
+    bool checkPositiveOnBBox(CcdImage const &ccdImage) const;
+
     friend std::ostream &operator<<(std::ostream &s, PhotometryModel const &model) {
         model.dump(s);
         return s;
@@ -163,6 +183,9 @@ class PhotometryModel {
     /// Return a pointer to the mapping associated with this ccdImage.
     virtual PhotometryMappingBase *findMapping(CcdImage const &ccdImage) const = 0;
 
+    /// lsst.logging instance, to be created by a subclass so that messages have consistent name.
+    LOG_LOGGER _log;
+
     // Pedestal on flux/magnitude error (percent of flux or delta magnitude)
     double errorPedestal;
 };

include/lsst/jointcal/PhotometryTransfo.h

@@ -206,7 +206,8 @@ class PhotometryTransfoChebyshev : public PhotometryTransfo {
     /**
      * Create a Chebyshev transfo with the specified coefficients.
      *
-     * The polynomial order is determined from the number of coefficients.
+     * The polynomial order is determined from the number of coefficients, taking only the
+     * anti-diagonal upper triangle portion of the passed-in coefficients
      *
      * @param      coefficients  The polynomial coefficients.
      * @param[in]  bbox          The bounding box it is valid within, to rescale it to [-1,1].
@@ -291,6 +292,9 @@ class FluxTransfoChebyshev : public PhotometryTransfoChebyshev {
     }
 };
 
+/**
+ * nth-order 2d Chebyshev photometry transfo, plus the input flux.
+ */
 class MagnitudeTransfoChebyshev : public PhotometryTransfoChebyshev {
 public:
     MagnitudeTransfoChebyshev(size_t order, afw::geom::Box2D const &bbox)

include/lsst/jointcal/SimplePhotometryModel.h

@@ -19,8 +19,8 @@ class Point;
 //! Photometric response model which has a single photometric factor per CcdImage.
 class SimplePhotometryModel : public PhotometryModel {
 public:
-    SimplePhotometryModel(CcdImageList const &ccdImageList, double errorPedestal_ = 0)
-            : PhotometryModel(errorPedestal_) {
+    SimplePhotometryModel(CcdImageList const &ccdImageList, LOG_LOGGER log, double errorPedestal_ = 0)
+            : PhotometryModel(log, errorPedestal_) {
         _myMap.reserve(ccdImageList.size());
     }
 

python/lsst/jointcal/astrometryModels.cc

@@ -50,6 +50,7 @@ void declareAstrometryModel(py::module &mod) {
     cls.def("getSky2TP", &AstrometryModel::getSky2TP);
     cls.def("makeSkyWcs", &AstrometryModel::makeSkyWcs);
     cls.def("getTotalParameters", &AstrometryModel::getTotalParameters);
+    cls.def("validate", &AstrometryModel::validate);
 }
 
 void declareSimpleAstrometryModel(py::module &mod) {

python/lsst/jointcal/fitter.cc

@@ -51,13 +51,17 @@ void declareAstrometryFit(py::module &mod) {
 
     cls.def(py::init<std::shared_ptr<Associations>, std::shared_ptr<AstrometryModel>, double>(),
             "associations"_a, "astrometryModel"_a, "posError"_a);
+
+    cls.def("getModel", &AstrometryFit::getModel, py::return_value_policy::reference_internal);
 }
 
 void declarePhotometryFit(py::module &mod) {
     py::class_<PhotometryFit, std::shared_ptr<PhotometryFit>, FitterBase> cls(mod, "PhotometryFit");
 
     cls.def(py::init<std::shared_ptr<Associations>, std::shared_ptr<PhotometryModel>>(), "associations"_a,
             "photometryModel"_a);
+
+    cls.def("getModel", &PhotometryFit::getModel, py::return_value_policy::reference_internal);
 }
 
 PYBIND11_MODULE(fitter, mod) {

python/lsst/jointcal/jointcal.py

@@ -585,6 +585,17 @@ def _check_star_lists(self, associations, name):
         if associations.refStarListSize() == 0:
             raise RuntimeError('No stars in the {} reference star list!'.format(name))
 
+    def _logChi2AndValidate(self, associations, fit, model, chi2Label="Model"):
+        """Compute chi2, log it, validate the model, and return chi2."""
+        chi2 = fit.computeChi2()
+        self.log.info("%s %s", chi2Label, chi2)
+        self._check_stars(associations)
+        if not np.isfinite(chi2.chi2):
+            raise FloatingPointError('%s chi2 is invalid: %s', chi2Label, chi2)
+        if not model.validate(associations.getCcdImageList()):
+            raise ValueError("Model is not valid: check log messages for warnings.")
+        return chi2
+
     def _fit_photometry(self, associations, dataName=None):
         """
         Fit the photometric data.
@@ -632,37 +643,32 @@ def _fit_photometry(self, associations, dataName=None):
             doLineSearch = False  # purely linear in model parameters, so no line search needed
 
         fit = lsst.jointcal.PhotometryFit(associations, model)
-        chi2 = fit.computeChi2()
+        self._logChi2AndValidate(associations, fit, model, "Initialized")
+
         # TODO DM-12446: turn this into a "butler save" somehow.
         # Save reference and measurement chi2 contributions for this data
         if self.config.writeChi2ContributionFiles:
             baseName = "photometry_initial_chi2-{}.csv".format(dataName)
             fit.saveChi2Contributions(baseName)
 
-        if not np.isfinite(chi2.chi2):
-            raise FloatingPointError('Initial chi2 is invalid: %s'%chi2)
-        self.log.info("Initialized: %s", str(chi2))
         # The constrained model needs the visit transfo fit first; the chip
         # transfo is initialized from the singleFrame PhotoCalib, so it's close.
         dumpMatrixFile = "photometry_preinit" if self.config.writeInitMatrix else ""
         if self.config.photometryModel.startswith("constrained"):
             # no line search: should be purely (or nearly) linear,
             # and we want a large step size to initialize with.
             fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)
-            chi2 = fit.computeChi2()
-            self.log.info(str(chi2))
+            self._logChi2AndValidate(associations, fit, model)
             dumpMatrixFile = ""  # so we don't redo the output on the next step
+
         fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
-        chi2 = fit.computeChi2()
-        self.log.info(str(chi2))
+        self._logChi2AndValidate(associations, fit, model)
+
         fit.minimize("Fluxes")  # no line search: always purely linear.
-        chi2 = fit.computeChi2()
-        self.log.info(str(chi2))
+        self._logChi2AndValidate(associations, fit, model)
+
         fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
-        chi2 = fit.computeChi2()
-        if not np.isfinite(chi2.chi2):
-            raise FloatingPointError('Pre-iteration chi2 is invalid: %s'%chi2)
-        self.log.info("Fit prepared with %s", str(chi2))
+        self._logChi2AndValidate(associations, fit, model, "Fit prepared")
 
         model.freezeErrorTransform()
         self.log.debug("Photometry error scales are frozen.")
@@ -725,36 +731,30 @@ def _fit_astrometry(self, associations, dataName=None):
                                                         order=self.config.astrometrySimpleOrder)
 
         fit = lsst.jointcal.AstrometryFit(associations, model, self.config.positionErrorPedestal)
-        chi2 = fit.computeChi2()
+        self._logChi2AndValidate(associations, fit, model, "Initial")
+
         # TODO DM-12446: turn this into a "butler save" somehow.
         # Save reference and measurement chi2 contributions for this data
         if self.config.writeChi2ContributionFiles:
             baseName = "astrometry_initial_chi2-{}.csv".format(dataName)
             fit.saveChi2Contributions(baseName)
 
-        if not np.isfinite(chi2.chi2):
-            raise FloatingPointError('Initial chi2 is invalid: %s'%chi2)
-        self.log.info("Initialized: %s", str(chi2))
         dumpMatrixFile = "astrometry_preinit" if self.config.writeInitMatrix else ""
         # The constrained model needs the visit transfo fit first; the chip
         # transfo is initialized from the detector's cameraGeom, so it's close.
         if self.config.astrometryModel == "constrained":
             fit.minimize("DistortionsVisit", dumpMatrixFile=dumpMatrixFile)
-            chi2 = fit.computeChi2()
-            self.log.info(str(chi2))
+            self._logChi2AndValidate(associations, fit, model)
             dumpMatrixFile = ""  # so we don't redo the output on the next step
+
         fit.minimize("Distortions", dumpMatrixFile=dumpMatrixFile)
-        chi2 = fit.computeChi2()
-        self.log.info(str(chi2))
+        self._logChi2AndValidate(associations, fit, model)
+
         fit.minimize("Positions")
-        chi2 = fit.computeChi2()
-        self.log.info(str(chi2))
+        self._logChi2AndValidate(associations, fit, model)
+
         fit.minimize("Distortions Positions")
-        chi2 = fit.computeChi2()
-        self.log.info(str(chi2))
-        if not np.isfinite(chi2.chi2):
-            raise FloatingPointError('Pre-iteration chi2 is invalid: %s'%chi2)
-        self.log.info("Fit prepared with %s", str(chi2))
+        self._logChi2AndValidate(associations, fit, model, "Fit prepared")
 
         chi2 = self._iterate_fit(associations,
                                  fit,
@@ -831,17 +831,15 @@ def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit,
                                      doLineSearch=doLineSearch,
                                      dumpMatrixFile=dumpMatrixFile)
             dumpMatrixFile = ""  # clear it so we don't write the matrix again.
-            chi2 = fitter.computeChi2()
-            self._check_stars(associations)
-            self.log.info(str(chi2))
+            chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel())
+
             if result == MinimizeResult.Converged:
                 if doRankUpdate:
                     self.log.debug("fit has converged - no more outliers - redo minimization "
                                    "one more time in case we have lost accuracy in rank update.")
                     # Redo minimization one more time in case we have lost accuracy in rank update
                     result = fitter.minimize(whatToFit, self.config.outlierRejectSigma)
-                chi2 = fitter.computeChi2()
-                self.log.info("Fit completed with: %s", str(chi2))
+                    chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(), "Fit completed")
 
                 # log a message for a large final chi2, TODO: DM-15247 for something better
                 if chi2.chi2/chi2.ndof >= 4.0:

python/lsst/jointcal/photometryModels.cc

@@ -58,6 +58,9 @@ void declarePhotometryModel(py::module &mod) {
     cls.def("getRefError", &PhotometryModel::getRefError);
     cls.def("computeRefResidual", &PhotometryModel::computeRefResidual);
 
+    cls.def("checkPositiveOnBBox", &PhotometryModel::checkPositiveOnBBox);
+    cls.def("validate", &PhotometryModel::validate);
+
     cls.def("getMappingIndices", &PhotometryModel::getMappingIndices);
     cls.def("computeParameterDerivatives",
             [](PhotometryModel const &self, MeasuredStar const &star, CcdImage const &ccdImage) {