Skip to content
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 energy dispersion RMF integration #731

Merged
merged 3 commits into from Nov 8, 2016
Merged

Add energy dispersion RMF integration #731

merged 3 commits into from Nov 8, 2016

Conversation

@JouvinLea
Copy link

@JouvinLea JouvinLea commented Sep 29, 2016

Add energy dispersion RMF integration

@joleroi
Copy link
Contributor

@joleroi joleroi commented Sep 30, 2016

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Oct 5, 2016

@joleroi
Hi now it's working for HAP the rmf even with 150bin in e_reco.
The test doesn't pass on Travis but I don't really understand the issue in test_hess_chains

@cdeil cdeil added the feature label Oct 5, 2016
@cdeil cdeil added this to the 0.5 milestone Oct 5, 2016
@cdeil cdeil changed the title RMF integral Add energy dispersion RMF integration Oct 5, 2016
@cdeil
Copy link
Member

@cdeil cdeil commented Oct 5, 2016

@JouvinLea - Thanks!

Can you please add a short description what this PR does in the first comment?
It currently only says "START PR".

And could you please squash the commits and put a good commit message?
These aren't very helpful: https://github.com/gammapy/gammapy/pull/731/commits
Maybe "Add energy dispersion RMF integration"?
That's what I put in the PR title for now, and we use the PR title for the changelog entries.

The implementation of the integration method is complex code that's difficult to verify for correctness via code review. Is it possible to extract the code in the for loop into a helper function or class? And add a docstring describing with one sentence what the integration method is. Is it trapezoidal integration?

@joleroi - I don't know if you have time this week, for now I've assigned this for you to review. If you don't have time, please assign back to me.

@joleroi
Copy link
Contributor

@joleroi joleroi commented Oct 6, 2016

I would propose the following API:

def integrate(self, offset, e_true, e_reco_min, e_reco_max):
   ""Integrate Energy dispersion over a given range in reco energy

    Description of the method
    """
    ...
    return integral

def get_response(self, offset, e_true, e_reco)
   """Detector response for a given true energy
    """"
    ....
    return response

Note that trapezoidal integration is implemented e.g. here

https://github.com/gammapy/gammapy/blob/master/gammapy/spectrum/utils.py#L222

I don't know if it's possible to reuse that code introduced by @adonath

@cdeil
Copy link
Member

@cdeil cdeil commented Oct 6, 2016

That code is specifically for log(energy)-log(flux) integration.

If I understand correctly, here we just have "normal" x-y integration of a PDF.

For that you should call
http://docs.scipy.org/doc/numpy/reference/generated/numpy.trapz.html
or if you think it's useful one of the methods from
http://docs.scipy.org/doc/scipy/reference/integrate.html

Although -- not that EDISP can be very noisy, and "fancy" adaptive integration methods can look fine for most analyses and then give crazy results or subtly incorrect results for some noisy IRFs we'll encounter later. So maybe just call the robust trapezoidal and assume the distributions are well-sampled?

@JouvinLea JouvinLea force-pushed the JouvinLea:rmf branch from a18a878 to 42a55ab Oct 7, 2016
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Oct 7, 2016

Hi,
I don't see the point to use a trapez integration whereas the edisp is only defined at the center of the migra bin this is why I use just a simple riemann sum

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Oct 10, 2016

@cdeil @joleroi
Do you agree with my previous comments? that since I get only the value of the edisp at the center of the migra bin there is no point to do a trapeze integration?

@cdeil
Copy link
Member

@cdeil cdeil commented Oct 11, 2016

My suggestion here would be that we use
http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html
which does interpolation and has an integrate method. And code-wise, it would probably be a good idea to split that 1-dimensional density distribution out into a separate class, in a similar way that we did for the PSF:
http://docs.gammapy.org/en/latest/api/gammapy.irf.TablePSF.html
There could still be a method on the existing class, but it would just dispatch to the method on the new class that represents the 1-dim PDF distribution (which could be cached or not, performance is not a bit concern here I think).

I expect that this will give a more flexible (in terms of interpolation / integration methods) and shorter implementation (compared to what's done here, rectangle integration and interpolation by hand), and it's more consistent wrt. what we do for the case of the PSF.

But it's quite a bit of work, maybe a day for a version that does what's done here, but with cleaner code and with tests added, and I don't have that time at the moment.

@joleroi
Copy link
Contributor

@joleroi joleroi commented Oct 12, 2016

I don't see the point to use a trapez integration whereas the edisp is only defined at the center of the migra bin this is why I use just a simple riemann sum

I though this gives you problems when you oversample the EDISP distrubution, e.g. when chosing the fine reco energy binning of HAP-HD

And code-wise, it would probably be a good idea to split that 1-dimensional density distribution out into a separate class

👍

I expect that this will give a more flexible (in terms of interpolation / integration methods) and shorter implementation (compared to what's done here, rectangle integration and interpolation by hand), and it's more consistent wrt. what we do for the case of the PSF.

The interpolation has been filtered out of the IRF classes in #625, and we agreed that this should also be done for the PSF. 👎 for adding the same 100 lines of boilter plate interpolation code to every class.

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Oct 26, 2016

@cdeil @joleroi @registerrier
Hi,
I don't know where we are with this PR but I really want this new calculation of the rmf to be included.
Do you agree on the fact that there is no interest for a trapezoidal interpolation here? Do you absolutely want to use the same kind of integration than what it is done for the PSF with the UnivariateSpline?

@joleroi joleroi mentioned this pull request Nov 8, 2016
@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 8, 2016

@JouvinLea Can you add a description in the docstring of get_response please? It still has the old description of what I have coded. It doesnt have to be long, just so we now what's there (I would do it myself, but I don't know what's there)

@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 8, 2016

It was agreed upon that we merge this now, and do the refactoring discussed above later.

@joleroi joleroi merged commit 824e351 into gammapy:master Nov 8, 2016
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@joleroi joleroi mentioned this pull request Nov 17, 2016
0 of 2 tasks complete
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants