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 support for observations with different energy binning in SpectrumFit #881

Merged
merged 4 commits into from Feb 14, 2017

Conversation

@lmohrmann
Copy link
Contributor

@lmohrmann lmohrmann commented Feb 10, 2017

This PR fixes the stat computation in SpectrumFit for the case that observation with different energy binning are given and add a regression test.

@cdeil cdeil requested a review from joleroi Feb 10, 2017
@cdeil cdeil added the bug label Feb 10, 2017
@cdeil cdeil added this to the 0.6 milestone Feb 10, 2017
Copy link
Contributor

@joleroi joleroi left a comment

Thanks for these fixes, @lmohrmann . I have left a comment for both changes you made. In summary, it would be very helpful to have a failing code snipped that shows what problem you encountered exactly. Both to fix it (obvioulsy - already done by this PR, I guess) but also to improve our regression tests that potentially did not cover a use case that led to the errors you saw.

@@ -204,8 +204,11 @@ def _predict_counts_helper(self, obs, model, forward_folded=True):
e_reco=binning)
counts = temp.data.data
else:
# TODO: This could also be part of calculate predicted counts
counts = model.integral(binning[:-1], binning[1:])
temp = calculate_predicted_counts(model=model,

This comment has been minimized.

@joleroi

joleroi Feb 10, 2017
Contributor

The forward_folding=False flag was meant for the case, when no IRFs are given at all (i.e. the model is directly in counts space). The way you changed it, now forsees the case that only an effective area is given, right? This should have worked before by setting forward_folding=True and edisp=None. Didn'it it work for you? If it worked before, we should probably just improve the documentation. If not, we have to discuss again.

This comment has been minimized.

@lmohrmann

lmohrmann Feb 10, 2017
Author Contributor

I assumed that model.integral(binning[:-1], binning[1:]) always returns a flux. If it returns counts the code is probably fine. The documentation should be improved in this case I guess.

This comment has been minimized.

@joleroi

joleroi Feb 10, 2017
Contributor

model.integral returns whatever the unit of the model amplitude is (times energy). Do you want to update the docstring of SpectrumFit in this PR, or should I do this in a separate one?

This comment has been minimized.

@lmohrmann

lmohrmann Feb 10, 2017
Author Contributor

I added a commit that reverts the change and includes an updated docstring

This comment has been minimized.

@joleroi

joleroi Feb 10, 2017
Contributor

Thanks 👍

@@ -88,7 +88,7 @@ def __init__(self, fit):
def _calc(self, data, model, *args, **kwargs):
self.fit.calc_statval()
# Sum likelihood over all observations
total_stat = np.sum(self.fit.statval, dtype=np.float64)
total_stat = np.sum([np.sum(v) for v in self.fit.statval], dtype=np.float64)

This comment has been minimized.

@joleroi

joleroi Feb 10, 2017
Contributor

Consider the following code example

In [7]: array1 = np.array([1,2,3])

In [8]: array2 = np.array([1,2,3])

In [9]: listofarrays = [array1, array2]

In [10]: np.sum(listofarrays)
Out[10]: 12

So np.sum should be able to deal with a list of arrays. Do you have an example where you got an wrong fit for more than one observation?

This comment has been minimized.

@lmohrmann

lmohrmann Feb 10, 2017
Author Contributor

This does not work if the two observations have a different number of bins, like so:

In [9]: array1 = np.array([1,2,3])

In [10]: array2 = np.array([1,2,3,4])

In [11]: listofarrays = [array1, array2]

In [12]: np.sum(listofarrays)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-12-e04d0a9e51c5> in <module>()
----> 1 np.sum(listofarrays)

/Users/larsm/Software/anaconda2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in sum(a, axis, dtype, out, keepdims)
   1842         except AttributeError:
   1843             return _methods._sum(a, axis=axis, dtype=dtype,
-> 1844                                 out=out, **kwargs)
   1845         return sum(axis=axis, dtype=dtype, out=out, **kwargs)
   1846     else:

/Users/larsm/Software/anaconda2/lib/python2.7/site-packages/numpy/core/_methods.pyc in _sum(a, axis, dtype, out, keepdims)
     30 
     31 def _sum(a, axis=None, dtype=None, out=None, keepdims=False):
---> 32     return umr_sum(a, axis, dtype, out, keepdims)
     33 
     34 def _prod(a, axis=None, dtype=None, out=None, keepdims=False):

ValueError: operands could not be broadcast together with shapes (3,) (4,) 

In principle it could be possible to have observations with different binnings...

This comment has been minimized.

@joleroi

joleroi Feb 10, 2017
Contributor

I see, good point (different binning will definitely happen when doing joint fit with different instruments). I think it would be good to add this scenario to the regression tests. So we would need to add something here. Since I know these tests quite well, it would be quick for me to do it. If you want to get into the gammapy test, however, you're very welcome to do so. I could give you some hints how to implement this test, if you wish.

This comment has been minimized.

@lmohrmann

lmohrmann Feb 10, 2017
Author Contributor

Please feel free to go ahead and add the test ;-)

This comment has been minimized.

@joleroi

joleroi Feb 10, 2017
Contributor

Alright!

@joleroi joleroi changed the title two fixes in spectrum module Support observations with different energy binning in SpectrumFit Feb 10, 2017
@joleroi
Copy link
Contributor

@joleroi joleroi commented Feb 10, 2017

I changed the title of the PR, because @cdeil uses them for the changelog (and he's very picky about the words a PR title may start with 😉 )

@cdeil
Copy link
Member

@cdeil cdeil commented Feb 10, 2017

I changed the title of the PR, because @cdeil uses them for the changelog (and he's very picky about the words a PR title may start with 😉 )

Really, that's all I care about. Any code changes are fine with me.

@joleroi
Copy link
Contributor

@joleroi joleroi commented Feb 10, 2017

I updates the description this PR and added a regression test. Let's merge if tests pass...

@joleroi
Copy link
Contributor

@joleroi joleroi commented Feb 14, 2017

Remaining test failures are related to gamma cat

@joleroi joleroi merged commit c91132a into gammapy:master Feb 14, 2017
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build could not complete due to an error
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@joleroi
Copy link
Contributor

@joleroi joleroi commented Feb 14, 2017

@lmohrmann Thanks a lot for your first contribution 🎆

@cdeil cdeil changed the title Support observations with different energy binning in SpectrumFit Add support for observations with different energy binning in SpectrumFit Feb 18, 2017
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