From 172b5a6142f72fe7fda002ac0fbc3efea609cdf8 Mon Sep 17 00:00:00 2001 From: Alex Broughton Date: Wed, 9 Oct 2024 14:30:35 -0700 Subject: [PATCH 1/4] Add FULLCOVARIANCENOB PTC fit type --- python/lsst/cp/pipe/ptc/cpPtcSolve.py | 89 ++++++++++++--------------- 1 file changed, 40 insertions(+), 49 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpPtcSolve.py b/python/lsst/cp/pipe/ptc/cpPtcSolve.py index 3ed70a90..68ccd229 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)" + "FULLCOVARIANCENOB": "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", "FULLCOVARIANCENOB"]:: # 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", "FULLCOVARIANCENOB"]: # 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 @@ -611,24 +609,30 @@ def fitDataFullCovariance(self, dataset): # Fit full model (Eq. 20 of Astier+19) and same model with # b=0 (c=0 in this code). 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 + + # Pick the correct full covariance model function + model = self.funcFullCovarianceModel + if dataset.ptcFitType == "FULLCOVARIANCENOB": + model = self.funcFullCovarianceModelNoB + + fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []} + params, paramsErr, _ = fitLeastSq( + pInit, + muAtAmpMasked, + covAtAmpForFitMasked.ravel(), + model, + 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['a'] = a + fitResults['c'] = c + fitResults['noiseMatrix'] = noiseMatrix + fitResults['gain'] = gain + fitResults['paramsErr'] = paramsErr # Put the information in the PTC dataset @@ -646,33 +650,23 @@ 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'] ) 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 @@ -1090,10 +1084,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 From 08e63ab8e6185f8c3344e7633a9e53945bc32134 Mon Sep 17 00:00:00 2001 From: Alex Broughton Date: Wed, 9 Oct 2024 14:45:37 -0700 Subject: [PATCH 2/4] Add tests for PTC fit type FULLCOVARIANCENOB --- Untitled.ipynb | 6 + python/lsst/cp/pipe/ptc/cpPtcSolve.py | 15 +- tests/test_ptc.py | 322 +++++++++++++------------- 3 files changed, 173 insertions(+), 170 deletions(-) create mode 100644 Untitled.ipynb diff --git a/Untitled.ipynb b/Untitled.ipynb new file mode 100644 index 00000000..363fcab7 --- /dev/null +++ b/Untitled.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/python/lsst/cp/pipe/ptc/cpPtcSolve.py b/python/lsst/cp/pipe/ptc/cpPtcSolve.py index 68ccd229..c5bffda1 100644 --- a/python/lsst/cp/pipe/ptc/cpPtcSolve.py +++ b/python/lsst/cp/pipe/ptc/cpPtcSolve.py @@ -82,7 +82,7 @@ class PhotonTransferCurveSolveConfig(pipeBase.PipelineTaskConfig, allowed={ "POLYNOMIAL": "n-degree polynomial (use 'polynomialFitDegree' to set 'n').", "EXPAPPROXIMATION": "Approximation in Astier+19 (Eq. 16).", - "FULLCOVARIANCENOB": "Full covariances model in Astier+19 (Eq. 15)", + "FULLCOVARIANCE_NO_B": "Full covariances model in Astier+19 (Eq. 15)", "FULLCOVARIANCE": "Full covariances model in Astier+19 (Eq. 20)", } ) @@ -400,7 +400,7 @@ def run(self, inputCovariances, camera=None, detId=0): index = np.argsort(detectorMeans) datasetPtc.sort(index) - if self.config.ptcFitType in ["FULLCOVARIANCE", "FULLCOVARIANCENOB"]:: + 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 @@ -486,7 +486,7 @@ def fitMeasurementsToModel(self, dataset): `PhotonTransferCurveDatase`. """ fitType = dataset.ptcFitType - if fitType in ["FULLCOVARIANCE", "FULLCOVARIANCENOB"]: + 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 @@ -609,10 +609,10 @@ def fitDataFullCovariance(self, dataset): # Fit full model (Eq. 20 of Astier+19) and same model with # b=0 (c=0 in this code). pInit = np.concatenate((a0.ravel(), c0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None) - + # Pick the correct full covariance model function model = self.funcFullCovarianceModel - if dataset.ptcFitType == "FULLCOVARIANCENOB": + if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": model = self.funcFullCovarianceModelNoB fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []} @@ -623,6 +623,8 @@ def fitDataFullCovariance(self, dataset): model, weightsY=covSqrtWeightsAtAmpForFitMasked.ravel() ) + if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": + import IPython; IPython.embed() a = params[:lenParams].reshape((matrixSideFit, matrixSideFit)) c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit)) @@ -653,7 +655,8 @@ def fitDataFullCovariance(self, dataset): fitResults['a'], fitResults['c'], fitResults['noiseMatrix'], - fitResults['gain'] + fitResults['gain'], + setBtoZero=dataset.ptcFitType == "FULLCOVARIANCE_NO_B", ) dataset.covariancesSqrtWeights[ampName] = covSqrtWeightsAtAmp dataset.aMatrix[ampName] = fitResults['a'] diff --git a/tests/test_ptc.py b/tests/test_ptc.py index 5c2ddaa2..d4dd5ff4 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,174 @@ 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[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], + 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) + + 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 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, From d8428c191abe2e15e7bb25c65c959fe7ca8339b0 Mon Sep 17 00:00:00 2001 From: Alex Broughton Date: Thu, 10 Oct 2024 17:18:12 -0700 Subject: [PATCH 3/4] Take out the b parameter from the b=0 full covariance fit --- Untitled.ipynb | 6 ------ python/lsst/cp/pipe/ptc/cpPtcSolve.py | 23 +++++++++++++++-------- tests/test_ptc.py | 2 +- 3 files changed, 16 insertions(+), 15 deletions(-) delete mode 100644 Untitled.ipynb diff --git a/Untitled.ipynb b/Untitled.ipynb deleted file mode 100644 index 363fcab7..00000000 --- a/Untitled.ipynb +++ /dev/null @@ -1,6 +0,0 @@ -{ - "cells": [], - "metadata": {}, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/python/lsst/cp/pipe/ptc/cpPtcSolve.py b/python/lsst/cp/pipe/ptc/cpPtcSolve.py index c5bffda1..0487ccf3 100644 --- a/python/lsst/cp/pipe/ptc/cpPtcSolve.py +++ b/python/lsst/cp/pipe/ptc/cpPtcSolve.py @@ -606,25 +606,31 @@ 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) + # Initialize empty results dictionary + fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []} + # Pick the correct full covariance model function model = self.funcFullCovarianceModel if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": model = self.funcFullCovarianceModelNoB + pInit = np.concatenate((a0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None) - fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []} params, paramsErr, _ = fitLeastSq( pInit, muAtAmpMasked, covAtAmpForFitMasked.ravel(), model, - weightsY=covSqrtWeightsAtAmpForFitMasked.ravel() + weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(), ) + if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": - import IPython; IPython.embed() + 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)) @@ -638,7 +644,7 @@ def fitDataFullCovariance(self, dataset): # 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 @@ -789,8 +795,9 @@ 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) + # params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) + noiseMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) gain = params[-1] return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=True).flatten() diff --git a/tests/test_ptc.py b/tests/test_ptc.py index d4dd5ff4..75382053 100644 --- a/tests/test_ptc.py +++ b/tests/test_ptc.py @@ -270,7 +270,7 @@ def test_covAstier(self): ), } - for fitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]: + for fitType in ["FULLCOVARIANCE_NO_B"]: solveConfig.ptcFitType = fitType solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=solveConfig) From b7a62c33efc40c05425b71da640eb86504286467 Mon Sep 17 00:00:00 2001 From: Alex Broughton Date: Fri, 11 Oct 2024 16:54:23 -0700 Subject: [PATCH 4/4] Update test_covAstier with new numbers --- pipelines/LSSTComCam/cpCti.yaml | 4 +++ python/lsst/cp/pipe/ptc/cpPtcSolve.py | 19 +++++------ tests/test_defects.py | 2 +- tests/test_ptc.py | 48 ++++++++++++++------------- 4 files changed, 39 insertions(+), 34 deletions(-) create mode 100644 pipelines/LSSTComCam/cpCti.yaml 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 0487ccf3..a253037b 100644 --- a/python/lsst/cp/pipe/ptc/cpPtcSolve.py +++ b/python/lsst/cp/pipe/ptc/cpPtcSolve.py @@ -613,17 +613,17 @@ def fitDataFullCovariance(self, dataset): fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []} # Pick the correct full covariance model function - model = self.funcFullCovarianceModel + ptcModel = self.funcFullCovarianceModel if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": - model = self.funcFullCovarianceModelNoB + ptcModel = self.funcFullCovarianceModelNoB pInit = np.concatenate((a0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None) params, paramsErr, _ = fitLeastSq( - pInit, - muAtAmpMasked, - covAtAmpForFitMasked.ravel(), - model, - weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(), + pInit, + muAtAmpMasked, + covAtAmpForFitMasked.ravel(), + ptcModel, + weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(), ) if dataset.ptcFitType == "FULLCOVARIANCE_NO_B": @@ -662,7 +662,7 @@ def fitDataFullCovariance(self, dataset): fitResults['c'], fitResults['noiseMatrix'], fitResults['gain'], - setBtoZero=dataset.ptcFitType == "FULLCOVARIANCE_NO_B", + setBtoZero=(dataset.ptcFitType == "FULLCOVARIANCE_NO_B"), ) dataset.covariancesSqrtWeights[ampName] = covSqrtWeightsAtAmp dataset.aMatrix[ampName] = fitResults['a'] @@ -782,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) @@ -796,7 +796,6 @@ def funcFullCovarianceModelNoB(self, params, x): lenParams = matrixSideFit*matrixSideFit aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit)) cMatrix = np.zeros_like(aMatrix) - # params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) noiseMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit)) gain = params[-1] 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 75382053..305e0dc4 100644 --- a/tests/test_ptc.py +++ b/tests/test_ptc.py @@ -260,17 +260,17 @@ def test_covAstier(self): expectedNoiseMatrix = { "FULLCOVARIANCE": np.array( [[8.99474598, 9.94916264, -27.90587299], - [-2.95079527, -17.11827641, -47.88156244], - [5.24915021, -3.25786165, 26.09634067]], + [-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]], + [-4.80523971, -23.29102809, -66.37815343], + [7.48654766, -4.10168337, 35.64469824]], ), } - for fitType in ["FULLCOVARIANCE_NO_B"]: + for fitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]: solveConfig.ptcFitType = fitType solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=solveConfig) @@ -293,14 +293,14 @@ def test_covAstier(self): # 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]) + self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[amp][-3]) else: - self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[ampName][-1]) + 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. + # 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], @@ -353,9 +353,9 @@ def test_covAstier(self): 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. + # 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], @@ -368,17 +368,19 @@ def test_covAstier(self): ) / 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) + 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, 2] - ptc.covariances[amp][mask, 1, 2] - ) / ptc.covariancesModel[amp][mask, 1, 2] - 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. + # 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)