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.
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-12445: Set appropriate default configs for CompareWarp Coadds #162
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Some changes required, but shouldn't be too hard.
"Otherwise, do not mask.", | ||
dtype=float, | ||
default=0.075 | ||
maxNumEpochs = pexConfig.Field( |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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)
?
@@ -1769,25 +1766,25 @@ def computeAltMaskList(self, tempExpRefList, maskSpanSets): | |||
|
|||
return altMaskList | |||
|
|||
def _filterArtifacts(self, spanSetList, epochCountImage, maxNumEpochs=None): | |||
def _filterArtifacts(self, spanSetList, epochCountImage): |
There was a problem hiding this comment.
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.
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)) | ||
idx = numpy.where((counts > 0) & (counts <= self.config.maxNumEpochs)) | ||
percentBelowThreshold = len(idx[0]) / len(counts) | ||
if percentBelowThreshold > self.config.spatialThreshold: | ||
maskSpanSetList.append(span) | ||
return maskSpanSetList | ||
|
||
def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd): |
There was a problem hiding this comment.
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.
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)) | ||
idx = numpy.where((counts > 0) & (counts <= self.config.maxNumEpochs)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You should avoid using numpy.where
whenever you can (it's expensive), and just use a boolean array. Then, len(idx[0])
below simply becomes array.sum()
.
# Warp comparison must use PSF-Matched Warps regardless of requested coadd warp type | ||
warpName = self.getTempExpDatasetName('psfMatched') | ||
# Warp comparison should use PSF-Matched Warps regardless of requested coadd warp type | ||
warpName = self.getTempExpDatasetName(self.config.assembleStaticSkyModel.warpType) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's not a good idea to go reaching through multiple config elements to get what you want: what if I retarget
the assembleStaticSkyModel
and name warpType
something else? You should only access what you control. That suggests putting warpType
into the config at this level, and passing it down to the assembleStaticSkyModel
as it needs it.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hmmm. How would you have done this? Passed the new config down to the subtask in validate
? I configure the defaults for assembleStaticSkyModel subtask (an AssembleCoaddTask
instance) in the CompareWarpAssembleConfig's setDefaults
, already requiring that it have a warpType. It's going to be 'psfMatched' in every scenario I can think of (but then again I lack imagination). I'm debating if its worth making a new config parameter or just leaving it hardcoded.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We could just leave it hardcoded for now: I also lack the imagination to think as to why anyone would use a different warpType
, since the algorithm for rejecting artifacts requires PSF-matched inputs.
Or put warpType
into the AssembleCoaddTask
's Config
, and then pass it to the AssembleStaticSkyModelTask
constructor when you makeSubtask
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I thought configs frozen by the time we got to making the subtask, but I'll give it try.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The config is definitely frozen within CompareWarpAssembleCoadd.__init__
(Pdb) kwargs['config'].assembleStaticSkyModel.warpType = staticSkyModelWarpType
*** FieldValidationError: ConfigurableField 'assembleStaticSkyModel' failed validation: Cannot modify a frozen Config
For more information read the Field definition at:
File pex/config/configurableField.py:271 (__deepcopy__)And the Config definition at:
File pipe/tasks/assembleCoadd.py:1439 (<module>)
To clarify: AssembleStaticSkyModelTask
is not a thing. It is a vanilla AssembleCoaddTask
used as a subtask within CompareWarpAssembleCoaddTask
.
Only solution I can think of at the moment is:
def validate(self):
AssembleCoaddConfig.validate(self)
self.assembleStaticSkyModel.warpType = self.staticSkyModelWarpType
which I'm not sure is any better because it smells like code duplication. staticSkyModelWarpType
is just a copied and pasted warpType
config from:
https://github.com/lsst/pipe_tasks/blob/master/python/lsst/pipe/tasks/assembleCoadd.py#L55
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Reverting to hardcoded 'psfMatched'. Jenkins running.
doc="Unitless fraction of pixels defining how much of the outlier region has to meet the " | ||
"temporal criteria", | ||
dtype=float, | ||
default=0.5 | ||
default=0.5, | ||
min=0., max=1. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Add inclusiveMax=True
?
self.assembleStaticSkyModel.warpType = 'psfMatched' | ||
self.assembleStaticSkyModel.statistic = 'MEDIAN' | ||
self.assembleStaticSkyModel.statistic = 'MEANCLIP' | ||
self.assembleStaticSkyModel.sigmaClip = 1.5 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why such harsh clipping?
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.)
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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?
9b0b32e
to
d0a5edd
Compare
Replace config parameter that controls maximum number of epochs an artifact candidate can appear in without being masking, from a fraction of N to an integer.
d0a5edd
to
db3de1f
Compare
No description provided.