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 MapDatasetOnOff #2481

Merged
merged 13 commits into from Oct 29, 2019
Merged

Add MapDatasetOnOff #2481

merged 13 commits into from Oct 29, 2019

Conversation

@luca-giunti
Copy link
Contributor

luca-giunti commented Oct 23, 2019

Implement a map dataset containing on-off counts and acceptance.

Some tests: https://github.com/luca-giunti/share/blob/master/ONOFF_test.ipynb

@luca-giunti luca-giunti requested a review from adonath Oct 23, 2019
@cdeil cdeil added the feature label Oct 23, 2019
@cdeil cdeil added this to the 0.15 milestone Oct 23, 2019
@cdeil cdeil added this to To do in gammapy.maps via automation Oct 23, 2019
@cdeil cdeil changed the title Introduce map dataset on off Add MapDatasetOnOff Oct 23, 2019
@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Oct 23, 2019

@luca-giunti @adonath @registerrier - Any chance to reduce the amount of code (e.g. by using extending methods in the parent classes, or using helper functions called from both map dataset classes)?

Just a thought. I'm unsubscribing now.

Copy link
Member

adonath left a comment

Thanks al lot @luca-giunti! I've left a few inline comments how to reduce code-duplication between MapDataset and MapDatasetOnOff. I think once we make the (preliminary) assumption MapDatasetOnOff.background_model=None we can re-use a lot of methods from the MapDataset.

I think in general it is likely, that we remove MapDataset.background_model anyway and just handle it as part of the "normal" model...

counts_off : `~gammapy.maps.WcsNDMap`
Ring-convolved counts cube
acceptance : `~gammapy.maps.WcsNDMap` or float
Background from the IRFs

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

Let's keep consistency here and maybe not call it background but just acceptance. While technically it is the same, it is not used as a background model...

acceptance : `~gammapy.maps.WcsNDMap` or float
Background from the IRFs
acceptance_off : `~gammapy.maps.WcsNDMap` or float
Ring-convolved background from the IRFs

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

See my comment above, just name it "Acceptance off"

self.mask_safe = mask_safe
self.gti = gti

def __str__(self):

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

It think once we set background_model=None the MapDatasetOnOff.__str__ is the same as MapDataset.__str__? In this case there is no need to re-implement it. For now I would propose to just call super().__str__ and attach a little information about the acceptance, alpha and counts_off at the end. See here, for what we did for the SpectrumDatasetOnOff. While this is certainly not an ideal solution, it is pragmatic for now and avoids most code duplication.

This comment has been minimized.

Copy link
@luca-giunti

luca-giunti Oct 24, 2019

Author Contributor

OK. Notice that this means that also the .npred() method from the MapDataset is called, which is fine if we set background_model=None on the MapDatasetOnOff. I guess I can also eliminate the .npred_sig() on the MapDatasetOnOff and just use the .npred() inherited from the MapDataset...

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

Yes, I think this would be preferred solution for now, just do avoid the code duplication.

psf_map = Map.from_geom(geom_rad, unit="sr-1")
psf = PSFMap(psf_map, exposure_irf)

return cls(

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

Is it possible to simplify .create() to just call .from_geoms() and avoid e.g. re-creation of the data of the EDispMap?

This comment has been minimized.

Copy link
@luca-giunti

luca-giunti Oct 24, 2019

Author Contributor

Sure! I mean I just duplicated from the MapDataset. We should probably what you suggest there as well?

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

The change in MapDataset is already part of #2479, so no need to adapt it here...

hdulist : `~astropy.io.fits.HDUList`
Map dataset list of HDUs.
"""
# TODO: what todo about the model and background model parameters?

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

Again you can call super().to_hdulist() here and just extend the hdulist with the "off" quantities, such as "acceptance", "acceptance_off", etc.

other: `~gammapy.cube.MapDatasetOnOff`
Dataset to be stacked with this one.
"""
if self.counts and other.counts:

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

I guess you can just call super().stack(other) here and just keep the code to stack the off quantities...

if self.gti and other.gti:
self.gti = self.gti.stack(other.gti).union()

def likelihood(self):

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

Methods implemented on the base class, are inherited and don't have to be re-implemented. Just remove .likelihood()...

This comment has been minimized.

Copy link
@luca-giunti

luca-giunti Oct 24, 2019

Author Contributor

I think we need to skip the MapDataset parent class and go back to the general Dataset class for this. This because the MapDataset.likelihood() invokes some cython functions that work under the assumption of cash (

def likelihood(self):
)

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

Sorry, of course you are right. Just leave as is...

"""Total likelihood given the current model parameters."""
return Dataset.likelihood(self)

def fake(self, random_state="random-seed"):

This comment has been minimized.

Copy link
@adonath

adonath Oct 24, 2019

Member

I guess .fake() should also simulate the counts_off. In SpectrumDatasetOnOff, we decided to just pass a background_model. Maybe we do the same here?

gammapy.maps automation moved this from To do to In progress Oct 24, 2019
@adonath

This comment has been minimized.

Copy link
Member

adonath commented Oct 24, 2019

Note that currently there is a test-fail, because of a missing import of WcsNDMap, see https://travis-ci.org/gammapy/gammapy/jobs/601875197#L3089

@luca-giunti

This comment has been minimized.

Copy link
Contributor Author

luca-giunti commented Oct 24, 2019

Thanks @adonath for the review!

I addressed most of your comments. However I get 2 errors which after spending sufficient time I don't understand, maybe you might have an idea... one occurs in the test for the .to_hdulist() method and the other in the test of the stacking.

The error I get in the stacking part is especially puzzling. The strange behavior I get is the following: https://github.com/luca-giunti/share/blob/master/Atest.ipynb

I guess I am not making a right use of super()...

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Oct 24, 2019

Thanks @luca-giunti, I'll take a look now...

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Oct 24, 2019

@luca-giunti The usage of super() was basically correct. In e347057 I adapted the I/O test. The special behaviour there was that the test date are actually images, so no _BANDS HDUs are written. Once I adapted the reference HDU names the tests passed.

Taking a look at the implementation of the .stack() method I'm not sure it's correct (at least it seems different from what we've done in SpectrumDatasetOnOff.stack() ). @registerrier Could you maybe check the implementations for consistency?

Copy link
Contributor

registerrier left a comment

There is indeed an issue in the stack method for the computation of the acceptance_off. The mask should be taken into account and the division by the total OFF seems to be missing.

self.alpha.data * self.counts_off.data
+ other.alpha.data * other.counts_off.data
)
self.acceptance.data = np.ones(self.data_shape)

This comment has been minimized.

Copy link
@registerrier

registerrier Oct 24, 2019

Contributor

This looks incorrect indeed.

  • mask should be applied to compute the average alpha
  • the average alpha has to be divided by the total OFF
@luca-giunti

This comment has been minimized.

Copy link
Contributor Author

luca-giunti commented Oct 25, 2019

Thanks @adonath and @registerrier !

Regarding the stacking, indeed there were mistakes. However, this issue that I encountered in stacking a dataset with an identical one still persists... and moreover, it is not related to the .stack() in the MapDatasetOnOff but with the one in the MapDataset.

For some reason, the other dataset is changed as well as the self...

Here I reproduced the issue, using gammapy-0.14 (it is there also in the last dev...): https://github.com/luca-giunti/share/blob/master/MapDataset_stacking.ipynb

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Oct 25, 2019

@luca-giunti Actually this behaviour is expected when you stack a dataset with itself, because the data is modified in place. If dataset.data and other.data refer to the same numpy data buffer in memory, then once you change dataset.data, other.data is changed as well. I think in practice this is rarely a problem, but in your specific case it appears, because of the way you create the test datasets

So what you need to do is an explicit copy of the first dataset, before stacking e.g.like so:

dataset = MapDataset()

stacked = dataset.copy()
stacked.stack(dataset)
@luca-giunti

This comment has been minimized.

Copy link
Contributor Author

luca-giunti commented Oct 25, 2019

Thanks @adonath. The stacking should be better now, and I also added a test for the stacking of cutouts

@luca-giunti luca-giunti force-pushed the luca-giunti:Introduce-MapDatasetOnOff branch from b1805e3 to 7845a37 Oct 28, 2019
luca-giunti and others added 4 commits Oct 28, 2019
@adonath adonath force-pushed the luca-giunti:Introduce-MapDatasetOnOff branch from e38e3fa to 6012953 Oct 28, 2019
adonath added 2 commits Oct 29, 2019
@adonath adonath merged commit a658474 into gammapy:master Oct 29, 2019
8 of 9 checks passed
8 of 9 checks passed
Codacy/PR Quality Review Not up to standards. This pull request quality could be better.
Details
Scrutinizer Analysis: 2 new issues, 25 updated code elements – Tests: passed
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
gammapy.gammapy Build #20191029.4 succeeded
Details
gammapy.gammapy (DevDocs) DevDocs succeeded
Details
gammapy.gammapy (Lint) Lint succeeded
Details
gammapy.gammapy (Test Python36) Test Python36 succeeded
Details
gammapy.gammapy (Test Windows36) Test Windows36 succeeded
Details
gammapy.gammapy (Test Windows37) Test Windows37 succeeded
Details
gammapy.maps automation moved this from In progress to Done Oct 29, 2019
@adonath

This comment has been minimized.

Copy link
Member

adonath commented Oct 29, 2019

Thanks @luca-giunti!

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