Coverage for python / lsst / cp / pipe / ptc / cpPtcSolve.py: 10%
592 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 09:34 +0000
1# This file is part of cp_pipe.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21#
22import copy
23import 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(
510 setDate=True,
511 camera=camera,
512 detector=detector,
513 filterName=inputCovariances[0].metadata['FILTER'],
514 )
516 for ampName in ampNames:
517 if ampName in datasetPtc.badAmps:
518 continue
520 datasetPtc.overscanMedian[ampName] = np.nanmedian(
521 datasetPtc.overscanMedianLevelList[ampName][datasetPtc.expIdMask[ampName]]
522 )
523 datasetPtc.overscanMedianSigma[ampName] = median_abs_deviation(
524 datasetPtc.overscanMedianLevelList[ampName][datasetPtc.expIdMask[ampName]],
525 scale="normal",
526 )
528 return pipeBase.Struct(
529 outputPtcDataset=datasetPtc,
530 )
532 def fitDataFullCovariance(self, dataset):
533 """Fit measured flat covariances to the full model in
534 Astier+19 (Eq. 20).
536 Parameters
537 ----------
538 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
539 The dataset containing information such as the means,
540 (co)variances, and exposure times.
542 Returns
543 -------
544 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
545 This is the same dataset as the input parameter, however,
546 it has been modified to include information such as the
547 fit vectors and the fit parameters. See the class
548 `PhotonTransferCurveDatase`.
550 Notes
551 -----
552 The parameters of the full model for C_ij(mu) ("C_ij" and "mu"
553 in adu^2 and adu, respectively) in Astier+19 (Eq. 20) are:
555 - "a" coefficients (r by r matrix), units: 1/electron
556 - "b" coefficients (r by r matrix), units: 1/electron
557 - noise matrix (r by r matrix), units: electron^2
558 - gain, units: electron/adu
560 "b" appears in Eq. 20 only through the "ab" combination, which
561 is defined in this code as "c=ab".
563 Total number of parameters: #entries(a) + #entries(c) + #entries(noise)
564 + 1. This is equivalent to r^2 + r^2 + r^2 + 1, where "r" is the
565 maximum lag considered for the covariances calculation, and the
566 extra "1" is the gain. If "b" is 0, then "c" is 0, and len(pInit) will
567 have r^2 fewer entries.
568 """
569 matrixSide = dataset.covMatrixSide
570 matrixSideFit = dataset.covMatrixSideFullCovFit
571 lenParams = matrixSideFit*matrixSideFit
573 for ampName in dataset.ampNames:
574 lenInputTimes = len(dataset.rawExpTimes[ampName])
575 # Not used when ptcFitType is 'FULLCOVARIANCE'
576 dataset.ptcFitPars[ampName] = np.array([np.nan])
577 dataset.ptcFitParsError[ampName] = np.array([np.nan])
578 dataset.ptcFitChiSq[ampName] = np.nan
580 if ampName in dataset.badAmps:
581 # Bad amp
582 # Entries need to have proper dimensions so read/write
583 # with astropy.Table works.
584 nanMatrixFit = np.full((matrixSideFit, matrixSideFit), np.nan)
585 listNanMatrix = np.full((lenInputTimes, matrixSide, matrixSide), np.nan)
586 listNanMatrixFit = np.full((lenInputTimes, matrixSideFit, matrixSideFit), np.nan)
587 dataset.covariancesModel[ampName] = listNanMatrixFit
588 dataset.covariancesSqrtWeights[ampName] = listNanMatrix
589 dataset.aMatrix[ampName] = nanMatrixFit
590 dataset.bMatrix[ampName] = nanMatrixFit
591 dataset.noiseMatrix[ampName] = nanMatrixFit
593 dataset.expIdMask[ampName] = np.repeat(False, lenInputTimes)
594 dataset.gain[ampName] = np.nan
595 dataset.gainUnadjusted[ampName] = np.nan
596 dataset.gainErr[ampName] = np.nan
597 dataset.noise[ampName] = np.nan
598 dataset.noiseErr[ampName] = np.nan
599 dataset.finalVars[ampName] = np.repeat(np.nan, lenInputTimes)
600 dataset.finalModelVars[ampName] = np.repeat(np.nan, lenInputTimes)
601 dataset.finalMeans[ampName] = np.repeat(np.nan, lenInputTimes)
602 continue
604 muAtAmp = dataset.rawMeans[ampName].copy()
605 maskAtAmp = dataset.expIdMask[ampName]
606 if len(maskAtAmp) == 0:
607 maskAtAmp = np.repeat(True, len(muAtAmp))
609 if np.sum(maskAtAmp) < matrixSide:
610 self.log.warning("Not enough good points to fit PTC for amp %s.", ampName)
612 dataset.badAmps.append(ampName)
613 nanMatrixFit = np.full((matrixSideFit, matrixSideFit), np.nan)
614 listNanMatrix = np.full((lenInputTimes, matrixSide, matrixSide), np.nan)
615 listNanMatrixFit = np.full((lenInputTimes, matrixSideFit, matrixSideFit), np.nan)
616 dataset.covariancesModel[ampName] = listNanMatrixFit
617 dataset.covariancesSqrtWeights[ampName] = listNanMatrix
618 dataset.aMatrix[ampName] = nanMatrixFit
619 dataset.bMatrix[ampName] = nanMatrixFit
620 dataset.noiseMatrix[ampName] = nanMatrixFit
622 dataset.expIdMask[ampName] = np.repeat(False, lenInputTimes)
623 dataset.gain[ampName] = np.nan
624 dataset.gainUnadjusted[ampName] = np.nan
625 dataset.gainErr[ampName] = np.nan
626 dataset.noise[ampName] = np.nan
627 dataset.noiseErr[ampName] = np.nan
628 dataset.finalVars[ampName] = np.repeat(np.nan, lenInputTimes)
629 dataset.finalModelVars[ampName] = np.repeat(np.nan, lenInputTimes)
630 dataset.finalMeans[ampName] = np.repeat(np.nan, lenInputTimes)
631 continue
633 muAtAmpMasked = muAtAmp[maskAtAmp]
634 covAtAmp = dataset.covariances[ampName]
635 covAtAmpMasked = np.nan_to_num(covAtAmp)[maskAtAmp]
636 covSqrtWeightsAtAmp = dataset.covariancesSqrtWeights[ampName]
637 covSqrtWeightsAtAmpMasked = np.nan_to_num(covSqrtWeightsAtAmp)[maskAtAmp]
639 # Subtract long-range covariances
640 if self.config.doSubtractLongRangeCovariances:
641 startLag = self.config.startLongRangeCovariances
642 covAtAmpMasked, covSqrtWeightsAtAmpMasked = self.subtractDistantOffset(
643 muAtAmpMasked, covAtAmpMasked,
644 covSqrtWeightsAtAmpMasked,
645 start=startLag,
646 degree=self.config.polyDegLongRangeCovariances)
648 # In principle, we could fit to a lag smaller than the measured
649 # covariances.
650 r = self.config.maximumRangeCovariancesAstierFullCovFit
651 covAtAmpForFitMasked = covAtAmpMasked[:, :r, :r]
652 covSqrtWeightsAtAmpForFitMasked = covSqrtWeightsAtAmpMasked[:, :r, :r]
654 # Initial fit, to approximate parameters, with c=0
655 a0, c0, noiseMatrix0, gain0 = self.initialFitFullCovariance(
656 muAtAmpMasked,
657 covAtAmpForFitMasked,
658 covSqrtWeightsAtAmpForFitMasked
659 )
661 # Fit full model (Eq. 20 of Astier+19)
662 pInit = np.concatenate((a0.ravel(), c0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None)
664 # Initialize empty results dictionary
665 fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []}
667 # Pick the correct full covariance model function
668 ptcModel = self.funcFullCovarianceModel
669 if dataset.ptcFitType == "FULLCOVARIANCE_NO_B":
670 ptcModel = self.funcFullCovarianceModelNoB
671 pInit = np.concatenate((a0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None)
673 params, paramsErr, _ = fitLeastSq(
674 pInit,
675 muAtAmpMasked,
676 covAtAmpForFitMasked.ravel(),
677 ptcModel,
678 weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(),
679 )
681 if dataset.ptcFitType == "FULLCOVARIANCE_NO_B":
682 zeros = np.zeros_like(params[:lenParams])
683 params = np.insert(params, lenParams, zeros)
684 paramsErr = np.insert(paramsErr, lenParams, zeros)
686 a = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
687 c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
688 noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
689 gain = params[-1]
691 fitResults['a'] = a
692 fitResults['c'] = c
693 fitResults['noiseMatrix'] = noiseMatrix
694 fitResults['gain'] = gain
695 fitResults['paramsErr'] = paramsErr
697 # Put the information in the PTC dataset
699 # Not used when ptcFitType is 'FULLCOVARIANCE*'
700 dataset.ptcFitPars[ampName] = np.array([np.nan])
701 dataset.ptcFitParsError[ampName] = np.array([np.nan])
702 dataset.ptcFitChiSq[ampName] = np.nan
704 # Save full covariances, covariances models, and their weights.
705 # dataset.expIdMask is already full, but needs to be
706 # converted to bool.
707 dataset.expIdMask[ampName] = np.array(dataset.expIdMask[ampName], dtype=bool)
708 dataset.covariances[ampName] = covAtAmp
709 # We evaluate the covariance model everywhere, even the
710 # masked amps.
711 dataset.covariancesModel[ampName] = self.evalCovModel(
712 muAtAmp,
713 fitResults['a'],
714 fitResults['c'],
715 fitResults['noiseMatrix'],
716 fitResults['gain'],
717 setBtoZero=(dataset.ptcFitType == "FULLCOVARIANCE_NO_B"),
718 )
719 dataset.covariancesSqrtWeights[ampName] = covSqrtWeightsAtAmp
720 dataset.aMatrix[ampName] = fitResults['a']
721 dataset.bMatrix[ampName] = fitResults['c']/fitResults['a']
722 dataset.gain[ampName] = fitResults['gain']
723 dataset.gainUnadjusted[ampName] = fitResults['gain']
724 dataset.gainErr[ampName] = fitResults['paramsErr'][-1]
725 readoutNoiseSquared = fitResults['noiseMatrix'][0][0]
726 readoutNoise = np.sqrt(np.fabs(readoutNoiseSquared))
727 dataset.noise[ampName] = readoutNoise
728 readoutNoiseSquaredSigma = fitResults['paramsErr'][2*lenParams]
729 noiseErr = 0.5*(readoutNoiseSquaredSigma/np.fabs(readoutNoiseSquared))*readoutNoise
730 dataset.noiseErr[ampName] = noiseErr
731 dataset.noiseMatrix[ampName] = fitResults['noiseMatrix']
733 dataset.finalVars[ampName] = covAtAmp[:, 0, 0].copy()
734 dataset.finalVars[ampName][~maskAtAmp] = np.nan
735 dataset.finalModelVars[ampName] = dataset.covariancesModel[ampName][:, 0, 0].copy()
736 dataset.finalModelVars[ampName][~maskAtAmp] = np.nan
737 dataset.finalMeans[ampName] = muAtAmp.copy()
738 dataset.finalMeans[ampName][~maskAtAmp] = np.nan
740 return dataset
742 def initialFitFullCovariance(self, mu, cov, sqrtW):
743 """ Performs a crude parabolic fit of the data in order to start
744 the full fit close to the solution, setting b=0 (c=0) in Eq. 20
745 of Astier+19.
747 Parameters
748 ----------
749 mu : `numpy.array`, (N,)
750 Signal `mu` (adu)
751 cov : `numpy.array`, (N, M, M)
752 Covariance arrays of size `(M, M)` (with
753 `M = config.maximumRangeCovariancesAstier`),
754 indexed by mean signal `mu`.
755 sqrtW : `numpy.array`, (N,)
756 Covariance weights, defined as 1./sqrt(Variances)
758 Returns
759 -------
760 a : `numpy.array`, (M, M)
761 "a" parameter per flux in Eq. 20 of Astier+19
762 (units: 1/electron).
763 c : `numpy.array`, (M, M)
764 "c"="ab" parameter per flux in Eq. 20 of Astier+19.
765 (units: 1/electron^2).
766 noiseMatrix : `numpy.array`, (M, M)
767 "noise" parameter per flux in Eq. 20 of Astier+19.
768 (units: electron^2)
769 gain : `float`
770 Amplifier gain (electron/adu)
771 """
772 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
774 # Initialize fit parameters
775 a = np.zeros((matrixSideFit, matrixSideFit))
776 c = np.zeros((matrixSideFit, matrixSideFit))
777 noiseMatrix = np.zeros((matrixSideFit, matrixSideFit))
778 gain = 1.
780 # iterate the fit to account for higher orders
781 # the chi2 does not necessarily go down, so one could
782 # stop when it increases
783 oldChi2 = 1e30
784 for _ in range(5):
785 model = np.nan_to_num(self.evalCovModel(mu, a, c, noiseMatrix, gain, setBtoZero=True))
786 # loop on lags
787 for i in range(matrixSideFit):
788 for j in range(matrixSideFit):
789 # fit a parabola for a given lag
790 parsFit = np.polyfit(mu, cov[:, i, j] - model[:, i, j],
791 2, w=sqrtW[:, i, j])
792 # model equation (Eq. 20) in Astier+19, with c=a*b=0:
793 a[i, j] += parsFit[0]
794 noiseMatrix[i, j] += parsFit[2]
795 if i + j == 0:
796 gain = 1./(1/gain+parsFit[1])
797 weightedRes = (model - cov)*sqrtW
798 chi2 = (weightedRes.flatten()**2).sum()
799 if chi2 > oldChi2:
800 break
801 oldChi2 = chi2
803 return a, c, noiseMatrix, gain
805 def funcFullCovarianceModel(self, params, x):
806 """Model to fit covariances from flat fields; Equation 20 of
807 Astier+19.
809 Parameters
810 ----------
811 params : `list`
812 Parameters of the model: aMatrix, CMatrix, noiseMatrix,
813 gain (e/adu).
814 x : `numpy.array`, (N,)
815 Signal `mu` (adu)
817 Returns
818 -------
819 y : `numpy.array`, (N,)
820 Covariance matrix.
821 """
822 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
823 lenParams = matrixSideFit*matrixSideFit
824 aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
825 cMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
826 noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
827 gain = params[-1]
829 return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain).flatten()
831 def funcFullCovarianceModelNoB(self, params, x):
832 """Model to fit covariances from flat fields; Equation 20 of
833 Astier+19, with b=0 (equivalent to c=a*b=0 in this code).
835 Parameters
836 ----------
837 params : `list`
838 Parameters of the model: aMatrix, noiseMatrix,
839 gain (e/adu).
840 x : `numpy.array`, (N,)
841 Signal mu (adu)
843 Returns
844 -------
845 y : `numpy.array`, (N,)
846 Covariance matrix.
847 """
848 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
849 lenParams = matrixSideFit*matrixSideFit
850 aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
851 cMatrix = np.zeros_like(aMatrix)
852 noiseMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
853 gain = params[-1]
855 return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=True).flatten()
857 def evalCovModel(self, mu, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=False):
858 """Computes full covariances model (Eq. 20 of Astier+19).
860 Parameters
861 ----------
862 mu : `numpy.array`, (N,)
863 List of mean signals (units: adu)
864 aMatrix : `numpy.array`, (M, M)
865 "a" parameter per flux in Eq. 20 of Astier+19
866 (units: 1/electron)
867 cMatrix : `numpy.array`, (M, M)
868 "c"="ab" parameter per flux in Eq. 20 of Astier+19
869 (units: 1/electron^2)
870 noiseMatrix : `numpy.array`, (M, M)
871 "noise" parameter per flux in Eq. 20 of Astier+19
872 (units: electron^2)
873 gain : `float`
874 Amplifier gain (e/adu)
875 setBtoZero=False : `bool`, optional
876 Set "b" parameter in full model (see Astier+19) to zero.
878 Returns
879 -------
880 covModel : `numpy.array`, (N, M, M)
881 Covariances model.
883 Notes
884 -----
885 By default, computes the covModel for the mu's stored(self.mu).
886 Returns cov[Nmu, M, M]. The variance for the PTC is
887 cov[:, 0, 0]. mu and cov are in adus and adus squared. To use
888 electrons for both, the gain should be set to 1. This routine
889 implements the model in Astier+19 (1905.08677).
890 The parameters of the full model for C_ij(mu) ("C_ij" and "mu"
891 in adu^2 and adu, respectively) in Astier+19 (Eq. 20) are:
893 - "a" coefficients (M by M matrix), units: 1/electron
894 - "b" coefficients (M by M matrix), units: 1/electron
895 - noise matrix (M by M matrix), units: electron^2
896 - gain, units: electron/adu
898 "b" appears in Eq. 20 only through the "ab" combination, which
899 is defined in this code as "c=ab".
900 """
901 matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
902 sa = (matrixSideFit, matrixSideFit)
903 # pad a with zeros and symmetrize
904 aEnlarged = np.zeros((int(sa[0]*1.5)+1, int(sa[1]*1.5)+1))
905 aEnlarged[0:sa[0], 0:sa[1]] = aMatrix
906 aSym = symmetrize(aEnlarged)
907 # pad c with zeros and symmetrize
908 cEnlarged = np.zeros((int(sa[0]*1.5)+1, int(sa[1]*1.5)+1))
909 cEnlarged[0:sa[0], 0:sa[1]] = cMatrix
910 cSym = symmetrize(cEnlarged)
911 a2 = fftconvolve(aSym, aSym, mode='same')
912 a3 = fftconvolve(a2, aSym, mode='same')
913 ac = fftconvolve(aSym, cSym, mode='same')
914 (xc, yc) = np.unravel_index(np.abs(aSym).argmax(), a2.shape)
916 a1 = aMatrix[np.newaxis, :, :]
917 a2 = a2[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
918 a3 = a3[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
919 ac = ac[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
920 c1 = cMatrix[np.newaxis, ::]
922 # assumes that mu is 1d
923 bigMu = mu[:, np.newaxis, np.newaxis]*gain
924 # c(=a*b in Astier+19) also has a contribution to the last
925 # term, that is absent for now.
926 if setBtoZero:
927 c1 = np.zeros_like(c1)
928 ac = np.zeros_like(ac)
929 covModel = (bigMu/(gain*gain)*(a1*bigMu+2./3.*(bigMu*bigMu)*(a2 + c1)
930 + (1./3.*a3 + 5./6.*ac)*(bigMu*bigMu*bigMu)) + noiseMatrix[np.newaxis, :, :]/gain**2)
931 # add the Poisson term, and the read out noise (variance)
932 covModel[:, 0, 0] += mu/gain
934 return covModel
936 def subtractDistantOffset(self, muAtAmpMasked, covAtAmpMasked, covSqrtWeightsAtAmpMasked,
937 start, degree=1):
938 """Subtract distant offset from the covariance matrices.
940 Parameters
941 ----------
942 muAtAmpMasked : `numpy.array`
943 Masked mean flux array for a particular amplifier.
944 covAtAmpMasked : `numpy.array`
945 Masked measured covariances for a particular amplifier.
946 covSqrtWeightsAtAmpMasked : `numpy.array`
947 Masked inverse covariance weights for a particular amplifier.
948 start : int, optional
949 The starting index to eliminate the core for the fit.
950 degree : int, optional
951 Degree of the polynomial fit.
953 Returns
954 -------
955 covAtAmpMasked : `numpy.array`
956 Subtracted measured covariances for a particular amplifier.
957 covSqrtWeightsAtAmpMasked : `numpy.array`
958 Masked inverse covariance weights for a particular amplifier.
960 Notes
961 -----
962 Ported from https://gitlab.in2p3.fr/astier/bfptc by P. Astier.
964 This function subtracts a distant offset from the
965 covariance matrices using polynomial fitting. The core
966 of the matrices is eliminated for the fit.
968 The function modifies the internal state of the object, updating the
969 covariance matrices and related attributes.
970 """
971 for k in range(len(muAtAmpMasked)):
972 # Make a copy because it will be altered
973 w = np.copy(covSqrtWeightsAtAmpMasked[k, ...])
974 wShape = w.shape
975 i, j = np.meshgrid(range(wShape[0]), range(wShape[1]), indexing='ij')
977 # Eliminate the core for the fit
978 w[:start, :start] = 0
980 poly = Pol2D(i, j, covAtAmpMasked[k, ...], degree, w=w)
981 back = poly.eval(i, j)
983 covAtAmpMasked[k, ...] -= back
985 return covAtAmpMasked, covSqrtWeightsAtAmpMasked
987 # EXPAPPROXIMATION fit method
988 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
989 # and legacy turnoff
990 @staticmethod
991 @deprecated(reason="POLYNOMIAL PTC fit is no longer supported. Will be removed after v30.",
992 version="v30.0", category=FutureWarning)
993 def _initialParsForPolynomial(order):
994 assert order >= 2
995 pars = np.zeros(order, dtype=float)
996 pars[0] = 10
997 pars[1] = 1
998 pars[2:] = 0.0001
999 return pars
1001 @staticmethod
1002 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
1003 # and legacy turnoff
1004 @deprecated(reason="POLYNOMIAL PTC fit is no longer supported. Will be removed after v30.",
1005 version="v30.0", category=FutureWarning)
1006 def _boundsForPolynomial(initialPars, lowers=[], uppers=[]):
1007 if not len(lowers):
1008 lowers = [-np.inf for p in initialPars]
1009 if not len(uppers):
1010 uppers = [np.inf for p in initialPars]
1011 lowers[1] = 0 # no negative gains
1012 return (lowers, uppers)
1014 @staticmethod
1015 def _boundsForAstier(initialPars, lowers=[], uppers=[]):
1016 if not len(lowers):
1017 lowers = [-np.inf for p in initialPars]
1018 if not len(uppers):
1019 uppers = [np.inf for p in initialPars]
1020 return (lowers, uppers)
1022 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
1023 # and legacy turnoff
1024 @staticmethod
1025 @deprecated(reason="This is only used by doLegacyTurnoffSelection, which is deprecated, "
1026 "so this deprecated too. Will be removed after v30.",
1027 version="v30.0", category=FutureWarning)
1028 def _getInitialGoodPoints(means, variances, minVarPivotSearch, consecutivePointsVarDecreases):
1029 """Return a boolean array to mask bad points.
1031 Parameters
1032 ----------
1033 means : `numpy.array`
1034 Input array with mean signal values.
1035 variances : `numpy.array`
1036 Input array with variances at each mean value.
1037 minVarPivotSearch : `float`
1038 The variance (in adu^2), above which, the point
1039 of decreasing variance should be sought.
1040 consecutivePointsVarDecreases : `int`
1041 Required number of consecutive points/fluxes
1042 in the PTC where the variance
1043 decreases in order to find a first
1044 estimate of the PTC turn-off.
1046 Returns
1047 ------
1048 goodPoints : `numpy.array` [`bool`]
1049 Boolean array to select good (`True`) and bad (`False`)
1050 points.
1052 Notes
1053 -----
1054 Eliminate points beyond which the variance decreases.
1055 """
1056 goodPoints = np.ones_like(means, dtype=bool)
1057 # Variances are sorted and should monotonically increase
1058 pivotList = np.where(np.array(np.diff(variances)) < 0)[0]
1059 if len(pivotList) > 0:
1060 # For small values, sometimes the variance decreases slightly
1061 # Only look when var > self.config.minVarPivotSearch
1062 pivotList = [p for p in pivotList if variances[p] > minVarPivotSearch]
1063 # Require that the varince decreases during
1064 # consecutivePointsVarDecreases
1065 # consecutive points. This will give a first
1066 # estimate of the PTC turn-off, which
1067 # may be updated (reduced) further in the code.
1068 if len(pivotList) > 1:
1069 # enumerate(pivotList) creates tuples (index, value), for
1070 # each value in pivotList. The lambda function subtracts
1071 # each value from the index.
1072 # groupby groups elements by equal key value.
1073 for k, g in groupby(enumerate(pivotList), lambda x: x[0]-x[1]):
1074 group = (map(itemgetter(1), g))
1075 # Form groups of consecute values from pivotList
1076 group = list(map(int, group))
1077 # values in pivotList are indices where np.diff(variances)
1078 # is negative, i.e., where the variance starts decreasing.
1079 # Find the first group of consecutive numbers when
1080 # variance decreases.
1081 if len(group) >= consecutivePointsVarDecreases:
1082 pivotIndex = np.min(group)
1083 goodPoints[pivotIndex+1:] = False
1084 break
1086 # Finally, we filter out any infinities or NaNs.
1087 goodPoints[(~np.isfinite(means)) | (~np.isfinite(variances))] = False
1089 return goodPoints
1091 def _makeZeroSafe(self, array, substituteValue=1e-9):
1092 """"""
1093 array = np.array(array)
1094 nBad = Counter(np.ravel(array))[0]
1095 if nBad == 0:
1096 return array
1098 index, = np.where(array == 0)
1099 if len(index):
1100 msg = f"Found {nBad} zeros in array at elements {index}"
1101 self.log.warning(msg)
1103 array[index] = substituteValue
1105 return array
1107 def fitPtcRolloff(self, dataset):
1108 """Fit the photon transfer curve to the
1109 Astier+19 approximation (Eq. 16) with a roll-off model to try
1110 and capture saturation effects.
1112 Parameters
1113 ----------
1114 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1115 The dataset containing the means, variances and
1116 exposure times.
1118 Returns
1119 -------
1120 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1121 This is the same dataset as the input parameter, however,
1122 it has been modified to include information such as the
1123 fit vectors and the fit parameters. See the class
1124 `PhotonTransferCurveDataset`.
1125 """
1126 # For FULLCOVARIANCE model fit
1127 matrixSideFit = dataset.covMatrixSideFullCovFit
1128 nanMatrixFit = np.empty((matrixSideFit, matrixSideFit))
1129 nanMatrixFit[:] = np.nan
1131 def errFunc(p, x, y):
1132 return ptcFunc(p, x) - y
1134 for i, ampName in enumerate(dataset.ampNames):
1135 meanVec = dataset.rawMeans[ampName].copy()
1136 varVec = dataset.rawVars[ampName].copy()
1137 varVec = self._makeZeroSafe(varVec)
1139 # Make sure we have this properly sorted.
1140 goodPoints = dataset.expIdMask[ampName].copy()
1142 if not np.isfinite(dataset.ptcTurnoff[ampName]):
1143 msg = (f"SERIOUS: Initial PTC turnoff is not finite. "
1144 f"Setting {ampName} to BAD.")
1145 self.log.warning(msg)
1146 # Fill entries with NaNs
1147 self.fillBadAmp(dataset, ampName)
1148 continue
1150 mask = goodPoints.copy()
1152 # Compute the extended mask
1153 preTurnoff = dataset.ptcTurnoff[ampName]
1154 turnoffIdx = np.argwhere(meanVec == preTurnoff)[0][0]
1155 varianceAtTurnoff = varVec[turnoffIdx]
1157 # Add points up to some threshold below the variance of the turnoff
1158 pointsToFit = (
1159 (meanVec >= preTurnoff)
1160 & (meanVec < preTurnoff+self.config.extendRollofSearchMaskSizeAdu)
1161 & (varVec <= varianceAtTurnoff)
1162 )
1164 # Retain the original mask below the turnoff
1165 pointsToFit = np.logical_or(mask, pointsToFit)
1166 if np.count_nonzero(pointsToFit) == np.count_nonzero(mask):
1167 self.log.warning("Expanding fit to include saturation, but no points detected above "
1168 f"initial computed PTC turnoff for amp {ampName}.")
1169 dataset.expIdRolloffMask[ampName] = pointsToFit
1170 dataset.ptcRolloff[ampName] = dataset.ptcTurnoff[ampName]
1171 dataset.ptcRolloffError[ampName] = np.nan
1172 dataset.ptcRolloffTau[ampName] = np.nan
1173 dataset.ptcRolloffTauError[ampName] = np.nan
1174 continue
1176 # Fit initialization
1177 ptcFunc = funcAstierWithRolloff
1178 parsIniPtc = dataset.ptcFitPars[ampName] # a00, gain, noise^2
1180 # Estimate initial parameters
1181 estimateTurnoff = preTurnoff
1182 estimateTau = -1200 # This turns out to be a generally good guess!
1183 parsIniPtc = np.append(
1184 parsIniPtc,
1185 [estimateTurnoff, estimateTau],
1186 )
1188 # Set initial bounds
1189 bounds = self._boundsForAstier(
1190 parsIniPtc,
1191 lowers=[-1e-4, 0.1, 0, min(preTurnoff - 20000, 0), -10000],
1192 uppers=[0, 10.0, 2000, preTurnoff, -100],
1193 )
1195 # Clip any input values to the bounds.
1196 # This is needed because we have been given parsIniPtc from an
1197 # unconstrained fit.
1198 for index in range(len(parsIniPtc)):
1199 parsIniPtc[index] = float(np.clip(parsIniPtc[index], bounds[0][index], bounds[1][index]))
1201 # Perform the fit
1202 res = least_squares(
1203 errFunc,
1204 parsIniPtc,
1205 bounds=bounds,
1206 args=(meanVec[pointsToFit], varVec[pointsToFit]),
1207 )
1208 pars = res.x
1209 originalModelPars = pars[:-2]
1210 J = res.jac
1211 cov = np.linalg.inv(J.T @ J)
1212 parErrors = np.sqrt(np.diag(cov))
1213 a00, gain, noiseSquared, rolloff, tau = pars
1214 a00Err, gainErr, noiseSquaredError, rolloffError, tauError = parErrors
1216 if not res.success:
1217 self.log.warning(
1218 "Fit with saturation roll off model did not succeed for amp %s.",
1219 ampName,
1220 )
1221 continue
1222 else:
1223 self.log.info("Fit with saturation roll off model returned estimates: "
1224 "rolloff=%f adu (%f el) and tau=%f",
1225 rolloff, rolloff * gain, tau)
1227 # The PTC turnoff is not immediate, and we can tolerate some of
1228 # the physics in the rolloff, so we define the turnoff as a
1229 # threshold of deviation between the model without the rolloff
1230 # and the model with the rolloff. The rolloff estimate is the
1231 # last acceptable point.
1232 modelWithoutRolloff = funcAstier(originalModelPars, meanVec[pointsToFit])
1233 modelWithRolloff = funcAstierWithRolloff(pars, meanVec[pointsToFit])
1234 residual = np.fabs(modelWithRolloff / modelWithoutRolloff - 1)
1235 acceptablePoints = np.argwhere(residual <= self.config.maxPtcRolloffDeviation)
1236 lastGoodIndex = acceptablePoints[-1]
1237 rolloff = meanVec[pointsToFit][lastGoodIndex][0]
1238 rolloff = min(rolloff, preTurnoff)
1240 # Set the mask to the new mask
1241 newMask = pointsToFit * (meanVec < rolloff)
1242 dataset.expIdMask[ampName] = newMask
1244 # Save the rolloff point
1245 dataset.expIdRolloffMask[ampName] = pointsToFit
1246 dataset.ptcRolloff[ampName] = rolloff
1247 dataset.ptcRolloffError[ampName] = rolloffError
1248 dataset.ptcRolloffTau[ampName] = tau
1249 dataset.ptcRolloffTauError[ampName] = tauError
1251 return dataset
1253 def fitPtc(self, dataset, computePtcTurnoff=True):
1254 """Fit the photon transfer curve to the
1255 Astier+19 approximation (Eq. 16).
1257 Fit the photon transfer curve using the exponential
1258 approximation in Astier+19.
1260 Sigma clipping is performed iteratively for the fit, as
1261 well as an initial clipping of data points that are more
1262 than `config.initialNonLinearityExclusionThreshold` away
1263 from lying on a straight line. This other step is necessary
1264 because the photon transfer curve turns over catastrophically
1265 at very high flux (because saturation
1266 drops the variance to ~0) and these far outliers cause the
1267 initial fit to fail, meaning the sigma cannot be calculated
1268 to perform the sigma-clipping.
1270 If `doModelPtcRolloff` is True, a roll-off model will be added
1271 to the initial fit of the PTC to try and capture saturation
1272 effects. This will only be applied if
1273 `ptcFitType=EXPAPPROXIMATION`.
1275 Parameters
1276 ----------
1277 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1278 The dataset containing the means, variances and
1279 exposure times.
1280 computePtcTurnoff : `bool`
1281 Compute and save the PTC turnoff and PTC turnoff sampling
1282 error. If False, this will preserve the attributes in the
1283 input `dataset`.
1285 Returns
1286 -------
1287 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1288 This is the same dataset as the input parameter, however,
1289 it has been modified to include information such as the
1290 fit vectors and the fit parameters. See the class
1291 `PhotonTransferCurveDatase`.
1293 Raises
1294 ------
1295 RuntimeError
1296 Raised if dataset.ptcFitType is None or empty.
1297 """
1298 # For FULLCOVARIANCE model fit
1299 matrixSideFit = dataset.covMatrixSideFullCovFit
1300 nanMatrixFit = np.empty((matrixSideFit, matrixSideFit))
1301 nanMatrixFit[:] = np.nan
1303 for amp in dataset.ampNames:
1304 lenInputTimes = len(dataset.rawExpTimes[amp])
1305 listNanMatrixFit = np.empty((lenInputTimes, matrixSideFit, matrixSideFit))
1306 listNanMatrixFit[:] = np.nan
1308 dataset.covariancesModel[amp] = listNanMatrixFit
1309 dataset.aMatrix[amp] = nanMatrixFit
1310 dataset.bMatrix[amp] = nanMatrixFit
1311 dataset.noiseMatrix[amp] = nanMatrixFit
1313 def errFunc(p, x, y):
1314 return ptcFunc(p, x) - y
1316 sigmaCutPtcOutliers = self.config.sigmaCutPtcOutliers
1317 maxIterationsPtcOutliers = self.config.maxIterationsPtcOutliers
1319 for i, ampName in enumerate(dataset.ampNames):
1320 # The PTC is sorted by detector means,
1321 # however this might not be the true sorted
1322 # order for any specific amplifier. We must
1323 # therefore make sure that the amplifier
1324 # itself is sorted by the amplifier means.
1325 # At the end we will re-order the mask to
1326 # match the sort by detector means.
1328 # These are sorted by detector means
1329 meanVecRaw = dataset.rawMeans[ampName].copy()
1330 varVecRaw = dataset.rawVars[ampName].copy()
1331 varVecRaw = self._makeZeroSafe(varVecRaw)
1333 # These are sorted by the amplifier means
1334 meanVecSort = np.argsort(meanVecRaw)
1335 meanVecSorted = meanVecRaw[meanVecSort]
1336 varVecSorted = varVecRaw[meanVecSort]
1338 # Check if all points are bad from the 'cpExtractPtcTask'
1339 goodPoints = dataset.expIdMask[ampName].copy()
1340 goodPoints = goodPoints[meanVecSort]
1342 if not goodPoints.any():
1343 msg = ("SERIOUS: All points in goodPoints are bad. "
1344 f"Setting {ampName} to BAD.")
1345 self.log.warning(msg)
1346 # Fill entries with NaNs
1347 self.fillBadAmp(dataset, ampName)
1348 continue
1350 mask = goodPoints.copy()
1352 # Set the fitting function and search region
1353 ptcFunc = funcAstier
1354 parsIniPtc = [-2e-6, 1.5, 25.] # a00, gain, noise^2
1355 if self.config.binSize > 1:
1356 bounds = self._boundsForAstier(parsIniPtc)
1357 else:
1358 bounds = self._boundsForAstier(
1359 parsIniPtc,
1360 lowers=[-1e-4, 0.1, 0],
1361 uppers=[0, 10.0, 2000],
1362 )
1364 # We perform an initial (unweighted) fit of variance vs signal
1365 # (after initial KS test or post-drop selection) to look for
1366 # outliers, particularly at the high-flux end. The initial fit
1367 # is performed only for points that are guaranteed to be below
1368 # the PTC turnoff and then extrapolated to ensure that high
1369 # flux points that have abnormal variance values can be properly
1370 # rejected in this phase without biasing the initial fit.
1371 # This algorithm was initially developed by Seth Digel for
1372 # the EO Testing pipeline.
1373 if self.config.scaleMaxSignalInitialPtcOutlierFit:
1374 approxGain = np.nanmedian(meanVecRaw/varVecRaw)
1375 maxADUInitialPtcOutlierFit = self.config.maxSignalInitialPtcOutlierFit/approxGain
1376 maxDeltaADUInitialPtcOutlierFit = self.config.maxDeltaInitialPtcOutlierFit/approxGain
1377 self.log.info(
1378 "Using approximate gain %.3f and ADU signal cutoff of %.1f "
1379 "and max delta %.1f for amplifier %s",
1380 approxGain,
1381 maxADUInitialPtcOutlierFit,
1382 maxDeltaADUInitialPtcOutlierFit,
1383 ampName,
1384 )
1385 else:
1386 maxADUInitialPtcOutlierFit = self.config.maxSignalInitialPtcOutlierFit
1387 maxDeltaADUInitialPtcOutlierFit = self.config.maxDeltaInitialPtcOutlierFit
1389 if maxIterationsPtcOutliers == 0:
1390 # We are not doing any outlier rejection here, but we do want
1391 # an initial fit.
1392 res = least_squares(
1393 errFunc,
1394 parsIniPtc,
1395 bounds=bounds,
1396 args=(meanVecSorted[mask], varVecSorted[mask]),
1397 )
1398 pars = res.x
1399 newMask = mask.copy()
1400 else:
1401 newMask = (
1402 mask
1403 & (meanVecSorted <= maxADUInitialPtcOutlierFit)
1404 & (varVecSorted <= maxADUInitialPtcOutlierFit)
1405 )
1407 converged = False
1408 count = 0
1409 lastMask = mask.copy()
1410 while count < maxIterationsPtcOutliers:
1411 res = least_squares(
1412 errFunc,
1413 parsIniPtc,
1414 bounds=bounds,
1415 args=(meanVecSorted[newMask], varVecSorted[newMask]),
1416 )
1417 pars = res.x
1419 sigResids = (varVecSorted - ptcFunc(pars, meanVecSorted))/np.sqrt(varVecSorted)
1420 # The new mask includes points where the residuals are
1421 # finite, are less than the cut, and include the original
1422 # mask of known points that should not be used.
1423 newMask = (
1424 np.isfinite(sigResids)
1425 & (np.abs(np.nan_to_num(sigResids)) < sigmaCutPtcOutliers)
1426 & mask
1427 )
1428 # Demand at least 2 points to continue.
1429 if np.count_nonzero(newMask) < 2:
1430 msg = (f"SERIOUS: All points after outlier rejection are bad. "
1431 f"Setting {ampName} to BAD.")
1432 self.log.warning(msg)
1433 # Fill entries with NaNs
1434 self.fillBadAmp(dataset, ampName)
1435 break
1437 self.log.debug(
1438 "Iteration %d: Removed %d points in total for %s.",
1439 count,
1440 np.count_nonzero(mask) - np.count_nonzero(newMask),
1441 ampName,
1442 )
1444 # Loop over all used (True) points. If one of them follows
1445 # a False point, then it must be within
1446 # maxDeltaADUInitialPtcOutlierFit of a True point. If it
1447 # is a large gap, everything above is marked False.
1448 useMask, = np.where(newMask)
1449 for useIndex, usePoint in enumerate(useMask):
1450 if useIndex == 0 or newMask[usePoint - 1]:
1451 # The previous point was good; continue.
1452 continue
1453 deltaADU = meanVecSorted[usePoint] - meanVecSorted[useMask[useIndex - 1]]
1454 if deltaADU < maxDeltaADUInitialPtcOutlierFit:
1455 # This jump is fine; continue.
1456 continue
1458 self.log.info(
1459 "Found a jump of %.2f for amp %s, greater than %.2f. Masking points higher "
1460 "than %.2f ADU.",
1461 deltaADU,
1462 ampName,
1463 maxDeltaADUInitialPtcOutlierFit,
1464 meanVecSorted[useMask[useIndex - 1]],
1465 )
1467 # Mark all further points bad.
1468 newMask[usePoint:] = False
1469 break
1471 # If the mask hasn't changed then break out.
1472 if np.all(newMask == lastMask):
1473 self.log.debug("Convergence at iteration %d; breaking loop for %s.", count, ampName)
1474 converged = True
1475 break
1477 lastMask = newMask.copy()
1479 count += 1
1481 if not converged and maxIterationsPtcOutliers > 0:
1482 self.log.warning(
1483 "Outlier detection was not converged prior to %d iteration for %s",
1484 count,
1485 ampName
1486 )
1488 # Set the mask to the new mask, and reset the sorting.
1489 mask = np.zeros(len(meanVecSort), dtype=np.bool_)
1490 mask[meanVecSort[newMask]] = True
1491 maskSorted = newMask.copy()
1493 if not mask.any():
1494 # We have already filled the bad amp above, so continue.
1495 continue
1497 dataset.expIdMask[ampName] = mask
1499 parsIniPtc = pars
1500 meanVecFinal = meanVecRaw[mask]
1501 varVecFinal = varVecRaw[mask]
1503 if Counter(mask)[False] > 0:
1504 self.log.info("Number of points discarded in PTC of amplifier %s:"
1505 " %d out of %d", ampName, Counter(mask)[False], len(meanVecSorted))
1507 if (len(meanVecFinal) < len(parsIniPtc)):
1508 msg = (f"SERIOUS: Not enough data points ({len(meanVecFinal)}) compared to the number of "
1509 f"parameters of the PTC model({len(parsIniPtc)}). Setting {ampName} to BAD.")
1510 self.log.warning(msg)
1511 # Fill entries with NaNs
1512 self.fillBadAmp(dataset, ampName)
1513 continue
1515 # Fit the PTC.
1516 # The variance of the variance is Var(v)=2*v^2/Npix. This is
1517 # already calculated in `makeCovArray` of CpPtcExtract.
1518 # dataset.covariancesSqrtWeights[ampName][:,0,0]
1519 # has 1/sqrt(Var(v)).
1520 weightsY = dataset.covariancesSqrtWeights[ampName][:, 0, 0][mask]
1521 if self.config.doFitBootstrap:
1522 parsFit, parsFitErr, reducedChiSqPtc = fitBootstrap(parsIniPtc, meanVecFinal,
1523 varVecFinal, ptcFunc,
1524 weightsY=weightsY)
1525 else:
1526 parsFit, parsFitErr, reducedChiSqPtc = fitLeastSq(parsIniPtc, meanVecFinal,
1527 varVecFinal, ptcFunc,
1528 weightsY=weightsY)
1530 # Determine PTC turnoff
1531 if computePtcTurnoff:
1532 dataset.ptcTurnoff[ampName] = np.max(meanVecFinal)
1533 # And compute the ptcTurnoffSamplingError as one half the
1534 # difference between the previous and next point.
1535 lastGoodIndex = np.where(maskSorted)[0][-1]
1536 ptcTurnoffLow = meanVecSorted[lastGoodIndex - 1]
1537 if lastGoodIndex == (len(meanVecSorted) - 1):
1538 # If it's the last index, just use the interval.
1539 ptcTurnoffSamplingError = dataset.ptcTurnoff[ampName] - ptcTurnoffLow
1540 elif not np.isfinite(meanVecSorted[lastGoodIndex + 1]):
1541 # If the next index is not finite, just use the interval.
1542 ptcTurnoffSamplingError = dataset.ptcTurnoff[ampName] - ptcTurnoffLow
1543 else:
1544 ptcTurnoffSamplingError = (meanVecSorted[lastGoodIndex + 1] - ptcTurnoffLow)/2.
1545 dataset.ptcTurnoffSamplingError[ampName] = ptcTurnoffSamplingError
1547 dataset.ptcFitPars[ampName] = parsFit
1548 dataset.ptcFitParsError[ampName] = parsFitErr
1549 dataset.ptcFitChiSq[ampName] = reducedChiSqPtc
1551 dataset.finalVars[ampName] = varVecRaw
1552 dataset.finalVars[ampName][~mask] = np.nan
1553 dataset.finalModelVars[ampName] = ptcFunc(parsFit, meanVecRaw)
1554 dataset.finalModelVars[ampName][~mask] = np.nan
1555 dataset.finalMeans[ampName] = meanVecRaw
1556 dataset.finalMeans[ampName][~mask] = np.nan
1558 ptcGain = parsFit[1]
1559 ptcGainErr = parsFitErr[1]
1560 ptcNoise = np.sqrt(np.fabs(parsFit[2]))
1561 ptcNoiseErr = 0.5*(parsFitErr[2]/np.fabs(parsFit[2]))*np.sqrt(np.fabs(parsFit[2]))
1563 # Save results
1564 dataset.gain[ampName] = ptcGain
1565 dataset.gainUnadjusted[ampName] = ptcGain
1566 dataset.gainErr[ampName] = ptcGainErr
1567 dataset.noise[ampName] = ptcNoise
1568 dataset.noiseErr[ampName] = ptcNoiseErr
1570 if len(dataset.badAmps) == 0:
1571 dataset.badAmps = []
1573 return dataset
1575 def fillBadAmp(self, dataset, ampName):
1576 """Fill the dataset with NaNs if there are not enough
1577 good points.
1579 Parameters
1580 ----------
1581 dataset : `lsst.ip.isr.ptcDataset.PhotonTransferCurveDataset`
1582 The dataset containing the means, variances and
1583 exposure times.
1584 ampName : `str`
1585 Amplifier name.
1586 """
1587 if ampName not in dataset.badAmps:
1588 dataset.badAmps.append(ampName)
1589 dataset.expIdMask[ampName] = np.repeat(False, len(dataset.rawExpTimes[ampName]))
1590 dataset.gain[ampName] = np.nan
1591 dataset.gainUnadjusted[ampName] = np.nan
1592 dataset.gainErr[ampName] = np.nan
1593 dataset.noise[ampName] = np.nan
1594 dataset.noiseErr[ampName] = np.nan
1595 dataset.ptcFitPars[ampName] = (np.repeat(np.nan, 3))
1596 dataset.ptcFitParsError[ampName] = (np.repeat(np.nan, 3))
1597 dataset.ptcFitChiSq[ampName] = np.nan
1598 dataset.ptcTurnoff[ampName] = np.nan
1599 dataset.finalVars[ampName] = np.repeat(np.nan, len(dataset.rawExpTimes[ampName]))
1600 dataset.finalModelVars[ampName] = np.repeat(np.nan, len(dataset.rawExpTimes[ampName]))
1601 dataset.finalMeans[ampName] = np.repeat(np.nan, len(dataset.rawExpTimes[ampName]))
1603 return