Coverage for python / lsst / cp / pipe / ptc / cpPtcAdjustGainRatios.py: 20%
141 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 09:34 +0000
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 = np.zeros(len(exp_ids_first), dtype=np.int64)
177 for amp_name in input_ptc.ampNames:
178 if amp_name in input_ptc.badAmps:
179 continue
181 finite = np.isfinite(input_ptc.finalMeans[amp_name])
182 avg[finite] += input_ptc.finalMeans[amp_name][finite]
183 n_amp[finite] += 1
185 avg /= n_amp
187 exp_ids = exp_ids_first[((avg > self.config.min_adu) & (avg < self.config.max_adu))]
188 if len(exp_ids) > self.config.n_flat:
189 exp_ids = rng.choice(exp_ids, size=self.config.n_flat, replace=False)
190 exp_ids = np.sort(exp_ids)
192 fixed_amp_index = _choose_reference_amplifier(
193 input_ptc,
194 self.config.max_noise_reference,
195 self.config.turnoff_percentile_reference,
196 )
198 if fixed_amp_index < 0:
199 return lsst.pipe.base.Struct(output_ptc=input_ptc, output_adjust_summary=Table())
201 self.log.info(
202 "Using amplifier %d (%s) as fixed reference amplifier.",
203 fixed_amp_index,
204 input_ptc.ampNames[fixed_amp_index],
205 )
207 gain_ratio_array = np.zeros((len(exp_ids), len(input_ptc.ampNames)))
209 for i, exp_id in enumerate(exp_ids):
210 for handle in exposures:
211 if exp_id == handle.dataId["exposure"]:
212 exposure = handle.get()
213 break
215 self.log.info("Fitting gain ratios on exposure %d.", exp_id)
216 if not exposure.metadata.get("LSST ISR BOOTSTRAP", True):
217 raise RuntimeError(
218 "PhotonTransferCurveAdjustGainRatiosTask can only be run on bootstrap exposures.",
219 )
221 binned = bin_flat(
222 input_ptc,
223 exposure,
224 bin_factor=self.config.bin_factor,
225 amp_boundary=self.config.amp_boundary,
226 apply_gains=True,
227 )
229 # Clip out non-finite and extreme values.
230 # We need to be careful about unmasked defects, so we take
231 # a tighter percentile and expand.
232 lo, hi = np.nanpercentile(binned["value"], [5.0, 95.0])
233 lo *= 0.8
234 hi *= 1.2
235 use = (np.isfinite(binned["value"]) & (binned["value"] >= lo) & (binned["value"] <= hi))
236 binned = binned[use]
238 gain_ratio_array[i, :] = _compute_gain_ratios(
239 exposure.getDetector(),
240 binned,
241 fixed_amp_index,
242 do_remove_radial_gradient=self.config.do_remove_radial_gradient,
243 radial_gradient_n_spline_nodes=self.config.radial_gradient_n_spline_nodes,
244 chebyshev_gradient_order=self.config.chebyshev_gradient_order,
245 max_fractional_gain_ratio=self.config.max_fractional_gain_ratio,
246 log=self.log,
247 )
249 output_ptc = copy.copy(input_ptc)
251 # Compute the summary table.
252 summary = Table()
253 summary.meta["fixed_amp_index"] = fixed_amp_index
254 summary.meta["fixed_amp_name"] = input_ptc.ampNames[fixed_amp_index]
256 summary["exp_id"] = exp_ids
258 summary["mean_adu"] = np.zeros(len(exp_ids))
260 a, b = esutil.numpy_util.match(exp_ids_first, exp_ids)
261 summary["mean_adu"][b] = avg[a]
263 corrections = np.zeros(len(input_ptc.ampNames))
264 for i, amp_name in enumerate(input_ptc.ampNames):
265 summary[f"{amp_name}_gain_ratio"] = gain_ratio_array[:, i]
267 corrections[i] = np.median(gain_ratio_array[:, i])
269 # Adjust the median correction to 1.0 so we do not change the
270 # gain of the detector on average.
271 # This is needed in case the reference amplifier is skewed in
272 # terms of offsets even though it has the median gain.
273 med_correction = np.median(corrections)
275 summary.meta["median_correction"] = med_correction
277 for i, amp_name in enumerate(output_ptc.ampNames):
278 correction = corrections[i] / med_correction
279 new_gain = output_ptc.gainUnadjusted[amp_name] / correction
280 self.log.info(
281 "Adjusting gain from amplifier %s by factor of %.5f (from %.5f to %.5f)",
282 amp_name,
283 correction,
284 output_ptc.gain[amp_name],
285 new_gain,
286 )
287 output_ptc.gain[amp_name] = new_gain
289 return lsst.pipe.base.Struct(output_ptc=output_ptc, output_adjust_summary=summary)
292def _choose_reference_amplifier(
293 ptc,
294 max_noise_reference,
295 turnoff_percentile_reference,
296 log=None,
297):
298 """Choose a reference amplifier from a PTC.
300 Parameters
301 ----------
302 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
303 Input PTC.
304 max_noise_reference : `float`
305 Maximum read noise (e-) in the PTC for an amp to be considered.
306 turnoff_percentile_reference : `float`
307 Percentile threshold for sorting PTC turnoff for an amp to be
308 considered.
309 log : `logging.logger`, optional
310 Log object.
312 Returns
313 -------
314 reference_amp_index : `int`
315 Index of the reference amplifier.
316 """
317 log = log if log else logging.getLogger(__name__)
319 gain_array = np.zeros(len(ptc.ampNames))
320 turnoff_array = np.zeros(len(ptc.ampNames))
321 for i, amp_name in enumerate(ptc.ampNames):
322 gain_array[i] = ptc.gain[amp_name]
323 turnoff_array[i] = ptc.ptcTurnoff[amp_name]
324 # Do not consider any amplifier with high noise
325 # as a reference amplifier.
326 if ptc.noise[amp_name] > max_noise_reference:
327 log.info(
328 "Excluding amplifier %s as a possible reference amplifier due to high noise (%.2f)",
329 amp_name,
330 ptc.noise[amp_name],
331 )
332 gain_array[i] = np.nan
334 turnoff_threshold = np.nanpercentile(turnoff_array, turnoff_percentile_reference)
335 gain_array[turnoff_array < turnoff_threshold] = np.nan
336 good, = np.where(np.isfinite(gain_array))
338 if len(good) <= 1:
339 log.warning("Insufficient good amplifiers for PTC gain adjustment.")
340 return -1
342 st = np.argsort(np.nan_to_num(gain_array[good]))
343 ref_amp_index = good[st[int(0.5*len(good))]]
345 return ref_amp_index
348def _compute_gain_ratios(
349 detector,
350 binned,
351 fixed_amp_index,
352 do_remove_radial_gradient=True,
353 radial_gradient_n_spline_nodes=20,
354 chebyshev_gradient_order=1,
355 max_fractional_gain_ratio=0.05,
356 nsig_clip=5.0,
357 log=None,
358):
359 """Compute the gain ratios from a given non-gain-corrected exposure.
361 Parameters
362 ----------
363 detector : `lsst.afw.cameraGeom.Detector`
364 Detector object.
365 binned : `astropy.table.Table`
366 Table of binned values. Will contain ``xd``, ``yd`` (detector
367 bin positions); ``value`` (binned value); ``amp_index``
368 (index of the amplifiers).
369 fixed_amp_index : `int`
370 Use this amp index as the fixed point (gain ratio == 1.0).
371 do_remove_radial_gradient : `bool`, optional
372 Remove radial gradient before fitting amp gain ratios?
373 radial_gradient_n_spline_nodes : `int`, optional
374 Number of radial spline nodes for radial gradient.
375 chebyshev_gradient_order : `int`, optional
376 Order of chebyshev x/y polynomials to remove additional gradients.
377 max_fractional_gain_ratio : `int`, optional
378 Maximum fractional gain ratio to consider per-exposure. Any
379 amps with larger offset will be excluded from the gradient fit
380 and will have no corrections computed.
381 nsig_clip : `float`, optional
382 Number of sigma on distribution of gain ratios to clip when fitting
383 out Chebyshev gradient.
384 log : `logging.logger`, optional
385 Log object.
387 Returns
388 -------
389 amp_gain_ratios : `np.ndarray`
390 Amp gain ratios, relative to fixed amp.
391 """
392 log = log if log else logging.getLogger(__name__)
394 n_amp = len(detector)
395 amp_names = [amp.getName() for amp in detector]
397 # We iterate up to 8x to look for any bad amps (that have a gain offset
398 # fit greater than max_fractional_gain_ratio) and reject them
399 # from the fits.
400 bad_amps_converged = False
401 n_iter = 0
402 value_uncorrected = binned["value"].copy()
403 while not bad_amps_converged and n_iter < n_amp // 2:
404 # If configured, fit the radial gradient.
405 if do_remove_radial_gradient:
406 transform = detector.getTransform(
407 lsst.afw.cameraGeom.PIXELS,
408 lsst.afw.cameraGeom.FOCAL_PLANE,
409 )
410 xy = np.vstack((binned["xd"], binned["yd"]))
411 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)
412 xf = xf.ravel()
413 yf = yf.ravel()
414 radius = np.sqrt(xf**2. + yf**2.)
416 nodes = np.linspace(
417 np.min(radius),
418 np.max(radius),
419 radial_gradient_n_spline_nodes,
420 )
422 # Put in a normalization for fitting.
423 norm = np.nanpercentile(binned["value"], 95.0)
425 fitter = FlatGradientFitter(
426 nodes,
427 xf,
428 yf,
429 binned["value"]/norm,
430 np.array([]),
431 constrain_zero=False,
432 fit_centroid=False,
433 fit_gradient=False,
434 fp_centroid_x=0.0,
435 fp_centroid_y=0.0,
436 )
437 p0 = fitter.compute_p0()
438 pars = fitter.fit(p0)
440 binned["value"] /= fitter.compute_model(pars)
442 # Fit gain ratios.
443 fitter = FlatGainRatioFitter(
444 detector.getBBox(),
445 chebyshev_gradient_order,
446 binned["xd"],
447 binned["yd"],
448 binned["amp_index"],
449 binned["value"],
450 np.arange(n_amp),
451 fixed_amp_index,
452 )
453 pars = fitter.fit(nsig_clip=nsig_clip)
454 amp_pars = pars[fitter.indices["amp_pars"]]
456 fractional_gain_ratio = np.abs(amp_pars - 1.0)
457 max_ratio_ind = np.argmax(fractional_gain_ratio)
459 if fractional_gain_ratio[max_ratio_ind] > max_fractional_gain_ratio:
460 log.warning(
461 "Found bad amp %s with offset parameter %.2f",
462 amp_names[max_ratio_ind],
463 amp_pars[max_ratio_ind],
464 )
465 good = (binned["amp_index"] != max_ratio_ind)
466 binned = binned[good]
467 binned["value"] = value_uncorrected[good]
469 value_uncorrected = binned["value"].copy()
470 else:
471 bad_amps_converged = True
473 n_iter += 1
475 return amp_pars