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

Spectrum energy grouping #709

Merged
merged 15 commits into from Sep 15, 2016

Conversation

Projects
None yet
2 participants
@cdeil
Member

cdeil commented Sep 13, 2016

This pull request

  • renames gammapy/spectrum/group.py to gammapy/spectrum/obs_group.py to make the distinction between grouping energy bins or observations clear.
  • adds ObservationStats.to_dict and SpectrumObservation.stats_table, which are useful to have to compute energy groupings.
  • add table_from_row_data utility function
  • add SpectrumStats class
  • introduces gammapy/spectrum/energy_group.py and moves the calculate_flux_point_binning function which was in gammapy/spectrum/utils.py there.
  • Adds a SpectrumEnergyGroupMaker class that takes a SpectrumObservation and has methods to compute energy groupings (energy bin range selection and grouping)
  • Under the hood, the SpectrumEnergyGroupMaker uses a SpectrumEnergyGroups object to do the bookkeeping and some of the calculations, which is a list of SpectrumEnergyGroup objects, which are simple containers to do book-keeping and some simple operations.
  • Adds some tests.

This is left to future PRs:

  • Changing DifferentialFluxPoints.compute and calculate_flux_point_binning to work with ENERGY_GROUP_ID index arrays instead of ENERGY_BOUNDS float arrays with the brittle float eps check will be done in a future pull request.
  • There's also many open TODOs and I'm not fully happy with the overall design.
  • Implement adaptive grouping
  • Add more tests (also degenerate cases that result in one or no groups)

@cdeil cdeil added the cleanup label Sep 13, 2016

@cdeil cdeil added this to the 0.5 milestone Sep 13, 2016

@cdeil cdeil self-assigned this Sep 13, 2016

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Sep 13, 2016

Member

@joleroi @registerrier @leajouvin - If you have time, could you please have a look at this and let me know if you have any thoughts on how energy bin grouping should be implemented and exposed in Gammapy?

(Alternatively, do you have examples how to use the energy bin grouping in Sherpa and do you think we should just use that or wrap it in Gammapy?)

Member

cdeil commented Sep 13, 2016

@joleroi @registerrier @leajouvin - If you have time, could you please have a look at this and let me know if you have any thoughts on how energy bin grouping should be implemented and exposed in Gammapy?

(Alternatively, do you have examples how to use the energy bin grouping in Sherpa and do you think we should just use that or wrap it in Gammapy?)

@cdeil cdeil changed the title from Add new module for spectrum energy grouping to Spectrum energy grouping Sep 13, 2016

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Sep 14, 2016

Member

@joleroi - I'll merge this tomorrow morning. But comments or even a full review would still be very welcome. I plan to continue improving and extending this in the coming days.

Member

cdeil commented Sep 14, 2016

@joleroi - I'll merge this tomorrow morning. But comments or even a full review would still be very welcome. I plan to continue improving and extending this in the coming days.

@cdeil cdeil merged commit a6cd2e2 into gammapy:master Sep 15, 2016

2 checks passed

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

This comment has been minimized.

Show comment
Hide comment
@joleroi

joleroi Sep 16, 2016

Contributor

I had a look at your notes. Here are my comments

  • Note that the current functionality was hacked together in a few hours in order to produce some flux points with gammapy. I tried to leave the code in isolated places (e.g. utils) to make refactoring it easier.
  • I agree that the flux point computation methods should only accept GROUP_ID vector and not fiddle with the energy binning
  • I like your new proposal. As input for the grouping class a stacked SpectrumObservations seems ideal for me. IMO the purpose if the SpectrumObservations.stats_table should be easier debug output not actual functionality. We left root TH1 behind, and it should stay like this
  • The examples you give seem straight forward 👍
Contributor

joleroi commented Sep 16, 2016

I had a look at your notes. Here are my comments

  • Note that the current functionality was hacked together in a few hours in order to produce some flux points with gammapy. I tried to leave the code in isolated places (e.g. utils) to make refactoring it easier.
  • I agree that the flux point computation methods should only accept GROUP_ID vector and not fiddle with the energy binning
  • I like your new proposal. As input for the grouping class a stacked SpectrumObservations seems ideal for me. IMO the purpose if the SpectrumObservations.stats_table should be easier debug output not actual functionality. We left root TH1 behind, and it should stay like this
  • The examples you give seem straight forward 👍
OVERFLOW_BIN_INDEX = -2
class SpectrumEnergyGroupMaker(object):

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

I think at the coding sprint we agreed to not use Makers in Gammapy.
Why not just SpectrumEnergyGrouping?

@joleroi

joleroi Sep 16, 2016

Contributor

I think at the coding sprint we agreed to not use Makers in Gammapy.
Why not just SpectrumEnergyGrouping?

Attributes
----------
obs : `~gammapy.spectrum.SpectrumObservation`
Spectrum observation data

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

Is this stacked or not? Looks like stacked, but it should be clear from the docstring

@joleroi

joleroi Sep 16, 2016

Contributor

Is this stacked or not? Looks like stacked, but it should be clear from the docstring

# Methods to compute groupings
def compute_groups_fixed(self, ebounds):

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

The stacked spectrum object does not have data outside the safe energy range. What does this method do if ebounds[0] < obs.lo_threshold?

@joleroi

joleroi Sep 16, 2016

Contributor

The stacked spectrum object does not have data outside the safe energy range. What does this method do if ebounds[0] < obs.lo_threshold?

def compute_groups_npoints(self, npoints):
"""TODO: document"""
emin, emax = self.get_safe_range()

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

I don't find the get_safe_range() method? Anyway, IMO it should be on SpectrumObservation.

@joleroi

joleroi Sep 16, 2016

Contributor

I don't find the get_safe_range() method? Anyway, IMO it should be on SpectrumObservation.

def compute_groups_npoints(self, npoints):
"""TODO: document"""
emin, emax = self.get_safe_range()
npoints = self.config['n_points']

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

What is the config attribute?

@joleroi

joleroi Sep 16, 2016

Contributor

What is the config attribute?

return ss
class EnergyRange(object):

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

I don't think this is the appropriate place for this class. It reproduces a lot of the functionality that is already in the EnergyBounds. See also
http://docs.gammapy.org/en/latest/utils/index.html#energy-handling-in-gammapy

@joleroi

joleroi Sep 16, 2016

Contributor

I don't think this is the appropriate place for this class. It reproduces a lot of the functionality that is already in the EnergyBounds. See also
http://docs.gammapy.org/en/latest/utils/index.html#energy-handling-in-gammapy

return str(self.to_dict())
class SpectrumEnergyGroups(list):

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

It's hard to see the purpose of this class without much documentation. At first sight it looks super-complicated. In the end we just want a list with integers to put each fine energy bin into one flux points group, this list then goes to the flux point computation method. Can you add some explanation what this is for?

@joleroi

joleroi Sep 16, 2016

Contributor

It's hard to see the purpose of this class without much documentation. At first sight it looks super-complicated. In the end we just want a list with integers to put each fine energy bin into one flux points group, this list then goes to the flux point computation method. Can you add some explanation what this is for?

'SpectrumObservation',
'SpectrumObservationList',
]
class SpectrumStats(ObservationStats):

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

Nice, I wanted to have this since a long time ago

@joleroi

joleroi Sep 16, 2016

Contributor

Nice, I wanted to have this since a long time ago

@@ -78,6 +99,11 @@ def alpha(self):
return self.on_vector.backscal / self.off_vector.backscal
@property
def ebounds(self):
"""Energy bounds array."""
return self.on_vector.energy.data.to('TeV')

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

Why would you have the ``to('TeV')`?

@joleroi

joleroi Sep 16, 2016

Contributor

Why would you have the ``to('TeV')`?

# in order to add spectrum specific information
kwargs = dict(
return SpectrumStats(
energy_min=self.ebounds[:-1],

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

This should be self.lo_threshold, no? That's at least what I would expect

@joleroi

joleroi Sep 16, 2016

Contributor

This should be self.lo_threshold, no? That's at least what I would expect

obs_id=self.obs_id,
livetime=self.livetime,
)
return ObservationStats(**kwargs)
def stats_table(self):

This comment has been minimized.

@joleroi

joleroi Sep 16, 2016

Contributor

👍

@joleroi

joleroi Sep 16, 2016

Contributor

👍

@joleroi

This comment has been minimized.

Show comment
Hide comment
@joleroi

joleroi Sep 20, 2016

Contributor

@cdeil Do you still plan to work on this?

Contributor

joleroi commented Sep 20, 2016

@cdeil Do you still plan to work on this?

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