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

Tickets/dm 9229 #54

Merged
merged 4 commits into from
Apr 3, 2017
Merged

Tickets/dm 9229 #54

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a7e9e79
Remove unnecessary comment
yalsayyad Mar 27, 2017
02a9d97
Remove config padPsf
yalsayyad Mar 27, 2017
4b3ab69
Pad science PSF if too small for matching kernel in ModelPsfMatchTask
yalsayyad Mar 28, 2017
edbf328
import numpy as np per RFC-34
yalsayyad Apr 3, 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
99 changes: 67 additions & 32 deletions python/lsst/ip/diffim/modelPsfMatch.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
from __future__ import absolute_import, division, print_function

from builtins import range
import numpy as num
import numpy as np

from . import diffimLib
import lsst.afw.geom as afwGeom
Expand All @@ -36,7 +36,14 @@
from . import utils as diUtils
import lsst.afw.display.ds9 as ds9

sigma2fwhm = 2. * num.sqrt(2. * num.log(2.))
__all__ = ("ModelPsfMatchTask", "ModelPsfMatchConfig")

sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))


def nextOddInteger(x):
nextInt = int(np.ceil(x))
return nextInt + 1 if nextInt % 2 == 0 else nextInt


class ModelPsfMatchConfig(pexConfig.Config):
Expand All @@ -49,11 +56,27 @@ class ModelPsfMatchConfig(pexConfig.Config):
),
default="AL",
)

padPsf = pexConfig.Field(
doAutoPadPsf = pexConfig.Field(
dtype=bool,
doc="Grow PSF dimensions to largest PSF dimension on input exposure; else clip to smallest",
default=False,
doc=("If too small, automatically pad the science Psf? "
"Pad to smallest dimensions appropriate for the matching kernel dimensions, "
"as specified by autoPadPsfTo. If false, pad by the padPsfBy config."),
default=True,
)
autoPadPsfTo = pexConfig.RangeField(
dtype=float,
doc=("Minimum Science Psf dimensions as a fraction of matching kernel dimensions. "
"If the dimensions of the Psf to be matched are less than the "
"matching kernel dimensions * autoPadPsfTo, pad Science Psf to this size. "
"Ignored if doAutoPadPsf=False."),
default=1.4,
min=1.0,
max=2.0
)
padPsfBy = pexConfig.Field(
dtype=int,
doc="Pixels (even) to pad Science Psf by before matching. Ignored if doAutoPadPsf=True",
default=0,
)

def setDefaults(self):
Expand Down Expand Up @@ -288,11 +311,11 @@ def run(self, exposure, referencePsfModel, kernelSum=1.0):
spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)

if psfMatchingKernel.isSpatiallyVarying():
sParameters = num.array(psfMatchingKernel.getSpatialParameters())
sParameters = np.array(psfMatchingKernel.getSpatialParameters())
sParameters[0][0] = kernelSum
psfMatchingKernel.setSpatialParameters(sParameters)
else:
kParameters = num.array(psfMatchingKernel.getKernelParameters())
kParameters = np.array(psfMatchingKernel.getKernelParameters())
kParameters[0] = kernelSum
psfMatchingKernel.setKernelParameters(kParameters)

Expand Down Expand Up @@ -344,13 +367,6 @@ def _buildCellSet(self, exposure, referencePsfModel):
dimenR = referencePsfModel.getLocalKernel().getDimensions()
psfWidth, psfHeight = dimenR

maxKernelSize = min(psfWidth, psfHeight) - 1
if maxKernelSize % 2 == 0:
maxKernelSize -= 1
if self.kConfig.kernelSize > maxKernelSize:
raise ValueError("Kernel size (%d) too big to match Psfs of size %d; reduce to at least %d" % (
self.kConfig.kernelSize, psfWidth, maxKernelSize))

regionSizeX, regionSizeY = scienceBBox.getDimensions()
scienceX0, scienceY0 = scienceBBox.getMin()

Expand All @@ -368,7 +384,6 @@ def _buildCellSet(self, exposure, referencePsfModel):

# Survey the PSF dimensions of the Spatial Cell Set
# to identify the minimum enclosed or maximum bounding square BBox.
# This loop can be made faster with a cheaper method to inspect PSF dimensions
widthList = []
heightList = []
for row in range(nCellY):
Expand All @@ -379,15 +394,42 @@ def _buildCellSet(self, exposure, referencePsfModel):
widthList.append(widthS)
heightList.append(heightS)

lenMax = max(max(heightList), max(widthList))
lenMin = min(min(heightList), min(widthList))
psfSize = max(max(heightList), max(widthList))

lenPsfScience = lenMax if self.config.padPsf else lenMin
dimenS = afwGeom.Extent2I(lenPsfScience, lenPsfScience)
if self.config.doAutoPadPsf:
minPsfSize = nextOddInteger(self.kConfig.kernelSize*self.config.autoPadPsfTo)
Copy link
Contributor

Choose a reason for hiding this comment

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

It may not be necessary, as this function is only called a small number of times, but it is possible to avoid all the overhead of calling a function, and the steps in the function by saying:
minPsfSize = int(np.ceil(self.KConfig.kernelSize*self.config.autoPadPsfTo))|1
This will return the number if it is odd, or the next integer if it is even

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I acknowledge two suggestions here: (1) to take the minPsfSize calculation out of a function and (2) to switch the last bit to 1 to get the next odd integer.
1: This function is not in a loop and only called once. My rationale for putting it in a function was to reduce code duplication with the unit test. If you argue that unit tests should calculate it differently, then I'll take it out of the function.
2: Flipping the last bit to 1 smells clever to me, because it doesn't read like english.

The function plus english-reading % 2 version also appears to be (counterintuitively) faster, but a microsecond (or even 100 milliseconds) won't make a difference in a long task like ModelPsfMatch:

import timeit

setupBitFlip = """
from random import random
import numpy as np
def nextOddInt(x):
    return int(np.ceil(x))|1

def genKRand():
    for i in range(1000):
        yield 10*random()
"""

setupIfMod = """
from random import random
import numpy as np

def genKRand():
    for i in range(1000):
        yield 10*random()

def nextOddInt(x):
    return x if x % 2 == 0 else x + 1
"""

print('%.02f us'% (1000*timeit.timeit('[nextOddInt(y) for y in genKRand()]', setup=setupIfMod, number=10)/10.))
# 0.63 us
print('%.02f us'% (1000*timeit.timeit('[nextOddInt(y) for y in genKRand()]', setup=setupBitFlip, number=10)/10.))
# 1.26 us
print('%.02f us'% (1000*timeit.timeit('[int(np.ceil(y))|1 for y in genKRand()]', setup=setupBitFlip, number=10)/10.))
# 1.07 us

# Single value
print('%.02f us'%(timeit.timeit('nextOddInt(3.4)', setup=setupIfMod, number=1000000)))
# 0.35 us
print('%.02f us'%(timeit.timeit('nextOddInt(3.4)', setup=setupBitFlip, number=1000000)))
# 0.99 us
print('%.02f us'%(timeit.timeit('int(np.ceil(3.4))|1', setup=setupBitFlip, number=1000000)))
# 0.82 us

Copy link
Contributor

Choose a reason for hiding this comment

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

The speed difference here comes from setupBitFlip using int(np.ceil( and setupIfMod not, it seems to just be doing the mod on whatever type that x is. Over all I don't have a strong feeling on this and reuse in the unit test is a good enough argument for me.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah, whoops. Let me update those last times for posterity:

# # Single value UPDATED
print('%.02f us'%(timeit.timeit('nextOddInt(3.4)', setup=setupIfMod, number=1000000)))
# 1.23 us
print('%.02f us'%(timeit.timeit('nextOddInt(3.4)', setup=setupBitFlip, number=1000000)))
# 1.02 us
print('%.02f us'%(timeit.timeit('int(np.ceil(3.4))|1', setup=setupBitFlip, number=1000000)))
# 0.87 us

paddingPix = max(0, minPsfSize - psfSize)
else:
if self.config.padPsfBy % 2 != 0:
raise ValueError("Config padPsfBy (%i pixels) must be even number." %
Copy link
Contributor

Choose a reason for hiding this comment

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

Is it worth just moving to the next even number rather than throwing an error? Or maybe being explicit is good, like you have it. If the algorithm must have an even to work, perhaps move to the next even number and put out a warning? What you have is fine, these are just some things for you to consider.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I considered that. My reasoning for raising an error is that when setting doAutoPadPsf=False, I presume the user has reasons for turning off the automatic tuning, and thus wants exact manual control over the configs and wants ModelPsfMatchTask to just do exactly what they ask of it or quit. The behavior is consistent: Kernel too big? Raise. Padding inappropriately odd? Raise. I suspect this config option will only be used when users want precise control such as for testing config parameters like you did for the HSC image

self.config.padPsfBy)
paddingPix = self.config.padPsfBy

if paddingPix > 0:
self.log.info("Padding Science PSF from (%s, %s) to (%s, %s) pixels" %
(psfSize, psfSize, paddingPix + psfSize, paddingPix + psfSize))
psfSize += paddingPix

# Check that PSF is larger than the matching kernel
maxKernelSize = psfSize - 1
if maxKernelSize % 2 == 0:
maxKernelSize -= 1
if self.kConfig.kernelSize > maxKernelSize:
message = """
Kernel size (%d) too big to match Psfs of size %d.
Please reconfigure by setting one of the following:
1) kernel size to <= %d
2) doAutoPadPsf=True
3) padPsfBy to >= %s
""" % (self.kConfig.kernelSize, psfSize,
maxKernelSize, self.kConfig.kernelSize - maxKernelSize)
raise ValueError(message)

dimenS = afwGeom.Extent2I(psfSize, psfSize)

if (dimenR != dimenS):
self.log.info("Adjusting dimensions of reference PSF model from %s to %s" % (dimenR, dimenS))
referencePsfModel = referencePsfModel.resized(lenPsfScience, lenPsfScience)
referencePsfModel = referencePsfModel.resized(psfSize, psfSize)
dimenR = dimenS

policy = pexConfig.makePolicy(self.kConfig)
Expand Down Expand Up @@ -417,17 +459,10 @@ def _buildCellSet(self, exposure, referencePsfModel):
else:
# make image of proper size
kernelImageS = afwImage.ImageF(dimenR)
if self.config.padPsf:
bboxToPlace = afwGeom.Box2I(afwGeom.Point2I((lenMax - rawKernel.getWidth())//2,
(lenMax - rawKernel.getHeight())//2),
rawKernel.getDimensions())
kernelImageS.assign(rawKernel, bboxToPlace)
else:
bboxToExtract = afwGeom.Box2I(afwGeom.Point2I((rawKernel.getWidth() - lenMin)//2,
(rawKernel.getHeight() - lenMin)//2),
dimenR)
rawKernel.setXY0(0, 0)
kernelImageS = rawKernel.Factory(rawKernel, bboxToExtract)
bboxToPlace = afwGeom.Box2I(afwGeom.Point2I((psfSize - rawKernel.getWidth())//2,
(psfSize - rawKernel.getHeight())//2),
rawKernel.getDimensions())
kernelImageS.assign(rawKernel, bboxToPlace)

kernelMaskS = afwImage.MaskU(dimenS)
kernelMaskS.set(0)
Expand Down
37 changes: 37 additions & 0 deletions tests/testModelPsfMatch.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
import lsst.ip.diffim as ipDiffim
import lsst.log.utils as logUtils
import lsst.meas.algorithms as measAlg
from lsst.ip.diffim.modelPsfMatch import nextOddInteger
Copy link
Contributor

Choose a reason for hiding this comment

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

If you changed to the binary operator in the other file, remember to change where the function is used here as well


logUtils.traceSetAt("ip.diffim", 4)

Expand All @@ -28,6 +29,10 @@ def setUp(self):
self.exp.setPsf(measAlg.DoubleGaussianPsf(self.ksize, self.ksize, self.sigma1))

def testTooBig(self):
"""Test that modelPsfMatchTask raises if kernel size is too big and
and automatic padding disabled
"""
self.config.doAutoPadPsf = False
self.subconfig.kernelSize = self.ksize
psf = measAlg.DoubleGaussianPsf(self.ksize, self.ksize, self.sigma2)
psfMatch = ipDiffim.ModelPsfMatchTask(config=self.config)
Expand Down Expand Up @@ -61,13 +66,45 @@ def testAdjustModelSize(self):
"""Test that modelPsfMatch correctly adjusts the model PSF dimensions to
match those of the science PSF.
"""
self.config.doAutoPadPsf = False
psfModel = measAlg.DoubleGaussianPsf(self.ksize + 2, self.ksize + 2, self.sigma2)
psfMatch = ipDiffim.ModelPsfMatchTask(config=self.config)
results = psfMatch.run(self.exp, psfModel)
self.assertEqual(results.psfMatchedExposure.getPsf().computeImage().getDimensions(),
self.exp.getPsf().computeImage().getDimensions())
self.assertEqual(results.psfMatchedExposure.getPsf().getSigma1(), self.sigma2)

def testPadPsf(self):
"""Test automatic and manual PSF padding

Compare expected Psf size, after padding, to the reference Psf size.
The reference Psf Size is proxy for the Sciencee Psf size here.
"""
psfModel = measAlg.DoubleGaussianPsf(self.ksize, self.ksize, self.sigma2)

# Test automatic padding (doAutoPadPsf is True by default)
autoPaddedKernel = nextOddInteger(self.subconfig.kernelSize*self.config.autoPadPsfTo)
psfMatch = ipDiffim.ModelPsfMatchTask(config=self.config)
results = psfMatch.run(self.exp, psfModel)
self.assertEqual(results.psfMatchedExposure.getPsf().computeImage().getWidth(), autoPaddedKernel)

# Test manual padding
self.config.doAutoPadPsf = False
PAD_EVEN_VALUES = [0, 2, 4]
for padPix in PAD_EVEN_VALUES:
self.config.padPsfBy = padPix
psfMatch = ipDiffim.ModelPsfMatchTask(config=self.config)
results = psfMatch.run(self.exp, psfModel)
self.assertEqual(results.psfMatchedExposure.getPsf().computeImage().getWidth(),
self.ksize + padPix)

PAD_ODD_VALUES = [1, 3, 5]
for padPix in PAD_ODD_VALUES:
self.config.padPsfBy = padPix
psfMatch = ipDiffim.ModelPsfMatchTask(config=self.config)
with self.assertRaises(ValueError):
results = psfMatch.run(self.exp, psfModel)

def tearDown(self):
del self.exp
del self.subconfig
Expand Down