Coverage for tests / test_doubleShapeletPsfApprox.py: 16%
262 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:12 +0000
1#
2# LSST Data Management System
3#
4# Copyright 2008-2016 AURA/LSST.
5#
6# This product includes software developed by the
7# LSST Project (http://www.lsst.org/).
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the LSST License Statement and
20# the GNU General Public License along with this program. If not,
21# see <https://www.lsstcorp.org/LegalNotices/>.
22#
23import os
24import unittest
25import numpy
26from io import StringIO
27import warnings
29import lsst.utils.tests
30import lsst.afw.detection
31import lsst.afw.image
32import lsst.geom
33import lsst.afw.geom
34import lsst.afw.geom.ellipses
35import lsst.utils.logging
36import lsst.meas.modelfit
37import lsst.meas.algorithms
39# Set trace to 0-5 to view debug messages. Level 5 enables all traces.
40lsst.utils.logging.trace_set_at("lsst.meas.modelfit.optimizer.Optimizer", -1)
41lsst.utils.logging.trace_set_at("lsst.meas.modelfit.optimizer.solveTrustRegion", -1)
44class DoubleShapeletPsfApproxTestMixin:
46 Algorithm = lsst.meas.modelfit.DoubleShapeletPsfApproxAlgorithm
48 def initialize(self, psf, ctrl=None, atol=1E-4, **kwds):
49 if not isinstance(psf, lsst.afw.detection.Psf):
50 kernel = lsst.afw.math.FixedKernel(psf)
51 psf = lsst.meas.algorithms.KernelPsf(kernel)
52 self.psf = psf
53 self.atol = atol
54 if ctrl is None:
55 ctrl = lsst.meas.modelfit.DoubleShapeletPsfApproxControl()
56 self.ctrl = ctrl
57 for name, value in kwds.items():
58 setattr(self.ctrl, name, value)
59 self.exposure = lsst.afw.image.ExposureF(1, 1)
60 scale = 5.0e-5 * lsst.geom.degrees
61 wcs = lsst.afw.geom.makeSkyWcs(crpix=lsst.geom.Point2D(0.0, 0.0),
62 crval=lsst.geom.SpherePoint(45, 45, lsst.geom.degrees),
63 cdMatrix=lsst.afw.geom.makeCdMatrix(scale=scale))
64 self.exposure.setWcs(wcs)
65 self.exposure.setPsf(self.psf)
67 def tearDown(self):
68 del self.exposure
69 del self.psf
70 del self.ctrl
71 del self.atol
73 def setupTaskConfig(self, config):
74 config.slots.shape = None
75 config.slots.psfFlux = None
76 config.slots.apFlux = None
77 config.slots.gaussianFlux = None
78 config.slots.modelFlux = None
79 config.slots.calibFlux = None
80 config.doReplaceWithNoise = False
81 config.plugins.names = ["modelfit_DoubleShapeletPsfApprox"]
82 config.plugins["modelfit_DoubleShapeletPsfApprox"].readControl(self.ctrl)
84 def checkBounds(self, msf):
85 """Check that the bounds specified in the control object are met by a MultiShapeletFunction.
87 These requirements must be true after a call to any fit method or measure().
88 """
89 pos = self.psf.getAveragePosition()
90 self.assertEqual(len(msf.getComponents()), 2)
91 self.assertEqual(
92 lsst.shapelet.computeSize(self.ctrl.innerOrder),
93 len(msf.getComponents()[0].getCoefficients())
94 )
95 self.assertEqual(
96 lsst.shapelet.computeSize(self.ctrl.outerOrder),
97 len(msf.getComponents()[1].getCoefficients())
98 )
99 self.assertGreater(
100 self.ctrl.maxRadiusBoxFraction * (self.psf.computeKernelImage(pos).getBBox().getArea())**0.5,
101 lsst.afw.geom.ellipses.Axes(msf.getComponents()[0].getEllipse().getCore()).getA()
102 )
103 self.assertGreater(
104 self.ctrl.maxRadiusBoxFraction * (self.psf.computeKernelImage(pos).getBBox().getArea())**0.5,
105 lsst.afw.geom.ellipses.Axes(msf.getComponents()[1].getEllipse().getCore()).getA()
106 )
107 self.assertLess(
108 self.ctrl.minRadius,
109 lsst.afw.geom.ellipses.Axes(msf.getComponents()[0].getEllipse().getCore()).getB()
110 )
111 self.assertLess(
112 self.ctrl.minRadius,
113 lsst.afw.geom.ellipses.Axes(msf.getComponents()[1].getEllipse().getCore()).getB()
114 )
115 self.assertLess(
116 self.ctrl.minRadiusDiff,
117 (msf.getComponents()[1].getEllipse().getCore().getDeterminantRadius()
118 - msf.getComponents()[0].getEllipse().getCore().getDeterminantRadius())
119 )
121 def checkRatios(self, msf):
122 """Check that the ratios specified in the control object are met by a MultiShapeletFunction.
124 These requirements must be true after initializeResult and fitMoments, but will are relaxed
125 in later stages of the fit.
126 """
127 inner = msf.getComponents()[0]
128 outer = msf.getComponents()[1]
129 position = msf.getComponents()[0].getEllipse().getCenter()
130 self.assertFloatsAlmostEqual(position.getX(), msf.getComponents()[1].getEllipse().getCenter().getX())
131 self.assertFloatsAlmostEqual(position.getY(), msf.getComponents()[1].getEllipse().getCenter().getY())
132 self.assertFloatsAlmostEqual(outer.evaluate()(position),
133 inner.evaluate()(position)*self.ctrl.peakRatio)
134 self.assertFloatsAlmostEqual(
135 outer.getEllipse().getCore().getDeterminantRadius(),
136 inner.getEllipse().getCore().getDeterminantRadius() * self.ctrl.radiusRatio
137 )
139 def makeImages(self, msf):
140 """Return an Image of the data and an Image of the model for comparison.
141 """
142 pos = self.exposure.getPsf().getAveragePosition()
143 dataImage = self.exposure.getPsf().computeKernelImage(pos)
144 modelImage = dataImage.Factory(dataImage.getBBox())
145 msf.evaluate().addToImage(modelImage)
146 return dataImage, modelImage
148 def checkFitQuality(self, msf):
149 """Check the quality of the fit by comparing to the PSF image.
150 """
151 dataImage, modelImage = self.makeImages(msf)
152 self.assertFloatsAlmostEqual(dataImage.getArray(), modelImage.getArray(), atol=self.atol,
153 plotOnFailure=True)
155 def testSingleFramePlugin(self):
156 """Run the algorithm as a single-frame plugin and check the quality of the fit.
157 """
158 with warnings.catch_warnings():
159 warnings.filterwarnings("ignore", message="ignoreSlotPluginChecks", category=FutureWarning)
160 config = lsst.meas.base.SingleFrameMeasurementTask.ConfigClass(ignoreSlotPluginChecks=True)
161 self.setupTaskConfig(config)
162 config.slots.centroid = "centroid"
163 schema = lsst.afw.table.SourceTable.makeMinimalSchema()
164 centroidKey = lsst.afw.table.Point2DKey.addFields(schema, "centroid", "centroid", "pixel")
165 task = lsst.meas.base.SingleFrameMeasurementTask(config=config, schema=schema)
166 measCat = lsst.afw.table.SourceCatalog(schema)
167 measRecord = measCat.addNew()
168 measRecord.set(centroidKey, lsst.geom.Point2D(0.0, 0.0))
169 task.run(measCat, self.exposure)
170 self.assertFalse(measRecord.get("modelfit_DoubleShapeletPsfApprox_flag"))
171 key = lsst.shapelet.MultiShapeletFunctionKey(schema["modelfit"]["DoubleShapeletPsfApprox"])
172 msf = measRecord.get(key)
173 self.checkBounds(msf)
174 self.checkFitQuality(msf)
176 def testForcedPlugin(self):
177 """Run the algorithm as a forced plugin and check the quality of the fit.
178 """
179 config = lsst.meas.base.ForcedMeasurementTask.ConfigClass()
180 config.copyColumns = {"id": "objectId", "parent": "parentObjectId"}
181 self.setupTaskConfig(config)
182 config.slots.centroid = "base_TransformedCentroid"
183 config.plugins.names |= ["base_TransformedCentroid"]
184 refSchema = lsst.afw.table.SourceTable.makeMinimalSchema()
185 refCentroidKey = lsst.afw.table.Point2DKey.addFields(refSchema, "centroid", "centroid", "pixel")
186 refSchema.getAliasMap().set("slot_Centroid", "centroid")
187 refCat = lsst.afw.table.SourceCatalog(refSchema)
188 refRecord = refCat.addNew()
189 refRecord.set(refCentroidKey, lsst.geom.Point2D(0.0, 0.0))
190 refWcs = self.exposure.getWcs() # same as measurement Wcs
191 task = lsst.meas.base.ForcedMeasurementTask(config=config, refSchema=refSchema)
192 measCat = task.generateMeasCat(self.exposure, refCat, refWcs)
193 task.run(measCat, self.exposure, refCat, refWcs)
194 measRecord = measCat[0]
195 self.assertFalse(measRecord.get("modelfit_DoubleShapeletPsfApprox_flag"))
196 measSchema = measCat.schema
197 key = lsst.shapelet.MultiShapeletFunctionKey(measSchema["modelfit"]["DoubleShapeletPsfApprox"])
198 msf = measRecord.get(key)
199 self.checkBounds(msf)
200 self.checkFitQuality(msf)
202 def testInitializeResult(self):
203 """Test that initializeResult() returns a unit-flux, unit-circle MultiShapeletFunction
204 with the right peakRatio and radiusRatio.
205 """
206 msf = self.Algorithm.initializeResult(self.ctrl)
207 self.assertFloatsAlmostEqual(msf.evaluate().integrate(), 1.0)
208 moments = msf.evaluate().computeMoments()
209 axes = lsst.afw.geom.ellipses.Axes(moments.getCore())
210 self.assertFloatsAlmostEqual(moments.getCenter().getX(), 0.0)
211 self.assertFloatsAlmostEqual(moments.getCenter().getY(), 0.0)
212 self.assertFloatsAlmostEqual(axes.getA(), 1.0)
213 self.assertFloatsAlmostEqual(axes.getB(), 1.0)
214 self.assertEqual(len(msf.getComponents()), 2)
215 self.checkRatios(msf)
217 def testFitMoments(self):
218 """Test that fitMoments() preserves peakRatio and radiusRatio while setting moments
219 correctly.
220 """
221 MOMENTS_RTOL = 1E-13
222 image = self.psf.computeKernelImage(self.psf.getAveragePosition())
223 array = image.getArray()
224 bbox = image.getBBox()
225 x, y = numpy.meshgrid(
226 numpy.arange(bbox.getBeginX(), bbox.getEndX()),
227 numpy.arange(bbox.getBeginY(), bbox.getEndY())
228 )
229 msf = self.Algorithm.initializeResult(self.ctrl)
230 self.Algorithm.fitMoments(msf, self.ctrl, image)
231 self.assertFloatsAlmostEqual(msf.evaluate().integrate(), array.sum(), rtol=MOMENTS_RTOL)
232 moments = msf.evaluate().computeMoments()
233 q = lsst.afw.geom.ellipses.Quadrupole(moments.getCore())
234 cx = (x*array).sum()/array.sum()
235 cy = (y*array).sum()/array.sum()
236 self.assertFloatsAlmostEqual(moments.getCenter().getX(), cx, rtol=MOMENTS_RTOL)
237 self.assertFloatsAlmostEqual(moments.getCenter().getY(), cy, rtol=MOMENTS_RTOL)
238 self.assertFloatsAlmostEqual(q.getIxx(), ((x - cx)**2 * array).sum()/array.sum(), rtol=MOMENTS_RTOL)
239 self.assertFloatsAlmostEqual(q.getIyy(), ((y - cy)**2 * array).sum()/array.sum(), rtol=MOMENTS_RTOL)
240 self.assertFloatsAlmostEqual(q.getIxy(), ((x - cx)*(y - cy)*array).sum()/array.sum(),
241 rtol=MOMENTS_RTOL)
242 self.assertEqual(len(msf.getComponents()), 2)
243 self.checkRatios(msf)
244 self.checkBounds(msf)
246 def testObjective(self):
247 """Test that model evaluation agrees with derivative evaluation in the objective object.
248 """
249 image = self.psf.computeKernelImage(self.psf.getAveragePosition())
250 msf = self.Algorithm.initializeResult(self.ctrl)
251 self.Algorithm.fitMoments(msf, self.ctrl, image)
252 moments = msf.evaluate().computeMoments()
253 r0 = moments.getCore().getDeterminantRadius()
254 objective = self.Algorithm.makeObjective(moments, self.ctrl, image)
255 image, model = self.makeImages(msf)
256 parameters = numpy.zeros(4, dtype=float)
257 parameters[0] = msf.getComponents()[0].getCoefficients()[0]
258 parameters[1] = msf.getComponents()[1].getCoefficients()[0]
259 parameters[2] = msf.getComponents()[0].getEllipse().getCore().getDeterminantRadius() / r0
260 parameters[3] = msf.getComponents()[1].getEllipse().getCore().getDeterminantRadius() / r0
261 residuals = numpy.zeros(image.getArray().size, dtype=float)
262 objective.computeResiduals(parameters, residuals)
263 self.assertFloatsAlmostEqual(
264 residuals.reshape(image.getHeight(), image.getWidth()),
265 image.getArray() - model.getArray()
266 )
267 step = 1E-6
268 derivatives = numpy.zeros((parameters.size, residuals.size), dtype=float).transpose()
269 objective.differentiateResiduals(parameters, derivatives)
270 for i in range(parameters.size):
271 original = parameters[i]
272 r1 = numpy.zeros(residuals.size, dtype=float)
273 r2 = numpy.zeros(residuals.size, dtype=float)
274 parameters[i] = original + step
275 objective.computeResiduals(parameters, r1)
276 parameters[i] = original - step
277 objective.computeResiduals(parameters, r2)
278 parameters[i] = original
279 d = (r1 - r2)/(2.0*step)
280 self.assertFloatsAlmostEqual(
281 d.reshape(image.getHeight(), image.getWidth()),
282 derivatives[:, i].reshape(image.getHeight(), image.getWidth()),
283 atol=1E-11
284 )
286 def testFitProfile(self):
287 """Test that fitProfile() does not modify the ellipticity, that it improves the fit, and
288 that small perturbations to the zeroth-order amplitudes and radii do not improve the fit.
289 """
290 image = self.psf.computeKernelImage(self.psf.getAveragePosition())
291 msf = self.Algorithm.initializeResult(self.ctrl)
292 self.Algorithm.fitMoments(msf, self.ctrl, image)
293 prev = lsst.shapelet.MultiShapeletFunction(msf)
294 self.Algorithm.fitProfile(msf, self.ctrl, image)
296 def getEllipticity(m, c):
297 s = lsst.afw.geom.ellipses.SeparableDistortionDeterminantRadius(
298 m.getComponents()[c].getEllipse().getCore()
299 )
300 return numpy.array([s.getE1(), s.getE2()])
301 self.assertFloatsAlmostEqual(getEllipticity(prev, 0), getEllipticity(msf, 0), rtol=1E-13)
302 self.assertFloatsAlmostEqual(getEllipticity(prev, 1), getEllipticity(msf, 1), rtol=1E-13)
304 def computeChiSq(m):
305 data, model = self.makeImages(m)
306 return numpy.sum((data.getArray() - model.getArray())**2)
307 bestChiSq = computeChiSq(msf)
308 self.assertLessEqual(bestChiSq, computeChiSq(prev))
309 step = 1E-4
310 for component in msf.getComponents():
311 # 0th-order amplitude perturbation
312 original = component.getCoefficients()[0]
313 component.getCoefficients()[0] = original + step
314 self.assertLessEqual(bestChiSq, computeChiSq(msf))
315 component.getCoefficients()[0] = original - step
316 self.assertLessEqual(bestChiSq, computeChiSq(msf))
317 component.getCoefficients()[0] = original
318 # Radius perturbation
319 original = component.getEllipse()
320 component.getEllipse().getCore().scale(1.0 + step)
321 self.assertLessEqual(bestChiSq, computeChiSq(msf))
322 component.setEllipse(original)
323 component.getEllipse().getCore().scale(1.0 - step)
324 self.assertLessEqual(bestChiSq, computeChiSq(msf))
325 component.setEllipse(original)
327 def testFitShapelets(self):
328 """Test that fitShapelets() does not modify the zeroth order coefficients or ellipse,
329 that it improves the fit, and that small perturbations to the higher-order coefficients
330 do not improve the fit.
331 """
332 image = self.psf.computeKernelImage(self.psf.getAveragePosition())
333 msf = self.Algorithm.initializeResult(self.ctrl)
334 self.Algorithm.fitMoments(msf, self.ctrl, image)
335 self.Algorithm.fitProfile(msf, self.ctrl, image)
336 prev = lsst.shapelet.MultiShapeletFunction(msf)
337 self.Algorithm.fitShapelets(msf, self.ctrl, image)
338 self.assertFloatsAlmostEqual(
339 prev.getComponents()[0].getEllipse().getParameterVector(),
340 msf.getComponents()[0].getEllipse().getParameterVector()
341 )
342 self.assertFloatsAlmostEqual(
343 prev.getComponents()[1].getEllipse().getParameterVector(),
344 msf.getComponents()[1].getEllipse().getParameterVector()
345 )
347 def computeChiSq(m):
348 data, model = self.makeImages(m)
349 return numpy.sum((data.getArray() - model.getArray())**2)
350 bestChiSq = computeChiSq(msf)
351 self.assertLessEqual(bestChiSq, computeChiSq(prev))
352 step = 1E-4
353 for component in msf.getComponents():
354 for i in range(1, len(component.getCoefficients())):
355 original = component.getCoefficients()[i]
356 component.getCoefficients()[i] = original + step
357 self.assertLessEqual(bestChiSq, computeChiSq(msf))
358 component.getCoefficients()[i] = original - step
359 self.assertLessEqual(bestChiSq, computeChiSq(msf))
360 component.getCoefficients()[i] = original
362 def testSingleFrameConfigIO(self):
363 config1 = lsst.meas.base.SingleFrameMeasurementTask.ConfigClass()
364 config2 = lsst.meas.base.SingleFrameMeasurementTask.ConfigClass()
365 self.setupTaskConfig(config1)
366 stream = StringIO()
367 config1.saveToStream(stream)
368 config2.loadFromString(stream.getvalue())
369 self.assertEqual(config1, config2)
372class SingleGaussianTestCase(DoubleShapeletPsfApproxTestMixin, lsst.utils.tests.TestCase):
374 def setUp(self):
375 numpy.random.seed(500)
376 DoubleShapeletPsfApproxTestMixin.initialize(
377 self, psf=lsst.afw.detection.GaussianPsf(25, 25, 2.0),
378 innerOrder=0, outerOrder=0, peakRatio=0.0
379 )
382class HigherOrderTestCase0(DoubleShapeletPsfApproxTestMixin, lsst.utils.tests.TestCase):
384 def setUp(self):
385 numpy.random.seed(500)
386 image = lsst.afw.image.ImageD(os.path.join(os.path.dirname(os.path.realpath(__file__)),
387 "data", "psfs/great3-0.fits"))
388 DoubleShapeletPsfApproxTestMixin.initialize(
389 self, psf=image,
390 innerOrder=3, outerOrder=2,
391 atol=0.0005
392 )
395class HigherOrderTestCase1(DoubleShapeletPsfApproxTestMixin, lsst.utils.tests.TestCase):
397 def setUp(self):
398 numpy.random.seed(500)
399 image = lsst.afw.image.ImageD(os.path.join(os.path.dirname(os.path.realpath(__file__)),
400 "data", "psfs/great3-1.fits"))
401 DoubleShapeletPsfApproxTestMixin.initialize(
402 self, psf=image,
403 innerOrder=2, outerOrder=1,
404 atol=0.002
405 )
408class TestMemory(lsst.utils.tests.MemoryTestCase):
409 pass
412def setup_module(module):
413 lsst.utils.tests.init()
416if __name__ == "__main__": 416 ↛ 417line 416 didn't jump to line 417 because the condition on line 416 was never true
417 lsst.utils.tests.init()
418 unittest.main()