Coverage for tests / test_psfFitter.py: 15%

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 os 

25import glob 

26import numpy 

27 

28import lsst.utils.tests 

29import lsst.shapelet 

30import lsst.geom 

31import lsst.utils.logging 

32import lsst.meas.modelfit 

33import lsst.meas.base 

34 

35numpy.random.seed(500) 

36 

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

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

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

40 

41ELLIPSE_PARAMETER_NAMES = ["eta1", "eta2", "logR", "x", "y"] 

42DATA_DIR = os.path.join(os.environ["MEAS_MODELFIT_DIR"], "tests", "data") 

43 

44 

45def computeMoments(image): 

46 """Helper function to compute moments of a postage stamp about its origin.""" 

47 maskedImage = lsst.afw.image.MaskedImageD(image) 

48 result = lsst.meas.base.SdssShapeAlgorithm.computeAdaptiveMoments( 

49 maskedImage, 

50 lsst.geom.Point2D(0.0, 0.0) 

51 ) 

52 return result.getShape() 

53 

54 

55class GeneralPsfFitterTestCase(lsst.utils.tests.TestCase): 

56 

57 def setUp(self): 

58 self.configs = {} 

59 self.configs['fixed'] = lsst.meas.modelfit.GeneralPsfFitterConfig() 

60 self.configs['fixed'].primary.ellipticityPriorSigma = 0.0 

61 self.configs['fixed'].primary.radiusPriorSigma = 0.0 

62 self.configs['fixed'].primary.positionPriorSigma = 0.0 

63 self.configs['fixed'].wings.ellipticityPriorSigma = 0.0 

64 self.configs['fixed'].wings.radiusPriorSigma = 0.0 

65 self.configs['fixed'].wings.positionPriorSigma = 0.0 

66 self.configs['ellipse'] = lsst.meas.modelfit.GeneralPsfFitterConfig() 

67 self.configs['ellipse'].primary.positionPriorSigma = 0.0 

68 self.configs['ellipse'].wings.positionPriorSigma = 0.0 

69 self.configs['full'] = lsst.meas.modelfit.GeneralPsfFitterConfig() 

70 self.configs['full'].inner.order = 0 

71 self.configs['full'].primary.order = 4 

72 self.configs['full'].wings.order = 4 

73 self.configs['full'].outer.order = 0 

74 self.configs['full'].inner.ellipticityPriorSigma = 0.3 

75 self.configs['full'].inner.radiusPriorSigma = 0.5 

76 self.configs['full'].inner.positionPriorSigma = 0.1 

77 self.configs['full'].primary.ellipticityPriorSigma = 0.3 

78 self.configs['full'].primary.radiusPriorSigma = 0.5 

79 self.configs['full'].primary.positionPriorSigma = 0.1 

80 self.configs['full'].wings.ellipticityPriorSigma = 0.3 

81 self.configs['full'].wings.radiusPriorSigma = 0.5 

82 self.configs['full'].wings.positionPriorSigma = 0.1 

83 self.configs['full'].outer.ellipticityPriorSigma = 0.3 

84 self.configs['full'].outer.radiusPriorSigma = 0.5 

85 self.configs['full'].outer.positionPriorSigma = 0.1 

86 

87 def tearDown(self): 

88 del self.configs 

89 

90 def testFixedModel(self): 

91 fitter = lsst.meas.modelfit.GeneralPsfFitter(self.configs['fixed'].makeControl()) 

92 model = fitter.getModel() 

93 

94 # check that we have the right numbers and names for parameters 

95 self.assertEqual(model.getNonlinearDim(), 0) 

96 self.assertEqual(model.getFixedDim(), 10) 

97 self.assertEqual(model.getAmplitudeDim(), 2) 

98 self.assertEqual(model.getBasisCount(), 2) 

99 self.assertEqual(list(model.getNonlinearNames()), []) 

100 self.assertEqual(list(model.getAmplitudeNames()), ["primary.alpha[0,0]", "wings.alpha[0,0]"]) 

101 self.assertEqual(list(model.getFixedNames()), 

102 [f"primary.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

103 + [f"wings.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES]) 

104 

105 # test that we can round-trip ellipses through the model, and that this agrees 

106 # with makeShapeletFunction 

107 ellipseParameters = numpy.array([[0.01, -0.01, 1.1, 0.03, -0.04], 

108 [0.02, -0.02, 0.9, 0.05, -0.06]]) 

109 ellipses1 = model.makeEllipseVector() 

110 for i in range(len(ellipses1)): 

111 ellipses1[i].setParameterVector(ellipseParameters[i]) 

112 nonlinear = numpy.zeros(model.getNonlinearDim(), dtype=lsst.meas.modelfit.Scalar) 

113 fixed = numpy.zeros(model.getFixedDim(), dtype=lsst.meas.modelfit.Scalar) 

114 amplitudes = numpy.array([1.0, 0.1], dtype=lsst.meas.modelfit.Scalar) 

115 model.readEllipses(ellipses1, nonlinear, fixed) 

116 self.assertFloatsAlmostEqual(fixed, ellipseParameters.ravel()) 

117 ellipses2 = model.writeEllipses(nonlinear, fixed) 

118 msf = model.makeShapeletFunction(nonlinear, amplitudes, fixed) 

119 self.assertFloatsAlmostEqual(len(msf.getComponents()), len(ellipses1)) 

120 ellipses3 = model.makeEllipseVector() 

121 for i in range(len(ellipses2)): 

122 self.assertFloatsAlmostEqual(ellipses1[i].getParameterVector(), ellipses2[i].getParameterVector()) 

123 # need to convert ellipse parametrization 

124 ellipses3[i].setCore(msf.getComponents()[i].getEllipse().getCore()) 

125 ellipses3[i].setCenter(msf.getComponents()[i].getEllipse().getCenter()) 

126 self.assertFloatsAlmostEqual(ellipses1[i].getParameterVector(), ellipses3[i].getParameterVector()) 

127 self.assertFloatsAlmostEqual(amplitudes[i:i+1], msf.getComponents()[i].getCoefficients()) 

128 

129 def testEllipseModel(self): 

130 fitter = lsst.meas.modelfit.GeneralPsfFitter(self.configs['ellipse'].makeControl()) 

131 model = fitter.getModel() 

132 

133 # check that we have the right numbers and names for parameters 

134 self.assertEqual(model.getNonlinearDim(), 6) 

135 self.assertEqual(model.getFixedDim(), 10) 

136 self.assertEqual(model.getAmplitudeDim(), 2) 

137 self.assertEqual(model.getBasisCount(), 2) 

138 self.assertEqual(list(model.getNonlinearNames()), 

139 [f"primary.{s}" for s in ELLIPSE_PARAMETER_NAMES[:3]] 

140 + [f"wings.{s}" for s in ELLIPSE_PARAMETER_NAMES[:3]] 

141 ) 

142 self.assertEqual(list(model.getAmplitudeNames()), ["primary.alpha[0,0]", "wings.alpha[0,0]"]) 

143 self.assertEqual(list(model.getFixedNames()), 

144 [f"primary.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

145 + [f"wings.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES]) 

146 

147 # test that we can round-trip ellipses through the model, and that this agrees 

148 # with makeShapeletFunction 

149 ellipseParameters = numpy.array([[0.01, -0.01, 1.1, 0.03, -0.04], 

150 [0.02, -0.02, 0.9, 0.05, -0.06]]) 

151 ellipses1 = model.makeEllipseVector() 

152 for i in range(len(ellipses1)): 

153 ellipses1[i].setParameterVector(ellipseParameters[i]) 

154 nonlinear = numpy.zeros(model.getNonlinearDim(), dtype=lsst.meas.modelfit.Scalar) 

155 fixed = numpy.zeros(model.getFixedDim(), dtype=lsst.meas.modelfit.Scalar) 

156 amplitudes = numpy.array([1.0, 0.1], dtype=lsst.meas.modelfit.Scalar) 

157 model.readEllipses(ellipses1, nonlinear, fixed) 

158 self.assertFloatsAlmostEqual(nonlinear, numpy.zeros(model.getNonlinearDim(), 

159 dtype=lsst.meas.modelfit.Scalar)) 

160 self.assertFloatsAlmostEqual(fixed, ellipseParameters.ravel()) 

161 ellipses2 = model.writeEllipses(nonlinear, fixed) 

162 msf = model.makeShapeletFunction(nonlinear, amplitudes, fixed) 

163 self.assertEqual(len(msf.getComponents()), len(ellipses1)) 

164 ellipses3 = model.makeEllipseVector() 

165 for i in range(len(ellipses2)): 

166 self.assertFloatsAlmostEqual(ellipses1[i].getParameterVector(), ellipses2[i].getParameterVector(), 

167 rtol=1E-8) 

168 # need to convert ellipse parametrization 

169 ellipses3[i].setCore(msf.getComponents()[i].getEllipse().getCore()) 

170 ellipses3[i].setCenter(msf.getComponents()[i].getEllipse().getCenter()) 

171 self.assertFloatsAlmostEqual(ellipses1[i].getParameterVector(), ellipses3[i].getParameterVector(), 

172 rtol=1E-8) 

173 self.assertFloatsAlmostEqual(amplitudes[i:i+1], msf.getComponents()[i].getCoefficients(), 

174 rtol=1E-8) 

175 

176 # test the ellipse round-tripping again, this time starting with nonzero nonlinear parameters: 

177 # this will be read back in by adding to the fixed parameters and zeroing the nonlinear parameters. 

178 nonlinear[:] = 0.5*ellipseParameters[:, :3].ravel() 

179 ellipses4 = model.writeEllipses(nonlinear, fixed) 

180 model.readEllipses(ellipses4, nonlinear, fixed) 

181 self.assertFloatsAlmostEqual(nonlinear, numpy.zeros(model.getNonlinearDim(), 

182 dtype=lsst.meas.modelfit.Scalar), 

183 rtol=1E-8) 

184 self.assertFloatsAlmostEqual(fixed.reshape(2, 5)[:, :3], 1.5*ellipseParameters[:, :3], rtol=1E-8) 

185 self.assertFloatsAlmostEqual(fixed.reshape(2, 5)[:, 3:], ellipseParameters[:, 3:], rtol=1E-8) 

186 

187 def testFullModel(self): 

188 fitter = lsst.meas.modelfit.GeneralPsfFitter(self.configs['full'].makeControl()) 

189 model = fitter.getModel() 

190 

191 # check that we have the right numbers and names for parameters 

192 self.assertEqual(model.getNonlinearDim(), 20) 

193 self.assertEqual(model.getFixedDim(), 20) 

194 self.assertEqual(model.getAmplitudeDim(), 2*(1 + lsst.shapelet.computeSize(4))) 

195 self.assertEqual(model.getBasisCount(), 4) 

196 self.assertEqual(list(model.getNonlinearNames()), 

197 [f"inner.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

198 + [f"primary.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

199 + [f"wings.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

200 + [f"outer.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

201 ) 

202 self.assertEqual(list(model.getAmplitudeNames()), 

203 ["inner.alpha[0,0]"] 

204 + [f"primary.alpha[{x},{y}]" 

205 for n, x, y in lsst.shapelet.HermiteIndexGenerator(4)] 

206 + [f"wings.alpha[{x},{y}]" 

207 for n, x, y in lsst.shapelet.HermiteIndexGenerator(4)] 

208 + ["outer.alpha[0,0]"]) 

209 self.assertEqual(list(model.getFixedNames()), 

210 [f"inner.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

211 + [f"primary.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

212 + [f"wings.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

213 + [f"outer.fiducial.{s}" for s in ELLIPSE_PARAMETER_NAMES] 

214 ) 

215 

216 # test that we can round-trip ellipses through the model, and that this agrees 

217 # with makeShapeletFunction 

218 ellipseParameters = numpy.array([[0.01, -0.01, 1.1, 0.03, -0.04], 

219 [0.015, -0.015, 1.0, 0.04, -0.05], 

220 [0.02, -0.02, 0.9, 0.05, -0.06], 

221 [0.025, -0.025, 0.8, 0.06, -0.07], 

222 ]) 

223 ellipses1 = model.makeEllipseVector() 

224 for i in range(len(ellipses1)): 

225 ellipses1[i].setParameterVector(ellipseParameters[i]) 

226 nonlinear = numpy.zeros(model.getNonlinearDim(), dtype=lsst.meas.modelfit.Scalar) 

227 fixed = numpy.zeros(model.getFixedDim(), dtype=lsst.meas.modelfit.Scalar) 

228 amplitudes = numpy.random.randn(model.getAmplitudeDim()) 

229 model.readEllipses(ellipses1, nonlinear, fixed) 

230 self.assertFloatsAlmostEqual(nonlinear, numpy.zeros(model.getNonlinearDim(), 

231 dtype=lsst.meas.modelfit.Scalar)) 

232 self.assertFloatsAlmostEqual(fixed, ellipseParameters.ravel()) 

233 ellipses2 = model.writeEllipses(nonlinear, fixed) 

234 msf = model.makeShapeletFunction(nonlinear, amplitudes, fixed) 

235 self.assertFloatsAlmostEqual(len(msf.getComponents()), len(ellipses1)) 

236 ellipses3 = model.makeEllipseVector() 

237 amplitudeOffset = 0 

238 for i in range(len(ellipses2)): 

239 self.assertFloatsAlmostEqual(ellipses1[i].getParameterVector(), ellipses2[i].getParameterVector(), 

240 rtol=1E-8) 

241 # need to convert ellipse parametrization 

242 ellipses3[i].setCore(msf.getComponents()[i].getEllipse().getCore()) 

243 ellipses3[i].setCenter(msf.getComponents()[i].getEllipse().getCenter()) 

244 amplitudeCount = len(msf.getComponents()[i].getCoefficients()) 

245 self.assertFloatsAlmostEqual(ellipses1[i].getParameterVector(), ellipses3[i].getParameterVector(), 

246 rtol=1E-8) 

247 self.assertFloatsAlmostEqual(amplitudes[amplitudeOffset:amplitudeOffset+amplitudeCount], 

248 msf.getComponents()[i].getCoefficients(), rtol=1E-8) 

249 amplitudeOffset += amplitudeCount 

250 

251 # test the ellipse round-tripping again, this time starting with nonzero nonlinear parameters: 

252 # this will be read back in by adding to the fixed parameters and zeroing the nonlinear parameters. 

253 nonlinear[:] = 0.5*ellipseParameters.ravel() 

254 ellipses4 = model.writeEllipses(nonlinear, fixed) 

255 model.readEllipses(ellipses4, nonlinear, fixed) 

256 self.assertFloatsAlmostEqual(nonlinear, numpy.zeros(model.getNonlinearDim(), 

257 dtype=lsst.meas.modelfit.Scalar)) 

258 self.assertFloatsAlmostEqual(fixed, 1.5*ellipseParameters.ravel()) 

259 

260 def testApply(self): 

261 tolerances = {"full": 3E-4, "ellipse": 8E-3, "fixed": 1E-2} 

262 for filename in glob.glob(os.path.join(DATA_DIR, "psfs", "great3*.fits")): 

263 kernelImage = lsst.afw.image.ImageD(filename) 

264 shape = computeMoments(kernelImage) 

265 for configKey in ["full", "ellipse", "fixed"]: 

266 fitter = lsst.meas.modelfit.GeneralPsfFitter(self.configs[configKey].makeControl()) 

267 multiShapeletFit = fitter.apply(kernelImage, shape, 0.01) 

268 modelImage = lsst.afw.image.ImageD(kernelImage.getBBox(lsst.afw.image.PARENT)) 

269 multiShapeletFit.evaluate().addToImage(modelImage) 

270 self.assertFloatsAlmostEqual(kernelImage.getArray(), modelImage.getArray(), 

271 atol=tolerances[configKey], 

272 plotOnFailure=True) 

273 

274 

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

276 pass 

277 

278 

279def setup_module(module): 

280 lsst.utils.tests.init() 

281 

282 

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

284 lsst.utils.tests.init() 

285 unittest.main()