Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-37951: Allow PTC analysis to split each amp into subregions #171

Open
wants to merge 8 commits into
base: main
Choose a base branch
from
Open

DM-37951: Allow PTC analysis to split each amp into subregions #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
cd36f61
Initial push
craiglagegit Feb 10, 2023
5259bc8
Fixing indexing problem
craiglagegit Feb 11, 2023
c711f34
Another indexing correction
craiglagegit Feb 11, 2023
a2a0daf
Correcting misunderstanding on ptcDatasetList
craiglagegit Feb 11, 2023
3d2dade
Lots of changes to deal with the larger arrays. Still not working
craiglagegit Feb 11, 2023
0bca008
Changes to deal with larger arrays. cpExtract now runs, but still fa…
craiglagegit Feb 12, 2023
5fd25ba
More fixes and clean-up. Seems to be working now
craiglagegit Feb 12, 2023
fef1f50
Pushing changes to getInitialGoodPoints
craiglagegit Feb 12, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 105 additions & 72 deletions python/lsst/cp/pipe/ptc/cpExtractPtcTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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']

Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's probably not worth the effort now, but this could have used the cpCombine.CpCombineTask._subBBoxIter static method to chunk the image and reuse code.

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]
Expand All @@ -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)
Expand Down
72 changes: 27 additions & 45 deletions python/lsst/cp/pipe/ptc/cpSolvePtcTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't this change the dimensionality? Previously the partial inputs were lists of tuples, appended into a list-of-list-of-tuples. Using += flattens it into a list-of-tuples again. I think this less of a problem below, as it looks like it was just selecting the first element of the input lists for those entries.

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:
Expand Down Expand Up @@ -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
------
Expand All @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These should probably be config options. Was 101 chosen for a particular reason, or is just "an odd integer about 100"?

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):
Expand Down