Skip to content

Commit

Permalink
Add fit initialization steps chi2 file output
Browse files Browse the repository at this point in the history
No outliers are rejected in these steps, but they move very far in chi2-space
because they are purely linear (or very nearly so) in the fitted parameters.

Added some tests (using mocks) to check that the filename prefixes are
reasonable.
  • Loading branch information
parejkoj committed Apr 2, 2019
1 parent 4972b01 commit 29627c8
Show file tree
Hide file tree
Showing 2 changed files with 81 additions and 21 deletions.
30 changes: 22 additions & 8 deletions python/lsst/jointcal/jointcal.py
Original file line number Diff line number Diff line change
Expand Up @@ -810,24 +810,31 @@ def _fit_photometry(self, associations, dataName=None):
baseName = None
self._logChi2AndValidate(associations, fit, model, "Initialized", writeChi2Name=baseName)

def getChi2Name(whatToFit):
if self.config.writeChi2FilesOuterLoop:
return f"photometry_init-%s_chi2-{dataName}" % whatToFit
else:
return None

# The constrained model needs the visit transform fit first; the chip
# transform is initialized from the singleFrame PhotoCalib, so it's close.
dumpMatrixFile = "photometry_preinit" if self.config.writeInitMatrix else ""
if self.config.photometryModel.startswith("constrained"):
# no line search: should be purely (or nearly) linear,
# and we want a large step size to initialize with.
fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)
self._logChi2AndValidate(associations, fit, model)
self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("ModelVisit"))
dumpMatrixFile = "" # so we don't redo the output on the next step

fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
self._logChi2AndValidate(associations, fit, model)
self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Model"))

fit.minimize("Fluxes") # no line search: always purely linear.
self._logChi2AndValidate(associations, fit, model)
self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Fluxes"))

fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
self._logChi2AndValidate(associations, fit, model, "Fit prepared")
self._logChi2AndValidate(associations, fit, model, "Fit prepared",
writeChi2Name=getChi2Name("ModelFluxes"))

model.freezeErrorTransform()
self.log.debug("Photometry error scales are frozen.")
Expand Down Expand Up @@ -898,22 +905,29 @@ def _fit_astrometry(self, associations, dataName=None):
baseName = None
self._logChi2AndValidate(associations, fit, model, "Initial", writeChi2Name=baseName)

def getChi2Name(whatToFit):
if self.config.writeChi2FilesOuterLoop:
return f"astrometry_init-%s_chi2-{dataName}" % whatToFit
else:
return None

dumpMatrixFile = "astrometry_preinit" if self.config.writeInitMatrix else ""
# The constrained model needs the visit transform fit first; the chip
# transform is initialized from the detector's cameraGeom, so it's close.
if self.config.astrometryModel == "constrained":
fit.minimize("DistortionsVisit", dumpMatrixFile=dumpMatrixFile)
self._logChi2AndValidate(associations, fit, model)
self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("DistortionsVisit"))
dumpMatrixFile = "" # so we don't redo the output on the next step

fit.minimize("Distortions", dumpMatrixFile=dumpMatrixFile)
self._logChi2AndValidate(associations, fit, model)
self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Distortions"))

fit.minimize("Positions")
self._logChi2AndValidate(associations, fit, model)
self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Positions"))

fit.minimize("Distortions Positions")
self._logChi2AndValidate(associations, fit, model, "Fit prepared")
self._logChi2AndValidate(associations, fit, model, "Fit prepared",
writeChi2Name=getChi2Name("DistortionsPositions"))

chi2 = self._iterate_fit(associations,
fit,
Expand Down
72 changes: 59 additions & 13 deletions tests/test_jointcal.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>.

import unittest
import unittest.mock
from unittest import mock

import numpy as np

Expand Down Expand Up @@ -77,28 +77,31 @@ def setUp(self):

self.maxSteps = 20
self.name = "testing"
self.dataName = "fake-fake"
self.dataName = "fake"
self.whatToFit = "" # unneeded, since we're mocking the fitter

# so the refObjLoaders have something to call `get()` on
# Mock a Butler so the refObjLoaders have something to call `get()` on.
self.butler = unittest.mock.Mock(spec=lsst.daf.persistence.Butler)
self.butler.get.return_value.indexer = DatasetConfig().indexer

# Mock the association manager and give it access to the ccd list above.
self.associations = mock.Mock(spec=lsst.jointcal.Associations)
self.associations.getCcdImageList.return_value = self.ccdImageList

# a default config to be modified by individual tests
self.config = lsst.jointcal.jointcal.JointcalConfig()


class TestJointcalIterateFit(JointcalTestBase, lsst.utils.tests.TestCase):
def setUp(self):
super().setUp()

# Mock the fitter, association manager, and model, so we can force particular
# Mock the fitter and model, so we can force particular
# return values/exceptions. Default to "good" return values.
self.fitter = unittest.mock.Mock(spec=lsst.jointcal.PhotometryFit)
self.fitter = mock.Mock(spec=lsst.jointcal.PhotometryFit)
self.fitter.computeChi2.return_value = self.goodChi2
self.fitter.minimize.return_value = MinimizeResult.Converged
self.associations = unittest.mock.Mock(spec=lsst.jointcal.Associations)
self.associations.getCcdImageList.return_value = self.ccdImageList
self.model = unittest.mock.Mock(spec=lsst.jointcal.SimpleFluxModel)
self.model = mock.Mock(spec=lsst.jointcal.SimpleFluxModel)

self.config = lsst.jointcal.jointcal.JointcalConfig()
self.jointcal = lsst.jointcal.JointcalTask(config=self.config, butler=self.butler)

def test_iterateFit_success(self):
Expand Down Expand Up @@ -127,7 +130,7 @@ def test_iterateFit_failed(self):
self.assertEqual(self.fitter.minimize.call_count, 1)

def test_iterateFit_badFinalChi2(self):
log = unittest.mock.Mock(spec=lsst.log.Log)
log = mock.Mock(spec=lsst.log.Log)
self.jointcal.log = log
self.fitter.computeChi2.return_value = self.badChi2

Expand All @@ -138,7 +141,7 @@ def test_iterateFit_badFinalChi2(self):
log.error.assert_called_with("Potentially bad fit: High chi-squared/ndof.")

def test_iterateFit_exceedMaxSteps(self):
log = unittest.mock.Mock(spec=lsst.log.Log)
log = mock.Mock(spec=lsst.log.Log)
self.jointcal.log = log
self.fitter.minimize.return_value = MinimizeResult.Chi2Increased
maxSteps = 3
Expand Down Expand Up @@ -178,7 +181,7 @@ def _make_fake_refcat(self):
fakeRefCat = make_fake_refcat(center, flux, filterName)
fluxField = getRefFluxField(fakeRefCat.schema, filterName)
returnStruct = lsst.pipe.base.Struct(refCat=fakeRefCat, fluxField=fluxField)
refObjLoader = unittest.mock.Mock(spec=LoadIndexedReferenceObjectsTask)
refObjLoader = mock.Mock(spec=LoadIndexedReferenceObjectsTask)
refObjLoader.loadSkyCircle.return_value = returnStruct

return refObjLoader, center, radius, filterName, fakeRefCat
Expand Down Expand Up @@ -221,6 +224,49 @@ def test_load_reference_catalog_subselect(self):
self.assertEqual(len(refCat), 0)


class TestJointcalFitModel(JointcalTestBase, lsst.utils.tests.TestCase):
def test_fit_photometry_writeChi2(self):
"""Test that we are calling saveChi2 with appropriate file prefixes."""
self.config.photometryModel = "constrainedFlux"
self.config.writeChi2FilesOuterLoop = True
jointcal = lsst.jointcal.JointcalTask(config=self.config, butler=self.butler)
jointcal.focalPlaneBBox = lsst.geom.Box2D()

# Mock the fitter, so we can pretend it found a good fit
with mock.patch("lsst.jointcal.PhotometryFit", autospect=True) as fitPatch:
fitPatch.return_value.computeChi2.return_value = self.goodChi2
fitPatch.return_value.minimize.return_value = MinimizeResult.Converged

expected = ["photometry_init-ModelVisit_chi2", "photometry_init-Model_chi2",
"photometry_init-Fluxes_chi2", "photometry_init-ModelFluxes_chi2"]
expected = [mock.call(x+"-fake{type}") for x in expected]
jointcal._fit_photometry(self.associations, dataName=self.dataName)
fitPatch.return_value.saveChi2Contributions.assert_has_calls(expected)

def test_fit_astrometry_writeChi2(self):
"""Test that we are calling saveChi2 with appropriate file prefixes."""
self.config.astrometryModel = "constrained"
self.config.writeChi2FilesOuterLoop = True
jointcal = lsst.jointcal.JointcalTask(config=self.config, butler=self.butler)
jointcal.focalPlaneBBox = lsst.geom.Box2D()

# Mock the fitter, so we can pretend it found a good fit
fitPatch = mock.patch("lsst.jointcal.AstrometryFit")
# Mock the projection handler so we don't segfault due to not-fully initialized ccdImages
projectorPatch = mock.patch("lsst.jointcal.OneTPPerVisitHandler")
with fitPatch as fit, projectorPatch as projector:
fit.return_value.computeChi2.return_value = self.goodChi2
fit.return_value.minimize.return_value = MinimizeResult.Converged
# return a real ProjectionHandler to keep ConstrainedAstrometryModel() happy
projector.return_value = lsst.jointcal.IdentityProjectionHandler()

expected = ["astrometry_init-DistortionsVisit_chi2", "astrometry_init-Distortions_chi2",
"astrometry_init-Positions_chi2", "astrometry_init-DistortionsPositions_chi2"]
expected = [mock.call(x+"-fake{type}") for x in expected]
jointcal._fit_astrometry(self.associations, dataName=self.dataName)
fit.return_value.saveChi2Contributions.assert_has_calls(expected)


class MemoryTester(lsst.utils.tests.MemoryTestCase):
pass

Expand Down

0 comments on commit 29627c8

Please sign in to comment.