Coverage for python / lsst / cp / pipe / ptc / cpPtcSolve.py: 10%
592 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-26 13:06 -0800
1# This file is part of cp_pipe.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21#
22import copy
23import warnings
24import numpy as np
25from collections import Counter
26from itertools import groupby
27from operator import itemgetter
28from scipy.signal import fftconvolve
29from scipy.optimize import least_squares
30from scipy.stats import median_abs_deviation
32import lsst.pex.config as pexConfig
33import lsst.pipe.base as pipeBase
34import lsst.pipe.base.connectionTypes as cT
35from lsst.ip.isr import PhotonTransferCurveDataset
36from lsst.cp.pipe.utils import (fitLeastSq, fitBootstrap,
37 funcAstier, funcAstierWithRolloff,
38 symmetrize, Pol2D, ampOffsetGainRatioFixup)
40from deprecated.sphinx import deprecated
42__all__ = ['PhotonTransferCurveSolveConfig', 'PhotonTransferCurveSolveTask']
45class PhotonTransferCurveSolveConnections(pipeBase.PipelineTaskConnections,
46 dimensions=("instrument", "detector")):
47 inputCovariances = cT.Input(
48 name="ptcCovariances",
49 doc="Tuple with measured covariances from flats.",
50 storageClass="PhotonTransferCurveDataset",
51 dimensions=("instrument", "exposure", "detector"),
52 isCalibration=True,
53 multiple=True,
54 )
55 camera = cT.PrerequisiteInput(
56 name="camera",
57 doc="Camera the input data comes from.",
58 storageClass="Camera",
59 dimensions=("instrument",),
60 isCalibration=True,
61 )
62 outputPtcDataset = cT.Output(
63 name="ptcDatsetProposal",
64 doc="Output proposed ptc dataset.",
65 storageClass="PhotonTransferCurveDataset",
66 dimensions=("instrument", "detector"),
67 multiple=False,
68 isCalibration=True,
69 )
72class PhotonTransferCurveSolveConfig(pipeBase.PipelineTaskConfig,
73 pipelineConnections=PhotonTransferCurveSolveConnections):
74 """Configuration for fitting measured covariances.
75 """
77 ptcFitType = pexConfig.ChoiceField(
78 dtype=str,
79 doc="Fit PTC to Eq. 16 or Eq. 20 in Astier+19.",
80 default="EXPAPPROXIMATION",
81 allowed={
82 "EXPAPPROXIMATION": "Approximation in Astier+19 (Eq. 16).",
83 "FULLCOVARIANCE_NO_B": "Full covariances model in Astier+19 (Eq. 15)",
84 "FULLCOVARIANCE": "Full covariances model in Astier+19 (Eq. 20)",
85 }
86 )
87 minMeanSignal = pexConfig.DictField(
88 keytype=str,
89 itemtype=float,
90 doc="Minimum values (inclusive) of mean signal (in adu) per amp to use."
91 " The same cut is applied to all amps if this parameter [`dict`] is passed as "
92 " {'ALL_AMPS': value}",
93 default={'ALL_AMPS': 0.0},
94 )
95 maxMeanSignal = pexConfig.DictField(
96 keytype=str,
97 itemtype=float,
98 doc="Maximum values (inclusive) of mean signal (in adu) below which to consider, per amp."
99 " The same cut is applied to all amps if this dictionary is of the form"
100 " {'ALL_AMPS': value}",
101 default={'ALL_AMPS': 1e6},
102 )
103 maximumRangeCovariancesAstier = pexConfig.Field(
104 dtype=int,
105 doc="Maximum range of measured covariances as in Astier+19",
106 default=10,
107 )
108 maximumRangeCovariancesAstierFullCovFit = pexConfig.Field(
109 dtype=int,
110 doc="Maximum range up to where to fit covariances as in Astier+19, "
111 "for the FULLCOVARIANCE model."
112 "This is different from maximumRangeCovariancesAstier."
113 "It should be less or equal than maximumRangeCovariancesAstier."
114 "The number of parameters for this model is "
115 "3*maximumRangeCovariancesAstierFullCovFit^2 + 1, so increase with care "
116 "so that the fit is not too slow.",
117 default=10,
118 )
119 doSubtractLongRangeCovariances = pexConfig.Field(
120 dtype=bool,
121 doc="Subtract long-range covariances before FULLCOVARIANCE fit, "
122 "beyond startLongRangeCovariances?",
123 default=False,
124 )
125 startLongRangeCovariances = pexConfig.Field(
126 dtype=int,
127 doc="If doSubtractLongRangeCovariances is True, subtract covariances "
128 "beyond this range. It should be less than maximumRangeCovariancesAstier. ",
129 default=4,
130 )
131 polyDegLongRangeCovariances = pexConfig.Field(
132 dtype=int,
133 doc="If doSubtractLongRangeCovariances is True, polynomial "
134 "degree to fit data beyond startLongRangeCovariances.",
135 default=1,
136 )
137 sigmaClipFullFitCovariancesAstier = pexConfig.Field(
138 dtype=float,
139 doc="sigma clip for full model fit for FULLCOVARIANCE ptcFitType ",
140 default=5.0,
141 )
142 maxIterFullFitCovariancesAstier = pexConfig.Field(
143 dtype=int,
144 doc="Maximum number of iterations in full model fit for FULLCOVARIANCE ptcFitType",
145 default=3,
146 )
147 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
148 # and legacy turnoff
149 polynomialFitDegree = pexConfig.Field(
150 dtype=int,
151 doc="Degree of polynomial to fit the PTC, when 'ptcFitType'=POLYNOMIAL.",
152 default=3,
153 deprecated="This field is no longer used. Will be removed after v30."
154 )
155 doModelPtcRolloff = pexConfig.Field(
156 dtype=bool,
157 doc="Model the roll-off in the PTC turnoff as a exponential decay.",
158 default=False,
159 )
160 maxPtcRolloffDeviation = pexConfig.Field(
161 dtype=float,
162 doc="Maximum percent difference between the model with saturation rolloff and the "
163 "model without to set the PTC turnoff. Only used if doModelPtcRolloff=True.",
164 default=0.01,
165 )
166 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
167 # and legacy turnoff
168 doLegacyTurnoffSelection = pexConfig.Field(
169 dtype=bool,
170 doc="Use 'legacy' computation for PTC turnoff selection. If set "
171 "to False, then the KS test p-value selection will be used instead.",
172 default=False,
173 deprecated="This option has been deprecated and will be removed after v30.",
174 )
175 sigmaCutPtcOutliers = pexConfig.Field(
176 dtype=float,
177 doc="Sigma cut for outlier rejection in PTC.",
178 default=5.0,
179 )
180 maxIterationsPtcOutliers = pexConfig.RangeField(
181 dtype=int,
182 doc="Maximum number of iterations for outlier rejection in PTC.",
183 default=2,
184 min=0
185 )
186 maxSignalInitialPtcOutlierFit = pexConfig.Field(
187 dtype=float,
188 doc="Maximum signal considered for intial outlier fit. This should be below "
189 "the PTC turnoff to ensure accurate outlier rejection. If "
190 "scaleMaxSignalInitialPtcOutlierFit=True then the units are electrons; "
191 "otherwise adu.",
192 default=50_000.,
193 )
194 maxDeltaInitialPtcOutlierFit = pexConfig.Field(
195 dtype=float,
196 doc="If there are any outliers in the initial fit that have mean greater than "
197 "maxSignalInitialPtcOutlierFit, then no points that have this delta "
198 "mean from the previous ``good`` point are allowed. If "
199 "scaleMaxSignalInitialPtcOutlierFit=True then the units are electrons; "
200 "otherwise adu.",
201 default=9_000.,
202 )
203 extendRollofSearchMaskSizeAdu = pexConfig.Field(
204 dtype=float,
205 doc="The amount to extend the mask beyond the initial PTC turnoff to fit "
206 "the rolloff. Only used if doModelPtcRolloff=True (units: adu). ",
207 default=5_000.,
208 )
209 scaleMaxSignalInitialPtcOutlierFit = pexConfig.Field(
210 dtype=bool,
211 doc="Scale maxSignalInitialPtcOutlierFit and maxDeltaInitialPtcOutlierFit "
212 "by approximate gain? If yes then "
213 "maxSignalInitialPtcOutlierFit and maxDeltaInitialPtcOutlierFit "
214 "are assumed to have units of electrons, otherwise adu.",
215 default=True,
216 )
217 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
218 # and legacy turnoff
219 minVarPivotSearch = pexConfig.Field(
220 dtype=float,
221 doc="The code looks for a pivot signal point after which the variance starts decreasing at high-flux"
222 " to exclude then from the PTC model fit. However, sometimes at low fluxes, the variance"
223 " decreases slightly. Set this variable for the variance value, in adu^2, after which the pivot "
224 " should be sought. Only used if doLegacyTurnoffSelection is True.",
225 default=10000,
226 deprecated="This option has been deprecated and will be removed after v30.",
227 )
228 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
229 # and legacy turnoff
230 consecutivePointsVarDecreases = pexConfig.RangeField(
231 dtype=int,
232 doc="Required number of consecutive points/fluxes in the PTC where the variance "
233 "decreases in order to find a first estimate of the PTC turn-off. "
234 "Only used if doLegacyTurnoffSelection is True.",
235 default=2,
236 min=2,
237 deprecated="This option has been deprecated and will be removed after v30.",
238 )
239 ksTestMinPvalue = pexConfig.Field(
240 dtype=float,
241 doc="Minimum value of the Gaussian histogram KS test p-value to be used in PTC fit. ",
242 default=0.01,
243 )
244 doFitBootstrap = pexConfig.Field(
245 dtype=bool,
246 doc="Use bootstrap for the PTC fit parameters and errors?.",
247 default=False,
248 )
249 binSize = pexConfig.Field(
250 dtype=int,
251 doc="Bin the image by this factor in both dimensions.",
252 default=1,
253 )
254 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
255 # and legacy turnoff
256 doAmpOffsetGainRatioFixup = pexConfig.Field(
257 dtype=bool,
258 doc="Do gain ratio fixup based on amp offsets?",
259 default=False,
260 deprecated="This option has been deprecated and will be removed after v30.",
261 )
262 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
263 # and legacy turnoff
264 ampOffsetGainRatioMinAdu = pexConfig.Field(
265 dtype=float,
266 doc="Minimum number of adu to use for amp offset gain ratio fixup.",
267 default=1000.0,
268 deprecated="This option has been deprecated and will be removed after v30.",
269 )
270 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
271 # and legacy turnoff
272 ampOffsetGainRatioMaxAdu = pexConfig.Field(
273 dtype=float,
274 doc="Maximum number of adu to use for amp offset gain ratio fixup.",
275 default=20000.0,
276 deprecated="This option has been deprecated and will be removed after v30.",
277 )
279 def validate(self):
280 super().validate()
281 fitMatrixSide = self.maximumRangeCovariancesAstierFullCovFit
282 measureMatrixSide = self.maximumRangeCovariancesAstier
283 if self.ptcFitType == "FULLCOVARIANCE":
284 if fitMatrixSide > measureMatrixSide:
285 raise RuntimeError("Covariance fit size %s is larger than"
286 "measurement size %s.",
287 fitMatrixSide, measureMatrixSide)
288 if self.doSubtractLongRangeCovariances:
289 startLag = self.startLongRangeCovariances
290 if measureMatrixSide < startLag:
291 raise RuntimeError("Covariance measure size %s is smaller than long"
292 "-range covariance starting point %s.",
293 measureMatrixSide, startLag)
296class PhotonTransferCurveSolveTask(pipeBase.PipelineTask):
297 """Task to fit the PTC from flat covariances.
299 The first task of the PTC measurement pipeline,
300 ``PhotonTransferCurveMeasureTask`` (and assumed to have been run
301 before this task), produced a list of
302 `~lsst.ip.isr.PhotonTransferCurveDataset` objects. Each dataset
303 contains the mean signal and covariances of the
304 difference image of the flat-field images taken at
305 the same exposure time. The list also contains dummy
306 datasets (with no measurements), whose purpose is to have
307 the input and output dimensions of ``PhotonTransferCurveMeasureTask``
308 match.
310 This task, ``PhotonTransferCurveSolveTask``, assembles the list
311 of individual PTC datasets produced
312 by ``PhotonTransferCurveMeasureTask`` into one single final PTC
313 dataset, discarding the dummy datset as appropiate.
314 The task fits the measured (co)variances to one of three models:
315 any of the models described in equations 16 and 20 of Astier+19 and
316 equation 20 with specifically fixed to 0. These options are referred
317 to as ``EXPAPPROXIMATION``, ``FULLCOVARIANCE``, and ``FULLCOVARIANCE_NO_B``
318 in the configuration options of the task, respectively).
319 Parameters of interest such as the gain and noise are derived
320 from the fits. The ``FULLCOVARIANCE`` model is fitted to the
321 full covariance data (as oppossed to the other two models, which
322 are fit to the variance vs mean measurements only).
324 Astier+19: "The Shape of the Photon Transfer Curve
325 of CCD sensors", arXiv:1905.08677
326 """
328 ConfigClass = PhotonTransferCurveSolveConfig
329 _DefaultName = 'cpPhotonTransferCurveSolve'
331 def runQuantum(self, butlerQC, inputRefs, outputRefs):
332 """Ensure that the input and output dimensions are passed along.
334 Parameters
335 ----------
336 butlerQC : `~lsst.daf.butler.QuantumContext`
337 Butler to operate on.
338 inputRefs : `~lsst.pipe.base.InputQuantizedConnection`
339 Input data refs to load.
340 ouptutRefs : `~lsst.pipe.base.OutputQuantizedConnection`
341 Output data refs to persist.
342 """
343 inputs = butlerQC.get(inputRefs)
344 detId = inputRefs.inputCovariances[0].dataId['detector']
345 outputs = self.run(inputCovariances=inputs['inputCovariances'], camera=inputs['camera'], detId=detId)
346 butlerQC.put(outputs, outputRefs)
348 def run(self, inputCovariances, camera=None, detId=0):
349 """Fit measured covariances to different models.
351 Parameters
352 ----------
353 inputCovariances : `list` [`lsst.ip.isr.PhotonTransferCurveDataset`]
354 List of lsst.ip.isr.PhotonTransferCurveDataset datasets.
355 camera : `lsst.afw.cameraGeom.Camera`, optional
356 Input camera.
357 detId : `int`
358 Detector ID to locate the detector in the camera and
359 populate the `lsst.ip.isr.PhotonTransferCurveDataset`
360 metadata.
361 Returns
362 -------
363 results : `lsst.pipe.base.Struct`
364 The resultins structure contains:
366 ``outputPtcDatset``
367 Final PTC dataset, containing information such as the
368 means, variances, and exposure times
369 (`lsst.ip.isr.PhotonTransferCurveDataset`).
370 """
371 # Find the ampNames from a non-dummy ptc.
372 ampNames = []
373 for partialPtcDataset in inputCovariances:
374 if partialPtcDataset.ptcFitType != 'DUMMY':
375 ampNames = partialPtcDataset.ampNames
376 break
378 # Each amp may have a different min and max adu signal
379 # specified in the config.
380 maxMeanSignalDict = {ampName: 1e6 for ampName in ampNames}
381 minMeanSignalDict = {ampName: 0.0 for ampName in ampNames}
382 for ampName in ampNames:
383 if 'ALL_AMPS' in self.config.maxMeanSignal:
384 maxMeanSignalDict[ampName] = self.config.maxMeanSignal['ALL_AMPS']
385 elif ampName in self.config.maxMeanSignal:
386 maxMeanSignalDict[ampName] = self.config.maxMeanSignal[ampName]
388 if 'ALL_AMPS' in self.config.minMeanSignal:
389 minMeanSignalDict[ampName] = self.config.minMeanSignal['ALL_AMPS']
390 elif ampName in self.config.minMeanSignal:
391 minMeanSignalDict[ampName] = self.config.minMeanSignal[ampName]
393 # Assemble individual PTC datasets into a single PTC dataset.
394 datasetPtc = PhotonTransferCurveDataset(
395 ampNames=ampNames,
396 ptcFitType=self.config.ptcFitType,
397 covMatrixSide=self.config.maximumRangeCovariancesAstier,
398 covMatrixSideFullCovFit=self.config.maximumRangeCovariancesAstierFullCovFit)
400 for partialPtcDataset in inputCovariances:
401 # Ignore dummy datasets
402 if partialPtcDataset.ptcFitType == 'DUMMY':
403 continue
405 # Apply min/max masking to the partial PTC.
406 for ampName in ampNames:
407 rawMean = partialPtcDataset.rawMeans[ampName][0]
408 rawVar = partialPtcDataset.rawVars[ampName][0]
409 if (rawMean <= minMeanSignalDict[ampName]) or (rawMean >= maxMeanSignalDict[ampName]) \
410 or not np.isfinite(rawMean) or not np.isfinite(rawVar):
411 partialPtcDataset.expIdMask[ampName][0] = False
413 kspValue = partialPtcDataset.kspValues[ampName][0]
414 if kspValue < self.config.ksTestMinPvalue:
415 partialPtcDataset.expIdMask[ampName][0] = False
417 # Append to the full PTC.
418 datasetPtc.appendPartialPtc(partialPtcDataset)
420 # Sort arrays that are filled so far in the final dataset by
421 # rawMeans index.
422 # First compute the mean across all the amps to make sure that they are
423 # all sorted the same way.
424 detectorMeans = np.zeros(len(datasetPtc.inputExpIdPairs[ampNames[0]]))
426 for i in range(len(detectorMeans)):
427 arr = np.array([datasetPtc.rawMeans[ampName][i] for ampName in ampNames])
428 good, = (np.isfinite(arr)).nonzero()
429 if good.size == 0:
430 detectorMeans[i] = np.nan
431 else:
432 detectorMeans[i] = np.mean(arr[good])
434 # Make sure this is sorted
435 index = np.argsort(detectorMeans)
436 datasetPtc.sort(index)
438 # Always start with an initial EXPAPPROXIMATION fit
439 # Fit the measured covariances vs mean signal to
440 # the Astier+19 full model (Eq. 20). Before that
441 # do a preliminary fit to the variance (C_00) vs mean
442 # signal (mu) curve using the EXPAPPROXIMATION model
443 # (Eq. 16 in Astier+19) in order to
444 # get the flat pairs that are masked. The
445 # points at these fluxes will also be masked when
446 # calculating the other elements of the covariance
447 # matrix, C_ij, i!=j).
448 tempDatasetPtc = copy.copy(datasetPtc)
449 tempDatasetPtc.ptcFitType = "EXPAPPROXIMATION"
450 tempDatasetPtc = self.fitPtc(tempDatasetPtc, computePtcTurnoff=True)
452 # Model the PTC rolloff
453 if self.config.doModelPtcRolloff:
454 tempDatasetPtc = self.fitPtcRolloff(tempDatasetPtc)
456 # Do the final fit
457 # The turnoff was already computed during
458 # the initial fit.
459 tempDatasetPtc.fitType = self.config.ptcFitType
460 if self.config.ptcFitType == "EXPAPPROXIMATION":
461 # Fit with the exponential approximation
462 datasetPtc = self.fitPtc(tempDatasetPtc, computePtcTurnoff=False)
463 elif self.config.ptcFitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]:
464 # Fit with the full covariance model (with/without b)
465 # This does not change the PTC turnoff.
466 tempDatasetPtc.ptcFitType = self.config.ptcFitType
467 datasetPtc = self.fitDataFullCovariance(tempDatasetPtc)
468 else:
469 # This fit type is not recognized
470 raise RuntimeError(
471 f"Fitting option {self.config.ptcFitType} not one of "
472 "'EXPAPPROXIMATION', 'FULLCOVARIANCE', or 'FULLCOVARIANCE_NO_B'"
473 )
475 # Initial validation of PTC fit.
476 for ampName in ampNames:
477 if ampName in datasetPtc.badAmps:
478 # Do not check bad amps!
479 continue
481 # These may be all nan (esp. in tests) and will be filtered
482 # as appropriate later.
483 with warnings.catch_warnings():
484 warnings.simplefilter("ignore")
485 overscanNoise = np.nanmedian(datasetPtc.noiseList[ampName]) # adu
486 overscanNoise *= datasetPtc.gain[ampName] # now in electron
487 noiseFitted = datasetPtc.noise[ampName] # electron
489 # Check if noise is close to noiseFitted
490 if not np.isclose(noiseFitted, overscanNoise, rtol=0.05, atol=0.0, equal_nan=True):
491 self.log.warning(f"Read noise from PTC fit ({noiseFitted}) is not consistent "
492 f"with read noise measured from overscan ({overscanNoise}) for "
493 f"amplifier {ampName}. Try adjusting the fit range.")
495 # Do amp-offset gain ratio fixup if configured.
496 if self.config.doAmpOffsetGainRatioFixup:
497 ampOffsetGainRatioFixup(
498 datasetPtc,
499 self.config.ampOffsetGainRatioMinAdu,
500 self.config.ampOffsetGainRatioMaxAdu,
501 log=self.log,
502 )
504 if camera:
505 detector = camera[detId]
506 else:
507 detector = None
508 datasetPtc.updateMetadataFromExposures(inputCovariances)
509 datasetPtc.updateMetadata(setDate=True, camera=camera, detector=detector)
511 for ampName in ampNames:
512 if ampName in datasetPtc.badAmps:
513 continue
515 datasetPtc.overscanMedian[ampName] = np.nanmedian(
516 datasetPtc.overscanMedianLevelList[ampName][datasetPtc.expIdMask[ampName]]
517 )
518 datasetPtc.overscanMedianSigma[ampName] = median_abs_deviation(
519 datasetPtc.overscanMedianLevelList[ampName][datasetPtc.expIdMask[ampName]],
520 scale="normal",
521 )
523 return pipeBase.Struct(
524 outputPtcDataset=datasetPtc,
525 )
527 def fitDataFullCovariance(self, dataset):
528 """Fit measured flat covariances to the full model in
529 Astier+19 (Eq. 20).
531 Parameters
532 ----------
533 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
534 The dataset containing information such as the means,
535 (co)variances, and exposure times.
537 Returns
538 -------
539 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
540 This is the same dataset as the input parameter, however,
541 it has been modified to include information such as the
542 fit vectors and the fit parameters. See the class
543 `PhotonTransferCurveDatase`.
545 Notes
546 -----
547 The parameters of the full model for C_ij(mu) ("C_ij" and "mu"
548 in adu^2 and adu, respectively) in Astier+19 (Eq. 20) are:
550 - "a" coefficients (r by r matrix), units: 1/electron
551 - "b" coefficients (r by r matrix), units: 1/electron
552 - noise matrix (r by r matrix), units: electron^2
553 - gain, units: electron/adu
555 "b" appears in Eq. 20 only through the "ab" combination, which
556 is defined in this code as "c=ab".
558 Total number of parameters: #entries(a) + #entries(c) + #entries(noise)
559 + 1. This is equivalent to r^2 + r^2 + r^2 + 1, where "r" is the
560 maximum lag considered for the covariances calculation, and the
561 extra "1" is the gain. If "b" is 0, then "c" is 0, and len(pInit) will
562 have r^2 fewer entries.
563 """
564 matrixSide = dataset.covMatrixSide
565 matrixSideFit = dataset.covMatrixSideFullCovFit
566 lenParams = matrixSideFit*matrixSideFit
568 for ampName in dataset.ampNames:
569 lenInputTimes = len(dataset.rawExpTimes[ampName])
570 # Not used when ptcFitType is 'FULLCOVARIANCE'
571 dataset.ptcFitPars[ampName] = np.array([np.nan])
572 dataset.ptcFitParsError[ampName] = np.array([np.nan])
573 dataset.ptcFitChiSq[ampName] = np.nan
575 if ampName in dataset.badAmps:
576 # Bad amp
577 # Entries need to have proper dimensions so read/write
578 # with astropy.Table works.
579 nanMatrixFit = np.full((matrixSideFit, matrixSideFit), np.nan)
580 listNanMatrix = np.full((lenInputTimes, matrixSide, matrixSide), np.nan)
581 listNanMatrixFit = np.full((lenInputTimes, matrixSideFit, matrixSideFit), np.nan)
582 dataset.covariancesModel[ampName] = listNanMatrixFit
583 dataset.covariancesSqrtWeights[ampName] = listNanMatrix
584 dataset.aMatrix[ampName] = nanMatrixFit
585 dataset.bMatrix[ampName] = nanMatrixFit
586 dataset.noiseMatrix[ampName] = nanMatrixFit
588 dataset.expIdMask[ampName] = np.repeat(False, lenInputTimes)
589 dataset.gain[ampName] = np.nan
590 dataset.gainUnadjusted[ampName] = np.nan
591 dataset.gainErr[ampName] = np.nan
592 dataset.noise[ampName] = np.nan
593 dataset.noiseErr[ampName] = np.nan
594 dataset.finalVars[ampName] = np.repeat(np.nan, lenInputTimes)
595 dataset.finalModelVars[ampName] = np.repeat(np.nan, lenInputTimes)
596 dataset.finalMeans[ampName] = np.repeat(np.nan, lenInputTimes)
597 continue
599 muAtAmp = dataset.rawMeans[ampName].copy()
600 maskAtAmp = dataset.expIdMask[ampName]
601 if len(maskAtAmp) == 0:
602 maskAtAmp = np.repeat(True, len(muAtAmp))
604 if np.sum(maskAtAmp) < matrixSide:
605 self.log.warning("Not enough good points to fit PTC for amp %s.", ampName)
607 dataset.badAmps.append(ampName)
608 nanMatrixFit = np.full((matrixSideFit, matrixSideFit), np.nan)
609 listNanMatrix = np.full((lenInputTimes, matrixSide, matrixSide), np.nan)
610 listNanMatrixFit = np.full((lenInputTimes, matrixSideFit, matrixSideFit), np.nan)
611 dataset.covariancesModel[ampName] = listNanMatrixFit
612 dataset.covariancesSqrtWeights[ampName] = listNanMatrix
613 dataset.aMatrix[ampName] = nanMatrixFit
614 dataset.bMatrix[ampName] = nanMatrixFit
615 dataset.noiseMatrix[ampName] = nanMatrixFit
617 dataset.expIdMask[ampName] = np.repeat(False, lenInputTimes)
618 dataset.gain[ampName] = np.nan
619 dataset.gainUnadjusted[ampName] = np.nan
620 dataset.gainErr[ampName] = np.nan
621 dataset.noise[ampName] = np.nan
622 dataset.noiseErr[ampName] = np.nan
623 dataset.finalVars[ampName] = np.repeat(np.nan, lenInputTimes)
624 dataset.finalModelVars[ampName] = np.repeat(np.nan, lenInputTimes)
625 dataset.finalMeans[ampName] = np.repeat(np.nan, lenInputTimes)
626 continue
628 muAtAmpMasked = muAtAmp[maskAtAmp]
629 covAtAmp = dataset.covariances[ampName]
630 covAtAmpMasked = np.nan_to_num(covAtAmp)[maskAtAmp]
631 covSqrtWeightsAtAmp = dataset.covariancesSqrtWeights[ampName]
632 covSqrtWeightsAtAmpMasked = np.nan_to_num(covSqrtWeightsAtAmp)[maskAtAmp]
634 # Subtract long-range covariances
635 if self.config.doSubtractLongRangeCovariances:
636 startLag = self.config.startLongRangeCovariances
637 covAtAmpMasked, covSqrtWeightsAtAmpMasked = self.subtractDistantOffset(
638 muAtAmpMasked, covAtAmpMasked,
639 covSqrtWeightsAtAmpMasked,
640 start=startLag,
641 degree=self.config.polyDegLongRangeCovariances)
643 # In principle, we could fit to a lag smaller than the measured
644 # covariances.
645 r = self.config.maximumRangeCovariancesAstierFullCovFit
646 covAtAmpForFitMasked = covAtAmpMasked[:, :r, :r]
647 covSqrtWeightsAtAmpForFitMasked = covSqrtWeightsAtAmpMasked[:, :r, :r]
649 # Initial fit, to approximate parameters, with c=0
650 a0, c0, noiseMatrix0, gain0 = self.initialFitFullCovariance(
651 muAtAmpMasked,
652 covAtAmpForFitMasked,
653 covSqrtWeightsAtAmpForFitMasked
654 )
656 # Fit full model (Eq. 20 of Astier+19)
657 pInit = np.concatenate((a0.ravel(), c0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None)
659 # Initialize empty results dictionary
660 fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []}
662 # Pick the correct full covariance model function
663 ptcModel = self.funcFullCovarianceModel
664 if dataset.ptcFitType == "FULLCOVARIANCE_NO_B":
665 ptcModel = self.funcFullCovarianceModelNoB
666 pInit = np.concatenate((a0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None)
668 params, paramsErr, _ = fitLeastSq(
669 pInit,
670 muAtAmpMasked,
671 covAtAmpForFitMasked.ravel(),
672 ptcModel,
673 weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(),
674 )
676 if dataset.ptcFitType == "FULLCOVARIANCE_NO_B":
677 zeros = np.zeros_like(params[:lenParams])
678 params = np.insert(params, lenParams, zeros)
679 paramsErr = np.insert(paramsErr, lenParams, zeros)
681 a = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
682 c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
683 noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
684 gain = params[-1]
686 fitResults['a'] = a
687 fitResults['c'] = c
688 fitResults['noiseMatrix'] = noiseMatrix
689 fitResults['gain'] = gain
690 fitResults['paramsErr'] = paramsErr
692 # Put the information in the PTC dataset
694 # Not used when ptcFitType is 'FULLCOVARIANCE*'
695 dataset.ptcFitPars[ampName] = np.array([np.nan])
696 dataset.ptcFitParsError[ampName] = np.array([np.nan])
697 dataset.ptcFitChiSq[ampName] = np.nan
699 # Save full covariances, covariances models, and their weights.
700 # dataset.expIdMask is already full, but needs to be
701 # converted to bool.
702 dataset.expIdMask[ampName] = np.array(dataset.expIdMask[ampName], dtype=bool)
703 dataset.covariances[ampName] = covAtAmp
704 # We evaluate the covariance model everywhere, even the
705 # masked amps.
706 dataset.covariancesModel[ampName] = self.evalCovModel(
707 muAtAmp,
708 fitResults['a'],
709 fitResults['c'],
710 fitResults['noiseMatrix'],
711 fitResults['gain'],
712 setBtoZero=(dataset.ptcFitType == "FULLCOVARIANCE_NO_B"),
713 )
714 dataset.covariancesSqrtWeights[ampName] = covSqrtWeightsAtAmp
715 dataset.aMatrix[ampName] = fitResults['a']
716 dataset.bMatrix[ampName] = fitResults['c']/fitResults['a']
717 dataset.gain[ampName] = fitResults['gain']
718 dataset.gainUnadjusted[ampName] = fitResults['gain']
719 dataset.gainErr[ampName] = fitResults['paramsErr'][-1]
720 readoutNoiseSquared = fitResults['noiseMatrix'][0][0]
721 readoutNoise = np.sqrt(np.fabs(readoutNoiseSquared))
722 dataset.noise[ampName] = readoutNoise
723 readoutNoiseSquaredSigma = fitResults['paramsErr'][2*lenParams]
724 noiseErr = 0.5*(readoutNoiseSquaredSigma/np.fabs(readoutNoiseSquared))*readoutNoise
725 dataset.noiseErr[ampName] = noiseErr
726 dataset.noiseMatrix[ampName] = fitResults['noiseMatrix']
728 dataset.finalVars[ampName] = covAtAmp[:, 0, 0].copy()
729 dataset.finalVars[ampName][~maskAtAmp] = np.nan
730 dataset.finalModelVars[ampName] = dataset.covariancesModel[ampName][:, 0, 0].copy()
731 dataset.finalModelVars[ampName][~maskAtAmp] = np.nan
732 dataset.finalMeans[ampName] = muAtAmp.copy()
733 dataset.finalMeans[ampName][~maskAtAmp] = np.nan
735 return dataset
737 def initialFitFullCovariance(self, mu, cov, sqrtW):
738 """ Performs a crude parabolic fit of the data in order to start
739 the full fit close to the solution, setting b=0 (c=0) in Eq. 20
740 of Astier+19.
742 Parameters
743 ----------
744 mu : `numpy.array`, (N,)
745 Signal `mu` (adu)
746 cov : `numpy.array`, (N, M, M)
747 Covariance arrays of size `(M, M)` (with
748 `M = config.maximumRangeCovariancesAstier`),
749 indexed by mean signal `mu`.
750 sqrtW : `numpy.array`, (N,)
751 Covariance weights, defined as 1./sqrt(Variances)
753 Returns
754 -------
755 a : `numpy.array`, (M, M)
756 "a" parameter per flux in Eq. 20 of Astier+19
757 (units: 1/electron).
758 c : `numpy.array`, (M, M)
759 "c"="ab" parameter per flux in Eq. 20 of Astier+19.
760 (units: 1/electron^2).
761 noiseMatrix : `numpy.array`, (M, M)
762 "noise" parameter per flux in Eq. 20 of Astier+19.
763 (units: electron^2)
764 gain : `float`
765 Amplifier gain (electron/adu)
766 """
767 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
769 # Initialize fit parameters
770 a = np.zeros((matrixSideFit, matrixSideFit))
771 c = np.zeros((matrixSideFit, matrixSideFit))
772 noiseMatrix = np.zeros((matrixSideFit, matrixSideFit))
773 gain = 1.
775 # iterate the fit to account for higher orders
776 # the chi2 does not necessarily go down, so one could
777 # stop when it increases
778 oldChi2 = 1e30
779 for _ in range(5):
780 model = np.nan_to_num(self.evalCovModel(mu, a, c, noiseMatrix, gain, setBtoZero=True))
781 # loop on lags
782 for i in range(matrixSideFit):
783 for j in range(matrixSideFit):
784 # fit a parabola for a given lag
785 parsFit = np.polyfit(mu, cov[:, i, j] - model[:, i, j],
786 2, w=sqrtW[:, i, j])
787 # model equation (Eq. 20) in Astier+19, with c=a*b=0:
788 a[i, j] += parsFit[0]
789 noiseMatrix[i, j] += parsFit[2]
790 if i + j == 0:
791 gain = 1./(1/gain+parsFit[1])
792 weightedRes = (model - cov)*sqrtW
793 chi2 = (weightedRes.flatten()**2).sum()
794 if chi2 > oldChi2:
795 break
796 oldChi2 = chi2
798 return a, c, noiseMatrix, gain
800 def funcFullCovarianceModel(self, params, x):
801 """Model to fit covariances from flat fields; Equation 20 of
802 Astier+19.
804 Parameters
805 ----------
806 params : `list`
807 Parameters of the model: aMatrix, CMatrix, noiseMatrix,
808 gain (e/adu).
809 x : `numpy.array`, (N,)
810 Signal `mu` (adu)
812 Returns
813 -------
814 y : `numpy.array`, (N,)
815 Covariance matrix.
816 """
817 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
818 lenParams = matrixSideFit*matrixSideFit
819 aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
820 cMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
821 noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
822 gain = params[-1]
824 return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain).flatten()
826 def funcFullCovarianceModelNoB(self, params, x):
827 """Model to fit covariances from flat fields; Equation 20 of
828 Astier+19, with b=0 (equivalent to c=a*b=0 in this code).
830 Parameters
831 ----------
832 params : `list`
833 Parameters of the model: aMatrix, noiseMatrix,
834 gain (e/adu).
835 x : `numpy.array`, (N,)
836 Signal mu (adu)
838 Returns
839 -------
840 y : `numpy.array`, (N,)
841 Covariance matrix.
842 """
843 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
844 lenParams = matrixSideFit*matrixSideFit
845 aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
846 cMatrix = np.zeros_like(aMatrix)
847 noiseMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
848 gain = params[-1]
850 return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=True).flatten()
852 def evalCovModel(self, mu, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=False):
853 """Computes full covariances model (Eq. 20 of Astier+19).
855 Parameters
856 ----------
857 mu : `numpy.array`, (N,)
858 List of mean signals (units: adu)
859 aMatrix : `numpy.array`, (M, M)
860 "a" parameter per flux in Eq. 20 of Astier+19
861 (units: 1/electron)
862 cMatrix : `numpy.array`, (M, M)
863 "c"="ab" parameter per flux in Eq. 20 of Astier+19
864 (units: 1/electron^2)
865 noiseMatrix : `numpy.array`, (M, M)
866 "noise" parameter per flux in Eq. 20 of Astier+19
867 (units: electron^2)
868 gain : `float`
869 Amplifier gain (e/adu)
870 setBtoZero=False : `bool`, optional
871 Set "b" parameter in full model (see Astier+19) to zero.
873 Returns
874 -------
875 covModel : `numpy.array`, (N, M, M)
876 Covariances model.
878 Notes
879 -----
880 By default, computes the covModel for the mu's stored(self.mu).
881 Returns cov[Nmu, M, M]. The variance for the PTC is
882 cov[:, 0, 0]. mu and cov are in adus and adus squared. To use
883 electrons for both, the gain should be set to 1. This routine
884 implements the model in Astier+19 (1905.08677).
885 The parameters of the full model for C_ij(mu) ("C_ij" and "mu"
886 in adu^2 and adu, respectively) in Astier+19 (Eq. 20) are:
888 - "a" coefficients (M by M matrix), units: 1/electron
889 - "b" coefficients (M by M matrix), units: 1/electron
890 - noise matrix (M by M matrix), units: electron^2
891 - gain, units: electron/adu
893 "b" appears in Eq. 20 only through the "ab" combination, which
894 is defined in this code as "c=ab".
895 """
896 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
897 sa = (matrixSideFit, matrixSideFit)
898 # pad a with zeros and symmetrize
899 aEnlarged = np.zeros((int(sa[0]*1.5)+1, int(sa[1]*1.5)+1))
900 aEnlarged[0:sa[0], 0:sa[1]] = aMatrix
901 aSym = symmetrize(aEnlarged)
902 # pad c with zeros and symmetrize
903 cEnlarged = np.zeros((int(sa[0]*1.5)+1, int(sa[1]*1.5)+1))
904 cEnlarged[0:sa[0], 0:sa[1]] = cMatrix
905 cSym = symmetrize(cEnlarged)
906 a2 = fftconvolve(aSym, aSym, mode='same')
907 a3 = fftconvolve(a2, aSym, mode='same')
908 ac = fftconvolve(aSym, cSym, mode='same')
909 (xc, yc) = np.unravel_index(np.abs(aSym).argmax(), a2.shape)
911 a1 = aMatrix[np.newaxis, :, :]
912 a2 = a2[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
913 a3 = a3[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
914 ac = ac[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
915 c1 = cMatrix[np.newaxis, ::]
917 # assumes that mu is 1d
918 bigMu = mu[:, np.newaxis, np.newaxis]*gain
919 # c(=a*b in Astier+19) also has a contribution to the last
920 # term, that is absent for now.
921 if setBtoZero:
922 c1 = np.zeros_like(c1)
923 ac = np.zeros_like(ac)
924 covModel = (bigMu/(gain*gain)*(a1*bigMu+2./3.*(bigMu*bigMu)*(a2 + c1)
925 + (1./3.*a3 + 5./6.*ac)*(bigMu*bigMu*bigMu)) + noiseMatrix[np.newaxis, :, :]/gain**2)
926 # add the Poisson term, and the read out noise (variance)
927 covModel[:, 0, 0] += mu/gain
929 return covModel
931 def subtractDistantOffset(self, muAtAmpMasked, covAtAmpMasked, covSqrtWeightsAtAmpMasked,
932 start, degree=1):
933 """Subtract distant offset from the covariance matrices.
935 Parameters
936 ----------
937 muAtAmpMasked : `numpy.array`
938 Masked mean flux array for a particular amplifier.
939 covAtAmpMasked : `numpy.array`
940 Masked measured covariances for a particular amplifier.
941 covSqrtWeightsAtAmpMasked : `numpy.array`
942 Masked inverse covariance weights for a particular amplifier.
943 start : int, optional
944 The starting index to eliminate the core for the fit.
945 degree : int, optional
946 Degree of the polynomial fit.
948 Returns
949 -------
950 covAtAmpMasked : `numpy.array`
951 Subtracted measured covariances for a particular amplifier.
952 covSqrtWeightsAtAmpMasked : `numpy.array`
953 Masked inverse covariance weights for a particular amplifier.
955 Notes
956 -----
957 Ported from https://gitlab.in2p3.fr/astier/bfptc by P. Astier.
959 This function subtracts a distant offset from the
960 covariance matrices using polynomial fitting. The core
961 of the matrices is eliminated for the fit.
963 The function modifies the internal state of the object, updating the
964 covariance matrices and related attributes.
965 """
966 for k in range(len(muAtAmpMasked)):
967 # Make a copy because it will be altered
968 w = np.copy(covSqrtWeightsAtAmpMasked[k, ...])
969 wShape = w.shape
970 i, j = np.meshgrid(range(wShape[0]), range(wShape[1]), indexing='ij')
972 # Eliminate the core for the fit
973 w[:start, :start] = 0
975 poly = Pol2D(i, j, covAtAmpMasked[k, ...], degree, w=w)
976 back = poly.eval(i, j)
978 covAtAmpMasked[k, ...] -= back
980 return covAtAmpMasked, covSqrtWeightsAtAmpMasked
982 # EXPAPPROXIMATION fit method
983 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
984 # and legacy turnoff
985 @staticmethod
986 @deprecated(reason="POLYNOMIAL PTC fit is no longer supported. Will be removed after v30.",
987 version="v30.0", category=FutureWarning)
988 def _initialParsForPolynomial(order):
989 assert order >= 2
990 pars = np.zeros(order, dtype=float)
991 pars[0] = 10
992 pars[1] = 1
993 pars[2:] = 0.0001
994 return pars
996 @staticmethod
997 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
998 # and legacy turnoff
999 @deprecated(reason="POLYNOMIAL PTC fit is no longer supported. Will be removed after v30.",
1000 version="v30.0", category=FutureWarning)
1001 def _boundsForPolynomial(initialPars, lowers=[], uppers=[]):
1002 if not len(lowers):
1003 lowers = [-np.inf for p in initialPars]
1004 if not len(uppers):
1005 uppers = [np.inf for p in initialPars]
1006 lowers[1] = 0 # no negative gains
1007 return (lowers, uppers)
1009 @staticmethod
1010 def _boundsForAstier(initialPars, lowers=[], uppers=[]):
1011 if not len(lowers):
1012 lowers = [-np.inf for p in initialPars]
1013 if not len(uppers):
1014 uppers = [np.inf for p in initialPars]
1015 return (lowers, uppers)
1017 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
1018 # and legacy turnoff
1019 @staticmethod
1020 @deprecated(reason="This is only used by doLegacyTurnoffSelection, which is deprecated, "
1021 "so this deprecated too. Will be removed after v30.",
1022 version="v30.0", category=FutureWarning)
1023 def _getInitialGoodPoints(means, variances, minVarPivotSearch, consecutivePointsVarDecreases):
1024 """Return a boolean array to mask bad points.
1026 Parameters
1027 ----------
1028 means : `numpy.array`
1029 Input array with mean signal values.
1030 variances : `numpy.array`
1031 Input array with variances at each mean value.
1032 minVarPivotSearch : `float`
1033 The variance (in adu^2), above which, the point
1034 of decreasing variance should be sought.
1035 consecutivePointsVarDecreases : `int`
1036 Required number of consecutive points/fluxes
1037 in the PTC where the variance
1038 decreases in order to find a first
1039 estimate of the PTC turn-off.
1041 Returns
1042 ------
1043 goodPoints : `numpy.array` [`bool`]
1044 Boolean array to select good (`True`) and bad (`False`)
1045 points.
1047 Notes
1048 -----
1049 Eliminate points beyond which the variance decreases.
1050 """
1051 goodPoints = np.ones_like(means, dtype=bool)
1052 # Variances are sorted and should monotonically increase
1053 pivotList = np.where(np.array(np.diff(variances)) < 0)[0]
1054 if len(pivotList) > 0:
1055 # For small values, sometimes the variance decreases slightly
1056 # Only look when var > self.config.minVarPivotSearch
1057 pivotList = [p for p in pivotList if variances[p] > minVarPivotSearch]
1058 # Require that the varince decreases during
1059 # consecutivePointsVarDecreases
1060 # consecutive points. This will give a first
1061 # estimate of the PTC turn-off, which
1062 # may be updated (reduced) further in the code.
1063 if len(pivotList) > 1:
1064 # enumerate(pivotList) creates tuples (index, value), for
1065 # each value in pivotList. The lambda function subtracts
1066 # each value from the index.
1067 # groupby groups elements by equal key value.
1068 for k, g in groupby(enumerate(pivotList), lambda x: x[0]-x[1]):
1069 group = (map(itemgetter(1), g))
1070 # Form groups of consecute values from pivotList
1071 group = list(map(int, group))
1072 # values in pivotList are indices where np.diff(variances)
1073 # is negative, i.e., where the variance starts decreasing.
1074 # Find the first group of consecutive numbers when
1075 # variance decreases.
1076 if len(group) >= consecutivePointsVarDecreases:
1077 pivotIndex = np.min(group)
1078 goodPoints[pivotIndex+1:] = False
1079 break
1081 # Finally, we filter out any infinities or NaNs.
1082 goodPoints[(~np.isfinite(means)) | (~np.isfinite(variances))] = False
1084 return goodPoints
1086 def _makeZeroSafe(self, array, substituteValue=1e-9):
1087 """"""
1088 array = np.array(array)
1089 nBad = Counter(np.ravel(array))[0]
1090 if nBad == 0:
1091 return array
1093 index, = np.where(array == 0)
1094 if len(index):
1095 msg = f"Found {nBad} zeros in array at elements {index}"
1096 self.log.warning(msg)
1098 array[index] = substituteValue
1100 return array
1102 def fitPtcRolloff(self, dataset):
1103 """Fit the photon transfer curve to the
1104 Astier+19 approximation (Eq. 16) with a roll-off model to try
1105 and capture saturation effects.
1107 Parameters
1108 ----------
1109 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1110 The dataset containing the means, variances and
1111 exposure times.
1113 Returns
1114 -------
1115 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1116 This is the same dataset as the input parameter, however,
1117 it has been modified to include information such as the
1118 fit vectors and the fit parameters. See the class
1119 `PhotonTransferCurveDataset`.
1120 """
1121 # For FULLCOVARIANCE model fit
1122 matrixSideFit = dataset.covMatrixSideFullCovFit
1123 nanMatrixFit = np.empty((matrixSideFit, matrixSideFit))
1124 nanMatrixFit[:] = np.nan
1126 def errFunc(p, x, y):
1127 return ptcFunc(p, x) - y
1129 for i, ampName in enumerate(dataset.ampNames):
1130 meanVec = dataset.rawMeans[ampName].copy()
1131 varVec = dataset.rawVars[ampName].copy()
1132 varVec = self._makeZeroSafe(varVec)
1134 # Make sure we have this properly sorted.
1135 goodPoints = dataset.expIdMask[ampName].copy()
1137 if not np.isfinite(dataset.ptcTurnoff[ampName]):
1138 msg = (f"SERIOUS: Initial PTC turnoff is not finite. "
1139 f"Setting {ampName} to BAD.")
1140 self.log.warning(msg)
1141 # Fill entries with NaNs
1142 self.fillBadAmp(dataset, ampName)
1143 continue
1145 mask = goodPoints.copy()
1147 # Compute the extended mask
1148 preTurnoff = dataset.ptcTurnoff[ampName]
1149 turnoffIdx = np.argwhere(meanVec == preTurnoff)[0][0]
1150 varianceAtTurnoff = varVec[turnoffIdx]
1152 # Add points up to some threshold below the variance of the turnoff
1153 pointsToFit = (
1154 (meanVec >= preTurnoff)
1155 & (meanVec < preTurnoff+self.config.extendRollofSearchMaskSizeAdu)
1156 & (varVec <= varianceAtTurnoff)
1157 )
1159 # Retain the original mask below the turnoff
1160 pointsToFit = np.logical_or(mask, pointsToFit)
1161 if np.count_nonzero(pointsToFit) == np.count_nonzero(mask):
1162 self.log.warning("Expanding fit to include saturation, but no points detected above "
1163 f"initial computed PTC turnoff for amp {ampName}.")
1164 dataset.expIdRolloffMask[ampName] = pointsToFit
1165 dataset.ptcRolloff[ampName] = dataset.ptcTurnoff[ampName]
1166 dataset.ptcRolloffError[ampName] = np.nan
1167 dataset.ptcRolloffTau[ampName] = np.nan
1168 dataset.ptcRolloffTauError[ampName] = np.nan
1169 continue
1171 # Fit initialization
1172 ptcFunc = funcAstierWithRolloff
1173 parsIniPtc = dataset.ptcFitPars[ampName] # a00, gain, noise^2
1175 # Estimate initial parameters
1176 estimateTurnoff = preTurnoff
1177 estimateTau = -1200 # This turns out to be a generally good guess!
1178 parsIniPtc = np.append(
1179 parsIniPtc,
1180 [estimateTurnoff, estimateTau],
1181 )
1183 # Set initial bounds
1184 bounds = self._boundsForAstier(
1185 parsIniPtc,
1186 lowers=[-1e-4, 0.1, 0, min(preTurnoff - 20000, 0), -10000],
1187 uppers=[0, 10.0, 2000, preTurnoff, -100],
1188 )
1190 # Clip any input values to the bounds.
1191 # This is needed because we have been given parsIniPtc from an
1192 # unconstrained fit.
1193 for index in range(len(parsIniPtc)):
1194 parsIniPtc[index] = float(np.clip(parsIniPtc[index], bounds[0][index], bounds[1][index]))
1196 # Perform the fit
1197 res = least_squares(
1198 errFunc,
1199 parsIniPtc,
1200 bounds=bounds,
1201 args=(meanVec[pointsToFit], varVec[pointsToFit]),
1202 )
1203 pars = res.x
1204 originalModelPars = pars[:-2]
1205 J = res.jac
1206 cov = np.linalg.inv(J.T @ J)
1207 parErrors = np.sqrt(np.diag(cov))
1208 a00, gain, noiseSquared, rolloff, tau = pars
1209 a00Err, gainErr, noiseSquaredError, rolloffError, tauError = parErrors
1211 if not res.success:
1212 self.log.warning(
1213 "Fit with saturation roll off model did not succeed for amp %s.",
1214 ampName,
1215 )
1216 continue
1217 else:
1218 self.log.info("Fit with saturation roll off model returned estimates: "
1219 "rolloff=%f adu (%f el) and tau=%f",
1220 rolloff, rolloff * gain, tau)
1222 # The PTC turnoff is not immediate, and we can tolerate some of
1223 # the physics in the rolloff, so we define the turnoff as a
1224 # threshold of deviation between the model without the rolloff
1225 # and the model with the rolloff. The rolloff estimate is the
1226 # last acceptable point.
1227 modelWithoutRolloff = funcAstier(originalModelPars, meanVec[pointsToFit])
1228 modelWithRolloff = funcAstierWithRolloff(pars, meanVec[pointsToFit])
1229 residual = np.fabs(modelWithRolloff / modelWithoutRolloff - 1)
1230 acceptablePoints = np.argwhere(residual <= self.config.maxPtcRolloffDeviation)
1231 lastGoodIndex = acceptablePoints[-1]
1232 rolloff = meanVec[pointsToFit][lastGoodIndex][0]
1233 rolloff = min(rolloff, preTurnoff)
1235 # Set the mask to the new mask
1236 newMask = pointsToFit * (meanVec < rolloff)
1237 dataset.expIdMask[ampName] = newMask
1239 # Save the rolloff point
1240 dataset.expIdRolloffMask[ampName] = pointsToFit
1241 dataset.ptcRolloff[ampName] = rolloff
1242 dataset.ptcRolloffError[ampName] = rolloffError
1243 dataset.ptcRolloffTau[ampName] = tau
1244 dataset.ptcRolloffTauError[ampName] = tauError
1246 return dataset
1248 def fitPtc(self, dataset, computePtcTurnoff=True):
1249 """Fit the photon transfer curve to the
1250 Astier+19 approximation (Eq. 16).
1252 Fit the photon transfer curve using the exponential
1253 approximation in Astier+19.
1255 Sigma clipping is performed iteratively for the fit, as
1256 well as an initial clipping of data points that are more
1257 than `config.initialNonLinearityExclusionThreshold` away
1258 from lying on a straight line. This other step is necessary
1259 because the photon transfer curve turns over catastrophically
1260 at very high flux (because saturation
1261 drops the variance to ~0) and these far outliers cause the
1262 initial fit to fail, meaning the sigma cannot be calculated
1263 to perform the sigma-clipping.
1265 If `doModelPtcRolloff` is True, a roll-off model will be added
1266 to the initial fit of the PTC to try and capture saturation
1267 effects. This will only be applied if
1268 `ptcFitType=EXPAPPROXIMATION`.
1270 Parameters
1271 ----------
1272 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1273 The dataset containing the means, variances and
1274 exposure times.
1275 computePtcTurnoff : `bool`
1276 Compute and save the PTC turnoff and PTC turnoff sampling
1277 error. If False, this will preserve the attributes in the
1278 input `dataset`.
1280 Returns
1281 -------
1282 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1283 This is the same dataset as the input parameter, however,
1284 it has been modified to include information such as the
1285 fit vectors and the fit parameters. See the class
1286 `PhotonTransferCurveDatase`.
1288 Raises
1289 ------
1290 RuntimeError
1291 Raised if dataset.ptcFitType is None or empty.
1292 """
1293 # For FULLCOVARIANCE model fit
1294 matrixSideFit = dataset.covMatrixSideFullCovFit
1295 nanMatrixFit = np.empty((matrixSideFit, matrixSideFit))
1296 nanMatrixFit[:] = np.nan
1298 for amp in dataset.ampNames:
1299 lenInputTimes = len(dataset.rawExpTimes[amp])
1300 listNanMatrixFit = np.empty((lenInputTimes, matrixSideFit, matrixSideFit))
1301 listNanMatrixFit[:] = np.nan
1303 dataset.covariancesModel[amp] = listNanMatrixFit
1304 dataset.aMatrix[amp] = nanMatrixFit
1305 dataset.bMatrix[amp] = nanMatrixFit
1306 dataset.noiseMatrix[amp] = nanMatrixFit
1308 def errFunc(p, x, y):
1309 return ptcFunc(p, x) - y
1311 sigmaCutPtcOutliers = self.config.sigmaCutPtcOutliers
1312 maxIterationsPtcOutliers = self.config.maxIterationsPtcOutliers
1314 for i, ampName in enumerate(dataset.ampNames):
1315 # The PTC is sorted by detector means,
1316 # however this might not be the true sorted
1317 # order for any specific amplifier. We must
1318 # therefore make sure that the amplifier
1319 # itself is sorted by the amplifier means.
1320 # At the end we will re-order the mask to
1321 # match the sort by detector means.
1323 # These are sorted by detector means
1324 meanVecRaw = dataset.rawMeans[ampName].copy()
1325 varVecRaw = dataset.rawVars[ampName].copy()
1326 varVecRaw = self._makeZeroSafe(varVecRaw)
1328 # These are sorted by the amplifier means
1329 meanVecSort = np.argsort(meanVecRaw)
1330 meanVecSorted = meanVecRaw[meanVecSort]
1331 varVecSorted = varVecRaw[meanVecSort]
1333 # Check if all points are bad from the 'cpExtractPtcTask'
1334 goodPoints = dataset.expIdMask[ampName].copy()
1335 goodPoints = goodPoints[meanVecSort]
1337 if not goodPoints.any():
1338 msg = ("SERIOUS: All points in goodPoints are bad. "
1339 f"Setting {ampName} to BAD.")
1340 self.log.warning(msg)
1341 # Fill entries with NaNs
1342 self.fillBadAmp(dataset, ampName)
1343 continue
1345 mask = goodPoints.copy()
1347 # Set the fitting function and search region
1348 ptcFunc = funcAstier
1349 parsIniPtc = [-2e-6, 1.5, 25.] # a00, gain, noise^2
1350 if self.config.binSize > 1:
1351 bounds = self._boundsForAstier(parsIniPtc)
1352 else:
1353 bounds = self._boundsForAstier(
1354 parsIniPtc,
1355 lowers=[-1e-4, 0.1, 0],
1356 uppers=[0, 10.0, 2000],
1357 )
1359 # We perform an initial (unweighted) fit of variance vs signal
1360 # (after initial KS test or post-drop selection) to look for
1361 # outliers, particularly at the high-flux end. The initial fit
1362 # is performed only for points that are guaranteed to be below
1363 # the PTC turnoff and then extrapolated to ensure that high
1364 # flux points that have abnormal variance values can be properly
1365 # rejected in this phase without biasing the initial fit.
1366 # This algorithm was initially developed by Seth Digel for
1367 # the EO Testing pipeline.
1368 if self.config.scaleMaxSignalInitialPtcOutlierFit:
1369 approxGain = np.nanmedian(meanVecRaw/varVecRaw)
1370 maxADUInitialPtcOutlierFit = self.config.maxSignalInitialPtcOutlierFit/approxGain
1371 maxDeltaADUInitialPtcOutlierFit = self.config.maxDeltaInitialPtcOutlierFit/approxGain
1372 self.log.info(
1373 "Using approximate gain %.3f and ADU signal cutoff of %.1f "
1374 "and max delta %.1f for amplifier %s",
1375 approxGain,
1376 maxADUInitialPtcOutlierFit,
1377 maxDeltaADUInitialPtcOutlierFit,
1378 ampName,
1379 )
1380 else:
1381 maxADUInitialPtcOutlierFit = self.config.maxSignalInitialPtcOutlierFit
1382 maxDeltaADUInitialPtcOutlierFit = self.config.maxDeltaInitialPtcOutlierFit
1384 if maxIterationsPtcOutliers == 0:
1385 # We are not doing any outlier rejection here, but we do want
1386 # an initial fit.
1387 res = least_squares(
1388 errFunc,
1389 parsIniPtc,
1390 bounds=bounds,
1391 args=(meanVecSorted[mask], varVecSorted[mask]),
1392 )
1393 pars = res.x
1394 newMask = mask.copy()
1395 else:
1396 newMask = (
1397 mask
1398 & (meanVecSorted <= maxADUInitialPtcOutlierFit)
1399 & (varVecSorted <= maxADUInitialPtcOutlierFit)
1400 )
1402 converged = False
1403 count = 0
1404 lastMask = mask.copy()
1405 while count < maxIterationsPtcOutliers:
1406 res = least_squares(
1407 errFunc,
1408 parsIniPtc,
1409 bounds=bounds,
1410 args=(meanVecSorted[newMask], varVecSorted[newMask]),
1411 )
1412 pars = res.x
1414 sigResids = (varVecSorted - ptcFunc(pars, meanVecSorted))/np.sqrt(varVecSorted)
1415 # The new mask includes points where the residuals are
1416 # finite, are less than the cut, and include the original
1417 # mask of known points that should not be used.
1418 newMask = (
1419 np.isfinite(sigResids)
1420 & (np.abs(np.nan_to_num(sigResids)) < sigmaCutPtcOutliers)
1421 & mask
1422 )
1423 # Demand at least 2 points to continue.
1424 if np.count_nonzero(newMask) < 2:
1425 msg = (f"SERIOUS: All points after outlier rejection are bad. "
1426 f"Setting {ampName} to BAD.")
1427 self.log.warning(msg)
1428 # Fill entries with NaNs
1429 self.fillBadAmp(dataset, ampName)
1430 break
1432 self.log.debug(
1433 "Iteration %d: Removed %d points in total for %s.",
1434 count,
1435 np.count_nonzero(mask) - np.count_nonzero(newMask),
1436 ampName,
1437 )
1439 # Loop over all used (True) points. If one of them follows
1440 # a False point, then it must be within
1441 # maxDeltaADUInitialPtcOutlierFit of a True point. If it
1442 # is a large gap, everything above is marked False.
1443 useMask, = np.where(newMask)
1444 for useIndex, usePoint in enumerate(useMask):
1445 if useIndex == 0 or newMask[usePoint - 1]:
1446 # The previous point was good; continue.
1447 continue
1448 deltaADU = meanVecSorted[usePoint] - meanVecSorted[useMask[useIndex - 1]]
1449 if deltaADU < maxDeltaADUInitialPtcOutlierFit:
1450 # This jump is fine; continue.
1451 continue
1453 self.log.info(
1454 "Found a jump of %.2f for amp %s, greater than %.2f. Masking points higher "
1455 "than %.2f ADU.",
1456 deltaADU,
1457 ampName,
1458 maxDeltaADUInitialPtcOutlierFit,
1459 meanVecSorted[useMask[useIndex - 1]],
1460 )
1462 # Mark all further points bad.
1463 newMask[usePoint:] = False
1464 break
1466 # If the mask hasn't changed then break out.
1467 if np.all(newMask == lastMask):
1468 self.log.debug("Convergence at iteration %d; breaking loop for %s.", count, ampName)
1469 converged = True
1470 break
1472 lastMask = newMask.copy()
1474 count += 1
1476 if not converged and maxIterationsPtcOutliers > 0:
1477 self.log.warning(
1478 "Outlier detection was not converged prior to %d iteration for %s",
1479 count,
1480 ampName
1481 )
1483 # Set the mask to the new mask, and reset the sorting.
1484 mask = np.zeros(len(meanVecSort), dtype=np.bool_)
1485 mask[meanVecSort[newMask]] = True
1486 maskSorted = newMask.copy()
1488 if not mask.any():
1489 # We have already filled the bad amp above, so continue.
1490 continue
1492 dataset.expIdMask[ampName] = mask
1494 parsIniPtc = pars
1495 meanVecFinal = meanVecRaw[mask]
1496 varVecFinal = varVecRaw[mask]
1498 if Counter(mask)[False] > 0:
1499 self.log.info("Number of points discarded in PTC of amplifier %s:"
1500 " %d out of %d", ampName, Counter(mask)[False], len(meanVecSorted))
1502 if (len(meanVecFinal) < len(parsIniPtc)):
1503 msg = (f"SERIOUS: Not enough data points ({len(meanVecFinal)}) compared to the number of "
1504 f"parameters of the PTC model({len(parsIniPtc)}). Setting {ampName} to BAD.")
1505 self.log.warning(msg)
1506 # Fill entries with NaNs
1507 self.fillBadAmp(dataset, ampName)
1508 continue
1510 # Fit the PTC.
1511 # The variance of the variance is Var(v)=2*v^2/Npix. This is
1512 # already calculated in `makeCovArray` of CpPtcExtract.
1513 # dataset.covariancesSqrtWeights[ampName][:,0,0]
1514 # has 1/sqrt(Var(v)).
1515 weightsY = dataset.covariancesSqrtWeights[ampName][:, 0, 0][mask]
1516 if self.config.doFitBootstrap:
1517 parsFit, parsFitErr, reducedChiSqPtc = fitBootstrap(parsIniPtc, meanVecFinal,
1518 varVecFinal, ptcFunc,
1519 weightsY=weightsY)
1520 else:
1521 parsFit, parsFitErr, reducedChiSqPtc = fitLeastSq(parsIniPtc, meanVecFinal,
1522 varVecFinal, ptcFunc,
1523 weightsY=weightsY)
1525 # Determine PTC turnoff
1526 if computePtcTurnoff:
1527 dataset.ptcTurnoff[ampName] = np.max(meanVecFinal)
1528 # And compute the ptcTurnoffSamplingError as one half the
1529 # difference between the previous and next point.
1530 lastGoodIndex = np.where(maskSorted)[0][-1]
1531 ptcTurnoffLow = meanVecSorted[lastGoodIndex - 1]
1532 if lastGoodIndex == (len(meanVecSorted) - 1):
1533 # If it's the last index, just use the interval.
1534 ptcTurnoffSamplingError = dataset.ptcTurnoff[ampName] - ptcTurnoffLow
1535 elif not np.isfinite(meanVecSorted[lastGoodIndex + 1]):
1536 # If the next index is not finite, just use the interval.
1537 ptcTurnoffSamplingError = dataset.ptcTurnoff[ampName] - ptcTurnoffLow
1538 else:
1539 ptcTurnoffSamplingError = (meanVecSorted[lastGoodIndex + 1] - ptcTurnoffLow)/2.
1540 dataset.ptcTurnoffSamplingError[ampName] = ptcTurnoffSamplingError
1542 dataset.ptcFitPars[ampName] = parsFit
1543 dataset.ptcFitParsError[ampName] = parsFitErr
1544 dataset.ptcFitChiSq[ampName] = reducedChiSqPtc
1546 dataset.finalVars[ampName] = varVecRaw
1547 dataset.finalVars[ampName][~mask] = np.nan
1548 dataset.finalModelVars[ampName] = ptcFunc(parsFit, meanVecRaw)
1549 dataset.finalModelVars[ampName][~mask] = np.nan
1550 dataset.finalMeans[ampName] = meanVecRaw
1551 dataset.finalMeans[ampName][~mask] = np.nan
1553 ptcGain = parsFit[1]
1554 ptcGainErr = parsFitErr[1]
1555 ptcNoise = np.sqrt(np.fabs(parsFit[2]))
1556 ptcNoiseErr = 0.5*(parsFitErr[2]/np.fabs(parsFit[2]))*np.sqrt(np.fabs(parsFit[2]))
1558 # Save results
1559 dataset.gain[ampName] = ptcGain
1560 dataset.gainUnadjusted[ampName] = ptcGain
1561 dataset.gainErr[ampName] = ptcGainErr
1562 dataset.noise[ampName] = ptcNoise
1563 dataset.noiseErr[ampName] = ptcNoiseErr
1565 if len(dataset.badAmps) == 0:
1566 dataset.badAmps = []
1568 return dataset
1570 def fillBadAmp(self, dataset, ampName):
1571 """Fill the dataset with NaNs if there are not enough
1572 good points.
1574 Parameters
1575 ----------
1576 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1577 The dataset containing the means, variances and
1578 exposure times.
1579 ampName : `str`
1580 Amplifier name.
1581 """
1582 if ampName not in dataset.badAmps:
1583 dataset.badAmps.append(ampName)
1584 dataset.expIdMask[ampName] = np.repeat(False, len(dataset.rawExpTimes[ampName]))
1585 dataset.gain[ampName] = np.nan
1586 dataset.gainUnadjusted[ampName] = np.nan
1587 dataset.gainErr[ampName] = np.nan
1588 dataset.noise[ampName] = np.nan
1589 dataset.noiseErr[ampName] = np.nan
1590 dataset.ptcFitPars[ampName] = (np.repeat(np.nan, 3))
1591 dataset.ptcFitParsError[ampName] = (np.repeat(np.nan, 3))
1592 dataset.ptcFitChiSq[ampName] = np.nan
1593 dataset.ptcTurnoff[ampName] = np.nan
1594 dataset.finalVars[ampName] = np.repeat(np.nan, len(dataset.rawExpTimes[ampName]))
1595 dataset.finalModelVars[ampName] = np.repeat(np.nan, len(dataset.rawExpTimes[ampName]))
1596 dataset.finalMeans[ampName] = np.repeat(np.nan, len(dataset.rawExpTimes[ampName]))
1598 return