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.

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-9615: Add dcrAssembleCoadd task #173

Merged
merged 109 commits into from Jul 20, 2018
Merged

DM-9615: Add dcrAssembleCoadd task #173

merged 109 commits into from Jul 20, 2018

Conversation

isullivan
Copy link
Contributor

This pull request is being created early to allow comments on the changes. It is not expected to be merged for a while.

__all__ = ["DcrAssembleCoaddTask"]


class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig):
Copy link
Contributor

@yalsayyad yalsayyad Dec 22, 2017

Choose a reason for hiding this comment

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

What extra methods/configs do you get from inheriting from CompareWarp vs vanilla Assemble? Oh I see, the templateCoadd.

I'm imagining you do want asteroids/CRs/ghosts clipped... but this skeleton overrides all the artifact-finding functionality of CompareWarp?

doc="Maximum number of iterations of forward modeling.",
default=0.001,
)
assembleStaticSkyModel = pexConfig.ConfigurableField(
Copy link
Contributor

Choose a reason for hiding this comment

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

Keep something like this if you inherit from Assemble, but you don't need it if you inherit from CompareWarp

" naive/first-iteration model of the static sky.",
)

def setDefaults(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

I assume this is here because you're expecting to set some Assemble defaults by the time its time to merge.



class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig):
dcrBufferSize = 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.

Probably don't need the word dcr in these configs.

CompareWarpAssembleCoaddTask.__init__(self, *args, **kwargs)

@pipeBase.timeMethod
def run(self, dataRef, selectDataList=[]):
Copy link
Contributor

Choose a reason for hiding this comment

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

I'm optimistic that we can change AssembleCoaddTask.run so that you don't need to override the entirety of run here.

self.log.info("Persisting dcrCoadd")
dataRef.put(retStruct.coaddExposure, "dcrCoadd", subfilter=retStruct.subFilter)

if coaddExp is None:
Copy link
Contributor

Choose a reason for hiding this comment

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

Oh, does this mean that this task will generate the the regular coadds too?

else:
mimage = coaddExp.getMaskedImage()
mimage += retStruct.coaddExposure.getMaskedImage()
if self.config.doWrite:
Copy link
Contributor

Choose a reason for hiding this comment

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

Indent more?


if self.config.doWrite:
self.log.info("Persisting dcrCoadd")
dataRef.put(retStruct.coaddExposure, "dcrCoadd", subfilter=retStruct.subFilter)
Copy link
Contributor

Choose a reason for hiding this comment

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

We could make a AssembleCoaddTask.write that you can override.

supplementaryData=supplementaryData,
)
coaddExp = None
for retStruct in retStructGen:
Copy link
Contributor

Choose a reason for hiding this comment

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

I think our ability to make a generic run() for both dcrAssemble and the other classes depends on how others feel about always iterating over the outputs of assemble. In most cases there will be one. It's fine with me. @r-owen? @TallJimbo?

Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

I'm not done, but I'm about to go to bed. Will pick up again later.

@@ -0,0 +1,26 @@
#!/Users/sullivan/LSST/code/lsstsw/miniconda/bin/python
Copy link
Contributor

Choose a reason for hiding this comment

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

Fix the path.


#
# LSST Data Management System
# Copyright 2008, 2009, 2010 LSST Corporation.
Copy link
Contributor

Choose a reason for hiding this comment

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

I think the new standard is to put the copyleft in a COPYRIGHT file and an AUTHORS file.
The copyright holder should be your institution, and the year should be updated.

@@ -471,6 +471,30 @@ def prepareInputs(self, refList):
return pipeBase.Struct(tempExpRefList=tempExpRefList, weightList=weightList,
imageScalerList=imageScalerList)

def prepareStats(self, mask=None):
"""!
\brief Prepare the statistics for coadding images.
Copy link
Contributor

Choose a reason for hiding this comment

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

You need to change all this to use numpydoc before merging.

Copy link
Contributor

Choose a reason for hiding this comment

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

This hasn't been converted to numpydoc yet?

bit = afwImage.Mask.getMaskPlane(plane)
statsCtrl.setMaskPropagationThreshold(bit, threshold)
statsFlags = afwMath.stringToStatisticsProperty(self.config.statistic)
return (statsCtrl, statsFlags)
Copy link
Contributor

Choose a reason for hiding this comment

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

Prefer a Struct to a tuple.


#
# LSST Data Management System
# Copyright 2008-2016 AURA/LSST.
Copy link
Contributor

Choose a reason for hiding this comment

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

Follow modern instructions for copyleft.

"""
retStruct = AssembleCoaddTask.run(self, dataRef, selectDataList=selectDataList)
for subfilter, coaddExposure in enumerate(retStruct.dcrCoadds):
if self.config.doWrite:
Copy link
Contributor

Choose a reason for hiding this comment

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

Put this before the for loop.

for subfilter, coaddExposure in enumerate(retStruct.dcrCoadds):
if self.config.doWrite:
self.writeDcrCoadd(dataRef, coaddExposure, subfilter)
return pipeBase.Struct(coaddExposure=retStruct.coaddExposure, nImage=retStruct.nImage)
Copy link
Contributor

Choose a reason for hiding this comment

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

You should return the DCR coadds too.

self.setBrightObjectMasks(coaddExposure, dataRef.dataId, brightObjectMasks)

self.log.info("Persisting dcrCoadd")
dataRef.put(coaddExposure, "dcrCoadd", subfilter=subfilter)
Copy link
Contributor

Choose a reason for hiding this comment

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

I think there's a problem with the data model: if you change the configuration of the reader to use a different number of sub-bands (why are they called "sub-bands" elsewhere and "subfilters" here?), then you would get the wrong data and wouldn't be aware of it. I think you either want to write these as data cubes (which I don't think we can do in afw) or include the total number of sub-bands in the data model (maybe the template should be ...%(subfilter)dof%(numSubfilters)d...).


if self.config.doMaskBrightObjects:
brightObjectMasks = self.readBrightObjectMasks(dataRef)
self.setBrightObjectMasks(coaddExposure, dataRef.dataId, brightObjectMasks)
Copy link
Contributor

Choose a reason for hiding this comment

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

The above operations (both interpolation and masking bright objects) don't belong in a method that writes data. They should be done before calling writeDcrCoadd so the results can be returned from the run method. It would be best if you didn't copy the implementations but just called a method to do them.


@param[in] dataRef: Data reference defining the patch for coaddition and the reference Warp
- [out] "dcrCoadd"
@param[in] coaddExposure: pipeBase.Struct with coaddExposure
Copy link
Contributor

Choose a reason for hiding this comment

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

Are you sure this is a Struct? I think it's actually an Exposure based on how you're treating it below.

Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

Still more to do...


Find artifacts and apply them to the warps' masks creating a list of alternative masks with a
new "CLIPPED" plane and updated "NO_DATA" plane.
Then pass these alternative masks to the base class's assemble method.
Copy link
Contributor

Choose a reason for hiding this comment

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

You've gotta update this description for DCR.

statsCtrl, statsFlags = self.prepareStats(mask=badPixelMask)

subregionSizeArr = self.config.subregionSize
subregionSize = afwGeom.Extent2I(subregionSizeArr[0], subregionSizeArr[1])
Copy link
Contributor

Choose a reason for hiding this comment

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

subregionSize = afwGeom.Extent2I(*self.config.subregionSize)

self.log.info("Initial convergence : %s", convergenceMetric)
convergenceList = [convergenceMetric]
convergenceCheck = 1.
while (convergenceCheck > self.config.convergenceThreshold) or (modelIter < self.config.minNIter):
Copy link
Contributor

Choose a reason for hiding this comment

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

Don't need parens here.

self.log.info("Initial convergence : %s", convergenceMetric)
convergenceList = [convergenceMetric]
convergenceCheck = 1.
while (convergenceCheck > self.config.convergenceThreshold) or (modelIter < self.config.minNIter):
Copy link
Contributor

Choose a reason for hiding this comment

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

You're not respecting the useConvergence configuration parameter. In fact, I don't see that used anywhere.

convergenceCheck = (convergenceList[-1] - convergenceMetric)/convergenceMetric
convergenceList.append(convergenceMetric)
except Exception as e:
self.log.warn("Error during iteration %s while computing coadd %s: %s", subBBox, e)
Copy link
Contributor

Choose a reason for hiding this comment

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

More % outputs than there are inputs, so this will fail.

What's supposed to happen after this error? It looks like you simply break out of the convergence loop and proceed on your merry way. Surely that's a mistake!

mask &= ~mask.getPlaneBitMask(maskPlane)
except Exception as e:
self.log.warn("Unable to remove mask plane %s: %s", maskPlane, e.message)
maskedImage -= templateImage
Copy link
Contributor

Choose a reason for hiding this comment

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

Much of the above looks like boilerplate copied from somewhere else. Can you factor it out?

try:
mask &= ~mask.getPlaneBitMask(maskPlane)
except Exception as e:
self.log.warn("Unable to remove mask plane %s: %s", maskPlane, e.message)
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't think you should do e.message. Simply casting e to str as the %s does should be sufficient.

except Exception as e:
self.log.warn("Unable to remove mask plane %s: %s", maskPlane, e.message)
maskedImage -= templateImage
obsWeight = self.calculateWeight(maskedImage, convergeMask)*1e3
Copy link
Contributor

Choose a reason for hiding this comment

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

Why the magic factor of 1000? If it's useful, shouldn't it be part of calculateWeight?

residual.setXY0(bboxGrow.getBegin())
newModel = self.clampModel(residual, oldModel, bboxGrow,
useNonNegative=self.config.useNonNegative,
clamp=self.config.clampModel)
Copy link
Contributor

Choose a reason for hiding this comment

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

Why do you need to pass config parameters that are already readily accessible to the method?

dcrSubModelOut.append(newModel)
self.regularizeModel(dcrSubModelOut, bboxGrow, baseMask,
nSigma=self.config.regularizeSigma,
clamp=self.config.clampFrequency)
Copy link
Contributor

Choose a reason for hiding this comment

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

Ditto on passing config parameters.

Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

Got through the main bulk of the code. Just the test left, I think.

"""! Calculate the weight of an exposure based on the goodness of fit of the matched template.

@param[in] residual: residual masked image after subtracting a DCR-matched template.
@param[in] goodMask: Mask plane to evaluate the goodness of fit of the residual over.
Copy link
Contributor

Choose a reason for hiding this comment

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

"Mask plane over which to evaluate..."

Also, it looks like it isn't a mask plane, but a bitmask.

@return weight: Goodness of fit metric of the residual image.
"""
residualVals = residual.getImage().getArray()
finiteInds = (numpy.isfinite(residualVals))
Copy link
Contributor

Choose a reason for hiding this comment

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

Outer parens not needed.


@return weight: Goodness of fit metric of the residual image.
"""
residualVals = residual.getImage().getArray()
Copy link
Contributor

Choose a reason for hiding this comment

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

residual.image.array

"""
residualVals = residual.getImage().getArray()
finiteInds = (numpy.isfinite(residualVals))
goodMaskInds = (residual.getMask().getArray() & goodMask) == goodMask
Copy link
Contributor

Choose a reason for hiding this comment

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

Why do you want == goodMask? Essentially it's saying that the pixels can only have DETECTED set, and nothing else. Why do you want to prohibit anything else? If you really do want to prohibit other mask planes from being set, why not specify a badMask value?

residualVals = residual.getImage().getArray()
finiteInds = (numpy.isfinite(residualVals))
goodMaskInds = (residual.getMask().getArray() & goodMask) == goodMask
weight = 1./numpy.std(residualVals[finiteInds*goodMaskInds])
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 write finiteInds & goodMaskInds since they're booleans, but that works.

"""
parAngle = visitInfo.getBoresightParAngle().asRadians()
cd = wcs.getCdMatrix()
cdAngle = (numpy.arctan2(-cd[0, 1], cd[0, 0]) + numpy.arctan2(cd[1, 0], cd[1, 1]))/2.
Copy link
Contributor

Choose a reason for hiding this comment

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

Can you still do this with SkyWcs? If so, that's a bug. We should be hiding the CD matrix from the user. Not sure of the Proper Way to do this.

rotAngle = afwGeom.Angle(cdAngle + parAngle)
return rotAngle

@staticmethod
Copy link
Contributor

Choose a reason for hiding this comment

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

Drop staticmethod.

bbox.grow(afwGeom.Extent2I(-bufferXSize, -bufferYSize))
bbox.shift(afwGeom.Extent2I(dx0 + x0, dy0 + y0))
retMask[bbox, afwImage.PARENT] |= mask[bboxBase, afwImage.PARENT]
return retMask
Copy link
Contributor

Choose a reason for hiding this comment

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

Why aren't you implementing this with our own warping functions from afw?



def wavelengthGenerator(lambdaEff, filterWidth, dcrNSubbands):
"""! tIterate over the wavelength endpoints of subfilters.
Copy link
Contributor

Choose a reason for hiding this comment

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

Typo: "tIterate"

from __future__ import absolute_import, division, print_function
#
# LSST Data Management System
# Copyright 2008-2017 AURA/LSST.
Copy link
Contributor

Choose a reason for hiding this comment

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

Update copyleft to modern standards.

Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

Hooray, finally made it all the way through!

# see <https://www.lsstcorp.org/LegalNotices/>.

from collections import namedtuple
from astropy import units as u
Copy link
Contributor

Choose a reason for hiding this comment

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

I really hate having something a module named simply u, but I understand it's standard practice...



class DcrAssembleCoaddTestTask(lsst.utils.tests.TestCase):
"""! A test case for the DCR-aware image coaddition algorithm.
Copy link
Contributor

Choose a reason for hiding this comment

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

Please use numpydoc (no leading !).

"""! A test case for the DCR-aware image coaddition algorithm.
"""

def setUp(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

There's a lot of code in here. Can you summarise in the method docstring what you're going to set up?

psfSize = 2
nSrc = 5
seed = 5
self.randGen = np.random
Copy link
Contributor

Choose a reason for hiding this comment

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

Why would you put the module into self? It's not to save on typing. If you want something in self, you could do self.randGen = np.random.RandomState(seed). I would rename randGen as rng (random number generator).

noiseLevel = 5
sourceSigma = 20.
fluxRange = 2.
edgeDist = self.config.bufferSize
Copy link
Contributor

Choose a reason for hiding this comment

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

Why edgeDist instead of bufferSize?

useInverse=False)
shift = (dcrShift[0].dy, dcrShift[0].dx)
refImage = scipy.ndimage.interpolation.shift(self.dcrModels[0].getImage().getArray(), shift)
refVariance = scipy.ndimage.interpolation.shift(self.dcrModels[0].getVariance().getArray(), shift)
Copy link
Contributor

Choose a reason for hiding this comment

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

All you're doing is rewriting the code in DcrAssembleCoaddTask.convolveDcrModelPlane --- that's not a real test. You should do a proper test: maybe centroid the sources on the model images to verify the shifts?

bboxClip.grow(afwGeom.Extent2I(-bufferXSize, -bufferYSize))

# Shifting the mask grows each mask plane by one pixel in the direction of the shift,
# so a shift followed by the reverse shift should be the same as a dilation by one pixel.
Copy link
Contributor

Choose a reason for hiding this comment

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

Only if the shift is small. You should do a large (> 1 pixel) shift as well to ensure that your code is capable of doing more than just growing a mask, but also actually shifting.

# Shifting the mask grows each mask plane by one pixel in the direction of the shift,
# so a shift followed by the reverse shift should be the same as a dilation by one pixel.
convolutionStruct = np.array([[True, True, True], [True, True, True], [True, True, True]])
maskRefCheck = dilate(model[bboxClip].getMask().getArray() == detectMask,
Copy link
Contributor

Choose a reason for hiding this comment

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

Again, why not use our own code (in this instance, lsst.afw.geom.SpanSet) instead of a third party's?

self.assertFloatsEqual(newMaskCheck, maskRefCheck)

def testRegularizationLargeClamp(self):
"""! Frequency regularization should leave the models unchanged if all the clamp factor is large.
Copy link
Contributor

Choose a reason for hiding this comment

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

"all the clamp factor"

imageDiffHigh = model.getImage().getArray() - (templateImage*clamp + noiseLevel)
self.assertLessEqual(np.max(imageDiffHigh), 0.)
imageDiffLow = model.getImage().getArray() - (templateImage/clamp - noiseLevel)
self.assertGreaterEqual(np.max(imageDiffLow), 0.)
Copy link
Contributor

Choose a reason for hiding this comment

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

This looks like a reimplementation of the code --- that's not a real test.

@isullivan isullivan force-pushed the tickets/DM-9615 branch 3 times, most recently from 5cf9cb2 to 6336f2d Compare April 6, 2018 20:59
@isullivan isullivan force-pushed the tickets/DM-9615 branch 2 times, most recently from 71c25b6 to 911ad6f Compare April 20, 2018 21:07
@isullivan isullivan force-pushed the tickets/DM-9615 branch 2 times, most recently from 2c22596 to e40ff4a Compare July 11, 2018 02:22
Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

Last major concern is I/O.

----------
coaddExposure : `lsst.afw.image.Exposure`
The coadded exposure to process.
dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
Copy link
Contributor

Choose a reason for hiding this comment

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

Please don't include the unnecessary butlerSubset.

Parameters
----------
mask : `int`, optional
Mask planes to exclude from coaddition.
Copy link
Contributor

Choose a reason for hiding this comment

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

"Mask planes" would have the type "list of str", but this is an int. Since you're plugging this into setAndMask, I think you mean the bit mask value.

Statistics structure with the following fields:

- ``statsCtrl``: Statistics control object for coadd
- ``statsFlags``: afwMath.Property object for statistic for coadd
Copy link
Contributor

Choose a reason for hiding this comment

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

I really would like to see this follow the example in the developer guide, where the type is included in parens and backticks. And the type should be fully-qualified (e.g., lsst.afw.math.Property).


Parameters
----------
maskedImage : `lsst.afw.image.maskedImage`
Copy link
Contributor

Choose a reason for hiding this comment

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

Typo: type is lsst.afw.image.MaskedImage.

------
subBBox : `lsst.afw.geom.Box2I`
Next sub-bounding box of size subregionSize or smaller; each subBBox
is contained within bbox, so it may be smaller than subregionSize at
Copy link
Contributor

Choose a reason for hiding this comment

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

Put parameter names in double-backticks, e.g.,

is contained within ``bbox``, so it may be smaller than ``subregionSize`` at


Returns
-------
`float`
Copy link
Contributor

Choose a reason for hiding this comment

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

Needs a name.


Returns
-------
`lsst.afw.image.ExposureF`
Copy link
Contributor

Choose a reason for hiding this comment

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

Needs a name.


Returns
-------
`list` of `lsst.afw.image.ExposureF`
Copy link
Contributor

Choose a reason for hiding this comment

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

Needs a name.

Best fit model of the true sky after correcting chromatic effects.
skyInfo : `lsst.pipe.base.Struct`
Patch geometry information, from getSkyInfo
tempExpRefList : `list` of `ButlerDataRef`
Copy link
Contributor

Choose a reason for hiding this comment

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

Use fully-qualified namespaces, lsst.daf.persistence.ButlerDataRef.

Best fit model of the true sky after correcting chromatic effects.
bbox : `lsst.afw.geom.box.Box2I`
Sub-region to coadd
tempExpRefList : `list` of `ButlerDataRef`
Copy link
Contributor

Choose a reason for hiding this comment

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

lsst.daf.persistence.ButlerDataRef.

@PaulPrice
Copy link
Contributor

@isullivan points out that saving the components of the DcrModel in separate FITS files makes it easy to do multiband measurements (by setting coaddName=dcr and using --id subfilter=0..2). I'm not perfectly happy about this, but it seems like a reasonable solution for the time being. In that case, the existing solution for I/O is reasonable, so I withdraw my objections.

@RobertLuptonTheGood
Copy link
Member

RobertLuptonTheGood commented Jul 19, 2018 via email

@isullivan isullivan merged commit 81851de into master Jul 20, 2018
@ktlim ktlim deleted the tickets/DM-9615 branch August 25, 2018 06:46
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants