From d6d6dab644b117e355f7a3c766c6cbc00c9b8c18 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 6 Dec 2023 14:35:49 +0100 Subject: [PATCH 1/5] avoid saving exposure cutout in evaluator Signed-off-by: Quentin Remy --- gammapy/datasets/evaluator.py | 34 +++++++++++++++++++++++----------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/gammapy/datasets/evaluator.py b/gammapy/datasets/evaluator.py index 184d698aa3..408aeac9c8 100644 --- a/gammapy/datasets/evaluator.py +++ b/gammapy/datasets/evaluator.py @@ -63,7 +63,7 @@ def __init__( ): self.model = model - self.exposure = exposure + self._exposure = exposure self.psf = psf self.edisp = edisp self.mask = mask @@ -72,6 +72,7 @@ def __init__( self._init_position = None self.contributes = True self.psf_containment = None + self._geom = None if evaluation_mode not in {"local", "global"}: raise ValueError(f"Invalid evaluation_mode: {evaluation_mode!r}") @@ -93,7 +94,8 @@ def __init__( self._cached_position = (0, 0) self._computation_cache = None self._spatial_oversampling_factor = 1 - if self.exposure is not None: + if self._exposure is not None: + self._geom = self._exposure.geom if not self.geom.is_region or self.geom.region is not None: self.update_spatial_oversampling_factor(self.geom) @@ -111,7 +113,7 @@ def reset_cache_properties(self): @property def geom(self): """True energy map geometry (`~gammapy.maps.Geom`).""" - return self.exposure.geom + return self._geom @property def _geom_reco(self): @@ -130,7 +132,7 @@ def needs_update(self): return False elif not self.contributes: return False - elif self.exposure is None: + elif self._exposure is None: return True elif self.geom.is_region: return False @@ -216,13 +218,10 @@ def update(self, exposure, psf, edisp, geom, mask): if self.evaluation_mode == "local": self.contributes = self.model.contributes(mask=mask, margin=self.psf_width) - - if self.contributes: - self.exposure = exposure.cutout( - position=self.model.position, width=self.cutout_width, odd_npix=True - ) - else: - self.exposure = exposure + self._geom = exposure.geom.cutout( + position=self.model.position, width=self.cutout_width, odd_npix=True + ) + self._exposure = exposure if self.contributes: if not self.geom.is_region or self.geom.region is not None: @@ -239,6 +238,19 @@ def _edisp_diagonal(self): energy_axis=self.edisp.axes["energy"], ) + @property + def exposure(self): + if ( + self._exposure is not None + and self.contributes + and self.model.evaluation_radius is not None + ): + return self._exposure.cutout( + position=self.model.position, width=self.cutout_width, odd_npix=True + ) + else: + return self._exposure + def update_spatial_oversampling_factor(self, geom): """Update spatial oversampling_factor for model evaluation.""" res_scale = self.model.evaluation_bin_size_min From 681a71a906f03e64d80cae02b0736296ca259fbe Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Thu, 14 Dec 2023 23:32:38 +0100 Subject: [PATCH 2/5] fix Signed-off-by: Quentin Remy --- gammapy/datasets/evaluator.py | 43 ++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/gammapy/datasets/evaluator.py b/gammapy/datasets/evaluator.py index 408aeac9c8..0c298b96a5 100644 --- a/gammapy/datasets/evaluator.py +++ b/gammapy/datasets/evaluator.py @@ -72,7 +72,6 @@ def __init__( self._init_position = None self.contributes = True self.psf_containment = None - self._geom = None if evaluation_mode not in {"local", "global"}: raise ValueError(f"Invalid evaluation_mode: {evaluation_mode!r}") @@ -94,10 +93,11 @@ def __init__( self._cached_position = (0, 0) self._computation_cache = None self._spatial_oversampling_factor = 1 - if self._exposure is not None: - self._geom = self._exposure.geom - if not self.geom.is_region or self.geom.region is not None: - self.update_spatial_oversampling_factor(self.geom) + if exposure is not None: + self._geom = exposure.geom + self.update_spatial_oversampling_factor(self.geom) + else: + self._geom = None def _repr_html_(self): try: @@ -216,16 +216,15 @@ def update(self, exposure, psf, edisp, geom, mask): max_radius=PSF_MAX_RADIUS, ) + self._exposure = exposure + self._geom = exposure.geom if self.evaluation_mode == "local": self.contributes = self.model.contributes(mask=mask, margin=self.psf_width) - self._geom = exposure.geom.cutout( - position=self.model.position, width=self.cutout_width, odd_npix=True - ) - self._exposure = exposure - - if self.contributes: - if not self.geom.is_region or self.geom.region is not None: - self.update_spatial_oversampling_factor(self.geom) + if self.contributes and not self.geom.is_region: + self._geom = exposure.geom.cutout( + position=self.model.position, width=self.cutout_width, odd_npix=True + ) + self.update_spatial_oversampling_factor(self.geom) self.reset_cache_properties() self._computation_cache = None @@ -242,8 +241,8 @@ def _edisp_diagonal(self): def exposure(self): if ( self._exposure is not None + and self.evaluation_mode == "local" and self.contributes - and self.model.evaluation_radius is not None ): return self._exposure.cutout( position=self.model.position, width=self.cutout_width, odd_npix=True @@ -253,15 +252,17 @@ def exposure(self): def update_spatial_oversampling_factor(self, geom): """Update spatial oversampling_factor for model evaluation.""" - res_scale = self.model.evaluation_bin_size_min - res_scale = res_scale.to_value("deg") if res_scale is not None else 0 + if self.contributes and (not geom.is_region or geom.region is not None): + res_scale = self.model.evaluation_bin_size_min + + res_scale = res_scale.to_value("deg") if res_scale is not None else 0 - if res_scale != 0: - if geom.is_region or geom.is_hpx: - geom = geom.to_wcs_geom() - factor = int(np.ceil(np.max(geom.pixel_scales.deg) / res_scale)) - self._spatial_oversampling_factor = factor + if res_scale != 0: + if geom.is_region or geom.is_hpx: + geom = geom.to_wcs_geom() + factor = int(np.ceil(np.max(geom.pixel_scales.deg) / res_scale)) + self._spatial_oversampling_factor = factor def compute_dnde(self): """Compute model differential flux at map pixel centers. From f779d07bb11ae943d2efde86a10f3548881c40bd Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Fri, 15 Dec 2023 12:22:33 +0100 Subject: [PATCH 3/5] lazy cutout width and position Signed-off-by: Quentin Remy --- gammapy/datasets/evaluator.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/gammapy/datasets/evaluator.py b/gammapy/datasets/evaluator.py index 0c298b96a5..fdc334518d 100644 --- a/gammapy/datasets/evaluator.py +++ b/gammapy/datasets/evaluator.py @@ -69,7 +69,6 @@ def __init__( self.mask = mask self.gti = gti self.use_cache = use_cache - self._init_position = None self.contributes = True self.psf_containment = None @@ -161,7 +160,12 @@ def use_psf_containment(self, geom): is_circle_region = isinstance(geom.region, CircleSkyRegion) return is_point_model & is_circle_region - @property + @lazyproperty + def position(self): + """Latest evaluation position.""" + return self.model.position + + @lazyproperty def cutout_width(self): """Cutout width for the model component.""" return self.psf_width + 2 * (self.model.evaluation_radius + CUTOUT_MARGIN) @@ -185,11 +189,14 @@ def update(self, exposure, psf, edisp, geom, mask): # TODO: simplify and clean up log.debug("Updating model evaluator") + del self.position + del self.cutout_width + # lookup edisp if edisp: energy_axis = geom.axes["energy"] self.edisp = edisp.get_edisp_kernel( - position=self.model.position, energy_axis=energy_axis + position=self.position, energy_axis=energy_axis ) del self._edisp_diagonal @@ -210,7 +217,7 @@ def update(self, exposure, psf, edisp, geom, mask): geom_psf = geom_psf.to_wcs_geom() self.psf = psf.get_psf_kernel( - position=self.model.position, + position=self.position, geom=geom_psf, containment=PSF_CONTAINMENT, max_radius=PSF_MAX_RADIUS, @@ -222,7 +229,7 @@ def update(self, exposure, psf, edisp, geom, mask): self.contributes = self.model.contributes(mask=mask, margin=self.psf_width) if self.contributes and not self.geom.is_region: self._geom = exposure.geom.cutout( - position=self.model.position, width=self.cutout_width, odd_npix=True + position=self.position, width=self.cutout_width, odd_npix=True ) self.update_spatial_oversampling_factor(self.geom) @@ -245,7 +252,7 @@ def exposure(self): and self.contributes ): return self._exposure.cutout( - position=self.model.position, width=self.cutout_width, odd_npix=True + position=self.position, width=self.cutout_width, odd_npix=True ) else: return self._exposure From 769a4acb9866fb0818da8a9eb8fa3e0c6401d44b Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Mon, 1 Jan 2024 16:24:14 +0100 Subject: [PATCH 4/5] cutout_view Signed-off-by: Quentin Remy --- gammapy/maps/wcs/ndmap.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/gammapy/maps/wcs/ndmap.py b/gammapy/maps/wcs/ndmap.py index 25323bc61c..4ff2051a93 100644 --- a/gammapy/maps/wcs/ndmap.py +++ b/gammapy/maps/wcs/ndmap.py @@ -994,6 +994,38 @@ def cutout(self, position, width, mode="trim", odd_npix=False): return self._init_copy(geom=geom_cutout, data=data) + def _cutout_view(self, position, width, odd_npix=False): + """ + Create a cutout around a given position without copy of the data. + + Parameters + ---------- + position : `~astropy.coordinates.SkyCoord` + Center position of the cutout region. + width : tuple of `~astropy.coordinates.Angle` + Angular sizes of the region in (lon, lat) in that specific order. + If only one value is passed, a square region is extracted. + odd_npix : bool, optional + Force width to odd number of pixels. + Default is False. + + Returns + ------- + cutout : `~gammapy.maps.WcsNDMap` + Cutout map. + """ + geom_cutout = self.geom.cutout( + position=position, width=width, mode="trim", odd_npix=odd_npix + ) + cutout_info = geom_cutout.cutout_slices(self.geom, mode="trim") + + slices = cutout_info["parent-slices"] + parent_slices = Ellipsis, slices[0], slices[1] + + return self.__class__.from_geom( + geom=geom_cutout, data=self.quantity[parent_slices] + ) + def stack(self, other, weights=None, nan_to_num=True): """Stack cutout into map. From 8f708ed4ce9a15b697668c5694a5d8e29bc888e8 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Mon, 1 Jan 2024 16:26:53 +0100 Subject: [PATCH 5/5] exposure view Signed-off-by: Quentin Remy --- gammapy/datasets/evaluator.py | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/gammapy/datasets/evaluator.py b/gammapy/datasets/evaluator.py index fdc334518d..b8416fa901 100644 --- a/gammapy/datasets/evaluator.py +++ b/gammapy/datasets/evaluator.py @@ -63,7 +63,7 @@ def __init__( ): self.model = model - self._exposure = exposure + self.exposure = exposure self.psf = psf self.edisp = edisp self.mask = mask @@ -93,10 +93,7 @@ def __init__( self._computation_cache = None self._spatial_oversampling_factor = 1 if exposure is not None: - self._geom = exposure.geom self.update_spatial_oversampling_factor(self.geom) - else: - self._geom = None def _repr_html_(self): try: @@ -112,7 +109,10 @@ def reset_cache_properties(self): @property def geom(self): """True energy map geometry (`~gammapy.maps.Geom`).""" - return self._geom + if self.exposure is not None: + return self.exposure.geom + else: + return None @property def _geom_reco(self): @@ -131,7 +131,7 @@ def needs_update(self): return False elif not self.contributes: return False - elif self._exposure is None: + elif self.exposure is None: return True elif self.geom.is_region: return False @@ -223,12 +223,11 @@ def update(self, exposure, psf, edisp, geom, mask): max_radius=PSF_MAX_RADIUS, ) - self._exposure = exposure - self._geom = exposure.geom + self.exposure = exposure if self.evaluation_mode == "local": self.contributes = self.model.contributes(mask=mask, margin=self.psf_width) if self.contributes and not self.geom.is_region: - self._geom = exposure.geom.cutout( + self.exposure = exposure._cutout_view( position=self.position, width=self.cutout_width, odd_npix=True ) self.update_spatial_oversampling_factor(self.geom) @@ -244,19 +243,6 @@ def _edisp_diagonal(self): energy_axis=self.edisp.axes["energy"], ) - @property - def exposure(self): - if ( - self._exposure is not None - and self.evaluation_mode == "local" - and self.contributes - ): - return self._exposure.cutout( - position=self.position, width=self.cutout_width, odd_npix=True - ) - else: - return self._exposure - def update_spatial_oversampling_factor(self, geom): """Update spatial oversampling_factor for model evaluation."""