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-12445: Set appropriate default configs for CompareWarp Coadds #162

Merged
merged 3 commits into from
Nov 18, 2017
Merged

DM-12445: Set appropriate default configs for CompareWarp Coadds #162

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7e191fa
Replace temporalThreshold with maxNumEpochs
yalsayyad Nov 18, 2017
eea84d2
Add a docstring and remove a duplicate function call
yalsayyad Nov 18, 2017
db3de1f
Set CompareWarp configuration defaults
yalsayyad Nov 13, 2017
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
60 changes: 40 additions & 20 deletions python/lsst/pipe/tasks/assembleCoadd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1446,28 +1446,31 @@ class CompareWarpAssembleCoaddConfig(AssembleCoaddConfig):
target=SourceDetectionTask,
doc="Detect outlier sources on difference between each psfMatched warp and static sky model"
)
temporalThreshold = pexConfig.Field(
doc="Unitless fraction of number of epochs to classify as an artifact/outlier source versus"
" a source that is intrinsically variable or difficult to subtract cleanly. "
"If outlier region in warp-diff Chi-image is mostly (defined by spatialThreshold) "
"an outlier in less than temporalThreshold * number of epochs, then mask. "
"Otherwise, do not mask.",
dtype=float,
default=0.075
maxNumEpochs = pexConfig.Field(
Copy link
Contributor

Choose a reason for hiding this comment

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

I would have thought that a fraction would be so much more generally useful than a hard number. I'm worried that a hard number isn't going to scale as you go from n=4 to n=400 inputs.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Just made a new ticket for this: DM-12692. I had talked briefly with @RobertLuptonTheGood about it. You're right maxNumEpochs won't scale, but neither will a fraction of the maximum number of epochs. I've seen some steep Nepoch gradients in a patch. Its also depends on how many repeat pointings line up the artifacts.

Copy link
Contributor

Choose a reason for hiding this comment

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

Yeah, it probably needs some combination of a number and a fraction, maybe max(maxNumEpochs, fracMaxEpochs*num)?

doc="Maximum number of epochs/visits in which an artifact candidate can appear and still be masked. "
"For each footprint detected on the image difference between the psfMatched warp and static sky "
"model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more "
"than maxNumEpochs, the artifact candidate is persistant rather than transient and not masked.",
dtype=int,
default=2
)
spatialThreshold = pexConfig.Field(
spatialThreshold = pexConfig.RangeField(
doc="Unitless fraction of pixels defining how much of the outlier region has to meet the "
"temporal criteria",
"temporal criteria. If 0, clip all. If 1, clip none.",
dtype=float,
default=0.5
default=0.5,
min=0., max=1.,
inclusiveMin=True, inclusiveMax=True
)

def setDefaults(self):
AssembleCoaddConfig.setDefaults(self)
self.statistic = 'MEAN'
self.assembleStaticSkyModel.warpType = 'psfMatched'
self.assembleStaticSkyModel.statistic = 'MEDIAN'
self.assembleStaticSkyModel.statistic = 'MEANCLIP'
self.assembleStaticSkyModel.sigmaClip = 1.5
Copy link
Contributor

Choose a reason for hiding this comment

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

Why such harsh clipping?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Borrowed from SafeClip. Running a grid of configs is still a to do: DM-12698

Copy link
Contributor

Choose a reason for hiding this comment

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

But harsh clipping like this is part of the SafeClip algorithm (we deliberately toss out anything that looks even remotely bad in order to make a very clean model of the sky), whereas I wouldn't have thought it would be for this.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This is the same idea. Make a very very clean and naive robust PSF-match coadd to serve as the naive model of the static sky.

Copy link
Contributor

Choose a reason for hiding this comment

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

Ah, I see --- this is for the assembleStaticSkyModel. Fair enough, then.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

For 2< N<=5 especially, we have to be strict to get the junk out of the naive template:

import  lsst.afw.math as afwMath
arr = [1,10,10,20]
sctrl = afwMath.StatisticsControl()
sctrl.setNumSigmaClip(1.5)
sctrl.setNumIter(3)
print afwMath.makeStatistics(arr, afwMath.MEANCLIP, sctrl).getValue()
print afwMath.makeStatistics(arr, afwMath.STDEV, sctrl).getValue()
> 10.0
> 7.76208734813

Copy link
Contributor

Choose a reason for hiding this comment

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

Perhaps that's due to a shortcoming in our statistics suite (no robust stdev), or because you're not using STDEVCLIP:

>>> ss = afwMath.makeStatistics(arr, afwMath.MEANCLIP | afwMath.STDEVCLIP | afwMath.NCLIPPED, sctrl)
>>> ss.getValue(afwMath.MEANCLIP)
10.0
>>> ss.getValue(afwMath.STDEVCLIP)
0.0
>>> ss.getValue(afwMath.NCLIPPED)
2.0

(NCLIPPED is a feature being implemented on DM-9953.)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That's cool. My point was that if I increase sctrl.setNumSigmaClip(1.5) to sctrl.setNumSigmaClip(2.6) the 2 outliers don't get clipped anymore.

>>> sctrl.setNumSigmaClip(2.6)
>>> print afwMath.makeStatistics(arr, afwMath.MEANCLIP, sctrl).getValue()
10.25

Copy link
Contributor

Choose a reason for hiding this comment

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

OK. How common is it that 40% of the inputs have outliers?

self.assembleStaticSkyModel.clipIter = 3
self.assembleStaticSkyModel.doWrite = False
self.statistic = 'MEAN'
self.detect.doTempLocalBackground = False
self.detect.reEstimateBackground = False
self.detect.returnOriginalFootprints = False
Expand Down Expand Up @@ -1721,7 +1724,6 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList):
spanSetArtifactList = []
spanSetNoDataMaskList = []

maxNumEpochs = int(max(1, self.config.temporalThreshold*len(tempExpRefList)))
for warpRef, imageScaler in zip(tempExpRefList, imageScalerList):
warpDiffExp = self._readAndComputeWarpDiff(warpRef, imageScaler, templateCoadd)
if warpDiffExp is not None:
Expand Down Expand Up @@ -1758,8 +1760,7 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList):

for i, spanSetList in enumerate(spanSetArtifactList):
if spanSetList:
filteredSpanSetList = self._filterArtifacts(spanSetList, epochCountImage,
maxNumEpochs=maxNumEpochs)
filteredSpanSetList = self._filterArtifacts(spanSetList, epochCountImage)
spanSetArtifactList[i] = filteredSpanSetList

return pipeBase.Struct(artifacts=spanSetArtifactList,
Expand Down Expand Up @@ -1794,25 +1795,44 @@ def computeAltMaskList(self, tempExpRefList, maskSpanSets):

return altMaskList

def _filterArtifacts(self, spanSetList, epochCountImage, maxNumEpochs=None):
def _filterArtifacts(self, spanSetList, epochCountImage):
Copy link
Contributor

Choose a reason for hiding this comment

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

This method needs a docstring.

"""!
\brief Filter artifact candidates

@param spanSetList: List of SpanSets representing artifact candidates
@param epochCountImage: Image of accumulated number of warpDiff detections

return List of SpanSets with artifacts
"""

maskSpanSetList = []
x0, y0 = epochCountImage.getXY0()
for i, span in enumerate(spanSetList):
y, x = span.indices()
counts = epochCountImage.array[[y1 - y0 for y1 in y], [x1 - x0 for x1 in x]]
idx = numpy.where((counts > 0) & (counts <= maxNumEpochs))
percentBelowThreshold = len(idx[0]) / len(counts)
nCountsBelowThreshold = numpy.count_nonzero((counts > 0) & (counts <= self.config.maxNumEpochs))
percentBelowThreshold = nCountsBelowThreshold / len(counts)
if percentBelowThreshold > self.config.spatialThreshold:
maskSpanSetList.append(span)
return maskSpanSetList

def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd):
Copy link
Contributor

Choose a reason for hiding this comment

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

This method needs a docstring.

"""!
\brief Fetch a warp from the butler and return a warpDiff

@param warpRef: `Butler dataRef` for the warp
@param imageScaler: `scaleZeroPoint.ImageScaler` object
@param templateCoadd: Exposure to be substracted from the scaled warp

return Exposure of the image difference between the warp and template
"""

# Warp comparison must use PSF-Matched Warps regardless of requested coadd warp type
warpName = self.getTempExpDatasetName('psfMatched')
if not warpRef.datasetExists(warpName):
self.log.warn("Could not find %s %s; skipping it", warpName, warpRef.dataId)
return None
warp = warpRef.get(self.getTempExpDatasetName('psfMatched'), immediate=True)
warp = warpRef.get(warpName, immediate=True)
# direct image scaler OK for PSF-matched Warp
imageScaler.scaleMaskedImage(warp.getMaskedImage())
mi = warp.getMaskedImage()
Expand Down