Merge pull request #182 from lsst/tickets/DM-30254

DM-30254: Fix jointcal crash when doing outlier rejection on only the model

include/lsst/jointcal/Associations.h

@@ -68,7 +68,8 @@ class Associations {
      */
     Associations()
             : _commonTangentPoint(Point(std::numeric_limits<double>::quiet_NaN(),
-                                        std::numeric_limits<double>::quiet_NaN())) {}
+                                        std::numeric_limits<double>::quiet_NaN())),
+              _maxMeasuredStars(0) {}
 
     /**
      * Create an Associations object from a pre-built list of ccdImages.
@@ -81,7 +82,8 @@ class Associations {
     Associations(CcdImageList const &imageList)
             : ccdImageList(imageList),
               _commonTangentPoint(Point(std::numeric_limits<double>::quiet_NaN(),
-                                        std::numeric_limits<double>::quiet_NaN())) {}
+                                        std::numeric_limits<double>::quiet_NaN())),
+              _maxMeasuredStars(0) {}
 
     /// No moves or copies: jointcal only ever needs one Associations object.
     Associations(Associations const &) = delete;
@@ -104,6 +106,8 @@ class Associations {
     //! can be used to project sidereal coordinates related to the image set on a plane.
     Point getCommonTangentPoint() const { return _commonTangentPoint; }
 
+    size_t getMaxMeasuredStars() const { return _maxMeasuredStars; }
+
     /**
      * @brief      Create a ccdImage from an exposure catalog and metadata, and add it to the list.
      *
@@ -119,9 +123,8 @@ class Associations {
      * @param[in]  control    The JointcalControl object
      */
     void createCcdImage(afw::table::SourceCatalog &catalog, std::shared_ptr<lsst::afw::geom::SkyWcs> wcs,
-                        std::shared_ptr<lsst::afw::image::VisitInfo> visitInfo,
-                        lsst::geom::Box2I const &bbox, std::string const &filter,
-                        std::shared_ptr<afw::image::PhotoCalib> photoCalib,
+                        std::shared_ptr<lsst::afw::image::VisitInfo> visitInfo, lsst::geom::Box2I const &bbox,
+                        std::string const &filter, std::shared_ptr<afw::image::PhotoCalib> photoCalib,
                         std::shared_ptr<afw::cameraGeom::Detector> detector, int visit, int ccd,
                         lsst::jointcal::JointcalControl const &control);
 
@@ -168,6 +171,16 @@ class Associations {
      */
     void prepareFittedStars(int minMeasurements);
 
+    /**
+     * Remove FittedStars that have no measured stars; this can happen after outlier rejection.
+     *
+     * Use this to perform e.g. position minimization with outlier rejection after model minimization has
+     * removed measuredStar outliers, to prevent the matrix from becoming non-positive definite.
+     * Once this has been called, prepareFittedStars() has to be called again if the full
+     * measuredStar->FittedStar relationship needs to be reconstructed.
+     */
+    void cleanFittedStars();
+
     CcdImageList const &getCcdImageList() const { return ccdImageList; }
 
     //! Number of different bands in the input image list. Not implemented so far
@@ -209,9 +222,14 @@ class Associations {
      * Only call after selectFittedStars() has been called: it assumes that each measuredStar points to a
      * fittedStar, and that the measurementCount for each fittedStar is correct.
      */
-    void normalizeFittedStars() const;
+    void normalizeFittedStars();
 
     Point _commonTangentPoint;
+
+    // The number of MeasuredStars at the start of fitting, before any outliers are removed.
+    // This is used to reserve space in vectors for e.g. outlier removal, but is not updated during outlier
+    // removal or cleanup, so should only be used as an upper bound on the number of MeasuredStars.
+    size_t _maxMeasuredStars;
 };
 
 }  // namespace jointcal

include/lsst/jointcal/AstrometryFit.h

@@ -142,8 +142,6 @@ class AstrometryFit : public FitterBase {
     double _JDRef;                    // average Julian date
 
     // counts in parameter subsets.
-    std::size_t _nParDistortions;
-    std::size_t _nParPositions;
     std::size_t _nParRefrac;
 
     double _posError;  // constant term on error on position (in pixel unit)
@@ -159,8 +157,7 @@ class AstrometryFit : public FitterBase {
 
     void accumulateStatRefStars(Chi2Accumulator &accum) const override;
 
-    void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar,
-                                  IndexVector &indices) const override;
+    void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const override;
 
     Point transformFittedStar(FittedStar const &fittedStar, AstrometryTransform const &sky2TP,
                               Point const &refractionVector, double refractionCoeff, double mjd) const;

include/lsst/jointcal/BaseStar.h

@@ -80,7 +80,8 @@ class BaseStar : public FatPoint {
     }
 
     virtual void print(std::ostream &out) const {
-        out << "x: " << x << " y: " << y << " flux: " << _flux << " fluxErr: " << _fluxErr;
+        FatPoint::print(out);
+        out << " flux: " << std::setprecision(9) << _flux << " fluxErr: " << _fluxErr;
     }
 
     BaseStar &operator=(Point const &point) {

include/lsst/jointcal/FatPoint.h

@@ -48,7 +48,7 @@ class FatPoint : public Point {
 
     void print(std::ostream& s = std::cout) const {
         Point::print(s);
-        s << " vxx,vyy,vxy " << vx << ' ' << vy << ' ' << vxy;
+        s << std::setprecision(8) << " vxx,vyy,vxy " << vx << ',' << vy << ',' << vxy;
     }
 };
 }  // namespace jointcal

include/lsst/jointcal/FitterBase.h

@@ -53,7 +53,12 @@ enum class MinimizeResult {
 class FitterBase {
 public:
     explicit FitterBase(std::shared_ptr<Associations> associations)
-            : _associations(associations), _whatToFit(""), _lastNTrip(0), _nParTot(0), _nMeasuredStars(0) {}
+            : _associations(associations),
+              _whatToFit(""),
+              _lastNTrip(0),
+              _nTotal(0),
+              _nModelParams(0),
+              _nStarParams(0) {}
 
     /// No copy or move: there is only ever one fitter of a given type.
     FitterBase(FitterBase const &) = delete;
@@ -156,9 +161,10 @@ class FitterBase {
     std::shared_ptr<Associations> _associations;
     std::string _whatToFit;
 
-    Eigen::Index _lastNTrip;  // last triplet count, used to speed up allocation
-    Eigen::Index _nParTot;
-    Eigen::Index _nMeasuredStars;
+    Eigen::Index _lastNTrip;     // last triplet count, used to speed up allocation
+    Eigen::Index _nTotal;        // Total number of parameters being fit.
+    Eigen::Index _nModelParams;  // Number of model parameters that are being fit.
+    Eigen::Index _nStarParams;   // Number of star positions/fluxes that are being fit.
 
     // lsst.logging instance, to be created by subclass so that messages have consistent name while fitting.
     LOG_LOGGER _log;
@@ -203,8 +209,7 @@ class FitterBase {
     void removeRefOutliers(FittedStarList &outliers);
 
     /// Set the indices of a measured star from the full matrix, for outlier removal.
-    virtual void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar,
-                                          IndexVector &indices) const = 0;
+    virtual void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const = 0;
 
     /// Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.
     virtual void accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const = 0;

include/lsst/jointcal/PhotometryFit.h

@@ -55,9 +55,7 @@ class PhotometryFit : public FitterBase {
             : FitterBase(associations),
               _fittingModel(false),
               _fittingFluxes(false),
-              _photometryModel(photometryModel),
-              _nParModel(0),
-              _nParFluxes(0) {
+              _photometryModel(photometryModel) {
         _log = LOG_GET("jointcal.PhotometryFit");
     }
 
@@ -98,16 +96,11 @@ class PhotometryFit : public FitterBase {
     bool _fittingModel, _fittingFluxes;
     std::shared_ptr<PhotometryModel> _photometryModel;
 
-    // counts in parameter subsets.
-    std::size_t _nParModel;
-    std::size_t _nParFluxes;
-
     void accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const override;
 
     void accumulateStatRefStars(Chi2Accumulator &accum) const override;
 
-    void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar,
-                                  IndexVector &indices) const override;
+    void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const override;
 
     void leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList,
                                            Eigen::VectorXd &grad,

include/lsst/jointcal/Point.h

@@ -25,6 +25,7 @@
 #ifndef LSST_JOINTCAL_POINT_H
 #define LSST_JOINTCAL_POINT_H
 #include <iostream>
+#include <iomanip>
 #include <cmath>
 
 namespace lsst {
@@ -62,10 +63,10 @@ class Point {
     //! Difference
     Point operator-(const Point& Right) const { return Point(x - Right.x, y - Right.y); }
 
-    //! utility
-    virtual void print(std::ostream& s = std::cout) const { s << "x: " << x << " y: " << y; }
+    virtual void print(std::ostream& s = std::cout) const {
+        s << std::setprecision(15) << "x: " << x << " y: " << y;
+    }
 
-    //! -
     friend std::ostream& operator<<(std::ostream& stream, const Point& point) {
         point.print(stream);
         return stream;

python/lsst/jointcal/associations.cc

@@ -55,6 +55,7 @@ void declareAssociations(py::module &mod) {
     cls.def("createCcdImage", &Associations::createCcdImage);
     cls.def("addCcdImage", &Associations::addCcdImage);
     cls.def("prepareFittedStars", &Associations::prepareFittedStars);
+    cls.def("cleanFittedStars", &Associations::cleanFittedStars);
 
     cls.def("getCcdImageList", &Associations::getCcdImageList, py::return_value_policy::reference_internal);
     cls.def_property_readonly("ccdImageList", &Associations::getCcdImageList,

python/lsst/jointcal/jointcal.py

@@ -451,7 +451,7 @@ class JointcalConfig(pipeBase.PipelineTaskConfig,
         default=20,
     )
     maxAstrometrySteps = pexConfig.Field(
-        doc="Maximum number of minimize iterations to take when fitting photometry.",
+        doc="Maximum number of minimize iterations to take when fitting astrometry.",
         dtype=int,
         default=20,
     )

src/Associations.cc

@@ -344,7 +344,7 @@ void Associations::selectFittedStars(int minMeasurements) {
     LOGLS_INFO(_log, "Total, valid number of Measured stars: " << totalMeasured << ", " << validMeasured);
 }
 
-void Associations::normalizeFittedStars() const {
+void Associations::normalizeFittedStars() {
     // Clear positions in order to take the average of the measuredStars.
     for (auto &fittedStar : fittedStarList) {
         fittedStar->x = 0.0;
@@ -370,15 +370,36 @@ void Associations::normalizeFittedStars() const {
         }
     }
 
+    _maxMeasuredStars = 0;
     for (auto &fi : fittedStarList) {
         auto measurementCount = fi->getMeasurementCount();
+        _maxMeasuredStars += measurementCount;
         fi->x /= measurementCount;
         fi->y /= measurementCount;
         fi->getFlux() /= measurementCount;
         fi->getMag() = utils::nanojanskyToABMagnitude(fi->getFlux());
     }
 }
 
+void Associations::cleanFittedStars() {
+    auto iter = fittedStarList.begin();
+    auto end = fittedStarList.end();
+    size_t count = 0;
+    while (iter != end) {
+        auto fittedStar = *iter;
+        if (fittedStar->getMeasurementCount() == 0) {
+            LOGLS_TRACE(_log, "Deleting FittedStar (has no measuredStars): " << *fittedStar);
+            iter = fittedStarList.erase(iter);
+            count++;
+        } else {
+            ++iter;
+        }
+    }
+    if (count > 0) {
+        LOGLS_INFO(_log, "Removed " << count << " fittedStars that had no associated measuredStar.");
+    }
+}
+
 void Associations::assignMags() {
     for (auto const &ccdImage : ccdImageList) {
         MeasuredStarList &catalog = ccdImage->getCatalogForFit();

src/AstrometryFit.cc

@@ -48,8 +48,6 @@ AstrometryFit::AstrometryFit(std::shared_ptr<Associations> associations,
         : FitterBase(associations),
           _astrometryModel(astrometryModel),
           _refractionCoefficient(0),
-          _nParDistortions(0),
-          _nParPositions(0),
           _nParRefrac(_associations->getNFilters()),
           _posError(posError) {
     _log = LOG_GET("jointcal.AstrometryFit");
@@ -428,8 +426,7 @@ void AstrometryFit::accumulateStatRefStars(Chi2Accumulator &accum) const {
 //! this routine is to be used only in the framework of outlier removal
 /*! it fills the array of indices of parameters that a Measured star
     constrains. Not really all of them if you check. */
-void AstrometryFit::getIndicesOfMeasuredStar(MeasuredStar const &measuredStar,
-                                             IndexVector &indices) const {
+void AstrometryFit::getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const {
     if (_fittingDistortions) {
         const AstrometryMapping *mapping = _astrometryModel->getMapping(measuredStar.getCcdImage());
         mapping->getMappingIndices(indices);
@@ -462,9 +459,9 @@ void AstrometryFit::assignIndices(std::string const &whatToFit) {
     _fittingPM = (_whatToFit.find("PM") != std::string::npos);
     // When entering here, we assume that whatToFit has already been interpreted.
 
-    _nParDistortions = 0;
-    if (_fittingDistortions) _nParDistortions = _astrometryModel->assignIndices(_whatToFit, 0);
-    std::size_t ipar = _nParDistortions;
+    _nModelParams = 0;
+    if (_fittingDistortions) _nModelParams = _astrometryModel->assignIndices(_whatToFit, 0);
+    std::size_t ipar = _nModelParams;
 
     if (_fittingPos) {
         FittedStarList &fittedStarList = _associations->fittedStarList;
@@ -478,16 +475,18 @@ void AstrometryFit::assignIndices(std::string const &whatToFit) {
             if ((_fittingPM)&fittedStar->mightMove) ipar += NPAR_PM;
         }
     }
-    _nParPositions = ipar - _nParDistortions;
+    _nStarParams = ipar - _nModelParams;
     if (_fittingRefrac) {
         _refracPosInMatrix = ipar;
         ipar += _nParRefrac;
     }
-    _nParTot = ipar;
+    _nTotal = ipar;
+    LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
+                                            << " values: " << _nStarParams);
 }
 
 void AstrometryFit::offsetParams(Eigen::VectorXd const &delta) {
-    if (delta.size() != _nParTot)
+    if (delta.size() != _nTotal)
         throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                           "AstrometryFit::offsetParams : the provided vector length is not compatible with "
                           "the current whatToFit setting");
@@ -529,13 +528,13 @@ void AstrometryFit::checkStuff() {
     for (unsigned k = 0; k < sizeof(what2fit) / sizeof(what2fit[0]); ++k) {
         assignIndices(what2fit[k]);
         TripletList tripletList(10000);
-        Eigen::VectorXd grad(_nParTot);
+        Eigen::VectorXd grad(_nTotal);
         grad.setZero();
         leastSquareDerivatives(tripletList, grad);
-        SparseMatrixD jacobian(_nParTot, tripletList.getNextFreeIndex());
+        SparseMatrixD jacobian(_nTotal, tripletList.getNextFreeIndex());
         jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
         SparseMatrixD hessian = jacobian * jacobian.transpose();
-        LOGLS_DEBUG(_log, "npar : " << _nParTot << ' ' << _nParDistortions);
+        LOGLS_DEBUG(_log, "npar : " << _nTotal << ' ' << _nModelParams);
     }
 }