New energy axis computation

baccmod/base_acceptance_map_creator.py

@@ -137,12 +137,16 @@ def __init__(self,
         self.exclude_regions = exclude_regions
 
         # Store energy axis computation information
-        self.energy_axis_computation = self.energy_axis if energy_axis_computation is None else energy_axis_computation
+        if energy_axis_computation is None:
+            self.base_energy_axis_computation = self.energy_axis
+        else:
+            self.base_energy_axis_computation = energy_axis_computation
+        self._energy_axis_computation = None
         self.dynamic_energy_axis = dynamic_energy_axis
         self.dynamic_energy_axis_target_statistics = dynamic_energy_axis_target_statistics
         self.dynamic_energy_axis_maximum_wideness_bin = dynamic_energy_axis_maximum_wideness_bin
         self.dynamic_energy_axis_merge_zeros_low_energy = dynamic_energy_axis_merge_zeros_low_energy
-        self. dynamic_energy_axis_merge_zeros_high_energy = dynamic_energy_axis_merge_zeros_high_energy
+        self.dynamic_energy_axis_merge_zeros_high_energy = dynamic_energy_axis_merge_zeros_high_energy
 
         # Calculate map parameter
         self.n_bins_map = 2 * int(np.rint((self.max_offset / spatial_resolution).to(u.dimensionless_unscaled)))
@@ -180,6 +184,13 @@ def __init__(self,
         self.use_mini_irf_computation = use_mini_irf_computation
         self.mini_irf_time_resolution = mini_irf_time_resolution
 
+    @property
+    def energy_axis_computation(self):
+        if self.dynamic_energy_axis:
+            return self._energy_axis_computation
+        else:
+            return self.base_energy_axis_computation
+
     @abstractmethod
     def create_model(self, observations: Observations) -> BackgroundIRF:
         """
@@ -570,29 +581,27 @@ def _compute_time_intervals_based_on_zenith_bin(self, obs: Observation, edge_zen
 
         return time_interval
 
-    def _compute_dynamic_energy_axis(self, base_energy_axis: MapAxis, data: np.ndarray, nb_spatial_bin: int) -> MapAxis:
+    def _compute_dynamic_energy_axis(self, data: np.ndarray, nb_spatial_bin: int):
         """
         Compute a new energy axis from the base one to better have a more uniform number of event in each bin
 
         Parameters
         ----------
-        base_energy_axis: gammapy.maps.MapAxis
-            the base energy axis
         data : np.ndarray
-            the count distribution only binned in energy (using base_energy_axis binning)
+            the count distribution only binned in energy (using base_energy_axis_computation binning)
         nb_spatial_bin : int
             the number of spatial bin covering the data_energy_distribution
 
-        Returns
+        Updates
         -------
-        final_energy_axis: gammapy.maps.MapAxis
+        self.energy_axis_computation: gammapy.maps.MapAxis
             the optimal energy axis
         """
 
         # Relative numerical tolerance for the energy bin width limit
         r_tol = 1 + 1e-7
 
-        edges_energy_axis = base_energy_axis.edges
+        edges_energy_axis = self.base_energy_axis_computation.edges
 
         cumsumdata = np.cumsum(data)/nb_spatial_bin
         rev_cumsumdata = np.cumsum(data[::-1])[::-1]/nb_spatial_bin
@@ -637,13 +646,14 @@ def _compute_dynamic_energy_axis(self, base_energy_axis: MapAxis, data: np.ndarr
         if self.dynamic_energy_axis_merge_zeros_high_energy and len(zeros_highE)>0:
             zeros_highE = np.array([i0], dtype=int)
         indexes_edges=np.sort(np.concatenate([zeros_highE, indexes_edges, zeros_lowE], dtype=int))
-        energy_axis = MapAxis.from_energy_edges(edges_energy_axis[np.array(indexes_edges, dtype=int)], name='energy')
-        logger.log(MOREINFO,'Dynamic energy binning : %s',  np.array_str(np.round(energy_axis.edges, 3)))
+        self._energy_axis_computation = MapAxis.from_energy_edges(
+            edges_energy_axis[np.array(indexes_edges, dtype=int)], name='energy')
+        logger.log(MOREINFO,'Dynamic energy binning : %s',
+                   np.array_str(np.round(self._energy_axis_computation.edges, 3)))
         logger.log(MOREINFO,'Number of bin limited by the maximum bin wideness : %d',  maximum_wideness_hit)
         logger.debug('Number of counts per spatial bin in each energy bin:\n%s', np.array_str(
             np.abs(np.append(np.diff(rev_cumsumdata[np.array(indexes_edges[:-1], dtype=int)]),
                    cumsumdata[-1]+rev_cumsumdata[indexes_edges[-2]]-rev_cumsumdata[0]))))
-        return energy_axis
 
     @staticmethod
     def _split_observations_azimuth(observations: Observations) -> Tuple[Observations, Observations, Dict[int, Dict[str, Any]]]:
@@ -1124,7 +1134,7 @@ def _apply_model_per_observation(self,
             acceptance_map = self._apply_model_all_run(observations=observations, model=model)
         return acceptance_map
 
-    def _interpolate_bkg_to_energy_axis(self, data_bkg: u.Quantity, energy_axis_computation: MapAxis):
+    def _interpolate_bkg_to_energy_axis(self, data_bkg: u.Quantity):
         """
             Compute the final background model from the provided data by interpolating the energy axis used for
             computation to the ones for the model
@@ -1133,8 +1143,6 @@ def _interpolate_bkg_to_energy_axis(self, data_bkg: u.Quantity, energy_axis_comp
             ----------
             data_bkg : u.Quantity
                 The data cube of the background model. The energy axis needs to be the first one.
-            energy_axis_computation : gammapy.maps.geom.MapAxis
-                The energy axis used for computation
 
             Returns
             -------
@@ -1143,7 +1151,8 @@ def _interpolate_bkg_to_energy_axis(self, data_bkg: u.Quantity, energy_axis_comp
         """
 
         # Return the provided bkg data if energy axis are matching
-        if energy_axis_computation.nbin == self.energy_axis.nbin and np.all(energy_axis_computation.edges == self.energy_axis.edges):
+        if self.energy_axis_computation.nbin == self.energy_axis.nbin and np.all(
+                self.energy_axis_computation.edges == self.energy_axis.edges):
             logger.info('Identical computation energy axis and model energy axis, no interpolation required')
             return data_bkg
 
@@ -1155,7 +1164,7 @@ def _interpolate_bkg_to_energy_axis(self, data_bkg: u.Quantity, energy_axis_comp
         min_value = np.min(data_bkg.value[~mask_zero_input])
         max_value = np.max(data_bkg.value[~mask_zero_input])
 
-        interp_func = interp1d(x=np.log10(energy_axis_computation.center.to_value(u.TeV)),
+        interp_func = interp1d(x=np.log10(self.energy_axis_computation.center.to_value(u.TeV)),
                                y=raw_log_data,
                                axis=0,
                                fill_value='extrapolate')

baccmod/fit_acceptance_map_creator.py

@@ -83,15 +83,15 @@ def create_model(self, observations) -> Background3D:
         (count_background,
          exp_map_background,
          exp_map_background_total,
-         livetime, energy_axis_computation) = self._create_base_computation_map(observations)
+         livetime) = self._create_base_computation_map(observations)
 
         # 2) downsample exposures
         exp_ds = exp_map_background.downsample(self.oversample_map, preserve_counts=True)
         exp_total_ds = exp_map_background_total.downsample(self.oversample_map, preserve_counts=True)
 
         # 3) Get spatial axis for the bkg model and the bins size
         extended_offset_axis_x, extended_offset_axis_y, solid_angle = self._get_ext_axis_and_solid_angle()
-        energy_bin_width = energy_axis_computation.bin_width
+        energy_bin_width = self.energy_axis_computation.bin_width
 
         # 4) fit function on counts → “corrected counts”
         predicted_counts = np.empty(count_background.shape)
@@ -103,7 +103,7 @@ def create_model(self, observations) -> Background3D:
         centers = self.offset_axis.center.to_value('deg')
         centers = np.concatenate((-np.flip(centers), centers), axis=None)
         if self.model_to_fit.n_inputs != 2:
-            coordinates.append(energy_axis_computation.center.to_value('TeV'))
+            coordinates.append(self.energy_axis_computation.center.to_value('TeV'))
         if self.model_to_fit.n_inputs > 1:
             coordinates.append(centers)
             coordinates.append(centers)
@@ -126,13 +126,13 @@ def create_model(self, observations) -> Background3D:
         elif self.model_to_fit.n_inputs == 2:
             bin_size = solid_angle.value
             logger.info(
-                "Fitting background per enery bin"
+                "Fitting background per energy bin"
                 )
             for e in range(count_background.shape[0]):
                 logger.log(MOREINFO,
                     "Fitting energy bin [%.2f, %.2f] TeV",
-                    energy_axis_computation.edges[e].to_value('TeV'),
-                    energy_axis_computation.edges[e + 1].to_value('TeV')
+                    self.energy_axis_computation.edges[e].to_value('TeV'),
+                    self.energy_axis_computation.edges[e + 1].to_value('TeV')
                 )
                 predicted_counts[e] = self._fit_background(
                     self.model_to_fit,
@@ -164,7 +164,7 @@ def create_model(self, observations) -> Background3D:
             raise RuntimeError(f"The provided model dimension is incorrect : {self.model_to_fit.n_inputs}")
 
         logger.info(
-            "Average event sqrt‐residuals per energy: %s, std = %s",
+            "Average event sqrt‐residuals per fit: %s, std = %s",
             np.array_str(np.round(self.sq_rel_residuals['mean'], 2)),
             np.array_str(np.round(self.sq_rel_residuals['std'], 2))
         )
@@ -173,10 +173,10 @@ def create_model(self, observations) -> Background3D:
         data_background = (
             predicted_counts
             / solid_angle[np.newaxis, :, :]
-            / energy_axis_computation.bin_width[:, np.newaxis, np.newaxis]
+            / self.energy_axis_computation.bin_width[:, np.newaxis, np.newaxis]
             / livetime
         )
-        data_background = self._interpolate_bkg_to_energy_axis(data_background, energy_axis_computation)
+        data_background = self._interpolate_bkg_to_energy_axis(data_background)
 
         # 6) instantiate Background3D
         acceptance_map = Background3D(axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y],

baccmod/fitting.py

@@ -144,6 +144,5 @@ def log_factorial(x: np.ndarray) -> np.ndarray:
 def log_poisson(mu: np.ndarray, x: np.ndarray, log_factorial_x: np.ndarray) -> np.ndarray:
     """ Poisson pdf in log scale. """
     if np.any(mu <= 0):
-        logger.warning('log poisson received a zero or negative value.')
         mu[mu <= 0] = np.finfo(mu.dtype).tiny
     return -mu + x * np.log(mu) - log_factorial_x

baccmod/grid3d_acceptance_map_creator.py

@@ -109,7 +109,7 @@ def create_model(self, observations: Observations) -> Background3D:
         """
 
         # Compute base data
-        count_background, exp_map_background, exp_map_background_total, livetime, energy_axis_computation = self._create_base_computation_map(
+        count_background, exp_map_background, exp_map_background_total, livetime = self._create_base_computation_map(
             observations)
 
         # Downsample map to background model resolution
@@ -125,10 +125,10 @@ def create_model(self, observations: Observations) -> Background3D:
                                                exp_map_background_downsample.data)
         data_background = (corrected_counts /
                            solid_angle[np.newaxis, :, :] /
-                           energy_axis_computation.bin_width[:, np.newaxis, np.newaxis] /
+                           self.energy_axis_computation.bin_width[:, np.newaxis, np.newaxis] /
                            livetime)
 
-        data_background = self._interpolate_bkg_to_energy_axis(data_background, energy_axis_computation)
+        data_background = self._interpolate_bkg_to_energy_axis(data_background)
 
         acceptance_map = Background3D(axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y],
                                       data=data_background.to(u.Unit('s-1 MeV-1 sr-1')),
@@ -163,21 +163,18 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[np.n
         offset_edges = self.offset_axis.edges
         offset_bins = np.round(np.concatenate((-np.flip(offset_edges), offset_edges[1:]), axis=None), 3)
 
-        energy_axis_computation = self.energy_axis_computation.copy()
         if self.dynamic_energy_axis:
-            data_energy_distribution = np.zeros(self.energy_axis_computation.nbin, dtype=np.int64)
+            data_energy_distribution = np.zeros(self.base_energy_axis_computation.nbin, dtype=np.int64)
             for obs in observations:
                 events_camera_frame = self._get_events_in_camera_frame(obs)
                 mask_event = np.logical_and(np.abs(events_camera_frame.lon) <= self.max_offset,
                                             np.abs(events_camera_frame.lat) <= self.max_offset)
-                distrib, _ = np.histogram(obs.events.energy[mask_event], energy_axis_computation.edges)
+                distrib, _ = np.histogram(obs.events.energy[mask_event], self.base_energy_axis_computation.edges)
                 data_energy_distribution += distrib
-            energy_axis_computation = self._compute_dynamic_energy_axis(energy_axis_computation,
-                                                                        data_energy_distribution,
-                                                                        len(offset_bins)**2)
+            self._compute_dynamic_energy_axis(data_energy_distribution, len(offset_bins)**2)
 
-        map_bins = (energy_axis_computation.edges, offset_bins, offset_bins)
-        geom = self._get_geom(energy_axis_computation)
+        map_bins = (self.energy_axis_computation.edges, offset_bins, offset_bins)
+        geom = self._get_geom(self.energy_axis_computation)
         count_background = np.zeros((len(map_bins[0]) - 1,
                                      len(map_bins[1]) - 1,
                                      len(map_bins[2]) - 1))
@@ -243,4 +240,4 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[np.n
                 exp_map_background_total.data += exp_map_obs_total.counts.data
                 livetime += obs.observation_live_time_duration
 
-        return count_background, exp_map_background, exp_map_background_total, livetime, energy_axis_computation
+        return count_background, exp_map_background, exp_map_background_total, livetime

baccmod/radial_acceptance_map_creator.py

@@ -78,11 +78,14 @@ def create_model(self, observations: Observations) -> Background2D:
         -------
         acceptance_map : Background2D
         """
-        count_map_background, exp_map_background, exp_map_background_total, livetime, energy_axis_computation = self._create_base_computation_map(
-            observations)
-        geom = self._get_geom(energy_axis_computation)
-
-        data_background = np.zeros((energy_axis_computation.nbin, self.offset_axis.nbin)) * u.Unit('s-1 MeV-1 sr-1')
+        (count_map_background,
+         exp_map_background,
+         exp_map_background_total,
+         livetime) = self._create_base_computation_map(observations)
+        geom = self._get_geom(self.energy_axis_computation)
+
+        data_background = np.zeros((self.energy_axis_computation.nbin, self.offset_axis.nbin)
+                                   ) * u.Unit('s-1 MeV-1 sr-1')
         for i in range(self.offset_axis.nbin):
             if np.isclose(0. * u.deg, self.offset_axis.edges[i]):
                 selection_region = CircleSkyRegion(center=self.center_map, radius=self.offset_axis.edges[i + 1])
@@ -91,18 +94,19 @@ def create_model(self, observations: Observations) -> Background2D:
                                                           inner_radius=self.offset_axis.edges[i],
                                                           outer_radius=self.offset_axis.edges[i + 1])
             selection_map = geom.to_image().region_mask([selection_region])
-            for j in range(energy_axis_computation.nbin):
+            for j in range(self.energy_axis_computation.nbin):
                 value = u.dimensionless_unscaled * np.sum(count_map_background.data[j, :, :] * selection_map)
                 value *= np.sum(exp_map_background_total.data[j, :, :] * selection_map) / np.sum(
                     exp_map_background.data[j, :, :] * selection_map)
 
-                value /= (energy_axis_computation.edges[j + 1] - energy_axis_computation.edges[j])
+                value /= (self.energy_axis_computation.edges[j + 1] - self.energy_axis_computation.edges[j])
                 value /= 2. * np.pi * (
                             np.cos(self.offset_axis.edges[i]) - np.cos(self.offset_axis.edges[i + 1])) * u.steradian
                 value /= livetime
                 data_background[j, i] = value
 
-        acceptance_map = Background2D(axes=[self.energy_axis, self.offset_axis], data=self._interpolate_bkg_to_energy_axis(data_background, energy_axis_computation))
+        acceptance_map = Background2D(axes=[self.energy_axis, self.offset_axis],
+                                      data=self._interpolate_bkg_to_energy_axis(data_background))
 
         return acceptance_map
 
@@ -130,24 +134,23 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[np.n
             The energy axis used for the computation
         """
 
-        energy_axis_computation = self.energy_axis_computation.copy()
         if self.dynamic_energy_axis:
-            data_energy_distribution = np.zeros(self.energy_axis_computation.nbin, dtype=np.int64)
+            data_energy_distribution = np.zeros(self.base_energy_axis_computation.nbin, dtype=np.int64)
             for obs in observations:
                 mask_event = obs.events.offset <= self.max_offset
-                distrib, _ = np.histogram(obs.events.energy[mask_event], energy_axis_computation.edges)
+                distrib, _ = np.histogram(obs.events.energy[mask_event], self.base_energy_axis_computation.edges)
                 data_energy_distribution += distrib
-            energy_axis_computation = self._compute_dynamic_energy_axis(energy_axis_computation, data_energy_distribution, self.offset_axis.nbin)
+            self._compute_dynamic_energy_axis(data_energy_distribution, self.offset_axis.nbin)
 
-        geom = self._get_geom(energy_axis_computation)
+        geom = self._get_geom(self.energy_axis_computation)
         count_map_background = WcsNDMap(geom=geom)
         exp_map_background = WcsNDMap(geom=geom, unit=u.s)
         exp_map_background_total = WcsNDMap(geom=geom, unit=u.s)
         livetime = 0. * u.s
 
         for obs in observations:
             geom = WcsGeom.create(skydir=obs.pointing.fixed_icrs, npix=(self.n_bins_map, self.n_bins_map),
-                                  binsz=self.spatial_bin_size, frame="icrs", axes=[energy_axis_computation])
+                                  binsz=self.spatial_bin_size, frame="icrs", axes=[self.energy_axis_computation])
             count_map_obs, exclusion_mask = self._create_map(obs, geom, self.exclude_regions, add_bkg=False)
 
             exp_map_obs = MapDataset.create(geom=count_map_obs.geoms['geom'])
@@ -164,4 +167,4 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[np.n
             exp_map_background_total.data += exp_map_obs_total.counts.data
             livetime += obs.observation_live_time_duration
 
-        return count_map_background, exp_map_background, exp_map_background_total, livetime, energy_axis_computation
+        return count_map_background, exp_map_background, exp_map_background_total, livetime