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 PSFMap class #1432

Merged
merged 22 commits into from Jul 10, 2018
Merged

Add PSFMap class #1432

merged 22 commits into from Jul 10, 2018

Conversation

@registerrier
Copy link
Contributor

@registerrier registerrier commented Jun 18, 2018

This PR adds a function to project a gammapy.irf.PSF3D on a gammapy.map.
It also adds a class to work with the resulting Map to extract an gammapy.irf.EnergyDependentTablePSF at any position or to produce a PSF kernel following the implementation of PR #1388.
The complete computation of a PSFMap for a collection of observation has to be performed on a higher level maker class. This will be added in a coming PR.

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

@cdeil cdeil left a comment

@registerrier - Thanks!

I've left some inline comments.

More generally, I think extending the class-level docstring a bit to describe what this is would be useful.
An example code snippet how to make one of those things in an Examples section is always gold for users and developers, because then they can just copy & paste that and start to learn / try out a new class in ipython / Jupyter. In this case it could be possible to read a CTA PSF IRF and make such a map with ~ 5 lines?

The other main comment I have is the name. You called the class PSFMap. But is is not a map. It has a map. Looking at how this is used then from callers, you get code like in your test:

pmap = PSFMap(psfmap)

It's easy to be confused in calling code what the "psf_map" and the "pmap" are, no?
As you've probably seen in
https://github.com/gammapy/gammapy/pull/1388/files#diff-cde80870ee7dc037d5c9a8863a870762R15
I put "kernel" in the class name, because for me a fundamental difference wrt the IRF classes is that those represent probability density functions, but this is a discretised probability mass function, i.e. what is commonly called a "kernel".
Suggest to avoid the use of "Map" in the class name here to avoid confusion.

@@ -0,0 +1,189 @@
import numpy as np

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

For new files, always add the two boilerplate lines at the top: license and future inits (copy from any other file).
Same for the new test file you're adding here.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

Done

'PSFMap'
]

def make_psf3d(psf_analytical, rad):

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Suggest to attach this converted as a to_psf3d method on EnergyDependentMultiGaussPSF.

As a standalone function it's hard to find, and the name make_psf3d doesn't capture well that this is really convert_energy_dependent_multi_gauss_psf_to_psf_3d.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

Good point.
I have moved make_psf3d to EnergyDependentMultiGaussPSF.to_psf3d.

rad : `~astropy.unit.Quantity` or `~astropy.coordinates.Angle`
the array of position errors (rad) on which the PSF3D will be defined
Return

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Returns with an s.

Please check once locally that python setup.py build_docs works and the resulting HTML docs formatting is OK for the new functions / classes you're adding.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

Done

----------
psf_analytical : `~gammapy.irf.EnergyDependentMultiGaussPSF`
the analytical PSF to be transformed to PSF3D
rad : `~astropy.unit.Quantity` or `~astropy.coordinates.Angle`

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

units with an s.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

Done

valid = np.where(separations < max_offset)

# Compute PSF values
psf_values = np.transpose(psf.evaluate(offset=separations[valid], energy=energy, rad=rad), axes=(2, 0, 1))

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Suggest to put the transpose call on the next line.
Very long lines are hard to read, and also harder to debug if one ever needs to enter a debugger here.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

Done

normalize : bool
normalize PSF per energy bin.
factor : int
oversampling factor to compute the PSF

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Add empty line before Returns section.

return EnergyDependentTablePSF(energy=self.energies,rad=self.rad,psf_value=psf_values.T)

def get_psf_kernel(self, position, npix, pixel_size, normalize=True, factor=5):
"""Returns a PSF kernel at the given position.

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Always add empty line after the first summary line in a docstring. Then the docs appear nicely in the method summary overview.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

Done

Parameters
----------
fraction : float
the containment fraction (a positive number <=1). Default 0.68.

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Add empty line before Returns section.

This comment has been minimized.

@registerrier

registerrier Jul 4, 2018
Author Contributor

Done

psfmap = make_psf_map(psf, pointing, geom, 3*u.deg)

pmap = PSFMap(psfmap)
# need to assert something...

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Yes, please execute each method (especially get_energy_dependent_table_psf) and add an assert on the output.

return psfmap


class PSFMap():

This comment has been minimized.

@cdeil

cdeil Jun 19, 2018
Member

Write class PSFMap(object).

If you omit the object you'll get an "old-style class" on Python 2. Usually you don't notice, but really those old-style classes should not be used since Python 2.1 or something like that, and then in some cases they do behave in weird ways.

This comment has been minimized.

@registerrier

registerrier Jun 19, 2018
Author Contributor

My bad...

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Jun 19, 2018

Regarding the name, indeed PSFMapis not very satisfactory since it is not a Map. Yet the only thing it contains is a Map. The class mostly provides convenience methods to interact with the map of PSFs, such as the extraction of an EnergyDependentTablePSFat any position. A possibility could be to have only functions taking the actual Map as argument. But this seemed an easier solution for users.

Making an inheritance on Map did not seem like a good idea since the Map could be an WcsNDMapor an HpxNDMap and it is expected to have specific axes.

This object does not contain kernels a priori, but can be used to create them. It could be called PSFKernelMaker or PSFKernelFactory or?

What we are doing here is really to create and handle maps of IRFs (here PSF but the approach would be nearly identical for Edisp). So PSFMapHandler could be a better name?

@cdeil
Copy link
Member

@cdeil cdeil commented Jun 19, 2018

@registerrier - You're right. I was confusing this PSFMap object with what I started to implement in #1388, a real PSF kernel array where I also used a Map object to hold it.

I'm afraid I'm out of good ideas concerning the name. Maybe leave as-is? I certainly didn't want to suggest changing to subclassing or something like that, I was just confused about what this is and was thinking if there is a better name.

@registerrier registerrier moved this from To do to In progress in gammapy.maps Jun 20, 2018
@registerrier registerrier force-pushed the registerrier:irf_map branch from 851459b to 79a6e51 Jul 5, 2018
registerrier added 2 commits Jul 5, 2018
raise ValueError("EnergyDependentTablePSF can be extracted at one single position only.")

# axes ordering fixed. Could be changed.
pix_ener = np.arange(self.psf_map.geom.axes[1].nbin)

This comment has been minimized.

@cdeil

cdeil Jul 5, 2018
Member

Suggest to pick fixed names and always access axes by name, not position index.

Copy link
Member

@cdeil cdeil left a comment

Some quick comments inline.
I'll try this out locally now ...

assert psf_map.geom.axes[0] == rad_axis
assert psf_map.geom.axes[1] == energy_axis

# Check unit

This comment has been minimized.

@cdeil

cdeil Jul 5, 2018
Member

Suggest to remove such comment lines as "check unit", that don't add any extra info.
Please look over the diff and remove such comment lines.

table_psf = self.get_energy_dependent_table_psf(position)
return PSFKernel.from_table_psf(table_psf, geom, max_radius, factor)

def containment_radius_map(self, fraction = 0.68):

This comment has been minimized.

@cdeil

cdeil Jul 5, 2018
Member

How much work is it to implement the containment radius map?
Would be very nice to have, no?
Alternatively, please open an issue for this and label it beginner friendly and we'll try to find soemone to do it next week.

@cdeil cdeil assigned cdeil and unassigned registerrier Jul 10, 2018
@cdeil
cdeil approved these changes Jul 10, 2018
@cdeil
Copy link
Member

@cdeil cdeil commented Jul 10, 2018

@registerrier - Thanks!

@cdeil cdeil merged commit 2f9afdf into gammapy:master Jul 10, 2018
0 of 2 checks passed
0 of 2 checks passed
continuous-integration/appveyor/pr Waiting for AppVeyor build to complete
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
gammapy.maps automation moved this from In progress to Done Jul 10, 2018
@cdeil cdeil changed the title Add functions and class to work with PSFs with gammapy.maps Add PSFMap class Aug 15, 2018
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

2 participants