Background estimation map

[facero]

Hi,
While preparing my notebooks for the DC1 gammapy/ctools comparison, I reran my notebook that I used for the PHYS meeting in Sept 2017.
The background map with the latest gammapy version (0.7.dev5223) produces a background map with higher (more diffuse/smoothed) background values than the map that I produced in september with 0.7.dev5025.
I reverted my current version to 0.7.dev5025 (67babf6) and the background map is similar to the one obtained in september.

There has been quite some modifications since september. Has any of you noticed this change ?
I'm not sure which background image is the more realistic but the change is quite significant and I was wondering whether someone had checked this.
In particular there is some kind of square grid pattern that is visible.

@bkhelifi when you implemented changes during the last coding sprint did you see changes in the background maps ?

0.7.dev5223

0.7.dev5025 (67babf6)

For a minimal example, you can use the notebook here:
#1143

and use this command:
images['background'].smooth(radius=0.1).plot(cmap='viridis',add_cbar=True)

[adonath]

Note that a few modifications have been made to the RingBackgroundEstimator in #1166, that probably explain the changes you see. Still the rectangular pattern looks suspicious to me. The ring background estimation is still not tested very well still does not take the instrument acceptance into account. I'll make an issue that the RingBackgroundEstimator class should be improved further...

[facero]

Thanks for the info. I produced the bkg map with 8bfae5c and it has no grid.
As you mention, the grid feature appears right after #1166 (I tested right before/after).
Here is a bkg map without mosaicing that shows the grid effect:

The grid seems to be centered on the ref_image defined by SkyImage.empty.
When changing the center of the empty image but keeping the same run list, the grid stays in the center of the image:

So the effect seems not FoV related but rather SkyImage related.

[cdeil]

Is this using FFT convolution? These artifacts look like something FFT convolution could introduce, no?

Probably writing a few tests to check if SkyImage.convolve works for FFT convolution (i.e. gives almost the same output) should be done first (including a test for behaviour at the edge). Then fixes / tests for ring background estimation can rely on convolve working properly and don't need to be duplicated for the FFT and direct convolve case.

Here's what we have for SkyImage.convolve at the moment:

gammapy/gammapy/image/core.py

Line 1175 in c421ba9

def convolve(self, kernel, use_fft=False, **kwargs):

Note that in Astropy a lot of time was spent over the past year to get the same output from direct and FFT convolution. Initially it was very buggy (i.e. different behaviour for NaN values, edges, kernel normalisation, ...), but I think by now it's really a common and nice / well-tested uniform interface for both FFT and direct convolve:
http://astropy.readthedocs.io/en/stable/convolution/index.html
I was against using it in Gammapy before, because it was buggy, and because we didn't really need FFT convolution and scipy.ndimage.convolve was simple / correct / fast.

So does anyone have time to tackly the convolve and then ring bg estimate situation and bring us towards a good, well-understood and well-tested situation? @adonath - I think for someone that knows about the methods it should be a day or two of work, and for someone new to Gammapy development / writing tests or even convolution it's a week of work?

[cdeil]

Just one more thing I forgot to mention: one concern I would have with switching to astropy.convolve is speed, at least for direct convolution. If possible, I would like to stick with scipy.ndimage.convolve for the direct convolution case which we've used for many years and I know is fast and has the edge options and N-dim we need, unless someone shows that astropy.convolve has similar performance. This is one case where performance really matters, a good number of analyses will be dominates in CPU time by speed of the convolve call.

[facero]

This doesn't seem related to FFT. The examples above were generated with use_fft_convolution=False. the problem seems to come from the convolve mode which was changed in #1166 from mode='reflect' to mode='constant'.
Reverting the option reflect fixes the issues but I guess @adonath changed it for a good reason.

Concerning the FFT, one interessting thing is that when keeping the code as is (mode='constant') and using FFT the problem disappears.
The problem happens only when using use_fft_convolution=False and mode='constant'.

For speed for 5 observations and an image of 200x200 pixels, fftconvolve runs in 3.7 sec and convolve in 8.4sec on my laptop.

[adonath]

@facero Now, I realize what the problem is. Could you increase the size of the reference sky image, so that the FoV is fully contained (and maybe include an extra margin of 0.5 deg) and see if that fixes the problem?

Actually the boundary handling of the convolution should not influence the result, because the code implicitly requires that the image is much larger than the FoV, so the part of the image, that is affected by the boundary handling is always outside the FoV. This is no the case in your example and we have to think about how to handle this. We could simply require a minimum size of the reference image. But it is probably best to adapt the algorithm and pad the image by half the kernel size prior to the ring convolution and crop it afterwards again. Does that make sense?

[facero]

I will try that.
Why would the problem go away when using FFT ?

I know that astropy.fftconvolve does padding (I don't remember if it is by default).

[cdeil]

Why is padding needed? We want to have the ring be a counting kernel. So https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.convolve.html with mode='constant' and the default cval=0 should be used for the direct convolution case, no?
Is it possible for the FFT convolution case also?

If we can avoid padding, it'll avoid having to introduce special code for handling different cases of where the background is in the image and different kernel / image sizes.

[adonath]

@cdeil The problem here is a different one: the FoV is not fully contained in the reference image and instead of using the full information for the background computation (i.e. the data in the FoV that is outside the reference image) the zero padding is used for the convolution, which causes the artefacts.

There are now different ways to handle this problem:

• require a minimal size for the reference image so that all runs are fully contained
• internally pad the reference image to the required size and crop again
• change the algorithm to always compute a correctly sized image for every run and stitch them together in the reference image

[cdeil]

@adonath - I don't think we should pad anywhere.

So I read the RingBackgroundEstimator code and I think this is a bug, no?

gammapy/gammapy/background/ring.py

Line 331 in 0d475fc

result['exposure_off'] = exposure_on_excluded.convolve(ring.array, mode='reflect',

You're convolving exposure_on with mode='reflect', but you should do the same thing on on_exposure and you do to counts, i.e. mode='constant', no?

Another suggestion: change the test case so that the run isn't fully contained in the image:

gammapy/gammapy/background/tests/test_ring.py

Line 13 in 0d475fc

def example_test_images():

and then add asserts on a pixel at the edge to establish what the algorithm does.
To debug it, probably a notebook showing images would be better than just running pytest and printout.

@facero or @adonath - Does one of you have time to start a pull request soon? Or should I do it?

[facero]

I tried @adonath suggestion and increasing the size does not fix the issue.

No FFT

With FFT

The problem is fixed when using the FFT option, this is an interesting fact and could help us to understand what is going on.

As @cdeil mentions, the problem seems related to the mode options used here.

[cdeil]

Here is another bug in the image computation (not really related to the feature discussed here, but could be fixed in the same PR if someone starts one with the goal to get correct images):

gammapy/gammapy/image/basic.py

Line 223 in 0d475fc

def _counts(self, observation):

Instead I think one has to process counts like this, to make sure the FOV offset cut is 100% consistent between exposure and counts, and the lines at the FOV edge with counts / exposure mismatch are gone for sure:
https://github.com/gammasky/cta-dc1-gps-analysis/blob/d91345301665cef41986f786c127578b03708339/make_cube.py#L64
Most of that code should be moved over to Gammapy, but especially a function like this should be called for images and cubes to process counts correctly:
https://github.com/gammasky/cta-dc1-gps-analysis/blob/d91345301665cef41986f786c127578b03708339/make_cube.py#L468

[adonath]

The problem @cdeil mentioned with mode='reflect' comes on top to the one I mentioned. When you change the line to mode='constant', the results between FFT and normal convolution should be consistent. But both still give incorrect results, when the FoV is not fully contained in the reference image and the question how to handle the latter is still open...

@cdeil Note that the data outside the FoV of the counts image is handled here:

gammapy/gammapy/image/basic.py

Line 293 in 0d475fc

not_has_exposure = ~(exposure.data > 0)

But I agree this is better handled inside _counts() using data[offset >= offset_max] = 0 but the result will be the same...

[cdeil]

but the result will be the same...

I think you are right! I just quickly read the code and thought you had a "select_cone" call on the events, but actually there is only a "select_energy" call:

gammapy/gammapy/image/basic.py

Line 233 in e167145

events = observation.events.select_energy((p['emin'], p['emax']))

[cdeil]

But both still give incorrect results, when the FoV is not fully contained in the reference image and the question how to handle the latter is still open...

This is what I don't understand: yes, usually the caller should compute images so that all data is fully contained. But if there are runs at the edge of the image, then counts / exposure / background computations should still work, just loose some data. You know the code better than me, maybe indeed there is another bug in the background computation code. But the solution would be to fix that without introducing padding. Suggest to first fix the obvious bug "reflect" -> "constant" and then look at an example image.

[lmohrmann]

Before having read this issue, I also found the bug that @cdeil discovered and opened a PR for this: #1204.

[cdeil]

I talked to @adonath about this offline:

• With Consistently use mode='constant' in convolutions of RingBackgroundEstimator #1204 merged results should be OK now.
• The discussion above about a second bug isn't really a bug, it's just the point that at the moment the per-run images aren't extended to cover the whole run FOV. This could stay as-is, but documentation should be improved to explain this (how run position selection and sky image region relate)
• As a concrete action item: the test for ring background estimation should be changed to have a run intersect the edge of the sky image, so that this effect appears in the reference results.

Keeping this issue open for now.

[cdeil]

@registerrier has started to work on new code to make maps that should resolve the issues discussed here. Regis - assigning this issue to you.

[cdeil]

I'm closing this old issue. I think it has been resolved in the old image estimator a few months ago, and anyways that code will be deleted very soon.

Please use gammapy.maps and the new MapMaker (see tutorial), and do report any issue you find there!