Merge pull request #91 from lsst/tickets/DM-14510

DM-14510 implement line search in jointcal

include/lsst/jointcal/CcdImage.h

@@ -25,7 +25,13 @@ typedef std::list<std::shared_ptr<CcdImage>> CcdImageList;
 typedef int VisitIdType;
 typedef int CcdIdType;
 /// For hashing a ccdImage: the pair of (visit, ccd) IDs should be unique to each ccdImage.
-typedef std::pair<VisitIdType, CcdIdType> CcdImageKey;
+struct CcdImageKey {
+    VisitIdType visit;
+    CcdIdType ccd;
+    bool operator!=(CcdImageKey const &right) const { return !(*this == right); }
+    bool operator==(CcdImageKey const &right) const { return (visit == right.visit) && (ccd == right.ccd); }
+};
+// typedef std::pair<VisitIdType, CcdIdType> CcdImageKey;
 std::ostream &operator<<(std::ostream &out, CcdImageKey const &key);
 
 /**
@@ -119,7 +125,7 @@ class CcdImage {
 
     std::shared_ptr<afw::cameraGeom::Detector> getDetector() const { return _detector; }
 
-    CcdImageKey getHashKey() const { return CcdImageKey(_visit, _ccdId); }
+    CcdImageKey getHashKey() const { return CcdImageKey{_visit, _ccdId}; }
 
     //!  Airmass
     double getAirMass() const { return _airMass; }
@@ -210,9 +216,8 @@ template <>
  * most-significant-bit, and the visitId in the least for a simple, unique, hash per ccdImage.
  */
 struct hash<lsst::jointcal::CcdImageKey> {
-    size_t operator()(lsst::jointcal::CcdImageKey const &ccdImage) const {
-        return hash<size_t>()(static_cast<size_t>(ccdImage.first) |
-                              (static_cast<size_t>(ccdImage.second) << 32));
+    size_t operator()(lsst::jointcal::CcdImageKey const &key) const {
+        return hash<size_t>()(static_cast<size_t>(key.visit) | (static_cast<size_t>(key.ccd) << 32));
     }
 };
 }  // namespace std

include/lsst/jointcal/FitterBase.h

@@ -41,18 +41,23 @@ class FitterBase {
      * Does a 1 step minimization, assuming a linear model.
      *
      * This is a complete Newton Raphson step. Compute first and second
-     * derivatives, solve for the step and apply it, without a line search.
+     * derivatives, solve for the step and apply it, with an optional line search.
      *
      * It calls assignIndices, leastSquareDerivatives, solves the linear system and calls
      * offsetParams, then removes outliers in a loop if requested.
-     * Relies on sparse linear algebra.
+     * Relies on sparse linear algebra via Eigen's CholmodSupport package.
      *
      * @param[in]  whatToFit  See child method assignIndices for valid string values.
      * @param[in]  nSigmaCut  How many sigma to reject outliers at. Outlier
      *                        rejection ignored for nSigmaCut=0.
      * @param[in]  doRankUpdate  Use CholmodSimplicialLDLT2.update() to do a fast rank update after outlier
      *                           removal; otherwise do a slower full recomputation of the matrix.
      *                           Only matters if nSigmaCut != 0.
+     * @param[in]  doLineSearch  Use boost's brent_find_minima to perform a line search after the gradient
+     *                           solution is found, and apply the scale factor to the computed offsets.
+     *                           The line search is done in the domain [-1, 2], but if the scale factor
+     *                           is far from 1.0, then the problem is likely in a significantly non-linear
+     *                           regime.
      * @param[in] dumpMatrixFile  Write the pre-fit Hessian matrix and gradient to the files with "-mat.txt"
      *                            and "-grad.txt". Be aware, this requires a large increase in memory usage
      *                            to create a dense matrix before writing it; the output file may be large.
@@ -74,7 +79,8 @@ class FitterBase {
      *         the fitted parameters, if the calculations and indices are defined correctly.
      */
     MinimizeResult minimize(std::string const &whatToFit, double const nSigmaCut = 0,
-                            bool const doRankUpdate = true, std::string const &dumpMatrixFile = "");
+                            bool const doRankUpdate = true, bool const doLineSearch = false,
+                            std::string const &dumpMatrixFile = "");
 
     /**
      * Returns the chi2 for the current state.
@@ -189,6 +195,19 @@ class FitterBase {
     /// Compute the derivatives of the reference terms
     virtual void leastSquareDerivativesReference(FittedStarList const &fittedStarList,
                                                  TripletList &tripletList, Eigen::VectorXd &grad) const = 0;
+
+private:
+    /**
+     * Performe a line search along vector delta, returning a scale factor for the minimum.
+     *
+     * Note that this offsets and restores the model during each iteration of the line search,
+     * as part of the minimization schema.
+     *
+     * @param delta The vector of offsets that is expected to reach the minimium value.
+     *
+     * @return The scale factor to apply to delta that gets it to the true minimum.
+     */
+    double _lineSearch(Eigen::VectorXd const &delta);
 };
 }  // namespace jointcal
 }  // namespace lsst

include/lsst/jointcal/SimpleAstrometryModel.h

@@ -71,7 +71,7 @@ class SimpleAstrometryModel : public AstrometryModel {
     ~SimpleAstrometryModel(){};
 
 private:
-    std::map<CcdImageKey, std::unique_ptr<SimpleGtransfoMapping>> _myMap;
+    std::unordered_map<CcdImageKey, std::unique_ptr<SimpleGtransfoMapping>> _myMap;
     const std::shared_ptr<ProjectionHandler const> _sky2TP;
 
     /// @copydoc AstrometryModel::findMapping

python/lsst/jointcal/fitter.cc

@@ -41,7 +41,7 @@ void declareFitterBase(py::module &mod) {
     py::class_<FitterBase, std::shared_ptr<FitterBase>> cls(mod, "FitterBase");
 
     cls.def("minimize", &FitterBase::minimize, "whatToFit"_a, "nSigRejCut"_a = 0, "doRankUpdate"_a = true,
-            "dumpMatrixFile"_a = "");
+            "doLineSearch"_a = false, "dumpMatrixFile"_a = "");
     cls.def("computeChi2", &FitterBase::computeChi2);
     cls.def("saveChi2Contributions", &FitterBase::saveChi2Contributions);
 }

python/lsst/jointcal/jointcal.py

@@ -152,6 +152,12 @@ class JointcalConfig(pexConfig.Config):
         dtype=int,
         default=30,
     )
+    allowLineSearch = pexConfig.Field(
+        doc="Allow a line search during minimization, if it is reasonable for the model"
+        " (models with a significant non-linear component, e.g. constrainedPhotometry).",
+        dtype=bool,
+        default=False
+    )
     astrometrySimpleOrder = pexConfig.Field(
         doc="Polynomial order for fitting the simple astrometry model.",
         dtype=int,
@@ -594,8 +600,11 @@ def _fit_photometry(self, associations, dataName=None):
             model = lsst.jointcal.ConstrainedPhotometryModel(associations.getCcdImageList(),
                                                              self.focalPlaneBBox,
                                                              visitOrder=self.config.photometryVisitOrder)
+            # potentially nonlinear problem, so we may need a line search to converge.
+            doLineSearch = self.config.allowLineSearch
         elif self.config.photometryModel == "simple":
             model = lsst.jointcal.SimplePhotometryModel(associations.getCcdImageList())
+            doLineSearch = False  # purely linear problem, so no line search needed
 
         fit = lsst.jointcal.PhotometryFit(associations, model)
         chi2 = fit.computeChi2()
@@ -613,17 +622,19 @@ def _fit_photometry(self, associations, dataName=None):
         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)
             chi2 = fit.computeChi2()
             self.log.info(str(chi2))
             dumpMatrixFile = ""  # so we don't redo the output on the next step
-        fit.minimize("Model", dumpMatrixFile=dumpMatrixFile)
+        fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
         chi2 = fit.computeChi2()
         self.log.info(str(chi2))
-        fit.minimize("Fluxes")
+        fit.minimize("Fluxes")  # no line search: always purely linear.
         chi2 = fit.computeChi2()
         self.log.info(str(chi2))
-        fit.minimize("Model Fluxes")
+        fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
         chi2 = fit.computeChi2()
         if not np.isfinite(chi2.chi2):
             raise FloatingPointError('Pre-iteration chi2 is invalid: %s'%chi2)
@@ -638,7 +649,8 @@ def _fit_photometry(self, associations, dataName=None):
                                  self.config.maxPhotometrySteps,
                                  "photometry",
                                  "Model Fluxes",
-                                 doRankUpdate=self.config.photometryDoRankUpdate)
+                                 doRankUpdate=self.config.photometryDoRankUpdate,
+                                 doLineSearch=doLineSearch)
 
         add_measurement(self.job, 'jointcal.photometry_final_chi2', chi2.chi2)
         add_measurement(self.job, 'jointcal.photometry_final_ndof', chi2.ndof)
@@ -746,7 +758,8 @@ def _check_stars(self, associations):
                 self.log.warn("ccdImage %s has only %s RefStars (desired %s)",
                               ccdImage.getName(), nRefStars, self.config.minRefStarsPerCcd)
 
-    def _iterate_fit(self, associations, fit, model, max_steps, name, whatToFit, doRankUpdate=True):
+    def _iterate_fit(self, associations, fit, model, max_steps, name, whatToFit, doRankUpdate=True,
+                     doLineSearch=False):
         """Run fit.minimize up to max_steps times, returning the final chi2."""
 
         dumpMatrixFile = "%s_postinit" % name if self.config.writeInitMatrix else ""
@@ -755,6 +768,7 @@ def _iterate_fit(self, associations, fit, model, max_steps, name, whatToFit, doR
             r = fit.minimize(whatToFit,
                              self.config.outlierRejectSigma,
                              doRankUpdate=doRankUpdate,
+                             doLineSearch=doLineSearch,
                              dumpMatrixFile=dumpMatrixFile)
             dumpMatrixFile = ""  # clear it so we don't write the matrix again.
             chi2 = fit.computeChi2()

src/CcdImage.cc

@@ -26,7 +26,7 @@ namespace lsst {
 namespace jointcal {
 
 std::ostream &operator<<(std::ostream &out, CcdImageKey const &key) {
-    out << "(visit: " << key.first << ", ccd: " << key.second << ")";
+    out << "(visit: " << key.visit << ", ccd: " << key.ccd << ")";
     return out;
 }
 

src/FitterBase.cc

@@ -1,6 +1,8 @@
 #include <vector>
 #include "Eigen/Core"
 
+#include <boost/math/tools/minima.hpp>
+
 #include "lsst/log/Log.h"
 
 #include "lsst/jointcal/Chi2.h"
@@ -141,7 +143,7 @@ void dumpMatrixAndGradient(SparseMatrixD const &matrix, Eigen::VectorXd const &g
 }  // namespace
 
 MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaCut, bool doRankUpdate,
-                                    std::string const &dumpMatrixFile) {
+                                    bool const doLineSearch, std::string const &dumpMatrixFile) {
     assignIndices(whatToFit);
 
     MinimizeResult returnCode = MinimizeResult::Converged;
@@ -151,6 +153,7 @@ MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaC
     TripletList tripletList(nTrip);
     Eigen::VectorXd grad(_nParTot);
     grad.setZero();
+    double scale = 1.0;
 
     // Fill the triplets
     leastSquareDerivatives(tripletList, grad);
@@ -186,7 +189,10 @@ MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaC
 
     while (true) {
         Eigen::VectorXd delta = chol.solve(grad);
-        offsetParams(delta);
+        if (doLineSearch) {
+            scale = _lineSearch(delta);
+        }
+        offsetParams(scale * delta);
         Chi2Statistic currentChi2(computeChi2());
         LOGLS_DEBUG(_log, currentChi2);
         if (currentChi2.chi2 > oldChi2 && totalMeasOutliers + totalRefOutliers != 0) {
@@ -299,5 +305,21 @@ void FitterBase::saveChi2Contributions(std::string const &baseName) const {
     saveChi2RefContributions(refTuple);
 }
 
+double FitterBase::_lineSearch(Eigen::VectorXd const &delta) {
+    auto func = [this, &delta](double scale) {
+        auto offset = scale * delta;
+        offsetParams(offset);
+        auto chi2 = computeChi2();
+        // reset the system to where it was before offsetting.
+        offsetParams(-offset);
+        return chi2.chi2;
+    };
+    // The maximum theoretical precision is half the number of bits in the mantissa (see boost docs).
+    auto bits = std::numeric_limits<double>::digits / 2;
+    auto result = boost::math::tools::brent_find_minima(func, -1.0, 2.0, bits);
+    LOGLS_DEBUG(_log, "Line search scale factor: " << result.first);
+    return result.first;
+}
+
 }  // namespace jointcal
 }  // namespace lsst

tests/test_jointcal_cfht.py

@@ -172,6 +172,16 @@ def test_jointcalTask_2_visits_constrainedPhotometry_no_rank_update(self):
 
         self._testJointcalTask(2, None, None, pa1, metrics=metrics)
 
+    def test_jointcalTask_2_visits_constrainedPhotometry_lineSearch(self):
+        """Activating the line search should only slightly change the chi2 in this case."""
+        pa1, metrics = self.setup_jointcalTask_2_visits_constrainedPhotometry()
+        self.config.allowLineSearch = True
+
+        # Only this value should differ from the metrics defined in setup above.
+        metrics['photometry_final_chi2'] = 2642.54
+
+        self._testJointcalTask(2, None, None, pa1, metrics=metrics)
+
     def test_jointcalTask_2_visits_constrainedPhotometry_flagged(self):
         """Test the use of the FlaggedSourceSelector."""
         pa1, metrics = self.setup_jointcalTask_2_visits_constrainedPhotometry()