DM-14738: Support DCR-matched templates for image differencing #93
Conversation
isullivan
force-pushed
the
tickets/DM-14738
branch
from
a35dbcc
to
f1a11b4
Compare
tests/test_dcrModel.py
Outdated
| @@ -0,0 +1,350 @@ | |||
| # This file is part of pipe_tasks. | |||
python/lsst/ip/diffim/dcrModel.py
Outdated
| @@ -0,0 +1,619 @@ | |||
| # This file is part of pipe_tasks. | |||
| Defaults to 1.0, which gives equal weight to both solutions. | ||
| """ | ||
| # Calculate weighted averages of the image and variance planes. | ||
| # Note that ``newModel *= gain`` would multiply the variance by ``gain**2`` |
There was a problem hiding this comment.
If it works, it works, but its is still not obvious WHY multiplying gain**2 and dividing by (1+ gain)**2 would be bad when you're averaging two images:
>>> var1 = 25
>>> var2 = 0.4
>>> s1 = np.random.normal(2, np.sqrt(var1), 10000)
>>> s2 = np.random.normal(3, np.sqrt(var2), 10000)
>>> (var1 + var2)/4.
6.35
>>> np.var(0.5*(s1+s2))
6.4071829805944898
There was a problem hiding this comment.
Using gain**2 would work fine as long as gain=1., but otherwise it biases the variance. And, since this algorithm is iterative that means that it gets progressively higher (or lower) with each iteration. The way it is written ensures that the variance does not drift high or low just from running for more iterations. I'll add more explanation to the comment.
python/lsst/ip/diffim/dcrModel.py
Outdated
| model.image.array[badPixels] = 0. | ||
| model.variance.array[badPixels] = 0. | ||
| model.image.array /= dcrNumSubfilters | ||
| model.variance.array /= dcrNumSubfilters |
There was a problem hiding this comment.
Really divide by dcrNumSubfilters not dcrNumSubfilters**2?
There was a problem hiding this comment.
I believe so, yes. Since the subfilter images are treated as independent, when we divide an image into N subfilters I think each should be treated as having less information, and sqrt(N) worse signal to noise. I will add unit tests in DM-15342 and DM-15343 that will test this and hopefull clarify the expected behavior.
python/lsst/ip/diffim/getTemplate.py
Outdated
| if not sensorRef.datasetExists(**patchArgDict): | ||
| self.log.warn("%(datasetType)s, tract=%(tract)s, patch=%(patch)s does not exist" | ||
| % patchArgDict) | ||
| continue | ||
|
|
||
| nPatchesFound += 1 |
There was a problem hiding this comment.
nPatchesFound was incremented if the patch was found by the butler. This way if there was no data to make a template, it would raise. I think you meant to put your block that reads or builds a patch above nPatchedFound += 0.
python/lsst/ip/diffim/getTemplate.py
Outdated
| self.log.info("Reading patch %s" % patchArgDict) | ||
| coaddPatch = sensorRef.get(**patchArgDict) | ||
| coaddExposure.getMaskedImage().assign(coaddPatch.getMaskedImage(), coaddPatch.getBBox()) | ||
| if self.config.coaddName == 'dcr': |
There was a problem hiding this comment.
I not crazy about the "if dcr do these three things else do these other three things." It's fine for now, but will you help me refactor it per RFC-352 soon?
There was a problem hiding this comment.
I quite agree that the "if dcr..." block is ugly, and I would strongly welcome any suggestions for how to change it. One possibility would be to modify the butler so that coaddPatch = sensorRef.get(**patchArgDict) would know to calculate and return the same thing as the DcrModel calculations, but I worry that would look surprising and magic.
There was a problem hiding this comment.
I am also happy to help implement RFC-352 if that is useful, though DcrModel here is not a Task, so I'm not sure what needs to be done.
There was a problem hiding this comment.
I meant that I plan to refactor GetCoaddAsTemplateTask. I might need to create a DcrModel with your fromImage interface rather than your fromDataRef interface, and might need advice.
There was a problem hiding this comment.
Got it. Yes, I will of course be happy to help with that refactor.
python/lsst/ip/diffim/getTemplate.py
Outdated
| continue | ||
| self.log.info("Reading patch %s" % patchArgDict) | ||
| coaddPatch = sensorRef.get(**patchArgDict) | ||
| coaddExposure.maskedImage.assign(coaddPatch.maskedImage, coaddPatch.getBBox()) |
There was a problem hiding this comment.
Is this equivalent to coaddExposure.maskedImage[coaddPatch.getBBox()] = coaddPatch.maskedImage? If so, that'd be clearer
There was a problem hiding this comment.
I believe it is equivalent in this instance, but I'd like to leave this as it is for consistency with assembleCoadd. If you'd like that changed, I (or someone else) could do that on a separate ticket for both ip_diffim and pipe_tasks, and possibly other packages as well.
isullivan
force-pushed
the
tickets/DM-14738
branch
from
b859a15
to
6d51d89
Compare
No description provided.