Add keepdim option for maps

[AtreyeeS]

This PR adds an option for squashing an axis into a corresponding axis with only one bin

[adonath]

Thanks, @AtreyeeS! I've left a few minor in-line comments. Would you like to have this PR in Gammapy 0.9?

[adonath]

Can you add asserts on the edge min and max values here?

[AtreyeeS]

Ok. Also added assert on the center.

[adonath]

Maybe move this line to an else block below? I think this improves the code structure a bit.

[AtreyeeS]

I think an else block is not necessary here.
But I have moved the data = np.nansum(self.data, axis=axis, keepdims=keepdims) to after the if block, which does make the code clearer.

[adonath]

Please adapt the docstring here and fix the typo.

[AtreyeeS]

Done. Also fixed typo in preceeding function

[AtreyeeS]

Thanks @adonath . Made the changes as you said. I dont think there is any special need to put this in v0.9, so if you are ready with the release already, this can wait.

[adonath]

@AtreyeeS Just one last comment: could you extend the test_wcsndmap_sum_over_axes() in test_wcsnd.py to cover the option keepdims=True and add an additional assert? Just to make sure we have test coverage for sum_over_axes(keepdims=True) as well. Thanks!

[adonath]

And what about #1935? Do you want to implement the option for MapMaker.make_images(keepdims=True) as well in this PR? Or a separate one?

[AtreyeeS]

And what about #1935? Do you want to implement the option for MapMaker.make_images(keepdims=True) as well in this PR? Or a separate one?

Thanks for this. I had added the option, but missed doing a git add on those files!

[cdeil]

Shouldn't axis.squash preserve the node_type?

[adonath]

@cdeil I've thought about this as well. But apart from #1952, there is also the question where to set the bin center if you sum over the axes. I'm not sure if there is a unique answer. Summing over the axes also mainly applies to integral quantities such as counts or flux, so I'm not sure if there is a use case for differential quantities at all.

But I agree we should not silently change the node_type. So maybe raise an error if node_type="center" and keepdims=True for now?

[cdeil]

Summing over the axes also mainly applies to integral quantities such as counts or flux, so I'm not sure if there is a use case for differential quantities at all.

I think the main use case (see #1935) is to allow getting a flux map with correct value via the MapMaker, and fitting a source on 2D maps to obtain a source flux, like what we did in HGPS.

Personally I think only supporting that kind of analysis by running 3D with multiple bins and fixing spectral paramters except for amplitude would be fine, because getting the 2D "right" is a headache. Probably we should use node_type="center" for exposure, and also the maps that represent the PSF? I remember this was extensively discussed before, but I don't know what the status of our code is.

But I agree we should not silently change the node_type. So maybe raise an error if node_type="center" and keepdims=True for now?

Or maybe keep as-is to allow using this, and add a comment at the top of the squash method?

# TODO: decide how to handle `node_type="center"`
# See https://github.com/gammapy/gammapy/issues/1952

[AtreyeeS]

Shouldn't axis.squash preserve the node_type?

I thought about this while coding, and reached #1952 . But then, I argues that this discrepancy does not matter too much for the purpose of MapMaker.

Personally I think only supporting that kind of analysis by running 3D with multiple bins and fixing spectral paramters except for amplitude would be fine

No, this is definitely not fine. For low flux sources, it is offensive to ask users to make too many bins in e_reco and then have zero counts in many bins and end up with wrong statistics.

Or maybe keep as-is to allow using this, and add a comment at the top of the squash method?

I think this a good idea for the present, and then, once we decide the fix for #1952, we fix the TODO

[cdeil]

No, this is definitely not fine. For low flux sources, it is offensive to ask users to make too many bins in e_reco and then have zero counts in many bins and end up with wrong statistics.

The cash / cstat likelihood that is used for map analyses doesn't have this issue, only WSTAT, which is only used for spectra and there we already support it.

But no need to argue - if you take care of figuring out the node_type and single bin case and make it work in Gammapy - that's great!

There is still the question if special-casing is needed for 2D maps, especially to support ring background estimation. We talked about this in Madrid and even did a bit of prototyping (#1850). Not sure if we reached a conclusion on what to do. If it's possible to support 2D analysis as you like by making it work with one MapMaker / MapFit, IMO that would be much better than to develop separate data reduction / fitting code just for the 2D case. Although if some algorithm is only used for 2D, it's fine to have a separate function or class for it - that's already the case e.g. for making TS maps.