Coverage for python/lsst/cp/pipe/ptc/cpPtcAdjustGainRatios.py: 20%
122 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-09 04:42 -0700
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 <https://www.gnu.org/licenses/>.
21#
22import copy
23import esutil
24import numpy as np
25from astropy.table import Table
27import lsst.pipe.base
28import lsst.pex.config
29import lsst.afw.cameraGeom
31from ..utils import FlatGradientFitter, FlatGainRatioFitter, bin_flat
34__all__ = [
35 "PhotonTransferCurveAdjustGainRatiosConfig",
36 "PhotonTransferCurveAdjustGainRatiosTask",
37]
40class PhotonTransferCurveAdjustGainRatiosConnections(
41 lsst.pipe.base.PipelineTaskConnections,
42 dimensions=("instrument", "detector"),
43):
44 exposures = lsst.pipe.base.connectionTypes.Input(
45 name="cpPtcIsrExp",
46 doc="Input exposures (from PTC ISR) for gain ratio adjustment.",
47 storageClass="Exposure",
48 dimensions=("instrument", "exposure", "detector"),
49 multiple=True,
50 deferLoad=True,
51 )
52 input_ptc = lsst.pipe.base.connectionTypes.Input(
53 name="ptcUnadjusted",
54 doc="Input PTC to have gain ratios adjusted.",
55 storageClass="PhotonTransferCurveDataset",
56 dimensions=("instrument", "detector"),
57 isCalibration=True,
58 )
59 output_ptc = lsst.pipe.base.connectionTypes.Output(
60 name="ptc",
61 doc="Output PTC after gain ratio adjustment.",
62 storageClass="PhotonTransferCurveDataset",
63 dimensions=("instrument", "detector"),
64 isCalibration=True,
65 )
66 output_adjust_summary = lsst.pipe.base.connectionTypes.Output(
67 name="ptc_adjustment_summary",
68 doc="Summary of gain adjustments.",
69 storageClass="ArrowAstropy",
70 dimensions=("instrument", "detector"),
71 )
74class PhotonTransferCurveAdjustGainRatiosConfig(
75 lsst.pipe.base.PipelineTaskConfig,
76 pipelineConnections=PhotonTransferCurveAdjustGainRatiosConnections,
77):
78 do_remove_radial_gradient = lsst.pex.config.Field(
79 dtype=bool,
80 doc="Remove radial gradient before fitting amp gain ratios?",
81 default=True,
82 )
83 radial_gradient_n_spline_nodes = lsst.pex.config.Field(
84 dtype=int,
85 doc="Number of radial spline nodes for radial gradient.",
86 default=20,
87 )
88 chebyshev_gradient_order = lsst.pex.config.Field(
89 dtype=int,
90 doc="Order of chebyshev x/y polynomials to remove additional gradients.",
91 default=1,
92 )
93 min_adu = lsst.pex.config.Field(
94 dtype=float,
95 doc="Minimum number of adu for an exposure to use in gain ratio calculation.",
96 default=1000.0,
97 )
98 max_adu = lsst.pex.config.Field(
99 dtype=float,
100 doc="Maximum number of adu for an exposure to use in gain ratio calculation.",
101 default=20000.0,
102 )
103 n_flat = lsst.pex.config.Field(
104 dtype=int,
105 doc="Number of flats (from min_adu to max_adu) to use in gain ratio calculation.",
106 default=50,
107 )
108 random_seed = lsst.pex.config.Field(
109 dtype=int,
110 doc="Random seed to use for down-sampling input flats.",
111 default=12345,
112 )
113 bin_factor = lsst.pex.config.Field(
114 dtype=int,
115 doc="Binning factor to compute gradients/gain ratios (pixels).",
116 default=8,
117 )
118 amp_boundary = lsst.pex.config.Field(
119 dtype=int,
120 doc="Amplifier boundary to ignore when computing gradients/gain ratios (pixels).",
121 default=20,
122 )
123 max_fractional_gain_ratio = lsst.pex.config.Field(
124 dtype=float,
125 doc="Maximum fractional gain ratio to consider per-exposure. Any amps with larger "
126 "offset will be excluded from the gradient fit and will have no corrections "
127 "applied.",
128 default=0.05,
129 )
132class PhotonTransferCurveAdjustGainRatiosTask(lsst.pipe.base.PipelineTask):
133 """Task to remove gradients to fit amp ratio gain adjustments.
134 """
135 ConfigClass = PhotonTransferCurveAdjustGainRatiosConfig
136 _DefaultName = "cpPhotonTransferCurveAdjustGainRatios"
138 def run(self, *, exposures, input_ptc):
139 """Run the gain adjustment task.
141 Parameters
142 ----------
143 exposures : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
144 Handles for input exposures.
145 input_ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
146 Input PTC to adjust.
148 Returns
149 -------
150 results : `lsst.pipe.base.Struct`
151 The output struct contains:
153 ``output_ptc``
154 The output modified PTC.
155 """
156 # Choose the exposures that will be used.
157 # For simplicity, we always use the first of a PTC pair.
158 rng = np.random.RandomState(seed=self.config.random_seed)
160 exp_ids_first = np.asarray(input_ptc.inputExpIdPairs[input_ptc.ampNames[0]])[:, 0]
162 avg = np.zeros(len(exp_ids_first))
163 n_amp = 0
164 for amp_name in input_ptc.ampNames:
165 if amp_name in input_ptc.badAmps:
166 continue
168 avg[:] += input_ptc.finalMeans[amp_name]
169 n_amp += 1
171 avg /= n_amp
173 exp_ids = exp_ids_first[((avg > self.config.min_adu) & (avg < self.config.max_adu))]
174 if len(exp_ids) > self.config.n_flat:
175 exp_ids = rng.choice(exp_ids, size=self.config.n_flat, replace=False)
176 exp_ids = np.sort(exp_ids)
178 # Figure out the reference amplifier.
179 gain_array = np.zeros(len(input_ptc.ampNames))
180 for i, amp_name in enumerate(input_ptc.ampNames):
181 gain_array[i] = input_ptc.gain[amp_name]
182 good, = np.where(np.isfinite(gain_array))
184 if len(good) <= 1:
185 self.log.warning("Insufficient good amplifiers for PTC gain adjustment.")
186 return lsst.pipe.base.Struct(output_ptc=input_ptc)
188 st = np.argsort(gain_array[good])
189 fixed_amp_index = good[st[int(0.5*len(good))]]
190 self.log.info(
191 "Using amplifier %d (%s) as fixed reference amplifier.",
192 fixed_amp_index,
193 input_ptc.ampNames[fixed_amp_index],
194 )
196 gain_ratio_array = np.zeros((len(exp_ids), len(input_ptc.ampNames)))
198 for i, exp_id in enumerate(exp_ids):
199 for handle in exposures:
200 if exp_id == handle.dataId["exposure"]:
201 exposure = handle.get()
202 break
204 self.log.info("Fitting gain ratios on exposure %d.", exp_id)
205 gain_ratio_array[i, :] = self._compute_gain_ratios(input_ptc, exposure, fixed_amp_index)
207 output_ptc = copy.copy(input_ptc)
209 # Compute the summary table.
210 summary = Table()
211 summary.meta["fixed_amp_index"] = fixed_amp_index
212 summary.meta["fixed_amp_name"] = input_ptc.ampNames[fixed_amp_index]
214 summary["exp_id"] = exp_ids
216 summary["mean_adu"] = np.zeros(len(exp_ids))
218 a, b = esutil.numpy_util.match(exp_ids_first, exp_ids)
219 summary["mean_adu"][b] = avg[a]
221 corrections = np.zeros(len(input_ptc.ampNames))
222 for i, amp_name in enumerate(input_ptc.ampNames):
223 summary[f"{amp_name}_gain_ratio"] = gain_ratio_array[:, i]
225 corrections[i] = np.median(gain_ratio_array[:, i])
227 # Adjust the median correction to 1.0 so we do not change the
228 # gain of the detector on average.
229 # This is needed in case the reference amplifier is skewed in
230 # terms of offsets even though it has the median gain.
231 med_correction = np.median(corrections)
233 summary.meta["median_correction"] = med_correction
235 for i, amp_name in enumerate(output_ptc.ampNames):
236 correction = corrections[i] / med_correction
237 new_gain = output_ptc.gainUnadjusted[amp_name] / correction
238 self.log.info(
239 "Adjusting gain from amplifier %s by factor of %.5f (from %.5f to %.5f)",
240 amp_name,
241 correction,
242 output_ptc.gain[amp_name],
243 new_gain,
244 )
245 output_ptc.gain[amp_name] = new_gain
247 return lsst.pipe.base.Struct(output_ptc=output_ptc, gain_adjust_summary=summary)
249 def _compute_gain_ratios(self, ptc, exposure, fixed_amp_index):
250 """Compute the gain ratios from a given non-gain-corrected exposure.
252 Parameters
253 ----------
254 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
255 PTC to correct exposure.
256 exposure : `lsst.afw.image.Exposure`
257 Exposure to measure gain ratios.
258 fixed_amp_index : `int`
259 Use this amp index as the fixed point (gain ratio == 1.0).
261 Returns
262 -------
263 amp_gain_ratios : `np.ndarray`
264 Amp gain ratios, relative to fixed amp.
265 """
266 if not exposure.metadata.get("LSST ISR BOOTSTRAP", True):
267 raise RuntimeError(
268 "PhotonTransferCurveAdjustGainRatiosTask can only be run on bootstrap exposures.",
269 )
271 detector = exposure.getDetector()
273 binned = bin_flat(
274 ptc,
275 exposure,
276 bin_factor=self.config.bin_factor,
277 amp_boundary=self.config.amp_boundary,
278 apply_gains=True,
279 )
281 # Clip out non-finite and extreme values.
282 # We need to be careful about unmasked defects, so we take
283 # a tighter percentile and expand.
284 lo, hi = np.nanpercentile(binned["value"], [5.0, 95.0])
285 lo *= 0.8
286 hi *= 1.2
287 use = (np.isfinite(binned["value"]) & (binned["value"] >= lo) & (binned["value"] <= hi))
288 binned = binned[use]
290 # We iterate up to 8x to look for any bad amps (that have a gain offset
291 # fit greater than config.max_fractional_gain_ratio) and reject them
292 # from the fits.
293 bad_amps_converged = False
294 n_iter = 0
295 value_uncorrected = binned["value"].copy()
296 while not bad_amps_converged and n_iter < len(ptc.ampNames) // 2:
297 # If configured, fit the radial gradient.
298 if self.config.do_remove_radial_gradient:
299 transform = detector.getTransform(
300 lsst.afw.cameraGeom.PIXELS,
301 lsst.afw.cameraGeom.FOCAL_PLANE,
302 )
303 xy = np.vstack((binned["xd"], binned["yd"]))
304 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)
305 xf = xf.ravel()
306 yf = yf.ravel()
307 radius = np.sqrt(xf**2. + yf**2.)
309 nodes = np.linspace(
310 np.min(radius),
311 np.max(radius),
312 self.config.radial_gradient_n_spline_nodes,
313 )
315 # Put in a normalization for fitting.
316 norm = np.nanpercentile(binned["value"], 95.0)
318 fitter = FlatGradientFitter(
319 nodes,
320 xf,
321 yf,
322 binned["value"]/norm,
323 np.array([]),
324 constrain_zero=False,
325 fit_centroid=False,
326 fit_gradient=False,
327 fp_centroid_x=0.0,
328 fp_centroid_y=0.0,
329 )
330 p0 = fitter.compute_p0()
331 pars = fitter.fit(p0)
333 binned["value"] /= fitter.compute_model(pars)
335 # Fit gain ratios.
336 fitter = FlatGainRatioFitter(
337 exposure.getBBox(),
338 self.config.chebyshev_gradient_order,
339 binned["xd"],
340 binned["yd"],
341 binned["amp_index"],
342 binned["value"],
343 np.arange(len(ptc.ampNames)),
344 fixed_amp_index,
345 )
346 pars = fitter.fit()
347 amp_pars = pars[fitter.indices["amp_pars"]]
349 fractional_gain_ratio = np.abs(amp_pars - 1.0)
350 max_ratio_ind = np.argmax(fractional_gain_ratio)
352 if fractional_gain_ratio[max_ratio_ind] > self.config.max_fractional_gain_ratio:
353 self.log.warning(
354 "Found bad amp %s with offset parameter %.2f",
355 ptc.ampNames[max_ratio_ind],
356 amp_pars[max_ratio_ind],
357 )
358 good = (binned["amp_index"] != max_ratio_ind)
359 binned = binned[good]
360 binned["value"] = value_uncorrected[good]
362 value_uncorrected = binned["value"].copy()
363 else:
364 bad_amps_converged = True
366 n_iter += 1
368 return amp_pars