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

Enable rate models in SpectrumSimulation #1337

Merged
merged 3 commits into from Mar 13, 2018

Conversation

Projects
None yet
2 participants
@joleroi
Copy link
Contributor

commented Mar 12, 2018

As discussed in #1330 the SpectrumSimulation class should be able to handle rate models. This PR adds this options. It is still not possible however, to have different kinds of models for source and background. This should be very easy to implement, however. I'll leave it for another PR. I don't think a review is needed but if someone want to give his 2 cents, you're welcome.

@joleroi joleroi added the feature label Mar 12, 2018

@joleroi joleroi self-assigned this Mar 12, 2018

@joleroi

This comment has been minimized.

Copy link
Contributor Author

commented Mar 12, 2018

Note that with this PR it's also necessary to specify a lifetime in the SpectumFit. This avoids disambiguities when fitting rate models (weather the lifetime should be multplied to the model before or during the fit). In the long run we should rewrite SpectrumFit anyway to accept only models with unit ( count / energy ). All the forward folding and livetime stuff should be done before on the model (like in Sherpa). Then this problem will go away.

@joleroi

This comment has been minimized.

Copy link
Contributor Author

commented Mar 13, 2018

Test failures are unrelated

@joleroi joleroi merged commit 7a1e501 into gammapy:master Mar 13, 2018

1 of 2 checks passed

continuous-integration/travis-ci/pr The Travis CI build failed
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details

@joleroi joleroi deleted the joleroi:rate_model branch Mar 13, 2018

e_true=e_true
)
sim.simulate_obs(seed=23, obs_id=23)
assert sim.obs.on_vector.total_counts == 10509 * u.ct

This comment has been minimized.

Copy link
@cdeil

cdeil Mar 15, 2018

Member

I thought we decided to not use u.ct or u.ph throughout Gammapy?
@joleroi - OK? Would you be willing to remove the u.ct here or in other places in gammapy.spectrum where it still exists?

This comment has been minimized.

Copy link
@joleroi

joleroi Mar 16, 2018

Author Contributor

I can't remember that we had a discussion about this, let alone a decision. Can you point me somewhere?

This comment has been minimized.

Copy link
@cdeil

cdeil Mar 16, 2018

Member

There is a discussion about pros and cons of having u.ph and / or u.ct or not here: #978
Do you dislike the solution to always omit them?
I think that is what we mostly have in Gammapy at the moment, and as I said in #978 I think it's the simplest solution.

This comment has been minimized.

Copy link
@joleroi

joleroi Mar 16, 2018

Author Contributor

I missunderstood, I though you don't like the way the astropy units package is used. I'm in favour of always omitting u.ph or u.ct but I'm afraid gammapy.spectrum is full of them. I can take care of it.

This comment has been minimized.

Copy link
@cdeil

cdeil Mar 16, 2018

Member

@joleroi - It's not so bad, I think it can be done in one pull request:

$ ack "u\.ct" gammapy --type=python
gammapy/spectrum/simulation.py
186:        self.on_vector.data.data += bkg_counts * u.ct

gammapy/spectrum/results.py
238:        data = self.npred_src * u.ct
248:            data = self.npred_bkg * u.ct

gammapy/spectrum/tests/test_core.py
16:        self.counts = [0, 0, 2, 5, 17, 3] * u.ct
26:        assert spec.data.data.unit.is_equivalent(u.ct)
32:        assert spec.data.data.unit.is_equivalent(u.ct)
66:        desired = [0, 7, 20] * u.ct
89:        assert spec.data.data.unit.is_equivalent(u.ct)
94:        assert spec.data.data.unit.is_equivalent(u.ct)

gammapy/spectrum/tests/test_simulation.py
44:        assert self.sim.obs.on_vector.total_counts == 160 * u.ct
51:        assert self.sim.obs.on_vector.total_counts == 530 * u.ct
52:        assert self.sim.obs.off_vector.total_counts == 1112 * u.ct
58:        assert self.sim.result[0].on_vector.total_counts == 158 * u.ct
59:        assert self.sim.result[1].on_vector.total_counts == 158 * u.ct
60:        assert self.sim.result[2].on_vector.total_counts == 161 * u.ct
61:        assert self.sim.result[3].on_vector.total_counts == 168 * u.ct
62:        assert self.sim.result[4].on_vector.total_counts == 186 * u.ct
70:        assert sim.obs.on_vector.total_counts == 161 * u.ct
85:        assert sim.obs.on_vector.total_counts == 10509 * u.ct

gammapy/spectrum/core.py
49:        data = np.arange(10) * u.ct
148:        self.data.data = binned_val * u.ct
267:        retval.data.data = counts_rebinned * u.ct
488:            data=counts_table['COUNTS'] * u.ct,
625:            kwargs['data'] = row['COUNTS'] * u.ct
$ ack "u\.ph" gammapy --type=python
<nothing>

Please add a changelog entry, as this likely is a change that will break a user script or two.

@cdeil

This comment has been minimized.

Copy link
Member

commented Mar 15, 2018

@joleroi - Please always put a milestone for new issues or PRs - that helps when filling the changelog.
For PRs it's always the next upcoming release, in this case v0.8, and then if the PR stalls, it gets moved to the next release or closed and put on the "never" milestone. I'll do it here now.

@cdeil cdeil added this to the 0.8 milestone Mar 15, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.