Merge pull request #96 from lsst/tickets/DM-15011

DM-15011 Add ModelVisit/ModelChip options for photometry

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);
 }

python/lsst/jointcal/testUtils.py

@@ -34,7 +34,7 @@
 import lsst.jointcal.star
 
 
-def createTwoFakeCcdImages(num1=4, num2=4, seed=100):
+def createTwoFakeCcdImages(num1=4, num2=4, seed=100, fakeCcdId=12):
     """Return two fake ccdImages built on CFHT Megacam metadata.
 
     If ``num1 == num2``, the catalogs will align on-sky so each source will
@@ -50,6 +50,10 @@ def createTwoFakeCcdImages(num1=4, num2=4, seed=100):
         a square, to have sqrt(num) centroids on a grid.
     seed : `int`, optional
         Seed value for np.random.
+    fakeCcdId : `int`, optional
+        Sensor identifier to use for both CcdImages. The wcs, bbox, calib, etc.
+        will still be drawn from the CFHT ccd=12 files, as that is the only
+        testdata that is included in this simple test dataset.
 
     Returns
     -------
@@ -68,7 +72,6 @@ def createTwoFakeCcdImages(num1=4, num2=4, seed=100):
 
     visit1 = 849375
     visit2 = 850587
-    ccdId = 12
     instFluxKeyName = "SomeFlux"
 
     # Load or fake the necessary metadata for each CcdImage
@@ -77,21 +80,21 @@ def createTwoFakeCcdImages(num1=4, num2=4, seed=100):
     butler = lsst.daf.persistence.Butler(inputDir)
 
     # so we can access parts of the camera later (e.g. focal plane)
-    camera = butler.get('camera', visit=visit1, ccd=ccdId)
+    camera = butler.get('camera', visit=visit1)
 
-    struct1 = createFakeCcdImage(butler, visit1, ccdId, num1, instFluxKeyName,
-                                 photoCalibMean=100.0, photoCalibErr=1.0)
-    struct2 = createFakeCcdImage(butler, visit2, ccdId, num2, instFluxKeyName,
-                                 photoCalibMean=120.0, photoCalibErr=5.0)
+    struct1 = createFakeCcdImage(butler, visit1, num1, instFluxKeyName,
+                                 photoCalibMean=100.0, photoCalibErr=1.0, fakeCcdId=fakeCcdId)
+    struct2 = createFakeCcdImage(butler, visit2, num2, instFluxKeyName,
+                                 photoCalibMean=120.0, photoCalibErr=5.0, fakeCcdId=fakeCcdId)
 
     return lsst.pipe.base.Struct(camera=camera,
                                  catalogs=[struct1.catalog, struct2.catalog],
                                  ccdImageList=[struct1.ccdImage, struct2.ccdImage],
                                  bbox=struct1.bbox)
 
 
-def createFakeCcdImage(butler, visit, ccdId, num, instFluxKeyName,
-                       photoCalibMean=100.0, photoCalibErr=1.0):
+def createFakeCcdImage(butler, visit, num, instFluxKeyName,
+                       photoCalibMean=100.0, photoCalibErr=1.0, fakeCcdId=12):
     """Create a fake CcdImage by making a fake catalog.
 
     Parameters
@@ -100,8 +103,6 @@ def createFakeCcdImage(butler, visit, ccdId, num, instFluxKeyName,
         Butler to load metadata from.
     visit : `int`
         Visit identifier to build a butler dataId.
-    ccdId : `int`
-        CCD identifier to build a butler dataId.
     num : `int`
         Number of sources to put in the catalogs. Should be
         a square, to have sqrt(num) centroids on a grid.
@@ -111,6 +112,8 @@ def createFakeCcdImage(butler, visit, ccdId, num, instFluxKeyName,
         Value to set for calibrationMean in the created PhotoCalib.
     photoCalibErr : `float`, optional
         Value to set for calibrationErr in the created PhotoCalib.
+    fakeCcdId : `int`, optional
+        Use this as the ccdId in the returned CcdImage.
 
     Returns
     -------
@@ -123,6 +126,8 @@ def createFakeCcdImage(butler, visit, ccdId, num, instFluxKeyName,
            (`lsst.jointcal.CcdImage`).
        - `bbox` : Bounding Box of the image (`lsst.afw.geom.Box2I`).
     """
+    ccdId = 12  # we only have data for ccd=12
+
     dataId = dict(visit=visit, ccd=ccdId)
     skyWcs = butler.get('calexp_wcs', dataId=dataId)
     visitInfo = butler.get('calexp_visitInfo', dataId=dataId)
@@ -133,7 +138,7 @@ def createFakeCcdImage(butler, visit, ccdId, num, instFluxKeyName,
 
     catalog = createFakeCatalog(num, bbox, instFluxKeyName, skyWcs=skyWcs)
     ccdImage = lsst.jointcal.ccdImage.CcdImage(catalog, skyWcs, visitInfo, bbox, filt, photoCalib,
-                                               detector, visit, ccdId, instFluxKeyName)
+                                               detector, visit, fakeCcdId, instFluxKeyName)
 
     return lsst.pipe.base.Struct(catalog=catalog, ccdImage=ccdImage, bbox=bbox)