Inheritance diagram for lsst.cp.pipe.ptc.cpExtractPtcTask.PhotonTransferCurveExtractTask:

Public Member Functions
def	runQuantum (self, butlerQC, inputRefs, outputRefs)

def	run (self, inputExp, inputDims)

def	measureMeanVarCov (self, exposure1, exposure2, region=None, covAstierRealSpace=False)

Static Public Attributes
	ConfigClass = PhotonTransferCurveExtractConfig

Detailed Description

Task to measure covariances from flat fields.
This task receives as input a list of flat-field images
(flats), and sorts these flats in pairs taken at the
same time (if there's a different number of flats,
those flats are discarded). The mean, variance, and
covariances are measured from the difference of the flat
pairs at a given time. The variance is calculated
via afwMath, and the covariance via the methods in Astier+19
(appendix A). In theory, var = covariance[0,0]. This should
be validated, and in the future, we may decide to just keep
one (covariance).

The measured covariances at a particular time (along with
other quantities such as the mean) are stored in a PTC dataset
object (`PhotonTransferCurveDataset`), which gets partially
filled. The number of partially-filled PTC dataset objects
will be less than the number of input exposures, but gen3
requires/assumes that the number of input dimensions matches
bijectively the number of output dimensions. Therefore, a
number of "dummy" PTC dataset are inserted in the output list
that has the partially-filled PTC datasets with the covariances.
This output list will be used as input of
`PhotonTransferCurveSolveTask`, which will assemble the multiple
`PhotonTransferCurveDataset`s into a single one in order to fit
the measured covariances as a function of flux to a particular
model.

Astier+19: "The Shape of the Photon Transfer Curve of CCD
sensors", arXiv:1905.08677.

Definition at line 156 of file cpExtractPtcTask.py.

Member Function Documentation

◆ measureMeanVarCov()

def lsst.cp.pipe.ptc.cpExtractPtcTask.PhotonTransferCurveExtractTask.measureMeanVarCov	(	self,
		exposure1,
		exposure2,
		region = `None`,
		covAstierRealSpace = `False`
	)

Calculate the mean of each of two exposures and the variance
and covariance of their difference. The variance is calculated
via afwMath, and the covariance via the methods in Astier+19
(appendix A). In theory, var = covariance[0,0]. This should
be validated, and in the future, we may decide to just keep
one (covariance).

Parameters
----------
exposure1 : `lsst.afw.image.exposure.exposure.ExposureF`
    First exposure of flat field pair.
exposure2 : `lsst.afw.image.exposure.exposure.ExposureF`
    Second exposure of flat field pair.
region : `lsst.geom.Box2I`, optional
    Region of each exposure where to perform the calculations (e.g, an amplifier).
covAstierRealSpace : `bool`, optional
    Should the covariannces in Astier+19 be calculated in real space or via FFT?
    See Appendix A of Astier+19.

Returns
-------
mu : `float` or `NaN`
    0.5*(mu1 + mu2), where mu1, and mu2 are the clipped means of the regions in
    both exposures. If either mu1 or m2 are NaN's, the returned value is NaN.
varDiff : `float` or `NaN`
    Half of the clipped variance of the difference of the regions inthe two input
    exposures. If either mu1 or m2 are NaN's, the returned value is NaN.
covDiffAstier : `list` or `NaN`
    List with tuples of the form (dx, dy, var, cov, npix), where:
        dx : `int`
            Lag in x
        dy : `int`
            Lag in y
        var : `float`
            Variance at (dx, dy).
        cov : `float`
            Covariance at (dx, dy).
        nPix : `int`
            Number of pixel pairs used to evaluate var and cov.
    If either mu1 or m2 are NaN's, the returned value is NaN.

Definition at line 361 of file cpExtractPtcTask.py.

     def measureMeanVarCov(self, exposure1, exposure2, region=None, covAstierRealSpace=False):
         """Calculate the mean of each of two exposures and the variance
         and covariance of their difference. The variance is calculated
         via afwMath, and the covariance via the methods in Astier+19
         (appendix A). In theory, var = covariance[0,0]. This should
         be validated, and in the future, we may decide to just keep
         one (covariance).
  
         Parameters
         ----------
         exposure1 : `lsst.afw.image.exposure.exposure.ExposureF`
             First exposure of flat field pair.
         exposure2 : `lsst.afw.image.exposure.exposure.ExposureF`
             Second exposure of flat field pair.
         region : `lsst.geom.Box2I`, optional
             Region of each exposure where to perform the calculations (e.g, an amplifier).
         covAstierRealSpace : `bool`, optional
             Should the covariannces in Astier+19 be calculated in real space or via FFT?
             See Appendix A of Astier+19.
  
         Returns
         -------
         mu : `float` or `NaN`
             0.5*(mu1 + mu2), where mu1, and mu2 are the clipped means of the regions in
             both exposures. If either mu1 or m2 are NaN's, the returned value is NaN.
         varDiff : `float` or `NaN`
             Half of the clipped variance of the difference of the regions inthe two input
             exposures. If either mu1 or m2 are NaN's, the returned value is NaN.
         covDiffAstier : `list` or `NaN`
             List with tuples of the form (dx, dy, var, cov, npix), where:
                 dx : `int`
                     Lag in x
                 dy : `int`
                     Lag in y
                 var : `float`
                     Variance at (dx, dy).
                 cov : `float`
                     Covariance at (dx, dy).
                 nPix : `int`
                     Number of pixel pairs used to evaluate var and cov.
             If either mu1 or m2 are NaN's, the returned value is NaN.
         """
  
         if region is not None:
             im1Area = exposure1.maskedImage[region]
             im2Area = exposure2.maskedImage[region]
         else:
             im1Area = exposure1.maskedImage
             im2Area = exposure2.maskedImage
  
         if self.config.binSize > 1:
             im1Area = afwMath.binImage(im1Area, self.config.binSize)
             im2Area = afwMath.binImage(im2Area, self.config.binSize)
  
         im1MaskVal = exposure1.getMask().getPlaneBitMask(self.config.maskNameList)
         im1StatsCtrl = afwMath.StatisticsControl(self.config.nSigmaClipPtc,
                                                  self.config.nIterSigmaClipPtc,
                                                  im1MaskVal)
         im1StatsCtrl.setNanSafe(True)
         im1StatsCtrl.setAndMask(im1MaskVal)
  
         im2MaskVal = exposure2.getMask().getPlaneBitMask(self.config.maskNameList)
         im2StatsCtrl = afwMath.StatisticsControl(self.config.nSigmaClipPtc,
                                                  self.config.nIterSigmaClipPtc,
                                                  im2MaskVal)
         im2StatsCtrl.setNanSafe(True)
         im2StatsCtrl.setAndMask(im2MaskVal)
  
         #  Clipped mean of images; then average of mean.
         mu1 = afwMath.makeStatistics(im1Area, afwMath.MEANCLIP, im1StatsCtrl).getValue()
         mu2 = afwMath.makeStatistics(im2Area, afwMath.MEANCLIP, im2StatsCtrl).getValue()
         if np.isnan(mu1) or np.isnan(mu2):
             self.log.warn(f"Mean of amp in image 1 or 2 is NaN: {mu1}, {mu2}.")
             return np.nan, np.nan, None
         mu = 0.5*(mu1 + mu2)
  
         # Take difference of pairs
         # symmetric formula: diff = (mu2*im1-mu1*im2)/(0.5*(mu1+mu2))
         temp = im2Area.clone()
         temp *= mu1
         diffIm = im1Area.clone()
         diffIm *= mu2
         diffIm -= temp
         diffIm /= mu
  
         diffImMaskVal = diffIm.getMask().getPlaneBitMask(self.config.maskNameList)
         diffImStatsCtrl = afwMath.StatisticsControl(self.config.nSigmaClipPtc,
                                                     self.config.nIterSigmaClipPtc,
                                                     diffImMaskVal)
         diffImStatsCtrl.setNanSafe(True)
         diffImStatsCtrl.setAndMask(diffImMaskVal)
  
         # Variance calculation via afwMath
         varDiff = 0.5*(afwMath.makeStatistics(diffIm, afwMath.VARIANCECLIP, diffImStatsCtrl).getValue())
  
         # Covariances calculations
         # Get the pixels that were not clipped
         varClip = afwMath.makeStatistics(diffIm, afwMath.VARIANCECLIP, diffImStatsCtrl).getValue()
         meanClip = afwMath.makeStatistics(diffIm, afwMath.MEANCLIP, diffImStatsCtrl).getValue()
         cut = meanClip + self.config.nSigmaClipPtc*np.sqrt(varClip)
         unmasked = np.where(np.fabs(diffIm.image.array) <= cut, 1, 0)
  
         # Get the pixels in the mask planes of teh differenc eimage that were ignored
         # by the clipping algorithm
         wDiff = np.where(diffIm.getMask().getArray() == 0, 1, 0)
         # Combine the two sets of pixels ('1': use; '0': don't use) into a final weight matrix
         # to be used in the covariance calculations below.
         w = unmasked*wDiff
  
         if np.sum(w) < self.config.minNumberGoodPixelsForCovariance:
             self.log.warn(f"Number of good points for covariance calculation ({np.sum(w)}) is less "
                           f"(than threshold {self.config.minNumberGoodPixelsForCovariance})")
             return np.nan, np.nan, None
  
         maxRangeCov = self.config.maximumRangeCovariancesAstier
         if covAstierRealSpace:
             # Calculate  covariances in real space.
             covDiffAstier = computeCovDirect(diffIm.image.array, w, maxRangeCov)
         else:
             # Calculate covariances via FFT (default).
             shapeDiff = np.array(diffIm.image.array.shape)
             # Calculate the sizes of FFT dimensions.
             s = shapeDiff + maxRangeCov
             tempSize = np.array(np.log(s)/np.log(2.)).astype(int)
             fftSize = np.array(2**(tempSize+1)).astype(int)
             fftShape = (fftSize[0], fftSize[1])
  
             c = CovFastFourierTransform(diffIm.image.array, w, fftShape, maxRangeCov)
             covDiffAstier = c.reportCovFastFourierTransform(maxRangeCov)
  
         # Compare Cov[0,0] and afwMath.VARIANCECLIP
         # covDiffAstier[0] is the Cov[0,0] element, [3] is the variance, and there's a factor of 0.5
         # difference with afwMath.VARIANCECLIP.
         thresholdPercentage = self.config.thresholdDiffAfwVarVsCov00
         fractionalDiff = 100*np.fabs(1 - varDiff/(covDiffAstier[0][3]*0.5))
         if fractionalDiff >= thresholdPercentage:
             self.log.warn("Absolute fractional difference between afwMatch.VARIANCECLIP and Cov[0,0] "
                           f"is more than {thresholdPercentage}%: {fractionalDiff}")
  
         return mu, varDiff, covDiffAstier

◆ run()

def lsst.cp.pipe.ptc.cpExtractPtcTask.PhotonTransferCurveExtractTask.run	(	self,
		inputExp,
		inputDims
	)

Measure covariances from difference of flat pairs

Parameters
----------
inputExp : `dict` [`float`,
                (`~lsst.afw.image.exposure.exposure.ExposureF`,
                `~lsst.afw.image.exposure.exposure.ExposureF`, ...,
                `~lsst.afw.image.exposure.exposure.ExposureF`)]
    Dictionary that groups flat-field exposures that have the same
    exposure time (seconds).

inputDims : `list`
    List of exposure IDs.

Definition at line 216 of file cpExtractPtcTask.py.

     def run(self, inputExp, inputDims):
         """Measure covariances from difference of flat pairs
  
         Parameters
         ----------
         inputExp : `dict` [`float`,
                         (`~lsst.afw.image.exposure.exposure.ExposureF`,
                         `~lsst.afw.image.exposure.exposure.ExposureF`, ...,
                         `~lsst.afw.image.exposure.exposure.ExposureF`)]
             Dictionary that groups flat-field exposures that have the same
             exposure time (seconds).
  
         inputDims : `list`
             List of exposure IDs.
         """
         # inputExp.values() returns a view, which we turn into a list. We then
         # access the first exposure-ID tuple to get the detector.
         detector = list(inputExp.values())[0][0][0].getDetector()
         detNum = detector.getId()
         amps = detector.getAmplifiers()
         ampNames = [amp.getName() for amp in amps]
  
         # Each amp may have a different  min and max ADU signal specified in the config.
         maxMeanSignalDict = {ampName: 1e6 for ampName in ampNames}
         minMeanSignalDict = {ampName: 0.0 for ampName in ampNames}
         for ampName in ampNames:
             if 'ALL_AMPS' in self.config.maxMeanSignal:
                 maxMeanSignalDict[ampName] = self.config.maxMeanSignal['ALL_AMPS']
             elif ampName in self.config.maxMeanSignal:
                 maxMeanSignalDict[ampName] = self.config.maxMeanSignal[ampName]
  
             if 'ALL_AMPS' in self.config.minMeanSignal:
                 minMeanSignalDict[ampName] = self.config.minMeanSignal['ALL_AMPS']
             elif ampName in self.config.minMeanSignal:
                 minMeanSignalDict[ampName] = self.config.minMeanSignal[ampName]
         # These are the column names for `tupleRows` below.
         tags = [('mu', '<f8'), ('afwVar', '<f8'), ('i', '<i8'), ('j', '<i8'), ('var', '<f8'),
                 ('cov', '<f8'), ('npix', '<i8'), ('ext', '<i8'), ('expTime', '<f8'), ('ampName', '<U3')]
         # Create a dummy ptcDataset
         dummyPtcDataset = PhotonTransferCurveDataset(ampNames, 'DUMMY',
                                                      self.config.maximumRangeCovariancesAstier)
         # Initialize amps of `dummyPtcDatset`.
         for ampName in ampNames:
             dummyPtcDataset.setAmpValues(ampName)
         # Output list with PTC datasets.
         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)):
             partialPtcDatasetList.append(dummyPtcDataset)
  
         if self.config.numEdgeSuspect > 0:
             isrTask = IsrTask()
             self.log.info(f"Masking {self.config.numEdgeSuspect} pixels from the edges "
                           "of all exposures as SUSPECT.")
  
         for expTime in inputExp:
             exposures = inputExp[expTime]
             if len(exposures) == 1:
                 self.log.warn(f"Only one exposure found at expTime {expTime}. Dropping exposure "
                               f"{exposures[0][1]}")
                 continue
             else:
                 # Only use the first two exposures at expTime. Each elements is a tuple (exposure, expId)
                 exp1, expId1 = exposures[0]
                 exp2, expId2 = exposures[1]
                 if len(exposures) > 2:
                     self.log.warn(f"Already found 2 exposures at expTime {expTime}. "
                                   "Ignoring exposures: "
                                   f"{i[1] for i in exposures[2:]}")
             # Mask pixels at the edge of the detector or of each amp
             if self.config.numEdgeSuspect > 0:
                 isrTask.maskEdges(exp1, numEdgePixels=self.config.numEdgeSuspect,
                                   maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
                 isrTask.maskEdges(exp2, numEdgePixels=self.config.numEdgeSuspect,
                                   maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
  
             nAmpsNan = 0
             partialPtcDataset = PhotonTransferCurveDataset(ampNames, '',
                                                            self.config.maximumRangeCovariancesAstier)
             for ampNumber, amp in enumerate(detector):
                 ampName = amp.getName()
                 # covAstier: [(i, j, var (cov[0,0]), cov, npix) for (i,j) in {maxLag, maxLag}^2]
                 doRealSpace = self.config.covAstierRealSpace
                 if self.config.detectorMeasurementRegion == 'AMP':
                     region = amp.getBBox()
                 elif self.config.detectorMeasurementRegion == 'FULL':
                     region = None
                 # `measureMeanVarCov` is the function that measures the variance and  covariances from
                 # the difference image of two flats at the same exposure time.
                 # The variable `covAstier` is of the form: [(i, j, var (cov[0,0]), cov, npix) for (i,j)
                 # in {maxLag, maxLag}^2]
                 muDiff, varDiff, covAstier = self.measureMeanVarCov(exp1, exp2, region=region,
                                                                     covAstierRealSpace=doRealSpace)
                 # Correction factor for 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
                 if np.isnan(muDiff) or np.isnan(varDiff) or (covAstier is None):
                     msg = (f"NaN mean or var, or None cov in amp {ampName} in exposure pair {expId1},"
                            f" {expId2} of detector {detNum}.")
                     self.log.warn(msg)
                     nAmpsNan += 1
                     expIdMask = False
                     covArray = np.full((1, self.config.maximumRangeCovariancesAstier,
                                         self.config.maximumRangeCovariancesAstier), np.nan)
                     covSqrtWeights = np.full_like(covArray, np.nan)
  
                 if (muDiff <= minMeanSignalDict[ampName]) or (muDiff >= maxMeanSignalDict[ampName]):
                     expIdMask = False
  
                 if covAstier is not None:
                     tupleRows = [(muDiff, varDiff) + covRow + (ampNumber, expTime,
                                                                ampName) for covRow in covAstier]
                     tempStructArray = np.array(tupleRows, dtype=tags)
                     covArray, vcov, _ = 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]
                 covArray[0, 0] /= varFactor
  
                 partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff],
                                                rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)],
                                                expIdMask=[expIdMask], covArray=covArray,
                                                covSqrtWeights=covSqrtWeights)
             # 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[datasetIndex] = partialPtcDataset
  
             if nAmpsNan == len(ampNames):
                 msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}."
                 self.log.warn(msg)
         return pipeBase.Struct(
             outputCovariances=partialPtcDatasetList,
         )
  

◆ runQuantum()

def lsst.cp.pipe.ptc.cpExtractPtcTask.PhotonTransferCurveExtractTask.runQuantum	(	self,
		butlerQC,
		inputRefs,
		outputRefs
	)

Ensure that the input and output dimensions are passed along.

Parameters
----------
butlerQC : `~lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
    Butler to operate on.
inputRefs : `~lsst.pipe.base.connections.InputQuantizedConnection`
    Input data refs to load.
ouptutRefs : `~lsst.pipe.base.connections.OutputQuantizedConnection`
    Output data refs to persist.

Definition at line 191 of file cpExtractPtcTask.py.

     def runQuantum(self, butlerQC, inputRefs, outputRefs):
         """Ensure that the input and output dimensions are passed along.
  
         Parameters
         ----------
         butlerQC : `~lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
             Butler to operate on.
         inputRefs : `~lsst.pipe.base.connections.InputQuantizedConnection`
             Input data refs to load.
         ouptutRefs : `~lsst.pipe.base.connections.OutputQuantizedConnection`
             Output data refs to persist.
         """
         inputs = butlerQC.get(inputRefs)
         # Ids of input list of exposures
         inputs['inputDims'] = [expId.dataId['exposure'] for expId in inputRefs.inputExp]
  
         # Dictionary, keyed by expTime, with tuples containing flat exposures and their IDs.
         if self.config.matchByExposureId:
             inputs['inputExp'] = arrangeFlatsByExpId(inputs['inputExp'], inputs['inputDims'])
         else:
             inputs['inputExp'] = arrangeFlatsByExpTime(inputs['inputExp'], inputs['inputDims'])
  
         outputs = self.run(**inputs)
         butlerQC.put(outputs, outputRefs)
  

Member Data Documentation

◆ ConfigClass

lsst.cp.pipe.ptc.cpExtractPtcTask.PhotonTransferCurveExtractTask.ConfigClass = PhotonTransferCurveExtractConfig

static

Definition at line 188 of file cpExtractPtcTask.py.

The documentation for this class was generated from the following file:

/j/snowflake/release/lsstsw/stack/lsst-scipipe-0.6.0/Linux64/cp_pipe/22.0.1-18-gb17765a+2264247a6b/python/lsst/cp/pipe/ptc/cpExtractPtcTask.py

Public Member Functions