Coverage for python / lsst / cp / pipe / cpLinearitySolve.py: 10%
872 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-15 07:22 -0300
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#
23__all__ = [
24 "LinearitySolveTask",
25 "LinearitySolveConfig",
26 "LinearityDoubleSplineSolveTask",
27 "LinearityDoubleSplineSolveConfig",
28]
30import copy
31import logging
32import numpy as np
33import esutil
34from scipy.stats import median_abs_deviation
35from scipy.interpolate import Akima1DInterpolator
37import lsst.afw.image as afwImage
38import lsst.pipe.base as pipeBase
39import lsst.pipe.base.connectionTypes as cT
40import lsst.pex.config as pexConfig
42from lsstDebug import getDebugFrame
43from lsst.ip.isr import (Linearizer, IsrProvenance)
45from .utils import (funcPolynomial, irlsFit, AstierSplineLinearityFitter,
46 extractCalibDate)
49def ptcLookup(datasetType, registry, quantumDataId, collections):
50 """Butler lookup function to allow PTC to be found.
52 Parameters
53 ----------
54 datasetType : `lsst.daf.butler.DatasetType`
55 Dataset type to look up.
56 registry : `lsst.daf.butler.Registry`
57 Registry for the data repository being searched.
58 quantumDataId : `lsst.daf.butler.DataCoordinate`
59 Data ID for the quantum of the task this dataset will be passed to.
60 This must include an "instrument" key, and should also include any
61 keys that are present in ``datasetType.dimensions``. If it has an
62 ``exposure`` or ``visit`` key, that's a sign that this function is
63 not actually needed, as those come with the temporal information that
64 would allow a real validity-range lookup.
65 collections : `lsst.daf.butler.registry.CollectionSearch`
66 Collections passed by the user when generating a QuantumGraph. Ignored
67 by this function (see notes below).
69 Returns
70 -------
71 refs : `list` [ `DatasetRef` ]
72 A zero- or single-element list containing the matching
73 dataset, if one was found.
75 Raises
76 ------
77 RuntimeError
78 Raised if more than one PTC reference is found.
79 """
80 refs = list(registry.queryDatasets(datasetType, dataId=quantumDataId, collections=collections,
81 findFirst=False))
82 if len(refs) >= 2:
83 RuntimeError("Too many PTC connections found. Incorrect collections supplied?")
85 return refs
88class LinearitySolveConnections(pipeBase.PipelineTaskConnections,
89 dimensions=("instrument", "detector")):
90 dummy = cT.Input(
91 name="raw",
92 doc="Dummy exposure.",
93 storageClass='Exposure',
94 dimensions=("instrument", "exposure", "detector"),
95 multiple=True,
96 deferLoad=True,
97 )
99 camera = cT.PrerequisiteInput(
100 name="camera",
101 doc="Camera Geometry definition.",
102 storageClass="Camera",
103 dimensions=("instrument", ),
104 isCalibration=True,
105 )
107 inputPtc = cT.PrerequisiteInput(
108 name="ptc",
109 doc="Input PTC dataset.",
110 storageClass="PhotonTransferCurveDataset",
111 dimensions=("instrument", "detector"),
112 isCalibration=True,
113 lookupFunction=ptcLookup,
114 )
116 inputLinearizerPtc = cT.Input(
117 name="linearizerPtc",
118 doc="Input linearizer PTC dataset.",
119 storageClass="PhotonTransferCurveDataset",
120 dimensions=("instrument", "detector"),
121 isCalibration=True,
122 )
124 inputPhotodiodeCorrection = cT.Input(
125 name="pdCorrection",
126 doc="Input photodiode correction.",
127 storageClass="IsrCalib",
128 dimensions=("instrument", ),
129 isCalibration=True,
130 )
132 inputNormalization = cT.Input(
133 name="cpLinearizerPtcNormalization",
134 doc="Focal-plane normalization table.",
135 storageClass="ArrowAstropy",
136 dimensions=("instrument",),
137 isCalibration=True,
138 )
140 outputLinearizer = cT.Output(
141 name="linearity",
142 doc="Output linearity measurements.",
143 storageClass="Linearizer",
144 dimensions=("instrument", "detector"),
145 isCalibration=True,
146 )
148 def __init__(self, *, config=None):
149 super().__init__(config=config)
151 if not config.applyPhotodiodeCorrection:
152 del self.inputPhotodiodeCorrection
154 if config.useLinearizerPtc:
155 del self.inputPtc
156 else:
157 del self.inputLinearizerPtc
159 if not config.useFocalPlaneNormalization:
160 del self.inputNormalization
163class LinearitySolveConfig(pipeBase.PipelineTaskConfig,
164 pipelineConnections=LinearitySolveConnections):
165 """Configuration for solving the linearity from PTC dataset.
166 """
167 linearityType = pexConfig.ChoiceField(
168 dtype=str,
169 doc="Type of linearizer to construct.",
170 default="Squared",
171 allowed={
172 "LookupTable": "Create a lookup table solution.",
173 "Polynomial": "Create an arbitrary polynomial solution.",
174 "Squared": "Create a single order squared solution.",
175 "Spline": "Create a spline based solution.",
176 "None": "Create a dummy solution.",
177 }
178 )
179 polynomialOrder = pexConfig.RangeField(
180 dtype=int,
181 doc="Degree of polynomial to fit. Must be at least 2.",
182 default=3,
183 min=2,
184 )
185 splineKnots = pexConfig.Field(
186 dtype=int,
187 doc="Number of spline knots to use in fit.",
188 default=10,
189 )
191 trimmedState = pexConfig.Field(
192 dtype=bool,
193 doc="Will this linearizer be used on trimmed data?",
194 default=True,
195 )
197 maxLookupTableAdu = pexConfig.Field(
198 dtype=int,
199 doc="Maximum DN value for a LookupTable linearizer.",
200 default=2**18,
201 )
202 maxLinearAdu = pexConfig.Field(
203 dtype=float,
204 doc="Maximum adu value to use to estimate linear term; not used with spline fits.",
205 default=20000.0,
206 )
207 minLinearAdu = pexConfig.Field(
208 dtype=float,
209 doc="Minimum adu value to use to estimate linear term.",
210 default=30.0,
211 )
212 nSigmaClipLinear = pexConfig.Field(
213 dtype=float,
214 doc="Maximum deviation from linear solution for Poissonian noise.",
215 default=5.0,
216 )
217 ignorePtcMask = pexConfig.Field(
218 dtype=bool,
219 doc="Ignore the expIdMask set by the PTC solver?",
220 default=False,
221 deprecated="This field is no longer used. Will be removed after v28.",
222 )
223 maxFracLinearityDeviation = pexConfig.Field(
224 dtype=float,
225 doc="Maximum fraction deviation from raw linearity to compute "
226 "linearityTurnoff and linearityMaxSignal.",
227 # TODO: DM-46811 investigate if this can be raised to 0.05.
228 default=0.01,
229 )
230 minSignalFitLinearityTurnoff = pexConfig.Field(
231 dtype=float,
232 doc="Minimum signal to compute raw linearity slope for linearityTurnoff.",
233 default=1000.0,
234 )
235 usePhotodiode = pexConfig.Field(
236 dtype=bool,
237 doc="Use the photodiode info instead of the raw expTimes?",
238 default=False,
239 )
240 applyPhotodiodeCorrection = pexConfig.Field(
241 dtype=bool,
242 doc="Calculate and apply a correction to the photodiode readings?",
243 default=False,
244 )
245 minPhotodiodeCurrent = pexConfig.Field(
246 dtype=float,
247 doc="Minimum value to trust photodiode signals.",
248 default=0.0,
249 )
250 doAutoGrouping = pexConfig.Field(
251 dtype=bool,
252 doc="Do automatic group detection? Cannot be True if splineGroupingColumn is also set. "
253 "The automatic group detection will use the ratio of signal to exposure time (if "
254 "autoGroupingUseExptime is True) or photodiode (if False) to determine which "
255 "flat pairs were taken with different illumination settings.",
256 default=False,
257 )
258 autoGroupingUseExptime = pexConfig.Field(
259 dtype=bool,
260 doc="Use exposure time to determine automatic grouping. Used if doAutoGrouping=True.",
261 default=True,
262 )
263 autoGroupingThreshold = pexConfig.Field(
264 dtype=float,
265 doc="Minimum relative jump from sorted conversion values to determine a group.",
266 default=0.1,
267 )
268 autoGroupingMaxSignalFraction = pexConfig.Field(
269 dtype=float,
270 doc="Only do auto-grouping when the signal is this fraction of the maximum signal. "
271 "All exposures with signal higher than this threshold will be put into the "
272 "largest signal group. This config is needed if the input PTC goes beyond "
273 "the linearity turnoff.",
274 default=0.9,
275 )
276 splineGroupingColumn = pexConfig.Field(
277 dtype=str,
278 doc="Column to use for grouping together points for Spline mode, to allow "
279 "for different proportionality constants. If None, then grouping will "
280 "only be done if doAutoGrouping is True.",
281 default=None,
282 optional=True,
283 )
284 splineGroupingMinPoints = pexConfig.Field(
285 dtype=int,
286 doc="Minimum number of linearity points to allow grouping together points "
287 "for Spline mode with splineGroupingColumn. This configuration is here "
288 "to prevent misuse of the Spline code to avoid over-fitting.",
289 default=100,
290 )
291 splineFitMinIter = pexConfig.Field(
292 dtype=int,
293 doc="Minimum number of iterations for spline fit.",
294 default=3,
295 )
296 splineFitMaxIter = pexConfig.Field(
297 dtype=int,
298 doc="Maximum number of iterations for spline fit.",
299 default=20,
300 )
301 splineFitMaxRejectionPerIteration = pexConfig.Field(
302 dtype=int,
303 doc="Maximum number of rejections per iteration for spline fit.",
304 default=5,
305 )
306 doSplineFitOffset = pexConfig.Field(
307 dtype=bool,
308 doc="Fit a scattered light offset in the spline fit.",
309 default=True,
310 )
311 doSplineFitWeights = pexConfig.Field(
312 dtype=bool,
313 doc="Fit linearity weight parameters in the spline fit.",
314 default=False,
315 )
316 splineFitWeightParsStart = pexConfig.ListField(
317 dtype=float,
318 doc="Starting parameters for weight fit, if doSplineFitWeights=True. "
319 "Parameters are such that sigma = sqrt(par[0]**2. + par[1]**2./mu)."
320 "If doSplineFitWeights=False then these are used as-is; otherwise "
321 "they are used as the initial values for fitting these parameters.",
322 length=2,
323 default=[1.0, 0.0],
324 )
325 doSplineFitTemperature = pexConfig.Field(
326 dtype=bool,
327 doc="Fit temperature coefficient in spline fit?",
328 default=False,
329 )
330 splineFitTemperatureColumn = pexConfig.Field(
331 dtype=str,
332 doc="Name of the temperature column to use when fitting temperature "
333 "coefficients in spline fit; this must not be None if "
334 "doSplineFitTemperature is True.",
335 default=None,
336 optional=True,
337 )
338 doSplineFitTemporal = pexConfig.Field(
339 dtype=bool,
340 doc="Fit a linear temporal parameter coefficient in spline fit?",
341 default=False,
342 )
343 useLinearizerPtc = pexConfig.Field(
344 dtype=bool,
345 doc="Use a linearizer ptc in a single pipeline?",
346 default=False,
347 )
348 useFocalPlaneNormalization = pexConfig.Field(
349 dtype=bool,
350 doc="Use focal-plane normalization in addition to/instead of photodiode? "
351 "(Only used with spline fitting).",
352 default=False,
353 )
355 def validate(self):
356 super().validate()
358 if self.doSplineFitTemperature and self.splineFitTemperatureColumn is None:
359 raise ValueError("Must set splineFitTemperatureColumn if doSplineFitTemperature is True.")
361 if self.doAutoGrouping and self.splineGroupingColumn is not None:
362 raise ValueError("Must not set doAutoGrouping and also splineGroupingColumn")
363 if self.doAutoGrouping:
364 if not self.autoGroupingUseExptime and not self.usePhotodiode:
365 raise ValueError("If doAutoGrouping is True and autoGroupingUseExptime is False, then "
366 "usePhotodiode must be True.")
369class LinearitySolveTask(pipeBase.PipelineTask):
370 """Fit the linearity from the PTC dataset.
371 """
373 ConfigClass = LinearitySolveConfig
374 _DefaultName = 'cpLinearitySolve'
376 def runQuantum(self, butlerQC, inputRefs, outputRefs):
377 """Ensure that the input and output dimensions are passed along.
379 Parameters
380 ----------
381 butlerQC : `lsst.daf.butler.QuantumContext`
382 Butler to operate on.
383 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
384 Input data refs to load.
385 ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection`
386 Output data refs to persist.
387 """
388 inputs = butlerQC.get(inputRefs)
390 # Use the dimensions to set calib/provenance information.
392 if self.config.useLinearizerPtc:
393 inputs["inputDims"] = dict(inputRefs.inputLinearizerPtc.dataId.required)
394 inputs["inputPtc"] = inputs["inputLinearizerPtc"]
395 del inputs["inputLinearizerPtc"]
396 else:
397 inputs["inputDims"] = dict(inputRefs.inputPtc.dataId.required)
399 # Add calibration provenance info to header.
400 kwargs = dict()
401 reference = getattr(inputRefs, "inputPtc", None)
403 if reference is not None and hasattr(reference, "run"):
404 runKey = "PTC_RUN"
405 runValue = reference.run
406 idKey = "PTC_UUID"
407 idValue = str(reference.id)
408 dateKey = "PTC_DATE"
409 calib = inputs.get("inputPtc", None)
410 dateValue = extractCalibDate(calib)
412 kwargs[runKey] = runValue
413 kwargs[idKey] = idValue
414 kwargs[dateKey] = dateValue
416 self.log.info("Using " + str(reference.run))
418 outputs = self.run(**inputs)
419 outputs.outputLinearizer.updateMetadata(setDate=False, **kwargs)
421 butlerQC.put(outputs, outputRefs)
423 def run(self, inputPtc, dummy, camera, inputDims,
424 inputPhotodiodeCorrection=None, inputNormalization=None):
425 """Fit non-linearity to PTC data, returning the correct Linearizer
426 object.
428 Parameters
429 ----------
430 inputPtc : `lsst.ip.isr.PtcDataset`
431 Pre-measured PTC dataset.
432 dummy : `lsst.afw.image.Exposure`
433 The exposure used to select the appropriate PTC dataset.
434 In almost all circumstances, one of the input exposures
435 used to generate the PTC dataset is the best option.
436 camera : `lsst.afw.cameraGeom.Camera`
437 Camera geometry.
438 inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
439 DataIds to use to populate the output calibration.
440 inputPhotodiodeCorrection :
441 `lsst.ip.isr.PhotodiodeCorrection`, optional
442 Pre-measured photodiode correction used in the case when
443 applyPhotodiodeCorrection=True.
444 inputNormalization : `astropy.table.Table`, optional
445 Focal plane normalization table to use if
446 useFocalPlaneNormalization is True.
448 Returns
449 -------
450 results : `lsst.pipe.base.Struct`
451 The results struct containing:
453 ``outputLinearizer``
454 Final linearizer calibration (`lsst.ip.isr.Linearizer`).
455 ``outputProvenance``
456 Provenance data for the new calibration
457 (`lsst.ip.isr.IsrProvenance`).
459 Notes
460 -----
461 This task currently fits only polynomial-defined corrections,
462 where the correction coefficients are defined such that:
463 :math:`corrImage = uncorrImage + \\sum_i c_i uncorrImage^(2 + i)`
464 These :math:`c_i` are defined in terms of the direct polynomial fit:
465 :math:`meanVector ~ P(x=timeVector) = \\sum_j k_j x^j`
466 such that :math:`c_(j-2) = -k_j/(k_1^j)` in units of DN^(1-j) (c.f.,
467 Eq. 37 of 2003.05978). The `config.polynomialOrder` or
468 `config.splineKnots` define the maximum order of :math:`x^j` to fit.
469 As :math:`k_0` and :math:`k_1` are degenerate with bias level and gain,
470 they are not included in the non-linearity correction.
471 """
472 if len(dummy) == 0:
473 self.log.warning("No dummy exposure found.")
475 detector = camera[inputDims['detector']]
476 if self.config.linearityType == 'LookupTable':
477 table = np.zeros((len(detector), self.config.maxLookupTableAdu), dtype=np.float32)
478 tableIndex = 0
479 else:
480 table = None
481 tableIndex = None # This will fail if we increment it.
483 # Initialize the linearizer.
484 linearizer = Linearizer(detector=detector, table=table, log=self.log)
485 linearizer.updateMetadataFromExposures([inputPtc])
486 if self.config.usePhotodiode and self.config.applyPhotodiodeCorrection:
487 abscissaCorrections = inputPhotodiodeCorrection.abscissaCorrections
489 groupingValues = _determineInputGroups(
490 inputPtc,
491 self.config.doAutoGrouping,
492 self.config.autoGroupingUseExptime,
493 self.config.autoGroupingMaxSignalFraction,
494 self.config.autoGroupingThreshold,
495 self.config.splineGroupingColumn,
496 self.config.minPhotodiodeCurrent,
497 )
499 if self.config.linearityType == 'Spline':
500 if self.config.doSplineFitTemperature:
501 if self.config.splineFitTemperatureColumn not in inputPtc.auxValues:
502 raise ValueError("Config requests fitting temperature coefficient for "
503 f"{self.config.splineFitTemperatureColumn} but this column "
504 "is not available in inputPtc.auxValues.")
505 temperatureValues = inputPtc.auxValues[self.config.splineFitTemperatureColumn]
506 else:
507 temperatureValues = None
509 # We set this to have a value to fill the bad amps.
510 fitOrder = self.config.splineKnots
511 else:
512 fitOrder = self.config.polynomialOrder
514 for i, amp in enumerate(detector):
515 ampName = amp.getName()
517 # Save the input gains
518 linearizer.inputGain[ampName] = inputPtc.gain[ampName]
520 if ampName in inputPtc.badAmps:
521 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
522 self.log.warning("Amp %s in detector %s has no usable PTC information. Skipping!",
523 ampName, detector.getName())
524 continue
526 # Check for too few points.
527 if self.config.linearityType == "Spline" \
528 and self.config.splineGroupingColumn is not None \
529 and len(inputPtc.inputExpIdPairs[ampName]) < self.config.splineGroupingMinPoints:
530 raise RuntimeError(
531 "The input PTC has too few points to reliably run with PD grouping. "
532 "The recommended course of action is to set splineGroupingColumn to None. "
533 "If you really know what you are doing, you may reduce "
534 "config.splineGroupingMinPoints.")
536 # We start with all finite values.
537 mask = np.isfinite(inputPtc.rawMeans[ampName])
539 if self.config.linearityType == "Spline" and temperatureValues is not None:
540 mask &= np.isfinite(temperatureValues)
542 if self.config.usePhotodiode:
543 modExpTimes = inputPtc.photoCharges[ampName].copy()
544 # Make sure any exposure pairs that do not have photodiode data
545 # are masked.
546 mask[~np.isfinite(modExpTimes)] = False
548 # Make sure any photodiode measurements below the configured
549 # minimum are masked.
550 mask[modExpTimes < self.config.minPhotodiodeCurrent] = False
552 # Get the photodiode correction.
553 if self.config.applyPhotodiodeCorrection:
554 for j, pair in enumerate(inputPtc.inputExpIdPairs[ampName]):
555 try:
556 correction = abscissaCorrections[str(pair)]
557 except KeyError:
558 correction = 0.0
559 modExpTimes[j] += correction
561 inputAbscissa = modExpTimes
562 else:
563 inputAbscissa = inputPtc.rawExpTimes[ampName].copy()
565 # Normalize if configured.
566 inputNorm = np.ones_like(inputAbscissa)
567 if self.config.useFocalPlaneNormalization:
568 exposures = np.asarray(inputPtc.inputExpIdPairs[ampName])[:, 0]
569 a, b = esutil.numpy_util.match(exposures, inputNormalization["exposure"])
570 inputNorm[a] = inputNormalization["normalization"][b]
571 inputAbscissa *= inputNorm
573 # Compute linearityTurnoff and linearitySignalMax.
574 turnoffMask = inputPtc.expIdMask[ampName].copy()
575 turnoffMask &= mask
577 turnoffIndex, turnoff, maxSignal, _ = _computeTurnoffAndMax(
578 inputAbscissa,
579 inputPtc.rawMeans[ampName],
580 turnoffMask,
581 groupingValues,
582 ampName=ampName,
583 minSignalFitLinearityTurnoff=self.config.minSignalFitLinearityTurnoff,
584 maxFracLinearityDeviation=self.config.maxFracLinearityDeviation,
585 log=self.log,
586 use_all_for_normalization=True,
587 )
589 if np.isnan(turnoff):
590 # This is a bad amp, with no linearizer.
591 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
592 self.log.warning("Amp %s in detector %s has no usable linearizer information. Skipping!",
593 ampName, detector.getName())
594 continue
596 linearizer.linearityTurnoff[ampName] = turnoff
597 linearizer.linearityMaxSignal[ampName] = maxSignal
599 inputOrdinate = inputPtc.rawMeans[ampName].copy()
601 linearizer.inputAbscissa[ampName] = inputAbscissa.copy()
602 linearizer.inputOrdinate[ampName] = inputOrdinate.copy()
603 linearizer.inputGroupingIndex[ampName] = groupingValues.copy()
604 linearizer.inputNormalization[ampName] = inputNorm.copy()
606 if self.config.linearityType != 'Spline':
607 mask &= (inputOrdinate < self.config.maxLinearAdu)
608 else:
609 # For spline fits, cut above the turnoff.
610 self.log.info("Using linearityTurnoff of %.4f adu for amplifier %s", turnoff, ampName)
611 extraMask = np.ones(len(inputOrdinate), dtype=bool)
612 extraMask[turnoffIndex + 1:] = False
613 mask &= extraMask
615 mask &= (inputOrdinate > self.config.minLinearAdu)
617 # Initial value for the mask.
618 linearizer.inputMask[ampName] = mask.copy()
620 if mask.sum() < 2:
621 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
622 self.log.warning("Amp %s in detector %s has not enough points for fit. Skipping!",
623 ampName, detector.getName())
624 continue
626 if self.config.linearityType != 'Spline':
627 linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 100.0], inputAbscissa[mask],
628 inputOrdinate[mask], funcPolynomial)
630 # Convert this proxy-to-flux fit into an expected linear flux
631 linearOrdinate = linearFit[0] + linearFit[1] * inputAbscissa
632 # Exclude low end outliers.
633 # This is compared to the original values.
634 threshold = self.config.nSigmaClipLinear * np.sqrt(abs(inputOrdinate))
636 mask[np.abs(inputOrdinate - linearOrdinate) >= threshold] = False
638 linearizer.inputMask[ampName] = mask.copy()
640 if mask.sum() < 2:
641 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
642 self.log.warning("Amp %s in detector %s has not enough points in linear ordinate. "
643 "Skipping!", ampName, detector.getName())
644 continue
646 self.debugFit('linearFit', inputAbscissa, inputOrdinate, linearOrdinate, mask, ampName)
648 # Do fits
649 if self.config.linearityType in ['Polynomial', 'Squared', 'LookupTable']:
650 polyFit = np.zeros(fitOrder + 1)
651 polyFit[1] = 1.0
652 polyFit, polyFitErr, chiSq, weights = irlsFit(polyFit, linearOrdinate[mask],
653 inputOrdinate[mask], funcPolynomial)
655 # Truncate the polynomial fit to the squared term.
656 k1 = polyFit[1]
657 linearityCoeffs = np.array(
658 [-coeff/(k1**order) for order, coeff in enumerate(polyFit)]
659 )[2:]
660 significant = np.where(np.abs(linearityCoeffs) > 1e-10)
661 self.log.info("Significant polynomial fits: %s", significant)
663 modelOrdinate = funcPolynomial(polyFit, linearOrdinate)
665 self.debugFit(
666 'polyFit',
667 inputAbscissa[mask],
668 inputOrdinate[mask],
669 modelOrdinate[mask],
670 None,
671 ampName,
672 )
674 if self.config.linearityType == 'Squared':
675 # The first term is the squared term.
676 linearityCoeffs = linearityCoeffs[0: 1]
677 elif self.config.linearityType == 'LookupTable':
678 # Use linear part to get time at which signal is
679 # maxAduForLookupTableLinearizer DN
680 tMax = (self.config.maxLookupTableAdu - polyFit[0])/polyFit[1]
681 timeRange = np.linspace(0, tMax, self.config.maxLookupTableAdu)
682 signalIdeal = polyFit[0] + polyFit[1]*timeRange
683 signalUncorrected = funcPolynomial(polyFit, timeRange)
684 lookupTableRow = signalIdeal - signalUncorrected # LinearizerLookupTable has correction
686 linearizer.tableData[tableIndex, :] = lookupTableRow
687 linearityCoeffs = np.array([tableIndex, 0])
688 tableIndex += 1
689 elif self.config.linearityType in ['Spline']:
690 # This is a spline fit with photodiode data based on a model
691 # from Pierre Astier.
692 # This model fits a spline with (optional) nuisance parameters
693 # to allow for different linearity coefficients with different
694 # photodiode settings. The minimization is a least-squares
695 # fit with the residual of
696 # Sum[(S(mu_i) + mu_i - O)/(k_j * D_i) - 1]**2, where S(mu_i)
697 # is an Akima Spline function of mu_i, the observed flat-pair
698 # mean; D_j is the photo-diode measurement corresponding to
699 # that flat-pair; and k_j is a constant of proportionality
700 # which is over index j as it is allowed to
701 # be different based on different photodiode settings (e.g.
702 # CCOBCURR); and O is a constant offset to allow for light
703 # leaks (and is only fit if doSplineFitOffset=True). In
704 # addition, if config.doSplineFitTemperature is True then
705 # the fit will adjust mu such that
706 # mu = mu_input*(1 + alpha*(T - T_ref))
707 # and T_ref is taken as the median temperature of the run.
708 # Finally, if config.doSplineFitTemporal is True then the
709 # fit will further adjust mu such that
710 # mu = mu_input*(1 + beta*(mjd - mjd_ref))
711 # and mjd_ref is taken as the median mjd of the run.
712 # Note that this fit is only valid if the input data
713 # was taken with a randomly shuffled order of exposure
714 # levels.
716 # The fit has additional constraints to ensure that the spline
717 # goes through the (0, 0) point, as well as a normalization
718 # condition so that the average of the spline over the full
719 # range is 0. The normalization ensures that the spline only
720 # fits deviations from linearity, rather than the linear
721 # function itself which is degenerate with the gain.
723 # We want to make sure the top node is above the top value
724 # to avoid edge issues with the top point.
725 nodes = np.linspace(0.0, np.max(inputOrdinate[mask]) + 1.0, self.config.splineKnots)
727 if temperatureValues is not None:
728 temperatureValuesScaled = temperatureValues - np.median(temperatureValues[mask])
729 else:
730 temperatureValuesScaled = None
732 if self.config.doSplineFitTemporal:
733 inputMjdScaled = inputPtc.inputExpPairMjdStartList[ampName].copy()
734 inputMjdScaled -= np.nanmedian(inputMjdScaled)
735 else:
736 inputMjdScaled = None
738 fitter = AstierSplineLinearityFitter(
739 nodes,
740 groupingValues,
741 inputAbscissa,
742 inputOrdinate,
743 mask=mask,
744 log=self.log,
745 fit_offset=self.config.doSplineFitOffset,
746 fit_weights=self.config.doSplineFitWeights,
747 weight_pars_start=self.config.splineFitWeightParsStart,
748 fit_temperature=self.config.doSplineFitTemperature,
749 temperature_scaled=temperatureValuesScaled,
750 max_signal_nearly_linear=inputPtc.ptcTurnoff[ampName],
751 fit_temporal=self.config.doSplineFitTemporal,
752 mjd_scaled=inputMjdScaled,
753 )
754 p0 = fitter.estimate_p0(use_all_for_normalization=True)
755 pars = fitter.fit(
756 p0,
757 min_iter=self.config.splineFitMinIter,
758 max_iter=self.config.splineFitMaxIter,
759 max_rejection_per_iteration=self.config.splineFitMaxRejectionPerIteration,
760 n_sigma_clip=self.config.nSigmaClipLinear,
761 )
763 # Confirm that the first parameter is 0, and set it to
764 # exactly zero.
765 if not np.isclose(pars[0], 0):
766 raise RuntimeError("Programmer error! First spline parameter must "
767 "be consistent with zero.")
768 pars[0] = 0.0
770 linearityChisq = fitter.compute_chisq_dof(pars)
772 linearityCoeffs = np.concatenate([nodes, pars[fitter.par_indices["values"]]])
773 linearFit = np.array([0.0, np.mean(pars[fitter.par_indices["groups"]])])
775 # We must modify the inputOrdinate according to the
776 # nuisance terms in the linearity fit for the residual
777 # computation code to work properly.
778 # The true mu (inputOrdinate) is given by
779 # mu = mu_in * (1 + alpha*t_scale) * (1 + beta*mjd_scale)
780 if self.config.doSplineFitTemperature:
781 inputOrdinate *= (1.0
782 + pars[fitter.par_indices["temperature_coeff"]]*temperatureValuesScaled)
783 if self.config.doSplineFitTemporal:
784 inputOrdinate *= (1.0
785 + pars[fitter.par_indices["temporal_coeff"]]*inputMjdScaled)
786 # We have to adjust the abscissa for the different groups.
787 # This is because we need a corrected abscissa to get a
788 # reasonable linear fit to look for residuals, particularly in
789 # the case of significantly different signal-vs-photodiode or
790 # signal-vs-exptime scalings for different groups. This then
791 # becomes a multiplication by the relative scaling of the
792 # different groups.
793 for j, group_index in enumerate(fitter.group_indices):
794 inputAbscissa[group_index] *= (pars[fitter.par_indices["groups"][j]] / linearFit[1])
795 # And remove the offset term.
796 if self.config.doSplineFitOffset:
797 inputOrdinate -= pars[fitter.par_indices["offset"]]
799 linearOrdinate = linearFit[1] * inputOrdinate
800 # For the spline fit, reuse the "polyFit -> fitParams"
801 # field to record the linear coefficients for the groups.
802 # We additionally append the offset and weight_pars;
803 # however these will be zero-length arrays if these were
804 # not configured to be fit.
805 polyFit = np.concatenate((
806 pars[fitter.par_indices["groups"]],
807 pars[fitter.par_indices["offset"]],
808 pars[fitter.par_indices["weight_pars"]],
809 pars[fitter.par_indices["temperature_coeff"]],
810 pars[fitter.par_indices["temporal_coeff"]],
811 ))
812 polyFitErr = np.zeros_like(polyFit)
813 chiSq = linearityChisq
815 # Update mask based on what the fitter rejected.
816 mask = fitter.mask
818 linearizer.inputMask[ampName] = mask.copy()
819 else:
820 polyFit = np.zeros(1)
821 polyFitErr = np.zeros(1)
822 chiSq = np.nan
823 linearityCoeffs = np.zeros(1)
825 linearizer.linearityType[ampName] = self.config.linearityType
826 linearizer.linearityCoeffs[ampName] = linearityCoeffs
827 if self.config.trimmedState:
828 linearizer.linearityBBox[ampName] = amp.getBBox()
829 else:
830 linearizer.linearityBBox[ampName] = amp.getRawBBox()
831 linearizer.fitParams[ampName] = polyFit
832 linearizer.fitParamsErr[ampName] = polyFitErr
833 linearizer.fitChiSq[ampName] = chiSq
834 linearizer.linearFit[ampName] = linearFit
836 image = afwImage.ImageF(len(inputOrdinate), 1)
837 image.array[:, :] = inputOrdinate
838 linearizeFunction = linearizer.getLinearityTypeByName(linearizer.linearityType[ampName])
839 linearizeFunction()(
840 image,
841 **{'coeffs': linearizer.linearityCoeffs[ampName],
842 'table': linearizer.tableData,
843 'log': linearizer.log}
844 )
845 linearizeModel = image.array[0, :]
847 # The residuals that we record are the final residuals compared to
848 # a linear model, after everything has been (properly?) linearized.
849 if mask.sum() < 2:
850 self.log.warning("Amp %s in detector %s has not enough points in linear ordinate "
851 "for residuals. Skipping!", ampName, detector.getName())
852 residuals = np.full_like(linearizeModel, np.nan)
853 residualsUnmasked = residuals.copy()
854 else:
855 postLinearFit, _, _, _ = irlsFit(
856 linearFit,
857 inputAbscissa[mask],
858 linearizeModel[mask],
859 funcPolynomial,
860 )
861 # When computing residuals, we only care about the slope of
862 # the postLinearFit and not the intercept. The intercept
863 # itself depends on a possibly unknown zero in the abscissa
864 # (often photodiode) which may have an arbitrary value.
865 residuals = linearizeModel - (postLinearFit[1] * inputAbscissa)
866 residualsUnmasked = residuals.copy()
867 # We set masked residuals to nan.
868 residuals[~mask] = np.nan
870 linearizer.fitResidualsUnmasked[ampName] = residualsUnmasked
871 linearizer.fitResiduals[ampName] = residuals
872 linearizer.fitResidualsModel[ampName] = linearizeModel.copy()
874 finite = np.isfinite(residuals)
875 if finite.sum() == 0:
876 sigmad = np.nan
877 else:
878 sigmad = median_abs_deviation(residuals[finite]/inputOrdinate[finite], scale="normal")
879 linearizer.fitResidualsSigmaMad[ampName] = sigmad
881 self.debugFit(
882 'solution',
883 inputOrdinate[mask],
884 linearOrdinate[mask],
885 linearizeModel[mask],
886 None,
887 ampName,
888 )
890 self.fixupBadAmps(linearizer)
892 linearizer.hasLinearity = True
893 linearizer.validate()
894 linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
895 linearizer.updateMetadata(setDate=True, setCalibId=True)
896 linearizer.updateMetadataFromExposures([inputPtc])
897 provenance = IsrProvenance(calibType='linearizer')
899 return pipeBase.Struct(
900 outputLinearizer=linearizer,
901 outputProvenance=provenance,
902 )
904 def fillBadAmp(self, linearizer, fitOrder, inputPtc, amp):
905 # Need to fill linearizer with empty values
906 # if the amp is non-functional
907 ampName = amp.getName()
908 nEntries = 1
909 pEntries = 1
910 if self.config.linearityType in ['Polynomial']:
911 # We discard the first 2 entries in the polynomial.
912 nEntries = fitOrder + 1 - 2
913 pEntries = fitOrder + 1 - 2
914 elif self.config.linearityType in ['Spline']:
915 nEntries = fitOrder * 2
916 elif self.config.linearityType in ['Squared', 'None']:
917 nEntries = 1
918 pEntries = fitOrder + 1
919 elif self.config.linearityType in ['LookupTable']:
920 nEntries = 2
921 pEntries = fitOrder + 1
923 nPair = len(inputPtc.inputExpIdPairs[ampName])
925 linearizer.linearityType[ampName] = "None"
926 linearizer.linearityCoeffs[ampName] = np.zeros(nEntries)
927 if self.config.trimmedState:
928 linearizer.linearityBBox[ampName] = amp.getBBox()
929 else:
930 linearizer.linearityBBox[ampName] = amp.getRawBBox()
931 linearizer.fitParams[ampName] = np.zeros(pEntries)
932 linearizer.fitParamsErr[ampName] = np.zeros(pEntries)
933 linearizer.fitChiSq[ampName] = np.nan
934 linearizer.fitResiduals[ampName] = np.zeros(nPair)
935 linearizer.fitResidualsSigmaMad[ampName] = np.nan
936 linearizer.fitResidualsUnmasked[ampName] = np.zeros(nPair)
937 linearizer.fitResidualsModel[ampName] = np.zeros(nPair)
938 linearizer.linearFit[ampName] = np.zeros(2)
939 linearizer.linearityTurnoff[ampName] = np.nan
940 linearizer.linearityMaxSignal[ampName] = np.nan
941 linearizer.inputMask[ampName] = np.zeros(nPair, dtype=np.bool_)
942 linearizer.inputAbscissa[ampName] = np.zeros(nPair)
943 linearizer.inputOrdinate[ampName] = np.zeros(nPair)
944 linearizer.inputGroupingIndex[ampName] = np.zeros(nPair, dtype=np.int64)
945 linearizer.inputNormalization[ampName] = np.ones(nPair)
947 return linearizer
949 def fixupBadAmps(self, linearizer):
950 """Fix nan padding in bad amplifiers.
952 Parameters
953 ----------
954 linearizer : `lsst.ip.isr.Linearizer`
955 """
956 fitParamsMaxLen = 0
957 for ampName in linearizer.ampNames:
958 if (length := len(linearizer.fitParams[ampName])) > fitParamsMaxLen:
959 fitParamsMaxLen = length
961 for ampName in linearizer.ampNames:
962 if linearizer.linearityType[ampName] == "None":
963 # Bad amplifier.
964 linearizer.fitParams[ampName] = np.zeros(fitParamsMaxLen)
965 linearizer.fitParamsErr[ampName] = np.zeros(fitParamsMaxLen)
966 elif len(linearizer.fitParams[ampName]) != fitParamsMaxLen:
967 raise RuntimeError("Linearity has mismatched fitParams; check code/data.")
969 def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName):
970 """Debug method for linearity fitting.
972 Parameters
973 ----------
974 stepname : `str`
975 A label to use to check if we care to debug at a given
976 line of code.
977 xVector : `numpy.array`, (N,)
978 The values to use as the independent variable in the
979 linearity fit.
980 yVector : `numpy.array`, (N,)
981 The values to use as the dependent variable in the
982 linearity fit.
983 yModel : `numpy.array`, (N,)
984 The values to use as the linearized result.
985 mask : `numpy.array` [`bool`], (N,) , optional
986 A mask to indicate which entries of ``xVector`` and
987 ``yVector`` to keep.
988 ampName : `str`
989 Amplifier name to lookup linearity correction values.
990 """
991 frame = getDebugFrame(self._display, stepname)
992 if frame:
993 import matplotlib.pyplot as plt
994 fig, axs = plt.subplots(2)
996 if mask is None:
997 mask = np.ones_like(xVector, dtype=bool)
999 fig.suptitle(f"{stepname} {ampName} {self.config.linearityType}")
1000 if stepname == 'linearFit':
1001 axs[0].set_xlabel("Input Abscissa (time or mondiode)")
1002 axs[0].set_ylabel("Input Ordinate (flux)")
1003 axs[1].set_xlabel("Linear Ordinate (linear flux)")
1004 axs[1].set_ylabel("Flux Difference: (input - linear)")
1005 elif stepname in ('polyFit', 'splineFit'):
1006 axs[0].set_xlabel("Linear Abscissa (linear flux)")
1007 axs[0].set_ylabel("Input Ordinate (flux)")
1008 axs[1].set_xlabel("Linear Ordinate (linear flux)")
1009 axs[1].set_ylabel("Flux Difference: (input - full model fit)")
1010 elif stepname == 'solution':
1011 axs[0].set_xlabel("Input Abscissa (time or mondiode)")
1012 axs[0].set_ylabel("Linear Ordinate (linear flux)")
1013 axs[1].set_xlabel("Model flux (linear flux)")
1014 axs[1].set_ylabel("Flux Difference: (linear - model)")
1016 axs[0].set_yscale('log')
1017 axs[0].set_xscale('log')
1018 axs[0].scatter(xVector, yVector)
1019 axs[0].scatter(xVector[~mask], yVector[~mask], c='red', marker='x')
1020 axs[1].set_xscale('log')
1022 axs[1].scatter(yModel, yVector[mask] - yModel)
1023 fig.tight_layout()
1024 fig.show()
1026 prompt = "Press Enter or c to continue [chpx]..."
1027 while True:
1028 ans = input(prompt).lower()
1029 if ans in ("", " ", "c",):
1030 break
1031 elif ans in ("p", ):
1032 import pdb
1033 pdb.set_trace()
1034 elif ans in ("h", ):
1035 print("[h]elp [c]ontinue [p]db")
1036 elif ans in ('x', ):
1037 exit()
1038 plt.close()
1041class LinearityDoubleSplineSolveConnections(
1042 pipeBase.PipelineTaskConnections,
1043 dimensions=("instrument", "detector"),
1044):
1045 dummy = cT.Input(
1046 name="raw",
1047 doc="Dummy exposure.",
1048 storageClass='Exposure',
1049 dimensions=("instrument", "exposure", "detector"),
1050 multiple=True,
1051 deferLoad=True,
1052 )
1053 camera = cT.PrerequisiteInput(
1054 name="camera",
1055 doc="Camera Geometry definition.",
1056 storageClass="Camera",
1057 dimensions=("instrument", ),
1058 isCalibration=True,
1059 )
1060 inputLinearizerPtc = cT.Input(
1061 name="linearizerPtc",
1062 doc="Input linearizer PTC dataset.",
1063 storageClass="PhotonTransferCurveDataset",
1064 dimensions=("instrument", "detector"),
1065 isCalibration=True,
1066 )
1067 inputNormalization = cT.Input(
1068 name="cpLinearizerPtcNormalization",
1069 doc="Focal-plane normalization table.",
1070 storageClass="ArrowAstropy",
1071 dimensions=("instrument",),
1072 isCalibration=True,
1073 )
1074 inputBinnedImagesHandles = cT.Input(
1075 name="cpPtcPairBinned",
1076 doc="Tabulated binned exposure pairs.",
1077 storageClass="ArrowAstropy",
1078 dimensions=("instrument", "exposure", "detector"),
1079 multiple=True,
1080 deferLoad=True,
1081 )
1082 outputLinearizer = cT.Output(
1083 name="linearity",
1084 doc="Output linearity measurements.",
1085 storageClass="Linearizer",
1086 dimensions=("instrument", "detector"),
1087 isCalibration=True,
1088 )
1090 def __init__(self, *, config=None):
1091 super().__init__(config=config)
1093 if not config.useFocalPlaneNormalization:
1094 del self.inputNormalization
1097class LinearityDoubleSplineSolveConfig(
1098 pipeBase.PipelineTaskConfig,
1099 pipelineConnections=LinearityDoubleSplineSolveConnections,
1100):
1101 maxFracLinearityDeviation = pexConfig.Field(
1102 dtype=float,
1103 doc="Maximum fraction deviation from raw linearity to compute "
1104 "linearityTurnoff and linearityMaxSignal.",
1105 # TODO: DM-46811 investigate if this can be raised to 0.05.
1106 default=0.01,
1107 )
1108 minSignalFitLinearityTurnoff = pexConfig.Field(
1109 dtype=float,
1110 doc="Minimum signal to compute raw linearity slope for linearityTurnoff.",
1111 default=1000.0,
1112 )
1113 maxLinearityTurnoffRelativeToPtcTurnoff = pexConfig.Field(
1114 dtype=float,
1115 doc="Maximum fractional allowed linearity turnoff relative to the PTC turnoff. Used "
1116 "to keep extra-high odd values from contaminating the fit.",
1117 default=1.3,
1118 )
1119 maxNoiseReference = pexConfig.Field(
1120 dtype=float,
1121 doc="Maximum read noise (e-) in the PTC for an amp to be considered as a reference.",
1122 default=12.0,
1123 )
1124 usePhotodiode = pexConfig.Field(
1125 dtype=bool,
1126 doc="Use the photodiode info instead of the raw expTimes?",
1127 default=False,
1128 )
1129 minPhotodiodeCurrent = pexConfig.Field(
1130 dtype=float,
1131 doc="Minimum value to trust photodiode signals.",
1132 default=0.0,
1133 )
1134 doAutoGrouping = pexConfig.Field(
1135 dtype=bool,
1136 doc="Do automatic group detection? Cannot be True if splineGroupingColumn is also set. "
1137 "The automatic group detection will use the ratio of signal to exposure time (if "
1138 "autoGroupingUseExptime is True) or photodiode (if False) to determine which "
1139 "flat pairs were taken with different illumination settings.",
1140 default=False,
1141 )
1142 autoGroupingUseExptime = pexConfig.Field(
1143 dtype=bool,
1144 doc="Use exposure time to determine automatic grouping. Used if doAutoGrouping=True.",
1145 default=True,
1146 )
1147 autoGroupingThreshold = pexConfig.Field(
1148 dtype=float,
1149 doc="Minimum relative jump from sorted conversion values to determine a group.",
1150 default=0.1,
1151 )
1152 autoGroupingMaxSignalFraction = pexConfig.Field(
1153 dtype=float,
1154 doc="Only do auto-grouping when the signal is this fraction of the maximum signal. "
1155 "All exposures with signal higher than this threshold will be put into the "
1156 "largest signal group. This config is needed if the input PTC goes beyond "
1157 "the linearity turnoff.",
1158 default=0.9,
1159 )
1160 groupingColumn = pexConfig.Field(
1161 dtype=str,
1162 doc="Column to use for grouping together points, to allow "
1163 "for different proportionality constants. If None, then grouping will "
1164 "only be done if doAutoGrouping is True.",
1165 default=None,
1166 optional=True,
1167 )
1168 absoluteSplineMinimumSignalNode = pexConfig.Field(
1169 dtype=float,
1170 doc="Smallest node (above 0) for absolute spline (adu).",
1171 default=0.0,
1172 )
1173 absoluteSplineLowThreshold = pexConfig.Field(
1174 dtype=float,
1175 doc="Threshold for the low-level linearity nodes for absolute spline (adu). "
1176 "If this is below ``absoluteSplineMinimumSignalNode`` then the low "
1177 "level checks will be skipped.",
1178 default=0.0,
1179 )
1180 absoluteSplineLowNodeSize = pexConfig.Field(
1181 dtype=float,
1182 doc="Minimum size for low-level linearity nodes for absolute spline (adu).",
1183 default=2000.0,
1184 )
1185 absoluteSplineNodeSize = pexConfig.Field(
1186 dtype=float,
1187 doc="Minimum size for linearity nodes for absolute spline above absoluteSplineLowThreshold e(adu); "
1188 "note that there will always be a node at the reference PTC turnoff.",
1189 default=3000.0,
1190 )
1191 absoluteSplineFitMinIter = pexConfig.Field(
1192 dtype=int,
1193 doc="Minimum number of iterations for absolute spline fit.",
1194 default=3,
1195 )
1196 absoluteSplineFitMaxIter = pexConfig.Field(
1197 dtype=int,
1198 doc="Maximum number of iterations for absolute spline fit.",
1199 default=20,
1200 )
1201 absoluteSplineFitMaxRejectionPerIteration = pexConfig.Field(
1202 dtype=int,
1203 doc="Maximum number of rejections per iteration for absolute spline fit.",
1204 default=5,
1205 )
1206 absoluteNSigmaClipLinear = pexConfig.Field(
1207 dtype=float,
1208 doc="Sigma-clipping for absolute spline solution.",
1209 default=5.0,
1210 )
1211 doAbsoluteSplineFitOffset = pexConfig.Field(
1212 dtype=bool,
1213 doc="Fit a scattered light offset in the spline fit.",
1214 default=True,
1215 )
1216 doAbsoluteSplineFitWeights = pexConfig.Field(
1217 dtype=bool,
1218 doc="Fit linearity weight parameters in the spline fit.",
1219 default=False,
1220 )
1221 absoluteSplineFitWeightParsStart = pexConfig.ListField(
1222 dtype=float,
1223 doc="Starting parameters for weight fit, if doSplineFitWeights=True. "
1224 "Parameters are such that sigma = sqrt(par[0]**2. + par[1]**2./mu)."
1225 "If doSplineFitWeights=False then these are used as-is; otherwise "
1226 "they are used as the initial values for fitting these parameters.",
1227 length=2,
1228 default=[1.0, 0.0],
1229 )
1230 doAbsoluteSplineFitTemperature = pexConfig.Field(
1231 dtype=bool,
1232 doc="Fit temperature coefficient in spline fit?",
1233 default=False,
1234 )
1235 absoluteSplineFitTemperatureColumn = pexConfig.Field(
1236 dtype=str,
1237 doc="Name of the temperature column to use when fitting temperature "
1238 "coefficients in spline fit; this must not be None if "
1239 "doSplineFitTemperature is True.",
1240 default=None,
1241 optional=True,
1242 )
1243 doAbsoluteSplineFitTemporal = pexConfig.Field(
1244 dtype=bool,
1245 doc="Fit a linear temporal parameter coefficient in spline fit?",
1246 default=False,
1247 )
1248 useFocalPlaneNormalization = pexConfig.Field(
1249 dtype=bool,
1250 doc="Use focal-plane normalization in addition to/instead of photodiode? "
1251 "(Only used with for absolute spline fitting).",
1252 default=False,
1253 )
1254 relativeSplineReferenceCounts = pexConfig.Field(
1255 dtype=float,
1256 doc="Number of target counts (adu) to select a reference image for "
1257 "relative spline solution.",
1258 default=10000.0,
1259 )
1260 relativeSplineMinimumSignalNode = pexConfig.Field(
1261 dtype=float,
1262 doc="Smallest node (above 0) for relative spline (adu).",
1263 default=100.0,
1264 )
1265 relativeSplineLowThreshold = pexConfig.Field(
1266 dtype=float,
1267 doc="Threshold for the low-level linearity nodes for relative spline (adu)."
1268 "If this is below ``relativeSplineMinimumSignalNode`` then the low "
1269 "level checks will be skipped.",
1270 default=5000.0,
1271 )
1272 relativeSplineLowNodeSize = pexConfig.Field(
1273 dtype=float,
1274 doc="Minimum size for low-level linearity nodes for relative spline (adu).",
1275 default=750.0,
1276 )
1277 relativeSplineMidNodeSize = pexConfig.Field(
1278 dtype=float,
1279 doc="Minimum size for mid-level linearity nodes for relative spline (adu); "
1280 "this applies to counts between relativeSplineLowThreshold and the "
1281 "PTC turnoff.",
1282 default=5000.0,
1283 )
1284 relativeSplineHighNodeSize = pexConfig.Field(
1285 dtype=float,
1286 doc="Minimum size for high-level linearity nodes for relative spline (adu); "
1287 "this applies to counts between the PTC and linearity turnoffs.",
1288 default=2000.0,
1289 )
1290 relativeSplineFitMinIter = pexConfig.Field(
1291 dtype=int,
1292 doc="Minimum number of iterations for relative spline fit.",
1293 default=3,
1294 )
1295 relativeSplineFitMaxIter = pexConfig.Field(
1296 dtype=int,
1297 doc="Maximum number of iterations for relative spline fit.",
1298 default=20,
1299 )
1300 relativeSplineFitMaxRejectionPerIteration = pexConfig.Field(
1301 dtype=int,
1302 doc="Maximum number of rejections per iteration for relative spline fit.",
1303 default=5,
1304 )
1305 relativeNSigmaClipLinear = pexConfig.Field(
1306 dtype=float,
1307 doc="Sigma-clipping for relative spline solution.",
1308 default=5.0,
1309 )
1311 def validate(self):
1312 super().validate()
1314 if self.doAbsoluteSplineFitTemperature and self.absoluteSplineFitTemperatureColumn is None:
1315 raise ValueError(
1316 "Must set absoluteSplineFitTemperatureColumn if doAbsoluteSplineFitTemperature is True.",
1317 )
1319 if self.doAutoGrouping and self.groupingColumn is not None:
1320 raise ValueError("Must not set doAutoGrouping and also groupingColumn")
1321 if self.doAutoGrouping:
1322 if not self.autoGroupingUseExptime and not self.usePhotodiode:
1323 raise ValueError("If doAutoGrouping is True and autoGroupingUseExptime is False, then "
1324 "usePhotodiode must be True.")
1327class LinearityDoubleSplineSolveTask(pipeBase.PipelineTask):
1328 ConfigClass = LinearityDoubleSplineSolveConfig
1329 _DefaultName = "cpLinearityDoubleSplineSolve"
1331 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1332 # docstring inherited
1333 inputs = butlerQC.get(inputRefs)
1335 if self.config.useFocalPlaneNormalization:
1336 inputNormalization = inputs["inputNormalization"]
1337 else:
1338 inputNormalization = None
1340 # Add calibration provenance info to header.
1341 kwargs = dict()
1342 reference = getattr(inputRefs, "inputLinearizerPtc", None)
1344 if reference is not None and hasattr(reference, "run"):
1345 runKey = "PTC_RUN"
1346 runValue = reference.run
1347 idKey = "PTC_UUID"
1348 idValue = str(reference.id)
1349 dateKey = "PTC_DATE"
1350 calib = inputs.get("inputPtc", None)
1351 dateValue = extractCalibDate(calib)
1353 kwargs[runKey] = runValue
1354 kwargs[idKey] = idValue
1355 kwargs[dateKey] = dateValue
1357 self.log.info("Using " + str(reference.run))
1359 outputs = self.run(
1360 inputPtc=inputs["inputLinearizerPtc"],
1361 camera=inputs["camera"],
1362 inputBinnedImagesHandles=inputs["inputBinnedImagesHandles"],
1363 inputNormalization=inputNormalization,
1364 )
1365 outputs.outputLinearizer.updateMetadata(setDate=False, **kwargs)
1367 butlerQC.put(outputs, outputRefs)
1369 def run(
1370 self,
1371 *,
1372 inputPtc,
1373 camera,
1374 inputBinnedImagesHandles,
1375 inputNormalization,
1376 ):
1377 """Fit the double-spline relative/absolute linearity correction.
1379 Parameters
1380 ----------
1381 inputPtc : `lsst.ip.isr.PtcDataset`
1382 Pre-measured PTC dataset.
1383 camera : `lsst.afw.cameraGeom.Camera`
1384 Camera geometry.
1385 inputBinnedImagesHandles : `list` [`DeferredDatasetHandle`]
1386 Handles for input binned image pairs.
1387 inputNormalization : `astropy.table.Table`, optional
1388 Focal plane normalization table to use if
1389 useFocalPlaneNormalization is True.
1391 Returns
1392 -------
1393 results : `lsst.pipe.base.Struct`
1394 The results struct containing:
1396 ``outputLinearizer``
1397 Final linearizer calibration (`lsst.ip.isr.Linearizer`).
1398 ``outputProvenance``
1399 Provenance data for the new calibration
1400 (`lsst.ip.isr.IsrProvenance`).
1401 """
1402 detector = camera[inputPtc.metadata["DETECTOR"]]
1404 binnedImagesHandleDict = {
1405 handle.dataId["exposure"]: handle for handle in inputBinnedImagesHandles
1406 }
1408 linearizer = Linearizer(detector=detector, log=self.log)
1409 linearizer.updateMetadataFromExposures([inputPtc])
1411 groupingValues = _determineInputGroups(
1412 inputPtc,
1413 self.config.doAutoGrouping,
1414 self.config.autoGroupingUseExptime,
1415 self.config.autoGroupingMaxSignalFraction,
1416 self.config.autoGroupingThreshold,
1417 self.config.groupingColumn,
1418 self.config.minPhotodiodeCurrent,
1419 )
1421 if self.config.doAbsoluteSplineFitTemperature:
1422 if self.config.absoluteSplineFitTemperatureColumn not in inputPtc.auxValues:
1423 raise ValueError(
1424 "Config requests fitting temperature coefficient for "
1425 f"{self.config.splineFitTemperatureColumn} but this column "
1426 "is not available in inputPtc.auxValues.",
1427 )
1428 temperatureValues = inputPtc.auxValues[self.config.splineFitTemperatureColumn]
1429 else:
1430 temperatureValues = None
1432 # Fill the linearizer with empty values.
1433 firstAmp = None
1434 for ampName in inputPtc.ampNames:
1435 if ampName not in inputPtc.badAmps:
1436 firstAmp = ampName
1437 break
1438 if firstAmp is None:
1439 raise pipeBase.NoWorkFound("No valid amps in input PTC.")
1440 nExp = len(inputPtc.inputExpIdPairs[firstAmp]) * 2
1441 nAmp = len(inputPtc.ampNames)
1443 for amp in detector:
1444 ampName = amp.getName()
1446 linearizer.inputGain[ampName] = inputPtc.gain[ampName]
1447 linearizer.linearityType[ampName] = "None"
1448 linearizer.linearityCoeffs[ampName] = np.zeros(1)
1449 # This is not used; kept for compatibility.
1450 linearizer.linearityBBox[ampName] = amp.getBBox()
1451 linearizer.fitParams[ampName] = np.zeros(1)
1452 linearizer.fitParamsErr[ampName] = np.zeros(1)
1453 linearizer.fitChiSq[ampName] = np.nan
1454 linearizer.fitResiduals[ampName] = np.zeros(nExp)
1455 linearizer.fitResidualsSigmaMad[ampName] = np.nan
1456 linearizer.fitResidualsUnmasked[ampName] = np.zeros(nExp)
1457 linearizer.fitResidualsModel[ampName] = np.zeros(nExp)
1458 linearizer.linearFit[ampName] = np.zeros(2)
1459 linearizer.linearityTurnoff[ampName] = np.nan
1460 linearizer.linearityMaxSignal[ampName] = np.nan
1461 linearizer.inputMask[ampName] = np.zeros(nExp, dtype=np.bool_)
1462 linearizer.inputAbscissa[ampName] = np.zeros(nExp)
1463 linearizer.inputOrdinate[ampName] = np.zeros(nExp)
1464 linearizer.inputGroupingIndex[ampName] = np.zeros(nExp, dtype=np.int64)
1465 linearizer.inputNormalization[ampName] = np.ones(nExp)
1467 linearizer.absoluteReferenceAmplifier = ""
1469 # Extract values in common, and per-amp.
1470 data = np.zeros(
1471 nExp,
1472 dtype=[
1473 ("exp_id", "i8"),
1474 ("exptime", "f8"),
1475 ("photocharge", "f8"),
1476 ("mjd", "f8"),
1477 ("raw_mean", ("f8", nAmp)),
1478 ("abscissa", "f8"),
1479 ("grouping", "i4"),
1480 # The following are computed in the relative scaling
1481 # measurements.
1482 ("ref_counts", "f8"),
1483 ("gain_ratio", ("f8", nAmp)),
1484 ],
1485 )
1487 data["exp_id"] = np.asarray(inputPtc.inputExpIdPairs[firstAmp]).ravel()
1488 data["exptime"] = np.repeat(inputPtc.rawExpTimes[firstAmp], 2)
1489 data["mjd"] = np.repeat(inputPtc.inputExpPairMjdStartList[firstAmp], 2)
1490 data["photocharge"] = np.repeat(inputPtc.photoCharges[firstAmp], 2)
1491 data["photocharge"][::2] -= inputPtc.photoChargeDeltas[firstAmp] / 2.
1492 data["photocharge"][1::2] += inputPtc.photoChargeDeltas[firstAmp] / 2.
1493 data["grouping"] = np.repeat(groupingValues, 2)
1495 for i, amp in enumerate(detector):
1496 ampName = amp.getName()
1498 data["raw_mean"][:, i] = np.repeat(inputPtc.rawMeans[ampName], 2)
1499 data["raw_mean"][::2, i] -= inputPtc.rawDeltas[ampName] / 2.
1500 data["raw_mean"][1::2, i] += inputPtc.rawDeltas[ampName] / 2.
1502 if self.config.usePhotodiode:
1503 data["abscissa"][:] = data["photocharge"]
1505 data["abscissa"][data["photocharge"] < self.config.minPhotodiodeCurrent] = np.nan
1506 else:
1507 data["abscissa"][:] = data["exptime"]
1509 # Normalize if configured.
1510 inputNorm = np.ones(nExp, dtype=np.float64)
1511 if self.config.useFocalPlaneNormalization:
1512 a, b = esutil.numpy_util.match(data["exp_id"], inputNormalization["exposure"])
1513 inputNorm[a] = inputNormalization["normalization"][b]
1514 data["abscissa"] *= inputNorm
1516 postTurnoffMasks = {}
1518 # Compute linearity turnoff for each amp.
1519 for i, amp in enumerate(detector):
1520 ampName = amp.getName()
1522 if ampName in inputPtc.badAmps:
1523 self.log.warning(
1524 "Amp %s in detector %s has no usable PTC information. Skipping!",
1525 ampName,
1526 detector.getName(),
1527 )
1528 continue
1530 mask = np.isfinite(data["raw_mean"][:, i])
1532 turnoffMask = np.repeat(inputPtc.expIdMask[ampName], 2)
1533 turnoffMask &= mask
1535 _, turnoff, maxSignal, goodPoints = _computeTurnoffAndMax(
1536 data["abscissa"],
1537 data["raw_mean"][:, i],
1538 turnoffMask,
1539 data["grouping"],
1540 ampName=ampName,
1541 minSignalFitLinearityTurnoff=self.config.minSignalFitLinearityTurnoff,
1542 maxFracLinearityDeviation=self.config.maxFracLinearityDeviation,
1543 log=self.log,
1544 maxTurnoff=inputPtc.ptcTurnoff[ampName] * self.config.maxLinearityTurnoffRelativeToPtcTurnoff,
1545 )
1547 # Use the goodPoints as an initial estimate of the mask
1548 # above the turnoff. But we only want to maintain the
1549 # "high end" outliers.
1550 postTurnoffMask = goodPoints
1551 postTurnoffMask[data["raw_mean"][:, i] < np.median(data["raw_mean"][goodPoints, i])] = True
1552 postTurnoffMasks[ampName] = postTurnoffMask
1554 if np.isnan(turnoff):
1555 # This is a bad amp, with no linearizer.
1556 self.log.warning(
1557 "Amp %s in detector %s has no usable linearizer information. Skipping!",
1558 ampName,
1559 detector.getName(),
1560 )
1561 continue
1563 linearizer.linearityTurnoff[ampName] = turnoff
1564 linearizer.linearityMaxSignal[ampName] = maxSignal
1566 self.log.info("Amplifier %s has a linearity turnoff of %.2f adu.", ampName, turnoff)
1568 # Choose the reference amplifier as the one with the largest
1569 # turnoff. This ensures that the absolute fit covers the full
1570 # range. We additionally confirm that the ptc turnoff is
1571 # finite for this amplifier.
1572 turnoffArray = np.asarray([linearizer.linearityTurnoff[ampName] for ampName in inputPtc.ampNames])
1573 # This is a possibly redundant check to make sure that a bad amp is
1574 # not chosen as a reference amp. We also check that a high noise
1575 # amp is not chosen as the reference amp.
1576 for i, ampName in enumerate(inputPtc.ampNames):
1577 if ampName in inputPtc.badAmps \
1578 or not np.isfinite(inputPtc.ptcTurnoff[ampName]) \
1579 or inputPtc.noise[ampName] > self.config.maxNoiseReference:
1580 turnoffArray[i] = np.nan
1582 if np.all(~np.isfinite(turnoffArray)):
1583 # Return the default blank linearizer.
1584 linearizer.hasLinearity = True
1585 linearizer.validate()
1586 linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
1587 linearizer.updateMetadata(setDate=True, setCalibId=True)
1588 linearizer.updateMetadataFromExposures([inputPtc])
1589 provenance = IsrProvenance(calibType='linearizer')
1591 return pipeBase.Struct(
1592 outputLinearizer=linearizer,
1593 outputProvenance=provenance,
1594 )
1596 refAmpIndex = np.argmax(np.nan_to_num(turnoffArray))
1597 refAmpName = inputPtc.ampNames[refAmpIndex]
1598 linearizer.absoluteReferenceAmplifier = refAmpName
1600 # Choose a reference image.
1601 refExpIndex = np.argmin(
1602 np.abs(
1603 np.nan_to_num(data["raw_mean"][:, refAmpIndex]) - self.config.relativeSplineReferenceCounts
1604 )
1605 )
1606 refExpId = data["exp_id"][refExpIndex]
1608 self.log.info(
1609 "Using exposure %d (%.2f adu in amp %s) as reference.",
1610 refExpId,
1611 data["raw_mean"][refExpIndex, refAmpIndex],
1612 refAmpName,
1613 )
1615 # We need to know if the reference exposure is the first or second
1616 # in the pair, because the binned are pairs.
1617 offset = refExpIndex % 2
1619 refBinned = binnedImagesHandleDict[data["exp_id"][refExpIndex - offset]].get()
1620 refBinned = copy.copy(refBinned)
1621 if offset == 0:
1622 refBinned["value"] = refBinned["value1"]
1623 else:
1624 refBinned["value"] = refBinned["value2"]
1626 # Scale reference according to reference amplifier.
1627 refScaling = np.median(refBinned["value"][refBinned["amp_index"] == refAmpIndex])
1628 refBinned["value"] /= refScaling
1630 # Get the invidual amp scalings.
1631 # These are the relative gains for the reference image.
1632 ampScalings = np.asarray(
1633 [
1634 np.median(refBinned["value"][refBinned["amp_index"] == ampIndex])
1635 for ampIndex in range(nAmp)
1636 ],
1637 )
1639 # Compute the gain ratios for every exposure.
1640 # The binned images are stored as pairs.
1641 self.log.info("Computing gain ratios for %d exposures.", len(data))
1642 for i in range(len(data)):
1643 expId = data["exp_id"][i]
1644 if (i % 2) == 0:
1645 binned = binnedImagesHandleDict[expId].get()
1646 binned["value"] = binned["value1"]
1647 else:
1648 binned["value"] = binned["value2"]
1650 binned["value"] /= refBinned["value"]
1652 gainRatios = np.asarray(
1653 [
1654 np.median(binned["value"][binned["amp_index"] == ampIndex])
1655 for ampIndex in range(nAmp)
1656 ]
1657 )
1658 ref_counts = gainRatios[refAmpIndex]
1659 gainRatios /= ref_counts
1661 data["ref_counts"][i] = ref_counts
1662 data["gain_ratio"][i, :] = gainRatios
1664 # We need to know which group has the largest size.
1665 groupAmplitudes = np.zeros(len(np.unique(data["grouping"])))
1666 for g in range(len(groupAmplitudes)):
1667 use = (data["grouping"] == g)
1668 groupAmplitudes[g] = np.nanmax(data["ref_counts"][use]) - np.nanmin(data["ref_counts"][use])
1669 maxAmplitudeGroup = np.argmax(groupAmplitudes)
1671 self.log.info("Illumination group %d has the largest signal amplitude.", maxAmplitudeGroup)
1673 # Compute relative linearization first.
1674 maxRelNodes = 0
1676 for i, amp in enumerate(detector):
1677 if i == refAmpIndex:
1678 continue
1680 ampName = amp.getName()
1682 ptcTurnoff = inputPtc.ptcTurnoff[ampName]
1683 linearityTurnoff = linearizer.linearityTurnoff[ampName]
1685 if not np.isfinite(ptcTurnoff) or not np.isfinite(linearityTurnoff):
1686 # This is a bad amp; skip it.
1687 continue
1689 if ptcTurnoff < self.config.relativeSplineLowThreshold:
1690 lowThreshold = 0.0
1691 else:
1692 lowThreshold = self.config.relativeSplineLowThreshold
1694 relAbscissa = data["ref_counts"] * ampScalings[i]
1695 relOrdinate = data["ref_counts"] * data["gain_ratio"][:, i] * ampScalings[i]
1697 # The mask here must exclude everything beyond the turnoff.
1698 # Note that we need to do this before we use the actual
1699 # turnoff to compute the nodes to avoid nodes going past the
1700 # data domain.
1701 relMask = (
1702 postTurnoffMasks[ampName]
1703 & np.isfinite(relAbscissa)
1704 & np.isfinite(relOrdinate)
1705 & (relOrdinate < linearityTurnoff)
1706 )
1708 # Make sure that the linearity turnoff used here does not
1709 # go beyond the max value of the relOrdinate
1710 relTurnoff = min(linearityTurnoff, np.max(relOrdinate[relMask]))
1712 relNodes = _noderator(
1713 lowThreshold,
1714 ptcTurnoff,
1715 relTurnoff,
1716 self.config.relativeSplineMinimumSignalNode,
1717 self.config.relativeSplineLowNodeSize,
1718 self.config.relativeSplineMidNodeSize,
1719 self.config.relativeSplineHighNodeSize,
1720 )
1722 self.log.info(
1723 "Relative linearity for amplifier %s using %d nodes from %.2f to %.2f counts.",
1724 ampName,
1725 len(relNodes),
1726 relNodes[0],
1727 relNodes[-1],
1728 )
1730 # Update the number of relative nodes to concatenation.
1731 if len(relNodes) > maxRelNodes:
1732 maxRelNodes = len(relNodes)
1734 linearizer.inputMask[ampName] = relMask.copy()
1735 linearizer.inputAbscissa[ampName] = relAbscissa.copy()
1736 linearizer.inputOrdinate[ampName] = relOrdinate.copy()
1737 linearizer.inputGroupingIndex[ampName] = data["grouping"].copy()
1738 linearizer.inputNormalization[ampName] = np.ones_like(relAbscissa)
1740 fitter = AstierSplineLinearityFitter(
1741 relNodes,
1742 data["grouping"],
1743 relAbscissa,
1744 relOrdinate,
1745 mask=relMask,
1746 fit_offset=False,
1747 fit_weights=False,
1748 fit_temperature=False,
1749 max_signal_nearly_linear=ptcTurnoff,
1750 fit_temporal=False,
1751 # Put a cap on the maximum correction in absolute value.
1752 max_frac_correction=np.inf,
1753 max_correction=10_000.0,
1754 )
1755 p0 = fitter.estimate_p0()
1756 pars = fitter.fit(
1757 p0,
1758 min_iter=self.config.relativeSplineFitMinIter,
1759 max_iter=self.config.relativeSplineFitMaxIter,
1760 max_rejection_per_iteration=self.config.relativeSplineFitMaxRejectionPerIteration,
1761 n_sigma_clip=self.config.relativeNSigmaClipLinear,
1762 )
1764 # Confirm that the first parameter is 0, and set it to
1765 # exactly zero.
1766 relValues = pars[fitter.par_indices["values"]]
1767 if not np.isclose(relValues[0], 0):
1768 raise RuntimeError("Programmer error! First spline parameter must "
1769 "be consistent with zero.")
1770 relValues[0] = 0.0
1772 if np.any(np.abs(pars[fitter.par_indices["values"]]) > 10_000.0):
1773 self.log.warning("Unconstrained nodes escaped containment; clipping.")
1774 lo = (pars[fitter.par_indices["values"]] < -10_000.0)
1775 if np.sum(lo) > 0:
1776 pars[fitter.par_indices["values"][lo]] = -10_000.0
1777 hi = (pars[fitter.par_indices["values"]] > 10_000.0)
1778 if np.sum(hi) > 0:
1779 pars[fitter.par_indices["values"][hi]] = 10_000.0
1781 # We adjust the node values according to the slope of the
1782 # group with the largest amplitude. This removes a degeneracy
1783 # in the normalization and ensures that the overall linearized
1784 # correction is as close to the reference as possible.
1785 relValues -= (1.0 - pars[fitter.par_indices["groups"][maxAmplitudeGroup]]) * relNodes
1787 relChisq = fitter.compute_chisq_dof(pars)
1789 # Our reference fit is always 1.0 slope.
1790 relLinearFit = np.array([0.0, 1.0])
1792 # Adjust the abscissa for different groups for residuals.
1793 for j, groupIndex in enumerate(fitter.group_indices):
1794 relAbscissa[groupIndex] *= (pars[fitter.par_indices["groups"][j]] / relLinearFit[1])
1796 relMask = fitter.mask
1798 # Record values in the linearizer.
1799 linearizer.linearityType[ampName] = "DoubleSpline"
1800 # Note that we have a placeholder for the number of nodes in
1801 # the absolute spline.
1802 linearizer.linearityCoeffs[ampName] = np.concatenate([[len(relNodes), 0], relNodes, relValues])
1803 linearizer.fitChiSq[ampName] = relChisq
1804 linearizer.linearFit[ampName] = relLinearFit
1806 # Compute residuals.
1807 spl = Akima1DInterpolator(relNodes, relValues, method="akima")
1808 relOffset = spl(np.clip(relOrdinate, relNodes[0], relNodes[-1]))
1809 relModel = relOrdinate - relOffset
1811 if relMask.sum() < 2:
1812 self.log.warning("Amp %s in detector %s has not enough points in linear ordinate "
1813 "for residuals. Skipping!", ampName, detector.getName())
1814 relResiduals = np.full_like(relModel, np.nan)
1815 relResidualsUnmasked = relResiduals.copy()
1816 else:
1817 postLinearFit, _, _, _ = irlsFit(
1818 relLinearFit,
1819 relAbscissa[relMask],
1820 relModel[relMask],
1821 funcPolynomial,
1822 )
1823 # When computing residuals, we only care about the slope of
1824 # the postLinearFit and not the intercept. The intercept
1825 # itself depends on a possibly unknown zero in the abscissa
1826 # (often photodiode) which may have an arbitrary value.
1827 relResiduals = relModel - (postLinearFit[1] * relAbscissa)
1828 relResidualsUnmasked = relResiduals.copy()
1829 # We set masked residuals to nan.
1830 relResiduals[~relMask] = np.nan
1832 linearizer.fitResidualsUnmasked[ampName] = relResidualsUnmasked
1833 linearizer.fitResiduals[ampName] = relResiduals
1834 linearizer.fitResidualsModel[ampName] = relModel.copy()
1836 finite = np.isfinite(relResiduals)
1837 if finite.sum() == 0:
1838 sigmad = np.nan
1839 else:
1840 sigmad = median_abs_deviation(relResiduals[finite]/relOrdinate[finite], scale="normal")
1841 linearizer.fitResidualsSigmaMad[ampName] = sigmad
1843 # Now compute absolute linearization.
1845 if temperatureValues is not None:
1846 temperatureValuesScaled = temperatureValues - np.median(temperatureValues)
1847 else:
1848 temperatureValuesScaled = None
1850 if self.config.doAbsoluteSplineFitTemporal:
1851 inputMjdScaled = data["mjd"].copy()
1852 inputMjdScaled -= np.nanmedian(inputMjdScaled)
1853 else:
1854 inputMjdScaled = None
1856 absAbscissa = data["abscissa"].copy()
1857 absOrdinate = data["ref_counts"].copy()
1859 # These are guaranteed to be finite (as checked previously).
1860 absPtcTurnoff = inputPtc.ptcTurnoff[refAmpName]
1861 absLinearityTurnoff = linearizer.linearityTurnoff[refAmpName]
1863 if absPtcTurnoff < self.config.absoluteSplineLowThreshold:
1864 lowThreshold = 0.0
1865 else:
1866 lowThreshold = self.config.absoluteSplineLowThreshold
1868 # The mask here must exclude everything beyond the turnoff.
1869 # Note that we need to do this before we use the actual
1870 # turnoff to compute the nodes to avoid nodes going past the
1871 # data domain.
1872 absMask = postTurnoffMasks[refAmpName] & np.isfinite(absAbscissa) & np.isfinite(absOrdinate)
1874 absLinearityTurnoff = min(absLinearityTurnoff, np.max(absOrdinate[absMask]))
1876 absNodes = _noderator(
1877 lowThreshold,
1878 absPtcTurnoff,
1879 absLinearityTurnoff,
1880 self.config.absoluteSplineMinimumSignalNode,
1881 self.config.absoluteSplineLowNodeSize,
1882 # The medium and high are matched for absolute spline.
1883 self.config.absoluteSplineNodeSize,
1884 self.config.absoluteSplineNodeSize,
1885 )
1887 self.log.info("Absolute linearity for using %d nodes.", len(absNodes))
1889 # We store the absolute residuals with the reference amplifier.
1890 linearizer.linearityType[refAmpName] = "DoubleSpline"
1891 linearizer.inputMask[refAmpName] = absMask.copy()
1892 linearizer.inputAbscissa[refAmpName] = absAbscissa.copy()
1893 linearizer.inputOrdinate[refAmpName] = absOrdinate.copy()
1894 linearizer.inputGroupingIndex[refAmpName] = data["grouping"].copy()
1895 linearizer.inputNormalization[refAmpName] = inputNorm.copy()
1897 fitter = AstierSplineLinearityFitter(
1898 absNodes,
1899 data["grouping"].copy(),
1900 absAbscissa,
1901 absOrdinate,
1902 mask=absMask,
1903 log=self.log,
1904 fit_offset=self.config.doAbsoluteSplineFitOffset,
1905 fit_weights=self.config.doAbsoluteSplineFitWeights,
1906 weight_pars_start=self.config.absoluteSplineFitWeightParsStart,
1907 fit_temperature=self.config.doAbsoluteSplineFitTemperature,
1908 temperature_scaled=temperatureValuesScaled,
1909 max_signal_nearly_linear=absPtcTurnoff,
1910 fit_temporal=self.config.doAbsoluteSplineFitTemporal,
1911 mjd_scaled=inputMjdScaled,
1912 )
1913 p0 = fitter.estimate_p0(use_all_for_normalization=True)
1914 pars = fitter.fit(
1915 p0,
1916 min_iter=self.config.absoluteSplineFitMinIter,
1917 max_iter=self.config.absoluteSplineFitMaxIter,
1918 max_rejection_per_iteration=self.config.absoluteSplineFitMaxRejectionPerIteration,
1919 n_sigma_clip=self.config.absoluteNSigmaClipLinear,
1920 )
1922 # Confirm that the first parameter is 0, and set it to
1923 # exactly zero.
1924 absValues = pars[fitter.par_indices["values"]]
1925 if not np.isclose(absValues[0], 0):
1926 raise RuntimeError("Programmer error! First spline parameter must "
1927 "be consistent with zero.")
1928 absValues[0] = 0.0
1930 # We need a place to store this.
1931 linearizer.fitChiSq[refAmpName] = fitter.compute_chisq_dof(pars)
1933 absLinearFit = np.array([0.0, np.mean(pars[fitter.par_indices["groups"]])])
1935 # We must modify the inputOrdinate according to the
1936 # nuisance terms in the linearity fit for the residual
1937 # computation code to work properly.
1938 # The true mu (inputOrdinate) is given by
1939 # mu = mu_in * (1 + alpha*t_scale) * (1 + beta*mjd_scale)
1940 if self.config.doAbsoluteSplineFitTemperature:
1941 absOrdinate *= (1.0
1942 + pars[fitter.par_indices["temperature_coeff"]]*temperatureValuesScaled)
1943 if self.config.doAbsoluteSplineFitTemporal:
1944 absOrdinate *= (1.0
1945 + pars[fitter.par_indices["temporal_coeff"]]*inputMjdScaled)
1947 # Adjust the abscissa for different groups for residuals.
1948 for j, groupIndex in enumerate(fitter.group_indices):
1949 absAbscissa[groupIndex] *= (pars[fitter.par_indices["groups"][j]] / absLinearFit[1])
1950 # And remove the offset term.
1951 if self.config.doAbsoluteSplineFitOffset:
1952 absOrdinate -= pars[fitter.par_indices["offset"]]
1954 absMask = fitter.mask
1956 # Compute residuals.
1957 spl = Akima1DInterpolator(absNodes, absValues, method="akima")
1958 absOffset = spl(np.clip(absOrdinate, absNodes[0], absNodes[-1]))
1959 absModel = absOrdinate - absOffset
1961 if absMask.sum() < 2:
1962 self.log.warning("Detector %s has not enough points in linear ordinate "
1963 "for residuals. Skipping!", detector.getName())
1964 # We have to KICK OUT HERE something is VERY wrong.
1965 absResiduals = np.full_like(absModel, np.nan)
1966 absResidualsUnmasked = relResiduals.copy()
1967 else:
1968 postLinearFit, _, _, _ = irlsFit(
1969 absLinearFit,
1970 absAbscissa[absMask],
1971 absModel[absMask],
1972 funcPolynomial,
1973 )
1974 # When computing residuals, we only care about the slope of
1975 # the postLinearFit and not the intercept. The intercept
1976 # itself depends on a possibly unknown zero in the abscissa
1977 # (often photodiode) which may have an arbitrary value.
1978 absResiduals = absModel - (postLinearFit[1] * absAbscissa)
1979 absResidualsUnmasked = absResiduals.copy()
1980 # We set masked residuals to nan.
1981 absResiduals[~absMask] = np.nan
1983 linearizer.fitResidualsUnmasked[refAmpName] = absResidualsUnmasked
1984 linearizer.fitResiduals[refAmpName] = absResiduals
1985 linearizer.fitResidualsModel[refAmpName] = absModel.copy()
1987 finite = np.isfinite(absResiduals)
1988 if finite.sum() == 0:
1989 sigmad = np.nan
1990 else:
1991 sigmad = median_abs_deviation(absResiduals[finite]/absOrdinate[finite], scale="normal")
1992 linearizer.fitResidualsSigmaMad[refAmpName] = sigmad
1994 # Record the absolute nodes and values in each individual amplifier,
1995 # along with extra padding for alignment.
1996 nAbsNodes = len(absNodes)
1997 for i, amp in enumerate(detector):
1998 ampName = amp.getName()
2000 coeffs = np.zeros(2 * nAbsNodes + 2 * maxRelNodes + 2)
2001 if ampName == refAmpName:
2002 # The reference amplifier only has the absolute spline.
2003 coeffs[1] = nAbsNodes
2004 coeffs[2: 2 + 2 * nAbsNodes] = np.concatenate([absNodes, absValues])
2005 else:
2006 nRelNodes = int(linearizer.linearityCoeffs[ampName][0])
2008 coeffs = np.zeros(2 * nAbsNodes + 2 * maxRelNodes + 2)
2009 coeffs[0] = nRelNodes
2010 coeffs[1] = nAbsNodes
2011 relStart = 2
2012 relEnd = relStart + 2 * nRelNodes
2013 coeffs[relStart: relEnd] = linearizer.linearityCoeffs[ampName][relStart: relEnd]
2014 absStart = relEnd
2015 absEnd = absStart + 2 * nAbsNodes
2016 coeffs[absStart: absEnd] = np.concatenate([absNodes, absValues])
2018 linearizer.linearityCoeffs[ampName] = coeffs
2020 linearizer.hasLinearity = True
2021 linearizer.validate()
2022 linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
2023 linearizer.updateMetadata(setDate=True, setCalibId=True)
2024 linearizer.updateMetadataFromExposures([inputPtc])
2025 provenance = IsrProvenance(calibType='linearizer')
2027 return pipeBase.Struct(
2028 outputLinearizer=linearizer,
2029 outputProvenance=provenance,
2030 )
2033def _determineInputGroups(
2034 ptc,
2035 doAutoGrouping,
2036 autoGroupingUseExptime,
2037 autoGroupingMaxSignalFraction,
2038 autoGroupingThreshold,
2039 groupingColumn,
2040 minPhotodiodeCurrent,
2041):
2042 """Determine input groups for linearity fit.
2044 If ``splineGroupingColumn`` is set, then grouping will be done
2045 based on this. Otherwise, if ``doAutoGrouping`` is False, then
2046 no grouping will be done. Finally, grouping will be done by measuring
2047 the ratio of signal to exposure time (if
2048 ``autoGroupingUseExptime`` is set; recommended) or photocharge.
2049 These are then clustered with a simple algorithm to split into groups.
2050 If the data was taking by varying exposure time at different
2051 illumination levels, this grouping is very robust as the clusters are
2052 very well separated.
2054 Parameters
2055 ----------
2056 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
2057 Input PTC to do grouping.
2058 doAutoGrouping : `bool`
2059 Do automatic grouping of pairs?
2060 autoGroupingUseExptime : `bool`
2061 Use exposure time for automatic grouping of pairs?
2062 autoGroupingMaxSignalFraction : `float`
2063 All exposures with signal higher than this threshold will
2064 be put into the largest signal group.
2065 autoGroupingThreshold : `float`
2066 Minimum relative jump from sorted values to determine a group.
2067 minPhotodiodeCurrent : `float`
2068 Minimum photodiode current if auto-grouping is used and
2069 autoGroupingUseExptime is False.
2070 splineGroupingColumn : `str` or `None`
2071 Column to be used for spline grouping (if doAutoGrouping is False).
2073 Returns
2074 -------
2075 groupingValues : `np.ndarray`
2076 Array of values that are unique for a given group.
2077 """
2078 nPair = np.asarray(ptc.inputExpIdPairs[ptc.ampNames[0]]).shape[0]
2079 groupingValues = np.zeros(nPair, dtype=np.int64)
2081 if not doAutoGrouping:
2082 if groupingColumn is not None:
2083 if groupingColumn not in ptc.auxValues:
2084 raise ValueError(f"Config requests grouping by {groupingColumn}, "
2085 "but this column is not available in ptc.auxValues.")
2087 uGroupValues = np.unique(ptc.auxValues[groupingColumn])
2088 for i, uGroupValue in enumerate(uGroupValues):
2089 groupingValues[ptc.auxValues[groupingColumn] == uGroupValue] = i
2090 else:
2091 means = np.zeros((nPair, len(ptc.ampNames)))
2092 exptimes = np.zeros_like(means)
2093 for i, ampName in enumerate(ptc.ampNames):
2094 means[:, i] = ptc.rawMeans[ampName] * ptc.gain[ampName]
2095 exptimes[:, i] = ptc.rawExpTimes[ampName]
2096 detMeans = np.nanmean(means, axis=1)
2097 detExptimes = np.nanmean(exptimes, axis=1)
2099 if autoGroupingUseExptime:
2100 abscissa = detExptimes
2101 else:
2102 abscissa = ptc.photoCharges[ptc.ampNames[0]].copy()
2103 # Set illegal photocharges to NaN.
2104 abscissa[abscissa < minPhotodiodeCurrent] = np.nan
2106 ratio = detMeans / abscissa
2107 ratio /= np.nanmedian(ratio)
2109 # Adjust those that are above threshold so they fall into the
2110 # largest group.
2111 above = (detMeans > autoGroupingMaxSignalFraction*np.nanmax(detMeans))
2112 maxIndex = np.argmax(detMeans[~above])
2113 ratio[above] = ratio[maxIndex]
2115 # The clustering of ratios into groups is performed with a simple
2116 # algorithm based on sorting and looking for the largest gaps.
2117 # See https://stackoverflow.com/a/18385795
2118 st = np.argsort(ratio)
2119 stratio = ratio[st]
2120 delta = stratio[1:] - stratio[0: -1]
2122 transitions, = np.where(delta > autoGroupingThreshold)
2123 # If there are no transitions then everything ends up in group 0.
2124 if len(transitions) > 0:
2125 ratioCuts = stratio[transitions] + autoGroupingThreshold/2.
2127 for i in range(len(transitions) + 1):
2128 if i == 0:
2129 inGroup, = np.where(ratio < ratioCuts[i])
2130 elif i == len(transitions):
2131 inGroup, = np.where(ratio > ratioCuts[i - 1])
2132 else:
2133 inGroup, = np.where((ratio > ratioCuts[i - 1]) & (ratio < ratioCuts[i]))
2134 groupingValues[inGroup] = i
2136 # Ensure out-of-range photoCharges/exptimes are in their own group.
2137 # These are masked later on.
2138 groupingValues[~np.isfinite(abscissa)] = -1
2139 # And put the high ones in the max group.
2140 groupingValues[above] = groupingValues[maxIndex]
2142 return groupingValues
2145def _computeTurnoffAndMax(
2146 abscissa,
2147 ordinate,
2148 initialMask,
2149 groupingValues,
2150 ampName="UNKNOWN",
2151 minSignalFitLinearityTurnoff=1000.0,
2152 maxFracLinearityDeviation=0.01,
2153 log=None,
2154 maxTurnoff=np.inf,
2155 use_all_for_normalization=False,
2156):
2157 """Compute the turnoff and max signal.
2159 Parameters
2160 ----------
2161 abscissa : `np.ndarray`
2162 Input x values, either photoCharges or exposure times.
2163 These should be cleaned of any non-finite values.
2164 ordinate : `np.ndarray`
2165 Input y values, the raw mean values for the amp.
2166 These should be cleaned of any non-finite values.
2167 initialMask : `np.ndarray`
2168 Mask to use for initial fit (usually from PTC).
2169 groupingValues : `np.ndarray`
2170 Array of values that are used to group different fits.
2171 ampName : `str`, optional
2172 Amplifier name (used for logging).
2173 minSignalFitLinearityTurnoff : `float`, optional
2174 Minimum signal to cmpute raw linearity slope for linearityTurnoff.
2175 maxFracLinearityDeviation : `float`, optional
2176 Maximum fraction deviation from raw linearity to compute turnoff.
2177 log : `logging.Logger`, optional
2178 Log object.
2179 maxTurnoff : `float`, optional
2180 Maximum turnoff allowed (will be set above PTC turnoff).
2181 use_all_for_normalization : `bool`, optional
2182 Use all the points (not just below turnoff) for normalization;
2183 this is for compatibility with the old linearizer fits.
2185 Returns
2186 -------
2187 turnoffIndex : `int`
2188 Fit turnoff index (keyed to raw input).
2189 turnoff : `float`
2190 Fit turnoff value.
2191 maxSignal : `float`
2192 Fit maximum signal value.
2193 goodPoints : `np.ndarray`
2194 Mask of good points used in turnoff fit.
2195 """
2196 if log is None:
2197 log = logging.getLogger(__name__)
2199 # Follow eo_pipe:
2200 # https://github.com/lsst-camera-dh/eo_pipe/blob/6afa546569f622b8d604921e248200481c445730/python/lsst/eo/pipe/linearityPlotsTask.py#L50
2201 # Replacing flux with abscissa, Ne with ordinate.
2203 # Fit a line with the y-intercept fixed to zero, using the
2204 # signal counts Ne as the variance in the chi-square, i.e.,
2205 # chi2 = sum( (ordinate - aa*abscissa)**2/ordinate )
2206 # Minimizing chi2 wrt aa, gives
2207 # aa = sum(abscissa) / sum(abscissa**2/ordinate)
2209 fitMask = initialMask.copy()
2210 fitMask[ordinate < minSignalFitLinearityTurnoff] = False
2211 fitMask[ordinate > maxTurnoff] = False
2212 fitMask[~np.isfinite(abscissa) | ~np.isfinite(ordinate)] = False
2213 goodPoints = fitMask.copy()
2215 gValues = np.unique(groupingValues)
2216 groupIndicesList = []
2217 for gValue in gValues:
2218 use, = np.where(groupingValues == gValue)
2219 groupIndicesList.append(use)
2221 found = False
2222 firstIteration = True
2223 while (fitMask.sum() >= 4) and not found:
2224 residuals = np.zeros_like(ordinate)
2226 abscissaMasked = abscissa.copy()
2227 abscissaMasked[~fitMask] = np.nan
2228 ordinateMasked = ordinate.copy()
2229 ordinateMasked[~fitMask] = np.nan
2231 for i, groupIndices in enumerate(groupIndicesList):
2232 num = np.nansum(abscissaMasked[groupIndices])
2233 denom = np.nansum(abscissaMasked[groupIndices]**2./ordinateMasked[groupIndices])
2235 if num == 0.0 or denom == 0.0:
2236 if firstIteration:
2237 log.info(
2238 "All points for %s were masked in linearity turnoff for group %d (first iteration).",
2239 ampName,
2240 i,
2241 )
2242 # We can try to recover this.
2243 nTry = min(10, len(groupIndices))
2244 num = np.nansum(abscissa[groupIndices][0: nTry])
2245 denom = np.nansum(abscissa[groupIndices][0: nTry]**2./ordinate[groupIndices][0: nTry])
2246 aa = num / denom
2247 else:
2248 log.warning("All points masked in linearity turnoff for group %d.", i)
2249 aa = np.nan
2250 else:
2251 aa = num / denom
2253 residuals[groupIndices] = (ordinate[groupIndices] - aa*abscissa[groupIndices]) / \
2254 ordinate[groupIndices]
2256 # Use the residuals to compute the turnoff.
2257 if use_all_for_normalization:
2258 residuals -= np.nanmedian(residuals)
2259 else:
2260 # Only subtract off the median from the previously estimated
2261 # fitMask.
2262 residuals -= np.nanmedian(residuals[fitMask])
2264 goodPoints = (np.abs(residuals) < maxFracLinearityDeviation) & (ordinate < maxTurnoff)
2266 if goodPoints.sum() > 4:
2267 # This was an adequate fit.
2268 found = True
2269 turnoff = np.max(ordinate[goodPoints])
2270 turnoffIndex = np.where(np.isclose(ordinate, turnoff))[0][0]
2271 else:
2272 # This was a bad fit; remove the largest outlier.
2273 badIndex = np.argmax(np.abs(residuals)[fitMask])
2274 fitIndices, = np.nonzero(fitMask)
2275 fitMask[fitIndices[badIndex]] = False
2277 firstIteration = False
2279 if not found:
2280 # Could not find any reasonable value.
2281 log.warning(
2282 "Could not find a reasonable initial linear fit to compute linearity turnoff for "
2283 "amplifier %s; may need finer sampling of input data?",
2284 ampName,
2285 )
2286 if np.all(~fitMask):
2287 return -1, np.nan, np.nan, goodPoints
2289 turnoff = np.max(ordinate[fitMask])
2290 turnoffIndex = np.where(np.isclose(ordinate, turnoff))[0][0]
2292 residuals = np.zeros(len(ordinate))
2294 # Fit the maximum signal.
2295 if turnoffIndex == (len(residuals) - 1):
2296 # This is the last point; we can't do a fit.
2297 # This is not a warning because we do not actually need this
2298 # value in practice.
2299 log.info(
2300 "No linearity turnoff detected for amplifier %s; try to increase the signal range.",
2301 ampName,
2302 )
2303 maxSignal = ordinate[turnoffIndex]
2304 else:
2305 maxSignalInitial = np.nanmax(ordinate)
2307 highFluxPoints = (np.nan_to_num(ordinate)
2308 > (1.0 - maxFracLinearityDeviation)*maxSignalInitial)
2309 maxSignal = np.median(ordinate[highFluxPoints])
2311 return turnoffIndex, turnoff, maxSignal, goodPoints
2314def _noderator(turnoff0, turnoff1, turnoff2, minNode, lowNodeSize, midNodeSize, highNodeSize):
2315 """The "noderator" node-finder.
2317 Parameters
2318 ----------
2319 turnoff0 : `float`
2320 Zeroth turnoff value (e.g. expectation of low-level
2321 non-linearity threshold) (adu).
2322 turnoff1 : `float`
2323 First turnoff value (e.g. ptc turnoff) (adu).
2324 turnoff2 : `float`
2325 Second turnoff value (e.g. linearity turnoff) (adu).
2326 minNode : `float`
2327 Location to place the first node after 0.0 (if this is <= 0.0
2328 it will be ignored) (adu).
2329 lowNodeSize : `float`
2330 Minimum node size in the low-level non-linearity regime
2331 (below turnoff0) (adu).
2332 midNodeSize : `float`
2333 Minimum node size in the mid-level non-linearity regime
2334 (between turnoff0 and turnoff1) (adu).
2335 highNodeSize : `float`
2336 Minimum node size in the high-level non-linearity regime
2337 (between turnoff1 and turnoff2) (adu).
2339 Returns
2340 -------
2341 nodes : `np.ndarray`
2342 Array of node values (adu).
2343 """
2344 if turnoff0 > minNode:
2345 # At least 2 nodes (edges) in the low signal regime.
2346 nNodesLow = np.clip(int(np.ceil((turnoff0 - minNode) / lowNodeSize)), 2, None)
2347 midStart = turnoff0
2348 else:
2349 nNodesLow = 0
2350 midStart = 0.0
2351 # At least 5 nodes (akima minimum) in the mid signal regime.
2352 nNodesMid = np.clip(int(np.ceil((turnoff1 - midStart) / midNodeSize)), 5, None)
2353 if turnoff2 > turnoff1:
2354 # At least 2 nodes (edges) in the high signal regime.
2355 nNodesHigh = np.clip(int(np.ceil((turnoff2 - turnoff1) / highNodeSize)), 3, None)
2356 else:
2357 nNodesHigh = 0
2358 nodesLow = np.linspace(minNode, turnoff0, nNodesLow)
2359 nodesMid = np.linspace(midStart, turnoff1, nNodesMid)
2360 nodesHigh = np.linspace(turnoff1, turnoff2, nNodesHigh)
2362 # Make sure we do not duplicate nodes when concatenating.
2363 nodeList = []
2364 if nNodesLow > 1:
2365 nodeList.append([0.0])
2366 nodeList.append(nodesLow[:-1])
2367 if nNodesMid > 1:
2368 nodeList.append(nodesMid)
2369 if nNodesHigh > 1:
2370 nodeList.append(nodesHigh[1:])
2371 return np.concatenate(nodeList)