Coverage for tests / test_cModel.py: 23%

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 unittest 

24import numpy 

25 

26import lsst.utils.tests 

27import lsst.shapelet 

28import lsst.afw.geom 

29import lsst.geom 

30import lsst.afw.image 

31import lsst.utils.logging 

32import lsst.meas.modelfit 

33import lsst.meas.base 

34 

35# Set trace to 0-5 to view debug messages. Level 5 enables all traces. 

36lsst.utils.logging.trace_set_at("lsst.meas.modelfit.optimizer.Optimizer", -1) 

37lsst.utils.logging.trace_set_at("lsst.meas.modelfit.optimizer.solveTrustRegion", -1) 

38 

39 

40def makeMultiShapeletCircularGaussian(sigma): 

41 s = lsst.shapelet.ShapeletFunction(0, lsst.shapelet.HERMITE, sigma) 

42 s.getCoefficients()[0] = 1.0 / lsst.shapelet.ShapeletFunction.FLUX_FACTOR 

43 m = lsst.shapelet.MultiShapeletFunction() 

44 m.addComponent(s) 

45 return m 

46 

47 

48def computePsfFlux(centroid, exposure): 

49 schema = lsst.afw.table.SourceTable.makeMinimalSchema() 

50 pointKey = lsst.afw.table.Point2DKey.addFields(schema, "centroid", "known input centroid", "pixel") 

51 schema.getAliasMap().set("slot_Centroid", "centroid") 

52 algorithm = lsst.meas.base.PsfFluxAlgorithm(lsst.meas.base.PsfFluxControl(), "base_PsfFlux", schema) 

53 table = lsst.afw.table.SourceTable.make(schema) 

54 record = table.makeRecord() 

55 record.set(pointKey, centroid) 

56 algorithm.measure(record, exposure) 

57 return record.get("base_PsfFlux_instFlux"), record.get("base_PsfFlux_instFluxErr") 

58 

59 

60class CModelTestCase(lsst.utils.tests.TestCase): 

61 

62 def setUp(self): 

63 # Setup test data: a single point source, initially with no noise. 

64 numpy.random.seed(500) 

65 crval = lsst.geom.SpherePoint(45.0, 45.0, lsst.geom.degrees) 

66 crpix = lsst.geom.Point2D(0.0, 0.0) 

67 scale = 0.2 * lsst.geom.arcseconds 

68 cdMatrix = lsst.afw.geom.makeCdMatrix(scale=scale, flipX=True) 

69 dataWcs = lsst.afw.geom.makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix) 

70 photoCalib = lsst.afw.image.PhotoCalib(4.0) 

71 self.xyPosition = lsst.geom.Point2D(1.1, -0.8) 

72 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-100, -100), lsst.geom.Point2I(100, 100)) 

73 self.exposure = lsst.afw.image.ExposureF(bbox) 

74 self.exposure.setWcs(dataWcs) 

75 self.exposure.setPhotoCalib(photoCalib) 

76 self.trueFlux = 65.0 

77 self.psfSigma = 2.0 

78 psf = lsst.afw.detection.GaussianPsf(25, 25, self.psfSigma) 

79 self.exposure.setPsf(psf) 

80 psfImage = psf.computeImage(self.xyPosition) 

81 psfImage.getArray()[:, :] *= self.trueFlux 

82 psfBBox = psfImage.getBBox(lsst.afw.image.PARENT) 

83 subImage = lsst.afw.image.ImageF(self.exposure.getMaskedImage().getImage(), psfBBox, 

84 lsst.afw.image.PARENT) 

85 subImage.getArray()[:, :] = psfImage.getArray() 

86 

87 def tearDown(self): 

88 del self.xyPosition 

89 del self.exposure 

90 del self.trueFlux 

91 del self.psfSigma 

92 

93 def testNoNoise(self): 

94 """Test that CModelAlgorithm.apply() works when applied to a postage-stamp 

95 containing only a point source with no noise. 

96 

97 We still have to pretend there is noise (i.e. have nonzero values in 

98 the variance plane) to allow it to compute a likelihood, though. 

99 """ 

100 ctrl = lsst.meas.modelfit.CModelControl() 

101 ctrl.initial.usePixelWeights = False 

102 algorithm = lsst.meas.modelfit.CModelAlgorithm(ctrl) 

103 var = 1E-16 

104 self.exposure.getMaskedImage().getVariance().getArray()[:, :] = var 

105 psfImage = self.exposure.getPsf().computeKernelImage(self.xyPosition).getArray() 

106 expectedFluxErr = var**0.5 * (psfImage**2).sum()**(-0.5) 

107 pos = self.exposure.getPsf().getAveragePosition() 

108 result = algorithm.apply( 

109 self.exposure, makeMultiShapeletCircularGaussian(self.psfSigma), 

110 self.xyPosition, self.exposure.getPsf().computeShape(pos) 

111 ) 

112 self.assertFalse(result.initial.flags[result.FAILED]) 

113 self.assertFloatsAlmostEqual(result.initial.instFlux, self.trueFlux, rtol=0.01) 

114 self.assertFloatsAlmostEqual(result.initial.instFluxErr, expectedFluxErr, rtol=0.01) 

115 self.assertLess(result.initial.ellipse.getDeterminantRadius(), 0.2) 

116 self.assertFalse(result.exp.flags[result.FAILED]) 

117 self.assertFloatsAlmostEqual(result.exp.instFlux, self.trueFlux, rtol=0.01) 

118 self.assertFloatsAlmostEqual(result.exp.instFluxErr, expectedFluxErr, rtol=0.01) 

119 self.assertLess(result.exp.ellipse.getDeterminantRadius(), 0.2) 

120 self.assertFalse(result.dev.flags[result.FAILED]) 

121 self.assertFloatsAlmostEqual(result.dev.instFlux, self.trueFlux, rtol=0.01) 

122 self.assertFloatsAlmostEqual(result.dev.instFluxErr, expectedFluxErr, rtol=0.01) 

123 self.assertLess(result.dev.ellipse.getDeterminantRadius(), 0.2) 

124 self.assertFalse(result.flags[result.FAILED]) 

125 self.assertFloatsAlmostEqual(result.instFlux, self.trueFlux, rtol=0.01) 

126 

127 def testVsPsfFlux(self): 

128 """Test that CModel produces results comparable to PsfFlux when run 

129 on point sources. 

130 """ 

131 noiseSigma = 1.0 

132 for fluxFactor in (1.0, 10.0, 100.0): 

133 exposure = self.exposure.Factory(self.exposure, True) 

134 exposure.getMaskedImage().getImage().getArray()[:] *= fluxFactor 

135 exposure.getMaskedImage().getVariance().getArray()[:] = noiseSigma**2 

136 exposure.getMaskedImage().getImage().getArray()[:] += \ 

137 noiseSigma*numpy.random.randn(exposure.getHeight(), exposure.getWidth()) 

138 ctrl = lsst.meas.modelfit.CModelControl() 

139 algorithm = lsst.meas.modelfit.CModelAlgorithm(ctrl) 

140 pos = self.exposure.getPsf().getAveragePosition() 

141 cmodel = algorithm.apply( 

142 exposure, makeMultiShapeletCircularGaussian(self.psfSigma), 

143 self.xyPosition, self.exposure.getPsf().computeShape(pos) 

144 ) 

145 psfFlux, psfFluxErr = computePsfFlux(self.xyPosition, exposure) 

146 self.assertFloatsAlmostEqual(psfFlux, cmodel.instFlux, rtol=0.1/fluxFactor**0.5) 

147 self.assertFloatsAlmostEqual(psfFluxErr, cmodel.instFluxErr, rtol=0.1/fluxFactor**0.5) 

148 

149 

150class TestMemory(lsst.utils.tests.MemoryTestCase): 

151 pass 

152 

153 

154def setup_module(module): 

155 lsst.utils.tests.init() 

156 

157 

158if __name__ == "__main__": 158 ↛ 159line 158 didn't jump to line 159 because the condition on line 158 was never true

159 lsst.utils.tests.init() 

160 unittest.main()