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

Jump to bottom

DM-42163: eo_pipe/cp_verify parity: ptc #228

Merged
merged 5 commits into from
Jan 25, 2024
Merged

DM-42163: eo_pipe/cp_verify parity: ptc #228

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9d2f93e
Add row means variance calculation.
Jan 17, 2024
84b8756
Add more modifications related to rowMeanVariance
Jan 17, 2024
6e9828b
Add a few tests
Jan 17, 2024
a50f319
Modify cpSolvePtc
Jan 17, 2024
c24e7b7
Fix check for np.nan
Jan 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
43 changes: 34 additions & 9 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 @@ -27,6 +27,7 @@
import lsst.afw.math as afwMath
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
from lsst.geom import (Box2I, Point2I, Extent2I)
from lsst.cp.pipe.utils import (arrangeFlatsByExpTime, arrangeFlatsByExpId,
arrangeFlatsByExpFlux, sigmaClipCorrection,
CovFastFourierTransform)
Expand Down Expand Up @@ -558,7 +559,11 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
# 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)
muDiff, varDiff, covAstier, rowMeanVariance = 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,
Expand All @@ -579,8 +584,10 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
# 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 "
if (np.isnan(muDiff) or np.isnan(varDiff) or (covAstier is None)
or (rowMeanVariance is np.nan)):
self.log.warning("NaN mean, var or rowmeanVariance, or None cov in amp %s "
Copy link
Contributor

Choose a reason for hiding this comment

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

"rowMeanVariance"? Just for consistency.

"in exposure pair %d, %d of "
"detector %d.", ampName, expId1, expId2, detNum)
nAmpsNan += 1
expIdMask = False
Expand Down Expand Up @@ -664,6 +671,7 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
histVar=histVar,
histChi2Dof=histChi2Dof,
kspValue=kspValue,
rowMeanVariance=rowMeanVariance,
)

partialPtcDataset.setAuxValuesPartialDataset(auxDict)
Expand Down Expand Up @@ -789,11 +797,11 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):

Returns
-------
mu : `float` or `NaN`
mu : `float`
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`
varDiff : `float`
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.
Expand All @@ -809,12 +817,15 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):
Covariance at (dx, dy).
nPix : `int`
Number of pixel pairs used to evaluate var and cov.
rowMeanVariance : `float`
Variance of the means of each row in the difference image.
Taken from `github.com/lsst-camera-dh/eo_pipe`.

If either mu1 or m2 are NaN's, the returned value is NaN.
"""
if np.isnan(mu1) or np.isnan(mu2):
self.log.warning("Mean of amp in image 1 or 2 is NaN: %f, %f.", mu1, mu2)
return np.nan, np.nan, None
return np.nan, np.nan, None, np.nan
Copy link
Contributor

Choose a reason for hiding this comment

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

The test above (L587) tests for is None, not np.nan. Same point for the other returns below.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Changed the test to check for np.nan

mu = 0.5*(mu1 + mu2)

# Take difference of pairs
Expand All @@ -829,6 +840,20 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):
if self.config.binSize > 1:
diffIm = afwMath.binImage(diffIm, self.config.binSize)

# Calculate the variance (in ADU^2) of the means of rows for diffIm.
# Taken from eo_pipe
region = diffIm.getBBox()
rowMeans = []
for row in range(region.minY, region.maxY):
regionRow = Box2I(Point2I(region.minX, row),
Extent2I(region.width, 1))
rowMeans.append(afwMath.makeStatistics(diffIm[regionRow],
afwMath.MEAN,
imStatsCtrl).getValue())
rowMeanVariance = afwMath.makeStatistics(
np.array(rowMeans), afwMath.VARIANCECLIP,
imStatsCtrl).getValue()

# Variance calculation via afwMath
varDiff = 0.5*(afwMath.makeStatistics(diffIm, afwMath.VARIANCECLIP, imStatsCtrl).getValue())

Expand All @@ -851,7 +876,7 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):
self.log.warning("Number of good points for covariance calculation (%s) is less "
"(than threshold %s)", np.sum(w),
self.config.minNumberGoodPixelsForCovariance/(self.config.binSize**2))
return np.nan, np.nan, None
return np.nan, np.nan, None, np.nan

maxRangeCov = self.config.maximumRangeCovariancesAstier

Expand All @@ -870,7 +895,7 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):
# This is raised if there are not enough pixels.
self.log.warning("Not enough pixels covering the requested covariance range in x/y (%d)",
self.config.maximumRangeCovariancesAstier)
return np.nan, np.nan, None
return np.nan, np.nan, None, np.nan

# Compare Cov[0,0] and afwMath.VARIANCECLIP covDiffAstier[0]
# is the Cov[0,0] element, [3] is the variance, and there's a
Expand All @@ -881,7 +906,7 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):
self.log.warning("Absolute fractional difference between afwMatch.VARIANCECLIP and Cov[0,0] "
"is more than %f%%: %f", thresholdPercentage, fractionalDiff)

return mu, varDiff, covDiffAstier
return mu, varDiff, covDiffAstier, rowMeanVariance

def getImageAreasMasksStats(self, exposure1, exposure2, region=None):
"""Get image areas in a region as well as masks and statistic objects.
Expand Down
3 changes: 3 additions & 0 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 @@ -350,6 +350,8 @@ def run(self, inputCovariances, camera=None, detId=0):
partialPtcDataset.rawMeans[ampName][0])
datasetPtc.rawVars[ampName] = np.append(datasetPtc.rawVars[ampName],
partialPtcDataset.rawVars[ampName][0])
datasetPtc.rowMeanVariance[ampName] = np.append(datasetPtc.rowMeanVariance[ampName],
partialPtcDataset.rowMeanVariance[ampName][0])
datasetPtc.photoCharges[ampName] = np.append(datasetPtc.photoCharges[ampName],
partialPtcDataset.photoCharges[ampName][0])
datasetPtc.histVars[ampName] = np.append(datasetPtc.histVars[ampName],
Expand Down Expand Up @@ -423,6 +425,7 @@ def run(self, inputCovariances, camera=None, detId=0):
datasetPtc.rawExpTimes[ampName] = datasetPtc.rawExpTimes[ampName][index]
datasetPtc.rawMeans[ampName] = datasetPtc.rawMeans[ampName][index]
datasetPtc.rawVars[ampName] = datasetPtc.rawVars[ampName][index]
datasetPtc.rowMeanVariance[ampName] = datasetPtc.rowMeanVariance[ampName][index]
datasetPtc.photoCharges[ampName] = datasetPtc.photoCharges[ampName][index]
datasetPtc.histVars[ampName] = datasetPtc.histVars[ampName][index]
datasetPtc.histChi2Dofs[ampName] = datasetPtc.histChi2Dofs[ampName][index]
Expand Down
15 changes: 9 additions & 6 deletions tests/test_ptc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@ def test_covAstier(self):
mu1,
mu2,
) = extractTask.getImageAreasMasksStats(mockExp1, mockExp2)
muDiff, varDiff, covAstier = extractTask.measureMeanVarCov(
muDiff, varDiff, covAstier, rowMeanVariance = extractTask.measureMeanVarCov(
im1Area, im2Area, imStatsCtrl, mu1, mu2
)
muStandard.setdefault(ampName, []).append(muDiff)
Expand Down Expand Up @@ -697,7 +697,7 @@ def test_meanVarMeasurement(self):
im1Area, im2Area, imStatsCtrl, mu1, mu2 = task.getImageAreasMasksStats(
self.flatExp1, self.flatExp2
)
mu, varDiff, _ = task.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2)
mu, varDiff, _, _ = task.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2)

self.assertLess(self.flatWidth - np.sqrt(varDiff), 1)
self.assertLess(self.flatMean - mu, 1)
Expand All @@ -714,7 +714,7 @@ def test_meanVarMeasurementWithNans(self):
im1Area, im2Area, imStatsCtrl, mu1, mu2 = task.getImageAreasMasksStats(
flatExp1, flatExp2
)
mu, varDiff, _ = task.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2)
mu, varDiff, _, _ = task.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2)

expectedMu1 = np.nanmean(flatExp1.image.array)
expectedMu2 = np.nanmean(flatExp2.image.array)
Expand Down Expand Up @@ -751,13 +751,14 @@ def test_meanVarMeasurementAllNan(self):
im1Area, im2Area, imStatsCtrl, mu1, mu2 = task.getImageAreasMasksStats(
flatExp1, flatExp2
)
mu, varDiff, covDiff = task.measureMeanVarCov(
mu, varDiff, covDiff, rowMeanVariance = task.measureMeanVarCov(
im1Area, im2Area, imStatsCtrl, mu1, mu2
)

self.assertTrue(np.isnan(mu))
self.assertTrue(np.isnan(varDiff))
self.assertTrue(covDiff is None)
self.assertTrue(np.isnan(rowMeanVariance))

def test_meanVarMeasurementTooFewPixels(self):
task = self.defaultTaskExtract
Expand All @@ -775,14 +776,15 @@ def test_meanVarMeasurementTooFewPixels(self):
flatExp1, flatExp2
)
with self.assertLogs(level=logging.WARNING) as cm:
mu, varDiff, covDiff = task.measureMeanVarCov(
mu, varDiff, covDiff, rowMeanVariance = task.measureMeanVarCov(
im1Area, im2Area, imStatsCtrl, mu1, mu2
)
self.assertIn("Number of good points", cm.output[0])

self.assertTrue(np.isnan(mu))
self.assertTrue(np.isnan(varDiff))
self.assertTrue(covDiff is None)
self.assertTrue(np.isnan(rowMeanVariance))

def test_meanVarMeasurementTooNarrowStrip(self):
# We need a new config to make sure the second covariance cut is
Expand All @@ -808,14 +810,15 @@ def test_meanVarMeasurementTooNarrowStrip(self):
flatExp1, flatExp2
)
with self.assertLogs(level=logging.WARNING) as cm:
mu, varDiff, covDiff = task.measureMeanVarCov(
mu, varDiff, covDiff, rowMeanVariance = task.measureMeanVarCov(
im1Area, im2Area, imStatsCtrl, mu1, mu2
)
self.assertIn("Not enough pixels", cm.output[0])

self.assertTrue(np.isnan(mu))
self.assertTrue(np.isnan(varDiff))
self.assertTrue(covDiff is None)
self.assertTrue(np.isnan(rowMeanVariance))

def test_makeZeroSafe(self):
noZerosArray = [1.0, 20, -35, 45578.98, 90.0, 897, 659.8]
Expand Down