From cd36f615b4f28c67987c4764208701c522273fa9 Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Fri, 10 Feb 2023 10:20:47 -0800 Subject: [PATCH 1/8] Initial push --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 222 +++++++++++--------- 1 file changed, 128 insertions(+), 94 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index 47602a02..f4e52e0a 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -32,6 +32,7 @@ from lsst.ip.isr import PhotonTransferCurveDataset from lsst.ip.isr import IsrTask +from lsst.geom import Box2I, Point2I __all__ = ['PhotonTransferCurveExtractConfig', 'PhotonTransferCurveExtractTask'] @@ -177,6 +178,16 @@ class PhotonTransferCurveExtractConfig(pipeBase.PipelineTaskConfig, 'FULL': 'Second order correction.' } ) + numSubregionsX = pexConfig.Field( + dtype=int, + doc="Number of subregions in each amp in the X direction.", + default=1, + ) + numSubregionsY = pexConfig.Field( + dtype=int, + doc="Number of subregions in each amp in the Y direction.", + default=1, + ) class PhotonTransferCurveExtractTask(pipeBase.PipelineTask): @@ -291,6 +302,11 @@ def run(self, inputExp, inputDims, taskMetadata): amps = detector.getAmplifiers() ampNames = [amp.getName() for amp in amps] + # We have the option to split each amp into a number of sub-regions + # These parameters determine how many subregions we split the amp into. + nSubX = self.config.numSubregionsX + nSubY = self.config.numSubregionsY + # Each amp may have a different min and max ADU signal # specified in the config. maxMeanSignalDict = {ampName: 1e6 for ampName in ampNames} @@ -326,7 +342,7 @@ def run(self, inputExp, inputDims, taskMetadata): partialPtcDatasetList = [] # The number of output references needs to match that of input # references: initialize outputlist with dummy PTC datasets. - for i in range(len(inputDims)): + for i in range(len(inputDims) * nSubX * nSubY): partialPtcDatasetList.append(dummyPtcDataset) if self.config.numEdgeSuspect > 0: @@ -371,99 +387,117 @@ def run(self, inputExp, inputDims, taskMetadata): elif self.config.detectorMeasurementRegion == 'FULL': region = None - # Get masked image regions, masking planes, statistic control - # objects, and clipped means. Calculate once to reuse in - # `measureMeanVarCov` and `getGainFromFlatPair`. - im1Area, im2Area, imStatsCtrl, mu1, mu2 = self.getImageAreasMasksStats(exp1, exp2, - region=region) - - # `measureMeanVarCov` is the function that measures - # the variance and covariances from a region of - # the difference image of two flats at the same - # exposure time. The variable `covAstier` that is - # returned is of the form: - # [(i, j, var (cov[0,0]), cov, npix) for (i,j) in - # {maxLag, maxLag}^2]. - muDiff, varDiff, covAstier = self.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2) - # Estimate the gain from the flat pair - if self.config.doGain: - gain = self.getGainFromFlatPair(im1Area, im2Area, imStatsCtrl, mu1, mu2, - correctionType=self.config.gainCorrectionType, - readNoise=readNoiseDict[ampName]) - else: - gain = np.nan - - # Correction factor for bias introduced by sigma - # clipping. - # Function returns 1/sqrt(varFactor), so it needs - # to be squared. varDiff is calculated via - # afwMath.VARIANCECLIP. - varFactor = sigmaClipCorrection(self.config.nSigmaClipPtc)**2 - varDiff *= varFactor - - expIdMask = True - # Mask data point at this mean signal level if - # the signal, variance, or covariance calculations - # from `measureMeanVarCov` resulted in NaNs. - if np.isnan(muDiff) or np.isnan(varDiff) or (covAstier is None): - self.log.warning("NaN mean or var, or None cov in amp %s in exposure pair %d, %d of " - "detector %d.", ampName, expId1, expId2, detNum) - nAmpsNan += 1 - expIdMask = False - covArray = np.full((1, self.config.maximumRangeCovariancesAstier, - self.config.maximumRangeCovariancesAstier), np.nan) - covSqrtWeights = np.full_like(covArray, np.nan) - - # Mask data point if it is outside of the - # specified mean signal range. - if (muDiff <= minMeanSignalDict[ampName]) or (muDiff >= maxMeanSignalDict[ampName]): - expIdMask = False - - if covAstier is not None: - # Turn the tuples with the measured information - # into covariance arrays. - # covrow: (i, j, var (cov[0,0]), cov, npix) - tupleRows = [(muDiff, varDiff) + covRow + (ampNumber, expTime, - ampName) for covRow in covAstier] - tempStructArray = np.array(tupleRows, dtype=tags) - covArray, vcov, _ = self.makeCovArray(tempStructArray, - self.config.maximumRangeCovariancesAstier) - covSqrtWeights = np.nan_to_num(1./np.sqrt(vcov)) - - # Correct covArray for sigma clipping: - # 1) Apply varFactor twice for the whole covariance matrix - covArray *= varFactor**2 - # 2) But, only once for the variance element of the - # matrix, covArray[0,0] (so divide one factor out). - covArray[0, 0] /= varFactor - - partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff], - rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)], - expIdMask=[expIdMask], covArray=covArray, - covSqrtWeights=covSqrtWeights, gain=gain, - noise=readNoiseDict[ampName]) - # Use location of exp1 to save PTC dataset from (exp1, exp2) pair. - # Below, np.where(expId1 == np.array(inputDims)) returns a tuple - # with a single-element array, so [0][0] - # is necessary to extract the required index. - datasetIndex = np.where(expId1 == np.array(inputDims))[0][0] - # `partialPtcDatasetList` is a list of - # `PhotonTransferCurveDataset` objects. Some of them - # will be dummy datasets (to match length of input - # and output references), and the rest will have - # datasets with the mean signal, variance, and - # covariance measurements at a given exposure - # time. The next ppart of the PTC-measurement - # pipeline, `solve`, will take this list as input, - # and assemble the measurements in the datasets - # in an addecuate manner for fitting a PTC - # model. - partialPtcDataset.updateMetadata(setDate=True, detector=detector) - partialPtcDatasetList[datasetIndex] = partialPtcDataset - - if nAmpsNan == len(ampNames): - msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." - self.log.warning(msg) + # Now split the region into subregions (if nSubX * nSubY > 1) + xStep = int(region.width / nSubX) + yStep = int(region.height / nSubY) + for iSub in range(nSubX): + xmin = region.minX + iSub * xStep + if iSub < (nSubX - 1): + xmax = region.minX + (iSub + 1) * xStep + else: + xmax = region.maxX + for jSub in range(nSubY): + ymin = region.minY + jSub * yStep + if jSub < (nSubY - 1): + ymax = region.minY + (jSub + 1) * yStep + else: + ymax = region.maxY + subRegion = Box2I(minimum=Point2I(xmin,ymin), maximum=Point2I(xmax,ymax)) + + # Get masked image regions, masking planes, statistic control + # objects, and clipped means. Calculate once to reuse in + # `measureMeanVarCov` and `getGainFromFlatPair`. + im1Area, im2Area, imStatsCtrl, mu1, mu2 = self.getImageAreasMasksStats(exp1, exp2, + region=subRegion) + + # `measureMeanVarCov` is the function that measures + # the variance and covariances from a region of + # the difference image of two flats at the same + # exposure time. The variable `covAstier` that is + # returned is of the form: + # [(i, j, var (cov[0,0]), cov, npix) for (i,j) in + # {maxLag, maxLag}^2]. + muDiff, varDiff, covAstier = self.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2) + # Estimate the gain from the flat pair + if self.config.doGain: + gain = self.getGainFromFlatPair(im1Area, im2Area, imStatsCtrl, mu1, mu2, + correctionType=self.config.gainCorrectionType, + readNoise=readNoiseDict[ampName]) + else: + gain = np.nan + + # Correction factor for bias introduced by sigma + # clipping. + # Function returns 1/sqrt(varFactor), so it needs + # to be squared. varDiff is calculated via + # afwMath.VARIANCECLIP. + varFactor = sigmaClipCorrection(self.config.nSigmaClipPtc)**2 + varDiff *= varFactor + + expIdMask = True + # Mask data point at this mean signal level if + # the signal, variance, or covariance calculations + # from `measureMeanVarCov` resulted in NaNs. + if np.isnan(muDiff) or np.isnan(varDiff) or (covAstier is None): + self.log.warning("NaN mean or var, or None cov in amp %s in exposure pair %d, %d of " + "detector %d.", ampName, expId1, expId2, detNum) + nAmpsNan += 1 + expIdMask = False + covArray = np.full((1, self.config.maximumRangeCovariancesAstier, + self.config.maximumRangeCovariancesAstier), np.nan) + covSqrtWeights = np.full_like(covArray, np.nan) + + # Mask data point if it is outside of the + # specified mean signal range. + if (muDiff <= minMeanSignalDict[ampName]) or (muDiff >= maxMeanSignalDict[ampName]): + expIdMask = False + + if covAstier is not None: + # Turn the tuples with the measured information + # into covariance arrays. + # covrow: (i, j, var (cov[0,0]), cov, npix) + tupleRows = [(muDiff, varDiff) + covRow + (ampNumber, expTime, + ampName) for covRow in covAstier] + tempStructArray = np.array(tupleRows, dtype=tags) + covArray, vcov, _ = self.makeCovArray(tempStructArray, + self.config.maximumRangeCovariancesAstier) + covSqrtWeights = np.nan_to_num(1./np.sqrt(vcov)) + + # Correct covArray for sigma clipping: + # 1) Apply varFactor twice for the whole covariance matrix + covArray *= varFactor**2 + # 2) But, only once for the variance element of the + # matrix, covArray[0,0] (so divide one factor out). + covArray[0, 0] /= varFactor + + partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff], + rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)], + expIdMask=[expIdMask], covArray=covArray, + covSqrtWeights=covSqrtWeights, gain=gain, + noise=readNoiseDict[ampName]) + # Use location of exp1 to save PTC dataset from (exp1, exp2) pair. + # Below, np.where(expId1 == np.array(inputDims)) returns a tuple + # with a single-element array, so [0][0] + # is necessary to extract the required index. + #datasetIndex = np.where(expId1 == np.array(inputDims))[0][0] + # `partialPtcDatasetList` is a list of + # `PhotonTransferCurveDataset` objects. Some of them + # will be dummy datasets (to match length of input + # and output references), and the rest will have + # datasets with the mean signal, variance, and + # covariance measurements at a given exposure + # time. The next ppart of the PTC-measurement + # pipeline, `solve`, will take this list as input, + # and assemble the measurements in the datasets + # in an addecuate manner for fitting a PTC + # model. + partialPtcDataset.updateMetadata(setDate=True, detector=detector) + #partialPtcDatasetList[datasetIndex] = partialPtcDataset + partialPtcDatasetList.append(partialPtcDataset) + + if nAmpsNan == len(ampNames): + msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." + self.log.warning(msg) return pipeBase.Struct( outputCovariances=partialPtcDatasetList, ) From 5259bc83c52500f80673639f5b209572a255542a Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sat, 11 Feb 2023 03:19:30 -0800 Subject: [PATCH 2/8] Fixing indexing problem --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index f4e52e0a..c2d66223 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -479,7 +479,10 @@ def run(self, inputExp, inputDims, taskMetadata): # Below, np.where(expId1 == np.array(inputDims)) returns a tuple # with a single-element array, so [0][0] # is necessary to extract the required index. - #datasetIndex = np.where(expId1 == np.array(inputDims))[0][0] + # After extracting the index, we add an integer to it + # to account for the multiple subregions + datasetIndex = np.where(expId1 == np.array(inputDims))[0][0] + datasetIndex += iSub * nSubY + jSub # `partialPtcDatasetList` is a list of # `PhotonTransferCurveDataset` objects. Some of them # will be dummy datasets (to match length of input @@ -489,11 +492,11 @@ def run(self, inputExp, inputDims, taskMetadata): # time. The next ppart of the PTC-measurement # pipeline, `solve`, will take this list as input, # and assemble the measurements in the datasets - # in an addecuate manner for fitting a PTC + # in an adecuate manner for fitting a PTC # model. partialPtcDataset.updateMetadata(setDate=True, detector=detector) - #partialPtcDatasetList[datasetIndex] = partialPtcDataset - partialPtcDatasetList.append(partialPtcDataset) + partialPtcDatasetList[datasetIndex] = partialPtcDataset + if nAmpsNan == len(ampNames): msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." From c711f3469b643b937f9c8ce63beef9456d11a70f Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sat, 11 Feb 2023 03:33:03 -0800 Subject: [PATCH 3/8] Another indexing correction --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index c2d66223..5a06f9d2 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -479,10 +479,10 @@ def run(self, inputExp, inputDims, taskMetadata): # Below, np.where(expId1 == np.array(inputDims)) returns a tuple # with a single-element array, so [0][0] # is necessary to extract the required index. - # After extracting the index, we add an integer to it + # After extracting the index, we modify it # to account for the multiple subregions - datasetIndex = np.where(expId1 == np.array(inputDims))[0][0] - datasetIndex += iSub * nSubY + jSub + expIndex = np.where(expId1 == np.array(inputDims))[0][0] + datasetIndex = expIndex * nSubX * nSubY + iSub * nSubY + jSub # `partialPtcDatasetList` is a list of # `PhotonTransferCurveDataset` objects. Some of them # will be dummy datasets (to match length of input From a2a0daf9dfbf93e7c5bff93a8b06da365879f5ed Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sat, 11 Feb 2023 04:17:25 -0800 Subject: [PATCH 4/8] Correcting misunderstanding on ptcDatasetList --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 53 ++++++++++----------- 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index 5a06f9d2..9517c489 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -268,6 +268,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs['inputExp'] = arrangeFlatsByExpId(inputs['inputExp'], inputs['inputDims']) outputs = self.run(**inputs) + print("In runQuantum", outputs) butlerQC.put(outputs, outputRefs) def run(self, inputExp, inputDims, taskMetadata): @@ -342,7 +343,7 @@ def run(self, inputExp, inputDims, taskMetadata): partialPtcDatasetList = [] # The number of output references needs to match that of input # references: initialize outputlist with dummy PTC datasets. - for i in range(len(inputDims) * nSubX * nSubY): + for i in range(len(inputDims)): partialPtcDatasetList.append(dummyPtcDataset) if self.config.numEdgeSuspect > 0: @@ -475,32 +476,30 @@ def run(self, inputExp, inputDims, taskMetadata): expIdMask=[expIdMask], covArray=covArray, covSqrtWeights=covSqrtWeights, gain=gain, noise=readNoiseDict[ampName]) - # Use location of exp1 to save PTC dataset from (exp1, exp2) pair. - # Below, np.where(expId1 == np.array(inputDims)) returns a tuple - # with a single-element array, so [0][0] - # is necessary to extract the required index. - # After extracting the index, we modify it - # to account for the multiple subregions - expIndex = np.where(expId1 == np.array(inputDims))[0][0] - datasetIndex = expIndex * nSubX * nSubY + iSub * nSubY + jSub - # `partialPtcDatasetList` is a list of - # `PhotonTransferCurveDataset` objects. Some of them - # will be dummy datasets (to match length of input - # and output references), and the rest will have - # datasets with the mean signal, variance, and - # covariance measurements at a given exposure - # time. The next ppart of the PTC-measurement - # pipeline, `solve`, will take this list as input, - # and assemble the measurements in the datasets - # in an adecuate manner for fitting a PTC - # model. - partialPtcDataset.updateMetadata(setDate=True, detector=detector) - partialPtcDatasetList[datasetIndex] = partialPtcDataset - - - if nAmpsNan == len(ampNames): - msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." - self.log.warning(msg) + # Use location of exp1 to save PTC dataset from (exp1, exp2) pair. + # Below, np.where(expId1 == np.array(inputDims)) returns a tuple + # with a single-element array, so [0][0] + # is necessary to extract the required index. + datasetIndex = np.where(expId1 == np.array(inputDims))[0][0] + # `partialPtcDatasetList` is a list of + # `PhotonTransferCurveDataset` objects. Some of them + # will be dummy datasets (to match length of input + # and output references), and the rest will have + # datasets with the mean signal, variance, and + # covariance measurements at a given exposure + # time. The next ppart of the PTC-measurement + # pipeline, `solve`, will take this list as input, + # and assemble the measurements in the datasets + # in an adecuate manner for fitting a PTC + # model. + partialPtcDataset.updateMetadata(setDate=True, detector=detector) + partialPtcDatasetList[datasetIndex] = partialPtcDataset + + + if nAmpsNan == len(ampNames): + msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." + self.log.warning(msg) + print("Finished covariances", len(inputDims), len(partialPtcDatasetList)) return pipeBase.Struct( outputCovariances=partialPtcDatasetList, ) From 3d2dadec5c2332446356a5141199d7fcce03d82e Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sat, 11 Feb 2023 10:17:46 -0800 Subject: [PATCH 5/8] Lots of changes to deal with the larger arrays. Still not working --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 5 +++-- python/lsst/cp/pipe/ptc/cpSolvePtcTask.py | 15 ++++++++------- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index 9517c489..5883e727 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -268,7 +268,6 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs['inputExp'] = arrangeFlatsByExpId(inputs['inputExp'], inputs['inputDims']) outputs = self.run(**inputs) - print("In runQuantum", outputs) butlerQC.put(outputs, outputRefs) def run(self, inputExp, inputDims, taskMetadata): @@ -470,7 +469,9 @@ def run(self, inputExp, inputDims, taskMetadata): # 2) But, only once for the variance element of the # matrix, covArray[0,0] (so divide one factor out). covArray[0, 0] /= varFactor - + print(f"Amp:{ampName}, iSub={iSub}, jSub={jSub}, loading a covariance set") + print(f"In cpExtract, type(covArray)={type(covArray)}, shape = {np.array(covArray).shape}") + print(f"In cpExtract, type(covSqrtWeights)={type(covSqrtWeights)}, shape = {np.array(covSqrtWeights).shape}") partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff], rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)], expIdMask=[expIdMask], covArray=covArray, diff --git a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py index 559340a8..2d806019 100644 --- a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py @@ -226,24 +226,25 @@ def run(self, inputCovariances, camera=None, detId=0): for ampName in ampNames: datasetPtc.inputExpIdPairs[ampName].append(partialPtcDataset.inputExpIdPairs[ampName]) if type(partialPtcDataset.rawExpTimes[ampName]) is list: - datasetPtc.rawExpTimes[ampName].append(partialPtcDataset.rawExpTimes[ampName][0]) + datasetPtc.rawExpTimes[ampName] += partialPtcDataset.rawExpTimes[ampName] else: datasetPtc.rawExpTimes[ampName].append(partialPtcDataset.rawExpTimes[ampName]) if type(partialPtcDataset.rawMeans[ampName]) is list: - datasetPtc.rawMeans[ampName].append(partialPtcDataset.rawMeans[ampName][0]) + print(f"Amp:{ampName}, rawMeans is a list of length {len(partialPtcDataset.rawMeans[ampName])}") + datasetPtc.rawMeans[ampName]+= partialPtcDataset.rawMeans[ampName] else: datasetPtc.rawMeans[ampName].append(partialPtcDataset.rawMeans[ampName]) if type(partialPtcDataset.rawVars[ampName]) is list: - datasetPtc.rawVars[ampName].append(partialPtcDataset.rawVars[ampName][0]) + datasetPtc.rawVars[ampName] += partialPtcDataset.rawVars[ampName] else: datasetPtc.rawVars[ampName].append(partialPtcDataset.rawVars[ampName]) if type(partialPtcDataset.expIdMask[ampName]) is list: - datasetPtc.expIdMask[ampName].append(partialPtcDataset.expIdMask[ampName][0]) + datasetPtc.expIdMask[ampName] += partialPtcDataset.expIdMask[ampName] else: datasetPtc.expIdMask[ampName].append(partialPtcDataset.expIdMask[ampName]) - datasetPtc.covariances[ampName].append(np.array(partialPtcDataset.covariances[ampName][0])) - datasetPtc.covariancesSqrtWeights[ampName].append( - np.array(partialPtcDataset.covariancesSqrtWeights[ampName][0])) + print(f"Amp:{ampName}, loading covariances") + datasetPtc.covariances[ampName] += partialPtcDataset.covariances[ampName] + datasetPtc.covariancesSqrtWeights[ampName] += partialPtcDataset.covariancesSqrtWeights[ampName] # Sort arrays that are filled so far in the final dataset by # rawMeans index for ampName in ampNames: From 0bca008c745d3a9971126da453f76c5bca5d78ec Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sun, 12 Feb 2023 03:51:04 -0800 Subject: [PATCH 6/8] Changes to deal with larger arrays. cpExtract now runs, but still fails in cpSolve --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 6 +++--- python/lsst/cp/pipe/ptc/cpSolvePtcTask.py | 11 ++++++++++- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index 5883e727..3dcc496d 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -469,9 +469,9 @@ def run(self, inputExp, inputDims, taskMetadata): # 2) But, only once for the variance element of the # matrix, covArray[0,0] (so divide one factor out). covArray[0, 0] /= varFactor - print(f"Amp:{ampName}, iSub={iSub}, jSub={jSub}, loading a covariance set") - print(f"In cpExtract, type(covArray)={type(covArray)}, shape = {np.array(covArray).shape}") - print(f"In cpExtract, type(covSqrtWeights)={type(covSqrtWeights)}, shape = {np.array(covSqrtWeights).shape}") + #print(f"Amp:{ampName}, iSub={iSub}, jSub={jSub}, loading a covariance set") + #print(f"In cpExtract, type(covArray)={type(covArray)}, shape = {np.array(covArray).shape}") + #print(f"In cpExtract, type(covSqrtWeights)={type(covSqrtWeights)}, shape = {np.array(covSqrtWeights).shape}") partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff], rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)], expIdMask=[expIdMask], covArray=covArray, diff --git a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py index 2d806019..3cafaf9a 100644 --- a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py @@ -224,7 +224,7 @@ def run(self, inputCovariances, camera=None, detId=0): if partialPtcDataset.ptcFitType == 'DUMMY': continue for ampName in ampNames: - datasetPtc.inputExpIdPairs[ampName].append(partialPtcDataset.inputExpIdPairs[ampName]) + datasetPtc.inputExpIdPairs[ampName] += partialPtcDataset.inputExpIdPairs[ampName] if type(partialPtcDataset.rawExpTimes[ampName]) is list: datasetPtc.rawExpTimes[ampName] += partialPtcDataset.rawExpTimes[ampName] else: @@ -248,7 +248,15 @@ def run(self, inputCovariances, camera=None, detId=0): # Sort arrays that are filled so far in the final dataset by # rawMeans index for ampName in ampNames: + print(f"amp:{ampName}, ExpIdPairs size = {np.array(datasetPtc.inputExpIdPairs[ampName]).shape}") + print(f"amp:{ampName}, rawMeans size = {np.array(datasetPtc.rawMeans[ampName]).shape}") + print(f"amp:{ampName}, covariances size = {np.array(datasetPtc.covariances[ampName]).shape}") + print(f"amp:{ampName}, ExpIdMask size = {np.array(datasetPtc.expIdMask[ampName]).shape}") + + + index = np.argsort(np.ravel(np.array(datasetPtc.rawMeans[ampName]))) + print(f"amp:{ampName}, index={index}") datasetPtc.inputExpIdPairs[ampName] = np.array(datasetPtc.inputExpIdPairs[ampName])[index] datasetPtc.rawExpTimes[ampName] = np.array(datasetPtc.rawExpTimes[ampName])[index] datasetPtc.rawMeans[ampName] = np.array(datasetPtc.rawMeans[ampName])[index] @@ -915,6 +923,7 @@ def errFunc(p, x, y): continue dataset.expIdMask[ampName] = np.array(dataset.expIdMask[ampName]) # store the final mask + print(f"Masking expIdMask shape = {dataset.expIdMask[ampName].shape}, lenMask={len(mask)}") if len(dataset.expIdMask[ampName]): dataset.expIdMask[ampName] &= mask # bitwise_and if there is already a mask else: From 5fd25ba86d2750059517105334272c84e80b46f8 Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sun, 12 Feb 2023 04:28:33 -0800 Subject: [PATCH 7/8] More fixes and clean-up. Seems to be working now --- python/lsst/cp/pipe/ptc/cpExtractPtcTask.py | 6 +----- python/lsst/cp/pipe/ptc/cpSolvePtcTask.py | 12 +----------- 2 files changed, 2 insertions(+), 16 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index 3dcc496d..059784c0 100644 --- a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py @@ -387,7 +387,7 @@ def run(self, inputExp, inputDims, taskMetadata): elif self.config.detectorMeasurementRegion == 'FULL': region = None - # Now split the region into subregions (if nSubX * nSubY > 1) + # Now split the region into subregions xStep = int(region.width / nSubX) yStep = int(region.height / nSubY) for iSub in range(nSubX): @@ -469,9 +469,6 @@ def run(self, inputExp, inputDims, taskMetadata): # 2) But, only once for the variance element of the # matrix, covArray[0,0] (so divide one factor out). covArray[0, 0] /= varFactor - #print(f"Amp:{ampName}, iSub={iSub}, jSub={jSub}, loading a covariance set") - #print(f"In cpExtract, type(covArray)={type(covArray)}, shape = {np.array(covArray).shape}") - #print(f"In cpExtract, type(covSqrtWeights)={type(covSqrtWeights)}, shape = {np.array(covSqrtWeights).shape}") partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff], rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)], expIdMask=[expIdMask], covArray=covArray, @@ -500,7 +497,6 @@ def run(self, inputExp, inputDims, taskMetadata): if nAmpsNan == len(ampNames): msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." self.log.warning(msg) - print("Finished covariances", len(inputDims), len(partialPtcDatasetList)) return pipeBase.Struct( outputCovariances=partialPtcDatasetList, ) diff --git a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py index 3cafaf9a..774c5d1a 100644 --- a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py @@ -223,6 +223,7 @@ def run(self, inputCovariances, camera=None, detId=0): # Ignore dummy datasets if partialPtcDataset.ptcFitType == 'DUMMY': continue + # Modifications here for the subregions to append new data to the old for ampName in ampNames: datasetPtc.inputExpIdPairs[ampName] += partialPtcDataset.inputExpIdPairs[ampName] if type(partialPtcDataset.rawExpTimes[ampName]) is list: @@ -230,7 +231,6 @@ def run(self, inputCovariances, camera=None, detId=0): else: datasetPtc.rawExpTimes[ampName].append(partialPtcDataset.rawExpTimes[ampName]) if type(partialPtcDataset.rawMeans[ampName]) is list: - print(f"Amp:{ampName}, rawMeans is a list of length {len(partialPtcDataset.rawMeans[ampName])}") datasetPtc.rawMeans[ampName]+= partialPtcDataset.rawMeans[ampName] else: datasetPtc.rawMeans[ampName].append(partialPtcDataset.rawMeans[ampName]) @@ -242,21 +242,12 @@ def run(self, inputCovariances, camera=None, detId=0): datasetPtc.expIdMask[ampName] += partialPtcDataset.expIdMask[ampName] else: datasetPtc.expIdMask[ampName].append(partialPtcDataset.expIdMask[ampName]) - print(f"Amp:{ampName}, loading covariances") datasetPtc.covariances[ampName] += partialPtcDataset.covariances[ampName] datasetPtc.covariancesSqrtWeights[ampName] += partialPtcDataset.covariancesSqrtWeights[ampName] # Sort arrays that are filled so far in the final dataset by # rawMeans index for ampName in ampNames: - print(f"amp:{ampName}, ExpIdPairs size = {np.array(datasetPtc.inputExpIdPairs[ampName]).shape}") - print(f"amp:{ampName}, rawMeans size = {np.array(datasetPtc.rawMeans[ampName]).shape}") - print(f"amp:{ampName}, covariances size = {np.array(datasetPtc.covariances[ampName]).shape}") - print(f"amp:{ampName}, ExpIdMask size = {np.array(datasetPtc.expIdMask[ampName]).shape}") - - - index = np.argsort(np.ravel(np.array(datasetPtc.rawMeans[ampName]))) - print(f"amp:{ampName}, index={index}") datasetPtc.inputExpIdPairs[ampName] = np.array(datasetPtc.inputExpIdPairs[ampName])[index] datasetPtc.rawExpTimes[ampName] = np.array(datasetPtc.rawExpTimes[ampName])[index] datasetPtc.rawMeans[ampName] = np.array(datasetPtc.rawMeans[ampName])[index] @@ -923,7 +914,6 @@ def errFunc(p, x, y): continue dataset.expIdMask[ampName] = np.array(dataset.expIdMask[ampName]) # store the final mask - print(f"Masking expIdMask shape = {dataset.expIdMask[ampName].shape}, lenMask={len(mask)}") if len(dataset.expIdMask[ampName]): dataset.expIdMask[ampName] &= mask # bitwise_and if there is already a mask else: From fef1f50a5f652ebf5b69f30df8440fa4855b279f Mon Sep 17 00:00:00 2001 From: craiglagegit Date: Sun, 12 Feb 2023 14:46:41 -0800 Subject: [PATCH 8/8] Pushing changes to getInitialGoodPoints --- python/lsst/cp/pipe/ptc/cpSolvePtcTask.py | 56 ++++++++--------------- 1 file changed, 19 insertions(+), 37 deletions(-) diff --git a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py index 774c5d1a..1cfb1783 100644 --- a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py +++ b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py @@ -28,6 +28,9 @@ from scipy.signal import fftconvolve from scipy.optimize import least_squares +from scipy.interpolate import interp1d +from scipy.signal import savgol_filter + from itertools import groupby from operator import itemgetter @@ -701,14 +704,6 @@ def _getInitialGoodPoints(means, variances, minVarPivotSearch, consecutivePoints Input array with mean signal values. variances : `numpy.array` Input array with variances at each mean value. - minVarPivotSearch : `float` - The variance (in ADU^2), above which, the point - of decreasing variance should be sought. - consecutivePointsVarDecreases : `int` - Required number of consecutive points/fluxes - in the PTC where the variance - decreases in order to find a first - estimate of the PTC turn-off. Returns ------ @@ -719,37 +714,24 @@ def _getInitialGoodPoints(means, variances, minVarPivotSearch, consecutivePoints Notes ----- Eliminate points beyond which the variance decreases. + This algorithm has been problematic, so I (Lage -12-Feb-23) + have rewritten it to smooth the (mean,var) points + and then find the max variance of the smoothed points """ goodPoints = np.ones_like(means, dtype=bool) - # Variances are sorted and should monotonically increase - pivotList = np.where(np.array(np.diff(variances)) < 0)[0] - if len(pivotList) > 0: - # For small values, sometimes the variance decreases slightly - # Only look when var > self.config.minVarPivotSearch - pivotList = [p for p in pivotList if variances[p] > minVarPivotSearch] - # Require that the varince decreases during - # consecutivePointsVarDecreases - # consecutive points. This will give a first - # estimate of the PTC turn-off, which - # may be updated (reduced) further in the code. - if len(pivotList) > 1: - # enumerate(pivotList) creates tuples (index, value), for - # each value in pivotList. The lambda function subtracts - # each value from the index. - # groupby groups elements by equal key value. - for k, g in groupby(enumerate(pivotList), lambda x: x[0]-x[1]): - group = (map(itemgetter(1), g)) - # Form groups of consecute values from pivotList - group = list(map(int, group)) - # values in pivotList are indices where np.diff(variances) - # is negative, i.e., where the variance starts decreasing. - # Find the first group of consecutive numbers when - # variance decreases. - if len(group) >= consecutivePointsVarDecreases: - pivotIndex = np.min(group) - goodPoints[pivotIndex+1:] = False - break - + x = np.array(means) + y = np.array(variances) + xx = np.linspace(x.min(),x.max(), 1000) + + # interpolate + smooth + itp = interp1d(x,y, kind='linear') + window_size, poly_order = 101, 3 + yy_sg = savgol_filter(itp(xx), window_size, poly_order) + max_var_index = np.argmax(yy_sg) + ptc_turnoff = xx[max_var_index] + for i, mean in enumerate(means): + if mean > ptc_turnoff: + goodPoints[i] = False return goodPoints def _makeZeroSafe(self, array, substituteValue=1e-9):