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