Coverage for python / lsst / ip / isr / ptcDataset.py: 4%
629 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:00 +0000
1#
2# LSST Data Management System
3# Copyright 2008-2017 AURA/LSST.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
22"""
23Define dataset class for MeasurePhotonTransferCurve task
24"""
26__all__ = ['PhotonTransferCurveDataset']
28import warnings
29import numbers
30import numpy as np
31import math
32from astropy.table import Table
33import numpy.polynomial.polynomial as poly
34from scipy.signal import fftconvolve
36from lsst.ip.isr import IsrCalib
39def symmetrize(inputArray):
40 """ Copy array over 4 quadrants prior to convolution.
42 Parameters
43 ----------
44 inputarray : `numpy.array`
45 Input array to symmetrize.
47 Returns
48 -------
49 aSym : `numpy.array`
50 Symmetrized array.
51 """
53 targetShape = list(inputArray.shape)
54 r1, r2 = inputArray.shape[-1], inputArray.shape[-2]
55 targetShape[-1] = 2*r1-1
56 targetShape[-2] = 2*r2-1
57 aSym = np.ndarray(tuple(targetShape))
58 aSym[..., r2-1:, r1-1:] = inputArray
59 aSym[..., r2-1:, r1-1::-1] = inputArray
60 aSym[..., r2-1::-1, r1-1::-1] = inputArray
61 aSym[..., r2-1::-1, r1-1:] = inputArray
63 return aSym
66class PhotonTransferCurveDataset(IsrCalib):
67 """A simple class to hold the output data from the PTC task.
69 The dataset is made up of a dictionary for each item, keyed by the
70 amplifiers' names, which much be supplied at construction time.
71 New items cannot be added to the class to save accidentally saving to the
72 wrong property, and the class can be frozen if desired.
73 inputExpIdPairs records the exposures used to produce the data.
74 When fitPtc() or fitCovariancesAstier() is run, a mask is built up, which
75 is by definition always the same length as inputExpIdPairs, rawExpTimes,
76 rawMeans and rawVars, and is a list of bools, which are incrementally set
77 to False as points are discarded from the fits.
79 Parameters
80 ----------
81 ampNames : `list`
82 List with the names of the amplifiers of the detector at hand.
83 ptcFitType : `str`, optional
84 Type of model fitted to the PTC: "EXPAPPROXIMATION",
85 or "FULLCOVARIANCE" or "FULLCOVARIANCE_NO_B".
86 covMatrixSide : `int`, optional
87 Maximum lag of measured covariances (size of square covariance
88 matrices).
89 covMatrixSideFullCovFit : `int`, optional
90 Maximum covariances lag for FULLCOVARIANCE fit. It should be less or
91 equal than covMatrixSide.
92 kwargs : `dict`, optional
93 Other keyword arguments to pass to the parent init.
95 Notes
96 -----
97 The stored attributes are:
99 badAmps : `list` [`str`]
100 List with bad amplifiers names.
101 inputExpIdPairs : `dict`, [`str`, `list`]
102 Dictionary keyed by amp names containing the input exposures IDs.
103 inputExpPairMjdStartList : `dict`, [`str`, `np.ndarray`]
104 Dictionary keyed by amp names containing the start mjd from
105 the first exposure in each flat pair.
106 expIdMask : `dict`, [`str`, `np.ndarray`]
107 Dictionary keyed by amp names containing the mask produced after
108 outlier rejection. The mask produced by the "FULLCOVARIANCE"
109 option may differ from the one produced in the other two PTC
110 fit types.
111 expIdRolloffMask : `dict`, [`str`, `np.ndarray`]
112 Dictionary keyed by amp names containing the mask produced after
113 outlier rejection and extended to fit the PTC rolloff.
114 rawExpTimes : `dict`, [`str`, `np.ndarray`]
115 Dictionary keyed by amp names containing the unmasked exposure times.
116 rawMeans : `dict`, [`str`, `np.ndarray`]
117 Dictionary keyed by amp names containing the unmasked average of the
118 means of the exposures in each flat pair (units: adu).
119 rawVars : `dict`, [`str`, `np.ndarray`]
120 Dictionary keyed by amp names containing the variance of the
121 difference image of the exposures in each flat pair (units: adu^2).
122 rawDeltas : `dict`, [`str`, `np.ndarray`]
123 Dictionary keyed by amp names containing the scaled unmasked delta of
124 the means of the exposures in each flat pair (mean2 - mean1)
125 (units: adu).
126 rowMeanVariance : `dict`, [`str`, `np.ndarray`]
127 Dictionary keyed by amp names containing the variance of the
128 means of the rows of the difference image of the exposures
129 in each flat pair (units: adu^2).
130 histVars : `dict`, [`str`, `np.ndarray`]
131 Dictionary keyed by amp names containing the variance of the
132 difference image of the exposures in each flat pair estimated
133 by fitting a Gaussian model.
134 histChi2Dofs : `dict`, [`str`, `np.ndarray`]
135 Dictionary keyed by amp names containing the chi-squared per degree
136 of freedom fitting the difference image to a Gaussian model.
137 kspValues : `dict`, [`str`, `np.ndarray`]
138 Dictionary keyed by amp names containing the KS test p-value from
139 fitting the difference image to a Gaussian model.
140 gain : `dict`, [`str`, `float`]
141 Dictionary keyed by amp names containing the fitted gains. May be
142 adjusted by amp-offset gain ratios if configured in PTC solver.
143 gainUnadjusted : `dict`, [`str`, `float`]
144 Dictionary keyed by amp names containing unadjusted (raw) fit gain
145 values. May be the same as gain values if amp-offset adjustment
146 is not turned on.
147 gainErr : `dict`, [`str`, `float`]
148 Dictionary keyed by amp names containing the errors on the
149 fitted gains.
150 gainList : `dict`, [`str`, `np.ndarray`]
151 Dictionary keyed by amp names containing the gain estimated from
152 each flat pair.
153 overscanMedianLevelList : `dict`, [`str`, `np.ndarray`]
154 Dictionary keyed by amp names containing the median overscan
155 level from each input flat pair (units: adu).
156 overscanMedian : `dict `, [`str`, `float`]
157 Dictionary keyed by amp names containing the median of
158 overscanMedianLevelList[expIdMask] (units: adu).
159 overscanMedianSigma : `dict `, [`str`, `float`]
160 Dictionary keyed by amp names containing the median absolute
161 deviation of overscanMedianLevelList[expIdMask] (units: adu).
162 noiseList : `dict`, [`str`, `np.ndarray`]
163 Dictionary keyed by amp names containing the mean overscan
164 standard deviation from each flat pair (units: adu).
165 noise : `dict`, [`str`, `float`]
166 Dictionary keyed by amp names containing the fitted noise
167 (units: electron).
168 noiseErr : `dict`, [`str`, `float`]
169 Dictionary keyed by amp names containing the errors on the fitted
170 noise (units: electron).
171 ampOffsets : `dict`, [`str`, `float`]
172 Dictionary keyed by amp names containing amp-to-amp offsets
173 (units: adu).
174 ptcFitPars : `dict`, [`str`, `np.ndarray`]
175 Dictionary keyed by amp names containing the raveled array fitted
176 parameters of the PTC model.
177 ptcFitParsError : `dict`, [`str`, `np.ndarray`]
178 Dictionary keyed by amp names containing the raveled array errors
179 on the fitted parameters of the PTC model for ptcFitType.
180 ptcFitChiSq : `dict`, [`str`, `float`]
181 Dictionary keyed by amp names containing the reduced chi squared
182 of the fit.
183 ptcTurnoff : `dict` [`str, `float`]
184 Flux value (in adu) where the variance of the PTC curve starts
185 decreasing consistently.
186 ptcTurnoffSamplingError : `dict` [`str`, `float`]
187 ``Sampling`` error on the ptcTurnoff, based on the flux sampling
188 of the input PTC (units: adu).
189 ptcRolloff : `dict` [`str, `float`]
190 Flux value (in adu) where the variance of the PTC curve begins starts
191 to move away from the PTC model by some threshold.
192 ptcRolloffError : `dict` [`str`, `float`]
193 Covariance error from the fit to the PTC rolloff (units: adu).
194 ptcRolloff : `dict` [`str, `float`]
195 Curvature parameter of the PTC rolloff.
196 ptcRolloffError : `dict` [`str`, `float`]
197 Fitting error of the curvature parameter of the PTC rolloff.
198 nPixelCovariances : `dict`, [`str`, `int`]
199 Dictionary keyed by amp names containing the number of pixels
200 that were used to measure the covariances.
201 covariances : `dict`, [`str`, `np.ndarray`]
202 Dictionary keyed by amp names containing a list of measured
203 covariances per mean flux (units: adu^2).
204 covariancesModel : `dict`, [`str`, `np.ndarray`]
205 Dictionary keyed by amp names containinging covariances model
206 (Eq. 20 of Astier+19) per mean flux (units: adu^2).
207 covariancesSqrtWeights : `dict`, [`str`, `np.ndarray`]
208 Dictionary keyed by amp names containinging sqrt. of covariances
209 weights (units: 1/adu).
210 aMatrix : `dict`, [`str`, `np.ndarray`]
211 Dictionary keyed by amp names containing the "a" parameters from
212 the model in Eq. 20 of Astier+19 (units: 1/electron).
213 bMatrix : `dict`, [`str`, `np.ndarray`]
214 Dictionary keyed by amp names containing the "b" parameters from
215 the model in Eq. 20 of Astier+19 (units: 1/electron).
216 noiseMatrix : `dict`, [`str`, `np.ndarray`]
217 Dictionary keyed by amp names containing the "noise" parameters from
218 the model in Eq. 20 of Astier+19 (units: electron^2).
219 covariancesModelNoB : `dict`, [`str`, `np.ndarray`]
220 Dictionary keyed by amp names containing covariances model
221 (with 'b'=0 in Eq. 20 of Astier+19) per mean flux (units:
222 adu^2). Will be removed after v29.
223 aMatrixNoB : `dict`, [`str`, `np.ndarray`]
224 Dictionary keyed by amp names containing the "a" parameters from the
225 model in Eq. 20 of Astier+19 (and 'b' = 0) (units: 1/electron).
226 Will be removed after v29.
227 noiseMatrixNoB : `dict`, [`str`, `np.ndarray`]
228 Dictionary keyed by amp names containing the "noise" parameters from
229 the model in Eq. 20 of Astier+19, with 'b' = 0 (units: electron^2).
230 Will be removed after v29.
231 finalVars : `dict`, [`str`, `np.ndarray`]
232 Dictionary keyed by amp names containing the masked variance of the
233 difference image of each flat
234 pair. If needed, each array will be right-padded with
235 np.nan to match the length of rawExpTimes.
236 finalModelVars : `dict`, [`str`, `np.ndarray`]
237 Dictionary keyed by amp names containing the masked modeled
238 variance of the difference image of each flat pair. If needed, each
239 array will be right-padded with np.nan to match the length of
240 rawExpTimes.
241 finalMeans : `dict`, [`str`, `np.ndarray`]
242 Dictionary keyed by amp names containing the masked average of the
243 means of the exposures in each flat pair. If needed, each array
244 will be right-padded with np.nan to match the length of
245 rawExpTimes.
246 photoCharges : `dict`, [`str`, `np.ndarray`]
247 Dictionary keyed by amp names containing the integrated photocharge
248 for linearity calibration.
249 photoChargeDeltas : `dict`, [`str`, `np.ndarray`]
250 Dictionary keyed by amp names containing the delta for the integrated
251 photocharge (photocharge2 - photocharge1).
252 auxValues : `dict`, [`str`, `np.ndarray`]
253 Dictionary of per-detector auxiliary header values that can be used
254 for PTC, linearity computation.
256 Version 1.1 adds the `ptcTurnoff` attribute.
257 Version 1.2 adds the `histVars`, `histChi2Dofs`, and `kspValues`
258 attributes.
259 Version 1.3 adds the `noiseMatrix` and `noiseMatrixNoB` attributes.
260 Version 1.4 adds the `auxValues` attribute.
261 Version 1.5 adds the `covMatrixSideFullCovFit` attribute.
262 Version 1.6 adds the `rowMeanVariance` attribute.
263 Version 1.7 adds the `noiseList` attribute.
264 Version 1.8 adds the `ptcTurnoffSamplingError` attribute.
265 Version 1.9 standardizes PTC noise units to electron.
266 Version 2.0 adds the `ampOffsets`, `gainUnadjusted`, and
267 `gainList` attributes.
268 Version 2.1 deprecates the `covariancesModelNoB`, `aMatrixNoB`, and
269 `noiseMatrixNoB` attributes.
270 Version 2.2 adds the `overscanMedianLevelList` and
271 `inputExpPairMjdStartList` attributes.
272 Version 2.3 adds the `overscanMedian` and
273 `overscanMedianSigma` attrbutes.
274 Version 2.4 adds the `nPixelCovariances` attribute.
275 Version 2.5 adds the `rawDeltas` and `photoChargeDeltas` attributes.
276 Version 2.6 adds the `expIdRolloffMask`, `ptcRolloff`, `ptcRolloffError`,
277 `ptcRolloffTau`, and `ptcRolloffTauError` attributes. Also
278 deprecates the POLYNOMIAL fit type.
279 """
281 _OBSTYPE = 'PTC'
282 _SCHEMA = 'Gen3 Photon Transfer Curve'
283 # When adding a new field to update the version, be sure to update the
284 # following methods:
285 # * __init__()
286 # * fromDict()
287 # * toDict()
288 # * fromTable()
289 # * toTable()
290 # * setAmpValuesPartialDataset()
291 # * appendPartialPtc()
292 # * sort()
293 # * _checkTypes() in test_ptcDataset.py
294 # * test_ptcDataset() in test_ptcDataset.py
295 # * test_ptcDatasetSort in test_ptcDataset.py
296 # * test_ptcDatasetAppend in test_ptcDataset.py
297 _VERSION = 2.6
299 def __init__(self, ampNames=[], ptcFitType=None, covMatrixSide=1,
300 covMatrixSideFullCovFit=None, **kwargs):
301 self.ptcFitType = ptcFitType
302 self.ampNames = ampNames
303 self.covMatrixSide = covMatrixSide
304 if covMatrixSideFullCovFit is None:
305 self.covMatrixSideFullCovFit = covMatrixSide
306 else:
307 self.covMatrixSideFullCovFit = covMatrixSideFullCovFit
309 self.badAmps = []
311 self.inputExpIdPairs = {ampName: [] for ampName in ampNames}
312 self.inputExpPairMjdStartList = {ampName: np.array([]) for ampName in ampNames}
313 self.expIdMask = {ampName: np.array([], dtype=bool) for ampName in ampNames}
314 self.expIdRolloffMask = {ampName: np.array([], dtype=bool) for ampName in ampNames}
315 self.rawExpTimes = {ampName: np.array([]) for ampName in ampNames}
316 self.rawMeans = {ampName: np.array([]) for ampName in ampNames}
317 self.rawVars = {ampName: np.array([]) for ampName in ampNames}
318 self.rawDeltas = {ampName: np.array([]) for ampName in ampNames}
319 self.rowMeanVariance = {ampName: np.array([]) for ampName in ampNames}
320 self.photoCharges = {ampName: np.array([]) for ampName in ampNames}
321 self.photoChargeDeltas = {ampName: np.array([]) for ampName in ampNames}
322 self.ampOffsets = {ampName: np.array([]) for ampName in ampNames}
324 self.gain = {ampName: np.nan for ampName in ampNames}
325 self.gainUnadjusted = {ampName: np.nan for ampName in ampNames}
326 self.gainErr = {ampName: np.nan for ampName in ampNames}
327 self.gainList = {ampName: np.array([]) for ampName in ampNames}
328 self.overscanMedianLevelList = {ampName: np.array([]) for ampName in ampNames}
329 self.overscanMedian = {ampName: np.nan for ampName in ampNames}
330 self.overscanMedianSigma = {ampName: np.nan for ampName in ampNames}
331 self.noiseList = {ampName: np.array([]) for ampName in ampNames}
332 self.noise = {ampName: np.nan for ampName in ampNames}
333 self.noiseErr = {ampName: np.nan for ampName in ampNames}
335 self.histVars = {ampName: np.array([]) for ampName in ampNames}
336 self.histChi2Dofs = {ampName: np.array([]) for ampName in ampNames}
337 self.kspValues = {ampName: np.array([]) for ampName in ampNames}
339 self.ptcFitPars = {ampName: np.array([]) for ampName in ampNames}
340 self.ptcFitParsError = {ampName: np.array([]) for ampName in ampNames}
341 self.ptcFitChiSq = {ampName: np.nan for ampName in ampNames}
342 self.ptcTurnoff = {ampName: np.nan for ampName in ampNames}
343 self.ptcTurnoffSamplingError = {ampName: np.nan for ampName in ampNames}
344 self.ptcRolloff = {ampName: np.nan for ampName in ampNames}
345 self.ptcRolloffError = {ampName: np.nan for ampName in ampNames}
346 self.ptcRolloffTau = {ampName: np.nan for ampName in ampNames}
347 self.ptcRolloffTauError = {ampName: np.nan for ampName in ampNames}
349 self.nPixelCovariances = {ampName: -1 for ampName in ampNames}
350 self.covariances = {ampName: np.array([]) for ampName in ampNames}
351 self.covariancesModel = {ampName: np.array([]) for ampName in ampNames}
352 self.covariancesSqrtWeights = {ampName: np.array([]) for ampName in ampNames}
353 self.aMatrix = {ampName: np.array([]) for ampName in ampNames}
354 self.bMatrix = {ampName: np.array([]) for ampName in ampNames}
355 self.noiseMatrix = {ampName: np.array([]) for ampName in ampNames}
357 self.finalVars = {ampName: np.array([]) for ampName in ampNames}
358 self.finalModelVars = {ampName: np.array([]) for ampName in ampNames}
359 self.finalMeans = {ampName: np.array([]) for ampName in ampNames}
361 # Try this as a dict of arrays.
362 self.auxValues = {}
364 super().__init__(**kwargs)
365 self.requiredAttributes.update(['badAmps', 'inputExpIdPairs', 'inputExpPairMjdStartList',
366 'expIdMask', 'expIdRolloffMask', 'rawExpTimes', 'rawMeans',
367 'rawVars', 'rowMeanVariance', 'gain', 'gainErr', 'gainList',
368 'noise', 'noiseErr', 'noiseList', 'overscanMedianLevelList',
369 'overscanMedian', 'overscanMedianSigma', 'ptcFitPars',
370 'ptcFitParsError', 'ptcFitChiSq', 'ptcTurnoff', 'covariances',
371 'covariancesModel', 'covariancesSqrtWeights', 'aMatrix',
372 'bMatrix', 'noiseMatrix', 'finalVars', 'finalModelVars',
373 'finalMeans', 'photoCharges', 'histVars', 'histChi2Dofs',
374 'kspValues', 'auxValues', 'ptcTurnoffSamplingError',
375 'ptcRolloff', 'ptcRolloffError', 'ptcRolloffTau',
376 'ptcRolloffTauError', 'ampOffsets', 'gainUnadjusted',
377 'nPixelCovariances', 'rawDeltas', 'photoChargeDeltas'])
379 self.updateMetadata(setCalibInfo=True, setCalibId=True, **kwargs)
380 self._validateCovarianceMatrixSizes()
382 def setAmpValuesPartialDataset(
383 self,
384 ampName,
385 inputExpIdPair=(-1, -1),
386 inputExpPairMjdStart=np.nan,
387 rawExpTime=np.nan,
388 rawMean=np.nan,
389 rawVar=np.nan,
390 rawDelta=np.nan,
391 rowMeanVariance=np.nan,
392 photoCharge=np.nan,
393 photoChargeDelta=np.nan,
394 ampOffset=np.nan,
395 expIdMask=False,
396 expIdRolloffMask=False,
397 nPixelCovariance=-1,
398 covariance=None,
399 covSqrtWeights=None,
400 gain=np.nan,
401 noise=np.nan,
402 overscanMedianLevel=np.nan,
403 histVar=np.nan,
404 histChi2Dof=np.nan,
405 kspValue=0.0,
406 ):
407 """
408 Set the amp values for a partial PTC Dataset (from cpExtractPtcTask).
410 Parameters
411 ----------
412 ampName : `str`
413 Name of the amp to set the values.
414 inputExpIdPair : `tuple` [`int`]
415 Exposure IDs of input pair.
416 inputExpPairMjdStart : `float`, optional
417 The start MJD of first exposure in the flat pair.
418 rawExpTime : `float`, optional
419 Exposure time for this exposure pair (units: sec).
420 rawMean : `float`, optional
421 Average of the means of the exposures in this pair
422 (units: adu).
423 rawVar : `float`, optional
424 Variance of the difference of the exposures in this pair
425 (units: adu^2).
426 rawDelta : `float`, optional
427 Delta of the means of the exposure in this pair
428 (units: adu).
429 rowMeanVariance : `float`, optional
430 Variance of the means of the rows in the difference image
431 of the exposures in this pair (units: adu^2).
432 photoCharge : `float`, optional
433 Integrated photocharge for flat pair for linearity calibration
434 (arbitrary units).
435 photoChargeDelta : `float`, optional
436 Delta between integrated photocharge for the flat pair
437 (arbitrary units).
438 ampOffset : `float`, optional
439 Amp offset for this amplifier.
440 expIdMask : `bool`, optional
441 Flag setting if this exposure pair should be used (True)
442 or not used (False).
443 expIdRolloffMask : `bool`, optional
444 Flag setting if this exposure pair should be used (True)
445 or not used for rolloff search (False).
446 nPixelCovariance : `int`, optional
447 Number of pixels that went into the covariance measurement.
448 covariance : `np.ndarray` or None, optional
449 Measured covariance for this exposure pair (units: adu^2).
450 covSqrtWeights : `np.ndarray` or None, optional
451 Measured sqrt of covariance weights in this exposure pair
452 (units: 1/adu).
453 gain : `float`, optional
454 Estimated gain for this exposure pair (units: electron/adu).
455 noise : `float`, optional
456 Estimated read noise for this exposure pair (units: electron).
457 overscanMedianLevel : `float`, optional
458 Average of the median overscan levels for this exposure pair.
459 (units: adu)
460 histVar : `float`, optional
461 Variance estimated from fitting a histogram with a Gaussian model
462 (units: adu).
463 histChi2Dof : `float`, optional
464 Chi-squared per degree of freedom from Gaussian histogram fit.
465 kspValue : `float`, optional
466 KS test p-value from the Gaussian histogram fit.
467 """
468 nanMatrix = np.full((self.covMatrixSide, self.covMatrixSide), np.nan)
469 nanMatrixFit = np.full((self.covMatrixSideFullCovFit,
470 self.covMatrixSideFullCovFit), np.nan)
471 if covariance is None:
472 covariance = nanMatrix
473 if covSqrtWeights is None:
474 covSqrtWeights = nanMatrix
476 self.inputExpIdPairs[ampName] = [inputExpIdPair]
477 self.inputExpPairMjdStartList[ampName] = np.array([inputExpPairMjdStart])
478 self.rawExpTimes[ampName] = np.array([rawExpTime])
479 self.rawMeans[ampName] = np.array([rawMean])
480 self.rawVars[ampName] = np.array([rawVar])
481 self.rawDeltas[ampName] = np.array([rawDelta])
482 self.rowMeanVariance[ampName] = np.array([rowMeanVariance])
483 self.photoCharges[ampName] = np.array([photoCharge])
484 self.photoChargeDeltas[ampName] = np.array([photoChargeDelta])
485 self.ampOffsets[ampName] = np.array([ampOffset])
486 self.expIdMask[ampName] = np.array([expIdMask])
487 self.expIdRolloffMask[ampName] = np.array([expIdRolloffMask])
488 self.nPixelCovariances[ampName] = nPixelCovariance
489 self.covariances[ampName] = np.array([covariance])
490 self.covariancesSqrtWeights[ampName] = np.array([covSqrtWeights])
491 self.gain[ampName] = gain
492 self.gainUnadjusted[ampName] = gain
493 self.gainList[ampName] = np.array([gain])
494 self.noise[ampName] = noise
495 self.noiseList[ampName] = np.array([noise])
496 self.overscanMedianLevelList[ampName] = np.array([overscanMedianLevel])
497 self.overscanMedian[ampName] = float(overscanMedianLevel)
498 self.overscanMedianSigma[ampName] = float(0.0)
499 self.histVars[ampName] = np.array([histVar])
500 self.histChi2Dofs[ampName] = np.array([histChi2Dof])
501 self.kspValues[ampName] = np.array([kspValue])
503 # From FULLCOVARIANCE model
504 self.covariancesModel[ampName] = np.array([nanMatrixFit])
505 self.aMatrix[ampName] = nanMatrixFit
506 self.bMatrix[ampName] = nanMatrixFit
507 self.noiseMatrix[ampName] = nanMatrixFit
509 # Filler values.
510 self.finalVars[ampName] = np.array([np.nan])
511 self.finalModelVars[ampName] = np.array([np.nan])
512 self.finalMeans[ampName] = np.array([np.nan])
514 def setAuxValuesPartialDataset(self, auxDict):
515 """
516 Set a dictionary of auxiliary values for a partial dataset.
518 Parameters
519 ----------
520 auxDict : `dict` [`str`, `float`]
521 Dictionary of float values.
522 """
523 for key, value in auxDict.items():
524 if isinstance(value, numbers.Integral):
525 self.auxValues[key] = np.atleast_1d(np.asarray(value).astype(np.int64))
526 elif isinstance(value, (str, np.str_, np.bytes_)):
527 self.auxValues[key] = np.atleast_1d(np.asarray(value))
528 else:
529 self.auxValues[key] = np.atleast_1d(np.array(value, dtype=np.float64))
531 def updateMetadata(self, **kwargs):
532 """Update calibration metadata.
533 This calls the base class's method after ensuring the required
534 calibration keywords will be saved.
536 Parameters
537 ----------
538 setDate : `bool`, optional
539 Update the CALIBDATE fields in the metadata to the current
540 time. Defaults to False.
541 kwargs :
542 Other keyword parameters to set in the metadata.
543 """
544 super().updateMetadata(PTC_FIT_TYPE=self.ptcFitType, **kwargs)
546 @classmethod
547 def fromDict(cls, dictionary):
548 """Construct a calibration from a dictionary of properties.
549 Must be implemented by the specific calibration subclasses.
551 Parameters
552 ----------
553 dictionary : `dict`
554 Dictionary of properties.
556 Returns
557 -------
558 calib : `lsst.ip.isr.PhotonTransferCurveDataset`
559 Constructed calibration.
561 Raises
562 ------
563 RuntimeError
564 Raised if the supplied dictionary is for a different
565 calibration.
566 """
567 calib = cls()
568 if calib._OBSTYPE != dictionary['metadata']['OBSTYPE']:
569 raise RuntimeError(f"Incorrect Photon Transfer Curve dataset supplied. "
570 f"Expected {calib._OBSTYPE}, found {dictionary['metadata']['OBSTYPE']}")
571 calib.setMetadata(dictionary['metadata'])
572 calib.ptcFitType = dictionary['ptcFitType']
573 calib.covMatrixSide = dictionary['covMatrixSide']
574 calib.covMatrixSideFullCovFit = dictionary['covMatrixSideFullCovFit']
575 calib.badAmps = np.array(dictionary['badAmps'], 'str').tolist()
576 calib.ampNames = []
578 # The cov matrices are square
579 covMatrixSide = calib.covMatrixSide
580 covMatrixSideFullCovFit = calib.covMatrixSideFullCovFit
581 # Number of final signal levels
582 covDimensionsProduct = len(np.array(list(dictionary['covariances'].values())[0]).ravel())
583 nSignalPoints = int(covDimensionsProduct/(covMatrixSide*covMatrixSide))
585 for ampName in dictionary['ampNames']:
586 calib.ampNames.append(ampName)
587 calib.inputExpIdPairs[ampName] = dictionary['inputExpIdPairs'][ampName]
588 calib.inputExpPairMjdStartList[ampName] = np.array(
589 dictionary['inputExpPairMjdStartList'][ampName],
590 )
591 calib.expIdMask[ampName] = np.array(dictionary['expIdMask'][ampName])
592 calib.expIdRolloffMask[ampName] = np.array(dictionary['expIdRolloffMask'][ampName])
593 calib.rawExpTimes[ampName] = np.array(dictionary['rawExpTimes'][ampName], dtype=np.float64)
594 calib.rawMeans[ampName] = np.array(dictionary['rawMeans'][ampName], dtype=np.float64)
595 calib.rawVars[ampName] = np.array(dictionary['rawVars'][ampName], dtype=np.float64)
596 calib.rawDeltas[ampName] = np.array(dictionary['rawDeltas'][ampName], dtype=np.float64)
597 calib.rowMeanVariance[ampName] = np.array(dictionary['rowMeanVariance'][ampName],
598 dtype=np.float64)
599 calib.gain[ampName] = float(dictionary['gain'][ampName])
600 calib.gainErr[ampName] = float(dictionary['gainErr'][ampName])
601 calib.gainUnadjusted[ampName] = float(dictionary['gainUnadjusted'][ampName])
602 calib.gainList[ampName] = np.array(dictionary['gainList'][ampName], dtype=np.float64)
603 calib.noiseList[ampName] = np.array(dictionary['noiseList'][ampName], dtype=np.float64)
604 calib.nPixelCovariances[ampName] = int(dictionary['nPixelCovariances'][ampName])
605 calib.overscanMedianLevelList[ampName] = np.array(
606 dictionary['overscanMedianLevelList'][ampName],
607 dtype=np.float64,
608 )
609 calib.overscanMedian[ampName] = float(dictionary['overscanMedian'][ampName])
610 calib.overscanMedianSigma[ampName] = float(dictionary['overscanMedianSigma'][ampName])
611 calib.noise[ampName] = float(dictionary['noise'][ampName])
612 calib.noiseErr[ampName] = float(dictionary['noiseErr'][ampName])
613 calib.histVars[ampName] = np.array(dictionary['histVars'][ampName], dtype=np.float64)
614 calib.histChi2Dofs[ampName] = np.array(dictionary['histChi2Dofs'][ampName], dtype=np.float64)
615 calib.kspValues[ampName] = np.array(dictionary['kspValues'][ampName], dtype=np.float64)
616 calib.ptcFitPars[ampName] = np.array(dictionary['ptcFitPars'][ampName], dtype=np.float64)
617 calib.ptcFitParsError[ampName] = np.array(dictionary['ptcFitParsError'][ampName],
618 dtype=np.float64)
619 calib.ptcFitChiSq[ampName] = float(dictionary['ptcFitChiSq'][ampName])
620 calib.ptcTurnoff[ampName] = float(dictionary['ptcTurnoff'][ampName])
621 calib.ptcTurnoffSamplingError[ampName] = float(dictionary['ptcTurnoffSamplingError'][ampName])
622 calib.ptcRolloff[ampName] = float(dictionary['ptcRolloff'][ampName])
623 calib.ptcRolloffError[ampName] = float(dictionary['ptcRolloffError'][ampName])
624 calib.ptcRolloffTau[ampName] = float(dictionary['ptcRolloff'][ampName])
625 calib.ptcRolloffTauError[ampName] = float(dictionary['ptcRolloffError'][ampName])
626 if nSignalPoints > 0:
627 # Regular dataset
628 calib.covariances[ampName] = np.array(dictionary['covariances'][ampName],
629 dtype=np.float64).reshape(
630 (nSignalPoints, covMatrixSide, covMatrixSide))
631 calib.covariancesModel[ampName] = np.array(
632 dictionary['covariancesModel'][ampName],
633 dtype=np.float64).reshape(
634 (nSignalPoints, covMatrixSideFullCovFit, covMatrixSideFullCovFit))
635 calib.covariancesSqrtWeights[ampName] = np.array(
636 dictionary['covariancesSqrtWeights'][ampName],
637 dtype=np.float64).reshape(
638 (nSignalPoints, covMatrixSide, covMatrixSide))
639 calib.aMatrix[ampName] = np.array(dictionary['aMatrix'][ampName],
640 dtype=np.float64).reshape(
641 (covMatrixSideFullCovFit, covMatrixSideFullCovFit))
642 calib.bMatrix[ampName] = np.array(dictionary['bMatrix'][ampName],
643 dtype=np.float64).reshape(
644 (covMatrixSideFullCovFit, covMatrixSideFullCovFit))
645 calib.noiseMatrix[ampName] = np.array(
646 dictionary['noiseMatrix'][ampName],
647 dtype=np.float64).reshape((covMatrixSideFullCovFit, covMatrixSideFullCovFit))
648 else:
649 # Empty dataset
650 calib.covariances[ampName] = np.array([], dtype=np.float64)
651 calib.covariancesModel[ampName] = np.array([], dtype=np.float64)
652 calib.covariancesSqrtWeights[ampName] = np.array([], dtype=np.float64)
653 calib.aMatrix[ampName] = np.array([], dtype=np.float64)
654 calib.bMatrix[ampName] = np.array([], dtype=np.float64)
655 calib.noiseMatrix[ampName] = np.array([], dtype=np.float64)
657 calib.finalVars[ampName] = np.array(dictionary['finalVars'][ampName], dtype=np.float64)
658 calib.finalModelVars[ampName] = np.array(dictionary['finalModelVars'][ampName], dtype=np.float64)
659 calib.finalMeans[ampName] = np.array(dictionary['finalMeans'][ampName], dtype=np.float64)
660 calib.photoCharges[ampName] = np.array(dictionary['photoCharges'][ampName], dtype=np.float64)
661 calib.photoChargeDeltas[ampName] = np.array(
662 dictionary['photoChargeDeltas'][ampName],
663 dtype=np.float64,
664 )
665 calib.ampOffsets[ampName] = np.array(dictionary['ampOffsets'][ampName], dtype=np.float64)
667 for key, value in dictionary['auxValues'].items():
668 if isinstance(value[0], numbers.Integral):
669 calib.auxValues[key] = np.atleast_1d(np.asarray(value).astype(np.int64))
670 elif isinstance(value[0], (str, np.str_, np.bytes_)):
671 calib.auxValues[key] = np.atleast_1d(np.asarray(value))
672 else:
673 calib.auxValues[key] = np.atleast_1d(np.array(value, dtype=np.float64))
675 calib.updateMetadata()
676 return calib
678 def toDict(self):
679 """Return a dictionary containing the calibration properties.
680 The dictionary should be able to be round-tripped through
681 `fromDict`.
683 Returns
684 -------
685 dictionary : `dict`
686 Dictionary of properties.
687 """
688 self.updateMetadata()
690 outDict = dict()
691 metadata = self.getMetadata()
692 outDict['metadata'] = metadata
694 def _dictOfArraysToDictOfLists(dictOfArrays):
695 dictOfLists = {}
696 for key, value in dictOfArrays.items():
697 dictOfLists[key] = value.ravel().tolist()
699 return dictOfLists
701 outDict['ptcFitType'] = self.ptcFitType
702 outDict['covMatrixSide'] = self.covMatrixSide
703 outDict['covMatrixSideFullCovFit'] = self.covMatrixSideFullCovFit
704 outDict['ampNames'] = self.ampNames
705 outDict['badAmps'] = self.badAmps
706 outDict['inputExpIdPairs'] = self.inputExpIdPairs
707 outDict['inputExpPairMjdStartList'] = _dictOfArraysToDictOfLists(self.inputExpPairMjdStartList)
708 outDict['expIdMask'] = _dictOfArraysToDictOfLists(self.expIdMask)
709 outDict['expIdRolloffMask'] = _dictOfArraysToDictOfLists(self.expIdRolloffMask)
710 outDict['rawExpTimes'] = _dictOfArraysToDictOfLists(self.rawExpTimes)
711 outDict['rawMeans'] = _dictOfArraysToDictOfLists(self.rawMeans)
712 outDict['rawVars'] = _dictOfArraysToDictOfLists(self.rawVars)
713 outDict['rawDeltas'] = _dictOfArraysToDictOfLists(self.rawDeltas)
714 outDict['rowMeanVariance'] = _dictOfArraysToDictOfLists(self.rowMeanVariance)
715 outDict['gain'] = self.gain
716 outDict['gainErr'] = self.gainErr
717 outDict['gainUnadjusted'] = self.gainUnadjusted
718 outDict['gainList'] = _dictOfArraysToDictOfLists(self.gainList)
719 outDict['noiseList'] = _dictOfArraysToDictOfLists(self.noiseList)
720 outDict['overscanMedianLevelList'] = _dictOfArraysToDictOfLists(self.overscanMedianLevelList)
721 outDict['overscanMedian'] = self.overscanMedian
722 outDict['overscanMedianSigma'] = self.overscanMedianSigma
723 outDict['noise'] = self.noise
724 outDict['noiseErr'] = self.noiseErr
725 outDict['histVars'] = self.histVars
726 outDict['histChi2Dofs'] = self.histChi2Dofs
727 outDict['kspValues'] = self.kspValues
728 outDict['ptcFitPars'] = _dictOfArraysToDictOfLists(self.ptcFitPars)
729 outDict['ptcFitParsError'] = _dictOfArraysToDictOfLists(self.ptcFitParsError)
730 outDict['ptcFitChiSq'] = self.ptcFitChiSq
731 outDict['ptcTurnoff'] = self.ptcTurnoff
732 outDict['ptcTurnoffSamplingError'] = self.ptcTurnoffSamplingError
733 outDict['ptcRolloff'] = self.ptcRolloff
734 outDict['ptcRolloffError'] = self.ptcRolloffError
735 outDict['ptcRolloffTau'] = self.ptcRolloff
736 outDict['ptcRolloffTauError'] = self.ptcRolloffTauError
737 outDict['nPixelCovariances'] = self.nPixelCovariances
738 outDict['covariances'] = _dictOfArraysToDictOfLists(self.covariances)
739 outDict['covariancesModel'] = _dictOfArraysToDictOfLists(self.covariancesModel)
740 outDict['covariancesSqrtWeights'] = _dictOfArraysToDictOfLists(self.covariancesSqrtWeights)
741 outDict['aMatrix'] = _dictOfArraysToDictOfLists(self.aMatrix)
742 outDict['bMatrix'] = _dictOfArraysToDictOfLists(self.bMatrix)
743 outDict['noiseMatrix'] = _dictOfArraysToDictOfLists(self.noiseMatrix)
744 outDict['finalVars'] = _dictOfArraysToDictOfLists(self.finalVars)
745 outDict['finalModelVars'] = _dictOfArraysToDictOfLists(self.finalModelVars)
746 outDict['finalMeans'] = _dictOfArraysToDictOfLists(self.finalMeans)
747 outDict['photoCharges'] = _dictOfArraysToDictOfLists(self.photoCharges)
748 outDict['photoChargeDeltas'] = _dictOfArraysToDictOfLists(self.photoChargeDeltas)
749 outDict['ampOffsets'] = _dictOfArraysToDictOfLists(self.ampOffsets)
750 outDict['auxValues'] = _dictOfArraysToDictOfLists(self.auxValues)
752 return outDict
754 @classmethod
755 def fromTable(cls, tableList):
756 """Construct calibration from a list of tables.
757 This method uses the `fromDict` method to create the
758 calibration, after constructing an appropriate dictionary from
759 the input tables.
761 Parameters
762 ----------
763 tableList : `list` [`lsst.afw.table.Table`]
764 List of tables to use to construct the datasetPtc.
766 Returns
767 -------
768 calib : `lsst.ip.isr.PhotonTransferCurveDataset`
769 The calibration defined in the tables.
770 """
771 ptcTable = tableList[0]
773 metadata = ptcTable.meta
774 inDict = dict()
775 inDict['metadata'] = metadata
776 inDict['ampNames'] = []
777 inDict['ptcFitType'] = []
778 inDict['covMatrixSide'] = []
779 inDict['covMatrixSideFullCovFit'] = []
780 inDict['inputExpIdPairs'] = dict()
781 inDict['inputExpPairMjdStartList'] = dict()
782 inDict['expIdMask'] = dict()
783 inDict['expIdRolloffMask'] = dict()
784 inDict['rawExpTimes'] = dict()
785 inDict['rawMeans'] = dict()
786 inDict['rawVars'] = dict()
787 inDict['rawDeltas'] = dict()
788 inDict['rowMeanVariance'] = dict()
789 inDict['gain'] = dict()
790 inDict['gainErr'] = dict()
791 inDict['gainUnadjusted'] = dict()
792 inDict['gainList'] = dict()
793 inDict['noiseList'] = dict()
794 inDict['overscanMedianLevelList'] = dict()
795 inDict['overscanMedian'] = dict()
796 inDict['overscanMedianSigma'] = dict()
797 inDict['noise'] = dict()
798 inDict['noiseErr'] = dict()
799 inDict['histVars'] = dict()
800 inDict['histChi2Dofs'] = dict()
801 inDict['kspValues'] = dict()
802 inDict['ptcFitPars'] = dict()
803 inDict['ptcFitParsError'] = dict()
804 inDict['ptcFitChiSq'] = dict()
805 inDict['ptcTurnoff'] = dict()
806 inDict['ptcTurnoffSamplingError'] = dict()
807 inDict['ptcRolloff'] = dict()
808 inDict['ptcRolloffError'] = dict()
809 inDict['ptcRolloffTau'] = dict()
810 inDict['ptcRolloffTauError'] = dict()
811 inDict['nPixelCovariances'] = dict()
812 inDict['covariances'] = dict()
813 inDict['covariancesModel'] = dict()
814 inDict['covariancesSqrtWeights'] = dict()
815 inDict['aMatrix'] = dict()
816 inDict['bMatrix'] = dict()
817 inDict['noiseMatrix'] = dict()
818 inDict['finalVars'] = dict()
819 inDict['finalModelVars'] = dict()
820 inDict['finalMeans'] = dict()
821 inDict['badAmps'] = []
822 inDict['photoCharges'] = dict()
823 inDict['photoChargeDeltas'] = dict()
824 inDict['ampOffsets'] = dict()
826 # TODO: DM-47610, remove after v29
827 inDict['noiseMatrixNoB'] = dict()
828 inDict['covariancesModelNoB'] = dict()
829 inDict['aMatrixNoB'] = dict()
831 calibVersion = metadata['PTC_VERSION']
832 if calibVersion == 1.0:
833 cls().log.warning(f"Previous version found for PTC dataset: {calibVersion}. "
834 f"Setting 'ptcTurnoff' in all amps to last value in 'finalMeans'.")
835 for record in ptcTable:
836 ampName = record['AMPLIFIER_NAME']
838 inDict['ptcFitType'] = record['PTC_FIT_TYPE']
839 inDict['covMatrixSide'] = record['COV_MATRIX_SIDE']
840 inDict['ampNames'].append(ampName)
841 inDict['inputExpIdPairs'][ampName] = record['INPUT_EXP_ID_PAIRS'].tolist()
842 inDict['expIdMask'][ampName] = record['EXP_ID_MASK']
843 inDict['rawExpTimes'][ampName] = record['RAW_EXP_TIMES']
844 inDict['rawMeans'][ampName] = record['RAW_MEANS']
845 inDict['rawVars'][ampName] = record['RAW_VARS']
846 inDict['gain'][ampName] = record['GAIN']
847 inDict['gainErr'][ampName] = record['GAIN_ERR']
848 inDict['noise'][ampName] = record['NOISE']
849 inDict['noiseErr'][ampName] = record['NOISE_ERR']
850 inDict['ptcFitPars'][ampName] = record['PTC_FIT_PARS']
851 inDict['ptcFitParsError'][ampName] = record['PTC_FIT_PARS_ERROR']
852 inDict['ptcFitChiSq'][ampName] = record['PTC_FIT_CHI_SQ']
853 inDict['covariances'][ampName] = record['COVARIANCES']
854 inDict['covariancesModel'][ampName] = record['COVARIANCES_MODEL']
855 inDict['covariancesSqrtWeights'][ampName] = record['COVARIANCES_SQRT_WEIGHTS']
856 inDict['aMatrix'][ampName] = record['A_MATRIX']
857 inDict['bMatrix'][ampName] = record['B_MATRIX']
858 inDict['finalVars'][ampName] = record['FINAL_VARS']
859 inDict['finalModelVars'][ampName] = record['FINAL_MODEL_VARS']
860 inDict['finalMeans'][ampName] = record['FINAL_MEANS']
861 inDict['badAmps'] = record['BAD_AMPS'].tolist()
862 inDict['photoCharges'][ampName] = record['PHOTO_CHARGE']
863 if calibVersion == 1.0:
864 mask = record['FINAL_MEANS'].mask
865 array = record['FINAL_MEANS'][~mask]
866 if len(array) > 0:
867 inDict['ptcTurnoff'][ampName] = record['FINAL_MEANS'][~mask][-1]
868 else:
869 inDict['ptcTurnoff'][ampName] = np.nan
870 else:
871 inDict['ptcTurnoff'][ampName] = record['PTC_TURNOFF']
872 if calibVersion < 1.2:
873 inDict['histVars'][ampName] = np.array([np.nan])
874 inDict['histChi2Dofs'][ampName] = np.array([np.nan])
875 inDict['kspValues'][ampName] = np.array([0.0])
876 else:
877 inDict['histVars'][ampName] = record['HIST_VARS']
878 inDict['histChi2Dofs'][ampName] = record['HIST_CHI2_DOFS']
879 inDict['kspValues'][ampName] = record['KS_PVALUES']
880 if calibVersion < 1.3:
881 nanMatrix = np.full_like(inDict['aMatrix'][ampName], np.nan)
882 inDict['noiseMatrix'][ampName] = nanMatrix
883 inDict['noiseMatrixNoB'][ampName] = nanMatrix
884 elif calibVersion >= 1.3 and calibVersion < 2.1:
885 inDict['noiseMatrix'][ampName] = record['NOISE_MATRIX']
886 inDict['noiseMatrixNoB'][ampName] = record['NOISE_MATRIX_NO_B']
887 else:
888 inDict['noiseMatrix'][ampName] = record['NOISE_MATRIX']
889 nanMatrix = np.full_like(inDict['aMatrix'][ampName], np.nan)
890 inDict['noiseMatrixNoB'][ampName] = nanMatrix
891 if calibVersion < 1.5:
892 # Matched to `COV_MATRIX_SIDE`. Same for all amps.
893 inDict['covMatrixSideFullCovFit'] = inDict['covMatrixSide']
894 else:
895 inDict['covMatrixSideFullCovFit'] = record['COV_MATRIX_SIDE_FULL_COV_FIT']
896 if calibVersion < 1.6:
897 inDict['rowMeanVariance'][ampName] = np.full((len(inDict['expIdMask'][ampName]),), np.nan)
898 else:
899 inDict['rowMeanVariance'][ampName] = record['ROW_MEAN_VARIANCE']
900 if calibVersion < 1.7:
901 inDict['noiseList'][ampName] = np.full_like(inDict['rawMeans'][ampName], np.nan)
902 else:
903 inDict['noiseList'][ampName] = record['NOISE_LIST']
904 if calibVersion < 1.8:
905 inDict['ptcTurnoffSamplingError'][ampName] = np.nan
906 else:
907 inDict['ptcTurnoffSamplingError'][ampName] = record['PTC_TURNOFF_SAMPLING_ERROR']
908 if calibVersion < 1.9 and inDict['ptcFitType'] == "FULLCOVARIANCE":
909 # Before version 1.9, the noise stored in the PTC was in
910 # units of electron^2 only if ptcFitType == FULLCOVARIANCE.
911 # After version 1.9, we standardized the
912 # PhotonTransferCurveDataset noise units to electron to fix
913 # this bug. If a user tries to use an earlier version of
914 # PTC with this fit type, we must be sure to do the
915 # calculations properly. More information about this noise
916 # issue can be found in DM-45976.
917 if ampName == inDict['ampNames'][0]:
918 cls().log.info(f"Input PTC VERSION ({calibVersion}) < 1.9 and"
919 " ptcFitType == FULLCOVARIANCE. Applying fix for"
920 f" the DM-45976 noise issue.")
921 # The noiseErr calculation was accidentally correct in the
922 # previous version, so we only need to upday the noise
923 # attribute.
924 inDict['noise'][ampName] = np.sqrt(record['noise'][ampName])
925 if calibVersion < 2.0:
926 inDict['ampOffsets'][ampName] = np.full_like(inDict['rawMeans'][ampName], np.nan)
927 inDict['gainUnadjusted'][ampName] = record['GAIN']
928 inDict['gainList'][ampName] = np.full_like(inDict['rawMeans'][ampName], np.nan)
929 else:
930 inDict['ampOffsets'][ampName] = record['AMP_OFFSETS']
931 inDict['gainUnadjusted'][ampName] = record['GAIN_UNADJUSTED']
932 inDict['gainList'][ampName] = record['GAIN_LIST']
933 if calibVersion < 2.1:
934 inDict['covariancesModelNoB'][ampName] = record['COVARIANCES_MODEL_NO_B']
935 inDict['aMatrixNoB'][ampName] = record['A_MATRIX_NO_B']
936 else:
937 nanMatrixList = np.full_like(inDict['covariances'][ampName], np.nan)
938 inDict['covariancesModelNoB'][ampName] = nanMatrixList
939 nanMatrix = np.full_like(inDict['aMatrix'][ampName], np.nan)
940 inDict['aMatrixNoB'][ampName] = nanMatrix
941 if calibVersion < 2.2:
942 inDict['inputExpPairMjdStartList'][ampName] = np.full_like(
943 inDict['rawMeans'][ampName],
944 np.nan,
945 )
946 inDict['overscanMedianLevelList'][ampName] = np.full_like(
947 inDict['rawMeans'][ampName],
948 np.nan,
949 )
950 else:
951 inDict['inputExpPairMjdStartList'][ampName] = record['INPUT_EXP_PAIR_MJD_START']
952 inDict['overscanMedianLevelList'][ampName] = record['OVERSCAN_MEDIAN_LIST']
953 if calibVersion < 2.3:
954 inDict['overscanMedian'][ampName] = np.nan
955 inDict['overscanMedianSigma'][ampName] = np.nan
956 else:
957 inDict['overscanMedian'][ampName] = record['OVERSCAN_MEDIAN']
958 inDict['overscanMedianSigma'][ampName] = record['OVERSCAN_MEDIAN_SIGMA']
959 if calibVersion < 2.4:
960 inDict['nPixelCovariances'][ampName] = -1
961 else:
962 inDict['nPixelCovariances'][ampName] = record['NPIXEL_COVARIANCES']
963 if calibVersion < 2.5:
964 inDict['rawDeltas'][ampName] = np.full_like(
965 inDict['rawMeans'][ampName],
966 np.nan,
967 )
968 inDict['photoChargeDeltas'][ampName] = np.full_like(
969 inDict['rawMeans'][ampName],
970 np.nan,
971 )
972 else:
973 inDict['rawDeltas'][ampName] = record['RAW_DELTAS']
974 inDict['photoChargeDeltas'][ampName] = record['PHOTO_CHARGE_DELTAS']
975 if calibVersion < 2.6:
976 inDict['ptcRolloff'][ampName] = np.nan
977 inDict['ptcRolloffError'][ampName] = np.nan
978 inDict['ptcRolloffTau'][ampName] = np.nan
979 inDict['ptcRolloffTauError'][ampName] = np.nan
980 inDict['expIdRolloffMask'][ampName] = np.full_like(
981 inDict['expIdMask'][ampName],
982 False,
983 )
984 else:
985 inDict['ptcRolloff'][ampName] = record['PTC_ROLLOFF']
986 inDict['ptcRolloffError'][ampName] = record['PTC_ROLLOFF_ERROR']
987 inDict['ptcRolloffTau'][ampName] = record['PTC_ROLLOFF_TAU']
988 inDict['ptcRolloffTauError'][ampName] = record['PTC_ROLLOFF_TAU_ERROR']
989 inDict['expIdRolloffMask'][ampName] = record['EXP_ID_ROLLOFF_MASK']
991 inDict['auxValues'] = {}
992 record = ptcTable[0]
993 for col in record.columns.keys():
994 if col.startswith('PTCAUX_'):
995 parts = col.split('PTCAUX_')
996 if isinstance(record[col][0], np.bytes_):
997 # Convert to a unicode string because astropy fits doesn't
998 # round-trip properly
999 inDict['auxValues'][parts[1]] = record[col].astype(np.str_)
1000 else:
1001 inDict['auxValues'][parts[1]] = record[col]
1003 return cls().fromDict(inDict)
1005 def toTable(self):
1006 """Construct a list of tables containing the information in this
1007 calibration.
1009 The list of tables should create an identical calibration
1010 after being passed to this class's fromTable method.
1012 Returns
1013 -------
1014 tableList : `list` [`astropy.table.Table`]
1015 List of tables containing the linearity calibration
1016 information.
1017 """
1018 tableList = []
1019 self.updateMetadata()
1021 badAmps = np.array(self.badAmps) if len(self.badAmps) else np.array([], dtype="U3")
1023 catalogList = []
1024 for ampName in self.ampNames:
1025 ampDict = {
1026 'AMPLIFIER_NAME': ampName,
1027 'PTC_FIT_TYPE': self.ptcFitType,
1028 'COV_MATRIX_SIDE': self.covMatrixSide,
1029 'COV_MATRIX_SIDE_FULL_COV_FIT': self.covMatrixSideFullCovFit,
1030 'INPUT_EXP_ID_PAIRS': self.inputExpIdPairs[ampName],
1031 'INPUT_EXP_PAIR_MJD_START': self.inputExpPairMjdStartList[ampName],
1032 'EXP_ID_MASK': self.expIdMask[ampName],
1033 'EXP_ID_ROLLOFF_MASK': self.expIdRolloffMask[ampName],
1034 'RAW_EXP_TIMES': self.rawExpTimes[ampName],
1035 'RAW_MEANS': self.rawMeans[ampName],
1036 'RAW_VARS': self.rawVars[ampName],
1037 'RAW_DELTAS': self.rawDeltas[ampName],
1038 'ROW_MEAN_VARIANCE': self.rowMeanVariance[ampName],
1039 'GAIN': self.gain[ampName],
1040 'GAIN_ERR': self.gainErr[ampName],
1041 'GAIN_UNADJUSTED': self.gainUnadjusted[ampName],
1042 'GAIN_LIST': self.gainList[ampName],
1043 'OVERSCAN_MEDIAN_LIST': self.overscanMedianLevelList[ampName],
1044 'OVERSCAN_MEDIAN': self.overscanMedian[ampName],
1045 'OVERSCAN_MEDIAN_SIGMA': self.overscanMedianSigma[ampName],
1046 'NOISE_LIST': self.noiseList[ampName],
1047 'NOISE': self.noise[ampName],
1048 'NOISE_ERR': self.noiseErr[ampName],
1049 'HIST_VARS': self.histVars[ampName],
1050 'HIST_CHI2_DOFS': self.histChi2Dofs[ampName],
1051 'KS_PVALUES': self.kspValues[ampName],
1052 'PTC_FIT_PARS': np.array(self.ptcFitPars[ampName]),
1053 'PTC_FIT_PARS_ERROR': np.array(self.ptcFitParsError[ampName]),
1054 'PTC_FIT_CHI_SQ': self.ptcFitChiSq[ampName],
1055 'PTC_TURNOFF': self.ptcTurnoff[ampName],
1056 'PTC_TURNOFF_SAMPLING_ERROR': self.ptcTurnoffSamplingError[ampName],
1057 'PTC_ROLLOFF': self.ptcRolloff[ampName],
1058 'PTC_ROLLOFF_ERROR': self.ptcRolloffError[ampName],
1059 'PTC_ROLLOFF_TAU': self.ptcRolloffTau[ampName],
1060 'PTC_ROLLOFF_TAU_ERROR': self.ptcRolloffTauError[ampName],
1061 'A_MATRIX': self.aMatrix[ampName].ravel(),
1062 'B_MATRIX': self.bMatrix[ampName].ravel(),
1063 'NOISE_MATRIX': self.noiseMatrix[ampName].ravel(),
1064 'BAD_AMPS': badAmps,
1065 'PHOTO_CHARGE': self.photoCharges[ampName],
1066 'PHOTO_CHARGE_DELTAS': self.photoChargeDeltas[ampName],
1067 'AMP_OFFSETS': self.ampOffsets[ampName],
1068 'NPIXEL_COVARIANCES': self.nPixelCovariances[ampName],
1069 'COVARIANCES': self.covariances[ampName].ravel(),
1070 'COVARIANCES_MODEL': self.covariancesModel[ampName].ravel(),
1071 'COVARIANCES_SQRT_WEIGHTS': self.covariancesSqrtWeights[ampName].ravel(),
1072 'FINAL_VARS': self.finalVars[ampName],
1073 'FINAL_MODEL_VARS': self.finalModelVars[ampName],
1074 'FINAL_MEANS': self.finalMeans[ampName],
1075 }
1077 if self.auxValues:
1078 for key, value in self.auxValues.items():
1079 ampDict[f"PTCAUX_{key}"] = value
1081 catalogList.append(ampDict)
1083 catalog = Table(catalogList)
1085 inMeta = self.getMetadata().toDict()
1086 outMeta = {k: v for k, v in inMeta.items() if v is not None}
1087 outMeta.update({k: "" for k, v in inMeta.items() if v is None})
1088 catalog.meta = outMeta
1089 tableList.append(catalog)
1091 return tableList
1093 def fromDetector(self, detector):
1094 """Read metadata parameters from a detector.
1096 Parameters
1097 ----------
1098 detector : `lsst.afw.cameraGeom.detector`
1099 Input detector with parameters to use.
1101 Returns
1102 -------
1103 calib : `lsst.ip.isr.PhotonTransferCurveDataset`
1104 The calibration constructed from the detector.
1105 """
1107 pass
1109 def appendPartialPtc(self, partialPtc):
1110 """Append a partial PTC dataset to this dataset.
1112 Parameters
1113 ----------
1114 partialPtc : `lsst.ip.isr.PhotonTransferCurveDataset`
1115 Partial PTC to append. Should only have one element.
1116 """
1117 if self.ampNames != partialPtc.ampNames:
1118 raise ValueError("partialPtc has mis-matched amps.")
1119 if len(partialPtc.rawMeans[self.ampNames[0]]) != 1 or partialPtc.ptcFitType != "PARTIAL":
1120 raise ValueError("partialPtc does not appear to be the correct format.")
1122 # Record the initial length of the PTC, for checking auxValues.
1123 initialLength = len(self.expIdMask[self.ampNames[0]])
1125 for key, value in partialPtc.auxValues.items():
1126 if key in self.auxValues:
1127 self.auxValues[key] = np.append(self.auxValues[key], value)
1128 elif initialLength == 0:
1129 # This is the first partial, so we can set the dict key.
1130 self.auxValues[key] = value
1131 else:
1132 raise ValueError(f"partialPtc has mismatched auxValue key {key}.")
1134 for ampName in self.ampNames:
1135 if initialLength == 0:
1136 # This is the first partial, so we can set the dict key.
1137 self.nPixelCovariances[ampName] = partialPtc.nPixelCovariances[ampName]
1138 elif partialPtc.nPixelCovariances[ampName] != self.nPixelCovariances[ampName]:
1139 raise ValueError(f"partialPtc has mismatched nPixelCovariances for amp {ampName}.")
1141 for ampName in self.ampNames:
1142 # The partial dataset consists of lists of values for each
1143 # quantity. In the case of the input exposure pairs and the
1144 # input exposure MJDs, this is a list of tuples. In all cases
1145 # we only want the first (and only) element of the list.
1146 self.inputExpIdPairs[ampName].append(partialPtc.inputExpIdPairs[ampName][0])
1147 self.inputExpPairMjdStartList[ampName] = np.append(
1148 self.inputExpPairMjdStartList[ampName],
1149 partialPtc.inputExpPairMjdStartList[ampName][0],
1150 )
1151 self.expIdMask[ampName] = np.append(self.expIdMask[ampName],
1152 partialPtc.expIdMask[ampName][0])
1153 self.expIdRolloffMask[ampName] = np.append(self.expIdRolloffMask[ampName],
1154 partialPtc.expIdRolloffMask[ampName][0])
1155 self.rawExpTimes[ampName] = np.append(self.rawExpTimes[ampName],
1156 partialPtc.rawExpTimes[ampName][0])
1157 self.rawMeans[ampName] = np.append(self.rawMeans[ampName],
1158 partialPtc.rawMeans[ampName][0])
1159 self.rawVars[ampName] = np.append(self.rawVars[ampName],
1160 partialPtc.rawVars[ampName][0])
1161 self.rawDeltas[ampName] = np.append(self.rawDeltas[ampName],
1162 partialPtc.rawDeltas[ampName][0])
1163 self.rowMeanVariance[ampName] = np.append(self.rowMeanVariance[ampName],
1164 partialPtc.rowMeanVariance[ampName][0])
1165 self.photoCharges[ampName] = np.append(self.photoCharges[ampName],
1166 partialPtc.photoCharges[ampName][0])
1167 self.photoChargeDeltas[ampName] = np.append(self.photoChargeDeltas[ampName],
1168 partialPtc.photoChargeDeltas[ampName][0])
1169 self.ampOffsets[ampName] = np.append(self.ampOffsets[ampName],
1170 partialPtc.ampOffsets[ampName][0])
1171 self.histVars[ampName] = np.append(self.histVars[ampName],
1172 partialPtc.histVars[ampName][0])
1173 self.histChi2Dofs[ampName] = np.append(self.histChi2Dofs[ampName],
1174 partialPtc.histChi2Dofs[ampName][0])
1175 self.kspValues[ampName] = np.append(self.kspValues[ampName],
1176 partialPtc.kspValues[ampName][0])
1177 self.gainList[ampName] = np.append(self.gainList[ampName],
1178 partialPtc.gain[ampName])
1179 self.overscanMedianLevelList[ampName] = np.append(
1180 self.overscanMedianLevelList[ampName],
1181 partialPtc.overscanMedianLevelList[ampName][0],
1182 )
1183 self.noiseList[ampName] = np.append(self.noiseList[ampName],
1184 partialPtc.noise[ampName])
1185 self.finalVars[ampName] = np.append(self.finalVars[ampName],
1186 partialPtc.finalVars[ampName][0])
1187 self.finalModelVars[ampName] = np.append(self.finalModelVars[ampName],
1188 partialPtc.finalModelVars[ampName][0])
1189 self.finalMeans[ampName] = np.append(self.finalMeans[ampName],
1190 partialPtc.finalMeans[ampName][0])
1191 self.covariances[ampName] = np.append(
1192 self.covariances[ampName].ravel(),
1193 partialPtc.covariances[ampName].ravel()
1194 ).reshape(
1195 (
1196 len(self.rawExpTimes[ampName]),
1197 self.covMatrixSide,
1198 self.covMatrixSide,
1199 )
1200 )
1201 self.covariancesSqrtWeights[ampName] = np.append(
1202 self.covariancesSqrtWeights[ampName].ravel(),
1203 partialPtc.covariancesSqrtWeights[ampName].ravel()
1204 ).reshape(
1205 (
1206 len(self.rawExpTimes[ampName]),
1207 self.covMatrixSide,
1208 self.covMatrixSide,
1209 )
1210 )
1211 self.covariancesModel[ampName] = np.append(
1212 self.covariancesModel[ampName].ravel(),
1213 partialPtc.covariancesModel[ampName].ravel()
1214 ).reshape(
1215 (
1216 len(self.rawExpTimes[ampName]),
1217 self.covMatrixSide,
1218 self.covMatrixSide,
1219 )
1220 )
1222 def sort(self, sortIndex):
1223 """Sort the components of the PTC by a given sort index.
1225 The PTC is sorted in-place.
1227 Parameters
1228 ----------
1229 sortIndex : `list` or `np.ndarray`
1230 The sorting index, which must be the same length as
1231 the number of elements of the PTC.
1232 """
1233 index = np.atleast_1d(sortIndex)
1235 # First confirm everything matches.
1236 for ampName in self.ampNames:
1237 if len(index) != len(self.rawExpTimes[ampName]):
1238 raise ValueError(
1239 f"Length of sortIndex ({len(index)}) does not match number of PTC "
1240 f"elements ({len(self.rawExpTimes[ampName])})",
1241 )
1243 # Note that gain, gainUnadjusted, gainErr, noise, noiseErr,
1244 # ptcTurnoff, ptcTurnoffSamplingError, and the full covariance fit
1245 # parameters are global and not sorted by input pair.
1247 for ampName in self.ampNames:
1248 self.inputExpIdPairs[ampName] = np.array(
1249 self.inputExpIdPairs[ampName]
1250 )[index].tolist()
1251 self.inputExpPairMjdStartList[ampName] = self.inputExpPairMjdStartList[ampName][index]
1253 self.expIdMask[ampName] = self.expIdMask[ampName][index]
1254 self.expIdRolloffMask[ampName] = self.expIdRolloffMask[ampName][index]
1255 self.rawExpTimes[ampName] = self.rawExpTimes[ampName][index]
1256 self.rawMeans[ampName] = self.rawMeans[ampName][index]
1257 self.rawVars[ampName] = self.rawVars[ampName][index]
1258 self.rawDeltas[ampName] = self.rawDeltas[ampName][index]
1259 self.rowMeanVariance[ampName] = self.rowMeanVariance[ampName][index]
1260 self.photoCharges[ampName] = self.photoCharges[ampName][index]
1261 self.photoChargeDeltas[ampName] = self.photoChargeDeltas[ampName][index]
1262 self.ampOffsets[ampName] = self.ampOffsets[ampName][index]
1264 self.gainList[ampName] = self.gainList[ampName][index]
1265 self.noiseList[ampName] = self.noiseList[ampName][index]
1267 self.overscanMedianLevelList[ampName] = self.overscanMedianLevelList[ampName][index]
1269 self.histVars[ampName] = self.histVars[ampName][index]
1270 self.histChi2Dofs[ampName] = self.histChi2Dofs[ampName][index]
1271 self.kspValues[ampName] = self.kspValues[ampName][index]
1273 self.covariances[ampName] = self.covariances[ampName][index]
1274 self.covariancesSqrtWeights[ampName] = self.covariancesSqrtWeights[ampName][index]
1275 self.covariancesModel[ampName] = self.covariancesModel[ampName][index]
1277 self.finalVars[ampName] = self.finalVars[ampName][index]
1278 self.finalModelVars[ampName] = self.finalModelVars[ampName][index]
1279 self.finalMeans[ampName] = self.finalMeans[ampName][index]
1281 # Sort the auxiliary values which are not stored per-amp.
1282 for key, value in self.auxValues.items():
1283 self.auxValues[key] = value[index]
1285 def getExpIdsUsed(self, ampName):
1286 """Get the exposures used, i.e. not discarded, for a given amp.
1287 If no mask has been created yet, all exposures are returned.
1289 Parameters
1290 ----------
1291 ampName : `str`
1293 Returns
1294 -------
1295 expIdsUsed : `list` [`tuple`]
1296 List of pairs of exposure ids used in PTC.
1297 """
1298 if len(self.expIdMask[ampName]) == 0:
1299 return self.inputExpIdPairs[ampName]
1301 # if the mask exists it had better be the same length as the expIdPairs
1302 assert len(self.expIdMask[ampName]) == len(self.inputExpIdPairs[ampName])
1304 pairs = self.inputExpIdPairs[ampName]
1305 mask = self.expIdMask[ampName]
1306 # cast to bool required because numpy
1307 try:
1308 expIdsUsed = [(exp1, exp2) for ((exp1, exp2), m) in zip(pairs, mask) if m]
1309 except ValueError:
1310 self.log.warning("The PTC file was written incorrectly; you should rerun the "
1311 "PTC solve task if possible.")
1312 expIdsUsed = []
1313 for pairList, m in zip(pairs, mask):
1314 if m:
1315 expIdsUsed.append(pairList[0])
1317 return expIdsUsed
1319 def getGoodAmps(self):
1320 """Get the good amps from this PTC."""
1321 return [amp for amp in self.ampNames if amp not in self.badAmps]
1323 def getGoodPoints(self, ampName):
1324 """Get the good points used for a given amp in the PTC.
1326 Parameters
1327 ----------
1328 ampName : `str`
1329 Amplifier's name.
1331 Returns
1332 -------
1333 goodPoints : `np.ndarray`
1334 Boolean array of good points used in PTC.
1335 """
1336 return self.expIdMask[ampName]
1338 def validateGainNoiseTurnoffValues(self, ampName, doWarn=False):
1339 """Ensure the gain, read noise, and PTC turnoff have
1340 sensible values.
1342 Parameters
1343 ----------
1344 ampName : `str`
1345 Amplifier's name.
1346 """
1348 gain = self.gain[ampName]
1349 noise = self.noise[ampName]
1350 ptcTurnoff = self.ptcTurnoff[ampName]
1352 # Check if gain is not positive or is np.nan
1353 if not (isinstance(gain, (int, float)) and gain > 0) or math.isnan(gain):
1354 if doWarn:
1355 self.log.warning(f"Invalid gain value {gain}"
1356 " Setting to default: Gain=1")
1357 gain = 1
1359 # Check if noise is not positive or is np.nan
1360 if not (isinstance(noise, (int, float)) and noise > 0) or math.isnan(noise):
1361 if doWarn:
1362 self.log.warning(f"Invalid noise value: {noise}"
1363 " Setting to default: Noise=1")
1364 noise = 1
1366 # Check if ptcTurnoff is not positive or is np.nan
1367 if not (isinstance(ptcTurnoff, (int, float)) and ptcTurnoff > 0) or math.isnan(ptcTurnoff):
1368 if doWarn:
1369 self.log.warning(f"Invalid PTC turnoff value: {ptcTurnoff}"
1370 " Setting to default: PTC Turnoff=2e19")
1371 ptcTurnoff = 2e19
1373 self.gain[ampName] = gain
1374 self.noise[ampName] = noise
1375 self.ptcTurnoff[ampName] = ptcTurnoff
1377 def evalPtcModel(self, mu):
1378 """Computes the covariance model at specific signal levels.
1380 Parameters
1381 ----------
1382 mu : `numpy.array`, (N,)
1383 List of mean signals in ADU.
1385 Raises
1386 ------
1387 RuntimeError
1388 Raised if ptcFitType is invalid.
1390 Returns
1391 -------
1392 covModel : `numpy.array`, (N, M, M)
1393 Covariances model at mu (in ADU^2).
1395 Notes
1396 -----
1397 Computes the covModel for all mu, and it returns
1398 cov[N, M, M], where the variance model is cov[:,0,0].
1399 Both mu and cov are in ADUs and ADUs squared. This
1400 routine evaluates the approximation in Eq. 16 of
1401 Astier+19 (1905.08677)
1402 if self.ptcFitType == EXPAPPROXIMATION, and Eq. 20 of
1403 Astier+19 if self.ptcFitType == FULLCOVARIANCE(_NO_B).
1405 The EXPAPPROXIMATION model (Eq. 16 of Astier+19) is
1406 only an approximation for the variance (cov[0,0]),
1407 so the function returns covModel of shape (N,),
1408 representing an array of [C_{00}]
1409 if self.ptcFitType == EXPAPPROXIMATION.
1410 """
1412 ampNames = self.ampNames
1413 covModel = {ampName: np.array([]) for ampName in ampNames}
1415 if self.ptcFitType == "POLYNOMIAL":
1416 warnings.warn("The `POLYNOMIAL` fit type is deprecated; it will "
1417 "be removed from the Rubin Observatory Science "
1418 "Pipelines after v30.",
1419 category=FutureWarning)
1420 # TODO: DM-52720 - remove deprecated POLYNOMIAL fit
1421 # and legacy turnoff
1422 pars = self.ptcFitPars
1424 for ampName in ampNames:
1425 c00 = poly.polyval(mu, [*pars[ampName]])
1426 covModel[ampName] = c00
1427 elif self.ptcFitType == "EXPAPPROXIMATION":
1428 pars = self.ptcFitPars
1430 for ampName in ampNames:
1431 a00, gain, noise = pars[ampName]
1432 f1 = 0.5/(a00*gain*gain)*(np.exp(2*a00*mu*gain)-1)
1433 f2 = noise/(gain*gain)
1434 c00 = f1 + f2
1435 covModel[ampName] = c00
1437 elif self.ptcFitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]:
1438 for ampName in ampNames:
1439 noiseMatrix = self.noiseMatrix[ampName]
1440 gain = self.gain[ampName]
1441 aMatrix = self.aMatrix[ampName]
1442 bMatrix = self.bMatrix[ampName]
1443 cMatrix = aMatrix*bMatrix
1445 matrixSideFit = self.covMatrixSideFullCovFit
1446 sa = (matrixSideFit, matrixSideFit)
1448 # pad a with zeros and symmetrize
1449 aEnlarged = np.zeros((int(sa[0]*1.5)+1, int(sa[1]*1.5)+1))
1450 aEnlarged[0:sa[0], 0:sa[1]] = aMatrix
1451 aSym = symmetrize(aEnlarged)
1453 # pad c with zeros and symmetrize
1454 cEnlarged = np.zeros((int(sa[0]*1.5)+1, int(sa[1]*1.5)+1))
1455 cEnlarged[0:sa[0], 0:sa[1]] = cMatrix
1457 cSym = symmetrize(cEnlarged)
1458 a2 = fftconvolve(aSym, aSym, mode='same')
1459 a3 = fftconvolve(a2, aSym, mode='same')
1460 ac = fftconvolve(aSym, cSym, mode='same')
1461 (xc, yc) = np.unravel_index(np.abs(aSym).argmax(), a2.shape)
1463 a1 = aMatrix[np.newaxis, :, :]
1464 a2 = a2[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
1465 a3 = a3[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
1466 ac = ac[np.newaxis, xc:xc + matrixSideFit, yc:yc + matrixSideFit]
1467 c1 = cMatrix[np.newaxis, ::]
1469 # assumes that mu is 1d
1470 bigMu = mu[:, np.newaxis, np.newaxis]*gain
1471 # c(=a*b in Astier+19) also has a contribution to the last
1472 # term, that is absent for now.
1474 covModel[ampName] = (bigMu/(gain*gain)*(a1*bigMu+2./3.*(bigMu*bigMu)*(a2 + c1)
1475 + (1./3.*a3 + 5./6.*ac)*(bigMu*bigMu*bigMu))
1476 + noiseMatrix[np.newaxis, :, :]/gain**2)
1478 # add the Poisson term, and the read out noise (variance)
1479 covModel[ampName][:, 0, 0] += mu/gain
1480 else:
1481 raise RuntimeError("Cannot compute PTC model for "
1482 "ptcFitType %s." % self.ptcFitType)
1484 return covModel
1486 def _validateCovarianceMatrixSizes(self):
1487 """Ensure covMatrixSideFullCovFit <= covMatrixSide."""
1488 if self.covMatrixSideFullCovFit > self.covMatrixSide:
1489 self.log.warning("covMatrixSideFullCovFit > covMatrixSide "
1490 f"({self.covMatrixSideFullCovFit} > {self.covMatrixSide})."
1491 "Setting the former to the latter.")
1492 self.covMatrixSideFullCovFit = self.covMatrixSide