DM-14786 implement ConstrainedMagnitudeModel

bin.src/plot_photoCalib.py

@@ -120,12 +120,14 @@ def makePhotoCalibImages(visit, butler, step=8, chips=[], tract=None,
 
         calibScaling = 1
         if singleCalib:
-            calib = butler.get('calexp_calib', dataId=dict(visit=int(visit), ccd=int(ccd), tract=0))
+            calib = butler.get('calexp_calib', dataId=dict(visit=int(visit), ccd=int(ccd), tract=tract))
             calibScaling = calib.getFluxMag0()[0]
             if verbose:
                 print('calib ccd %s: %s' % (ccd, calibScaling))
 
-        photoCalib = butler.get('jointcal_photoCalib', dataId=dict(visit=int(visit), ccd=int(ccd), tract=0))
+        photoCalib = butler.get('jointcal_photoCalib', dataId=dict(visit=int(visit),
+                                                                   ccd=int(ccd),
+                                                                   tract=tract))
         if verbose:
             print("photoCalib mean ccd %s: %s" % (ccd, photoCalib.getCalibrationMean()))
         scaled = photoCalib.computeScaledCalibration()

include/lsst/jointcal/AstrometryModel.h

@@ -30,9 +30,13 @@ class AstrometryModel {
     //! assigned index.
     virtual unsigned assignIndices(std::string const &whatToFit, unsigned firstIndex) = 0;
 
-    //! Offset the parameters by the provided amounts.
-    /*! The shifts are applied according to the indices given in
-        AssignIndices. */
+    /**
+     * Offset the parameters by the provided amounts (by -delta).
+     *
+     * The shifts are applied according to the indices given in assignIndices.
+     *
+     * @param[in]  delta  vector of offsets to apply
+     */
     virtual void offsetParams(Eigen::VectorXd const &delta) = 0;
 
     //! The transformation used to project the positions of FittedStars.

include/lsst/jointcal/ConstrainedPhotometryModel.h

@@ -33,7 +33,10 @@ class ConstrainedPhotometryModel : public PhotometryModel {
      * @param[in]  visitOrder    The order of the visit polynomial.
      */
     explicit ConstrainedPhotometryModel(CcdImageList const &ccdImageList,
-                                        afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7);
+                                        afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7)
+            : _fittingChips(false), _fittingVisits(false) {
+        _chipVisitMap.reserve(ccdImageList.size());
+    }
 
     /// No copy or move: there is only ever one instance of a given model (i.e. per ccd+visit)
     ConstrainedPhotometryModel(ConstrainedPhotometryModel const &) = delete;
@@ -47,18 +50,6 @@ class ConstrainedPhotometryModel : public PhotometryModel {
     /// @copydoc PhotometryModel::offsetParams
     void offsetParams(Eigen::VectorXd const &delta) override;
 
-    /// @copydoc PhotometryModel::offsetFittedStar
-    void offsetFittedStar(FittedStar &fittedStar, double delta) const { fittedStar.getFlux() -= delta; }
-
-    /// @copydoc PhotometryModel::computeResidual
-    double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
-
-    /// @copydoc PhotometryModel::transform
-    double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
-
-    /// @copydoc PhotometryModel::transformError
-    double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
-
     /// @copydoc PhotometryModel::freezeErrorTransform
     void freezeErrorTransform() override;
 
@@ -72,27 +63,12 @@ class ConstrainedPhotometryModel : public PhotometryModel {
     void computeParameterDerivatives(MeasuredStar const &measuredStar, CcdImage const &ccdImage,
                                      Eigen::VectorXd &derivatives) const override;
 
-    /// @copydoc PhotometryModel::getRefError
-    double getRefError(RefStar const &refStar) const override { return refStar.getFluxErr(); }
-
-    /// @copydoc PhotometryModel::computeRefResidual
-    double computeRefResidual(FittedStar const &fittedStar, RefStar const &refStar) const override {
-        return fittedStar.getFlux() - refStar.getFlux();
-    };
-
-    /// @copydoc PhotometryModel::toPhotoCalib
-    std::shared_ptr<afw::image::PhotoCalib> toPhotoCalib(CcdImage const &ccdImage) const override;
-
     /// @copydoc PhotometryModel::dump
     void dump(std::ostream &stream = std::cout) const override;
 
-private:
+protected:
     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).
@@ -106,6 +82,102 @@ class ConstrainedPhotometryModel : public PhotometryModel {
     // The per-chip transforms that go into _chipVisitMap.
     typedef std::map<CcdIdType, std::shared_ptr<PhotometryMapping>> ChipMapType;
     ChipMapType _chipMap;
+
+    /**
+     * Initialize the chip, visit, and chipVisit mappings by creating appropriate transfos and mappings.
+     */
+    template <class ChipTransfo, class VisitTransfo, class ChipVisitMapping>
+    void initialize(CcdImageList const &ccdImageList, afw::geom::Box2D const &focalPlaneBBox, int visitOrder);
+
+    /// Return the initial calibration to use from this photoCalib.
+    virtual double initialChipCalibration(std::shared_ptr<afw::image::PhotoCalib const> photoCalib) = 0;
+
+private:
+    // Which components of the model are we fitting currently?
+    bool _fittingChips;
+    bool _fittingVisits;
+};
+
+class ConstrainedFluxModel : public ConstrainedPhotometryModel {
+public:
+    explicit ConstrainedFluxModel(CcdImageList const &ccdImageList, afw::geom::Box2D const &focalPlaneBBox,
+                                  int visitOrder = 7)
+            : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder) {
+        initialize<FluxTransfoSpatiallyInvariant, FluxTransfoChebyshev, ChipVisitFluxMapping>(
+                ccdImageList, focalPlaneBBox, visitOrder);
+    }
+
+    /// @copydoc PhotometryModel::offsetFittedStar
+    void offsetFittedStar(FittedStar &fittedStar, double delta) const override {
+        fittedStar.getFlux() -= delta;
+    }
+
+    /// @copydoc PhotometryModel::computeResidual
+    double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
+
+    /// @copydoc PhotometryModel::computeRefResidual
+    double computeRefResidual(FittedStar const &fittedStar, RefStar const &refStar) const override {
+        return fittedStar.getFlux() - refStar.getFlux();
+    };
+
+    /// @copydoc PhotometryModel::getRefError
+    double getRefError(RefStar const &refStar) const override { return refStar.getFluxErr(); }
+
+    /// @copydoc PhotometryModel::transform
+    double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
+
+    /// @copydoc PhotometryModel::transformError
+    double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
+
+    /// @copydoc PhotometryModel::toPhotoCalib
+    std::shared_ptr<afw::image::PhotoCalib> toPhotoCalib(CcdImage const &ccdImage) const override;
+
+protected:
+    /// @copydoc ConstrainedPhotometryModel::initialChipCalibration
+    double initialChipCalibration(std::shared_ptr<afw::image::PhotoCalib const> photoCalib) {
+        return photoCalib->getCalibrationMean();
+    }
+};
+
+class ConstrainedMagnitudeModel : public ConstrainedPhotometryModel {
+public:
+    explicit ConstrainedMagnitudeModel(CcdImageList const &ccdImageList,
+                                       afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7)
+            : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder) {
+        initialize<MagnitudeTransfoSpatiallyInvariant, MagnitudeTransfoChebyshev, ChipVisitMagnitudeMapping>(
+                ccdImageList, focalPlaneBBox, visitOrder);
+    }
+
+    /// @copydoc PhotometryModel::offsetFittedStar
+    void offsetFittedStar(FittedStar &fittedStar, double delta) const override {
+        fittedStar.getMag() -= delta;
+    }
+
+    /// @copydoc PhotometryModel::computeResidual
+    double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
+
+    /// @copydoc PhotometryModel::computeRefResidual
+    double computeRefResidual(FittedStar const &fittedStar, RefStar const &refStar) const override {
+        return fittedStar.getMag() - refStar.getMag();
+    };
+
+    /// @copydoc PhotometryModel::getRefError
+    double getRefError(RefStar const &refStar) const override { return refStar.getMagErr(); }
+
+    /// @copydoc PhotometryModel::transform
+    double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
+
+    /// @copydoc PhotometryModel::transformError
+    double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;
+
+    /// @copydoc PhotometryModel::toPhotoCalib
+    std::shared_ptr<afw::image::PhotoCalib> toPhotoCalib(CcdImage const &ccdImage) const override;
+
+protected:
+    /// @copydoc ConstrainedPhotometryModel::initialChipCalibration
+    double initialChipCalibration(std::shared_ptr<afw::image::PhotoCalib const> photoCalib) {
+        return magFromFlux(photoCalib->getCalibrationMean());
+    }
 };
 
 }  // namespace jointcal

include/lsst/jointcal/FitterBase.h

@@ -17,7 +17,8 @@ namespace jointcal {
 enum class MinimizeResult {
     Converged,      // fit has converged - no more outliers
     Chi2Increased,  // still some ouliers but chi2 increases
-    Failed          // factorization failed
+    Failed,         // factorization failed
+    NonFinite       // non-finite chi2 statistic
 };
 
 /**

include/lsst/jointcal/PhotometryMapping.h

@@ -36,21 +36,21 @@ class PhotometryMappingBase {
      * Return the on-sky transformed flux for measuredStar on ccdImage.
      *
      * @param[in]  measuredStar  The measured star position to transform.
-     * @param[in]  instFlux      The instrument flux to transform.
+     * @param[in]  value         The instrument flux or magnitude to transform.
      *
-     * @return     The on-sky flux transformed from instFlux at measuredStar's position.
+     * @return     The on-sky value () transformed from value at measuredStar's position.
      */
-    virtual double transform(MeasuredStar const &measuredStar, double instFlux) const = 0;
+    virtual double transform(MeasuredStar const &measuredStar, double value) const = 0;
 
     /**
      * Return the on-sky transformed flux uncertainty for measuredStar on ccdImage.
      * Matches the underlying PhotometryTransfo's `transformError()` until freezeErrorTransform() is called.
      *
      * @param[in]  measuredStar  The measured star position to transform.
-     * @param[in]  value  The value to transform.
-     * @param[in]  valueErr  The value uncertainty to transform.
+     * @param[in]  value  The flux or magnitude to transform.
+     * @param[in]  valueErr  The flux or magnitude uncertainty to transform.
      *
-     * @return     The on-sky flux transformed from instFlux at measuredStar's position.
+     * @return     The on-sky value transformed from value at measuredStar's position.
      */
     virtual double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const = 0;
 
@@ -67,10 +67,10 @@ class PhotometryMappingBase {
      * Compute the derivatives with respect to the parameters (i.e. the coefficients).
      *
      * @param[in]  measuredStar The measured star position to transform.
-     * @param[in]  instFlux     The instrument flux to compute the derivative at.
+     * @param[in]  value        The instrument flux or magnitude to compute the derivative at.
      * @param[out] derivatives  The computed derivatives, in the same order as the deltas in offsetParams.
      */
-    virtual void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
+    virtual void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
                                              Eigen::Ref<Eigen::VectorXd> derivatives) const = 0;
 
     /// Make this mapping's parameters fixed (i.e. not varied during fitting).
@@ -192,35 +192,27 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     ChipVisitPhotometryMapping(std::shared_ptr<PhotometryMapping> chipMapping,
                                std::shared_ptr<PhotometryMapping> visitMapping)
             : PhotometryMappingBase(),
-              _nParChips(0),
-              _nParVisits(0),
+              _nParChip(0),
+              _nParVisit(0),
               _chipMapping(std::move(chipMapping)),
               _visitMapping(std::move(visitMapping)) {}
 
     /// @copydoc PhotometryMappingBase::getNpar
-    unsigned getNpar() const override { return _nParChips + _nParVisits; }
+    unsigned getNpar() const override { return _nParChip + _nParVisit; }
 
     /// @copydoc PhotometryMappingBase::transform
-    double transform(MeasuredStar const &measuredStar, double instFlux) const override {
-        double tempFlux = _chipMapping->getTransfo()->transform(measuredStar.x, measuredStar.y, instFlux);
+    double transform(MeasuredStar const &measuredStar, double value) const override {
+        double temp = _chipMapping->getTransfo()->transform(measuredStar.x, measuredStar.y, value);
         return _visitMapping->getTransfo()->transform(measuredStar.getXFocal(), measuredStar.getYFocal(),
-                                                      tempFlux);
+                                                      temp);
     }
 
-    /// @copydoc PhotometryMappingBase::transformError
-    double transformError(MeasuredStar const &measuredStar, double instFlux,
-                          double instFluxErr) const override;
-
     /// @copydoc PhotometryMappingBase::freezeErrorTransform
     void freezeErrorTransform() override {
         _chipMapping->freezeErrorTransform();
         _visitMapping->freezeErrorTransform();
     }
 
-    /// @copydoc PhotometryMappingBase::computeParameterDerivatives
-    void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
-                                     Eigen::Ref<Eigen::VectorXd> derivatives) const override;
-
     /// @copydoc PhotometryMappingBase::getParameters
     Eigen::VectorXd getParameters() override {
         Eigen::VectorXd joined(getNpar());
@@ -253,15 +245,51 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
     std::shared_ptr<PhotometryMapping> getChipMapping() const { return _chipMapping; }
     std::shared_ptr<PhotometryMapping> getVisitMapping() const { return _visitMapping; }
 
-private:
+    unsigned getNParChip() const { return _nParChip; }
+    unsigned getNParVisit() const { return _nParVisit; }
+
+protected:
     // These are either transfo.getNpar() or 0, depending on whether we are fitting that component or not.
-    unsigned _nParChips, _nParVisits;
+    unsigned _nParChip, _nParVisit;
 
     // the actual transformation to be fit
     std::shared_ptr<PhotometryMapping> _chipMapping;
     std::shared_ptr<PhotometryMapping> _visitMapping;
 };
 
+class ChipVisitFluxMapping : public ChipVisitPhotometryMapping {
+public:
+    ChipVisitFluxMapping(std::shared_ptr<PhotometryMapping> chipMapping,
+                         std::shared_ptr<PhotometryMapping> visitMapping)
+            : ChipVisitPhotometryMapping(chipMapping, visitMapping) {}
+
+    /// @copydoc PhotometryMappingBase::transformError
+    double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override;
+
+    /// @copydoc PhotometryMappingBase::computeParameterDerivatives
+    void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
+                                     Eigen::Ref<Eigen::VectorXd> derivatives) const override;
+};
+
+class ChipVisitMagnitudeMapping : public ChipVisitPhotometryMapping {
+public:
+    ChipVisitMagnitudeMapping(std::shared_ptr<PhotometryMapping> chipMapping,
+                              std::shared_ptr<PhotometryMapping> visitMapping)
+            : ChipVisitPhotometryMapping(chipMapping, visitMapping) {}
+
+    /**
+     * @copydoc PhotometryMappingBase::transformError
+     *
+     * @note This method takes instFlux and instFluxErr: the error calculation has to
+     * use fluxes to get the math right.
+     */
+    double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override;
+
+    /// @copydoc PhotometryMappingBase::computeParameterDerivatives
+    void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
+                                     Eigen::Ref<Eigen::VectorXd> derivatives) const override;
+};
+
 }  // namespace jointcal
 }  // namespace lsst
 

include/lsst/jointcal/PhotometryModel.h

@@ -34,7 +34,7 @@ class PhotometryModel {
     virtual unsigned assignIndices(std::string const &whatToFit, unsigned firstIndex) = 0;
 
     /**
-     * Offset the parameters by the provided amounts.
+     * Offset the parameters by the provided amounts (by -delta).
      *
      * The shifts are applied according to the indices given in assignIndices.
      *
@@ -43,13 +43,25 @@ class PhotometryModel {
     virtual void offsetParams(Eigen::VectorXd const &delta) = 0;
 
     /**
-     * Offset the appropriate flux or magnitude.
+     * Offset the appropriate flux or magnitude (by -delta).
      *
      * @param fittedStar The star to update.
      * @param delta The amount to update by.
      */
     virtual void offsetFittedStar(FittedStar &fittedStar, double delta) const = 0;
 
+    /**
+     * Compute the residual between the model applied to a star and its associated fittedStar.
+     *
+     * @f[
+     *     residual = Model(measuredStar) - fittedStar
+     * @f]
+     *
+     * @param ccdImage The ccdImage where measuredStar resides.
+     * @param measuredStar The measured star position to compute the residual of.
+     *
+     * @return The residual.
+     */
     virtual double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const = 0;
 
     /**