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