347 inputExpIdPair=(-1, -1),
348 inputExpPairMjdStart=np.nan,
352 rowMeanVariance=np.nan,
360 overscanMedianLevel=np.nan,
366 Set the amp values for a partial PTC Dataset (from cpExtractPtcTask).
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
381 rawVar : `float`, optional
382 Variance of the difference of the exposures in this pair
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
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)
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
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.
407 histVar : `float`, optional
408 Variance estimated from fitting a histogram with a Gaussian model
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.
418 if covariance
is None:
419 covariance = nanMatrix
420 if covSqrtWeights
is None:
421 covSqrtWeights = nanMatrix
426 self.
rawMeans[ampName] = np.array([rawMean])
427 self.
rawVars[ampName] = np.array([rawVar])
430 self.
ampOffsets[ampName] = np.array([ampOffset])
431 self.
expIdMask[ampName] = np.array([expIdMask])
434 self.
gain[ampName] = gain
436 self.
gainList[ampName] = np.array([gain])
437 self.
noise[ampName] = noise
438 self.
noiseList[ampName] = np.array([noise])
440 self.
histVars[ampName] = np.array([histVar])
442 self.
kspValues[ampName] = np.array([kspValue])
446 self.
aMatrix[ampName] = nanMatrixFit
447 self.
bMatrix[ampName] = nanMatrixFit
450 self.
finalVars[ampName] = np.array([np.nan])
493 """Construct a calibration from a dictionary of properties.
494 Must be implemented by the specific calibration subclasses.
499 Dictionary of properties.
503 calib : `lsst.ip.isr.PhotonTransferCurveDataset`
504 Constructed calibration.
509 Raised if the supplied dictionary is for a different
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()
523 calibVersion = dictionary[
'metadata'][
'PTC_VERSION']
526 covMatrixSide = calib.covMatrixSide
527 covMatrixSideFullCovFit = calib.covMatrixSideFullCovFit
529 covDimensionsProduct = len(np.array(list(dictionary[
'covariances'].values())[0]).ravel())
530 nSignalPoints = int(covDimensionsProduct/(covMatrixSide*covMatrixSide))
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],
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],
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],
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],
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:
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))
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))
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)
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)
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))
621 calib.auxValues[key] = np.atleast_1d(np.array(value, dtype=np.float64))
623 calib.updateMetadata()
627 """Return a dictionary containing the calibration properties.
628 The dictionary should be able to be round-tripped through
634 Dictionary of properties.
640 outDict[
'metadata'] = metadata
642 def _dictOfArraysToDictOfLists(dictOfArrays):
644 for key, value
in dictOfArrays.items():
645 dictOfLists[key] = value.ravel().tolist()
653 outDict[
'badAmps'] = self.
badAmps
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
664 outDict[
'gainList'] = _dictOfArraysToDictOfLists(self.
gainList)
665 outDict[
'noiseList'] = _dictOfArraysToDictOfLists(self.
noiseList)
667 outDict[
'noise'] = self.
noise
672 outDict[
'ptcFitPars'] = _dictOfArraysToDictOfLists(self.
ptcFitPars)
673 outDict[
'ptcFitParsError'] = _dictOfArraysToDictOfLists(self.
ptcFitParsError)
677 outDict[
'covariances'] = _dictOfArraysToDictOfLists(self.
covariances)
678 outDict[
'covariancesModel'] = _dictOfArraysToDictOfLists(self.
covariancesModel)
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)
991 """Append a partial PTC dataset to this dataset.
995 partialPtc : `lsst.ip.isr.PhotonTransferCurveDataset`
996 Partial PTC to append. Should only have one element.
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.")
1004 initialLength = len(self.expIdMask[self.ampNames[0]])
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:
1011 self.auxValues[key] = value
1013 raise ValueError(f
"partialPtc has mismatched auxValue key {key}.")
1015 for ampName
in self.ampNames:
1020 self.inputExpIdPairs[ampName].append(partialPtc.inputExpIdPairs[ampName][0])
1021 self.inputExpPairMjdStartList[ampName] = np.append(
1022 self.inputExpPairMjdStartList[ampName],
1023 partialPtc.inputExpPairMjdStartList[ampName][0],
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],
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])
1060 self.covariances[ampName] = np.append(
1061 self.covariances[ampName].ravel(),
1062 partialPtc.covariances[ampName].ravel()
1065 len(self.rawExpTimes[ampName]),
1070 self.covariancesSqrtWeights[ampName] = np.append(
1071 self.covariancesSqrtWeights[ampName].ravel(),
1072 partialPtc.covariancesSqrtWeights[ampName].ravel()
1075 len(self.rawExpTimes[ampName]),
1080 self.covariancesModel[ampName] = np.append(
1081 self.covariancesModel[ampName].ravel(),
1082 partialPtc.covariancesModel[ampName].ravel()
1085 len(self.rawExpTimes[ampName]),
1244 """Computes the covariance model at specific signal levels.
1248 mu : `numpy.array`, (N,)
1249 List of mean signals in ADU.
1254 Raised if ptcFitType is invalid.
1258 covModel : `numpy.array`, (N, M, M)
1259 Covariances model at mu (in ADU^2).
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.
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.
1281 covModel = {ampName: np.array([])
for ampName
in ampNames}
1286 for ampName
in ampNames:
1287 c00 = poly.polyval(mu, [*pars[ampName]])
1288 covModel[ampName] = c00
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)
1298 covModel[ampName] = c00
1300 elif self.
ptcFitType in [
"FULLCOVARIANCE",
"FULLCOVARIANCE_NO_B"]:
1301 for ampName
in ampNames:
1303 gain = self.
gain[ampName]
1304 aMatrix = self.
aMatrix[ampName]
1305 bMatrix = self.
bMatrix[ampName]
1306 cMatrix = aMatrix*bMatrix
1309 sa = (matrixSideFit, matrixSideFit)
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
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
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)
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, ::]
1333 bigMu = mu[:, np.newaxis, np.newaxis]*gain
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)
1342 covModel[ampName][:, 0, 0] += mu/gain
1344 raise RuntimeError(
"Cannot compute PTC model for "