diff --git a/pipelines/LSSTComCam/cpCti.yaml b/pipelines/LSSTComCam/cpCti.yaml new file mode 100644 index 00000000..b27bcc16 --- /dev/null +++ b/pipelines/LSSTComCam/cpCti.yaml @@ -0,0 +1,4 @@ +description: LSSTCam cti calibration construction +instrument: lsst.obs.lsst.LsstComCam +imports: + - location: $CP_PIPE_DIR/pipelines/_ingredients/cpCtiLSST.yaml diff --git a/python/lsst/cp/pipe/ptc/cpPtcSolve.py b/python/lsst/cp/pipe/ptc/cpPtcSolve.py index 3ed70a90..a253037b 100644 --- a/python/lsst/cp/pipe/ptc/cpPtcSolve.py +++ b/python/lsst/cp/pipe/ptc/cpPtcSolve.py @@ -82,7 +82,8 @@ class PhotonTransferCurveSolveConfig(pipeBase.PipelineTaskConfig, allowed={ "POLYNOMIAL": "n-degree polynomial (use 'polynomialFitDegree' to set 'n').", "EXPAPPROXIMATION": "Approximation in Astier+19 (Eq. 16).", - "FULLCOVARIANCE": "Full covariances model in Astier+19 (Eq. 20)" + "FULLCOVARIANCE_NO_B": "Full covariances model in Astier+19 (Eq. 15)", + "FULLCOVARIANCE": "Full covariances model in Astier+19 (Eq. 20)", } ) minMeanSignal = pexConfig.DictField( @@ -399,7 +400,7 @@ def run(self, inputCovariances, camera=None, detId=0): index = np.argsort(detectorMeans) datasetPtc.sort(index) - if self.config.ptcFitType == "FULLCOVARIANCE": + if self.config.ptcFitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]: # Fit the measured covariances vs mean signal to # the Astier+19 full model (Eq. 20). Before that # do a preliminary fit to the variance (C_00) vs mean @@ -419,7 +420,7 @@ def run(self, inputCovariances, camera=None, detId=0): # previous fit. for ampName in datasetPtc.ampNames: datasetPtc.expIdMask[ampName] = tempDatasetPtc.expIdMask[ampName] - datasetPtc.fitType = "FULLCOVARIANCE" + datasetPtc.fitType = self.config.ptcFitType datasetPtc = self.fitMeasurementsToModel(datasetPtc) # The other options are: self.config.ptcFitType in # ("EXPAPPROXIMATION", "POLYNOMIAL") @@ -485,7 +486,7 @@ def fitMeasurementsToModel(self, dataset): `PhotonTransferCurveDatase`. """ fitType = dataset.ptcFitType - if fitType in ["FULLCOVARIANCE", ]: + if fitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]: # This model uses the full covariance matrix in the fit. # The PTC is technically defined as variance vs signal, # with variance = Cov_00 @@ -560,10 +561,7 @@ def fitDataFullCovariance(self, dataset): dataset.covariancesSqrtWeights[ampName] = listNanMatrix dataset.aMatrix[ampName] = nanMatrixFit dataset.bMatrix[ampName] = nanMatrixFit - dataset.covariancesModelNoB[ampName] = listNanMatrixFit - dataset.aMatrixNoB[ampName] = nanMatrixFit dataset.noiseMatrix[ampName] = nanMatrixFit - dataset.noiseMatrixNoB[ampName] = nanMatrixFit dataset.expIdMask[ampName] = np.repeat(False, lenInputTimes) dataset.gain[ampName] = np.nan @@ -608,31 +606,45 @@ def fitDataFullCovariance(self, dataset): covSqrtWeightsAtAmpForFitMasked ) - # Fit full model (Eq. 20 of Astier+19) and same model with - # b=0 (c=0 in this code). + # Fit full model (Eq. 20 of Astier+19) pInit = np.concatenate((a0.ravel(), c0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None) - functionsDict = {'fullModel': self.funcFullCovarianceModel, - 'fullModelNoB': self.funcFullCovarianceModelNoB} - fitResults = {'fullModel': {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []}, - 'fullModelNoB': {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []}} - for key in functionsDict: - params, paramsErr, _ = fitLeastSq(pInit, muAtAmpMasked, - covAtAmpForFitMasked.ravel(), functionsDict[key], - weightsY=covSqrtWeightsAtAmpForFitMasked.ravel()) - a = params[:lenParams].reshape((matrixSideFit, matrixSideFit)) - c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) - noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit)) - gain = params[-1] - - fitResults[key]['a'] = a - fitResults[key]['c'] = c - fitResults[key]['noiseMatrix'] = noiseMatrix - fitResults[key]['gain'] = gain - fitResults[key]['paramsErr'] = paramsErr + + # Initialize empty results dictionary + fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []} + + # Pick the correct full covariance model function + ptcModel = self.funcFullCovarianceModel + if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": + ptcModel = self.funcFullCovarianceModelNoB + pInit = np.concatenate((a0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None) + + params, paramsErr, _ = fitLeastSq( + pInit, + muAtAmpMasked, + covAtAmpForFitMasked.ravel(), + ptcModel, + weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(), + ) + + if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": + zeros = np.zeros_like(params[:lenParams]) + params = np.insert(params, lenParams, zeros) + paramsErr = np.insert(paramsErr, lenParams, zeros) + + a = params[:lenParams].reshape((matrixSideFit, matrixSideFit)) + c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) + noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit)) + gain = params[-1] + + fitResults['a'] = a + fitResults['c'] = c + fitResults['noiseMatrix'] = noiseMatrix + fitResults['gain'] = gain + fitResults['paramsErr'] = paramsErr # Put the information in the PTC dataset - # Not used when ptcFitType is 'FULLCOVARIANCE' + # Not used when ptcFitType is 'FULLCOVARIANCE*' dataset.ptcFitPars[ampName] = np.array([np.nan]) dataset.ptcFitParsError[ampName] = np.array([np.nan]) dataset.ptcFitChiSq[ampName] = np.nan @@ -646,33 +658,24 @@ def fitDataFullCovariance(self, dataset): # masked amps. dataset.covariancesModel[ampName] = self.evalCovModel( muAtAmp, - fitResults['fullModel']['a'], - fitResults['fullModel']['c'], - fitResults['fullModel']['noiseMatrix'], - fitResults['fullModel']['gain'] + fitResults['a'], + fitResults['c'], + fitResults['noiseMatrix'], + fitResults['gain'], + setBtoZero=(dataset.ptcFitType == "FULLCOVARIANCE_NO_B"), ) dataset.covariancesSqrtWeights[ampName] = covSqrtWeightsAtAmp - dataset.aMatrix[ampName] = fitResults['fullModel']['a'] - dataset.bMatrix[ampName] = fitResults['fullModel']['c']/fitResults['fullModel']['a'] - dataset.covariancesModelNoB[ampName] = self.evalCovModel( - muAtAmp, - fitResults['fullModelNoB']['a'], - fitResults['fullModelNoB']['c'], - fitResults['fullModelNoB']['noiseMatrix'], - fitResults['fullModelNoB']['gain'], - setBtoZero=True - ) - dataset.aMatrixNoB[ampName] = fitResults['fullModelNoB']['a'] - dataset.gain[ampName] = fitResults['fullModel']['gain'] - dataset.gainErr[ampName] = fitResults['fullModel']['paramsErr'][-1] - readoutNoiseSquared = fitResults['fullModel']['noiseMatrix'][0][0] + dataset.aMatrix[ampName] = fitResults['a'] + dataset.bMatrix[ampName] = fitResults['c']/fitResults['a'] + dataset.gain[ampName] = fitResults['gain'] + dataset.gainErr[ampName] = fitResults['paramsErr'][-1] + readoutNoiseSquared = fitResults['noiseMatrix'][0][0] readoutNoise = np.sqrt(np.fabs(readoutNoiseSquared)) dataset.noise[ampName] = readoutNoise - readoutNoiseSquaredSigma = fitResults['fullModel']['paramsErr'][2*lenParams] + readoutNoiseSquaredSigma = fitResults['paramsErr'][2*lenParams] noiseErr = 0.5*(readoutNoiseSquaredSigma/np.fabs(readoutNoiseSquared))*readoutNoise dataset.noiseErr[ampName] = noiseErr - dataset.noiseMatrix[ampName] = fitResults['fullModel']['noiseMatrix'] - dataset.noiseMatrixNoB[ampName] = fitResults['fullModelNoB']['noiseMatrix'] + dataset.noiseMatrix[ampName] = fitResults['noiseMatrix'] dataset.finalVars[ampName] = covAtAmp[:, 0, 0].copy() dataset.finalVars[ampName][~maskAtAmp] = np.nan @@ -779,7 +782,7 @@ def funcFullCovarianceModelNoB(self, params, x): Parameters ---------- params : `list` - Parameters of the model: aMatrix, CMatrix, noiseMatrix, + Parameters of the model: aMatrix, noiseMatrix, gain (e/adu). x : `numpy.array`, (N,) Signal mu (adu) @@ -792,8 +795,8 @@ def funcFullCovarianceModelNoB(self, params, x): matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit lenParams = matrixSideFit*matrixSideFit aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit)) - cMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) - noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit)) + cMatrix = np.zeros_like(aMatrix) + noiseMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) gain = params[-1] return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=True).flatten() @@ -1090,10 +1093,7 @@ def fitPtc(self, dataset): dataset.covariancesModel[amp] = listNanMatrixFit dataset.aMatrix[amp] = nanMatrixFit dataset.bMatrix[amp] = nanMatrixFit - dataset.covariancesModelNoB[amp] = listNanMatrixFit - dataset.aMatrixNoB[amp] = nanMatrixFit dataset.noiseMatrix[amp] = nanMatrixFit - dataset.noiseMatrixNoB[amp] = nanMatrixFit def errFunc(p, x, y): return ptcFunc(p, x) - y diff --git a/tests/test_defects.py b/tests/test_defects.py index da17243d..7de2419c 100644 --- a/tests/test_defects.py +++ b/tests/test_defects.py @@ -549,7 +549,7 @@ def test_maskBlocks_y_out_of_order_dm38103(self): """A test for maskBlocksIfIntermitentBadPixelsInColumn, y out of order. This test is a variant of - notest_maskBlocks_every_other_pixel_bad_greater_than_threshold with + test_maskBlocks_every_other_pixel_bad_greater_than_threshold with an extra out-of-y-order bad pixel to trigger DM-38103. """ expectedDefects = [Box2I(corner=Point2I(50, 110), dimensions=Extent2I(1, 31))] diff --git a/tests/test_ptc.py b/tests/test_ptc.py index 5c2ddaa2..305e0dc4 100644 --- a/tests/test_ptc.py +++ b/tests/test_ptc.py @@ -149,8 +149,8 @@ def test_covAstier(self): extractConfig.auxiliaryHeaderKeys = ["CCOBCURR", "CCDTEMP"] extractTask = cpPipe.ptc.PhotonTransferCurveExtractTask(config=extractConfig) + # Create solve task config solveConfig = self.defaultConfigSolve - solveConfig.ptcFitType = "FULLCOVARIANCE" # Cut off the low-flux point which is a bad fit, and this # also exercises this functionality and makes the tests # run a lot faster. @@ -162,7 +162,6 @@ def test_covAstier(self): # reliably fit out to a 3x3 covariance matrix. # Improvements will be investigated on DM-46131. solveConfig.maximumRangeCovariancesAstierFullCovFit = 3 - solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=solveConfig) inputGain = self.gain @@ -256,179 +255,176 @@ def test_covAstier(self): # out sorted afterwards. outputCovariancesRev = resultsExtract.outputCovariances[::-1] - # Ensure no warnings re: read noise mismatches are logged. - with self.assertNoLogs(level=logging.WARNING): - resultsSolve = solveTask.run( - outputCovariancesRev, camera=FakeCamera([self.flatExp1.getDetector()]) - ) - - ptc = resultsSolve.outputPtcDataset - # Some expected values for noise matrix, just to check that # it was calculated. - noiseMatrixExpected = np.array( - [[8.99474598, 9.94916264, -27.90587299], - [-2.95079527, -17.11827641, -47.88156244], - [5.24915021, -3.25786165, 26.09634067]], - ) - noiseMatrixNoBExpected = np.array( - [[8.71049338, 12.48584043, -37.06585088], - [-4.80523971, -23.29102809, -66.37815343], - [7.48654766, -4.10168337, 35.64469824]], - ) + expectedNoiseMatrix = { + "FULLCOVARIANCE": np.array( + [[8.99474598, 9.94916264, -27.90587299], + [-2.95079527, -17.11827641, -47.88156244], + [5.24915021, -3.25786165, 26.09634067]], + ), + "FULLCOVARIANCE_NO_B": np.array( + [[8.71049338, 12.48584043, -37.06585088], + [-4.80523971, -23.29102809, -66.37815343], + [7.48654766, -4.10168337, 35.64469824]], + ), + } - for amp in self.ampNames: - self.assertAlmostEqual(ptc.gain[amp], inputGain, places=2) - self.assertFloatsAlmostEqual( - np.asarray(varStandard[amp])[ptc.expIdMask[amp]] / ptc.finalVars[amp][ptc.expIdMask[amp]], - 1.0, - rtol=1e-4, - ) + for fitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]: + solveConfig.ptcFitType = fitType + solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=solveConfig) - # Check that the PTC turnoff is correctly computed. - # This will be different for the C:0,0 amp. - if amp == "C:0,0": - self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[ampName][-3]) - else: - self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[ampName][-1]) - - # Test that all the quantities are correctly ordered and have - # not accidentally been masked. We check every other output ([::2]) - # because these datasets are in pairs of [real, dummy] to - # match the inputs to the extract task. - for i, extractPtc in enumerate(resultsExtract.outputCovariances[::2]): - self.assertFloatsAlmostEqual( - extractPtc.rawExpTimes[ampName][0], - ptc.rawExpTimes[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.rawMeans[ampName][0], - ptc.rawMeans[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.rawVars[ampName][0], - ptc.rawVars[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.photoCharges[ampName][0], - ptc.photoCharges[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.histVars[ampName][0], - ptc.histVars[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.histChi2Dofs[ampName][0], - ptc.histChi2Dofs[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.kspValues[ampName][0], - ptc.kspValues[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.covariances[ampName][0], - ptc.covariances[ampName][i], - ) - self.assertFloatsAlmostEqual( - extractPtc.covariancesSqrtWeights[ampName][0], - ptc.covariancesSqrtWeights[ampName][i], - ) - self.assertFloatsAlmostEqual( - ptc.noiseMatrix[ampName], noiseMatrixExpected, atol=1e-8, rtol=None - ) - self.assertFloatsAlmostEqual( - ptc.noiseMatrixNoB[ampName], - noiseMatrixNoBExpected, - atol=1e-8, - rtol=None, - ) - self.assertFloatsAlmostEqual( - ptc.ampOffsets[ampName], - 0.0, + # Ensure no warnings re: read noise mismatches are logged. + with self.assertNoLogs(level=logging.WARNING): + resultsSolve = solveTask.run( + outputCovariancesRev, camera=FakeCamera([self.flatExp1.getDetector()]) ) + + ptc = resultsSolve.outputPtcDataset + + for amp in self.ampNames: + self.assertAlmostEqual(ptc.gain[amp], inputGain, places=2) self.assertFloatsAlmostEqual( - ptc.noise[ampName], - np.nanmedian(ptc.noiseList[ampName]) * ptc.gain[ampName], - rtol=0.05, - ) - # If the noise error is greater than the noise, something - # is seriously wrong. Possibly some kind of gain application - # mismatch. - self.assertLess( - ptc.noiseErr[ampName], - ptc.noise[ampName], + np.asarray(varStandard[amp])[ptc.expIdMask[amp]] / ptc.finalVars[amp][ptc.expIdMask[amp]], + 1.0, + rtol=1e-4, ) - mask = ptc.getGoodPoints(amp) - - values = ( - ptc.covariancesModel[amp][mask, 0, 0] - ptc.covariances[amp][mask, 0, 0] - ) / ptc.covariancesModel[amp][mask, 0, 0] - np.testing.assert_array_less(np.abs(values), 2e-3) - - values = ( - ptc.covariancesModel[amp][mask, 1, 1] - ptc.covariances[amp][mask, 1, 1] - ) / ptc.covariancesModel[amp][mask, 1, 1] - np.testing.assert_array_less(np.abs(values), 0.3) - - values = ( - ptc.covariancesModel[amp][mask, 1, 2] - ptc.covariances[amp][mask, 1, 2] - ) / ptc.covariancesModel[amp][mask, 1, 2] - np.testing.assert_array_less(np.abs(values), 0.3) - - # And test that the auxiliary values are there and correctly ordered. - self.assertIn('CCOBCURR', ptc.auxValues) - self.assertIn('CCDTEMP', ptc.auxValues) - firstExpIds = np.array([i for i, _ in ptc.inputExpIdPairs['C:0,0']], dtype=np.float64) - self.assertFloatsAlmostEqual(ptc.auxValues['CCOBCURR'], firstExpIds) - self.assertFloatsAlmostEqual(ptc.auxValues['CCDTEMP'], firstExpIds + 1) - - expIdsUsed = ptc.getExpIdsUsed("C:0,0") - # Check that these are the same as the inputs, paired up, with the - # first two (low flux) and final four (outliers, nans) removed. - self.assertTrue( - np.all(expIdsUsed == np.array(expIds).reshape(len(expIds) // 2, 2)[1:-2]) - ) + # Check that the PTC turnoff is correctly computed. + # This will be different for the C:0,0 amp. + if amp == "C:0,0": + self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[amp][-3]) + else: + self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[amp][-1]) + + # Test that all the quantities are correctly ordered and + # have not accidentally been masked. We check every other + # output ([::2]) because these datasets are in pairs of + # [real, dummy] to match the inputs to the extract task. + for i, extractPtc in enumerate(resultsExtract.outputCovariances[::2]): + self.assertFloatsAlmostEqual( + extractPtc.rawExpTimes[ampName][0], + ptc.rawExpTimes[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.rawMeans[ampName][0], + ptc.rawMeans[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.rawVars[ampName][0], + ptc.rawVars[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.photoCharges[ampName][0], + ptc.photoCharges[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.histVars[ampName][0], + ptc.histVars[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.histChi2Dofs[ampName][0], + ptc.histChi2Dofs[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.kspValues[ampName][0], + ptc.kspValues[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.covariances[ampName][0], + ptc.covariances[ampName][i], + ) + self.assertFloatsAlmostEqual( + extractPtc.covariancesSqrtWeights[ampName][0], + ptc.covariancesSqrtWeights[ampName][i], + ) + self.assertFloatsAlmostEqual( + ptc.noiseMatrix[ampName], + expectedNoiseMatrix[fitType], + atol=1e-8, + rtol=None, + ) + self.assertFloatsAlmostEqual( + ptc.ampOffsets[ampName], + 0.0, + ) + self.assertFloatsAlmostEqual( + ptc.noise[ampName], + np.nanmedian(ptc.noiseList[ampName]) * ptc.gain[ampName], + rtol=0.05, + ) + # If the noise error is greater than the noise, + # something is seriously wrong. Possibly some + # kind of gain application mismatch. + self.assertLess( + ptc.noiseErr[ampName], + ptc.noise[ampName], + ) - goodAmps = ptc.getGoodAmps() - self.assertEqual(goodAmps, self.ampNames) + mask = ptc.getGoodPoints(amp) + + values = ( + ptc.covariancesModel[amp][mask, 0, 0] - ptc.covariances[amp][mask, 0, 0] + ) / ptc.covariancesModel[amp][mask, 0, 0] + np.testing.assert_array_less(np.abs(values), 2e-3) + + if ptc.ptcFitType == "FULLCOVARIANCE": + values = ( + ptc.covariancesModel[amp][mask, 0, 1] - ptc.covariances[amp][mask, 0, 1] + ) / ptc.covariancesModel[amp][mask, 0, 1] + np.testing.assert_array_less(np.abs(values), 0.3) + + values = ( + ptc.covariancesModel[amp][mask, 1, 0] - ptc.covariances[amp][mask, 1, 0] + ) / ptc.covariancesModel[amp][mask, 1, 0] + np.testing.assert_array_less(np.abs(values), 0.3) + + # And test that the auxiliary values are there and + # correctly ordered. + self.assertIn('CCOBCURR', ptc.auxValues) + self.assertIn('CCDTEMP', ptc.auxValues) + firstExpIds = np.array([i for i, _ in ptc.inputExpIdPairs['C:0,0']], dtype=np.float64) + self.assertFloatsAlmostEqual(ptc.auxValues['CCOBCURR'], firstExpIds) + self.assertFloatsAlmostEqual(ptc.auxValues['CCDTEMP'], firstExpIds + 1) + + expIdsUsed = ptc.getExpIdsUsed("C:0,0") + # Check that these are the same as the inputs, paired up, with the + # first two (low flux) and final four (outliers, nans) removed. + self.assertTrue( + np.all(expIdsUsed == np.array(expIds).reshape(len(expIds) // 2, 2)[1:-2]) + ) - # Check that every possibly modified field has the same length. - covShape = None - covSqrtShape = None - covModelShape = None - covModelNoBShape = None + goodAmps = ptc.getGoodAmps() + self.assertEqual(goodAmps, self.ampNames) - for ampName in self.ampNames: - if covShape is None: - covShape = ptc.covariances[ampName].shape - covSqrtShape = ptc.covariancesSqrtWeights[ampName].shape - covModelShape = ptc.covariancesModel[ampName].shape - covModelNoBShape = ptc.covariancesModelNoB[ampName].shape - else: - self.assertEqual(ptc.covariances[ampName].shape, covShape) - self.assertEqual( - ptc.covariancesSqrtWeights[ampName].shape, covSqrtShape - ) - self.assertEqual(ptc.covariancesModel[ampName].shape, covModelShape) - self.assertEqual( - ptc.covariancesModelNoB[ampName].shape, covModelNoBShape - ) - # Check if evalPtcModel produces expected values - nanMask = ~np.isnan(ptc.finalMeans[ampName]) - means = ptc.finalMeans[ampName][nanMask] - covModel = ptc.covariancesModel[ampName][nanMask] - covariancesModel = ptc.evalPtcModel(means, setBtoZero=False)[ampName] - self.assertFloatsAlmostEqual(covariancesModel, covModel, atol=1e-12) - # Note that the PhotoTransferCurveDataset does not store the gain - # fit parameter for FULLCOVARIANCE with b=0, so we can't compare - # exactly. - - # And check that this is serializable - with tempfile.NamedTemporaryFile(suffix=".fits") as f: - usedFilename = ptc.writeFits(f.name) - fromFits = PhotonTransferCurveDataset.readFits(usedFilename) - self.assertEqual(fromFits, ptc) + # Check that every possibly modified field has the same length. + covShape = None + covSqrtShape = None + covModelShape = None + + for ampName in self.ampNames: + if covShape is None: + covShape = ptc.covariances[ampName].shape + covSqrtShape = ptc.covariancesSqrtWeights[ampName].shape + covModelShape = ptc.covariancesModel[ampName].shape + else: + self.assertEqual(ptc.covariances[ampName].shape, covShape) + self.assertEqual( + ptc.covariancesSqrtWeights[ampName].shape, covSqrtShape + ) + self.assertEqual(ptc.covariancesModel[ampName].shape, covModelShape) + # Check if evalPtcModel produces expected values + nanMask = ~np.isnan(ptc.finalMeans[ampName]) + means = ptc.finalMeans[ampName][nanMask] + covModel = ptc.covariancesModel[ampName][nanMask] + covariancesModel = ptc.evalPtcModel(means)[ampName] + self.assertFloatsAlmostEqual(covariancesModel, covModel, atol=1e-12) + + # And check that this is serializable + with tempfile.NamedTemporaryFile(suffix=".fits") as f: + usedFilename = ptc.writeFits(f.name) + fromFits = PhotonTransferCurveDataset.readFits(usedFilename) + self.assertEqual(fromFits, ptc) def ptcFitAndCheckPtc( self,