Coverage for tests / test_fit_turbulence.py: 19%

1# This file is part of drp_tasks. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

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 GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21import unittest 

22 

23import astropy.units as u 

24import numpy as np 

25import treegp 

26from astropy.table import Table, vstack 

27 

28import lsst.afw.table as afwTable 

29import lsst.utils.tests 

30from lsst.drp.tasks.fit_turbulence import ( 

31 GaussianProcessesTurbulenceFitConfig, 

32 GaussianProcessesTurbulenceFitTask, 

33) 

34from lsst.obs.base import createInitialSkyWcsFromBoresight 

35from lsst.obs.base.instrument_tests import DummyCam 

36 

37 

38class FitTurbulenceTestCase(lsst.utils.tests.TestCase): 

39 

40 def setUp(self): 

41 

42 visit = 1234 

43 self.noise = 0.01 

44 self.npoints = 250 

45 self.rng = np.random.RandomState(12345) 

46 

47 config = GaussianProcessesTurbulenceFitConfig() 

48 # Set the max separation lower, since we only have two detectors. 

49 config.correlationSeparationMax = 0.1 

50 config.initKernel = "15**2 * AnisotropicVonKarman(invLam=array([[1.0/0.1**2,0],[0,1.0/0.1**2]]))" 

51 config.splineNNodes = 60 

52 self.task = GaussianProcessesTurbulenceFitTask(config=config) 

53 

54 kernel_skl = treegp.eval_kernel(self.task.config.initKernel) 

55 

56 # Make a simulated WCS catalog and position catalog. 

57 self.wcsCatalog = self._makeWcsCatalog(visit) 

58 self.positions = self._makeResidualCatalog(self.wcsCatalog, visit, kernel_skl) 

59 

60 # Run the Gaussian Processes and predict at the input positions. 

61 self.gpx, self.gpy, *_ = self.task.runGP(self.wcsCatalog, self.positions) 

62 

63 sourceCatalog = Table( 

64 { 

65 "x": self.positions["xpix"], 

66 "y": self.positions["ypix"], 

67 "detector": self.positions["deviceName"], 

68 } 

69 ) 

70 self.predictCat = self.task.predict(self.gpx, self.gpy, self.wcsCatalog, sourceCatalog) 

71 

72 def _makeWcsCatalog(self, visit): 

73 """Make a simulated WCS catalog.""" 

74 camera = DummyCam().getCamera() 

75 schema = afwTable.ExposureTable.makeMinimalSchema() 

76 schema.addField("visit", type="L", doc="Visit number") 

77 wcsCatalog = afwTable.ExposureCatalog(schema) 

78 

79 boresight = lsst.geom.SpherePoint(150 * lsst.geom.degrees, 2.5 * lsst.geom.degrees) 

80 

81 for detectorId in [0, 1]: 

82 detector = camera[detectorId] 

83 wcs = createInitialSkyWcsFromBoresight(boresight, 0 * lsst.geom.degrees, detector) 

84 record = wcsCatalog.addNew() 

85 record.setId(detectorId) 

86 record.setWcs(wcs) 

87 record.setBBox(detector.getBBox()) 

88 wcsCatalog["visit"] = visit 

89 return wcsCatalog 

90 

91 def _makeResidualCatalog(self, wcsCatalog, visit, kernel): 

92 """Make a simulated catalog of source positions and residuals.""" 

93 visitPositions = [] 

94 for detector in [0, 1]: 

95 bbox = wcsCatalog[detector].getBBox() 

96 pixToTP = wcsCatalog[detector].wcs.getFrameDict().getMapping("PIXELS", "IWC") 

97 

98 x = self.rng.uniform(bbox.getBeginX(), bbox.getEndX(), self.npoints) 

99 y = self.rng.uniform(bbox.getBeginY(), bbox.getEndY(), self.npoints) 

100 xTP, yTP = pixToTP.applyForward(np.array([x, y])) 

101 xSky, ySky = wcsCatalog[detector].wcs.pixelToSkyArray(x, y, degrees=True) 

102 detectorPositions = Table( 

103 { 

104 "xpix": x, 

105 "ypix": y, 

106 "xworld": xTP, 

107 "yworld": yTP, 

108 "xsky": xSky, 

109 "ysky": ySky, 

110 } 

111 ) 

112 detectorPositions["deviceName"] = detector 

113 visitPositions.append(detectorPositions) 

114 visitPositions = vstack(visitPositions) 

115 totPoints = len(visitPositions) 

116 

117 K = kernel.__call__(np.array([visitPositions["xworld"], visitPositions["yworld"]]).T) 

118 dx = self.rng.multivariate_normal(np.zeros(totPoints), K) 

119 dy = self.rng.multivariate_normal(np.zeros(totPoints), K) 

120 

121 dx += np.random.normal(scale=self.noise, size=totPoints) 

122 dy += np.random.normal(scale=self.noise, size=totPoints) 

123 

124 residError = np.ones(totPoints) * self.noise 

125 visitPositions["xresw"] = dx 

126 visitPositions["yresw"] = dy 

127 visitPositions["covTotalW_00"] = residError**2 

128 visitPositions["covTotalW_11"] = residError**2 

129 visitPositions["exposureName"] = str(visit) 

130 

131 return visitPositions 

132 

133 def test_runGP(self): 

134 """Check that the predicted position is about the same as the simulated 

135 dx/dy values, modulo the scatter in the data.""" 

136 dRA = self.positions["xsky"] - self.predictCat["coord_ra"] 

137 dDec = self.positions["ysky"] - self.predictCat["coord_dec"] 

138 

139 np.testing.assert_allclose(dRA, (self.positions["xresw"] * u.mas).to(u.degree).value, rtol=1e-1) 

140 np.testing.assert_allclose(dDec, (self.positions["yresw"] * u.mas).to(u.degree).value, rtol=1e-1) 

141 

142 def test_makeWcs(self): 

143 """Check that the output WCS maps pixel positions to the same RA and 

144 Dec as the Gaussian Processes prediction. 

145 """ 

146 

147 wcsWithGP = self.task.addGPToWcs(self.gpx, self.gpy, self.wcsCatalog) 

148 for detectorRow in self.wcsCatalog: 

149 detector = detectorRow.getId() 

150 newWcs = wcsWithGP.find(detector).wcs 

151 

152 sources = self.predictCat[self.predictCat["detector"] == detector] 

153 ra, dec = newWcs.pixelToSkyArray(sources["x"], sources["y"], degrees=True) 

154 

155 np.testing.assert_allclose(ra, sources["coord_ra"]) 

156 np.testing.assert_allclose(dec, sources["coord_dec"], rtol=1e-6) 

157 

158 

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

160 lsst.utils.tests.init() 

161 unittest.main()