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

Add a PSFKernel to perform PSF convolution on Maps #1447

Merged
merged 16 commits into from Jul 4, 2018

Conversation

@registerrier
Copy link
Contributor

@registerrier registerrier commented Jun 28, 2018

This is a first commit to discuss.

Two function that take a TablePSF or an EnergyDependentTablePSF as well as a MapGeom are provided and are used to build the PSF kernel map.

The creators on PSFKernel ensure that there is an odd number of pixels for convolution.

The PSFKernel.apply()method is performing the PSF kernel convolution on a Map that has the same geometry.

A number of tests are present already and should be completed.

@registerrier registerrier requested review from cdeil and adonath Jun 28, 2018
@registerrier registerrier self-assigned this Jun 28, 2018
@registerrier registerrier added this to To do in gammapy.maps via automation Jun 28, 2018
@registerrier registerrier added this to the 0.8 milestone Jun 28, 2018
Copy link
Member

@cdeil cdeil left a comment

@registerrier - Thanks!

See my inline comments. Almost all of them are to reduce the amount of code and code duplication. At the moment there's ~ 20 lines of non-trivial code here that's duplicated for the 2D and 3D case for no good reason, and IMO the method to compute the kernel is unnecessarily complex (e.g. the solid angle multiplication), given that a caller has little control e.g. over PSF placement anyways (at least you don't show how to do it by preparing a geom object that's centered on a source or expose a source position option at the moment in the constructor methods).

self._psf_kernel_map = psf_kernel_map

@classmethod
def from_map(cls, psf_kernel_map):

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

This from_map is the same as init -> useless -> remove.

from ..irf import TablePSF

__all__ = [
'table_psf_to_kernel_map',

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

The table_psf_to_kernel_map and energy_dependent_table_psf_to_kernel_map are exposed twice, once as functions and once as class methods on the PSFKernel class.

I'd suggest to remove the standalone functions and paste their code in the classmethods.
They have lots of comments that say exactly what the next line of code says, if you remove most of those comments they aren't that long. (and you get to delete the duplicated long docstring!)
If you still find it too long, I'd suggest to rename the functions to _table_psf_to_kernel_map and _energy_dependent_table_psf_to_kernel_map and not expose them in the docs to users, i.e. have them as private helper functions only.

Thinking about it: a lot of the code is duplicated. Shouldn't the energy-dependent one be a loop over energy bins, and then just call the to_table_psf() method on the e-dep PSF and call the compute kernel for TablePSF?


energy_unit = u.Unit(geom.axes[energy_idx].unit)

# TODO: change the logic to support non-regular geometry. This would allow energy dependent sizes for the kernel.

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Wrap this long comment line into two lines.
Suggest to move the geom check to the top.

This comment has been minimized.

@registerrier

registerrier Jun 29, 2018
Author Contributor

Done

factor : int
the oversample factor to compute the PSF
"""
# Find energy axis in geom

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

I see, the way to access axes is not very nice at the moment.
As a short-term solution, suggest to move this code into a method on geom.get_axis and here just call:

energy_axis = geom.get_axis(type='energy')

that behaves like dict.get. No need for special error handling here, if the get fails, the error is raised in geom.get_axis (I think either KeyError or ValueError would be OK here, not sure which is better).

# get coordinates
map_c = upsampled_image_geom.get_coord()

# get solid angles

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Is it really useful to do this complicated PSF evaluation and solid angle weighting?
My suggestion would have been to just use a pixel scale, like we did so far in Gammapy Kernel computations.
Example: https://github.com/gammapy/gammapy/pull/1388/files#diff-cde80870ee7dc037d5c9a8863a870762R76

If you want to put this computation, please don't duplicate it twice. do it just in one energy band, and from the energy-dependent method, just loop over energy bands and call the other method.

This comment has been minimized.

@registerrier

registerrier Jun 29, 2018
Author Contributor

I think this computation is not really complex in terms of code, although likely heavier than the previous methods used for kernel computation. The computing time is far from being a bottleneck anyway.
The interesting point to me is that using the ~gammapy.maps framework potentially offers real power with simple code. e.g. if we want to have energy dependent kernel size, the change should be minimal.

factor : int
the oversample factor to compute the PSF
"""
if max_radius is not None:

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Avoid the code duplication: you have the same if max_radius and code in from_energy_dependent_table_psf and from_table_psf. Move it to a helper function and call that twice.

This comment has been minimized.

@registerrier

registerrier Jul 2, 2018
Author Contributor

Done

# maximum at the center of map?
ind = np.unravel_index(np.argmax(kernel.data, axis=None), kernel.data.shape)
# absolute tolerance at 0.5 because of even number of pixel here
assert_allclose(ind, geom.center_pix, atol=0.5)

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Those two asserts are pretty indirect. Suggest this instead:

assert kernel.data.shape == (some, value)
assert kernel.data[center, pixel] == kernel.data.max()

I.e. put the shape values and center pixel index values, isntead of this atol=0.5

Just computing a kernel that's 3 x 3 and putting the 3x3 float values in the assert might also be an option.


kernel = PSFKernel.from_gauss(testmap.geom, sigma)

conv_map = kernel.apply(testmap)

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Does this FFT convolve work properly?
Suggest to do a e.g. (10 x 5) map (different nx / ny) and put a value of 1 or so in the center and at the edges or corners, and to assert on the resulting pixels at the same positions that come out.
Something that really establishes what the PSF apply does, also at the edge of the image.

This is a container class to store a PSF kernel
that can be used to convolve `~gammapy.maps.WcsNDMap` objects.
It is usually computed from an EnergyDependentTablePSF

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Suggest to either just remove this "It is usually computed from an EnergyDependentTablePSF", and mention the classmethods instead that are usually used to create a kernel. Or add an Examples section showing the user how to make a kernel. Having a minimal example one can copy & paste is really nice to learn a new class.
See e.g. http://docs.gammapy.org/dev/maps/index.html#constructing-with-factory-methods with nice examples how to make maps for testing and learning quickly.



class PSFKernel(object):
"""PSF kernel for `~gammapy.maps.Map`.

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Does this class only support 2D and 3D with energy axis?
Or are extra axes also tolerated and work to some extent?
Suggest to clarify in the docstring (but not spend any time to make extra axes work if they don't already).

This comment has been minimized.

@registerrier

registerrier Jun 28, 2018
Author Contributor

Extra axes are also tolerated but filled the same way.
Docstring clarified.

@registerrier registerrier moved this from To do to In progress in gammapy.maps Jun 29, 2018

# loop over images and fill map
for img, idx in kernel_map.iter_by_image():
img += vals.value

This comment has been minimized.

@adonath

adonath Jun 29, 2018
Member

Normalization of the kernels could be added here to avoid the second Python loop at the end.

return axis
raise ValueError("Cannot find axis named {}".format(name))

def get_axes_by_type(self, type):

This comment has been minimized.

@cdeil

cdeil Jul 2, 2018
Member

Is it really useful to try and support maps with multiple axes of a given type?
I would suggest to change to get_axis_by_type to get one, for symmetry with get_axis_by_name and also it's what you need in the one caller for this functionality.

Note that filtering the axes list by type is still easily possible with a one-liner, should we ever need this.

axes = [ax for ax in m.geom.axes if ax.type == 'energy']
registerrier added 6 commits Jul 2, 2018
@cdeil cdeil assigned cdeil and unassigned registerrier Jul 4, 2018
def read(cls, *args, **kwargs):
"""Read kernel Map from file."""

psf_kernel_map = WcsNDMap.read(*args, **kwargs)

This comment has been minimized.

@cdeil

cdeil Jul 4, 2018
Member

In PyCharm I see: Unresolved reference WcsNDMap
You're missing the import, and this will give a NameError if someone tries to read.
It also means this method isn't covered by a test.
Can you please add one?

This comment has been minimized.

@registerrier

registerrier Jul 4, 2018
Author Contributor

OK. Probably best to use Map.read to be more general. I'll do the change now.

This comment has been minimized.

@cdeil

cdeil Jul 4, 2018
Member

FYI: If you want a coverage report locally to see which methods / lines are covered by your tests, use this command:

python setup.py test -V --coverage -t gammapy/cube/tests/test_psf_kernel.py
open htmlcov/index.html
# click on `cube/psf_kernel.py` to see the coverage report for that file.

@classmethod
def from_gauss(cls, geom, sigma, max_radius=None, containment_fraction=0.99,
normalize=True, factor=4):

This comment has been minimized.

@cdeil

cdeil Jul 4, 2018
Member

PyCharm shows: Unused options normalize and factor.
I think normalize should be removed (also from the docstring), and factor passed along.

geom = WcsGeom.create(binsz=0.02*u.deg, width=10.0*u.deg, axes=[energy_axis])
# Extract EnergyDependentTablePSF from PSF IRF (here a PSF3D)
table_psf = psf.to_energy_dependent_table_psf(theta=0.5*u.deg)

This comment has been minimized.

@cdeil

cdeil Jul 4, 2018
Member

Can you fix this and make a minimal working example?

In [1]: %paste
        import numpy as np
        from gammapy.maps import WcsGeom, MapAxis
        from astropy import units as u

        # Define energy axis
        energy_axis = MapAxis.from_edges(np.logspace(-1., 1., 4), unit='TeV', name='energy')

        # Create WcsGeom
        geom = WcsGeom.create(binsz=0.02*u.deg, width=10.0*u.deg, axes=[energy_axis])

        # Extract EnergyDependentTablePSF from PSF IRF (here a PSF3D)
        table_psf = psf.to_energy_dependent_table_psf(theta=0.5*u.deg)

        psf_kernel = PSFKernel.from_table_psf(table_psf,geom, max_radius=1*u.deg)

        some_map_convolved = psf_kernel.apply(some_map)
## -- End pasted text --
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-0b05c82cb4ce> in <module>()
     10 
     11 # Extract EnergyDependentTablePSF from PSF IRF (here a PSF3D)
---> 12 table_psf = psf.to_energy_dependent_table_psf(theta=0.5*u.deg)
     13 
     14 psf_kernel = PSFKernel.from_table_psf(table_psf,geom, max_radius=1*u.deg)

NameError: name 'psf' is not defined

This comment has been minimized.

@registerrier

registerrier Jul 4, 2018
Author Contributor

Done. It works now.

@cdeil
cdeil approved these changes Jul 4, 2018
@cdeil cdeil merged commit 2437fa5 into gammapy:master Jul 4, 2018
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
gammapy.maps automation moved this from In progress to Done Jul 4, 2018
@cdeil
Copy link
Member

@cdeil cdeil commented Jul 4, 2018

@registerrier - Thank you very much!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
gammapy.maps
  
Done
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants