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.

include/lsst/jointcal/ConstrainedPhotometryModel.h

@@ -77,6 +77,10 @@ class ConstrainedPhotometryModel : public PhotometryModel {
 private:
     PhotometryMappingBase *findMapping(CcdImage const &ccdImage) const override;
 
+    // Which components of the model are we fitting?
+    bool _fittingChips;
+    bool _fittingVisits;
+
     typedef std::unordered_map<CcdImageKey, std::unique_ptr<ChipVisitPhotometryMapping>> MapType;
     MapType _myMap;
 

include/lsst/jointcal/PhotometryMapping.h

@@ -198,7 +198,7 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
               _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,10 +224,19 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
                                      Eigen::Ref<Eigen::VectorXd> derivatives) const override;
 
-    /// @copydoc PhotometryMappingBase::offsetParams
+    /**
+     * @copydoc PhotometryMappingBase::offsetParams
+     *
+     * NOTE: this is not used during fitting (the model manages the two mappings separately),
+     * but it can be helpful for debugging.
+     */
     void offsetParams(Eigen::VectorXd const &delta) override {
-        _chipMapping->offsetParams(delta.segment(0, _chipMapping->getNpar()));
-        _visitMapping->offsetParams(delta.segment(_chipMapping->getNpar(), _visitMapping->getNpar()));
+        if (_nParChips > 0) {
+            _chipMapping->offsetParams(delta.segment(0, _nParChips));
+        }
+        if (_nParVisits > 0) {
+            _visitMapping->offsetParams(delta.segment(_nParChips, _nParVisits));
+        }
     }
 
     /// @copydoc PhotometryMappingBase::getParameters
@@ -240,6 +249,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 +272,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;

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

@@ -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

@@ -76,33 +76,56 @@ ConstrainedPhotometryModel::ConstrainedPhotometryModel(CcdImageList const &ccdIm
 }
 
 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();
+    if (whatToFit.find("Model") == std::string::npos) {
+        LOGLS_ERROR(_log, "assignIndices was called and Model is *not* in whatToFit");
+        return 0;
     }
-    for (auto &i : _visitMap) {
-        auto mapping = i.second.get();
-        mapping->setIndex(index);
-        index += mapping->getNpar();
+
+    // 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 &i : _chipMap) {
+            auto mapping = i.second.get();
+            // Don't assign indices for fixed parameters.
+            if (mapping->isFixed()) continue;
+            mapping->setIndex(index);
+            index += mapping->getNpar();
+        }
+    }
+    if (_fittingVisits) {
+        for (auto &i : _visitMap) {
+            auto mapping = i.second.get();
+            mapping->setIndex(index);
+            index += mapping->getNpar();
+        }
+    }
+    for (auto &i : _myMap) {
+        i.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 &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()));
+        }
     }
-    for (auto &i : _visitMap) {
-        auto mapping = i.second.get();
-        mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
+    if (_fittingVisits) {
+        for (auto &i : _visitMap) {
+            auto mapping = i.second.get();
+            mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
+        }
     }
 }
 

src/PhotometryMapping.cc

@@ -23,30 +23,49 @@ void ChipVisitPhotometryMapping::computeParameterDerivatives(MeasuredStar const
 
     // NOTE: the chip derivatives start at 0, and the visit derivatives start at the last chip derivative.
     // This is independent of the full-fit indices.
-    Eigen::Ref<Eigen::VectorXd> chipBlock = derivatives.segment(0, _chipMapping->getNpar());
-    Eigen::Ref<Eigen::VectorXd> visitBlock =
-            derivatives.segment(_chipMapping->getNpar(), _visitMapping->getNpar());
 
     // NOTE: chipBlock is the product of the chip derivatives and the visit transforms, and vice versa.
     // NOTE: See DMTN-036 for the math behind this.
-    if (!_chipMapping->isFixed()) {
+    if (_nParChips > 0 && !_chipMapping->isFixed()) {
+        Eigen::Ref<Eigen::VectorXd> chipBlock = derivatives.segment(0, _nParChips);
         _chipMapping->getTransfo()->computeParameterDerivatives(measuredStar.x, measuredStar.y, instFlux,
                                                                 chipBlock);
         chipBlock *= visitScale;
     }
-    _visitMapping->getTransfo()->computeParameterDerivatives(measuredStar.getXFocal(),
-                                                             measuredStar.getYFocal(), instFlux, visitBlock);
-    visitBlock *= chipScale;
+    if (_nParVisits > 0) {
+        Eigen::Ref<Eigen::VectorXd> visitBlock = derivatives.segment(_nParChips, _nParVisits);
+        _visitMapping->getTransfo()->computeParameterDerivatives(
+                measuredStar.getXFocal(), measuredStar.getYFocal(), instFlux, visitBlock);
+        visitBlock *= chipScale;
+    }
 }
 
 void ChipVisitPhotometryMapping::getMappingIndices(std::vector<unsigned> &indices) const {
     if (indices.size() < getNpar()) indices.resize(getNpar());
-    _chipMapping->getMappingIndices(indices);
-    // TODO DM-12169: there is probably a better way to feed a subpart of a std::vector (a view or iterators?)
-    std::vector<unsigned> tempIndices(_visitMapping->getNpar());
-    _visitMapping->getMappingIndices(tempIndices);
-    for (unsigned k = 0; k < _visitMapping->getNpar(); ++k) {
-        indices.at(k + _chipMapping->getNpar()) = tempIndices.at(k);
+    if (_nParChips > 0) {
+        _chipMapping->getMappingIndices(indices);
+    }
+    if (_nParVisits > 0) {
+        // TODO DM-12169: there is probably a better way to feed a subpart of a std::vector
+        // (maybe a view or iterators?)
+        std::vector<unsigned> tempIndices(_visitMapping->getNpar());
+        _visitMapping->getMappingIndices(tempIndices);
+        for (unsigned k = 0; k < _visitMapping->getNpar(); ++k) {
+            indices.at(k + _nParChips) = tempIndices.at(k);
+        }
+    }
+}
+
+void ChipVisitPhotometryMapping::setWhatToFit(bool const fittingChips, bool const fittingVisits) {
+    if (fittingChips) {
+        _nParChips = _chipMapping->getNpar();
+    } else {
+        _nParChips = 0;
+    }
+    if (fittingVisits) {
+        _nParVisits = _visitMapping->getNpar();
+    } else {
+        _nParVisits = 0;
     }
 }
 

tests/test_jointcal_cfht.py

@@ -178,7 +178,7 @@ def test_jointcalTask_2_visits_constrainedPhotometry_lineSearch(self):
         self.config.allowLineSearch = True
 
         # Only this value should differ from the metrics defined in setup above.
-        metrics['photometry_final_chi2'] = 2642.54
+        metrics['photometry_final_chi2'] = 2642.47
 
         self._testJointcalTask(2, None, None, pa1, metrics=metrics)
 

tests/test_photometryMapping.py

@@ -94,6 +94,7 @@ def setUp(self):
         visitMapping.setIndex(self.visitIndex)
         self.mappingInvariants = lsst.jointcal.photometryMappings.ChipVisitPhotometryMapping(chipMapping,
                                                                                              visitMapping)
+        self.mappingInvariants.setWhatToFit(True, True)  # default to fitting both
 
         # Have to make a new chipMapping, as it stores shared_ptr to the transfo.
         chipTransfo = lsst.jointcal.photometryTransfo.PhotometryTransfoSpatiallyInvariant(self.chipScale)
@@ -105,6 +106,7 @@ def setUp(self):
         visitMapping2.setIndex(self.visitIndex)
         self.mappingCheby = lsst.jointcal.photometryMappings.ChipVisitPhotometryMapping(chipMapping,
                                                                                         visitMapping2)
+        self.mappingCheby.setWhatToFit(True, True)  # default to fitting both
 
     def test_getNpar(self):
         result = self.mappingInvariants.getNpar()
@@ -164,23 +166,33 @@ def test_offsetParams(self):
         self.mappingCheby.offsetParams(delta)
         self.assertFloatsAlmostEqual(self.mappingCheby.getParameters(), expect)
 
-    def test_computeParameterDerivatives(self):
-        result = self.mappingInvariants.computeParameterDerivatives(self.star0, self.instFlux)
-        expect = np.array([self.instFlux*self.visitScale, self.instFlux*self.chipScale])
-        self.assertFloatsAlmostEqual(result, expect)
-
+    def _computeChebyshevDerivative(self, star):
+        """Return the derivatives w.r.t. the Chebyshev components."""
         cx = (self.bbox.getMinX() + self.bbox.getMaxX())/2.0
         cy = (self.bbox.getMinY() + self.bbox.getMaxY())/2.0
         sx = 2.0 / self.bbox.getWidth()
         sy = 2.0 / self.bbox.getHeight()
-        Tx = np.array([CHEBYSHEV_T[i](sx*(self.star1.getXFocal() - cx))
+        Tx = np.array([CHEBYSHEV_T[i](sx*(star.getXFocal() - cx))
                       for i in range(self.order+1)], dtype=float)
-        Ty = np.array([CHEBYSHEV_T[i](sy*(self.star1.getYFocal() - cy))
+        Ty = np.array([CHEBYSHEV_T[i](sy*(star.getYFocal() - cy))
                       for i in range(self.order+1)], dtype=float)
-        expect = [self.instFlux*self._evaluate_chebyshev(self.star1.getXFocal(), self.star1.getYFocal()), ]
+        expect = []
         for j in range(len(Ty)):
             for i in range(0, self.order-j+1):
                 expect.append(Ty[j]*Tx[i]*self.instFlux*self.chipScale)
+        return expect
+
+    def _computeChipDerivative(self, star):
+        """Return the derivative w.r.t. the chip component."""
+        return self.instFlux*self._evaluate_chebyshev(star.getXFocal(), star.getYFocal())
+
+    def test_computeParameterDerivatives(self):
+        result = self.mappingInvariants.computeParameterDerivatives(self.star0, self.instFlux)
+        expect = np.array([self.instFlux*self.visitScale, self.instFlux*self.chipScale])
+        self.assertFloatsAlmostEqual(result, expect)
+
+        expect = [self._computeChipDerivative(self.star1), ]
+        expect.extend(self._computeChebyshevDerivative(self.star1))
         result = self.mappingCheby.computeParameterDerivatives(self.star1, self.instFlux)
         self.assertFloatsAlmostEqual(np.array(expect), result)
 
@@ -196,6 +208,35 @@ def test_getMappingIndices(self):
         result = self.mappingCheby.getMappingIndices()
         self.assertEqual(result, expect)
 
+    def test_setWhatToFit(self):
+        """Test that mapping methods behave correctly when chip and/or visit
+        fitting is disabled.
+
+        The fit both case (True, True) is tested by all of the above tests.
+        """
+        # Using mappingCheby so getNpar() will distinguish chips (1 param) from visits (3 params).
+
+        # fit nothing means 0 parameters and no indices
+        self.mappingCheby.setWhatToFit(False, False)
+        self.assertEqual(self.mappingCheby.getNpar(), 0)
+        self.assertEqual(self.mappingCheby.getMappingIndices(), [])
+        result = self.mappingCheby.computeParameterDerivatives(self.star1, self.instFlux)
+        self.assertFloatsEqual(result, np.array([]))
+
+        # fit just chips means 1 parameters and self.chipIndex
+        self.mappingCheby.setWhatToFit(True, False)
+        self.assertEqual(self.mappingCheby.getNpar(), 1)
+        self.assertEqual(self.mappingCheby.getMappingIndices(), [self.chipIndex])
+        result = self.mappingCheby.computeParameterDerivatives(self.star1, self.instFlux)
+        self.assertFloatsAlmostEqual(result, [self._computeChipDerivative(self.star1)])
+
+        # fit just visits means 3 parameters (order 1) and self.visitIndex
+        self.mappingCheby.setWhatToFit(False, True)
+        self.assertEqual(self.mappingCheby.getNpar(), 3)
+        self.assertEqual(self.mappingCheby.getMappingIndices(), [1000, 1001, 1002])
+        result = self.mappingCheby.computeParameterDerivatives(self.star1, self.instFlux)
+        self.assertFloatsAlmostEqual(result, [self._computeChebyshevDerivative(self.star1)])
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass

tests/test_photometryModel.py

@@ -106,7 +106,7 @@ def setUp(self):
                                                                                self.focalPlaneBBox,
                                                                                self.visitOrder)
         # have to call this once to let offsetParams work.
-        self.model.assignIndices("", self.firstIndex)
+        self.model.assignIndices("Model", self.firstIndex)
         # tweak to get more than just a constant field for the second ccdImage
         self.delta = np.arange(20, dtype=float)*-0.2 + 1
         # but keep the first ccdImage constant, to help distinguish test failures.
@@ -135,6 +135,40 @@ def test_toPhotoCalib(self):
         self._toPhotoCalib(self.ccdImageList[0])
         self._toPhotoCalib(self.ccdImageList[1])
 
+    def test_assignIndices(self):
+        """Test that the correct number of indices were assigned.
+        Does not check that the internal mappings are assigned the correct
+        indices.
+        """
+        # need at least two sensors to distinguish "Model" from "ModelVisit"
+        # NOTE: createTwoFakeCcdImages() always uses the same two visitIds,
+        # so there will be 2 visits total here.
+        struct1 = lsst.jointcal.testUtils.createTwoFakeCcdImages(100, 100, seed=100, fakeCcdId=12)
+        ccdImageList = struct1.ccdImageList
+        struct2 = lsst.jointcal.testUtils.createTwoFakeCcdImages(100, 100, seed=101, fakeCcdId=13)
+        ccdImageList.extend(struct2.ccdImageList)
+        camera = struct1.camera  # the camera is the same in both structs
+        visitOrder = 3
+        focalPlaneBBox = camera.getFpBBox()
+        model = lsst.jointcal.photometryModels.ConstrainedPhotometryModel(ccdImageList,
+                                                                          focalPlaneBBox,
+                                                                          visitOrder)
+
+        # one polynomial per visit, plus one fitted scale for the second chip.
+        expect = 2 * getNParametersPolynomial(self.visitOrder) + 1
+        index = model.assignIndices("Model", self.firstIndex)
+        self.assertEqual(index, expect)
+
+        # one polynomial per visit
+        expect = 2 * getNParametersPolynomial(self.visitOrder)
+        index = model.assignIndices("ModelVisit", self.firstIndex)
+        self.assertEqual(index, expect)
+
+        # one fitted chip
+        expect = 1
+        index = model.assignIndices("ModelChip", self.firstIndex)
+        self.assertEqual(index, expect)
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass