Implement photometry ModelVisit/ModelChip option, with tests.

refactor test_photometryMapping.py so that I can compute the chip and visit
derivatives separately.

The cfht ConstrainedPhotometryModel chi2 changed very slightly, which I
can believe, as it starts with a different initial step.

include/lsst/jointcal/ConstrainedPhotometryModel.h

@@ -77,11 +77,21 @@ class ConstrainedPhotometryModel : public PhotometryModel {
 private:
     PhotometryMappingBase *findMapping(CcdImage const &ccdImage) const override;
 
+    // Which components of the model are we fitting currently?
+    bool _fittingChips;
+    bool _fittingVisits;
+
+    /* The per-ccdImage transforms, each of which is a composition of a chip and visit transform.
+     * Not all pairs of _visitMap[visit] and _chipMap[chip] are guaranteed to have an entry in
+     * _chipVisitMap (for example, one ccd in one visit might fail to produce a catalog).
+     */
     typedef std::unordered_map<CcdImageKey, std::unique_ptr<ChipVisitPhotometryMapping>> MapType;
-    MapType _myMap;
+    MapType _chipVisitMap;
 
+    // The per-visit transforms that go into _chipVisitMap.
     typedef std::map<VisitIdType, std::shared_ptr<PhotometryMapping>> VisitMapType;
     VisitMapType _visitMap;
+    // The per-chip transforms that go into _chipVisitMap.
     typedef std::map<CcdIdType, std::shared_ptr<PhotometryMapping>> ChipMapType;
     ChipMapType _chipMap;
 };

include/lsst/jointcal/PhotometryMapping.h

@@ -72,14 +72,6 @@ class PhotometryMappingBase {
     virtual void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
                                              Eigen::Ref<Eigen::VectorXd> derivatives) const = 0;
 
-    /**
-     * Offset the transfo parameters by delta.
-     *
-     * @param[in]   delta vector to offset transfo parameters. Same ordering as derivatives in
-     *              computeParameterDerivatives();
-     */
-    virtual void offsetParams(Eigen::VectorXd const &delta) = 0;
-
     /// Make this mapping's parameters fixed (i.e. not varied during fitting).
     void setFixed(bool _fixed) { fixed = _fixed; }
     bool isFixed() { return fixed; }
@@ -155,8 +147,13 @@ class PhotometryMapping : public PhotometryMappingBase {
         }
     }
 
-    /// @copydoc PhotometryMappingBase::offsetParams
-    void offsetParams(Eigen::VectorXd const &delta) override { _transfo->offsetParams(delta); }
+    /**
+     * Offset the transfo parameters by delta.
+     *
+     * @param[in]   delta vector to offset transfo parameters. Same ordering as derivatives in
+     *              computeParameterDerivatives();
+     */
+    void offsetParams(Eigen::VectorXd const &delta) { _transfo->offsetParams(delta); }
 
     /// @copydoc PhotometryMappingBase::getParameters
     Eigen::VectorXd getParameters() override { return _transfo->getParameters(); }
@@ -194,11 +191,13 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     ChipVisitPhotometryMapping(std::shared_ptr<PhotometryMapping> chipMapping,
                                std::shared_ptr<PhotometryMapping> visitMapping)
             : PhotometryMappingBase(),
+              _nParChips(0),
+              _nParVisits(0),
               _chipMapping(std::move(chipMapping)),
               _visitMapping(std::move(visitMapping)) {}
 
     /// @copydoc PhotometryMappingBase::getNpar
-    unsigned getNpar() const override { return _chipMapping->getNpar() + _visitMapping->getNpar(); }
+    unsigned getNpar() const override { return _nParChips + _nParVisits; }
 
     /// @copydoc PhotometryMappingBase::transform
     double transform(MeasuredStar const &measuredStar, double instFlux) const override {
@@ -224,12 +223,6 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
                                      Eigen::Ref<Eigen::VectorXd> derivatives) const override;
 
-    /// @copydoc PhotometryMappingBase::offsetParams
-    void offsetParams(Eigen::VectorXd const &delta) override {
-        _chipMapping->offsetParams(delta.segment(0, _chipMapping->getNpar()));
-        _visitMapping->offsetParams(delta.segment(_chipMapping->getNpar(), _visitMapping->getNpar()));
-    }
-
     /// @copydoc PhotometryMappingBase::getParameters
     Eigen::VectorXd getParameters() override {
         Eigen::VectorXd joined(getNpar());
@@ -240,6 +233,17 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     /// @copydoc PhotometryMappingBase::getMappingIndices
     void getMappingIndices(std::vector<unsigned> &indices) const override;
 
+    /**
+     * Set whether to fit chips or visits.
+     *
+     * This must be called before anything that depends on knowing the number of parameters in the fit,
+     * such as offsetParams(), getParameters(), or computeParameterDerivatives().
+     *
+     * @param fittingChips Fit the chip transform.
+     * @param fittingVisits Fit the visit transform.
+     */
+    void setWhatToFit(bool const fittingChips, bool const fittingVisits);
+
     /// @copydoc PhotometryMappingBase::dump
     void dump(std::ostream &stream = std::cout) const override {
         stream << "index: " << index << " chipMapping: ";
@@ -252,6 +256,9 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     std::shared_ptr<PhotometryMapping> getVisitMapping() const { return _visitMapping; }
 
 private:
+    // These are either transfo.getNpar() or 0, depending on whether we are fitting that component or not.
+    unsigned _nParChips, _nParVisits;
+
     // the actual transformation to be fit
     std::shared_ptr<PhotometryMapping> _chipMapping;
     std::shared_ptr<PhotometryMapping> _visitMapping;

include/lsst/jointcal/PhotometryModel.h

@@ -21,7 +21,7 @@ class MeasuredStar;
 class PhotometryModel {
 public:
     /**
-     * Assign indices to parameters involved in mappings, starting at firstIndex.
+     * Assign indices in the full matrix to the parameters being fit in the mappings, starting at firstIndex.
      *
      * @param[in]  whatToFit   String containing parameters to fit.
      * @param[in]  firstIndex  Index to start assigning at.

python/lsst/jointcal/jointcal.py

@@ -621,7 +621,6 @@ def _fit_photometry(self, associations, dataName=None):
         # transfo is initialized from the singleFrame PhotoCalib, so it's close.
         dumpMatrixFile = "photometry_preinit" if self.config.writeInitMatrix else ""
         if self.config.photometryModel == "constrained":
-            # TODO: (related to DM-8046): implement Visit/Chip choice
             # no line search: should be purely (or nearly) linear,
             # and we want a large step size to initialize with.
             fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)

python/lsst/jointcal/photometryMappings.cc

@@ -51,7 +51,6 @@ void declarePhotometryMappingBase(py::module &mod) {
                 return derivatives;
             });
 
-    cls.def("offsetParams", &PhotometryMappingBase::offsetParams);
     cls.def("getParameters", &PhotometryMappingBase::getParameters);
 
     cls.def("getMappingIndices", [](PhotometryMappingBase const &self) {
@@ -69,6 +68,7 @@ void declarePhotometryMapping(py::module &mod) {
             mod, "PhotometryMapping");
     cls.def(py::init<std::shared_ptr<PhotometryTransfo>>(), "transfo"_a);
 
+    cls.def("offsetParams", &PhotometryMapping::offsetParams);
     cls.def("getTransfo", &PhotometryMapping::getTransfo);
     cls.def("getTransfoErrors", &PhotometryMapping::getTransfoErrors);
 }
@@ -80,6 +80,8 @@ void declareChipVisitPhotometryMapping(py::module &mod) {
     cls.def(py::init<std::shared_ptr<PhotometryMapping>, std::shared_ptr<PhotometryMapping>>(),
             "chipMapping"_a, "visitMapping"_a);
 
+    cls.def("setWhatToFit", &ChipVisitPhotometryMapping::setWhatToFit);
+
     cls.def("getChipMapping", &ChipVisitPhotometryMapping::getChipMapping);
     cls.def("getVisitMapping", &ChipVisitPhotometryMapping::getVisitMapping);
 }

src/ConstrainedPhotometryModel.cc

@@ -20,8 +20,8 @@ namespace lsst {
 namespace jointcal {
 
 ConstrainedPhotometryModel::ConstrainedPhotometryModel(CcdImageList const &ccdImageList,
-                                                       afw::geom::Box2D const &focalPlaneBBox,
-                                                       int visitOrder) {
+                                                       afw::geom::Box2D const &focalPlaneBBox, int visitOrder)
+        : _fittingChips(false), _fittingVisits(false) {
     // keep track of which chip we want to constrain (the one closest to the middle of the focal plane)
     double minRadius2 = std::numeric_limits<double>::infinity();
     CcdIdType constrainedChip = -1;
@@ -45,64 +45,86 @@ ConstrainedPhotometryModel::ConstrainedPhotometryModel(CcdImageList const &ccdIm
             // Use the single-frame processing calibration from the PhotoCalib as the default.
             auto chipTransfo =
                     std::make_shared<PhotometryTransfoSpatiallyInvariant>(photoCalib->getCalibrationMean());
-            _chipMap[chip] =
-                    std::unique_ptr<PhotometryMapping>(new PhotometryMapping(std::move(chipTransfo)));
+            _chipMap[chip] = std::make_unique<PhotometryMapping>(std::move(chipTransfo));
         }
         // If the visit is not in the map, add it, otherwise continue.
         if (visitPair == _visitMap.end()) {
             auto visitTransfo = std::make_shared<PhotometryTransfoChebyshev>(visitOrder, focalPlaneBBox);
-            _visitMap[visit] =
-                    std::unique_ptr<PhotometryMapping>(new PhotometryMapping(std::move(visitTransfo)));
+            _visitMap[visit] = std::make_unique<PhotometryMapping>(std::move(visitTransfo));
         }
     }
 
     // Fix one chip mapping, to remove the degeneracy from the system.
     _chipMap.at(constrainedChip)->setFixed(true);
 
     // Now create the ccdImage mappings, which are combinations of the chip/visit mappings above.
-    _myMap.reserve(ccdImageList.size());  // we know how big it will be, so pre-allocate space.
+    _chipVisitMap.reserve(
+            ccdImageList.size());  // we know how chipVisitbig it will be, so pre-allocate space.
     for (auto const &ccdImage : ccdImageList) {
         auto visit = ccdImage->getVisit();
         auto chip = ccdImage->getCcdId();
-        _myMap.emplace(ccdImage->getHashKey(),
-                       std::unique_ptr<ChipVisitPhotometryMapping>(
-                               new ChipVisitPhotometryMapping(_chipMap[chip], _visitMap[visit])));
+        _chipVisitMap.emplace(ccdImage->getHashKey(),
+                              std::make_unique<ChipVisitPhotometryMapping>(_chipMap[chip], _visitMap[visit]));
     }
     LOGLS_INFO(_log, "Got " << _chipMap.size() << " chip mappings and " << _visitMap.size()
                             << " visit mappings; holding chip " << constrainedChip << " fixed ("
                             << getTotalParameters() << " total parameters).");
-    LOGLS_DEBUG(_log, "CcdImage map has " << _myMap.size() << " mappings, with " << _myMap.bucket_count()
-                                          << " buckets and a load factor of " << _myMap.load_factor());
+    LOGLS_DEBUG(_log, "CcdImage map has " << _chipVisitMap.size() << " mappings, with "
+                                          << _chipVisitMap.bucket_count() << " buckets and a load factor of "
+                                          << _chipVisitMap.load_factor());
 }
 
 unsigned ConstrainedPhotometryModel::assignIndices(std::string const &whatToFit, unsigned firstIndex) {
-    // TODO DM-8046: currently ignoring whatToFit: eventually implement configurability.
     unsigned index = firstIndex;
-    for (auto &i : _chipMap) {
-        auto mapping = i.second.get();
-        // Don't assign indices for fixed parameters.
-        if (mapping->isFixed()) continue;
-        mapping->setIndex(index);
-        index += mapping->getNpar();
-    }
-    for (auto &i : _visitMap) {
-        auto mapping = i.second.get();
-        mapping->setIndex(index);
-        index += mapping->getNpar();
+    if (whatToFit.find("Model") == std::string::npos) {
+        LOGLS_WARN(_log, "assignIndices was called and Model is *not* in whatToFit");
+        return index;
+    }
+
+    // If we got here, "Model" is definitely in whatToFit.
+    _fittingChips = (whatToFit.find("ModelChip") != std::string::npos);
+    _fittingVisits = (whatToFit.find("ModelVisit") != std::string::npos);
+    // If nothing more than "Model" is specified, it means fit everything.
+    if ((!_fittingChips) && (!_fittingVisits)) {
+        _fittingChips = _fittingVisits = true;
+    }
+
+    if (_fittingChips) {
+        for (auto &idMapping : _chipMap) {
+            auto mapping = idMapping.second.get();
+            // Don't assign indices for fixed parameters.
+            if (mapping->isFixed()) continue;
+            mapping->setIndex(index);
+            index += mapping->getNpar();
+        }
+    }
+    if (_fittingVisits) {
+        for (auto &idMapping : _visitMap) {
+            auto mapping = idMapping.second.get();
+            mapping->setIndex(index);
+            index += mapping->getNpar();
+        }
+    }
+    for (auto &idMapping : _chipVisitMap) {
+        idMapping.second->setWhatToFit(_fittingChips, _fittingVisits);
     }
     return index;
 }
 
 void ConstrainedPhotometryModel::offsetParams(Eigen::VectorXd const &delta) {
-    for (auto &i : _chipMap) {
-        auto mapping = i.second.get();
-        // Don't offset indices for fixed parameters.
-        if (mapping->isFixed()) continue;
-        mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
+    if (_fittingChips) {
+        for (auto &idMapping : _chipMap) {
+            auto mapping = idMapping.second.get();
+            // Don't offset indices for fixed parameters.
+            if (mapping->isFixed()) continue;
+            mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
+        }
     }
-    for (auto &i : _visitMap) {
-        auto mapping = i.second.get();
-        mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
+    if (_fittingVisits) {
+        for (auto &idMapping : _visitMap) {
+            auto mapping = idMapping.second.get();
+            mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
+        }
     }
 }
 
@@ -119,11 +141,11 @@ double ConstrainedPhotometryModel::transformError(CcdImage const &ccdImage, Meas
 }
 
 void ConstrainedPhotometryModel::freezeErrorTransform() {
-    for (auto &i : _chipMap) {
-        i.second.get()->freezeErrorTransform();
+    for (auto &idMapping : _chipMap) {
+        idMapping.second.get()->freezeErrorTransform();
     }
-    for (auto &i : _visitMap) {
-        i.second.get()->freezeErrorTransform();
+    for (auto &idMapping : _visitMap) {
+        idMapping.second.get()->freezeErrorTransform();
     }
 }
 
@@ -135,11 +157,11 @@ void ConstrainedPhotometryModel::getMappingIndices(CcdImage const &ccdImage,
 
 int ConstrainedPhotometryModel::getTotalParameters() const {
     int total = 0;
-    for (auto &i : _chipMap) {
-        total += i.second->getNpar();
+    for (auto &idMapping : _chipMap) {
+        total += idMapping.second->getNpar();
     }
-    for (auto &i : _visitMap) {
-        total += i.second->getNpar();
+    for (auto &idMapping : _visitMap) {
+        total += idMapping.second->getNpar();
     }
     return total;
 }
@@ -211,23 +233,23 @@ std::shared_ptr<afw::image::PhotoCalib> ConstrainedPhotometryModel::toPhotoCalib
 }
 
 void ConstrainedPhotometryModel::dump(std::ostream &stream) const {
-    for (auto &i : _chipMap) {
-        i.second->dump(stream);
+    for (auto &idMapping : _chipMap) {
+        idMapping.second->dump(stream);
         stream << std::endl;
     }
     stream << std::endl;
-    for (auto &i : _visitMap) {
-        i.second->dump(stream);
+    for (auto &idMapping : _visitMap) {
+        idMapping.second->dump(stream);
         stream << std::endl;
     }
 }
 
 PhotometryMappingBase *ConstrainedPhotometryModel::findMapping(CcdImage const &ccdImage) const {
-    auto i = _myMap.find(ccdImage.getHashKey());
-    if (i == _myMap.end())
+    auto idMapping = _chipVisitMap.find(ccdImage.getHashKey());
+    if (idMapping == _chipVisitMap.end())
         throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                           "ConstrainedPhotometryModel cannot find CcdImage " + ccdImage.getName());
-    return i->second.get();
+    return idMapping->second.get();
 }
 
 }  // namespace jointcal