diff --git a/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py index 47602a02d..059784c05 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} @@ -371,77 +387,93 @@ 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]) + # Now split the region into subregions + 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] @@ -456,11 +488,12 @@ 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 + if nAmpsNan == len(ampNames): msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}." self.log.warning(msg) diff --git a/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py b/python/lsst/cp/pipe/ptc/cpSolvePtcTask.py index 559340a8e..1cfb17837 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 @@ -223,27 +226,27 @@ 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].append(partialPtcDataset.inputExpIdPairs[ampName]) + datasetPtc.inputExpIdPairs[ampName] += 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]) + 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])) + 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: @@ -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):