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 EffectiveAreaTable exporter to EffectiveAreaTable2D #276

Merged
merged 4 commits into from Jul 15, 2015

Conversation

Projects
None yet
3 participants
@joleroi
Contributor

joleroi commented Jun 2, 2015

TODO:

  • add test
  • add example

I dont know why the entire EffArea2D class seems to be "new" probably I messed something up with git

@joleroi joleroi force-pushed the joleroi:pha_export branch 2 times, most recently from bda2c1e to 5345cd2 Jun 2, 2015

@joleroi joleroi changed the title from Pha export to PHA export Jun 2, 2015

@joleroi joleroi changed the title from PHA export to WIP: PHA export Jun 2, 2015

@coveralls

This comment has been minimized.

coveralls commented Jun 2, 2015

Coverage Status

Coverage remained the same at 43.12% when pulling bda2c1e on kingj90:pha_export into 0883871 on gammapy:master.

@coveralls

This comment has been minimized.

coveralls commented Jun 2, 2015

Coverage Status

Coverage decreased (-0.06%) to 43.06% when pulling 5345cd2 on kingj90:pha_export into 0883871 on gammapy:master.

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jun 2, 2015

the test doesn't really test much, just writes to disk. Don't know what to test in order to verify that everything is ok

@coveralls

This comment has been minimized.

coveralls commented Jun 2, 2015

Coverage Status

Coverage increased (+0.04%) to 43.17% when pulling 51b29cb on kingj90:pha_export into 0883871 on gammapy:master.

@@ -444,6 +446,67 @@ def read(filename):
hdu_list = fits.open(filename)
return EffectiveAreaTable2D.from_fits(hdu_list)
def to_effectiveareatable(self, offset, energ_lo = None, energ_hi = None):

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

to_effectiveareatable -> to_effective_area_table

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

Please run PEP8 on this file. There's some whitespace things that should be changed, e.g. it should be energ_lo=None without spaces instead of energy_lo = None with spaces for keyword arguments.

If the effective area table is intended to be used for spectral analysis,
the final energy binning should be given at this point, since the
effective area table class is not able to perform interpolation at
this point.

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

This text is a bit cryptic and it's easy to forget to change / remove this comment when the EffectiveAreaTable class is updated. Can you at least leave a TODO here or if you have a concrete suggestion how to improve things, leave a TODO at that place?

This comment has been minimized.

@joleroi

joleroi Jun 3, 2015

Contributor

I don't think its necessary to update the EffectiveAreaTable class (in this regard), if it is meant to be a wrapper around ARF files.

Returns
-------
eff_area_table : `EffectiveAreaTable`
Effective area table class.

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

Remove "class".
Python is very dynamic and you can return a class.
But here you return an instance, which is pretty much what we do all the time, so no need to add "object" or "instance" at the end of this line.

if not isinstance(energ_lo, Quantity) or not isinstance(energ_hi, Quantity):
raise ValueError("Energy must be a Quantity object.")
if len(energ_lo) != len(energ_hi):
raise ValueError("Energy Vectors must have same dimension")

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

dimension -> length

return EffectiveAreaTable(energ_lo, energ_hi, area)
def write_ARF(self, filename, offset, energ_lo = None, energ_hi = None):

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

write_ARF -> write_arf

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

Do you think it's worth adding this 2-line convenience method?
I'm leaning very slightly towards no, but if you like it, keep it!

With it you'd write

effective_area_2d.write_arf('arf.fits', Angle(2, 'deg'))

Without it

effective_area_2d.to_effective_area(Angle(2, 'deg')).write('arf.fits', format='arf'))

This comment has been minimized.

@joleroi

joleroi Jun 3, 2015

Contributor

I agree - I will remove it

-------
arf : `EnergyDependentARF`
ARF object.
"""

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

This docstring gives a Sphinx warning, which makes the travis-ci build fail (i.e. this PR can't be merged):
https://travis-ci.org/gammapy/gammapy/jobs/65093013#L1843

I'm not sure what the issue is. Maybe the missing space between filename and colon in filename:?

Debugging / fixing Sphinx warnings and errors is tedious (or I'm simply not aware of a good method).
What I do is usually remove things one by one until the error goes away and then when I know where the issue is I can usually guess or Google how to fix it.

This comment has been minimized.

@joleroi

joleroi Jun 3, 2015

Contributor

problem solved :)

@@ -130,3 +130,13 @@ def test_EffectiveAreaTable2D(method):
actual = effarea.evaluate(offset=offset, energy=energy).shape
desired = np.zeros([nbinso, nbinse]).shape
assert_equal(actual, desired)
#Test ARF export

This comment has been minimized.

@cdeil

cdeil Jun 2, 2015

Member

The test that writes ARF files should be for the EffectiveAreaTable class.
There you could write it, read it back into a second object and assert that some things in it's content are the same before and after writing.

Here you should just test to_effective_area. You could e.g. assert that evaluating the 1D and 2D effective area gives identical or close results if you use the same offset.

@joleroi joleroi force-pushed the joleroi:pha_export branch from 51b29cb to 4eb5d79 Jun 10, 2015

>>> aeff2D = EffectiveAreaTable2D.from_fits(load_aeff2D_fits_table())
>>> offset = Angle(0.43, 'degree')
>>> nbins = 50
>>> energy = Quantity(np.logspace(0, 1, nbins+1), 'TeV')

This comment has been minimized.

@cdeil

cdeil Jun 10, 2015

Member

Use the energy_bounds_equal_log_spacing utility function here?
(adding an Examples section there with one or two examples would be nice, I guess it's not obvious how to use it).

This comment has been minimized.

@joleroi

joleroi Jun 16, 2015

Contributor

Thanks for pointing this out. Why is this method take a Quantity tuple as input, this is indeed a little counter intuitive, I think. I would it expect to work like
energy_bound...._spacing(emin, emax, nbins)
Then you could also to stuff like
e.._spacing(Quantity(100,'keV'),Quantity(10,'TeV'), 200)

What do you think?

This comment has been minimized.

@cdeil

cdeil Jun 16, 2015

Member

In some places in Gammapy we currently use an energy_band tuple, in some separate energy_min and energy_max variables / parameters.

I'm 50:50 on this, but we probably should choose one and be more consistent.

@mapazarr, @adonath - What do you think?

This comment has been minimized.

@joleroi

joleroi Jun 16, 2015

Contributor

the same holds for
energy_axis vs. energ_lo + energ_hi
it would be nicer to be consistent there

My opinion: Get rid of the energy_band tuple, I guess Quantities have few advantages here could be one
Get rid of energ_lo and energ_hi, just makes more effort to create (I guess the reason we have this is because its like that in OGIP standard)

This comment has been minimized.

@cdeil

cdeil Jun 16, 2015

Member

@kingj90 – I don't understand what your preference is. You say "get rid of energy_band tuple" and "get rid of energy_lo and energy_hi", but as I understand it these are the two options?

I don't think Quantity or not is an issue ... in Gammapy energies should always be quantities and then when reading / writing to file we just have to make sure we store the unit correctly (or convert to an agreed unit if the storage format doesn't support storing the unit).

@kingj90 - it might make sense to just do here whatever you like best and split this discussion out into a separate issue / pull request. As done for times here we should have a "How to handle energies" section in the docs that also explains how to handle energy bounds, and the PR where this is written could then also adjust the API to be more consistent.

This comment has been minimized.

@joleroi

joleroi Jun 16, 2015

Contributor

ok sorry, I was not clear (but I agree this is not the place to discuss this)
Just for completeness: I was talking about two different things

  1. energy bound (nbins from emin to emax): used to generate an energy axis. I would prefer the API (emin, emax, nbins)
  2. energy axis (nparray holding the axis): when passed to a function (e.g. effective area table). I would make more sense to me to pass the axis with lenght (nbins+1) instead of passing two arrays (energ_lo, enerh_hi) with upper and lower bound. If we would pass only one axis this could be generated with the energy_bounds... function, otherwise there are unnecessary lines like energ_lo = energy_axis[:-1] and ener_hi = energy_axis[1:]

This comment has been minimized.

@cdeil

cdeil Jun 16, 2015

Member

@kingj90 Sounds good to me, please proceed like this here, or even clean up other code in Gammapy in this PR or a separate one to follow this style.

Although, I haven't thought about it much / looked at which parts or the code would change in which way, so I reserve the right to say something else in the future. One thing that might make sense in the future is an EnergyBinning class that is a wrapper or sub-class of the n+1 length quantity corresponding the the EBOUNDS extension, with a way to read / write this to a HDU and convenience attributes or functions to get e.g. the array of lower or upper energy bin edges. I only vaguely remember, but I think in the HESS software we have this as an EnergyAxis class (although I'm not sure "axis" is the best term for this).

@@ -130,22 +130,26 @@ def to_fits(self, header=None, **kwargs):
For more info on the ARF FITS file format see:
http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/summary/cal_gen_92_002_summary.html
Recommended units for ARF tables are keV and cm^2, but TeV and m^2 are chosen here
as the more natural units for IACTs
Recommended units for ARF tables are keV and cm^2, but TeV and m^2 are chosen here by default, as they are the more natural units for IACTs

This comment has been minimized.

@cdeil

cdeil Jun 10, 2015

Member

Wrap line.

@@ -111,7 +111,7 @@ def __init__(self, energy_lo, energy_hi, effective_area,
self.energy_thresh_lo = energy_thresh_lo.to('TeV')
self.energy_thresh_hi = energy_thresh_hi.to('TeV')
def to_fits(self, header=None, **kwargs):
def to_fits(self, header=None, units=['TeV','m2'], **kwargs):

This comment has been minimized.

@cdeil

cdeil Jun 10, 2015

Member

You should never use mutable objects like lists as default arguments in Python as explained here.
Also this units option is currently not listed in Parameters in the docstring.

I haven't thought about this and don't have a good concept for default units here and in other places in Gammapy.
There's many options to handle this, e.g.

  • Separate keyword arguments
  • No option to change units at this point

Do what you think is best and it's good enough for now.
Other than lists as default keyword arguments, that is :-)

Returns
-------
eff_area_table : `EffectiveAreaTable`
Effective area table instance.

This comment has been minimized.

@cdeil

cdeil Jun 10, 2015

Member

Remove instance ... it's always an instance = an object and there's no point in spelling that out everywhere or sometimes.

energ_hi = self.energ_hi
elif energ_lo is None or energ_hi is None:
raise ValueError("Only 1 energy vector given, need 2")
if not isinstance(offset, Angle):

This comment has been minimized.

@cdeil

cdeil Jun 10, 2015

Member

Move this offset input validation below all the energy input checks (and separate by a blank line)

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 10, 2015

Thanks for finishing this up. I've left some more inline comments. 😄

When you addressed those, please change the issue title to something like "Implement EffectiveAreaTable2D.to_effective_area_table" and add an entry to the changelog

Then this is ready to be merged.

@cdeil cdeil added the feature label Jun 10, 2015

@cdeil cdeil added this to the 0.3 milestone Jun 10, 2015

@joleroi joleroi changed the title from WIP: PHA export to Implement EffectiveAreaTable2D.to_effective_area_table Jun 16, 2015

@joleroi joleroi force-pushed the joleroi:pha_export branch from d6d4ab7 to e9b0ee7 Jun 16, 2015

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jun 17, 2015

Is this something I should resolve? If so, how :) ?

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 17, 2015

@kingj90 - Yes, I'd prefer you resolve the merge conflict yourself. One advantage is less work for me, another is that travis-ci will re-run the tests.

You have to rebase on master. How to do it a bit depends on how you work with the upstream master branch locally, i.e. what exactly your git workflow is.

The recommended procedure for Astropy is here

I use this:

git checkout master
git status # should be clean
git pull # get latest version from upstream
git checkout mybranch
git rebase master --interactive
# resolve conflict
git push -f kingj90 mybranch

You might want to make a copy of your gammapy folder before ... it's not unheard of that people git in a pickle on their first few rebases. :-)

@joleroi joleroi force-pushed the joleroi:pha_export branch from e9b0ee7 to 676ced9 Jun 17, 2015

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jun 24, 2015

anything else to do here?

----------
offset : `~astropy.coordinates.Angle`
offset
energy : `~astropy.units.Quantity`

This comment has been minimized.

@cdeil

cdeil Jun 24, 2015

Member

Parameter docs (energy) is inconsistent with actual parameters (energ_lo and energhi).

I don't like energ ... can you use energy in Gammapy for variable names even if some FITS format uses energ?

@@ -111,7 +111,7 @@ def __init__(self, energy_lo, energy_hi, effective_area,
self.energy_thresh_lo = energy_thresh_lo.to('TeV')
self.energy_thresh_hi = energy_thresh_hi.to('TeV')
def to_fits(self, header=None, **kwargs):
def to_fits(self, header=None, energy_u='TeV', effarea_u='m2', **kwargs):

This comment has been minimized.

@cdeil

cdeil Jun 24, 2015

Member

I'd prefer more explicit: unit instead of u in the parameter name.

@@ -130,22 +136,28 @@ def to_fits(self, header=None, **kwargs):
For more info on the ARF FITS file format see:
http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/summary/cal_gen_92_002_summary.html
Recommended units for ARF tables are keV and cm^2, but TeV and m^2 are chosen here
as the more natural units for IACTs
Recommended units for ARF tables are keV and cm^2,

This comment has been minimized.

@cdeil

cdeil Jun 24, 2015

Member

"Recommended units" -> "Recommended units in X-ray astronomy"

@@ -120,6 +120,12 @@ def to_fits(self, header=None, **kwargs):
header : `~astropy.io.fits.header.Header`
Header to be written in the fits file.
energy_u : str
Unit in which the energy is written in the fits file

This comment has been minimized.

@cdeil

cdeil Jun 24, 2015

Member

No need for empty lines between parameters in docstrings.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 24, 2015

Almost ready to be merged ... left one last round of inline comments.

@joleroi joleroi force-pushed the joleroi:pha_export branch from 676ced9 to 4357848 Jun 24, 2015

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 15, 2015

@kingj90 – I think you addressed my last comment in 4357848 on Jun 24.
Github doesn't send notifications when new commits are made, so I didn't notice.
You have to leave a comment saying "OK to merge" or something ...

By now, there's a merge conflict, so you have to rebase on master so that travis-ci runs the tests again and so that the green merge button appears for me to click.

@joleroi joleroi force-pushed the joleroi:pha_export branch from 4357848 to 9cf89b0 Jul 15, 2015

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jul 15, 2015

OK to merge ;)

@cdeil cdeil changed the title from Implement EffectiveAreaTable2D.to_effective_area_table to Add EffectiveAreaTable exporter to EffectiveAreaTable2D Jul 15, 2015

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

Merge pull request #276 from kingj90/pha_export
Add EffectiveAreaTable exporter to EffectiveAreaTable2D

@cdeil cdeil merged commit 650707d into gammapy:master Jul 15, 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 15, 2015

Thanks!

@joleroi joleroi deleted the joleroi:pha_export branch Jul 15, 2015

>>> aeff2D = EffectiveAreaTable2D.from_fits(load_aeff2D_fits_table())
>>> offset = Angle(0.43, 'degree')
>>> nbins = 50
>>> energy = energy_bounds_equal_log_spacing(Quantity((1,10), 'TeV', nbins)

This comment has been minimized.

@cdeil

cdeil Jul 16, 2015

Member

I just tried this example locally.
There's a missing closing parens on this line.

If I add it, I get this error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-55b7c69c7675> in <module>()
      8 offset = Angle(0.43, 'degree')
      9 nbins = 50
---> 10 energy = energy_bounds_equal_log_spacing(Quantity((1,10), 'TeV', nbins))
     11 energ_lo = energy[:-1]
     12 energ_hi = energy[1:]

/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev12994-py3.4-macosx-10.10-x86_64.egg/astropy/units/quantity.py in __new__(cls, value, unit, dtype, copy, order, subok, ndmin)
    247 
    248         value = np.array(value, dtype=dtype, copy=copy, order=order,
--> 249                          subok=False, ndmin=ndmin)
    250 
    251         # check that array contains numbers or long int objects

TypeError: data type not understood

@kingj90 – Can you fix this example as a separate commit in your next PR?

Also: code examples that don't show output don't have to have the >>> on every line, you can remove that.

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