Skip to content

Integration in the fitting procedure #87

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
gabemery merged 2 commits into from
Jun 18, 2026
+46 −34
Merged

Integration in the fitting procedure #87

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f4b2f36
Implement center point integration in the fitting procedure. Fix supp…
gabemery Apr 28, 2026
2968298
Refactor spatial axis manipulation in Grid3DAcceptanceMapCreator
gabemery Jun 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion baccmod/base_acceptance_map_creator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1135,7 +1135,7 @@ 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 len(energy_axis_computation.edges) == len(self.energy_axis.edges) and np.all(energy_axis_computation.edges == self.energy_axis.edges):
if energy_axis_computation.nbin == self.energy_axis.nbin and np.all(energy_axis_computation.edges == self.energy_axis.edges):
logger.info('Identical computation energy axis and model energy axis, no interpolation required')
return data_bkg

Expand Down
45 changes: 23 additions & 22 deletions baccmod/fit_acceptance_map_creator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ class FitAcceptanceMapCreator(Grid3DAcceptanceMapCreator):

def __init__(
self,
energy_axis: MapAxis,
offset_axis: MapAxis,
model_to_fit: FittableModel = Gaussian2D(),
list_name_normalisation_parameter: List[str] = ['amplitude'],
maxiter: int = 1000,
Expand All @@ -56,12 +58,10 @@ def __init__(
"""

# Call the “stack”‐only constructor in Grid3D, to set up geometry, offset axes, etc.
super().__init__(**kwargs)
super().__init__(energy_axis, offset_axis, **kwargs)

self.sq_rel_residuals = {"mean": [], "std": []}
self.model_to_fit = model_to_fit
if model_to_fit.n_inputs !=2: #TODO remove once 1D and 3D models are tested
Comment thread
mdebony marked this conversation as resolved.
logger.warning('Only 2D models, fitted per energy bin, were thoroughly tested.')
self.list_name_normalisation_parameter = list_name_normalisation_parameter
self.maxiter = maxiter

Expand Down Expand Up @@ -89,8 +89,11 @@ def create_model(self, observations) -> Background3D:
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

# 3) fit function on counts → “corrected counts”
# 4) fit function on counts → “corrected counts”
predicted_counts = np.empty(count_background.shape)
self.sq_rel_residuals = {"mean": [], "std": []}

Expand All @@ -108,6 +111,7 @@ def create_model(self, observations) -> Background3D:

# perform the fit, looping over non fitted axes
if self.model_to_fit.n_inputs == 3:
bin_size = (solid_angle[np.newaxis, :, :] * energy_bin_width[:, np.newaxis, np.newaxis]).value
logger.info(
"Fitting background with a 3D model."
)
Expand All @@ -117,8 +121,10 @@ def create_model(self, observations) -> Background3D:
count_map=count_background.astype(int),
exp_map_total=exp_total_ds.data,
exp_map=exp_ds.data,
bin_size=bin_size
)
elif self.model_to_fit.n_inputs == 2:
bin_size = solid_angle.value
logger.info(
"Fitting background per enery bin"
)
Expand All @@ -134,8 +140,10 @@ def create_model(self, observations) -> Background3D:
count_map=count_background[e].astype(int),
exp_map_total=exp_total_ds.data[e],
exp_map=exp_ds.data[e],
bin_size=bin_size
)
elif self.model_to_fit.n_inputs == 1:
bin_size = energy_bin_width.value
logger.info(
"Fitting background per spatial bin"
)
Expand All @@ -150,6 +158,7 @@ def create_model(self, observations) -> Background3D:
count_map=count_background[:,x,y].astype(int),
exp_map_total=exp_total_ds.data[:,x,y],
exp_map=exp_ds.data[:,x,y],
bin_size=bin_size
)
else:
raise RuntimeError(f"The provided model dimension is incorrect : {self.model_to_fit.n_inputs}")
Expand All @@ -160,16 +169,6 @@ def create_model(self, observations) -> Background3D:
np.array_str(np.round(self.sq_rel_residuals['std'], 2))
)

# 4) build final offset axes (matching Grid3DAcceptanceMapCreator)
edges = self.offset_axis.edges
extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None)
extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon")
extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat")

bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1)
bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0)
solid_angle = 4.0 * (np.sin(bin_width_x / 2) * np.sin(bin_width_y / 2)) * u.steradian

# 5) normalize to flux units
data_background = (
predicted_counts
Expand All @@ -186,7 +185,7 @@ def create_model(self, observations) -> Background3D:
return acceptance_map

def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarray,
exp_map_total: np.ndarray, exp_map: np.ndarray) -> np.ndarray:
exp_map_total: np.ndarray, exp_map: np.ndarray, bin_size: np.ndarray) -> np.ndarray:
"""
Perform a Poisson fit on the given map (could be 1D, 2D or 3D), given
the fine‐binned exposure (exp_map) and total‐exposure (exp_map_total),
Expand All @@ -199,16 +198,18 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr
coords : np.ndarray
the list of coordinates
count_map : np.ndarray
counts in each pixel, integer
counts in each bin, integer
exp_map_total : np.ndarray
exposure WITHOUT exclusion regions
exp_map : np.ndarray
exposure CORRECTED for exclusion regions
bin_size : np.ndarray
integral size of each bin

Returns
-------
fitted_counts_model : 2D np.ndarray
The best‐fit model of counts (on the same pixel grid as count_map).
fitted_counts_model : np.ndarray
The best‐fit model of counts (on the same bin grid as count_map).
"""
total_counts = np.sum(count_map)
# Skip the fit if there are no events
Expand All @@ -230,7 +231,7 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr
# Correct normalisation of the model
if self.list_name_normalisation_parameter is not None and len(self.list_name_normalisation_parameter) > 0:
# Compute correction
init_count_model = np.sum(model_init(*coords)*exp_correction)
init_count_model = np.sum(model_init(*coords)*exp_correction*bin_size)
correction_norm = total_counts/init_count_model

# build list of free parameters
Expand All @@ -248,10 +249,10 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr

# Fit the model
best_model = poisson_fitter(model_init, *coords, data=count_map, exposure_correction=exp_correction,
mask=mask, maxiter=self.maxiter)
mask=mask, maxiter=self.maxiter, bin_size=bin_size)
# Collect results & (optionally) track residuals
if logger.getEffectiveLevel() <= logging.INFO:
fitted_model = best_model(*coords) * exp_correction
fitted_model = best_model(*coords) * exp_correction * bin_size
fitted_model[exp_map == 0] = 1.0
rel_resid = 100 * (count_map - fitted_model) / fitted_model
sq_rel_resid = (count_map - fitted_model) / np.sqrt(fitted_model)
Expand All @@ -263,4 +264,4 @@ def _fit_background(self, model: Model, *coords: np.ndarray, count_map: np.ndarr
"Avg rel residual: %.1f, Std = %.2f\n", np.mean(rel_resid), np.std(rel_resid)
)

return best_model(*coords)
return best_model(*coords) * bin_size
10 changes: 7 additions & 3 deletions baccmod/fitting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ def poisson_fitter(model: Model,
data: np.ndarray,
exposure_correction: np.ndarray = None,
mask: np.ndarray = None,
maxiter: int = 1000) -> Model:
maxiter: int = 1000,
bin_size: np.ndarray = 1) -> Model:
"""
Fit the model to the data following poisson likelihood statistics and using iminuit

Expand All @@ -40,6 +41,8 @@ def poisson_fitter(model: Model,
maxiter : int
maximum number of iteration for fitting, as the fitting is performed in two steps,
could be the double of this value in practice
bin_size : np.ndarray
integral size of each bin

Returns
-------
Expand Down Expand Up @@ -88,7 +91,7 @@ def apply_params_and_tied(pars):
# negative log‑likelihood
def neg_logL(**pars):
apply_params_and_tied(pars)
mu = model_copy(*coords) * exposure_correction
mu = model_copy(*coords) * exposure_correction * bin_size
return -np.sum(log_poisson(mu[mask], data[mask], log_fact[mask]))

# wrapper to accept either positional arguments or keyword arguments (exclusives)
Expand Down Expand Up @@ -139,7 +142,8 @@ 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.min(mu[mu > 0])/2
mu[mu <= 0] = np.finfo(mu.dtype).tiny
return -mu + x * np.log(mu) - log_factorial_x
23 changes: 15 additions & 8 deletions baccmod/grid3d_acceptance_map_creator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,19 @@ def __init__(self,
exclude_regions=exclude_regions,
**kwargs)

def _get_ext_axis_and_solid_angle(self):
""" Compute the solid angle of spatial bins """
edges = self.offset_axis.edges
extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None)
extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon")
extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat")

bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1)
bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0)
solid_angle = 4.0 * (np.sin(bin_width_x / 2) * np.sin(bin_width_y / 2)) * u.steradian

return extended_offset_axis_x, extended_offset_axis_y, solid_angle

def create_model(self, observations: Observations) -> Background3D:
"""
Calculate a 3D grid acceptance map
Expand All @@ -104,18 +117,12 @@ def create_model(self, observations: Observations) -> Background3D:
exp_map_background_total_downsample = exp_map_background_total.downsample(self.oversample_map,
preserve_counts=True)

# Create axis for bkg model
edges = self.offset_axis.edges
extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None)
extended_offset_axis_x = MapAxis.from_edges(extended_edges, name='fov_lon')
bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1)
extended_offset_axis_y = MapAxis.from_edges(extended_edges, name='fov_lat')
bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0)
# Create spatial axis for the bkg model and the bin solid angle
extended_offset_axis_x, extended_offset_axis_y, solid_angle = self._get_ext_axis_and_solid_angle()

# Compute acceptance_map
corrected_counts = count_background * (exp_map_background_total_downsample.data /
exp_map_background_downsample.data)
solid_angle = 4. * (np.sin(bin_width_x / 2.) * np.sin(bin_width_y / 2.)) * u.steradian
data_background = (corrected_counts /
solid_angle[np.newaxis, :, :] /
energy_axis_computation.bin_width[:, np.newaxis, np.newaxis] /
Expand Down
Loading