Coverage for python / lsst / cp / pipe / cpGainCorrection.py: 30%

1# This file is part of cp_pipe. 

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 <http://www.gnu.org/licenses/>. 

21import numpy as np 

22from scipy.stats import median_abs_deviation 

23 

24import lsst.pipe.base 

25from lsst.ip.isr import GainCorrection 

26import lsst.pex.config 

27from lsst.utils.plotting import make_figure 

28 

29from .ptc.cpPtcAdjustGainRatios import _choose_reference_amplifier, _compute_gain_ratios 

30from .utils import bin_flat 

31 

32__all__ = ["CpMeasureGainCorrectionTask", "CpMeasureGainCorrectionConfig"] 

33 

34 

35class CpMeasureGainCorrectionConnections( 

36 lsst.pipe.base.PipelineTaskConnections, 

37 dimensions=("instrument", "detector"), 

38): 

39 input_reference_flat = lsst.pipe.base.connectionTypes.PrerequisiteInput( 

40 name="flat", 

41 doc="Input reference flat (certified).", 

42 storageClass="Exposure", 

43 dimensions=("instrument", "detector", "physical_filter"), 

44 isCalibration=True, 

45 ) 

46 input_reference_ptc = lsst.pipe.base.connectionTypes.PrerequisiteInput( 

47 name="ptc", 

48 doc="Input reference ptc (certified).", 

49 storageClass="PhotonTransferCurveDataset", 

50 dimensions=("instrument", "detector"), 

51 isCalibration=True, 

52 ) 

53 input_flat = lsst.pipe.base.connectionTypes.Input( 

54 name="flat_gain_correction", 

55 doc="Input flat to derive gain correction.", 

56 storageClass="Exposure", 

57 dimensions=("instrument", "detector", "physical_filter"), 

58 isCalibration=True, 

59 ) 

60 output_gain_correction = lsst.pipe.base.connectionTypes.Output( 

61 name="gain_correction", 

62 doc="Output gain correction calibration.", 

63 storageClass="IsrCalib", 

64 dimensions=("instrument", "detector"), 

65 isCalibration=True, 

66 ) 

67 output_flat_ratio_plot = lsst.pipe.base.connectionTypes.Output( 

68 name="gain_correction_flat_ratio_plot", 

69 doc="Gain correction flat ratio plot.", 

70 storageClass="Plot", 

71 dimensions=("instrument", "detector"), 

72 ) 

73 

74 def adjustQuantum(self, inputs, outputs, label, dataId): 

75 # We need to remove reference flats that do not match 

76 # the physical filter of the input flat. 

77 input_flat_ref = inputs["input_flat"][1][0] 

78 physical_filter = input_flat_ref.dataId["physical_filter"] 

79 

80 input_reference_flat_refs = inputs["input_reference_flat"][1] 

81 input_reference_flat_ref = None 

82 for ref in input_reference_flat_refs: 

83 if ref.dataId["physical_filter"] == physical_filter: 

84 input_reference_flat_ref = ref 

85 break 

86 

87 if input_reference_flat_ref is None: 

88 raise lsst.pipe.base.NoWorkFound( 

89 f"No matched input flat with physical filter {physical_filter}.", 

90 ) 

91 

92 inputs["input_reference_flat"] = (inputs["input_reference_flat"][0], (input_reference_flat_ref,)) 

93 

94 return inputs, outputs 

95 

96 

97class CpMeasureGainCorrectionConfig( 

98 lsst.pipe.base.PipelineTaskConfig, 

99 pipelineConnections=CpMeasureGainCorrectionConnections, 

100): 

101 max_noise_reference = lsst.pex.config.Field( 

102 dtype=float, 

103 doc="Maximum read noise (e-) in the PTC for an amp to be considered as a reference", 

104 default=12.0, 

105 ) 

106 turnoff_percentile_reference = lsst.pex.config.Field( 

107 dtype=float, 

108 doc="Percentile threshold for sorting PTC turnoff for an amp to be considered as a reference", 

109 default=25.0, 

110 ) 

111 do_remove_radial_gradient = lsst.pex.config.Field( 

112 dtype=bool, 

113 doc="Remove radial gradient before measureting amp gain ratios?", 

114 default=True, 

115 ) 

116 radial_gradient_n_spline_nodes = lsst.pex.config.Field( 

117 dtype=int, 

118 doc="Number of radial spline nodes for radial gradient.", 

119 default=20, 

120 ) 

121 chebyshev_gradient_order = lsst.pex.config.Field( 

122 dtype=int, 

123 doc="Order of chebyshev x/y polynomials to remove additional gradients.", 

124 default=1, 

125 ) 

126 chebyshev_gradient_nsig_clip = lsst.pex.config.Field( 

127 dtype=float, 

128 doc="Number of sigma to clip gain ratios when fitting the Chebyshev gradient.", 

129 default=5.0, 

130 ) 

131 bin_factor = lsst.pex.config.Field( 

132 dtype=int, 

133 doc="Binning factor to compute gradients/gain ratios (pixels).", 

134 default=8, 

135 ) 

136 amp_boundary = lsst.pex.config.Field( 

137 dtype=int, 

138 doc="Amplifier boundary to ignore when computing gradients/gain ratios (pixels).", 

139 default=20, 

140 ) 

141 max_fractional_gain_ratio = lsst.pex.config.Field( 

142 dtype=float, 

143 doc="Maximum fractional gain ratio to consider. Any amps with larger " 

144 "offset will be excluded from the measurement and will have no corrections " 

145 "computed.", 

146 default=0.05, 

147 ) 

148 ratio_outlier_nsig_clip = lsst.pex.config.Field( 

149 dtype=float, 

150 doc="Number of sigma to clip in ratio image.", 

151 default=5.0, 

152 ) 

153 

154 

155class CpMeasureGainCorrectionTask(lsst.pipe.base.PipelineTask): 

156 """Task to measure gain corrections.""" 

157 

158 ConfigClass = CpMeasureGainCorrectionConfig 

159 _DefaultName = "cpMeasureGainCorrection" 

160 

161 def run(self, *, input_reference_flat, input_reference_ptc, input_flat): 

162 """ 

163 

164 Parameters 

165 ---------- 

166 input_reference_flat : `lsst.afw.image.Exposure` 

167 Input reference flat (typically a certified calibration). 

168 input_reference_ptc : `lsst.ip.isr.PhotonTransferCurveDataset` 

169 Input reference PTC (typically a certified calibration). 

170 input_flat : `lsst.afw.image.Exposure` 

171 Input flat to derive gain correction (relative to reference flat). 

172 

173 Returns 

174 ------- 

175 result : `lsst.pipe.base.Struct` 

176 Result struct containing: 

177 ``output_gain_correction`` : `lsst.ip.isr.GainCorrection` 

178 ``output_flat_ratio_plot`` : `matplotlib.Figure` 

179 """ 

180 gain_correction = GainCorrection( 

181 ampNames=input_reference_ptc.ampNames, 

182 gainAdjustments=np.ones(len(input_reference_ptc.ampNames)), 

183 ) 

184 

185 fixed_amp_index = _choose_reference_amplifier( 

186 input_reference_ptc, 

187 self.config.max_noise_reference, 

188 self.config.turnoff_percentile_reference, 

189 ) 

190 

191 if fixed_amp_index < 0: 

192 return lsst.pipe.base.Struct(output_gain_correction=gain_correction) 

193 

194 ratio = input_flat.clone() 

195 ratio.image /= input_reference_flat.image 

196 

197 binned_ref = bin_flat(input_reference_ptc, input_reference_flat, apply_gains=False) 

198 binned_flat = bin_flat(input_reference_ptc, input_flat, apply_gains=False) 

199 binned_ratio = bin_flat(input_reference_ptc, ratio, apply_gains=False) 

200 

201 # Clip out non-finite and extreme values from both the input flat 

202 # and the reference flat. 

203 ok = ( 

204 np.isfinite(binned_ref["value"]) 

205 & np.isfinite(binned_flat["value"]) 

206 & np.isfinite(binned_ratio["value"]) 

207 ) 

208 lo_ref, hi_ref = np.percentile(binned_ref["value"][ok], [5.0, 95.0]) 

209 lo_ref *= 0.8 

210 hi_ref *= 1.2 

211 lo_flat, hi_flat = np.percentile(binned_flat["value"][ok], [5.0, 95.0]) 

212 lo_flat *= 0.8 

213 hi_flat *= 1.2 

214 med_ratio = np.median(binned_ratio["value"][ok]) 

215 # Make sure this is non-zero for noise-free tests in particular. 

216 sig_ratio = np.clip( 

217 median_abs_deviation(binned_ratio["value"][ok], scale="normal"), 

218 0.01, 

219 None, 

220 ) 

221 lo_ratio = med_ratio - self.config.ratio_outlier_nsig_clip * sig_ratio 

222 hi_ratio = med_ratio + self.config.ratio_outlier_nsig_clip * sig_ratio 

223 use = ( 

224 np.isfinite(binned_ratio["value"]) 

225 & (binned_ref["value"] >= lo_ref) 

226 & (binned_ref["value"] <= hi_ref) 

227 & (binned_flat["value"] >= lo_flat) 

228 & (binned_flat["value"] <= hi_flat) 

229 & (binned_ratio["value"] >= lo_ratio) 

230 & (binned_ratio["value"] <= hi_ratio) 

231 ) 

232 binned = binned_ratio[use] 

233 

234 self.log.info( 

235 "Using amplifier %d (%s) as fixed reference amplifier in fit.", 

236 fixed_amp_index, 

237 input_reference_ptc.ampNames[fixed_amp_index], 

238 ) 

239 

240 gain_ratios = _compute_gain_ratios( 

241 input_flat.getDetector(), 

242 binned, 

243 fixed_amp_index, 

244 do_remove_radial_gradient=self.config.do_remove_radial_gradient, 

245 radial_gradient_n_spline_nodes=self.config.radial_gradient_n_spline_nodes, 

246 chebyshev_gradient_order=self.config.chebyshev_gradient_order, 

247 max_fractional_gain_ratio=self.config.max_fractional_gain_ratio, 

248 nsig_clip=self.config.chebyshev_gradient_nsig_clip, 

249 log=self.log, 

250 ) 

251 

252 # Sort ratios to find the approximate midpoint. The goal is 

253 # to have one amplifier exactly 1.0 correction, and if there 

254 # are an even number of amps simply taking the median will 

255 # interpolate between two ratios. In addition, this allows 

256 # us to explicitly log the new reference amplifier. 

257 st = np.argsort(gain_ratios) 

258 fixed_amp_index2 = st[int(0.5 * len(st))] 

259 

260 gain_ratios /= gain_ratios[fixed_amp_index2] 

261 

262 self.log.info( 

263 "Adjusting to amplifier %d (%s) as fixed reference amplifier after fit.", 

264 fixed_amp_index2, 

265 input_reference_ptc.ampNames[fixed_amp_index2], 

266 ) 

267 

268 gain_correction.gainAdjustments[:] = gain_ratios 

269 

270 value2 = binned["value"].copy() 

271 for ind in range(len(input_reference_ptc.ampNames)): 

272 value2[binned["amp_index"] == ind] /= gain_ratios[ind] 

273 

274 vmin, vmax = np.nanpercentile(binned["value"], [10, 90]) 

275 

276 fig = make_figure(figsize=(16, 6)) 

277 

278 ax1 = fig.add_subplot(121) 

279 

280 im1 = ax1.hexbin(binned["xd"], binned["yd"], C=binned["value"], vmin=vmin, vmax=vmax) 

281 

282 ax1.set_xlabel("Detector x (pix)") 

283 ax1.set_ylabel("Detector y (pix)") 

284 ax1.set_aspect("equal") 

285 ax1.set_title("Flat Ratio (uncorrected)") 

286 

287 fig.colorbar(im1, ax=ax1) 

288 

289 ax2 = fig.add_subplot(122) 

290 

291 im2 = ax2.hexbin(binned["xd"], binned["yd"], C=value2, vmin=vmin, vmax=vmax) 

292 

293 ax2.set_xlabel("Detector x (pix)") 

294 ax2.set_ylabel("Detector y (pix)") 

295 ax2.set_aspect("equal") 

296 ax2.set_title("Flat Ratio (corrected)") 

297 

298 fig.colorbar(im2, ax=ax2) 

299 

300 return lsst.pipe.base.Struct( 

301 output_gain_correction=gain_correction, 

302 output_flat_ratio_plot=fig, 

303 )