Implement HpxNDMap.convolve() #3323

LauraOlivera · 2021-04-27T16:44:07Z

This pull requests implements several methods required for the convolution of a HealPix map with a kernel. There are several options to do this: the kernel being a WCS or Hpx map, the convolution happening in a small cutout or in the full sky map.

For now, I've implemented convolve_wcs() which is the "simplest" of these approaches in my opinion. For a "small map" (I set a limit of 10 degree width, but we can also refine that), it uses the HpxToWcsMapping to project the initial HpxNDMap into a WcsNDMap, convolves with a kernel using the existing WcsNDMap.convolve() and makes use of the mapping to go back to the inital HpxGeom. To make this work this I implemented a method HpxToWcsMapping.fill_hpx_map_from_wcs_data() which just does the inverse as the existing HpxToWcsMapping.fill_wcs_map_from_hpx_data(). Testing for this is not implemented yet. (TO DO)

The caveat is that the geometry of the kernel has to be compatible with the intermediate WcsGeom, and so the kernel has to be extracted using hpx_geom.to_wcs_geom() where hpx_geom is the inital map geometry. That means, that if this was used for PSF convolution in a MapDataset with Healpix geometry, the geometry of the PSFKernel (and Map?) would need to be different than the rest of components.

TO DO: The next step would be to, from an inital kernel in an unspecified geometry, extract the radial profile of the kernel in angular space (assumed symmetric), convert to a beam function using healpy.sphtfunc.beam2bl() and then make use of healpy.sphtfunc.smoothing. This happens in a full-sky map, which would require a cutout to be stacked in a full-sky map as an initial step. I suspect this will be much slower than the case above, but I have to check. For a kernel also in an HpxMap, there might be faster options.

codecov · 2021-04-27T16:56:44Z

Codecov Report

Merging #3323 (ac6cdac) into master (4571366) will decrease coverage by 0.05%.
The diff coverage is 91.76%.

@@            Coverage Diff             @@
##           master    #3323      +/-   ##
==========================================
- Coverage   93.83%   93.78%   -0.06%     
==========================================
  Files         144      144              
  Lines       18050    18198     +148     
==========================================
+ Hits        16938    17067     +129     
- Misses       1112     1131      +19

Impacted Files	Coverage Δ
gammapy/maps/hpx.py	`88.53% <86.66%> (+0.08%)`	⬆️
gammapy/maps/hpxnd.py	`94.57% <92.85%> (-0.14%)`	⬇️
gammapy/modeling/scipy.py	`87.32% <0.00%> (-12.68%)`	⬇️
gammapy/makers/safe.py	`91.02% <0.00%> (-2.22%)`	⬇️
gammapy/modeling/fit.py	`95.33% <0.00%> (-0.41%)`	⬇️
gammapy/modeling/models/spectral.py	`96.52% <0.00%> (-0.22%)`	⬇️
gammapy/maps/wcsmap.py	`96.29% <0.00%> (-0.11%)`	⬇️
gammapy/irf/effective_area.py	`89.36% <0.00%> (+0.11%)`	⬆️
gammapy/estimators/flux_point.py	`93.80% <0.00%> (+0.13%)`	⬆️
... and 5 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4571366...ac6cdac. Read the comment docs.

adonath

Thanks lot @LauraOlivera, I have left some first comments, I'll check out the functionality locally now...

adonath · 2021-04-29T08:02:44Z

gammapy/maps/hpxnd.py

+        map : `HpxNDMap`
+            Convolved map.
+        """
+        if self.geom.width < 10*u.deg:


I think the number "10 deg" should be introduce as a class level variable like:

class HpxNDMap: convolve_wcs_max_width

This allows users in general to overwrite it, in case they know what they are doing...

Or maybe even better: alternatively we could just say all partial maps are supported and we make a statement in the docstring on the accuracy for larger maps >~10 deg

I agree. My idea was that ideally we would have this method (ok for small maps) and then a full-sky map one that involves stacking a partial map into a full-sky map if needed (which I expect to be slower, but if it isn't then I need to rethink things), and if you try to use convolve_wcs on a partial map too large, you just get a warning suggesting that the other one is more suited or something. But it needs some more thought/work.

adonath · 2021-04-29T08:07:43Z

gammapy/maps/hpxnd.py

+        """
+        if self.geom.width < 10*u.deg:
+            # Project to WCS and convolve
+            wcs_map  = self.to_wcs(fill_nan = False)


I think .to_wcs() takes a pre-computed mapping as well. There might be a performance improvement if the HpxToWcsMapping is created outside and re-used for both cases i.e. .to_wcs() as well as the later inverse operation.

I think this even offers the possibly to adjust the WCS geometry used for the conversion and convolution, to the WCS geometry of the kernel, that is passed in. This eliminates the need to compute a matching kernel outside. I haven't thought this through completely, but I think it works...

oh, you are right! that is a very good point! I'll try that out

It works but I am seeing a lot more issues with artifacts that only go away if the wcs_geom binning is very fine. The nice thing about using .to_wcs() is that it guaranteed the oversampling... I guess it could be fixed by requiring a minimum wcs pixel size ratio to hpx pixel size...

%%time energy = MapAxis.from_bounds(1,100, unit='TeV', nbin=2, name='energy') nside = 1024 hpx_geom = HpxGeom.create(nside=nside, axes=[energy], region='DISK(0,0,2.5)', ) hpx_map = Map.from_geom(hpx_geom) hpx_map.set_by_coord((0,0, [2,90]), 1) wcs_geom = WcsGeom.create(width=5, binsz = 0.05, axes = [energy]) kernel = PSFKernel.from_gauss(wcs_geom, 0.5*u.deg) convolved_map = hpx_map.convolve_wcs(kernel) convolved_map.plot_interactive(add_cbar=True)

Thanks for check! Indeed we should require the oversampling to avoid the artefacts. As an educated guess, I'd say an oversampling of >2 should be sufficient. So +1 to just include the check and document it.

LauraOlivera · 2021-05-04T12:48:02Z

I added another method convolve_full (not convinced about the name because it also works for not full-sky maps so suggestions welcome). What this does is extract the radial profile from the input kernel (assuming radial symmetry), convert it into a window beam (like in PR #3314) and convolve using healpy.sphtfunc.smoothing. If the input map is not a full-sky map, it first inserts the smaller map into a full sky map, and then extracts the same region once the convolution is over.

Because the convolution happens in a full sky map, the performance is worse if the pixelization is fine. That is why, even if the input map is "small", for a usual analysis pixel size (nside = 1024 is the standard for HAWC, which corresponds to a 0.025deg pixel size), this method is much slower than convolve_wcs even if the intermediate WCS geometry has smaller pixels by a factor >2:

energy = MapAxis.from_bounds(1,100, unit='TeV', nbin=2, name='energy_true')
nside = 1024
hpx_geom = HpxGeom.create(nside=nside,
                         axes=[energy],
                         region='DISK(0,0,3.5)',
                         )
hpx_map = Map.from_geom(hpx_geom)
hpx_map.set_by_coord((0,0, [2,90]), 1)
wcs_geom = WcsGeom.create(width=5,
                          binsz = 0.01, 
                         axes = [energy])
psf = PSFMap.from_gauss(
    energy_axis_true=energy, sigma=[0.3, 0.2] * u.deg
)
kernel = psf.get_psf_kernel(geom=wcs_geom, max_radius=1*u.deg)

#### Using the WCS projection
%%time 
convolved_map_wcs = hpx_map.convolve_wcs(kernel)
CPU times: user 36.7 ms, sys: 30 µs, total: 36.8 ms
Wall time: 35.8 ms

#### Using the full-sky convolution
%%time 
convolved_map = hpx_map.convolve_full(kernel)
CPU times: user 1min 56s, sys: 3.22 s, total: 2min
Wall time: 10.9 s

So my suggestion is to implement now the method HpxNDMap.convolve() which essentially decides which of the two are called based on the size of the input map.

I also compared the maps produced above with one another and this is the result of doing convolved_map - convolved_map_wcs

adonath · 2021-05-04T13:19:52Z

Thanks a lot for investigating @LauraOlivera! From your results it seems, there is hardly any situation where .convolve_all() would be useful, except if it is really an all sky map to be convolved. For the npred evaluation we basically only need the local convolution, because source models are evaluated on small cutouts only. I'll proceed with a short code review and then we discuss in there afternoon call, how to support the different convolve methods (i.e. whether to use a single .convolve() method etc.)

adonath

Thanks @LauraOlivera, for now I just left two minor comments. Let's discuss in the next co-working meeting on the details...

gammapy/maps/hpxnd.py

adonath · 2021-05-04T18:22:16Z

@LauraOlivera To keep the scope of this PR limited maybe just implement .convolve() with the WCS method here? Than we could go ahead already and make MapDataset.npred() work. The .convolve_all() could be added in a follow up PR, that can be used to evaluate the other method and provide better information on the validity of the WCS approximation...

adonath

Thanks a lot @LauraOlivera! I think the code could need some polishing, however it's not a priority right now. Let's merge the PR and I'll prepare a follow up PR to do the clean up, while you and @vikasj78 can continue to work on the npred computation for HEALPix...

LauraOlivera added 4 commits April 27, 2021 17:15

expose fill_nan in to_wcs

3f631a3

create fill_hpx_map_from_wcs_data

75a36ff

create convolve_wcs

c475c06

add convolve_wcs test

70b3a35

adonath self-assigned this Apr 27, 2021

adonath added this to the v0.19 milestone Apr 27, 2021

adonath added this to To do in gammapy.maps via automation Apr 27, 2021

adonath added the feature label Apr 27, 2021

adonath requested changes Apr 29, 2021

View reviewed changes

gammapy.maps automation moved this from To do to In progress Apr 29, 2021

LauraOlivera added 3 commits May 4, 2021 12:10

add warning, take geometry from kernel

898178d

add convolve_full

6443e59

add test

a2424e3

adonath requested changes May 4, 2021

View reviewed changes

gammapy/maps/hpxnd.py Outdated Show resolved Hide resolved

gammapy/maps/hpxnd.py Show resolved Hide resolved

LauraOlivera added 3 commits May 5, 2021 11:37

create convolve

596b5c5

fix test

269baba

fix typo

ac6cdac

LauraOlivera marked this pull request as ready for review May 5, 2021 09:39

adonath approved these changes May 5, 2021

View reviewed changes

adonath merged commit 9c53c3f into gammapy:master May 5, 2021

gammapy.maps automation moved this from In progress to Done May 5, 2021

adonath mentioned this pull request May 5, 2021

Cleanup HpxNDMap.convolve_wcs #3341

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement HpxNDMap.convolve() #3323

Implement HpxNDMap.convolve() #3323

LauraOlivera commented Apr 27, 2021

codecov bot commented Apr 27, 2021 •

edited

Loading

adonath left a comment •

edited

Loading

adonath Apr 29, 2021

adonath Apr 29, 2021

LauraOlivera Apr 29, 2021

adonath Apr 29, 2021

adonath Apr 29, 2021

LauraOlivera Apr 29, 2021

LauraOlivera Apr 29, 2021

adonath Apr 29, 2021

LauraOlivera commented May 4, 2021

adonath commented May 4, 2021

adonath left a comment

adonath commented May 4, 2021

adonath left a comment

Implement HpxNDMap.convolve() #3323

Implement HpxNDMap.convolve() #3323

Conversation

LauraOlivera commented Apr 27, 2021

codecov bot commented Apr 27, 2021 • edited Loading

Codecov Report

adonath left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LauraOlivera commented May 4, 2021

adonath commented May 4, 2021

adonath left a comment

Choose a reason for hiding this comment

adonath commented May 4, 2021

adonath left a comment

Choose a reason for hiding this comment

codecov bot commented Apr 27, 2021 •

edited

Loading

adonath left a comment •

edited

Loading