Coverage for python/lsst/cp/pipe/cpLinearitySolve.py: 11%
843 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 02:31 -0700
1# This file is part of cp_pipe.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21#
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()
516 if ampName in inputPtc.badAmps:
517 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
518 self.log.warning("Amp %s in detector %s has no usable PTC information. Skipping!",
519 ampName, detector.getName())
520 continue
522 # Check for too few points.
523 if self.config.linearityType == "Spline" \
524 and self.config.splineGroupingColumn is not None \
525 and len(inputPtc.inputExpIdPairs[ampName]) < self.config.splineGroupingMinPoints:
526 raise RuntimeError(
527 "The input PTC has too few points to reliably run with PD grouping. "
528 "The recommended course of action is to set splineGroupingColumn to None. "
529 "If you really know what you are doing, you may reduce "
530 "config.splineGroupingMinPoints.")
532 # We start with all finite values.
533 mask = np.isfinite(inputPtc.rawMeans[ampName])
535 if self.config.linearityType == "Spline" and temperatureValues is not None:
536 mask &= np.isfinite(temperatureValues)
538 if self.config.usePhotodiode:
539 modExpTimes = inputPtc.photoCharges[ampName].copy()
540 # Make sure any exposure pairs that do not have photodiode data
541 # are masked.
542 mask[~np.isfinite(modExpTimes)] = False
544 # Make sure any photodiode measurements below the configured
545 # minimum are masked.
546 mask[modExpTimes < self.config.minPhotodiodeCurrent] = False
548 # Get the photodiode correction.
549 if self.config.applyPhotodiodeCorrection:
550 for j, pair in enumerate(inputPtc.inputExpIdPairs[ampName]):
551 try:
552 correction = abscissaCorrections[str(pair)]
553 except KeyError:
554 correction = 0.0
555 modExpTimes[j] += correction
557 inputAbscissa = modExpTimes
558 else:
559 inputAbscissa = inputPtc.rawExpTimes[ampName].copy()
561 # Normalize if configured.
562 inputNorm = np.ones_like(inputAbscissa)
563 if self.config.useFocalPlaneNormalization:
564 exposures = np.asarray(inputPtc.inputExpIdPairs[ampName])[:, 0]
565 a, b = esutil.numpy_util.match(exposures, inputNormalization["exposure"])
566 inputNorm[a] = inputNormalization["normalization"][b]
567 inputAbscissa *= inputNorm
569 # Compute linearityTurnoff and linearitySignalMax.
570 turnoffMask = inputPtc.expIdMask[ampName].copy()
571 turnoffMask &= mask
573 turnoffIndex, turnoff, maxSignal, _ = _computeTurnoffAndMax(
574 inputAbscissa,
575 inputPtc.rawMeans[ampName],
576 turnoffMask,
577 groupingValues,
578 ampName=ampName,
579 minSignalFitLinearityTurnoff=self.config.minSignalFitLinearityTurnoff,
580 maxFracLinearityDeviation=self.config.maxFracLinearityDeviation,
581 log=self.log,
582 )
584 if np.isnan(turnoff):
585 # This is a bad amp, with no linearizer.
586 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
587 self.log.warning("Amp %s in detector %s has no usable linearizer information. Skipping!",
588 ampName, detector.getName())
589 continue
591 linearizer.linearityTurnoff[ampName] = turnoff
592 linearizer.linearityMaxSignal[ampName] = maxSignal
594 inputOrdinate = inputPtc.rawMeans[ampName].copy()
596 linearizer.inputAbscissa[ampName] = inputAbscissa.copy()
597 linearizer.inputOrdinate[ampName] = inputOrdinate.copy()
598 linearizer.inputGroupingIndex[ampName] = groupingValues.copy()
599 linearizer.inputNormalization[ampName] = inputNorm.copy()
601 if self.config.linearityType != 'Spline':
602 mask &= (inputOrdinate < self.config.maxLinearAdu)
603 else:
604 # For spline fits, cut above the turnoff.
605 self.log.info("Using linearityTurnoff of %.4f adu for amplifier %s", turnoff, ampName)
606 extraMask = np.ones(len(inputOrdinate), dtype=bool)
607 extraMask[turnoffIndex + 1:] = False
608 mask &= extraMask
610 mask &= (inputOrdinate > self.config.minLinearAdu)
612 # Initial value for the mask.
613 linearizer.inputMask[ampName] = mask.copy()
615 if mask.sum() < 2:
616 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
617 self.log.warning("Amp %s in detector %s has not enough points for fit. Skipping!",
618 ampName, detector.getName())
619 continue
621 if self.config.linearityType != 'Spline':
622 linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 100.0], inputAbscissa[mask],
623 inputOrdinate[mask], funcPolynomial)
625 # Convert this proxy-to-flux fit into an expected linear flux
626 linearOrdinate = linearFit[0] + linearFit[1] * inputAbscissa
627 # Exclude low end outliers.
628 # This is compared to the original values.
629 threshold = self.config.nSigmaClipLinear * np.sqrt(abs(inputOrdinate))
631 mask[np.abs(inputOrdinate - linearOrdinate) >= threshold] = False
633 linearizer.inputMask[ampName] = mask.copy()
635 if mask.sum() < 2:
636 linearizer = self.fillBadAmp(linearizer, fitOrder, inputPtc, amp)
637 self.log.warning("Amp %s in detector %s has not enough points in linear ordinate. "
638 "Skipping!", ampName, detector.getName())
639 continue
641 self.debugFit('linearFit', inputAbscissa, inputOrdinate, linearOrdinate, mask, ampName)
643 # Do fits
644 if self.config.linearityType in ['Polynomial', 'Squared', 'LookupTable']:
645 polyFit = np.zeros(fitOrder + 1)
646 polyFit[1] = 1.0
647 polyFit, polyFitErr, chiSq, weights = irlsFit(polyFit, linearOrdinate[mask],
648 inputOrdinate[mask], funcPolynomial)
650 # Truncate the polynomial fit to the squared term.
651 k1 = polyFit[1]
652 linearityCoeffs = np.array(
653 [-coeff/(k1**order) for order, coeff in enumerate(polyFit)]
654 )[2:]
655 significant = np.where(np.abs(linearityCoeffs) > 1e-10)
656 self.log.info("Significant polynomial fits: %s", significant)
658 modelOrdinate = funcPolynomial(polyFit, linearOrdinate)
660 self.debugFit(
661 'polyFit',
662 inputAbscissa[mask],
663 inputOrdinate[mask],
664 modelOrdinate[mask],
665 None,
666 ampName,
667 )
669 if self.config.linearityType == 'Squared':
670 # The first term is the squared term.
671 linearityCoeffs = linearityCoeffs[0: 1]
672 elif self.config.linearityType == 'LookupTable':
673 # Use linear part to get time at which signal is
674 # maxAduForLookupTableLinearizer DN
675 tMax = (self.config.maxLookupTableAdu - polyFit[0])/polyFit[1]
676 timeRange = np.linspace(0, tMax, self.config.maxLookupTableAdu)
677 signalIdeal = polyFit[0] + polyFit[1]*timeRange
678 signalUncorrected = funcPolynomial(polyFit, timeRange)
679 lookupTableRow = signalIdeal - signalUncorrected # LinearizerLookupTable has correction
681 linearizer.tableData[tableIndex, :] = lookupTableRow
682 linearityCoeffs = np.array([tableIndex, 0])
683 tableIndex += 1
684 elif self.config.linearityType in ['Spline']:
685 # This is a spline fit with photodiode data based on a model
686 # from Pierre Astier.
687 # This model fits a spline with (optional) nuisance parameters
688 # to allow for different linearity coefficients with different
689 # photodiode settings. The minimization is a least-squares
690 # fit with the residual of
691 # Sum[(S(mu_i) + mu_i - O)/(k_j * D_i) - 1]**2, where S(mu_i)
692 # is an Akima Spline function of mu_i, the observed flat-pair
693 # mean; D_j is the photo-diode measurement corresponding to
694 # that flat-pair; and k_j is a constant of proportionality
695 # which is over index j as it is allowed to
696 # be different based on different photodiode settings (e.g.
697 # CCOBCURR); and O is a constant offset to allow for light
698 # leaks (and is only fit if doSplineFitOffset=True). In
699 # addition, if config.doSplineFitTemperature is True then
700 # the fit will adjust mu such that
701 # mu = mu_input*(1 + alpha*(T - T_ref))
702 # and T_ref is taken as the median temperature of the run.
703 # Finally, if config.doSplineFitTemporal is True then the
704 # fit will further adjust mu such that
705 # mu = mu_input*(1 + beta*(mjd - mjd_ref))
706 # and mjd_ref is taken as the median mjd of the run.
707 # Note that this fit is only valid if the input data
708 # was taken with a randomly shuffled order of exposure
709 # levels.
711 # The fit has additional constraints to ensure that the spline
712 # goes through the (0, 0) point, as well as a normalization
713 # condition so that the average of the spline over the full
714 # range is 0. The normalization ensures that the spline only
715 # fits deviations from linearity, rather than the linear
716 # function itself which is degenerate with the gain.
718 # We want to make sure the top node is above the top value
719 # to avoid edge issues with the top point.
720 nodes = np.linspace(0.0, np.max(inputOrdinate[mask]) + 1.0, self.config.splineKnots)
722 if temperatureValues is not None:
723 temperatureValuesScaled = temperatureValues - np.median(temperatureValues[mask])
724 else:
725 temperatureValuesScaled = None
727 if self.config.doSplineFitTemporal:
728 inputMjdScaled = inputPtc.inputExpPairMjdStartList[ampName].copy()
729 inputMjdScaled -= np.nanmedian(inputMjdScaled)
730 else:
731 inputMjdScaled = None
733 fitter = AstierSplineLinearityFitter(
734 nodes,
735 groupingValues,
736 inputAbscissa,
737 inputOrdinate,
738 mask=mask,
739 log=self.log,
740 fit_offset=self.config.doSplineFitOffset,
741 fit_weights=self.config.doSplineFitWeights,
742 weight_pars_start=self.config.splineFitWeightParsStart,
743 fit_temperature=self.config.doSplineFitTemperature,
744 temperature_scaled=temperatureValuesScaled,
745 max_signal_nearly_linear=inputPtc.ptcTurnoff[ampName],
746 fit_temporal=self.config.doSplineFitTemporal,
747 mjd_scaled=inputMjdScaled,
748 )
749 p0 = fitter.estimate_p0()
750 pars = fitter.fit(
751 p0,
752 min_iter=self.config.splineFitMinIter,
753 max_iter=self.config.splineFitMaxIter,
754 max_rejection_per_iteration=self.config.splineFitMaxRejectionPerIteration,
755 n_sigma_clip=self.config.nSigmaClipLinear,
756 )
758 # Confirm that the first parameter is 0, and set it to
759 # exactly zero.
760 if not np.isclose(pars[0], 0):
761 raise RuntimeError("Programmer error! First spline parameter must "
762 "be consistent with zero.")
763 pars[0] = 0.0
765 linearityChisq = fitter.compute_chisq_dof(pars)
767 linearityCoeffs = np.concatenate([nodes, pars[fitter.par_indices["values"]]])
768 linearFit = np.array([0.0, np.mean(pars[fitter.par_indices["groups"]])])
770 # We must modify the inputOrdinate according to the
771 # nuisance terms in the linearity fit for the residual
772 # computation code to work properly.
773 # The true mu (inputOrdinate) is given by
774 # mu = mu_in * (1 + alpha*t_scale) * (1 + beta*mjd_scale)
775 if self.config.doSplineFitTemperature:
776 inputOrdinate *= (1.0
777 + pars[fitter.par_indices["temperature_coeff"]]*temperatureValuesScaled)
778 if self.config.doSplineFitTemporal:
779 inputOrdinate *= (1.0
780 + pars[fitter.par_indices["temporal_coeff"]]*inputMjdScaled)
781 # We have to adjust the abscissa for the different groups.
782 # This is because we need a corrected abscissa to get a
783 # reasonable linear fit to look for residuals, particularly in
784 # the case of significantly different signal-vs-photodiode or
785 # signal-vs-exptime scalings for different groups. This then
786 # becomes a multiplication by the relative scaling of the
787 # different groups.
788 for j, group_index in enumerate(fitter.group_indices):
789 inputAbscissa[group_index] *= (pars[fitter.par_indices["groups"][j]] / linearFit[1])
790 # And remove the offset term.
791 if self.config.doSplineFitOffset:
792 inputOrdinate -= pars[fitter.par_indices["offset"]]
794 linearOrdinate = linearFit[1] * inputOrdinate
795 # For the spline fit, reuse the "polyFit -> fitParams"
796 # field to record the linear coefficients for the groups.
797 # We additionally append the offset and weight_pars;
798 # however these will be zero-length arrays if these were
799 # not configured to be fit.
800 polyFit = np.concatenate((
801 pars[fitter.par_indices["groups"]],
802 pars[fitter.par_indices["offset"]],
803 pars[fitter.par_indices["weight_pars"]],
804 pars[fitter.par_indices["temperature_coeff"]],
805 pars[fitter.par_indices["temporal_coeff"]],
806 ))
807 polyFitErr = np.zeros_like(polyFit)
808 chiSq = linearityChisq
810 # Update mask based on what the fitter rejected.
811 mask = fitter.mask
813 linearizer.inputMask[ampName] = mask.copy()
814 else:
815 polyFit = np.zeros(1)
816 polyFitErr = np.zeros(1)
817 chiSq = np.nan
818 linearityCoeffs = np.zeros(1)
820 linearizer.linearityType[ampName] = self.config.linearityType
821 linearizer.linearityCoeffs[ampName] = linearityCoeffs
822 if self.config.trimmedState:
823 linearizer.linearityBBox[ampName] = amp.getBBox()
824 else:
825 linearizer.linearityBBox[ampName] = amp.getRawBBox()
826 linearizer.fitParams[ampName] = polyFit
827 linearizer.fitParamsErr[ampName] = polyFitErr
828 linearizer.fitChiSq[ampName] = chiSq
829 linearizer.linearFit[ampName] = linearFit
831 image = afwImage.ImageF(len(inputOrdinate), 1)
832 image.array[:, :] = inputOrdinate
833 linearizeFunction = linearizer.getLinearityTypeByName(linearizer.linearityType[ampName])
834 linearizeFunction()(
835 image,
836 **{'coeffs': linearizer.linearityCoeffs[ampName],
837 'table': linearizer.tableData,
838 'log': linearizer.log}
839 )
840 linearizeModel = image.array[0, :]
842 # The residuals that we record are the final residuals compared to
843 # a linear model, after everything has been (properly?) linearized.
844 if mask.sum() < 2:
845 self.log.warning("Amp %s in detector %s has not enough points in linear ordinate "
846 "for residuals. Skipping!", ampName, detector.getName())
847 residuals = np.full_like(linearizeModel, np.nan)
848 residualsUnmasked = residuals.copy()
849 else:
850 postLinearFit, _, _, _ = irlsFit(
851 linearFit,
852 inputAbscissa[mask],
853 linearizeModel[mask],
854 funcPolynomial,
855 )
856 # When computing residuals, we only care about the slope of
857 # the postLinearFit and not the intercept. The intercept
858 # itself depends on a possibly unknown zero in the abscissa
859 # (often photodiode) which may have an arbitrary value.
860 residuals = linearizeModel - (postLinearFit[1] * inputAbscissa)
861 residualsUnmasked = residuals.copy()
862 # We set masked residuals to nan.
863 residuals[~mask] = np.nan
865 linearizer.fitResidualsUnmasked[ampName] = residualsUnmasked
866 linearizer.fitResiduals[ampName] = residuals
867 linearizer.fitResidualsModel[ampName] = linearizeModel.copy()
869 finite = np.isfinite(residuals)
870 if finite.sum() == 0:
871 sigmad = np.nan
872 else:
873 sigmad = median_abs_deviation(residuals[finite]/inputOrdinate[finite], scale="normal")
874 linearizer.fitResidualsSigmaMad[ampName] = sigmad
876 self.debugFit(
877 'solution',
878 inputOrdinate[mask],
879 linearOrdinate[mask],
880 linearizeModel[mask],
881 None,
882 ampName,
883 )
885 self.fixupBadAmps(linearizer)
887 linearizer.hasLinearity = True
888 linearizer.validate()
889 linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
890 linearizer.updateMetadata(setDate=True, setCalibId=True)
891 linearizer.updateMetadataFromExposures([inputPtc])
892 provenance = IsrProvenance(calibType='linearizer')
894 return pipeBase.Struct(
895 outputLinearizer=linearizer,
896 outputProvenance=provenance,
897 )
899 def fillBadAmp(self, linearizer, fitOrder, inputPtc, amp):
900 # Need to fill linearizer with empty values
901 # if the amp is non-functional
902 ampName = amp.getName()
903 nEntries = 1
904 pEntries = 1
905 if self.config.linearityType in ['Polynomial']:
906 # We discard the first 2 entries in the polynomial.
907 nEntries = fitOrder + 1 - 2
908 pEntries = fitOrder + 1 - 2
909 elif self.config.linearityType in ['Spline']:
910 nEntries = fitOrder * 2
911 elif self.config.linearityType in ['Squared', 'None']:
912 nEntries = 1
913 pEntries = fitOrder + 1
914 elif self.config.linearityType in ['LookupTable']:
915 nEntries = 2
916 pEntries = fitOrder + 1
918 nPair = len(inputPtc.inputExpIdPairs[ampName])
920 linearizer.linearityType[ampName] = "None"
921 linearizer.linearityCoeffs[ampName] = np.zeros(nEntries)
922 if self.config.trimmedState:
923 linearizer.linearityBBox[ampName] = amp.getBBox()
924 else:
925 linearizer.linearityBBox[ampName] = amp.getRawBBox()
926 linearizer.fitParams[ampName] = np.zeros(pEntries)
927 linearizer.fitParamsErr[ampName] = np.zeros(pEntries)
928 linearizer.fitChiSq[ampName] = np.nan
929 linearizer.fitResiduals[ampName] = np.zeros(nPair)
930 linearizer.fitResidualsSigmaMad[ampName] = np.nan
931 linearizer.fitResidualsUnmasked[ampName] = np.zeros(nPair)
932 linearizer.fitResidualsModel[ampName] = np.zeros(nPair)
933 linearizer.linearFit[ampName] = np.zeros(2)
934 linearizer.linearityTurnoff[ampName] = np.nan
935 linearizer.linearityMaxSignal[ampName] = np.nan
936 linearizer.inputMask[ampName] = np.zeros(nPair, dtype=np.bool_)
937 linearizer.inputAbscissa[ampName] = np.zeros(nPair)
938 linearizer.inputOrdinate[ampName] = np.zeros(nPair)
939 linearizer.inputGroupingIndex[ampName] = np.zeros(nPair, dtype=np.int64)
940 linearizer.inputNormalization[ampName] = np.ones(nPair)
942 return linearizer
944 def fixupBadAmps(self, linearizer):
945 """Fix nan padding in bad amplifiers.
947 Parameters
948 ----------
949 linearizer : `lsst.ip.isr.Linearizer`
950 """
951 fitParamsMaxLen = 0
952 for ampName in linearizer.ampNames:
953 if (length := len(linearizer.fitParams[ampName])) > fitParamsMaxLen:
954 fitParamsMaxLen = length
956 for ampName in linearizer.ampNames:
957 if linearizer.linearityType[ampName] == "None":
958 # Bad amplifier.
959 linearizer.fitParams[ampName] = np.zeros(fitParamsMaxLen)
960 linearizer.fitParamsErr[ampName] = np.zeros(fitParamsMaxLen)
961 elif len(linearizer.fitParams[ampName]) != fitParamsMaxLen:
962 raise RuntimeError("Linearity has mismatched fitParams; check code/data.")
964 def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName):
965 """Debug method for linearity fitting.
967 Parameters
968 ----------
969 stepname : `str`
970 A label to use to check if we care to debug at a given
971 line of code.
972 xVector : `numpy.array`, (N,)
973 The values to use as the independent variable in the
974 linearity fit.
975 yVector : `numpy.array`, (N,)
976 The values to use as the dependent variable in the
977 linearity fit.
978 yModel : `numpy.array`, (N,)
979 The values to use as the linearized result.
980 mask : `numpy.array` [`bool`], (N,) , optional
981 A mask to indicate which entries of ``xVector`` and
982 ``yVector`` to keep.
983 ampName : `str`
984 Amplifier name to lookup linearity correction values.
985 """
986 frame = getDebugFrame(self._display, stepname)
987 if frame:
988 import matplotlib.pyplot as plt
989 fig, axs = plt.subplots(2)
991 if mask is None:
992 mask = np.ones_like(xVector, dtype=bool)
994 fig.suptitle(f"{stepname} {ampName} {self.config.linearityType}")
995 if stepname == 'linearFit':
996 axs[0].set_xlabel("Input Abscissa (time or mondiode)")
997 axs[0].set_ylabel("Input Ordinate (flux)")
998 axs[1].set_xlabel("Linear Ordinate (linear flux)")
999 axs[1].set_ylabel("Flux Difference: (input - linear)")
1000 elif stepname in ('polyFit', 'splineFit'):
1001 axs[0].set_xlabel("Linear Abscissa (linear flux)")
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 - full model fit)")
1005 elif stepname == 'solution':
1006 axs[0].set_xlabel("Input Abscissa (time or mondiode)")
1007 axs[0].set_ylabel("Linear Ordinate (linear flux)")
1008 axs[1].set_xlabel("Model flux (linear flux)")
1009 axs[1].set_ylabel("Flux Difference: (linear - model)")
1011 axs[0].set_yscale('log')
1012 axs[0].set_xscale('log')
1013 axs[0].scatter(xVector, yVector)
1014 axs[0].scatter(xVector[~mask], yVector[~mask], c='red', marker='x')
1015 axs[1].set_xscale('log')
1017 axs[1].scatter(yModel, yVector[mask] - yModel)
1018 fig.tight_layout()
1019 fig.show()
1021 prompt = "Press Enter or c to continue [chpx]..."
1022 while True:
1023 ans = input(prompt).lower()
1024 if ans in ("", " ", "c",):
1025 break
1026 elif ans in ("p", ):
1027 import pdb
1028 pdb.set_trace()
1029 elif ans in ("h", ):
1030 print("[h]elp [c]ontinue [p]db")
1031 elif ans in ('x', ):
1032 exit()
1033 plt.close()
1036class LinearityDoubleSplineSolveConnections(
1037 pipeBase.PipelineTaskConnections,
1038 dimensions=("instrument", "detector"),
1039):
1040 dummy = cT.Input(
1041 name="raw",
1042 doc="Dummy exposure.",
1043 storageClass='Exposure',
1044 dimensions=("instrument", "exposure", "detector"),
1045 multiple=True,
1046 deferLoad=True,
1047 )
1048 camera = cT.PrerequisiteInput(
1049 name="camera",
1050 doc="Camera Geometry definition.",
1051 storageClass="Camera",
1052 dimensions=("instrument", ),
1053 isCalibration=True,
1054 )
1055 inputLinearizerPtc = cT.Input(
1056 name="linearizerPtc",
1057 doc="Input linearizer PTC dataset.",
1058 storageClass="PhotonTransferCurveDataset",
1059 dimensions=("instrument", "detector"),
1060 isCalibration=True,
1061 )
1062 inputNormalization = cT.Input(
1063 name="cpLinearizerPtcNormalization",
1064 doc="Focal-plane normalization table.",
1065 storageClass="ArrowAstropy",
1066 dimensions=("instrument",),
1067 isCalibration=True,
1068 )
1069 inputBinnedImagesHandles = cT.Input(
1070 name="cpPtcPairBinned",
1071 doc="Tabulated binned exposure pairs.",
1072 storageClass="ArrowAstropy",
1073 dimensions=("instrument", "exposure", "detector"),
1074 multiple=True,
1075 deferLoad=True,
1076 )
1077 outputLinearizer = cT.Output(
1078 name="linearity",
1079 doc="Output linearity measurements.",
1080 storageClass="Linearizer",
1081 dimensions=("instrument", "detector"),
1082 isCalibration=True,
1083 )
1085 def __init__(self, *, config=None):
1086 super().__init__(config=config)
1088 if not config.useFocalPlaneNormalization:
1089 del self.inputNormalization
1092class LinearityDoubleSplineSolveConfig(
1093 pipeBase.PipelineTaskConfig,
1094 pipelineConnections=LinearityDoubleSplineSolveConnections,
1095):
1096 maxFracLinearityDeviation = pexConfig.Field(
1097 dtype=float,
1098 doc="Maximum fraction deviation from raw linearity to compute "
1099 "linearityTurnoff and linearityMaxSignal.",
1100 # TODO: DM-46811 investigate if this can be raised to 0.05.
1101 default=0.01,
1102 )
1103 minSignalFitLinearityTurnoff = pexConfig.Field(
1104 dtype=float,
1105 doc="Minimum signal to compute raw linearity slope for linearityTurnoff.",
1106 default=1000.0,
1107 )
1108 usePhotodiode = pexConfig.Field(
1109 dtype=bool,
1110 doc="Use the photodiode info instead of the raw expTimes?",
1111 default=False,
1112 )
1113 minPhotodiodeCurrent = pexConfig.Field(
1114 dtype=float,
1115 doc="Minimum value to trust photodiode signals.",
1116 default=0.0,
1117 )
1118 doAutoGrouping = pexConfig.Field(
1119 dtype=bool,
1120 doc="Do automatic group detection? Cannot be True if splineGroupingColumn is also set. "
1121 "The automatic group detection will use the ratio of signal to exposure time (if "
1122 "autoGroupingUseExptime is True) or photodiode (if False) to determine which "
1123 "flat pairs were taken with different illumination settings.",
1124 default=False,
1125 )
1126 autoGroupingUseExptime = pexConfig.Field(
1127 dtype=bool,
1128 doc="Use exposure time to determine automatic grouping. Used if doAutoGrouping=True.",
1129 default=True,
1130 )
1131 autoGroupingThreshold = pexConfig.Field(
1132 dtype=float,
1133 doc="Minimum relative jump from sorted conversion values to determine a group.",
1134 default=0.1,
1135 )
1136 autoGroupingMaxSignalFraction = pexConfig.Field(
1137 dtype=float,
1138 doc="Only do auto-grouping when the signal is this fraction of the maximum signal. "
1139 "All exposures with signal higher than this threshold will be put into the "
1140 "largest signal group. This config is needed if the input PTC goes beyond "
1141 "the linearity turnoff.",
1142 default=0.9,
1143 )
1144 groupingColumn = pexConfig.Field(
1145 dtype=str,
1146 doc="Column to use for grouping together points, to allow "
1147 "for different proportionality constants. If None, then grouping will "
1148 "only be done if doAutoGrouping is True.",
1149 default=None,
1150 optional=True,
1151 )
1152 absoluteSplineMinimumSignalNode = pexConfig.Field(
1153 dtype=float,
1154 doc="Smallest node (above 0) for absolute spline (adu).",
1155 default=0.0,
1156 )
1157 absoluteSplineLowThreshold = pexConfig.Field(
1158 dtype=float,
1159 doc="Threshold for the low-level linearity nodes for absolute spline (adu). "
1160 "If this is below ``absoluteSplineMinimumSignalNode`` then the low "
1161 "level checks will be skipped.",
1162 default=0.0,
1163 )
1164 absoluteSplineLowNodeSize = pexConfig.Field(
1165 dtype=float,
1166 doc="Minimum size for low-level linearity nodes for absolute spline (adu).",
1167 default=2000.0,
1168 )
1169 absoluteSplineNodeSize = pexConfig.Field(
1170 dtype=float,
1171 doc="Minimum size for linearity nodes for absolute spline above absoluteSplineLowThreshold e(adu); "
1172 "note that there will always be a node at the reference PTC turnoff.",
1173 default=3000.0,
1174 )
1175 absoluteSplineFitMinIter = pexConfig.Field(
1176 dtype=int,
1177 doc="Minimum number of iterations for absolute spline fit.",
1178 default=3,
1179 )
1180 absoluteSplineFitMaxIter = pexConfig.Field(
1181 dtype=int,
1182 doc="Maximum number of iterations for absolute spline fit.",
1183 default=20,
1184 )
1185 absoluteSplineFitMaxRejectionPerIteration = pexConfig.Field(
1186 dtype=int,
1187 doc="Maximum number of rejections per iteration for absolute spline fit.",
1188 default=5,
1189 )
1190 absoluteNSigmaClipLinear = pexConfig.Field(
1191 dtype=float,
1192 doc="Sigma-clipping for absolute spline solution.",
1193 default=5.0,
1194 )
1195 doAbsoluteSplineFitOffset = pexConfig.Field(
1196 dtype=bool,
1197 doc="Fit a scattered light offset in the spline fit.",
1198 default=True,
1199 )
1200 doAbsoluteSplineFitWeights = pexConfig.Field(
1201 dtype=bool,
1202 doc="Fit linearity weight parameters in the spline fit.",
1203 default=False,
1204 )
1205 absoluteSplineFitWeightParsStart = pexConfig.ListField(
1206 dtype=float,
1207 doc="Starting parameters for weight fit, if doSplineFitWeights=True. "
1208 "Parameters are such that sigma = sqrt(par[0]**2. + par[1]**2./mu)."
1209 "If doSplineFitWeights=False then these are used as-is; otherwise "
1210 "they are used as the initial values for fitting these parameters.",
1211 length=2,
1212 default=[1.0, 0.0],
1213 )
1214 doAbsoluteSplineFitTemperature = pexConfig.Field(
1215 dtype=bool,
1216 doc="Fit temperature coefficient in spline fit?",
1217 default=False,
1218 )
1219 absoluteSplineFitTemperatureColumn = pexConfig.Field(
1220 dtype=str,
1221 doc="Name of the temperature column to use when fitting temperature "
1222 "coefficients in spline fit; this must not be None if "
1223 "doSplineFitTemperature is True.",
1224 default=None,
1225 optional=True,
1226 )
1227 doAbsoluteSplineFitTemporal = pexConfig.Field(
1228 dtype=bool,
1229 doc="Fit a linear temporal parameter coefficient in spline fit?",
1230 default=False,
1231 )
1232 useFocalPlaneNormalization = pexConfig.Field(
1233 dtype=bool,
1234 doc="Use focal-plane normalization in addition to/instead of photodiode? "
1235 "(Only used with for absolute spline fitting).",
1236 default=False,
1237 )
1238 relativeSplineReferenceCounts = pexConfig.Field(
1239 dtype=float,
1240 doc="Number of target counts (adu) to select a reference image for "
1241 "relative spline solution.",
1242 default=10000.0,
1243 )
1244 relativeSplineMinimumSignalNode = pexConfig.Field(
1245 dtype=float,
1246 doc="Smallest node (above 0) for relative spline (adu).",
1247 default=100.0,
1248 )
1249 relativeSplineLowThreshold = pexConfig.Field(
1250 dtype=float,
1251 doc="Threshold for the low-level linearity nodes for relative spline (adu)."
1252 "If this is below ``relativeSplineMinimumSignalNode`` then the low "
1253 "level checks will be skipped.",
1254 default=5000.0,
1255 )
1256 relativeSplineLowNodeSize = pexConfig.Field(
1257 dtype=float,
1258 doc="Minimum size for low-level linearity nodes for relative spline (adu).",
1259 default=750.0,
1260 )
1261 relativeSplineMidNodeSize = pexConfig.Field(
1262 dtype=float,
1263 doc="Minimum size for mid-level linearity nodes for relative spline (adu); "
1264 "this applies to counts between relativeSplineLowThreshold and the "
1265 "PTC turnoff.",
1266 default=5000.0,
1267 )
1268 relativeSplineHighNodeSize = pexConfig.Field(
1269 dtype=float,
1270 doc="Minimum size for high-level linearity nodes for relative spline (adu); "
1271 "this applies to counts between the PTC and linearity turnoffs.",
1272 default=2000.0,
1273 )
1274 relativeSplineFitMinIter = pexConfig.Field(
1275 dtype=int,
1276 doc="Minimum number of iterations for relative spline fit.",
1277 default=3,
1278 )
1279 relativeSplineFitMaxIter = pexConfig.Field(
1280 dtype=int,
1281 doc="Maximum number of iterations for relative spline fit.",
1282 default=20,
1283 )
1284 relativeSplineFitMaxRejectionPerIteration = pexConfig.Field(
1285 dtype=int,
1286 doc="Maximum number of rejections per iteration for relative spline fit.",
1287 default=5,
1288 )
1289 relativeNSigmaClipLinear = pexConfig.Field(
1290 dtype=float,
1291 doc="Sigma-clipping for relative spline solution.",
1292 default=5.0,
1293 )
1295 def validate(self):
1296 super().validate()
1298 if self.doAbsoluteSplineFitTemperature and self.absoluteSplineFitTemperatureColumn is None:
1299 raise ValueError(
1300 "Must set absoluteSplineFitTemperatureColumn if doAbsoluteSplineFitTemperature is True.",
1301 )
1303 if self.doAutoGrouping and self.groupingColumn is not None:
1304 raise ValueError("Must not set doAutoGrouping and also groupingColumn")
1305 if self.doAutoGrouping:
1306 if not self.autoGroupingUseExptime and not self.usePhotodiode:
1307 raise ValueError("If doAutoGrouping is True and autoGroupingUseExptime is False, then "
1308 "usePhotodiode must be True.")
1311class LinearityDoubleSplineSolveTask(pipeBase.PipelineTask):
1312 ConfigClass = LinearityDoubleSplineSolveConfig
1313 _DefaultName = "cpLinearityDoubleSplineSolve"
1315 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1316 # docstring inherited
1317 inputs = butlerQC.get(inputRefs)
1319 if self.config.useFocalPlaneNormalization:
1320 inputNormalization = inputs["inputNormalization"]
1321 else:
1322 inputNormalization = None
1324 # Add calibration provenance info to header.
1325 kwargs = dict()
1326 reference = getattr(inputRefs, "inputLinearizerPtc", None)
1328 if reference is not None and hasattr(reference, "run"):
1329 runKey = "PTC_RUN"
1330 runValue = reference.run
1331 idKey = "PTC_UUID"
1332 idValue = str(reference.id)
1333 dateKey = "PTC_DATE"
1334 calib = inputs.get("inputPtc", None)
1335 dateValue = extractCalibDate(calib)
1337 kwargs[runKey] = runValue
1338 kwargs[idKey] = idValue
1339 kwargs[dateKey] = dateValue
1341 self.log.info("Using " + str(reference.run))
1343 outputs = self.run(
1344 inputPtc=inputs["inputLinearizerPtc"],
1345 camera=inputs["camera"],
1346 inputBinnedImagesHandles=inputs["inputBinnedImagesHandles"],
1347 inputNormalization=inputNormalization,
1348 )
1349 outputs.outputLinearizer.updateMetadata(setDate=False, **kwargs)
1351 butlerQC.put(outputs, outputRefs)
1353 def run(
1354 self,
1355 *,
1356 inputPtc,
1357 camera,
1358 inputBinnedImagesHandles,
1359 inputNormalization,
1360 ):
1361 """Fit the double-spline relative/absolute linearity correction.
1363 Parameters
1364 ----------
1365 inputPtc : `lsst.ip.isr.PtcDataset`
1366 Pre-measured PTC dataset.
1367 camera : `lsst.afw.cameraGeom.Camera`
1368 Camera geometry.
1369 inputBinnedImagesHandles : `list` [`DeferredDatasetHandle`]
1370 Handles for input binned image pairs.
1371 inputNormalization : `astropy.table.Table`, optional
1372 Focal plane normalization table to use if
1373 useFocalPlaneNormalization is True.
1375 Returns
1376 -------
1377 results : `lsst.pipe.base.Struct`
1378 The results struct containing:
1380 ``outputLinearizer``
1381 Final linearizer calibration (`lsst.ip.isr.Linearizer`).
1382 ``outputProvenance``
1383 Provenance data for the new calibration
1384 (`lsst.ip.isr.IsrProvenance`).
1385 """
1386 detector = camera[inputPtc.metadata["DETECTOR"]]
1388 binnedImagesHandleDict = {
1389 handle.dataId["exposure"]: handle for handle in inputBinnedImagesHandles
1390 }
1392 linearizer = Linearizer(detector=detector, log=self.log)
1393 linearizer.updateMetadataFromExposures([inputPtc])
1395 groupingValues = _determineInputGroups(
1396 inputPtc,
1397 self.config.doAutoGrouping,
1398 self.config.autoGroupingUseExptime,
1399 self.config.autoGroupingMaxSignalFraction,
1400 self.config.autoGroupingThreshold,
1401 self.config.groupingColumn,
1402 self.config.minPhotodiodeCurrent,
1403 )
1405 if self.config.doAbsoluteSplineFitTemperature:
1406 if self.config.absoluteSplineFitTemperatureColumn not in inputPtc.auxValues:
1407 raise ValueError(
1408 "Config requests fitting temperature coefficient for "
1409 f"{self.config.splineFitTemperatureColumn} but this column "
1410 "is not available in inputPtc.auxValues.",
1411 )
1412 temperatureValues = inputPtc.auxValues[self.config.splineFitTemperatureColumn]
1413 else:
1414 temperatureValues = None
1416 # Fill the linearizer with empty values.
1417 firstAmp = None
1418 for ampName in inputPtc.ampNames:
1419 if ampName not in inputPtc.badAmps:
1420 firstAmp = ampName
1421 break
1422 if firstAmp is None:
1423 raise pipeBase.NoWorkFound("No valid amps in input PTC.")
1424 nExp = len(inputPtc.inputExpIdPairs[firstAmp]) * 2
1425 nAmp = len(inputPtc.ampNames)
1427 for amp in detector:
1428 ampName = amp.getName()
1430 linearizer.linearityType[ampName] = "None"
1431 linearizer.linearityCoeffs[ampName] = np.zeros(1)
1432 # This is not used; kept for compatibility.
1433 linearizer.linearityBBox[ampName] = amp.getBBox()
1434 linearizer.fitParams[ampName] = np.zeros(1)
1435 linearizer.fitParamsErr[ampName] = np.zeros(1)
1436 linearizer.fitChiSq[ampName] = np.nan
1437 linearizer.fitResiduals[ampName] = np.zeros(nExp)
1438 linearizer.fitResidualsSigmaMad[ampName] = np.nan
1439 linearizer.fitResidualsUnmasked[ampName] = np.zeros(nExp)
1440 linearizer.fitResidualsModel[ampName] = np.zeros(nExp)
1441 linearizer.linearFit[ampName] = np.zeros(2)
1442 linearizer.linearityTurnoff[ampName] = np.nan
1443 linearizer.linearityMaxSignal[ampName] = np.nan
1444 linearizer.inputMask[ampName] = np.zeros(nExp, dtype=np.bool_)
1445 linearizer.inputAbscissa[ampName] = np.zeros(nExp)
1446 linearizer.inputOrdinate[ampName] = np.zeros(nExp)
1447 linearizer.inputGroupingIndex[ampName] = np.zeros(nExp, dtype=np.int64)
1448 linearizer.inputNormalization[ampName] = np.ones(nExp)
1450 linearizer.absoluteReferenceAmplifier = ""
1452 # Extract values in common, and per-amp.
1453 data = np.zeros(
1454 nExp,
1455 dtype=[
1456 ("exp_id", "i8"),
1457 ("exptime", "f8"),
1458 ("photocharge", "f8"),
1459 ("mjd", "f8"),
1460 ("raw_mean", ("f8", nAmp)),
1461 ("abscissa", "f8"),
1462 ("grouping", "i4"),
1463 # The following are computed in the relative scaling
1464 # measurements.
1465 ("ref_counts", "f8"),
1466 ("gain_ratio", ("f8", nAmp)),
1467 ],
1468 )
1470 data["exp_id"] = np.asarray(inputPtc.inputExpIdPairs[firstAmp]).ravel()
1471 data["exptime"] = np.repeat(inputPtc.rawExpTimes[firstAmp], 2)
1472 data["mjd"] = np.repeat(inputPtc.inputExpPairMjdStartList[firstAmp], 2)
1473 data["photocharge"] = np.repeat(inputPtc.photoCharges[firstAmp], 2)
1474 data["photocharge"][::2] -= inputPtc.photoChargeDeltas[firstAmp] / 2.
1475 data["photocharge"][1::2] += inputPtc.photoChargeDeltas[firstAmp] / 2.
1476 data["grouping"] = np.repeat(groupingValues, 2)
1478 for i, amp in enumerate(detector):
1479 ampName = amp.getName()
1481 data["raw_mean"][:, i] = np.repeat(inputPtc.rawMeans[ampName], 2)
1482 data["raw_mean"][::2, i] -= inputPtc.rawDeltas[ampName] / 2.
1483 data["raw_mean"][1::2, i] += inputPtc.rawDeltas[ampName] / 2.
1485 if self.config.usePhotodiode:
1486 data["abscissa"][:] = data["photocharge"]
1488 data["abscissa"][data["photocharge"] < self.config.minPhotodiodeCurrent] = np.nan
1489 else:
1490 data["abscissa"][:] = data["exptime"]
1492 # Normalize if configured.
1493 inputNorm = np.ones(nExp, dtype=np.float64)
1494 if self.config.useFocalPlaneNormalization:
1495 a, b = esutil.numpy_util.match(data["exp_id"], inputNormalization["exposure"])
1496 inputNorm[a] = inputNormalization["normalization"][b]
1497 data["abscissa"] *= inputNorm
1499 postTurnoffMasks = {}
1501 # Compute linearity turnoff for each amp.
1502 for i, amp in enumerate(detector):
1503 ampName = amp.getName()
1505 if ampName in inputPtc.badAmps:
1506 self.log.warning(
1507 "Amp %s in detector %s has no usable PTC information. Skipping!",
1508 ampName,
1509 detector.getName(),
1510 )
1511 continue
1513 mask = np.isfinite(data["raw_mean"][:, i])
1515 turnoffMask = np.repeat(inputPtc.expIdMask[ampName], 2)
1516 turnoffMask &= mask
1518 _, turnoff, maxSignal, goodPoints = _computeTurnoffAndMax(
1519 data["abscissa"],
1520 data["raw_mean"][:, i],
1521 turnoffMask,
1522 data["grouping"],
1523 ampName=ampName,
1524 minSignalFitLinearityTurnoff=self.config.minSignalFitLinearityTurnoff,
1525 maxFracLinearityDeviation=self.config.maxFracLinearityDeviation,
1526 log=self.log,
1527 )
1529 # Use the goodPoints as an initial estimate of the mask
1530 # above the turnoff. But we only want to maintain the
1531 # "high end" outliers.
1532 postTurnoffMask = goodPoints
1533 postTurnoffMask[data["raw_mean"][:, i] < np.median(data["raw_mean"][:, i])] = True
1534 postTurnoffMasks[ampName] = postTurnoffMask
1536 if np.isnan(turnoff):
1537 # This is a bad amp, with no linearizer.
1538 self.log.warning(
1539 "Amp %s in detector %s has no usable linearizer information. Skipping!",
1540 ampName,
1541 detector.getName(),
1542 )
1543 continue
1545 linearizer.linearityTurnoff[ampName] = turnoff
1546 linearizer.linearityMaxSignal[ampName] = maxSignal
1548 self.log.info("Amplifier %s has a linearity turnoff of %.2f adu.", ampName, turnoff)
1550 # Choose the reference amplifier as the one with the largest
1551 # turnoff. This ensures that the absolute fit covers the full
1552 # range. We additionally confirm that the ptc turnoff is
1553 # finite for this amplifier.
1554 turnoffArray = np.asarray([linearizer.linearityTurnoff[ampName] for ampName in inputPtc.ampNames])
1555 # This is a possibly redundant check to make sure that a bad amp is
1556 # not chosen as a reference amp.
1557 for i, ampName in enumerate(inputPtc.ampNames):
1558 if ampName in inputPtc.badAmps or not np.isfinite(inputPtc.ptcTurnoff[ampName]):
1559 turnoffArray[i] = np.nan
1561 if np.all(~np.isfinite(turnoffArray)):
1562 # Return the default blank linearizer.
1563 linearizer.hasLinearity = True
1564 linearizer.validate()
1565 linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
1566 linearizer.updateMetadata(setDate=True, setCalibId=True)
1567 linearizer.updateMetadataFromExposures([inputPtc])
1568 provenance = IsrProvenance(calibType='linearizer')
1570 return pipeBase.Struct(
1571 outputLinearizer=linearizer,
1572 outputProvenance=provenance,
1573 )
1575 refAmpIndex = np.argmax(np.nan_to_num(turnoffArray))
1576 refAmpName = inputPtc.ampNames[refAmpIndex]
1577 linearizer.absoluteReferenceAmplifier = refAmpName
1579 # Choose a reference image.
1580 refExpIndex = np.argmin(
1581 np.abs(
1582 np.nan_to_num(data["raw_mean"][:, refAmpIndex]) - self.config.relativeSplineReferenceCounts
1583 )
1584 )
1585 refExpId = data["exp_id"][refExpIndex]
1587 self.log.info(
1588 "Using exposure %d (%.2f adu in amp %s) as reference.",
1589 refExpId,
1590 data["raw_mean"][refExpIndex, refAmpIndex],
1591 refAmpName,
1592 )
1594 # We need to know if the reference exposure is the first or second
1595 # in the pair, because the binned are pairs.
1596 offset = refExpIndex % 2
1598 refBinned = binnedImagesHandleDict[data["exp_id"][refExpIndex - offset]].get()
1599 refBinned = copy.copy(refBinned)
1600 if offset == 0:
1601 refBinned["value"] = refBinned["value1"]
1602 else:
1603 refBinned["value"] = refBinned["value2"]
1605 # Scale reference according to reference amplifier.
1606 refScaling = np.median(refBinned["value"][refBinned["amp_index"] == refAmpIndex])
1607 refBinned["value"] /= refScaling
1609 # Get the invidual amp scalings.
1610 # These are the relative gains for the reference image.
1611 ampScalings = np.asarray(
1612 [
1613 np.median(refBinned["value"][refBinned["amp_index"] == ampIndex])
1614 for ampIndex in range(nAmp)
1615 ],
1616 )
1618 # Compute the gain ratios for every exposure.
1619 # The binned images are stored as pairs.
1620 self.log.info("Computing gain ratios for %d exposures.", len(data))
1621 for i in range(len(data)):
1622 expId = data["exp_id"][i]
1623 if (i % 2) == 0:
1624 binned = binnedImagesHandleDict[expId].get()
1625 binned["value"] = binned["value1"]
1626 else:
1627 binned["value"] = binned["value2"]
1629 binned["value"] /= refBinned["value"]
1631 gainRatios = np.asarray(
1632 [
1633 np.median(binned["value"][binned["amp_index"] == ampIndex])
1634 for ampIndex in range(nAmp)
1635 ]
1636 )
1637 ref_counts = gainRatios[refAmpIndex]
1638 gainRatios /= ref_counts
1640 data["ref_counts"][i] = ref_counts
1641 data["gain_ratio"][i, :] = gainRatios
1643 # We need to know which group has the largest size.
1644 groupAmplitudes = np.zeros(len(np.unique(data["grouping"])))
1645 for g in range(len(groupAmplitudes)):
1646 use = (data["grouping"] == g)
1647 groupAmplitudes[g] = np.nanmax(data["ref_counts"][use]) - np.nanmin(data["ref_counts"][use])
1648 maxAmplitudeGroup = np.argmax(groupAmplitudes)
1650 self.log.info("Illumination group %d has the largest signal amplitude.", maxAmplitudeGroup)
1652 # Compute relative linearization first.
1653 maxRelNodes = 0
1655 for i, amp in enumerate(detector):
1656 if i == refAmpIndex:
1657 continue
1659 ampName = amp.getName()
1661 ptcTurnoff = inputPtc.ptcTurnoff[ampName]
1662 linearityTurnoff = linearizer.linearityTurnoff[ampName]
1664 if not np.isfinite(ptcTurnoff) or not np.isfinite(linearityTurnoff):
1665 # This is a bad amp; skip it.
1666 continue
1668 if ptcTurnoff < self.config.relativeSplineLowThreshold:
1669 lowThreshold = 0.0
1670 else:
1671 lowThreshold = self.config.relativeSplineLowThreshold
1673 # Need to check what happens if these are the same (corner dets).
1674 relNodes = _noderator(
1675 lowThreshold,
1676 ptcTurnoff,
1677 linearityTurnoff,
1678 self.config.relativeSplineMinimumSignalNode,
1679 self.config.relativeSplineLowNodeSize,
1680 self.config.relativeSplineMidNodeSize,
1681 self.config.relativeSplineHighNodeSize,
1682 )
1684 self.log.info("Relative linearity for amplifier %s using %d nodes.", ampName, len(relNodes))
1686 if len(relNodes) > maxRelNodes:
1687 maxRelNodes = len(relNodes)
1689 relAbscissa = data["ref_counts"] * ampScalings[i]
1690 relOrdinate = data["ref_counts"] * data["gain_ratio"][:, i] * ampScalings[i]
1692 # The mask here must exclude everything beyond the turnoff.
1693 relMask = postTurnoffMasks[ampName] & np.isfinite(relAbscissa) & np.isfinite(relOrdinate)
1695 linearizer.inputMask[ampName] = relMask.copy()
1696 linearizer.inputAbscissa[ampName] = relAbscissa.copy()
1697 linearizer.inputOrdinate[ampName] = relOrdinate.copy()
1698 linearizer.inputGroupingIndex[ampName] = data["grouping"].copy()
1699 linearizer.inputNormalization[ampName] = np.ones_like(relAbscissa)
1701 fitter = AstierSplineLinearityFitter(
1702 relNodes,
1703 data["grouping"],
1704 relAbscissa,
1705 relOrdinate,
1706 mask=relMask,
1707 fit_offset=False,
1708 fit_weights=False,
1709 fit_temperature=False,
1710 max_signal_nearly_linear=ptcTurnoff,
1711 fit_temporal=False,
1712 # Turn off max correction clipping.
1713 max_correction=np.inf,
1714 )
1715 p0 = fitter.estimate_p0()
1716 pars = fitter.fit(
1717 p0,
1718 min_iter=self.config.relativeSplineFitMinIter,
1719 max_iter=self.config.relativeSplineFitMaxIter,
1720 max_rejection_per_iteration=self.config.relativeSplineFitMaxRejectionPerIteration,
1721 n_sigma_clip=self.config.relativeNSigmaClipLinear,
1722 )
1724 # Confirm that the first parameter is 0, and set it to
1725 # exactly zero.
1726 relValues = pars[fitter.par_indices["values"]]
1727 if not np.isclose(relValues[0], 0):
1728 raise RuntimeError("Programmer error! First spline parameter must "
1729 "be consistent with zero.")
1730 relValues[0] = 0.0
1732 # We adjust the node values according to the slope of the
1733 # group with the largest amplitude. This removes a degeneracy
1734 # in the normalization and ensures that the overall linearized
1735 # correction is as close to the reference as possible.
1736 relValues -= (1.0 - pars[fitter.par_indices["groups"][maxAmplitudeGroup]]) * relNodes
1738 relChisq = fitter.compute_chisq_dof(pars)
1740 # Our reference fit is always 1.0 slope.
1741 relLinearFit = np.array([0.0, 1.0])
1743 # Adjust the abscissa for different groups for residuals.
1744 for j, groupIndex in enumerate(fitter.group_indices):
1745 relAbscissa[groupIndex] *= (pars[fitter.par_indices["groups"][j]] / relLinearFit[1])
1747 relMask = fitter.mask
1749 # Record values in the linearizer.
1750 linearizer.linearityType[ampName] = "DoubleSpline"
1751 # Note that we have a placeholder for the number of nodes in
1752 # the absolute spline.
1753 linearizer.linearityCoeffs[ampName] = np.concatenate([[len(relNodes), 0], relNodes, relValues])
1754 linearizer.fitChiSq[ampName] = relChisq
1755 linearizer.linearFit[ampName] = relLinearFit
1757 # Compute residuals.
1758 spl = Akima1DInterpolator(relNodes, relValues, method="akima")
1759 relOffset = spl(np.clip(relOrdinate, relNodes[0], relNodes[-1]))
1760 relModel = relOrdinate - relOffset
1762 if relMask.sum() < 2:
1763 self.log.warning("Amp %s in detector %s has not enough points in linear ordinate "
1764 "for residuals. Skipping!", ampName, detector.getName())
1765 relResiduals = np.full_like(relModel, np.nan)
1766 relResidualsUnmasked = relResiduals.copy()
1767 else:
1768 postLinearFit, _, _, _ = irlsFit(
1769 relLinearFit,
1770 relAbscissa[relMask],
1771 relModel[relMask],
1772 funcPolynomial,
1773 )
1774 # When computing residuals, we only care about the slope of
1775 # the postLinearFit and not the intercept. The intercept
1776 # itself depends on a possibly unknown zero in the abscissa
1777 # (often photodiode) which may have an arbitrary value.
1778 relResiduals = relModel - (postLinearFit[1] * relAbscissa)
1779 relResidualsUnmasked = relResiduals.copy()
1780 # We set masked residuals to nan.
1781 relResiduals[~relMask] = np.nan
1783 linearizer.fitResidualsUnmasked[ampName] = relResidualsUnmasked
1784 linearizer.fitResiduals[ampName] = relResiduals
1785 linearizer.fitResidualsModel[ampName] = relModel.copy()
1787 finite = np.isfinite(relResiduals)
1788 if finite.sum() == 0:
1789 sigmad = np.nan
1790 else:
1791 sigmad = median_abs_deviation(relResiduals[finite]/relOrdinate[finite], scale="normal")
1792 linearizer.fitResidualsSigmaMad[ampName] = sigmad
1794 # Now compute absolute linearization.
1796 if temperatureValues is not None:
1797 temperatureValuesScaled = temperatureValues - np.median(temperatureValues)
1798 else:
1799 temperatureValuesScaled = None
1801 if self.config.doAbsoluteSplineFitTemporal:
1802 inputMjdScaled = data["mjd"].copy()
1803 inputMjdScaled -= np.nanmedian(inputMjdScaled)
1804 else:
1805 inputMjdScaled = None
1807 # These are guaranteed to be finite (as checked previously).
1808 absPtcTurnoff = inputPtc.ptcTurnoff[refAmpName]
1809 absLinearityTurnoff = linearizer.linearityTurnoff[refAmpName]
1811 if absPtcTurnoff < self.config.absoluteSplineLowThreshold:
1812 lowThreshold = 0.0
1813 else:
1814 lowThreshold = self.config.absoluteSplineLowThreshold
1816 absNodes = _noderator(
1817 lowThreshold,
1818 absPtcTurnoff,
1819 absLinearityTurnoff,
1820 self.config.absoluteSplineMinimumSignalNode,
1821 self.config.absoluteSplineLowNodeSize,
1822 # The medium and high are matched for absolute spline.
1823 self.config.absoluteSplineNodeSize,
1824 self.config.absoluteSplineNodeSize,
1825 )
1827 self.log.info("Absolute linearity for using %d nodes.", len(absNodes))
1829 absAbscissa = data["abscissa"].copy()
1830 absOrdinate = data["ref_counts"].copy()
1831 # The mask here must exclude everything beyond the turnoff.
1832 absMask = postTurnoffMasks[refAmpName] & np.isfinite(absAbscissa) & np.isfinite(absOrdinate)
1834 # We store the absolute residuals with the reference amplifier.
1835 linearizer.linearityType[refAmpName] = "DoubleSpline"
1836 linearizer.inputMask[refAmpName] = absMask.copy()
1837 linearizer.inputAbscissa[refAmpName] = absAbscissa.copy()
1838 linearizer.inputOrdinate[refAmpName] = absOrdinate.copy()
1839 linearizer.inputGroupingIndex[refAmpName] = data["grouping"].copy()
1840 linearizer.inputNormalization[refAmpName] = inputNorm.copy()
1842 fitter = AstierSplineLinearityFitter(
1843 absNodes,
1844 data["grouping"].copy(),
1845 absAbscissa,
1846 absOrdinate,
1847 mask=absMask,
1848 log=self.log,
1849 fit_offset=self.config.doAbsoluteSplineFitOffset,
1850 fit_weights=self.config.doAbsoluteSplineFitWeights,
1851 weight_pars_start=self.config.absoluteSplineFitWeightParsStart,
1852 fit_temperature=self.config.doAbsoluteSplineFitTemperature,
1853 temperature_scaled=temperatureValuesScaled,
1854 max_signal_nearly_linear=absPtcTurnoff,
1855 fit_temporal=self.config.doAbsoluteSplineFitTemporal,
1856 mjd_scaled=inputMjdScaled,
1857 )
1858 p0 = fitter.estimate_p0()
1859 pars = fitter.fit(
1860 p0,
1861 min_iter=self.config.absoluteSplineFitMinIter,
1862 max_iter=self.config.absoluteSplineFitMaxIter,
1863 max_rejection_per_iteration=self.config.absoluteSplineFitMaxRejectionPerIteration,
1864 n_sigma_clip=self.config.absoluteNSigmaClipLinear,
1865 )
1867 # Confirm that the first parameter is 0, and set it to
1868 # exactly zero.
1869 absValues = pars[fitter.par_indices["values"]]
1870 if not np.isclose(absValues[0], 0):
1871 raise RuntimeError("Programmer error! First spline parameter must "
1872 "be consistent with zero.")
1873 absValues[0] = 0.0
1875 # We need a place to store this.
1876 linearizer.fitChiSq[refAmpName] = fitter.compute_chisq_dof(pars)
1878 absLinearFit = np.array([0.0, np.mean(pars[fitter.par_indices["groups"]])])
1880 # We must modify the inputOrdinate according to the
1881 # nuisance terms in the linearity fit for the residual
1882 # computation code to work properly.
1883 # The true mu (inputOrdinate) is given by
1884 # mu = mu_in * (1 + alpha*t_scale) * (1 + beta*mjd_scale)
1885 if self.config.doAbsoluteSplineFitTemperature:
1886 absOrdinate *= (1.0
1887 + pars[fitter.par_indices["temperature_coeff"]]*temperatureValuesScaled)
1888 if self.config.doAbsoluteSplineFitTemporal:
1889 absOrdinate *= (1.0
1890 + pars[fitter.par_indices["temporal_coeff"]]*inputMjdScaled)
1892 # Adjust the abscissa for different groups for residuals.
1893 for j, groupIndex in enumerate(fitter.group_indices):
1894 absAbscissa[groupIndex] *= (pars[fitter.par_indices["groups"][j]] / absLinearFit[1])
1895 # And remove the offset term.
1896 if self.config.doAbsoluteSplineFitOffset:
1897 absOrdinate -= pars[fitter.par_indices["offset"]]
1899 absMask = fitter.mask
1901 # Compute residuals.
1902 spl = Akima1DInterpolator(absNodes, absValues, method="akima")
1903 absOffset = spl(np.clip(absOrdinate, absNodes[0], absNodes[-1]))
1904 absModel = absOrdinate - absOffset
1906 if absMask.sum() < 2:
1907 self.log.warning("Detector %s has not enough points in linear ordinate "
1908 "for residuals. Skipping!", detector.getName())
1909 # We have to KICK OUT HERE something is VERY wrong.
1910 absResiduals = np.full_like(absModel, np.nan)
1911 absResidualsUnmasked = relResiduals.copy()
1912 else:
1913 postLinearFit, _, _, _ = irlsFit(
1914 absLinearFit,
1915 absAbscissa[absMask],
1916 absModel[absMask],
1917 funcPolynomial,
1918 )
1919 # When computing residuals, we only care about the slope of
1920 # the postLinearFit and not the intercept. The intercept
1921 # itself depends on a possibly unknown zero in the abscissa
1922 # (often photodiode) which may have an arbitrary value.
1923 absResiduals = absModel - (postLinearFit[1] * absAbscissa)
1924 absResidualsUnmasked = absResiduals.copy()
1925 # We set masked residuals to nan.
1926 absResiduals[~absMask] = np.nan
1928 linearizer.fitResidualsUnmasked[refAmpName] = absResidualsUnmasked
1929 linearizer.fitResiduals[refAmpName] = absResiduals
1930 linearizer.fitResidualsModel[refAmpName] = absModel.copy()
1932 finite = np.isfinite(absResiduals)
1933 if finite.sum() == 0:
1934 sigmad = np.nan
1935 else:
1936 sigmad = median_abs_deviation(absResiduals[finite]/absOrdinate[finite], scale="normal")
1937 linearizer.fitResidualsSigmaMad[refAmpName] = sigmad
1939 # Record the absolute nodes and values in each individual amplifier,
1940 # along with extra padding for alignment.
1941 nAbsNodes = len(absNodes)
1942 for i, amp in enumerate(detector):
1943 ampName = amp.getName()
1945 coeffs = np.zeros(2 * nAbsNodes + 2 * maxRelNodes + 2)
1946 if ampName == refAmpName:
1947 # The reference amplifier only has the absolute spline.
1948 coeffs[1] = nAbsNodes
1949 coeffs[2: 2 + 2 * nAbsNodes] = np.concatenate([absNodes, absValues])
1950 else:
1951 nRelNodes = int(linearizer.linearityCoeffs[ampName][0])
1953 coeffs = np.zeros(2 * nAbsNodes + 2 * maxRelNodes + 2)
1954 coeffs[0] = nRelNodes
1955 coeffs[1] = nAbsNodes
1956 relStart = 2
1957 relEnd = relStart + 2 * nRelNodes
1958 coeffs[relStart: relEnd] = linearizer.linearityCoeffs[ampName][relStart: relEnd]
1959 absStart = relEnd
1960 absEnd = absStart + 2 * nAbsNodes
1961 coeffs[absStart: absEnd] = np.concatenate([absNodes, absValues])
1963 linearizer.linearityCoeffs[ampName] = coeffs
1965 linearizer.hasLinearity = True
1966 linearizer.validate()
1967 linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
1968 linearizer.updateMetadata(setDate=True, setCalibId=True)
1969 linearizer.updateMetadataFromExposures([inputPtc])
1970 provenance = IsrProvenance(calibType='linearizer')
1972 return pipeBase.Struct(
1973 outputLinearizer=linearizer,
1974 outputProvenance=provenance,
1975 )
1978def _determineInputGroups(
1979 ptc,
1980 doAutoGrouping,
1981 autoGroupingUseExptime,
1982 autoGroupingMaxSignalFraction,
1983 autoGroupingThreshold,
1984 groupingColumn,
1985 minPhotodiodeCurrent,
1986):
1987 """Determine input groups for linearity fit.
1989 If ``splineGroupingColumn`` is set, then grouping will be done
1990 based on this. Otherwise, if ``doAutoGrouping`` is False, then
1991 no grouping will be done. Finally, grouping will be done by measuring
1992 the ratio of signal to exposure time (if
1993 ``autoGroupingUseExptime`` is set; recommended) or photocharge.
1994 These are then clustered with a simple algorithm to split into groups.
1995 If the data was taking by varying exposure time at different
1996 illumination levels, this grouping is very robust as the clusters are
1997 very well separated.
1999 Parameters
2000 ----------
2001 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
2002 Input PTC to do grouping.
2003 doAutoGrouping : `bool`
2004 Do automatic grouping of pairs?
2005 autoGroupingUseExptime : `bool`
2006 Use exposure time for automatic grouping of pairs?
2007 autoGroupingMaxSignalFraction : `float`
2008 All exposures with signal higher than this threshold will
2009 be put into the largest signal group.
2010 autoGroupingThreshold : `float`
2011 Minimum relative jump from sorted values to determine a group.
2012 minPhotodiodeCurrent : `float`
2013 Minimum photodiode current if auto-grouping is used and
2014 autoGroupingUseExptime is False.
2015 splineGroupingColumn : `str` or `None`
2016 Column to be used for spline grouping (if doAutoGrouping is False).
2018 Returns
2019 -------
2020 groupingValues : `np.ndarray`
2021 Array of values that are unique for a given group.
2022 """
2023 nPair = np.asarray(ptc.inputExpIdPairs[ptc.ampNames[0]]).shape[0]
2024 groupingValues = np.zeros(nPair, dtype=np.int64)
2026 if not doAutoGrouping:
2027 if groupingColumn is not None:
2028 if groupingColumn not in ptc.auxValues:
2029 raise ValueError(f"Config requests grouping by {groupingColumn}, "
2030 "but this column is not available in ptc.auxValues.")
2032 uGroupValues = np.unique(ptc.auxValues[groupingColumn])
2033 for i, uGroupValue in enumerate(uGroupValues):
2034 groupingValues[ptc.auxValues[groupingColumn] == uGroupValue] = i
2035 else:
2036 means = np.zeros((nPair, len(ptc.ampNames)))
2037 exptimes = np.zeros_like(means)
2038 for i, ampName in enumerate(ptc.ampNames):
2039 means[:, i] = ptc.rawMeans[ampName] * ptc.gain[ampName]
2040 exptimes[:, i] = ptc.rawExpTimes[ampName]
2041 detMeans = np.nanmean(means, axis=1)
2042 detExptimes = np.nanmean(exptimes, axis=1)
2044 if autoGroupingUseExptime:
2045 abscissa = detExptimes
2046 else:
2047 abscissa = ptc.photoCharges[ptc.ampNames[0]].copy()
2048 # Set illegal photocharges to NaN.
2049 abscissa[abscissa < minPhotodiodeCurrent] = np.nan
2051 ratio = detMeans / abscissa
2052 ratio /= np.nanmedian(ratio)
2054 # Adjust those that are above threshold so they fall into the
2055 # largest group.
2056 above = (detMeans > autoGroupingMaxSignalFraction*np.nanmax(detMeans))
2057 maxIndex = np.argmax(detMeans[~above])
2058 ratio[above] = ratio[maxIndex]
2060 # The clustering of ratios into groups is performed with a simple
2061 # algorithm based on sorting and looking for the largest gaps.
2062 # See https://stackoverflow.com/a/18385795
2063 st = np.argsort(ratio)
2064 stratio = ratio[st]
2065 delta = stratio[1:] - stratio[0: -1]
2067 transitions, = np.where(delta > autoGroupingThreshold)
2068 # If there are no transitions then everything ends up in group 0.
2069 if len(transitions) > 0:
2070 ratioCuts = stratio[transitions] + autoGroupingThreshold/2.
2072 for i in range(len(transitions) + 1):
2073 if i == 0:
2074 inGroup, = np.where(ratio < ratioCuts[i])
2075 elif i == len(transitions):
2076 inGroup, = np.where(ratio > ratioCuts[i - 1])
2077 else:
2078 inGroup, = np.where((ratio > ratioCuts[i - 1]) & (ratio < ratioCuts[i]))
2079 groupingValues[inGroup] = i
2081 # Ensure out-of-range photoCharges/exptimes are in their own group.
2082 # These are masked later on.
2083 groupingValues[~np.isfinite(abscissa)] = -1
2084 # And put the high ones in the max group.
2085 groupingValues[above] = groupingValues[maxIndex]
2087 return groupingValues
2090def _computeTurnoffAndMax(
2091 abscissa,
2092 ordinate,
2093 initialMask,
2094 groupingValues,
2095 ampName="UNKNOWN",
2096 minSignalFitLinearityTurnoff=1000.0,
2097 maxFracLinearityDeviation=0.01,
2098 log=None,
2099):
2100 """Compute the turnoff and max signal.
2102 Parameters
2103 ----------
2104 abscissa : `np.ndarray`
2105 Input x values, either photoCharges or exposure times.
2106 These should be cleaned of any non-finite values.
2107 ordinate : `np.ndarray`
2108 Input y values, the raw mean values for the amp.
2109 These should be cleaned of any non-finite values.
2110 initialMask : `np.ndarray`
2111 Mask to use for initial fit (usually from PTC).
2112 groupingValues : `np.ndarray`
2113 Array of values that are used to group different fits.
2114 ampName : `str`, optional
2115 Amplifier name (used for logging).
2116 minSignalFitLinearityTurnoff : `float`, optional
2117 Minimum signal to cmpute raw linearity slope for linearityTurnoff.
2118 maxFracLinearityDeviation : `float`, optional
2119 Maximum fraction deviation from raw linearity to compute turnoff.
2120 log : `logging.Logger`, optional
2121 Log object.
2123 Returns
2124 -------
2125 turnoffIndex : `int`
2126 Fit turnoff index (keyed to raw input).
2127 turnoff : `float`
2128 Fit turnoff value.
2129 maxSignal : `float`
2130 Fit maximum signal value.
2131 goodPoints : `np.ndarray`
2132 Mask of good points used in turnoff fit.
2133 """
2134 if log is None:
2135 log = logging.getLogger(__name__)
2137 # Follow eo_pipe:
2138 # https://github.com/lsst-camera-dh/eo_pipe/blob/6afa546569f622b8d604921e248200481c445730/python/lsst/eo/pipe/linearityPlotsTask.py#L50
2139 # Replacing flux with abscissa, Ne with ordinate.
2141 # Fit a line with the y-intercept fixed to zero, using the
2142 # signal counts Ne as the variance in the chi-square, i.e.,
2143 # chi2 = sum( (ordinate - aa*abscissa)**2/ordinate )
2144 # Minimizing chi2 wrt aa, gives
2145 # aa = sum(abscissa) / sum(abscissa**2/ordinate)
2147 fitMask = initialMask.copy()
2148 fitMask[ordinate < minSignalFitLinearityTurnoff] = False
2149 goodPoints = fitMask.copy()
2151 gValues = np.unique(groupingValues)
2152 groupIndicesList = []
2153 for gValue in gValues:
2154 use, = np.where(groupingValues == gValue)
2155 groupIndicesList.append(use)
2157 found = False
2158 while (fitMask.sum() >= 4) and not found:
2159 residuals = np.zeros_like(ordinate)
2161 abscissaMasked = abscissa.copy()
2162 abscissaMasked[~fitMask] = np.nan
2163 ordinateMasked = ordinate.copy()
2164 ordinateMasked[~fitMask] = np.nan
2166 for groupIndices in groupIndicesList:
2167 num = np.nansum(abscissaMasked[groupIndices])
2168 denom = np.nansum(abscissaMasked[groupIndices]**2./ordinateMasked[groupIndices])
2169 aa = num / denom
2171 residuals[groupIndices] = (ordinate[groupIndices] - aa*abscissa[groupIndices]) / \
2172 ordinate[groupIndices]
2174 # Use the residuals to compute the turnoff.
2175 residuals -= np.nanmedian(residuals)
2177 goodPoints = np.abs(residuals) < maxFracLinearityDeviation
2179 if goodPoints.sum() > 4:
2180 # This was an adequate fit.
2181 found = True
2182 turnoff = np.max(ordinate[goodPoints])
2183 turnoffIndex = np.where(np.isclose(ordinate, turnoff))[0][0]
2184 else:
2185 # This was a bad fit; remove the largest outlier.
2186 badIndex = np.argmax(np.abs(residuals)[fitMask])
2187 fitIndices, = np.nonzero(fitMask)
2188 fitMask[fitIndices[badIndex]] = False
2190 if not found:
2191 # Could not find any reasonable value.
2192 log.warning(
2193 "Could not find a reasonable initial linear fit to compute linearity turnoff for "
2194 "amplifier %s; may need finer sampling of input data?",
2195 ampName,
2196 )
2197 if np.all(~fitMask):
2198 return -1, np.nan, np.nan, goodPoints
2200 turnoff = np.max(ordinate[fitMask])
2201 turnoffIndex = np.where(np.isclose(ordinate, turnoff))[0][0]
2203 residuals = np.zeros(len(ordinate))
2205 # Fit the maximum signal.
2206 if turnoffIndex == (len(residuals) - 1):
2207 # This is the last point; we can't do a fit.
2208 # This is not a warning because we do not actually need this
2209 # value in practice.
2210 log.info(
2211 "No linearity turnoff detected for amplifier %s; try to increase the signal range.",
2212 ampName,
2213 )
2214 maxSignal = ordinate[turnoffIndex]
2215 else:
2216 maxSignalInitial = np.nanmax(ordinate)
2218 highFluxPoints = (np.nan_to_num(ordinate)
2219 > (1.0 - maxFracLinearityDeviation)*maxSignalInitial)
2220 maxSignal = np.median(ordinate[highFluxPoints])
2222 return turnoffIndex, turnoff, maxSignal, goodPoints
2225def _noderator(turnoff0, turnoff1, turnoff2, minNode, lowNodeSize, midNodeSize, highNodeSize):
2226 """The "noderator" node-finder.
2228 Parameters
2229 ----------
2230 turnoff0 : `float`
2231 Zeroth turnoff value (e.g. expectation of low-level
2232 non-linearity threshold) (adu).
2233 turnoff1 : `float`
2234 First turnoff value (e.g. ptc turnoff) (adu).
2235 turnoff2 : `float`
2236 Second turnoff value (e.g. linearity turnoff) (adu).
2237 minNode : `float`
2238 Location to place the first node after 0.0 (if this is <= 0.0
2239 it will be ignored) (adu).
2240 lowNodeSize : `float`
2241 Minimum node size in the low-level non-linearity regime
2242 (below turnoff0) (adu).
2243 midNodeSize : `float`
2244 Minimum node size in the mid-level non-linearity regime
2245 (between turnoff0 and turnoff1) (adu).
2246 highNodeSize : `float`
2247 Minimum node size in the high-level non-linearity regime
2248 (between turnoff1 and turnoff2) (adu).
2250 Returns
2251 -------
2252 nodes : `np.ndarray`
2253 Array of node values (adu).
2254 """
2255 if turnoff0 > minNode:
2256 # At least 2 nodes (edges) in the low signal regime.
2257 nNodesLow = np.clip(int(np.ceil((turnoff0 - minNode) / lowNodeSize)), 2, None)
2258 midStart = turnoff0
2259 else:
2260 nNodesLow = 0
2261 midStart = 0.0
2262 # At least 5 nodes (akima minimum) in the mid signal regime.
2263 nNodesMid = np.clip(int(np.ceil((turnoff1 - midStart) / midNodeSize)), 5, None)
2264 if turnoff2 > turnoff1:
2265 # At least 2 nodes (edges) in the high signal regime.
2266 nNodesHigh = np.clip(int(np.ceil((turnoff2 - turnoff1) / highNodeSize)), 2, None)
2267 else:
2268 nNodesHigh = 0
2269 nodesLow = np.linspace(minNode, turnoff0, nNodesLow)
2270 nodesMid = np.linspace(midStart, turnoff1, nNodesMid)
2271 nodesHigh = np.linspace(turnoff1, turnoff2, nNodesHigh)
2273 # Make sure we do not duplicate nodes when concatenating.
2274 nodeList = []
2275 if nNodesLow > 1:
2276 nodeList.append([0.0])
2277 nodeList.append(nodesLow[:-1])
2278 if nNodesMid > 1:
2279 nodeList.append(nodesMid)
2280 if nNodesHigh > 1:
2281 nodeList.append(nodesHigh[1:])
2282 return np.concatenate(nodeList)