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 function to fill acceptance image from curve #248

Merged
merged 4 commits into from Apr 9, 2015

Conversation

Projects
None yet
2 participants
@mapazarr
Member

mapazarr commented Mar 27, 2015

This pull request addresses issue #5.

The method is slow (loop over pixels for calculating offset and interpolate the acceptance value), so maybe it could still be improved.

I still am thinking about a good test for the method. While developping the method I used a script that plots the radial acceptance image that I can see, but I am not sure how to best implement an automatic test on a method working on images. In addition, since the mthod is slow, the test might take longer than desirable. Suggestions are welcome!

@mapazarr

This comment has been minimized.

Member

mapazarr commented Mar 27, 2015

I checked the output from travis: the failed status seems to come from an error in the connection to the repository.

def fill_acceptance_image(image, center, offset, acceptance):
"""Generate a 2D image of a radial acceptance curve given as

This comment has been minimized.

@cdeil

cdeil Mar 27, 2015

Member

Try to have a one-line description ... it then shows up nicer in the overview docs.
In this case.

Fill acceptance image from radial acceptance curve.

... detailed description if needed ...

This comment has been minimized.

@mapazarr

mapazarr Apr 1, 2015

Member

Fixed.

Parameters
----------
image : ImageHDU

This comment has been minimized.

@cdeil

cdeil Mar 27, 2015

Member

Try to write the data types like this so that in the HTML docs links are generated:

image : `~astropy.io.fits.ImageHDU`

This comment has been minimized.

@mapazarr

mapazarr Apr 1, 2015

Member

Fixed.

----------
image : ImageHDU
Empty image to fill.
center : SkyCoord

This comment has been minimized.

@cdeil

cdeil Mar 27, 2015

Member

Same here:

center : `~astropy.coordinates.SkyCoord`
Empty image to fill.
center : SkyCoord
Coordinate of the center of the image.
offset : `~numpy.ndarray`

This comment has been minimized.

@cdeil

cdeil Mar 27, 2015

Member

We should use Astropy quantity objects as much as possible.
FYI: this is not the case in Astropy core because most of that code was written before astropy.units existed, but the guideline is that new code should use quantities where possible / useful.

In this case offset should be an Angle object, which is a Quantity sub-class.
Example: https://gammapy.readthedocs.org/en/latest/api/gammapy.irf.EnergyDependentTablePSF.html

This comment has been minimized.

@mapazarr

mapazarr Apr 1, 2015

Member

Fixed.

Returns
-------
image : ImageHDU

This comment has been minimized.

@cdeil

cdeil Mar 27, 2015

Member
image : `~astropy.io.fits.ImageHDU`

This comment has been minimized.

@mapazarr

mapazarr Apr 1, 2015

Member

Fixed.

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 27, 2015

@mapazarr - Thanks!

As you say, a test needs to be added (do this first) and the function needs to be re-written so that it's fast (do this second).

Some suggestions for the test:

  • Use make_empty_image to get an input image.
  • Create some offset and acceptance objects from scratch, or add a utliity function to make it easy to create them for testing (see example here where I've done something similar).

To make it fast, you have to avoid the for loop in Python and instead use numpy arrays.
Use coordinates to get 2D coordinate images and then from that compute a 2D offset image and then evaluate the acceptance curve. For that last part you have to use interpolation, I think, have a look here.

Some info about pyfact.py: it's old spagetti code with useful functionality, but it has to be re-written or re-factored completely and tests / docs added ... this will be a large part of your work ... in most cases it will be simpler if you start from scratch and we just remove pyfact.py at some point when the equivalent functionality is in place. Here's an issue listing what remains to be done to integrate the old pyfact code in to gammapy: #78

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 27, 2015

If you have time to learn more about how to write array-oriented numpy code that avoids for loops in Python, have a look here or here (there's many more resources, but these are pretty good and with excercises I think)

@mapazarr

This comment has been minimized.

Member

mapazarr commented Mar 29, 2015

About the test. I already had implemented what you say. I uploaded the script that I'm using to an external github repository:
https://github.com/mapazarr/astropy_scripts/blob/master/astropy_scripts/test_fill_acceptance_image.py
My question is: since the test should run automatically, I need to define some "assert", but I don't know the best strategy for it. Should I for instance recompute the interpolation of the acceptance in the test script and check a few bins of the image? (it would use the same implementation as in the method to test, so I'm not sure if the test would be meaningful.)
Another possibility would be to fit the image and check that it delivers a radial function with the same parameters as the one used for the radial acceptance passed as input. This one, with the 2D fit, seems a bit too complicated.

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 29, 2015

@mapazarr – In my experience, writing the asserts for tests of analysis code like most of the code in Gammapy requires roughly as much (sometimes even more) work / thought as implementing the algorithms in the first place.

As you say, duplicating the method implementation in the test doesn't make sense.
Here's a few thoughts / options:

  • Even a test without assert is useful ... if it's in place, travis-ci makes sure that the code runs without raising an exception with various Python / Numpy / Astropy versions.
  • You could assert on the sum of the returned image. Then you have a regression test that makes sure the result is the same when you re-factor the code (e.g. introduce a RadialAcceptanceCurve class that stores the offset and acceptance arrays and has methods to read / write it to a FITS file).
  • In addition to the sum, you could assert on the values of two or three pixels, or on the maximum and minimum pixel position and value.
  • You could use photutils.detection.find_peaks with subpixel=True to get the maximum pixel position in the image with subpixel precision ... it should be very close to the input position and this can catch subtle off-by-one-pixel or off-by-0.5-pixel bugs.
  • You could try to break your function e.g. by passing an image that is larger than the acceptance curve, or a position completely outside your image, and check that your function behaves reasonably.

For now, my recommendation to you would be that you don't implement a ton of tests. This function is most likely not the final API we want (e.g. I think it makes sense to have a RadialAcceptanceCurve class), and you'll have to re-factor your tests if you re-factor your functions. I'm willing to merge this if you have a single assert on the sum of the image, doing more would be nice, but is not required at this point.

One more thing: please try to keep test datasets small. Think how gib the test image / curve needs to be so that you can test everything you want. Personally here I'd have chosen 5 x 5 or 10 x 10 pixel (although probably 3 x 3 would be fine as well), and not 100 x 100 or some big nice image. The advantage is that for small datasets the tests run faster, and if one gets an error, the printed arrays are easier to read / understand.

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 29, 2015

I forgot to mention one very good option to implement tests: check that the output matches some other code. E.g. here you could check against the HESS software. But I'd recommend against it because it takes too much time to set up a small test case in such a way that you can feed it into the HESS software ... which is one reason I started / like Gammapy ... it's simpler to use and more flexible.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Apr 1, 2015

@cdeil:
I improved the method by avoiding the loop over the pixels. Now it's very fast even for 100x100 pixel images.
The implemented test fails: I committed it for the purpose of asking for assistance. I have a test that tries to compare the integral of the radial acceptance curve (defined as a Gaussian), to the sum of the image, weighted by the pixel areas, which is in practice a discretization of the integral. Im principle, both results should be equivalent (within a certain small epsilon). But in practice, the difference is quite large and I'm not sure where does it come from.
A debug version of the test is posted in:
https://github.com/mapazarr/astropy_scripts/blob/master/astropy_scripts/test_fill_acceptance_image.py
BTW: I am still using 100x100 pixel images for the test because:

  • I need to correctly sample the radial acceptance (Gaussian) curve in order to have precise results to compare the integral and the sum of the image weighted by the pixel areas.
  • With the new method it is fast enough, so it doesn't take too much time to run the test.

Any input is welcome!

@mapazarr

This comment has been minimized.

Member

mapazarr commented Apr 1, 2015

On another note: could you tell me which is the convention to adopt for the origin in images? Is the origin pixel "1" as for FITS or "0"?
I'll use "0" as for photutils, and most of the code in gammapy (see issue #251).

ypix_coord_grid[x, 0:nx] = np.arange(0, ny)
#calculate pixel coordinates (in world coordinates)
coord = pixel_to_skycoord(xpix_coord_grid, ypix_coord_grid, w, 1)

This comment has been minimized.

@mapazarr

mapazarr Apr 1, 2015

Member

Here for instance, in the call to pixel_to_skycoord, I need to explicitly say if the origin is 0 or 1. I used 1 in this case.

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Changed to "0".

print(ypix_coord_grid) #debug
#calculate pixel offset from center (in world coordinates)
coord = pixel_to_skycoord(xpix_coord_grid, ypix_coord_grid, w, 1)

This comment has been minimized.

@mapazarr

mapazarr Apr 1, 2015

Member

Here also (pixel_to_skycoord).

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Changed to "0".

Empty image to fill.
center : `~astropy.coordinates.SkyCoord`
Coordinate of the center of the image.
offset : `~numpy.ndarray` of `~astropy.coordinates.Angle`

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

Simply write:

offset : `~astropy.coordinates.Angle`

(I think it's clear that it should be a 1D array, not a scalar or higher-dimensional object)

This comment has been minimized.

@mapazarr

mapazarr Apr 4, 2015

Member

Fixed

center : `~astropy.coordinates.SkyCoord`
Coordinate of the center of the image.
offset : `~numpy.ndarray` of `~astropy.coordinates.Angle`
Array of offset values in degrees where acceptance is defined.

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

No need to require "in degrees".
Users should be allowed to pass in Angle in any units they like (usually radian or deg).
If you really need some units internally, e.g. deg, you can convert to those units.

This comment has been minimized.

@mapazarr

mapazarr Apr 4, 2015

Member

Fixed

image : `~astropy.io.fits.ImageHDU`
Image filled with radial acceptance.
"""
#initialize WCS to the header of the image

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

There should be a space after #, i.e. the comment is

# initialize WCS

Discussing code formatting during code review is tedious for both parties.
That's why the Python programmers agreed on a style guide in PEP 8 and developed tools (pep8 and autopep8) to check and auto-format their code.

Please do this for all the code you write (after a while it'll be correctly formatted without you thinking about it).
But avoid mixing formatting fixes for the whole file with your code changes.
I do this by selecting the function I wrote in Pycharm and selecting "format code".
I'm not sure which tool it runs in the background, but it does a good job.

This comment has been minimized.

@mapazarr

mapazarr Apr 4, 2015

Member

Fixed.

pix_acc = f(pix_off)
#fill value in image
image.data = pix_acc

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

I would have put += here (and adapt the docstring to say this) to make it a bit simpler to fill many curves into one image. But it's not important, if you prefer the current version, fine.

This comment has been minimized.

@mapazarr

mapazarr Apr 4, 2015

Member

@cdeil
I'd rather leave it the way it is, or make a flag to switch between either overwriting the image content or add the acceptance to its content as you suggest. Which option would you prefer?

This comment has been minimized.

@cdeil

cdeil Apr 4, 2015

Member

Leave as is, then.
As I said this code will be wrapped or re-factored anyways when you introduce an AcceptanceCurve class, no need to add bells and whistles now.

@@ -1,2 +1,110 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import print_function, division
from gammapy.background import fill_acceptance_image

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

Please always sort imports like this:

  • Python stdlib
  • Numpy
  • scipy
  • Astropy
  • gammapy

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

Also for Gammapy, please use relative imports like everywhere else.

This comment has been minimized.

@mapazarr

mapazarr Apr 4, 2015

Member

Fixed order of imports.
@cdeil
What do you mean by relative import? Use:
from .gammapy.image import make_empty_image
instead of:
from gammapy.image import make_empty_image
for instance?
Should I use it only for gammapy or also for astropy?

This comment has been minimized.

@cdeil

cdeil Apr 4, 2015

Member

Browse Gammapy code base a bit and you'll quickly see the pattern.

Here's an example of a relative import in gammapy (https://github.com/gammapy/gammapy/blob/master/gammapy/background/tests/test_models.py#L8):

from ...background.models import GaussianBand2D

To import something from astropy in gammapy you have to use an absolute import.
Try doing a relative import for something in Astropy from gammapy ... you'll get an error saying that relative imports of things outside a given package (gammapy in this case) is not allowed.

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Ok, fixed.

from astropy.units.quantity import Quantity
import numpy as np
from scipy.integrate import quad

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

scipy is an optional dependency:
http://astropy.readthedocs.org/en/latest/development/testguide.html#tests-requiring-optional-dependencies
(search for scipy in test files if you want an example)

If you don't do this, test collection will fail for users that don't have scipy, which is not what we want.
We can discuss if scipy should become a required dependency for Gammapy in a separate issue, but for now, please fix this (or travis-ci will fail and the PR can't be merged)

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Fixed.

import numpy as np
from scipy.integrate import quad
def radial_gaussian_2D(r, sig):

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

I think this function is superfluous and should be removed.
Either call astropy.modeling.models.Gaussian1D(offset) or astropy.modeling.models.Gaussian2D(offset, 0) (I didn't look in detail what you implement here or need).

If you absolutely want a model that takes Angle objects as input and the Astropy models don't work the way you want for some reason, OK to write a wrapper, but internally please call the Astropy model and don't re-implement the formula. (But I think it's not needed, especially if you only use this in a test.)

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Fixed: using astropy.modeling.models.Gaussian1D, thanks for the heads up!

def test_fill_acceptance_image():
#create empty image

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

Comments should only be used where they add extra value and things can't be expressed clearly in code.
So remove this #create empty image comment because it doesn't add extra info to the image = make_empty_image() code.

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Fixed.

y = image.header['CRVAL2']
#sanity checks: am I using galactic coordinates in degrees?
assert image.header['CTYPE1'] == 'GLON-CAR', "Error: x coordinate is not in Galactic coordinates!"

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

These asserts don't belong here. If you want to add tests for make_empty_image, they should be added to gammapy/image/tests/test_utils.py in the test_make_empty_image function, but my suggestion would be to simply delete them in this PR. Here you should assume that make_empty_image works OK.

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Fixed: removed asserts.
@cdeil the idea was to have them in case the defaults of make_empty_image would change in the future, in which case it would be easier to find the reason for the test failure.

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Instead of the asserts, I now give the parameters to make_empty_image explicitly.

#fill acceptance in the image
image = fill_acceptance_image(image, center, offset, acceptance)
#test: sum of the image (weigthed by the pixel areas) should be equal to the integral of the radial acceptance

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

Wrap line to respect 80 character line lenght limit.

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Fixed.

#initialize WCS to the header of the image
w = WCS(image.header)
#define grids of pixel coorinates

This comment has been minimized.

@cdeil

This comment has been minimized.

@mapazarr

mapazarr Apr 5, 2015

Member

Yes! Thanks for the heads up. Fixed.

#integrate acceptance (i.e. gaussian function)
#1st integrate in r
acc_int = Quantity(quad(radial_gaussian_2D, Angle(0.*u.degree).to(u.radian).value, Angle(10.*u.degree).to(u.radian).value, args=(sigma.to(u.radian).value))*u.radian)

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

I can't read this line ... it's to long.
Use "temp variables" to split the expression over multiple lines.
(really there are no variables in Python, symbols are just references to objects on the stack, i.e. there's no computational overhead to use "temp variables" in Python. if this doesn't make sense and you don't have a Python book I can try to find a good article on the web that explains it).

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Fixed.

#check sum of the image:
# int ~= sum (bin content * bin size)
epsilon = 1.e-4
assert abs(image_int.to(u.rad**2).value - acc_int_value.to(u.rad**2).value) < epsilon, "image integral not compatible with radial acceptance integral"

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

Same here ... line is too long to read (on Github, but also in my editor).

Also, for such asserts please use np.testing.assert_almost_equal like in most other places in Gammapy.

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Fixed.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 1, 2015

I've pointed out some trivial stuff to improve in inline comments.
About your two core questions:

  1. why is the integral not correct? -> I have to sleep ... I'll try to find time tomorrow morning to have another look.
  2. pixel convention with offset 1 or 0 -> I'd prefer to explain via Skype than by typing lots of text and examples here ... let's Skype briefly sometime next week, email me which day would work for you.
from astropy.wcs import WCS
from astropy.wcs.utils import pixel_to_skycoord
import numpy as np
from scipy.interpolate import interp1d

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

This will fail on travis-ci. scipy is an optional dependency, so you should import it in the function where you use it.
The rule in Astropy is that users with missing optional dependencies should be able to import any Astropy / Gammapy module, and then only get the ImportError for the optional dependency when the call the function that needs it.

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Fixed.

#initialize WCS to the header of the image
w = WCS(image.header)
#define grids of pixel coorinates

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

There's a utility function for that. To generate 2d arrays of "pixel coordinates", there's a few numpy functions like np.indices or np.meshgrid that you can use. Look at how coordinates in gammapy.image is implemented.

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

But I think here you want to have 2d arrays of sky lon / lat arrays, i.e. just call gammapy.image.coordinates.

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Yes! Thanks for the heads up!

#interpolate
f = interp1d(offset, acceptance, kind='cubic')
pix_acc = f(pix_off)

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

I don't know if f can take 2d arrays as input here.
If it can't, you have to flatten the pix_off array here and then reshape it back to 2d (or use some other interpolation function from scipy).
This is the meat of your function and I'd have to look at the scipy docs / play with this a bit in an Ipython notebook myself to get it right.

#define radial acceptance and offset angles
offset = Angle(np.arange(0., 30., 0.1), unit=u.degree)
acceptance = np.zeros_like(offset)
sigma = Angle(1.0, unit=u.degree) #gaussian width

This comment has been minimized.

@cdeil

cdeil Apr 1, 2015

Member

I usually write

a = Angle(42, 'degree')

instead of

a = Angle(42 * u.degree)

If you don't mind, could you please also change to write Quantities by putting the units as strings?

There's two small advantages:

  1. I find "degree" more readable than "u.degree"
  2. Pycharm and other Python IDEs do static code analysis. But all the symbols in astropy.units are dynamically generated on import, so I get a warning displayed for a missing symbol "u.degree", which is distracting, because really all is OK.

This comment has been minimized.

@mapazarr

mapazarr Apr 8, 2015

Member

Fixed.

@cdeil cdeil changed the title from Added method for filling image with radial acceptance in FoV bg. to Add function to fill acceptance image from curve Apr 8, 2015

@cdeil cdeil added the feature label Apr 8, 2015

@cdeil cdeil added this to the 0.2 milestone Apr 8, 2015

@mapazarr

This comment has been minimized.

Member

mapazarr commented Apr 8, 2015

The test is now working. I followed your suggestion and I test on a few values, where the acceptance is defined. I had some trouble figuring out how the binning of the images work but now it's fine.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 9, 2015

Merging this now. Thanks!

cdeil added a commit that referenced this pull request Apr 9, 2015

Merge pull request #248 from mapazarr/fill_acc_image
Add function to fill acceptance image from curve

@cdeil cdeil merged commit 5d8447e into gammapy:master Apr 9, 2015

1 check passed

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

This comment has been minimized.

Member

cdeil commented Apr 9, 2015

The fill_acceptance_image function is now listed in the online docs. I hadn't checked the html formatting locally for this PR ... looks good!

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