DM-14786 implement ConstrainedMagnitudeModel #102

parejkoj · 2018-07-30T21:59:11Z

No description provided.

taranu

My main comment is that there is a lot of copy+pasted code between the FluxModel and MagnitudeModel versions of classes that appears to be easily avoidable in Python and somewhat less so in C++. I don't see a compelling reason for the repeated code as it is. Everything else is either minor or a request for clarification of something that I've likely misunderstood.

taranu · 2018-08-06T15:18:15Z

include/lsst/jointcal/ConstrainedPhotometryModel.h

+
+    /// @copydoc PhotometryModel::offsetFittedStar
+    void offsetFittedStar(FittedStar &fittedStar, double delta) const override {
+        fittedStar.getFlux() -= delta;


I realize that this is the same behaviour as in the non-constrained version but subtracting delta instead of adding it is counterintuitive to me. I would expect a method offsetting x by delta to add delta and nothing in the copied docs suggests that the offset is negative, so perhaps the docs should be updated.

Good catch. I have some notes about it in tests, but it doesn't say it in any of the model docs. I've updated the docs.

taranu · 2018-08-06T15:29:41Z

include/lsst/jointcal/ConstrainedPhotometryModel.h

+
+    /// @copydoc PhotometryModel::computeResidual
+    double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override {
+        return transform(ccdImage, measuredStar) - measuredStar.getFittedStar()->getFlux();


Although it's a one-liner, this seems like more than a simple accessor and shouldn't be inlined.

taranu · 2018-08-06T16:17:46Z

include/lsst/jointcal/PhotometryMapping.h

+     * @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;


You could rename the value and valueErr arguments to instFlux and instFluxErr if it only works properly that way. I don't think the style guide forbids this, although I guess you wouldn't be able to copy the docs then.

I think I'd rather avoid duplicating docs when possible and clarifying as necessary (as above). The SimpleMagnitudeModel does take mag and magErr, whereas the two FluxModels take flux and fluxErr, hence my choice of value and valueErr.

taranu · 2018-08-06T16:46:28Z

python/lsst/jointcal/jointcal.py

@@ -624,7 +635,7 @@ def _fit_photometry(self, associations, dataName=None):
        # The constrained model needs the visit transfo fit first; the chip
        # transfo is initialized from the singleFrame PhotoCalib, so it's close.
        dumpMatrixFile = "photometry_preinit" if self.config.writeInitMatrix else ""
-        if self.config.photometryModel == "constrained":
+        if "constrained" in self.config.photometryModel:


Consider self.config.photometryModel.startswith("constrained")

taranu · 2018-08-07T15:29:03Z

src/ConstrainedPhotometryModel.cc

+                                                     afw::geom::Box2D const &focalPlaneBBox, int visitOrder)
+        : ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder) {
+    // keep track of which chip we want to constrain (the one closest to the middle of the focal plane)
+    double minRadius2 = std::numeric_limits<double>::infinity();


This code is almost identical to the ConstrainedFluxModel constructor; the only differences are the type of chipTransfo and a few shared vs unique pointers. Can you not find a better way to share the code? For example:

template the entire class, e.g. with a bool isFluxModel or similar (probably too much work),

split it into a function that takes a visitMap, chipMap, chipVisitMap and is templated by the transform type

I've figured out a templated version that seems to do the trick.

taranu · 2018-08-07T18:34:42Z

tests/test_photometryModel.py

        self.assertEqual(index, expect)


+class ConstrainedFluxModelTestCase(ConstrainedPhotometryModelTestCase, lsst.utils.tests.TestCase):
+    def setUp(self):


Again, this should definitely be shared between Flux and Magnitude, and it seems very easy to do since there's already a useMagnitude property in the MagnitudeModelTestCase.

taranu · 2018-08-07T18:53:21Z

tests/test_photometryTransfo.py

+        expect1 = np.zeros((self.order1+1, self.order1+1), dtype=float)
+
+        expect2 = np.zeros((self.order1+1, self.order1+1), dtype=float)
+        expect2[0, 0] = -1


Besides the usual suggestion to reuse code as above, it would be helpful to add a comment explaining why expect1[0,0] and expect2[0,0] are 0, -1 in the Magnitude case.

I've restructured it to make it obvious.

taranu · 2018-08-07T18:56:06Z

python/lsst/jointcal/jointcal.py

-                    r = fit.minimize(whatToFit, 5)  # outliers removal at 5 sigma.
-                chi2 = fit.computeChi2()
+                    r = fitter.minimize(whatToFit, 5)  # outliers removal at 5 sigma.
+                chi2 = fitter.computeChi2()
                self.log.info("Fit completed with: %s", str(chi2))
                break
            elif r == MinimizeResult.Chi2Increased:
                self.log.warn("still some ouliers but chi2 increases - retry")


ouliers -> outliers

taranu · 2018-08-07T18:58:12Z

python/lsst/jointcal/jointcal.py

+        RuntimeError
+            Raised if the fitter fails for some other reason;
+            log messages will provide further details.
+        """

        dumpMatrixFile = "%s_postinit" % name if self.config.writeInitMatrix else ""
        for i in range(max_steps):
            # outlier removal at 5 sigma.


Isn't it outlier rejection at self.config.outlierRejectSigma sigma, which may not be 5? I see that it's set to 4 for CFHT and DECam.

I've removed the comment: it's unnecessary.

taranu · 2018-08-07T19:06:56Z

python/lsst/jointcal/jointcal.py

            self.log.info(str(chi2))
            if r == MinimizeResult.Converged:
                if doRankUpdate:
                    self.log.debug("fit has converged - no more outliers - redo minimization "
                                   "one more time in case we have lost accuracy in rank update.")
                    # Redo minimization one more time in case we have lost accuracy in rank update
-                    r = fit.minimize(whatToFit, 5)  # outliers removal at 5 sigma.
-                chi2 = fit.computeChi2()
+                    r = fitter.minimize(whatToFit, 5)  # outliers removal at 5 sigma.


Is this really meant to be hardcoded to 5? I would have expected it to use self.config.outlierRejectSigma, otherwise that should be named something like outlierRejectSigmaInitial.

Several classes were split up to make this possible: * ConstrainedPhotometryModel * ChipVisitPhotometryMapping * PhotometryTransfoChebyshev * tests were refactored to test the new classes. * minimized duplicate code between the flux/magnitude versions of the above. Fix constrained model initialization step to apply to both flux and magnitude. Revert chi2 value in cfht lineSearch test that came from me not fixing this in DM-14574 when I changed "constrained"->"constrainedFlux". Add missing overrides to SimplePhotometryModel. Cleanup some docstrings and variable names for consistence between flux and magnitude-based models.

Nonfinite chi2 is now discovered in minimize() and we abort outlier rejection and dump the chi2 meas/ref files to aid in debugging. Add docstring to jointcal._iterate_fit()

This tests for an error in construction that I introduced while cleaning up duplicate code: the magnitude and flux versions are almost identical, except for values of the chip transfos: one is just the photoCalib mean, the other is the log() of that.

taranu reviewed Aug 7, 2018

View reviewed changes

parejkoj force-pushed the tickets/DM-14786 branch 4 times, most recently from 3893824 to 7353d3d Compare August 22, 2018 20:53

parejkoj added 5 commits August 22, 2018 13:53

Remove unneeded lines from testUtils and test_photometryModel

4554dca

Improve handling of nonfinite chi2 during minimization

da9704e

Nonfinite chi2 is now discovered in minimize() and we abort outlier rejection and dump the chi2 meas/ref files to aid in debugging. Add docstring to jointcal._iterate_fit()

Bugfix in plot_photoCalib (hardcoded tract id)

bd290b2

parejkoj force-pushed the tickets/DM-14786 branch from 7353d3d to b240722 Compare August 22, 2018 20:54

parejkoj merged commit 46b432b into master Aug 23, 2018

ktlim deleted the tickets/DM-14786 branch August 25, 2018 06:44

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

DM-14786 implement ConstrainedMagnitudeModel #102

DM-14786 implement ConstrainedMagnitudeModel #102

parejkoj commented Jul 30, 2018

taranu left a comment

taranu Aug 6, 2018

parejkoj Aug 8, 2018

taranu Aug 6, 2018

taranu Aug 6, 2018

parejkoj Aug 8, 2018

taranu Aug 6, 2018

taranu Aug 7, 2018

parejkoj Aug 16, 2018

taranu Aug 7, 2018

taranu Aug 7, 2018

parejkoj Aug 17, 2018

taranu Aug 7, 2018

taranu Aug 7, 2018

parejkoj Aug 17, 2018

taranu Aug 7, 2018

DM-14786 implement ConstrainedMagnitudeModel #102

DM-14786 implement ConstrainedMagnitudeModel #102

Conversation

parejkoj commented Jul 30, 2018

taranu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment