Implement PSFKernelMap #3689

LauraOlivera · 2021-12-16T11:26:51Z

This pull request addresses issue #3681 to implement a PSFKernelMap.

So far, I have created the object and implemented a from_gauss() method to easily create Gaussian kernel maps. Since one of the uses of this implementation are asymmetric PSF kernels, I added an AsymmetricGauss2DPDF to gammapy.utils.gauss to easily create asymmetric Gaussian kernels.

I marked the PR as draft because I still have a bunch to do. An incomplete list below:

Implement from_geom method
Implement get_psf_kernel() ?
tests for everything

Things that will need to be done in a follow-up PR:

Containment radius methods ?
Orientation of the asymmetric Gaussian ?
Resolve name conflict between PSFKernelMap and PSFKernel.psf_kernel_map

codecov · 2021-12-16T11:48:37Z

Codecov Report

Merging #3689 (978157b) into master (ba7d377) will increase coverage by 0.00%.
The diff coverage is 90.90%.

@@           Coverage Diff            @@
##           master    #3689    +/-   ##
========================================
  Coverage   93.76%   93.76%            
========================================
  Files         162      162            
  Lines       19512    19630   +118     
========================================
+ Hits        18295    18406   +111     
- Misses       1217     1224     +7

Impacted Files	Coverage Δ
gammapy/irf/psf/map.py	`95.89% <90.56%> (-2.02%)`	⬇️
gammapy/utils/gauss.py	`82.29% <92.30%> (+1.56%)`	⬆️
gammapy/makers/background/fov.py	`95.23% <0.00%> (-2.17%)`	⬇️
gammapy/datasets/evaluator.py	`93.36% <0.00%> (-0.03%)`	⬇️
gammapy/datasets/spectrum.py	`99.09% <0.00%> (ø)`
gammapy/datasets/map.py	`93.23% <0.00%> (+<0.01%)`	⬆️
gammapy/makers/background/reflected.py	`98.14% <0.00%> (+0.05%)`	⬆️
gammapy/makers/map.py	`96.07% <0.00%> (+0.20%)`	⬆️
gammapy/makers/utils.py	`98.80% <0.00%> (+0.35%)`	⬆️
gammapy/maps/utils.py	`81.66% <0.00%> (+1.66%)`	⬆️

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 ba7d377...978157b. Read the comment docs.

LauraOlivera · 2022-01-20T19:07:07Z

I feel like I wasn't very creative when making the tests... Please let me know if you can think of something a better than what I implemented!

LauraOlivera · 2022-01-21T14:53:22Z

I made the geom input optional, as we discussed in the dev call.

I also realised I could construct the PSFKernelMap.from_gauss() using PSFKernel.from_spatial_model() using a GaussianSpatialModel with e>0 for the asymmetric case. From a quick test I think that option would be slower, but maybe it is nicer in terms of less code to maintain? It also would guarantee the kernels are normalized, although they get normalized when doing .get_psf_kernel() anyway...
Let me know what you think @adonath

registerrier

Thanks @LauraOlivera .

This is a really useful addition. With the coming of non-isotropic PSF, this type of PSFMap will be needed.

I have the following concerns:

I don't think this class represents a map of kernels but simply a map of non-isotropic PSF. kernel implies a geometry and a projection (and contain a pdf per bin and not per sr). This cannot be described by the lon and lat axes. So this is really a generalized PSFMap. It could even possibly inherit from it and the get_psf_kernel should follow more closely the implementation of PSFMap.get_psf_kernel.
A priori, the stored PSF is in physical coordinates. I assume it should be the same frame as the IRFMap geom, right? This would need to be documented somewhere.

registerrier · 2022-02-11T15:13:25Z

gammapy/irf/psf/map.py

+        LONLAT_AXIS_DEFAULT = MapAxis.from_bounds(-1, 1, nbin=11, unit='deg')
+
+        if not all([psf_lon_axis, psf_lat_axis]):
+            psf_lon_axis = LONLAT_AXIS_DEFAULT.copy()


It is probably better not to change the name attribute directly. You can rather set it during the copy:
LONLAT_AXIS_DEFAULT.copy(name="psf_lon")

registerrier · 2022-02-11T15:15:24Z

gammapy/irf/psf/map.py

+        sigma : `~astropy.coordinates.Angle`
+            Gaussian width.
+        geom : `Geom`
+            Image geometry. By default an allsky geometry is created.


Shouldn't the exposure be argument as well?

registerrier · 2022-02-11T15:17:25Z

gammapy/irf/psf/map.py

+
+        geom_exposure = geom.squash(axis_name='psf_lon').squash(axis_name='psf_lat')
+        exposure_map = Map.from_geom(
+            geom=geom_exposure, unit="m2 s", data=1.0


Is default exposure 0 not more appropriate?

registerrier · 2022-02-11T15:19:39Z

gammapy/irf/psf/map.py

+        sigma : `~astropy.coordinates.Angle`
+            Gaussian width.
+        geom : `Geom`
+            Image geometry. By default an allsky geometry is created.


Maybe specify the dimension of the default geometry?

registerrier · 2022-02-11T15:30:47Z

gammapy/irf/psf/map.py

+
+
+    @classmethod
+    def from_geom(cls, geom):


Since this is really an extension of the PSFMap this probably should follow the logic of the PSFMap.from_geom and create an empty Map. In the long run, if we have anisotropic PSF stored as IRFs we will have to use this class in MapDataset.create

registerrier · 2022-02-11T16:21:05Z

gammapy/irf/psf/map.py

+        else:
+            geom_center = geom.center_skydir
+
+        kernel_geom = WcsGeom.create(skydir=geom_center,


What about the projection to use? Is the default "CAR" OK in general?
Maybe don't allow for empty geom?

registerrier · 2022-02-11T16:22:22Z

gammapy/irf/psf/map.py

+
+        kernel_data = self.psf_kernel_map.to_region_nd_map(region=position)
+
+        width = 2*np.max(np.abs([kernel_data.geom.axes['psf_lon'].edges.min().to_value('deg'),


I would assume that the geometry with is defined from the input geom no?

registerrier · 2022-02-11T16:24:25Z

gammapy/irf/psf/map.py

+
+        position = self._get_nearest_valid_position(position)
+
+        kernel_data = self.psf_kernel_map.to_region_nd_map(region=position)


Shouldn't one compute the lon and lat of each pixel in the input geom and then evaluate the IRFMap at the given coords?

registerrier · 2022-02-11T16:26:06Z

gammapy/irf/psf/map.py

+
+        return cls.from_gauss(energy_axis_true, psf_lon_axis, psf_lat_axis, sigma=0.1 * u.deg, geom=geom.to_image())
+
+    def get_psf_kernel(


I would expect this method to follow much more closely PSFMap.get_psf_kernel except that the map evaluation is performed with lon and lat separation from the center of the target geometry

registerrier · 2022-02-11T16:37:47Z

gammapy/utils/gauss.py

@@ -327,3 +327,49 @@ def gauss_convolve(self, sigma, norm=1):
            norms.append(self.norms[ii] * norm)

        return MultiGauss2D(sigmas, norms)
+
+class AsymmetricGauss2DPDF:


It might prove useful to add a rotation angle no?

adonath · 2022-02-11T18:15:18Z

I don't think this class represents a map of kernels but simply a map of non-isotropic PSF.

@registerrier Your description would imply that the .get_psf_kernel() does the reprojection on the WCS grid, is this correct? I never thought of it this way. I think one huge advantage of the map of kernels is that the computation of the kernel is moved to the data reduction step, which solves one bottleneck of the fitting, where one could just extract the kernel instead of recomputing it, which takes time. The main question for me is whether we would ever need both. A map that really contains kernels (for performance) and one for the asymmetric PSF. My guess is no, because at the data reduction step the geometries are already fixed and don't change after. So one could directly compute the kernel at the data reduction step, just as we do for the EDispKernelMap.

On the other hand this makes the PSFKernelMap really bound to the geometry of e.g. the exposure that is used for the model evaluation, which limits the flexibility.

adonath · 2022-02-11T18:19:11Z

@registerrier Following your understanding of the asymmetric PSF map, we would maybe not need any PSFKernelMap at all, right? Just making sure the current PSFMap does not hard-code the dependency on rad should be sufficient...

adonath · 2022-02-11T18:35:11Z

Just taking the thought a bit further, maybe it still make sense to have both, a PSFKernelMap and a PSFMap like so:

class PSFMap(IRFMap):
    """Implements all common functionality, should be N-dimensional"""
    
    def to_psf_kernel_map(self, geom_kernels):
        """Compute PSFKernelMap"""
        return PSFKernelMap(..., geom_kernels)


# This corresponds to our current PSFMap
class RadiallySymmetricPSFMap(PSFMap):
    """"""
    required_axes = ["rad", "energy_true"]


# This is what Régis thought Laura implemented 
class NonRadiallySymmetricPSFMap(PSFMap):
    """"""
    required_axes = ["psf_lon", "psf_lat", "energy_true"]


# This is what I thought Laura implemented 
class PSFKernelMap(IRFMap):
    """Contains precomputed PSFKernel"""
    def __init__(map, exposure, geom_kernels):
        pass

But I'm just brain-storming for now...

The difference between the PSFKernelMap and the NonRadiallySymmetricPSFMap is the coordinate system the PSF is given in. While for the PSFKernelMap it should be the tangential system to the coordinate system of the exposure Map, given at the required position. For the NonRadiallySymmetricPSFMap it is some spatial coordinate system, independent of the sky position and the kernel is reprojected "on the" fly to sky coordinates.

bkhelifi · 2023-09-26T15:19:03Z

@LauraOlivera your seems to be updated from the main...

LauraOlivera added 4 commits December 9, 2021 19:58

get started

7ad6174

asymmetric gaussian

ea6d56c

asymmetric gauss but well done

ed07229

from gauss

824d2c4

adonath self-assigned this Dec 16, 2021

adonath added the feature label Dec 16, 2021

adonath added this to the 1.0 milestone Dec 16, 2021

adonath added this to To do in gammapy.irf via automation Dec 16, 2021

LauraOlivera added 7 commits December 16, 2021 15:28

from geom

37cf3c0

start get_psf_kernel

3b887a0

get psf kernel

840136d

test from gauss

c5f31c0

test from geom

0c6ff80

test get psf kernel

2e50e8c

typos and docstring

83c5c9b

forgot test gauss

a41a857

LauraOlivera marked this pull request as ready for review January 20, 2022 19:12

LauraOlivera added 2 commits January 21, 2022 14:55

geom optional

6997344

different inputs in test

978157b

registerrier linked an issue Jan 26, 2022 that may be closed by this pull request

Implement PSFKernelMap #3681

Open

registerrier reviewed Feb 11, 2022

View reviewed changes

adonath marked this pull request as draft April 22, 2022 09:00

adonath modified the milestones: 1.0rc, 1.0 Apr 22, 2022

adonath modified the milestones: 1.0, 1.1 May 11, 2022

registerrier modified the milestones: 1.1, 1.2 Apr 24, 2023

registerrier mentioned this pull request Dec 7, 2023

radially asymmetric IRFs #4952

Open

registerrier modified the milestones: 1.2, 1.3 Jan 24, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement PSFKernelMap #3689

Implement PSFKernelMap #3689

LauraOlivera commented Dec 16, 2021 •

edited

codecov bot commented Dec 16, 2021 •

edited

LauraOlivera commented Jan 20, 2022

LauraOlivera commented Jan 21, 2022

registerrier left a comment

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

registerrier Feb 11, 2022

adonath commented Feb 11, 2022 •

edited

adonath commented Feb 11, 2022

adonath commented Feb 11, 2022 •

edited

bkhelifi commented Sep 26, 2023


		kernel_data = self.psf_kernel_map.to_region_nd_map(region=position)

		width = 2*np.max(np.abs([kernel_data.geom.axes['psf_lon'].edges.min().to_value('deg'),


		position = self._get_nearest_valid_position(position)

		kernel_data = self.psf_kernel_map.to_region_nd_map(region=position)


		return cls.from_gauss(energy_axis_true, psf_lon_axis, psf_lat_axis, sigma=0.1 * u.deg, geom=geom.to_image())

		def get_psf_kernel(

Implement PSFKernelMap #3689

Are you sure you want to change the base?

Implement PSFKernelMap #3689

Conversation

LauraOlivera commented Dec 16, 2021 • edited

codecov bot commented Dec 16, 2021 • edited

Codecov Report

LauraOlivera commented Jan 20, 2022

LauraOlivera commented Jan 21, 2022

registerrier left a 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

adonath commented Feb 11, 2022 • edited

adonath commented Feb 11, 2022

adonath commented Feb 11, 2022 • edited

bkhelifi commented Sep 26, 2023

LauraOlivera commented Dec 16, 2021 •

edited

codecov bot commented Dec 16, 2021 •

edited

adonath commented Feb 11, 2022 •

edited

adonath commented Feb 11, 2022 •

edited