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 map geom pixel solid angle computation #1357

Merged
merged 6 commits into from Apr 11, 2018

Conversation

Projects
None yet
2 participants
@registerrier
Contributor

registerrier commented Apr 4, 2018

This adds a Geom.solid_angle() method to compute the solid angles of the map pixels.
For healpix geometry this simply returns a 1 element Quantity.
For Wcs Geom it returns a Quantity array with the same dimension as the Geom itself.

@registerrier registerrier added the feature label Apr 4, 2018

@registerrier registerrier added this to the 0.8 milestone Apr 4, 2018

@registerrier registerrier requested review from cdeil and woodmd Apr 4, 2018

@cdeil cdeil self-assigned this Apr 4, 2018

@cdeil cdeil changed the title from add solid angle computation for Geom objects to Add map geom pixel solid angle computation Apr 4, 2018

@cdeil

@registerrier - Please add at least one test (e.g. with 3x2 pix) case that exercises the new code, to show that it's working. Ideally there would even be a second test case where you have an energy or some extra axis, to show if the code still works and if it then returns a 2D or 3D array in that case (add an assert on shape of return value).

Probably http://docs.gammapy.org/dev/_modules/gammapy/image/core.html#SkyImage.solid_angle already has such a test case that you could copy and adapt here.

I'm not sure if it's useful to have the stub in the base class, I think you could also remove it there, as you like.

For all docstrings, you should give a one-line description so that it shows up nicely in the API summary lists. For cases like this one I'd suggestion to put a single line instead of a Returns section as explained here:
http://docs.gammapy.org/dev/development/howto.html#functions-or-class-methods-that-return-a-single-object

I see that the current
http://docs.gammapy.org/dev/_modules/gammapy/maps/wcs.html#WcsGeom.get_coord
doesn't expose mode='edge' and that _get_pix_coords is non-ideal. Here's what we had in the old class:
http://docs.gammapy.org/dev/api/gammapy.image.SkyImage.html#gammapy.image.SkyImage.coordinates

@registerrier - Presumably you don't want to improve the functionality to generate coordinate grids in this PR? If so, please change the comments in your implementation a bit to make it clearer that there's a TODO and what should be done. I'd suggest this:

# TODO: the following lines should be improved to be more clean and efficient.
# E.g. we could expose `mode='edge'` directly on `WcsGeom.get_coord`

Efficiency is not really a concern, because coordinate and solid angle arrays are only computed rarely. But IMO we should avoid re-implementing things, so I'm -1 on your proposal in the current code to call np.mgrid and self._wcs.all_pix2world directly from solid_angle.

@registerrier

This comment has been minimized.

Contributor

registerrier commented Apr 4, 2018

Some test added and some cleaning as suggested.

Note that solid_angle is called at least for each background map processed. This can be a significant cost if inefficient.

@cdeil

@registerrier - I've done a second round of review.

To significant suggestions:

  1. return a 2D array only
  2. change test case to AITOFF with large pixels, and check that something sensible happens for a pixel that is outside or partially outside the image.
@@ -1296,6 +1296,17 @@ def upsample(self, factor):
"""
pass
@abc.abstractmethod
def solid_angle(self):
"""

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

Adjust docstring here as well. Which docstring ends up in the derived classes in the HTML docs?

@@ -1712,6 +1713,17 @@ def skydir_to_pix(self, skydir):
return self.pix_to_coord((lat, lon))
def solid_angle(self):
"""

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

Adjust docstring here as well (always add a one-line summary).

solid_array = geom.solid_angle()
# Check array size
assert solid_array.shape == (9,3,100,100)

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

Maybe make the solid angle array 2D to avoid wasting memory by duplicating the same values for the extra axes?
I would do that and then leave broadcasting against extra axes (if needed) to callers using the solid angle arrays.

I think there is a method on geom to get the geom for the spatial part only, so should be easy to implement, no?

This comment has been minimized.

@registerrier

registerrier Apr 5, 2018

Contributor

Well, I initially did this but it seems more straightforward to do:
geom.to_image().solid_angle()
if one needs only the 2D data only rather than have to deal with the resize & broadcast oneself.

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

I think this could be implemented either way, i.e. only return 2D or ND or add a keyword argument so that the caller can get either version.

But in any case, mentioning geom.to_image().solid_angle() in the docstring here would be appropriate, as well as how to work with the memory-efficient way, i.e. get a 2D solid angle array, and then multiply it against an ND map using numpy broadcasting or some other way to get multiplication of 2D array with ND map to work, i.e. what one needs for background, model evaluation, ... to multiply with this solid angle array.

Can you please decide on the implementation, but add a few lines to the docstring to explain how to work with this?

This comment has been minimized.

@registerrier

registerrier Apr 5, 2018

Contributor

A possibility could be to pass an option to geom.solid_angle(mode='full') to have the full ND array and geom.solid_angle(mode='2D') to have only the 2D array.
Does it sound better like this?

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

I was going to suggest to always return 2D and only explain in the docstring how to multiply that with an ND map (or if really needed, convert it to ND with one or two lines of Numpy/ Maps code).

But now I realise that I'm forgetting about the case of variable map pixel npix size.
I think those multi-resolution maps already exist, see the last example in this section:
http://docs.gammapy.org/dev/maps/index.html#constructing-with-factory-methods
and they require an ND (or at least > 2D) solid angle array.

Does that case work in your current implementation?
If yes, can you add an extra test case to show that it does?
E.g. two bins that differ by factor 2, and sum of solid angles in each plane should be almost the same, or assert just on one pixel from each plane, showing that it's roughly factor 4 different in solid angle.

@registerrier - I realise that my comments are always towards extending implementation and tests, but really we're better served by small and quick pull requests to have the basic 3D working, and then to extend later. So please finish up the PR as you think best, possibly leaving TODO or other comments for us with ideas what to improve at a later time.

npix[0] += 1
npix[1] += 1
npix[0][0] += 1
npix[1][0] += 1

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

What changed here? Was it a bug? Add a regression test?

This comment has been minimized.

@registerrier

registerrier Apr 5, 2018

Contributor

geom.npix is a tuple of arrays. The current implementation does not work.
Maybe a better implementation would be npix[0] = npix[0] + 1

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

If npix[0] is a Numpy array, then npix[0][0] += 1 and npix[0] = npix[0] + 1 are the same, no? So pick either implementation. If you have time, please add a short regression test. Or I can also add it before I merge this in.

This comment has been minimized.

@registerrier

registerrier Apr 5, 2018

Contributor

This does not work either (tuple is immutable). Did this instead:

 for pix_num in npix[:2]:
                pix_num += 1
@@ -594,6 +596,30 @@ def upsample(self, factor):
def to_slice(self, slices):
raise NotImplementedError
def solid_angle(self):
""" Returns solid angle array (in `sr`) as `~astropy.units.Quantity`

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

Suggest to mention shape of returned array in an extra line. I.e. docs shoud say if this is 2D or higher dimensional if there are extra axes.

def test_wcsgeom_solid_angle():
binsz=1.0
npix=100

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

Why make so many pixels?
In my experience test with large arrays (like 100 x 100 here) are just harder to debug when things go wrong (e.g. when you see these arrays printed out in the tests reports on travis-ci), and, less importantly, are a bit slower to run.
I'd suggest to just make a few (and increase binsz to e.g. 10 or 20 if you want to go high-latitude.

assert solid_array.shape == (9,3,100,100)
# Test bin size at b=0
solid_lat0 = binsz*np.pi/180 * (np.sin(binsz * np.pi / 180.)) * u.sr

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

I'm -1 to have such code, that re-implements the CAR projection, in our tests. If you look around in Gammapy, you will find 1000s of lines of such complex tests, that are hard to understand. In my experience, this is more of a problem than an asset.

I would just hard-code the reference values of the pixels in the test.
There could still be a comment saying if the reference values just establish the current result, or if they can be understood. In this case, I think a value of very roughly binsz ** 2 is expected, so that could be mentioned in a comment.

One more thing that might be worth testing is what happens for pixels "outside the image", e.g. at the edge for AIT projection. Ideally the solid angle function would return NaN there and not error out. For partially outside pixels with e.g. one corner outside, ideally we'd get NaN also or some useful value, and not e.g. a crazy large or small value because the algorithm fails and then that might cause problems with analysis in those cases.

This comment has been minimized.

@registerrier

registerrier Apr 5, 2018

Contributor

I beg to differ. Giving directly the numerical value does not give any information on how it was extracted.

OK for the test with AIT external pixels.

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

In case you haven't seen it: when we implemented SkyImage, we already spent quite some time trying to get useful test cases and to validate the coordinates and solid angle methods:

class _TestImage:

My suggestion would have been to simply copy and adapt those to the new Map API, instead of trying to develop something new.
Although by now, I have changed my opinion a bit on how tests should be structured, preferring as simple and concrete tests over e.g. the way you're writing the tests here now (adding lines that re-implement CAR projection) or like we did in many other places in Gammapy, e.g. the complex way the test parametrisation was done via test_image for SkyImage; I already find that hard to read and maintain and would prefer a little bit of code duplication and having the reference values directly in the test function, i.e. to have two or three test functions each with one test case here.

pix = self._get_pix_coords(mode='edge')
# Note also that pix_to_coord is already called in _get_pix_coords. This should be made more efficient.
coords = self.pix_to_coord(pix)
lon=coords[0] * np.pi/180.

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

We shouldn't litter our code with np.pi/180 a lot.
Ideally coords would have methods / properties that allow you to do something like lon = coords.lon.to('rad').
Does it or is this easy to add?

This comment has been minimized.

@registerrier

registerrier Apr 5, 2018

Contributor

Well coords is simply a tuple of arrays, so one needs to create a SkyCoord. But then you have to know which coordinate system you are working with to do either coords.lor coords.ra.
The function ~gammapy.maps.geom.skycoord_to_lonlatdoes it already but it does not return a Quantityor an Angle but directly an array in degrees.
A possibility then is to add an option to this function to return radians instead of degrees.

It all boils down to a better integration of units into gammapy.maps.

This comment has been minimized.

@cdeil

cdeil Apr 5, 2018

Member

OK, then let's defer this point to another day / PR.
As you know, I'm +1 to use MapCoords objects more instead of raw tuples.

@cdeil

cdeil approved these changes Apr 11, 2018

@registerrier - I've done some edits in eaaaad9 .

I'm OK to merge as-is now, although it would be nice to make the docstring / test clearer on what shape array is returned for the HPX case.

It seems to be scalar, but I think you said it's array?

@registerrier

This comment has been minimized.

Contributor

registerrier commented on gammapy/maps/hpx.py in eaaaad9 Apr 11, 2018

Here indeed we return a simple number with a unit unless the maps have different nside in which case we return an array.
This might break the common usage with the WcsGeom.solid_angle().

@registerrier

This comment has been minimized.

Contributor

registerrier commented on eaaaad9 Apr 11, 2018

I am OK to merge.

@cdeil cdeil merged commit e1ee3bd into gammapy:master Apr 11, 2018

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Apr 11, 2018

@registerrier - Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment