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 cube background model class #299

Merged
merged 41 commits into from Jul 21, 2015

Conversation

Projects
None yet
2 participants
@mapazarr
Member

mapazarr commented Jul 11, 2015

This PR is a continuation of PR #292.
Task list:

  • the write method (and a test for it)
  • porting the test file to ~gammapy-extra
  • make the test that use external files work
  • implement read_image and write_image with a WCS object and a HDU table for the energy bin edges; tests
  • separate the plot functions into do_one_plot and do_many_plots functions
  • clean the code and autopep it
  • the sphinx doc; add link to examples/plot_background_model.py; add examples of many plots (mosaic of images, stack of spectra)
  • code review and final revision of inline comments

Further methods like histograming of event lists or interpolation and smoothing of the model will come on a later PR.

Sorry for the mess, I renamed the branch, which automatically closes the pull request. I thought the PR would be transparent to the rename.

@mapazarr mapazarr referenced this pull request Jul 11, 2015

Closed

Cube bg model container class: CubeBackgroundModel #292

5 of 8 tasks complete
@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 11, 2015

#292 (comment)
Ok, I moved the WCS functions to gammapy.utils.wcs. I named it wcs instead of the suggested wcsutils because being it already in utils makes the name redundant. In addition, none of the modules in utils are called somethingutils.

@mapazarr mapazarr force-pushed the mapazarr:cube-bg-container-class branch 2 times, most recently from 99431d4 to da60ebd Jul 11, 2015

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
You're right: the data members are already Quantity (or Angle) objects. But here's the thing.
This is the trick for storing the data in the FITS table: I need an array of an array, so that the table has only one line. I can either use Quantity, or np.array. The latter won't store the units, so I'd like to keep the Quantity, I think it's the most elegant and simple solution in order to store the data in the desired format, with the units that are already in the data members. Otherwise I'll need the lines that you wisely detected as not needed here:
#292 (comment)

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
Yes, the bg models (in FITS files) from @mimayer do have a E_THRES header keyword defined. In his scripts, he defines it as the energy of the lowest boundary of the 1st energy bin (as I defined them). I followed his workflow there. I don't think an extra data member is needed for it, but I don't have a strong preference: if you still want to have it, please let me know.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
Ok, the update method worked! Thanks for the tip.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
OK, but it might get confusing if we call both formats (astropy table and HDU bin table) the same (table).

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
I have to strip the units and put them back on later: np.log10 won't accept a quantity with units.
I leave it as it is, and put a TODO, that when PR #290 is done, it should be reviewed.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
#292 (comment)

Ok, we discussed this issue already here:

#292 (comment)
#292 (comment)

and I already simplified the methods. They are still long because I add 2 checks for a valid entry from the user. The method itself to find the bin is 4 lines on each case (X, Y and E). I'm not sure I can simplify it more. (We could drop the checks, if you think they don't belong there.)

In the case of the det bin (X, Y), I don't like using the wcs object since it's not a data member (we already discusse this: it makes more sense that the bins are data members instead of the wcs object).

In the case of the energy, I'm not sure what PR #290 does, but one can rewrite it once it is finished, if it makes sense. As I said, finding the bin is only 4 lines: actually it's 2 lines, plus extra 2 lines to find the edges (which I find nice/handy to have).

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
OK, I removed the print statement from here (the message can be added to the exception, as you say).

You say it is no good practice to add print statements to the library code. I think sometimes it can be useful to have some printouts, especially if the code takes long to run (which is not the case), but also to get some important information. An example would be:
https://github.com/gammapy/gammapy/pull/292/files#diff-82cdfbacfcc0119d2ccf6621da08de7bR719
The user requests the plot for a particular energy, and it gets to know the boundaries of the slice used for the plot. This is maybe not the most useful case, but serves as an example for my question: how would you issue such an information?

I removed in any case all print statements of the class.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

About the plots:

#292 (comment)
I do need it for the color axis. So it stays.

#292 (comment)
#292 (comment)
I removed the save options.

#292 (comment)
OK, I was using for instance the image object to get the data in the tests of the plotting functions. But I learnt how to get them from the ax object, so no need for them. Now the functions return only the ax object.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

#292 (comment)
I know. Ok, I adopted your naming, though I find more explicative the old denomination of the functions.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 12, 2015

OK I just pushed a commit addressing all open comments in the old PR #292. I commented above on some of them (the non-trivial ones). Please have a look at my replies, and if applicable please answer to them on a new comment.

For future comments, please use this PR (#299).

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 14, 2015

@mapazarr – I'm short on time this week, so I'll just do a final review tomorrow during the hangout.

Please update the task list and ask if you have any specific questions before you can finish this up.

Concerning the E_THRES header keyword ... if it's set to the energy array lower edge it's duplicated info and useless. But maybe in the future we do want to set it to indicate a valid range, or Gammalib requires it, so OK to keep if you like.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 15, 2015

Ok, I'm pretty much done with this PR. I think I have still one or 2 minor TODO's for you to review/comment on, on top of the questions I posted acouple of days ago above, and of course the final code review.

The link to the docs is here:
https://googledrive.com/host/0BzXHZQx0oCLBfmcyUDdxMU9sLVduMXM1QXk5RTNYU3dXYmthVG9OZ18wWE1uR3lYaEsxRU0/html_doc/docs/_build/html

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 15, 2015

This example doesn't work for me:

$ python examples/plot_background_model.py 
Traceback (most recent call last):
  File "examples/plot_background_model.py", line 23, in <module>
    write_kwargs=dict(clobber=True))
  File "/Users/deil/code/gammapy/gammapy/background/models.py", line 475, in write
    self.to_fits_image().writeto(outfile, **write_kwargs)
  File "/Users/deil/code/gammapy/gammapy/background/models.py", line 432, in to_fits_image
    imhdu.header.update(wcs_header)
  File "/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev12994-py3.4-macosx-10.10-x86_64.egg/astropy/io/fits/header.py", line 1167, in update
    update_from_dict(k, other[k])
  File "/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev12994-py3.4-macosx-10.10-x86_64.egg/astropy/io/fits/header.py", line 1149, in update_from_dict
    card = Card(k, v)
  File "/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev12994-py3.4-macosx-10.10-x86_64.egg/astropy/io/fits/card.py", line 448, in __init__
    self.value = value
  File "/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev12994-py3.4-macosx-10.10-x86_64.egg/astropy/io/fits/card.py", line 570, in value
    raise ValueError('Illegal value: %r.' % value)
ValueError: Illegal value: .

The problem are the '' and 'COMMENT' entries:

In [1]: dict(wcs_header)
Out[1]: 
{'': ,
 'CDELT1': 0.546914274042,
 'CDELT2': 0.546914274042,
 'COMMENT':  WCS header keyrecords produced by WCSLIB 5.8,
 'CRPIX1': 0.5,
 'CRPIX2': 0.5,
 'CRVAL1': -6.01605701447,
 'CRVAL2': -6.01605701447,
 'CTYPE1': 'DETX',
 'CTYPE2': 'DETY',
 'CUNIT1': 'deg',
 'CUNIT2': 'deg',
 'LATPOLE': 90.0,
 'WCSAXES': 2}

Merging the headers like this works for me, so that's what I'd suggest as a workaround:

header = wcs_header.update(imhdu.header)

The reason this error doesn't show up for you locally and on travis-ci is probably that you're not using the very latest Astropy master.
I'll make an issue that travis-ci docs builds with Python 3 and with the latest astropy master should be added.

@cdeil cdeil referenced this pull request Jul 15, 2015

Closed

Improve travis-ci setup / build matrix #300

4 of 4 tasks complete
@@ -27,11 +27,53 @@ TODO
Hello World.
.. _bg_models:
Background Models

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

I'd like the top-level gammapy.background page to remain short, just giving an overview of the available functionality and maybe one or two quick code examples / plots.

So please move this whole section to a sub page.
An example how to do this is here: https://gammapy.readthedocs.org/en/latest/image/plotting.html , i.e. plotting is a subpage of the image docs.

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

PS: the docs you put online on Google drive don't work (I get 404 page not found errors when I click on gammapy.background).
Contrary to what I said before: please don't put docs online any more! I'll simply check out your branch, built the docs locally and make inline comments on what I think should be fixed, using screenshots where needed. This works nicely.

This comment has been minimized.

@mapazarr
The file comes originally from the `~GammaLib` repository but has
been slightly modified.
TODO: the scripts produce some white canvases that I need to remove!

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Resolve this TODO.

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

It was resolved already. I forgot to remove the TODO.
Done.

# read
filename = '../test_datasets/background/bg_cube_model_test.fits'

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Remove this line.
Also the #read, #plot, #write comment lines can be removed IMO ... you've implemented a beautiful, clean API and there's no need to repeat in comments what the code is clearly expressing. :-)

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

Fixed.

:include-source:
More complex plots can be easily produced with a few lines of code:
TODO: mosaic/stack plots examples!!!

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Resolve this TODO.

I don't see the need for three plotting scripts, my suggestion would be to:

  • Remove the plotting from the first code example ... just have that read and write.
  • Combine the other two plotting scripts in one. Show three panels in one row (using plt.subplots), which show two images (one at very low and one at very high energy, so that they look different), and for the spectrum loop over a handful positions and draw the curves in the third panel.

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

I currently get this:

image

Why isn't the background rate a power-law, with lower rate at higher energy?

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

The TODO was resolved already. I forgot to remove it.
Done.

Changed the plot scheme to have only 1 plot with 3 pads.

The plots are correct. I still need to change the test file to have a more meaningful example.

@cdeil cdeil referenced this pull request Jul 15, 2015

Merged

Simplified attribute docstrings #301

2 of 2 tasks complete
table containing the bg cube
"""
# data arrays
a_detx_lo = Quantity([self.detx_bins[:-1]])

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Sorry if this has been discussed already, I don't remember: Why call Quantity here? self.detx_bins is already a quantity, so this is useless, no?

This comment has been minimized.

@mapazarr

mapazarr Jul 15, 2015

Member

As it is, it is the simplest way.
See this: #299 (comment)

This comment has been minimized.

@mapazarr

mapazarr Jul 15, 2015

Member

I can add an explanation in the code as a comment if you like.

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

So you basically want to reshape the array here?
Then please call reshape or maybe numpy.expand_dims to make the code clearer.
(I hope it works with Quantity objects).

Also, in this case the temp variables aren't really needed, you can just assign to the table columns below and it'll still nicely fit on one line, no?

table.meta['E_THRES'] = a_energy_lo.flatten()[0].value
# name

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Remove this comment.

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

Fixed.

tbhdu : `~astropy.io.fits.BinTableHDU`
table containing the bg cube
"""
# build astropy table and convert it to fits bin table

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Remove this comment ... it says exactly the same thing as the code on the next line.

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

Fixed.

enhdu = _table_to_fits_bin_table(energy_table)
hdu_list = fits.HDUList([fits.PrimaryHDU(), imhdu, enhdu])

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

Why put an empty PrimaryHDU?
Please make a PrimaryHDU directly instead of an ImageHDU, and then there's no need for an empty PrimaryHDU, right?

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

The issue with the primary HDU is that it's not possible to give it a name (in this case 'BACKGROUND'), which makes it a bit less easy to identify. Still, in order to have an image HDU on a fits file, a primary HDU is necessary.

I have seen other FITS images which use an empty PrimaryHDU, so I followed the trend.

I also was trying to follow the same pattern as in SpectralCube, which uses an image HDU for the cube, and a table HDU for the energy binning:
http://gammapy.readthedocs.org/en/latest/_modules/gammapy/data/spectral_cube.html#SpectralCube.to_fits
This BTW won't work in the writeto function: in order to store an image HDU a primary HDU has to be set as well. Since there is no test covering the writeto method, this problem passed unnoticed.

I don't have a strong opinion on this, but if you also don't, I'd like to keep it as it is now. Please let me know if you still want me to change it.

Once we agree on one convention I can apply the same to the SpectralCube, and do a simple test to check that the write method does not crash.
Reported in issue #304.

This comment has been minimized.

@cdeil

cdeil Jul 16, 2015

Member

As I said, I prefer to have the cube in the primary HDU.
Please change this here for the background cube and I think just fixing the SpectralCube write method in a commit here is better, no need for a separate PR.

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

OK. Done. It was also a bit of work to make the SpectralCube.writeto method to work.
Please see the TODO that I left inside.
This solves issue #304 .

return hdu_list
def write(self, outfile, format='table', write_kwargs=None):

This comment has been minimized.

@cdeil

cdeil Jul 15, 2015

Member

I'd prefer you use **kwargs here instead of write_kwargs=None.
That's the common pattern and it should work.

I remember this has been discussed before, and there was some issue.
But please, let's try to make this work and re-discuss if needed.
One thing you have to know is that **kwargs always has to be at the end.

This comment has been minimized.

@mapazarr

mapazarr Jul 16, 2015

Member

Done, though I have seen the other pattern in other functions as well.

dety_edges_centers : `~astropy.coordinates.Angle`
array of image bin centers (Y coord)
"""
detx_edges_centers = self.detx_bins[:-1] + np.diff(self.detx_bins)/2.

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

Thinking about this some more, I think this is a tick more readable:

bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[-1:])

and is more similar to the expression for log bin centers (geometric mean):

bin_centers = np.sqrt(bin_edges[:-1] * bin_edges[1:])

So please use that.

It's also faster! 🚤 🐎
:-)

This comment has been minimized.

@mapazarr

mapazarr Jul 20, 2015

Member

Done. It is how I had done it before, only I had temp vars for bin_edges[:-1] and bin_edges[1:] for readability purposes...

detx_center = (detx_max + detx_min)/2.
dety_center = (dety_max + dety_min)/2.
mask = (detx_points > detx_center) & (dety_points > dety_center)
#import IPython; IPython.embed()

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

Remove these comments

This comment has been minimized.

@mapazarr
bg_cube_model = make_test_bg_cube_model(apply_mask=True)
# test that values with (x, y) > (0, 0) are zero
x_points = Angle(np.arange(5), 'degree') + Angle(0.01, 'degree')

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

Slightly shorter:

x_points = Angle(np.arange(5) + 0.01, 'degree')

This comment has been minimized.

@mapazarr
assert_allclose(actual, desired, *args, **kwargs)
def assert_angle(actual, desired, *args, **kwargs):

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

I don't think this is really useful in addition to the quantity assert helper function?

This comment has been minimized.

@mapazarr

mapazarr Jul 20, 2015

Member

Ok. I removed it. The other one should fall soon also :-)

Parameters
----------
name_x : `~string`

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

string -> str

This comment has been minimized.

@mapazarr
w.wcs.cunit = [unit_x, unit_y]
w.wcs.cdelt = [delta_x.to(unit_x).value, delta_y.to(unit_y).value]
# ref as lower left corner (start of (X, Y) bin coordinates)
# coordinate start empiricaly determined at pix = 0.5: why 0.5?

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

I think we discussed this and it's understood, so remove the "why 0.5?" comment?

This comment has been minimized.

@mapazarr

mapazarr Jul 20, 2015

Member

Sorry I forgot. Done.

raise ValueError("The Y bins are not linear! Diff = {}".format(np.diff(bin_edges_y.value)))
# Create a new WCS object. The number of axes must be set from the start
w = WCS(naxis=2)

This comment has been minimized.

@cdeil

cdeil Jul 20, 2015

Member

w -> wcs (single-letter variables are a bit evil)

This comment has been minimized.

@mapazarr

mapazarr Jul 20, 2015

Member

I have changed it already. Maybe you're looking at an old version.

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 20, 2015

Left some more minor inline comments.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 20, 2015

Ok I implemented the last comments. I also fixed some problems in SpectralCube nad the tutorials, that where preventing the sphinx job in travis-ci to work.

Can we merge this PR already? Some of the minor comments are turning around on issues discussed a few times already, and undoing some of the proposed changes.

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 21, 2015

Merging this now.
Thanks for fixing all the issues that came up in this PR or following up in other Astropy / Gammapy issues!

I looked through the coverage report

time python setup.py test --remote -V --coverage
open htmlcov/gammapy_background_models.html

and generally coverage is good.
The only thing that's not tested is the code under except ValueError in _parse_bg_units
Why is that needed at all? Maybe because the background cubes from Michael or Gammalib don't have correct units per the FITS standard? If so we should tell them to fix it and make a note that this unit handling helper function should be removed as soon as it's done.

cdeil added a commit that referenced this pull request Jul 21, 2015

Merge pull request #299 from mapazarr/cube-bg-container-class
Add cube background model class

@cdeil cdeil merged commit f3ed37a into gammapy:master Jul 21, 2015

1 check passed

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

This comment has been minimized.

Member

cdeil commented Jul 21, 2015

🎆

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 21, 2015

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 21, 2015

Hi Christoph. Indeed the existence of _parse_bg_units is because of the units in @mimayer's files. I checked and they actually comply the FITS standard, but they are not understood by ~astropy.units.Units.

Please see the text at the examples section of: https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/general/ogip_93_001/ogip_93_001.html

'the number of spaces between individual unit sub-strings is left a matter for personal preference'

@mimayer doesn't use any spaces, i.e. 1/s/MeV/sr, while Quantity expects spaces and can't parse this unit, i.e. 1 / s / MeV / sr.

According to the FITS standards both are valid.

This could be another issue for astropy.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 21, 2015

Thanks for the merging! :-)

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 21, 2015

I didn't get any coveralls report.

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 21, 2015

You have to manually go to https://coveralls.io/github/gammapy/gammapy (I always do it via the link here) and find the right build.
There is the possibility to configure Coveralls to post comments on PRs with links to the right build.
But so far we weren't using that much and it was annoying.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Jul 21, 2015

BTW, I think the test file from gammalib also was having the issue with the units.

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Aug 30, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Aug 30, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Aug 30, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Aug 30, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Sep 25, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Sep 25, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Sep 25, 2015

JonathanDHarris added a commit to JonathanDHarris/gammapy that referenced this pull request Sep 25, 2015

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