Coverage for python/lsst/cp/pipe/ptc/cpPtcAdjustGainRatios.py: 20%
140 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-20 02:37 -0800
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 logging
25import numpy as np
26from astropy.table import Table
28import lsst.pipe.base
29import lsst.pex.config
30import lsst.afw.cameraGeom
32from ..utils import FlatGradientFitter, FlatGainRatioFitter, bin_flat
35__all__ = [
36 "PhotonTransferCurveAdjustGainRatiosConfig",
37 "PhotonTransferCurveAdjustGainRatiosTask",
38]
41class PhotonTransferCurveAdjustGainRatiosConnections(
42 lsst.pipe.base.PipelineTaskConnections,
43 dimensions=("instrument", "detector"),
44):
45 exposures = lsst.pipe.base.connectionTypes.Input(
46 name="cpPtcIsrExp",
47 doc="Input exposures (from PTC ISR) for gain ratio adjustment.",
48 storageClass="Exposure",
49 dimensions=("instrument", "exposure", "detector"),
50 multiple=True,
51 deferLoad=True,
52 )
53 input_ptc = lsst.pipe.base.connectionTypes.Input(
54 name="ptcUnadjusted",
55 doc="Input PTC to have gain ratios adjusted.",
56 storageClass="PhotonTransferCurveDataset",
57 dimensions=("instrument", "detector"),
58 isCalibration=True,
59 )
60 output_ptc = lsst.pipe.base.connectionTypes.Output(
61 name="ptc",
62 doc="Output PTC after gain ratio adjustment.",
63 storageClass="PhotonTransferCurveDataset",
64 dimensions=("instrument", "detector"),
65 isCalibration=True,
66 )
67 output_adjust_summary = lsst.pipe.base.connectionTypes.Output(
68 name="ptc_adjustment_summary",
69 doc="Summary of gain adjustments.",
70 storageClass="ArrowAstropy",
71 dimensions=("instrument", "detector"),
72 )
75class PhotonTransferCurveAdjustGainRatiosConfig(
76 lsst.pipe.base.PipelineTaskConfig,
77 pipelineConnections=PhotonTransferCurveAdjustGainRatiosConnections,
78):
79 do_remove_radial_gradient = lsst.pex.config.Field(
80 dtype=bool,
81 doc="Remove radial gradient before fitting amp gain ratios?",
82 default=True,
83 )
84 radial_gradient_n_spline_nodes = lsst.pex.config.Field(
85 dtype=int,
86 doc="Number of radial spline nodes for radial gradient.",
87 default=20,
88 )
89 chebyshev_gradient_order = lsst.pex.config.Field(
90 dtype=int,
91 doc="Order of chebyshev x/y polynomials to remove additional gradients.",
92 default=1,
93 )
94 min_adu = lsst.pex.config.Field(
95 dtype=float,
96 doc="Minimum number of adu for an exposure to use in gain ratio calculation.",
97 default=1000.0,
98 )
99 max_adu = lsst.pex.config.Field(
100 dtype=float,
101 doc="Maximum number of adu for an exposure to use in gain ratio calculation.",
102 default=20000.0,
103 )
104 max_noise_reference = lsst.pex.config.Field(
105 dtype=float,
106 doc="Maximum read noise (e-) in the PTC for an amp to be considered as a reference",
107 default=12.0,
108 )
109 turnoff_percentile_reference = lsst.pex.config.Field(
110 dtype=float,
111 doc="Percentile threshold for sorting PTC turnoff for an amp to be considered as a reference",
112 default=25.0,
113 )
114 n_flat = lsst.pex.config.Field(
115 dtype=int,
116 doc="Number of flats (from min_adu to max_adu) to use in gain ratio calculation.",
117 default=50,
118 )
119 random_seed = lsst.pex.config.Field(
120 dtype=int,
121 doc="Random seed to use for down-sampling input flats.",
122 default=12345,
123 )
124 bin_factor = lsst.pex.config.Field(
125 dtype=int,
126 doc="Binning factor to compute gradients/gain ratios (pixels).",
127 default=8,
128 )
129 amp_boundary = lsst.pex.config.Field(
130 dtype=int,
131 doc="Amplifier boundary to ignore when computing gradients/gain ratios (pixels).",
132 default=20,
133 )
134 max_fractional_gain_ratio = lsst.pex.config.Field(
135 dtype=float,
136 doc="Maximum fractional gain ratio to consider per-exposure. Any amps with larger "
137 "offset will be excluded from the gradient fit and will have no corrections "
138 "applied.",
139 default=0.05,
140 )
143class PhotonTransferCurveAdjustGainRatiosTask(lsst.pipe.base.PipelineTask):
144 """Task to remove gradients to fit amp ratio gain adjustments.
145 """
146 ConfigClass = PhotonTransferCurveAdjustGainRatiosConfig
147 _DefaultName = "cpPhotonTransferCurveAdjustGainRatios"
149 def run(self, *, exposures, input_ptc):
150 """Run the gain adjustment task.
152 Parameters
153 ----------
154 exposures : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
155 Handles for input exposures.
156 input_ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
157 Input PTC to adjust.
159 Returns
160 -------
161 results : `lsst.pipe.base.Struct`
162 The output struct contains:
164 ``output_ptc``
165 The output modified PTC.
166 ``output_adjust_summary``
167 The summary of adjustments.
168 """
169 # Choose the exposures that will be used.
170 # For simplicity, we always use the first of a PTC pair.
171 rng = np.random.RandomState(seed=self.config.random_seed)
173 exp_ids_first = np.asarray(input_ptc.inputExpIdPairs[input_ptc.ampNames[0]])[:, 0]
175 avg = np.zeros(len(exp_ids_first))
176 n_amp = 0
177 for amp_name in input_ptc.ampNames:
178 if amp_name in input_ptc.badAmps:
179 continue
181 avg[:] += input_ptc.finalMeans[amp_name]
182 n_amp += 1
184 avg /= n_amp
186 exp_ids = exp_ids_first[((avg > self.config.min_adu) & (avg < self.config.max_adu))]
187 if len(exp_ids) > self.config.n_flat:
188 exp_ids = rng.choice(exp_ids, size=self.config.n_flat, replace=False)
189 exp_ids = np.sort(exp_ids)
191 fixed_amp_index = _choose_reference_amplifier(
192 input_ptc,
193 self.config.max_noise_reference,
194 self.config.turnoff_percentile_reference,
195 )
197 if fixed_amp_index < 0:
198 return lsst.pipe.base.Struct(output_ptc=input_ptc, output_adjust_summary=Table())
200 self.log.info(
201 "Using amplifier %d (%s) as fixed reference amplifier.",
202 fixed_amp_index,
203 input_ptc.ampNames[fixed_amp_index],
204 )
206 gain_ratio_array = np.zeros((len(exp_ids), len(input_ptc.ampNames)))
208 for i, exp_id in enumerate(exp_ids):
209 for handle in exposures:
210 if exp_id == handle.dataId["exposure"]:
211 exposure = handle.get()
212 break
214 self.log.info("Fitting gain ratios on exposure %d.", exp_id)
215 if not exposure.metadata.get("LSST ISR BOOTSTRAP", True):
216 raise RuntimeError(
217 "PhotonTransferCurveAdjustGainRatiosTask can only be run on bootstrap exposures.",
218 )
220 binned = bin_flat(
221 input_ptc,
222 exposure,
223 bin_factor=self.config.bin_factor,
224 amp_boundary=self.config.amp_boundary,
225 apply_gains=True,
226 )
228 # Clip out non-finite and extreme values.
229 # We need to be careful about unmasked defects, so we take
230 # a tighter percentile and expand.
231 lo, hi = np.nanpercentile(binned["value"], [5.0, 95.0])
232 lo *= 0.8
233 hi *= 1.2
234 use = (np.isfinite(binned["value"]) & (binned["value"] >= lo) & (binned["value"] <= hi))
235 binned = binned[use]
237 gain_ratio_array[i, :] = _compute_gain_ratios(
238 exposure.getDetector(),
239 binned,
240 fixed_amp_index,
241 do_remove_radial_gradient=self.config.do_remove_radial_gradient,
242 radial_gradient_n_spline_nodes=self.config.radial_gradient_n_spline_nodes,
243 chebyshev_gradient_order=self.config.chebyshev_gradient_order,
244 max_fractional_gain_ratio=self.config.max_fractional_gain_ratio,
245 log=self.log,
246 )
248 output_ptc = copy.copy(input_ptc)
250 # Compute the summary table.
251 summary = Table()
252 summary.meta["fixed_amp_index"] = fixed_amp_index
253 summary.meta["fixed_amp_name"] = input_ptc.ampNames[fixed_amp_index]
255 summary["exp_id"] = exp_ids
257 summary["mean_adu"] = np.zeros(len(exp_ids))
259 a, b = esutil.numpy_util.match(exp_ids_first, exp_ids)
260 summary["mean_adu"][b] = avg[a]
262 corrections = np.zeros(len(input_ptc.ampNames))
263 for i, amp_name in enumerate(input_ptc.ampNames):
264 summary[f"{amp_name}_gain_ratio"] = gain_ratio_array[:, i]
266 corrections[i] = np.median(gain_ratio_array[:, i])
268 # Adjust the median correction to 1.0 so we do not change the
269 # gain of the detector on average.
270 # This is needed in case the reference amplifier is skewed in
271 # terms of offsets even though it has the median gain.
272 med_correction = np.median(corrections)
274 summary.meta["median_correction"] = med_correction
276 for i, amp_name in enumerate(output_ptc.ampNames):
277 correction = corrections[i] / med_correction
278 new_gain = output_ptc.gainUnadjusted[amp_name] / correction
279 self.log.info(
280 "Adjusting gain from amplifier %s by factor of %.5f (from %.5f to %.5f)",
281 amp_name,
282 correction,
283 output_ptc.gain[amp_name],
284 new_gain,
285 )
286 output_ptc.gain[amp_name] = new_gain
288 return lsst.pipe.base.Struct(output_ptc=output_ptc, output_adjust_summary=summary)
291def _choose_reference_amplifier(
292 ptc,
293 max_noise_reference,
294 turnoff_percentile_reference,
295 log=None,
296):
297 """Choose a reference amplifier from a PTC.
299 Parameters
300 ----------
301 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
302 Input PTC.
303 max_noise_reference : `float`
304 Maximum read noise (e-) in the PTC for an amp to be considered.
305 turnoff_percentile_reference : `float`
306 Percentile threshold for sorting PTC turnoff for an amp to be
307 considered.
308 log : `logging.logger`, optional
309 Log object.
311 Returns
312 -------
313 reference_amp_index : `int`
314 Index of the reference amplifier.
315 """
316 log = log if log else logging.getLogger(__name__)
318 gain_array = np.zeros(len(ptc.ampNames))
319 turnoff_array = np.zeros(len(ptc.ampNames))
320 for i, amp_name in enumerate(ptc.ampNames):
321 gain_array[i] = ptc.gain[amp_name]
322 turnoff_array[i] = ptc.ptcTurnoff[amp_name]
323 # Do not consider any amplifier with high noise
324 # as a reference amplifier.
325 if ptc.noise[amp_name] > max_noise_reference:
326 log.info(
327 "Excluding amplifier %s as a possible reference amplifier due to high noise (%.2f)",
328 amp_name,
329 ptc.noise[amp_name],
330 )
331 gain_array[i] = np.nan
333 turnoff_threshold = np.nanpercentile(turnoff_array, turnoff_percentile_reference)
334 gain_array[turnoff_array < turnoff_threshold] = np.nan
335 good, = np.where(np.isfinite(gain_array))
337 if len(good) <= 1:
338 log.warning("Insufficient good amplifiers for PTC gain adjustment.")
339 return -1
341 st = np.argsort(np.nan_to_num(gain_array[good]))
342 ref_amp_index = good[st[int(0.5*len(good))]]
344 return ref_amp_index
347def _compute_gain_ratios(
348 detector,
349 binned,
350 fixed_amp_index,
351 do_remove_radial_gradient=True,
352 radial_gradient_n_spline_nodes=20,
353 chebyshev_gradient_order=1,
354 max_fractional_gain_ratio=0.05,
355 nsig_clip=5.0,
356 log=None,
357):
358 """Compute the gain ratios from a given non-gain-corrected exposure.
360 Parameters
361 ----------
362 detector : `lsst.afw.cameraGeom.Detector`
363 Detector object.
364 binned : `astropy.table.Table`
365 Table of binned values. Will contain ``xd``, ``yd`` (detector
366 bin positions); ``value`` (binned value); ``amp_index``
367 (index of the amplifiers).
368 fixed_amp_index : `int`
369 Use this amp index as the fixed point (gain ratio == 1.0).
370 do_remove_radial_gradient : `bool`, optional
371 Remove radial gradient before fitting amp gain ratios?
372 radial_gradient_n_spline_nodes : `int`, optional
373 Number of radial spline nodes for radial gradient.
374 chebyshev_gradient_order : `int`, optional
375 Order of chebyshev x/y polynomials to remove additional gradients.
376 max_fractional_gain_ratio : `int`, optional
377 Maximum fractional gain ratio to consider per-exposure. Any
378 amps with larger offset will be excluded from the gradient fit
379 and will have no corrections computed.
380 nsig_clip : `float`, optional
381 Number of sigma on distribution of gain ratios to clip when fitting
382 out Chebyshev gradient.
383 log : `logging.logger`, optional
384 Log object.
386 Returns
387 -------
388 amp_gain_ratios : `np.ndarray`
389 Amp gain ratios, relative to fixed amp.
390 """
391 log = log if log else logging.getLogger(__name__)
393 n_amp = len(detector)
394 amp_names = [amp.getName() for amp in detector]
396 # We iterate up to 8x to look for any bad amps (that have a gain offset
397 # fit greater than max_fractional_gain_ratio) and reject them
398 # from the fits.
399 bad_amps_converged = False
400 n_iter = 0
401 value_uncorrected = binned["value"].copy()
402 while not bad_amps_converged and n_iter < n_amp // 2:
403 # If configured, fit the radial gradient.
404 if do_remove_radial_gradient:
405 transform = detector.getTransform(
406 lsst.afw.cameraGeom.PIXELS,
407 lsst.afw.cameraGeom.FOCAL_PLANE,
408 )
409 xy = np.vstack((binned["xd"], binned["yd"]))
410 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)
411 xf = xf.ravel()
412 yf = yf.ravel()
413 radius = np.sqrt(xf**2. + yf**2.)
415 nodes = np.linspace(
416 np.min(radius),
417 np.max(radius),
418 radial_gradient_n_spline_nodes,
419 )
421 # Put in a normalization for fitting.
422 norm = np.nanpercentile(binned["value"], 95.0)
424 fitter = FlatGradientFitter(
425 nodes,
426 xf,
427 yf,
428 binned["value"]/norm,
429 np.array([]),
430 constrain_zero=False,
431 fit_centroid=False,
432 fit_gradient=False,
433 fp_centroid_x=0.0,
434 fp_centroid_y=0.0,
435 )
436 p0 = fitter.compute_p0()
437 pars = fitter.fit(p0)
439 binned["value"] /= fitter.compute_model(pars)
441 # Fit gain ratios.
442 fitter = FlatGainRatioFitter(
443 detector.getBBox(),
444 chebyshev_gradient_order,
445 binned["xd"],
446 binned["yd"],
447 binned["amp_index"],
448 binned["value"],
449 np.arange(n_amp),
450 fixed_amp_index,
451 )
452 pars = fitter.fit(nsig_clip=nsig_clip)
453 amp_pars = pars[fitter.indices["amp_pars"]]
455 fractional_gain_ratio = np.abs(amp_pars - 1.0)
456 max_ratio_ind = np.argmax(fractional_gain_ratio)
458 if fractional_gain_ratio[max_ratio_ind] > max_fractional_gain_ratio:
459 log.warning(
460 "Found bad amp %s with offset parameter %.2f",
461 amp_names[max_ratio_ind],
462 amp_pars[max_ratio_ind],
463 )
464 good = (binned["amp_index"] != max_ratio_ind)
465 binned = binned[good]
466 binned["value"] = value_uncorrected[good]
468 value_uncorrected = binned["value"].copy()
469 else:
470 bad_amps_converged = True
472 n_iter += 1
474 return amp_pars