From 41386a3787e41be6218d1663b8bf09c6d8485be7 Mon Sep 17 00:00:00 2001
From: Ramiro <rtorrese@proton.me>
Date: Fri, 10 Jan 2025 21:01:27 +0800
Subject: [PATCH 1/7] adding files from the master upstream

---
 pyproject.toml | 101 +++++++++++++++++++++++++++++++++++++++++++++++++
 ruff.toml      |  15 ++++++++
 2 files changed, 116 insertions(+)
 create mode 100644 pyproject.toml
 create mode 100644 ruff.toml

diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..97491fb
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,101 @@
+[build-system]
+requires = ["poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"
+
+[tool.poetry]
+name = "hawc_hal"
+version = "1.1"
+authors = ["Giacomo Vianello <giacomov@stanford.edu>"]
+maintainers = [
+  "Xiaojie Wang <xwang32@mtu.edu>",
+  "Ramiro Torres-Escobedo <rtorrese@proton.me>",
+]
+description = "HAWC Accelerated Likelihood; Read and handle HAWC data"
+license = "BSD-3-Clause"
+packages = [{ include = "hawc_hal" }, { include = "scripts" }]
+
+[tool.poetry.urls]
+homepage = "https://threeml.readthedocs.io/en/stable/index.html"
+repository = "https://github.com/threeML/hawc_hal"
+documentation = "https://threeml.readthedocs.io/en/stable/notebooks/hal_example.html"
+"Bug Tracker" = "https://github.com/threeML/hawc_hal/issues"
+
+[tool.poetry.dependencies]
+python = "^3.8"
+numpy = ">=1.14"
+healpy = "*"
+threeml = "*"
+astromodels = "*"
+pandas = "*"
+pytest = "*"
+six = "*"
+astropy = "*"
+scipy = "*"
+matplotlib = "*"
+numba = "*"
+reproject = "*"
+tqdm = "*"
+uproot = "*"
+awkward = "*"
+mplhep = "*"
+hist = "*"
+ruff = "*"
+pyright = "*"
+
+
+[tool.poetry.scripts]
+hdf5tofits = "scripts.hal_hdf5_to_fits:main"
+halfitpointsrc = "scripts.hal_fit_point_source:main"
+
+[tool.ruff]
+include = [
+  "pyproject.toml",
+  "hawc_hal/**/*.py",
+  "scripts/**/*.py",
+  "notebooks/**/*.ipynb",
+]
+exclude = [
+  "codecov.yml",
+  "data/*",
+  "ci/*",
+  "notebooks",
+  "tests",
+  ".eggs",
+  ".git",
+  ".git-rewrite",
+  ".ipynb_checkpoints",
+  ".mypy_cache",
+  ".nox",
+  ".pants.d",
+  ".pyenv",
+  ".pytest_cache",
+  ".pytype",
+  ".ruff_cache",
+  ".venv",
+  ".vscode",
+  "__pypackages__",
+  "_build",
+  "buck-out",
+  "build",
+  "dist",
+  "node_modules",
+  "site-packages",
+  "venv",
+]
+line-length = 88
+indent-width = 4
+docstring-code-format = false
+
+
+[tool.ruff.format]
+# Like Black, use double quotes for strings.
+quote-style = "double"
+
+# Like Black, indent with spaces, rather than tabs.
+indent-style = "space"
+
+# Like Black, respect magic trailing commas.
+skip-magic-trailing-comma = false
+
+# Like Black, automatically detect the appropriate line ending.
+line-ending = "auto"
diff --git a/ruff.toml b/ruff.toml
new file mode 100644
index 0000000..eb00e5a
--- /dev/null
+++ b/ruff.toml
@@ -0,0 +1,15 @@
+line-length = 90
+target-version = "py310"
+[lint]
+select = ["F"]
+ignore = ["E501"]
+
+# [lint.isort]
+# line-after-imports = 2
+
+[format]
+quote-style = "double"
+docstring-code-format = true
+docstring-code-line-length = 60
+skip-magic-trailing-comma = true
+# quote-style = "single"

From efa331c225c18098d0706adef61a99a1340c4cb1 Mon Sep 17 00:00:00 2001
From: Ramiro Torres-Escobedo <rtorrese@proton.me>
Date: Mon, 20 Jan 2025 13:49:27 +0800
Subject: [PATCH 2/7] getting changes from test branch

---
 hawc_hal/HAL.py | 324 +++++++++++++++++++++++++-----------------------
 pyproject.toml  |   5 +-
 2 files changed, 175 insertions(+), 154 deletions(-)

diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py
index d83fb15..d12861e 100644
--- a/hawc_hal/HAL.py
+++ b/hawc_hal/HAL.py
@@ -4,6 +4,7 @@
 import contextlib
 import copy
 from builtins import range, str
+from typing import TypeAlias as T
 from typing import Union
 
 import astromodels
@@ -14,6 +15,8 @@
 from astromodels import Parameter
 from astropy.convolution import Gaussian2DKernel
 from astropy.convolution import convolve_fft as convolve
+from matplotlib.figure import Figure
+from numpy.typing import NDArray
 from past.utils import old_div
 from scipy.stats import poisson
 from threeML.io.logging import setup_logger
@@ -22,12 +25,17 @@
 from threeML.utils.statistics.gammaln import logfactorial
 from tqdm.auto import tqdm
 
-from hawc_hal.convolved_source import (ConvolvedExtendedSource2D,
-                                       ConvolvedExtendedSource3D,
-                                       ConvolvedPointSource,
-                                       ConvolvedSourcesContainer)
-from hawc_hal.healpix_handling import (FlatSkyToHealpixTransform,
-                                       SparseHealpix, get_gnomonic_projection)
+from hawc_hal.convolved_source import (
+    ConvolvedExtendedSource2D,
+    ConvolvedExtendedSource3D,
+    ConvolvedPointSource,
+    ConvolvedSourcesContainer,
+)
+from hawc_hal.healpix_handling import (
+    FlatSkyToHealpixTransform,
+    SparseHealpix,
+    get_gnomonic_projection,
+)
 from hawc_hal.log_likelihood import log_likelihood
 from hawc_hal.maptree import map_tree_factory
 from hawc_hal.maptree.data_analysis_bin import DataAnalysisBin
@@ -38,6 +46,7 @@
 
 log = setup_logger(__name__)
 log.propagate = False
+ndarray: T = NDArray[np.float64]
 
 
 class HAL(PluginPrototype):
@@ -124,7 +133,7 @@ def __init__(
         self._all_planes = list(self._maptree.analysis_bins_labels)
 
         # The active planes list always contains the list of *indexes* of the active planes
-        self._active_planes = None
+        self._active_planes: list[str] | None = None
 
         # Set up the transformations from the flat-sky projection to Healpix, as well as the list of active pixels
         # (one for each energy/nHit bin). We make a separate transformation because different energy bins might have
@@ -269,9 +278,9 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non
 
         # Check for legal input
         if bin_id_min is not None:
-            assert (
-                bin_id_max is not None
-            ), "If you provide a minimum bin, you also need to provide a maximum bin."
+            assert bin_id_max is not None, (
+                "If you provide a minimum bin, you also need to provide a maximum bin."
+            )
 
             # Make sure they are integers
             bin_id_min = int(bin_id_min)
@@ -281,14 +290,16 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non
             for this_bin in range(bin_id_min, bin_id_max + 1):
                 this_bin = str(this_bin)
                 if this_bin not in self._all_planes:
-                    raise ValueError(f"Bin {this_bin} is not contained in this maptree.")
+                    raise ValueError(
+                        f"Bin {this_bin} is not contained in this maptree."
+                    )
 
                 self._active_planes.append(this_bin)
 
         else:
-            assert (
-                bin_id_max is None
-            ), "If you provie a maximum bin, you also need to provide a minimum bin."
+            assert bin_id_max is None, (
+                "If you provie a maximum bin, you also need to provide a minimum bin."
+            )
 
             assert bin_list is not None
 
@@ -297,7 +308,9 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non
             for this_bin in bin_list:
                 # if not this_bin in self._all_planes:
                 if this_bin not in self._all_planes:
-                    raise ValueError(f"Bin {this_bin} is not contained in this maptree.")
+                    raise ValueError(
+                        f"Bin {this_bin} is not contained in this maptree."
+                    )
 
                 self._active_planes.append(this_bin)
 
@@ -385,27 +398,26 @@ def set_model(self, likelihood_model_instance):
 
                 self._convolved_ext_sources.append(this_convolved_ext_source)
 
-    def get_excess_background(self, ra: float, dec: float, radius: float):
-        """Calculates excess (data-bkg), background, and model counts at
-        different radial distances from origin of radial profile.
-
-
-        Parameters
-        ----------
-        ra : float
-            RA of origin of radial profile
-        dec : float
-           Dec of origin of radial profile
-        radius : float
-           distance from origin of radial profile
-
-        Returns
-        -------
-           returns a tuple of numpy arrays with info of areas (steradian) and
-           signal excess, background, and model in units of counts to be used
-           in the get_radial_profile method.
+    def _get_excess_background(
+        self, ra: float, dec: float, radius: float
+    ) -> tuple[ndarray, ...]:
+        """Compute the signal excess, background, and model counts at different radial
+        distances from the origin of the radial profile.
+
+        :param ra: right ascension of the origin of the radial profile
+        :param dec: declination of the origin of the radial profile
+        :param radius: distance from the origin of the radial profile
+        :raises ValueError: if no active planes have been set
+        :return: a tuple of numpy arrays with info of areas (steradian) and
+        signal excess, background, and model in units of counts to be used in
+        the `get_radial_profile` method
         """
 
+        if self._active_planes is None:
+            raise ValueError(
+                "No active planes have been set. Please use set_active_measurements()"
+            )
+
         radius_radians = np.deg2rad(radius)
 
         total_counts = np.zeros(len(self._active_planes), dtype=float)
@@ -415,47 +427,46 @@ def get_excess_background(self, ra: float, dec: float, radius: float):
         signal = np.zeros_like(total_counts)
         area = np.zeros_like(total_counts)
 
-        n_point_sources = self._likelihood_model.get_number_of_point_sources()
-        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
+        n_point_sources = self._likelihood_model.get_number_of_point_sources()  # type: ignore
+        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()  # type: ignore
 
         longitude = ra_to_longitude(ra)
         latitude = dec
         center = hp.ang2vec(longitude, latitude, lonlat=True)
 
-        for i, energy_id in enumerate(self._active_planes):
-            data_analysis_bin = self._maptree[energy_id]
-            this_nside = data_analysis_bin.observation_map.nside
+        # NOTE: the nside is the same for all bins, so it is okay to read it from the first bin
+        this_nside = self._maptree[self._active_planes[0]].nside
+        radial_bin_pixels = hp.query_disc(
+            this_nside, center, radius_radians, inclusive=False
+        )
 
-            radial_bin_pixels = hp.query_disc(
-                this_nside, center, radius_radians, inclusive=False
-            )
+        # select the pixels that are only within the radial bin
+        pixels_within_rad_bin = np.isin(
+            self._active_pixels[self._active_planes[0]], radial_bin_pixels
+        )
 
-            # calculate the areas per bin by the product
-            # of pixel area by the number of pixels at each radial bin
-            area[i] = hp.nside2pixarea(this_nside) * radial_bin_pixels.shape[0]
+        # calculate the areas per bin by the product
+        # of pixel area by the number of pixels at each radial bin
+        this_area = hp.nside2pixarea(this_nside) * radial_bin_pixels.shape[0]
+        area = np.full(this_area, len(self._active_planes))
 
-            # NOTE: select active pixels according to each radial bin
-            bin_active_pixel_indexes = np.intersect1d(
-                self._active_pixels[energy_id], radial_bin_pixels, return_indices=True
-            )[1]
+        for i, energy_id in enumerate(self._active_planes):
+            data_analysis_bin: DataAnalysisBin = self._maptree[energy_id]
 
             # obtain the excess, background, and expected excess at
             # each radial bin
-            data = data_analysis_bin.observation_map.as_partial()
-            bkg = data_analysis_bin.background_map.as_partial()
-            mdl = self._get_model_map(
+            data: ndarray = data_analysis_bin.observation_map.as_partial()
+            bkg: ndarray = data_analysis_bin.background_map.as_partial()
+            mdl: ndarray = self._get_model_map(
                 energy_id, n_point_sources, n_ext_sources
             ).as_partial()
 
-            # select counts only from the pixels within specifid distance from
-            # origin of radial profile
-            bin_data = np.array([data[i] for i in bin_active_pixel_indexes])
-            bin_bkg = np.array([bkg[i] for i in bin_active_pixel_indexes])
-            bin_model = np.array([mdl[i] for i in bin_active_pixel_indexes])
+            # select the information only from the pixels that are within the radial
+            # bin from origin of radial profile
 
-            this_data_tot = np.sum(bin_data)
-            this_bkg_tot = np.sum(bin_bkg)
-            this_model_tot = np.sum(bin_model)
+            this_data_tot = data[pixels_within_rad_bin].sum()
+            this_bkg_tot = bkg[pixels_within_rad_bin].sum()
+            this_model_tot = mdl[pixels_within_rad_bin].sum()
 
             background[i] = this_bkg_tot
             observation[i] = this_data_tot
@@ -464,37 +475,28 @@ def get_excess_background(self, ra: float, dec: float, radius: float):
 
         return area, signal, background, model
 
-    def get_radial_profile(
+    def _get_radial_profile(
         self,
         ra: float,
         dec: float,
-        active_planes: list = None,
+        active_planes: list[str] | None = None,
         max_radius: float = 3.0,
         n_radial_bins: int = 30,
-        model_to_subtract: astromodels.Model = None,
+        model_to_subtract: astromodels.Model | None = None,
         subtract_model_from_model: bool = False,
-    ):
-        """Calculates radial profiles for a source in units of excess counts
-           per steradian
-
-        Args:
-            ra (float): RA of origin of radial profile
-            dec (float): Declincation of origin of radial profile
-            active_planes (np.ndarray, optional): List of active planes over
-            which to average. Defaults to None.
-            max_radius (float, optional): Radius up to which evaluate the
-            radial profile. Defaults to 3.0.
-            n_radial_bins (int, optional): Number of radial bins to use for
-            the profile. Defaults to 30.
-            model_to_subtract (astromodels.model, optional): Another model to
-            subtract from the data excess. Defaults to None.
-            subtract_model_from_model (bool, optional): If True, and
-            model_to_subtract is not None,
-            subtract model from model too. Defaults to False.
-
-        Returns:
-            tuple(np.ndarray): returns list of radial distances, excess expected
-            counts, excess counts, counts uncertainty, and list of sorted active_planes
+    ) -> tuple[ndarray, ...]:
+        """Cacluate radial prfiles for a source in units of excess counts per steradian
+
+        :param ra: RA of origin of radial profile
+        :param dec: declination of origin of radial profile
+        :param active_planes: list of active planes over which to average, defaults to None
+        :param max_radius: radius up to which evaluate the radial profile, defaults to 3.0
+        :param n_radial_bins: number of radial bins to use for the profile, defaults to 30
+        :param model_to_subtract: another model to subtract from the data excess, defaults to None
+        :param subtract_model_from_model: if True, and model_to_subtract is not None,
+        subtract model from model too, defaults to False
+        :return: returns list of radial distances, excess expected counts, excess counts,
+        counts uncertainty, and list of sorted active_planes
         """
         # default is to use all active bins
         if active_planes is None:
@@ -512,7 +514,10 @@ def get_radial_profile(
         # The area of each ring is then given by the difference between two
         # subsequent circe areas.
         area = np.array(
-            [self.get_excess_background(ra, dec, r + offset * delta_r)[0] for r in radii]
+            [
+                self._get_excess_background(ra, dec, r + offset * delta_r)[0]
+                for r in radii
+            ]
         )
 
         temp = area[1:] - area[:-1]
@@ -520,7 +525,10 @@ def get_radial_profile(
 
         # signals
         signal = np.array(
-            [self.get_excess_background(ra, dec, r + offset * delta_r)[1] for r in radii]
+            [
+                self._get_excess_background(ra, dec, r + offset * delta_r)[1]
+                for r in radii
+            ]
         )
 
         temp = signal[1:] - signal[:-1]
@@ -528,7 +536,10 @@ def get_radial_profile(
 
         # backgrounds
         bkg = np.array(
-            [self.get_excess_background(ra, dec, r + offset * delta_r)[2] for r in radii]
+            [
+                self._get_excess_background(ra, dec, r + offset * delta_r)[2]
+                for r in radii
+            ]
         )
 
         temp = bkg[1:] - bkg[:-1]
@@ -539,7 +550,10 @@ def get_radial_profile(
         # model
         # convert 'top hat' excess into 'ring' excesses.
         model = np.array(
-            [self.get_excess_background(ra, dec, r + offset * delta_r)[3] for r in radii]
+            [
+                self._get_excess_background(ra, dec, r + offset * delta_r)[3]
+                for r in radii
+            ]
         )
 
         temp = model[1:] - model[:-1]
@@ -551,7 +565,7 @@ def get_radial_profile(
 
             model_subtract = np.array(
                 [
-                    self.get_excess_background(ra, dec, r + offset * delta_r)[3]
+                    self._get_excess_background(ra, dec, r + offset * delta_r)[3]
                     for r in radii
                 ]
             )
@@ -573,17 +587,14 @@ def get_radial_profile(
         # them to the data later. Weight is normalized (sum of weights over
         # the bins = 1).
 
-        np.array(self.get_excess_background(ra, dec, max_radius)[1])[good_planes]
+        # TODO: check if this is the correct way to calculate the weights
+        # np.array(self._get_excess_background(ra, dec, max_radius)[1])[good_planes]
 
-        total_bkg = np.array(self.get_excess_background(ra, dec, max_radius)[2])[
-            good_planes
-        ]
-
-        total_model = np.array(self.get_excess_background(ra, dec, max_radius)[3])[
-            good_planes
-        ]
+        total_bkg = self._get_excess_background(ra, dec, max_radius)[2][good_planes]
+        total_model = self._get_excess_background(ra, dec, max_radius)[3][good_planes]
 
-        w = np.divide(total_model, total_bkg)
+        # w = np.divide(total_model, total_bkg)
+        w = total_model / total_bkg
         weight = np.array([w / np.sum(w) for _ in radii])
 
         # restrict profiles to the user-specified analysis bins
@@ -598,65 +609,56 @@ def get_radial_profile(
         excess_error = np.sqrt(np.sum(counts * weight * weight / (area * area), axis=1))
         excess_model = np.average(model / area, weights=weight, axis=1)
 
-        return radii, excess_model, excess_data, excess_error, sorted(plane_ids)
+        return (
+            radii,
+            excess_model,
+            excess_data,
+            excess_error,
+            sorted(plane_ids),  # pyright:ignore
+            w[0, :],
+        )
 
     def plot_radial_profile(
         self,
         ra: float,
         dec: float,
-        active_planes: list = None,
+        active_planes: list[str] | None = None,
         max_radius: float = 3.0,
         n_radial_bins: int = 30,
-        model_to_subtract: astromodels.Model = None,
+        model_to_subtract: astromodels.Model | None = None,
         subtract_model_from_model: bool = False,
-    ):
-        """Plots radial profiles of data-background & model
-
-        Args:
-            ra (float): RA of origin of radial profile
-            dec (float): Declination of origin of radial profile.
-            active_planes (np.ndarray, optional): List of analysis bins over
-            which to average.
-            Defaults to None.
-            max_radius (float, optional): Radius up to which the radial profile
-            is evaluate; also used as the radius for the disk to calculate the
-            gamma/hadron weights. Defaults to 3.0.
-            n_radial_bins (int, optional): number of radial bins used for ring
-            calculation. Defaults to 30.
-            model_to_subtract (astromodels.model, optional): Another model that
-            is to be subtracted from the data excess. Defaults to None.
-            subtract_model_from_model (bool, optional): If True and
-            model_to_subtract is not None, subtract from model too.
-            Defaults to False.
-
-        Returns:
-            tuple(matplotlib.pyplot.Figure, pd.DataFrame): plot of data-background
-            & model radial profile for source and a dataframe with all
-            values for easy retrieval
+    ) -> tuple[Figure, pd.DataFrame]:
+        """Plot radial prfiles for a source in units of excess counts per steradian
+
+        :param ra: RA of origin of radial profile (J2000)
+        :param dec: declination of origin of radial profile (J2000)
+        :param active_planes: list of active planes over which to average, defaults to None
+        :param max_radius: radius up to which evaluate the radial profile, defaults to 3.0
+        :param n_radial_bins: number of radial bins to use for the profile, defaults to 30
+        :param model_to_subtract: another model to subtract from the data excess, defaults to None
+        :param subtract_model_from_model: if True, and model_to_subtract is not None,
+        subtract model from model too, defaults to False
+        :return: radial profile figure and a dataframe with all values for easy retrieval
         """
 
-        (
-            radii,
-            excess_model,
-            excess_data,
-            excess_error,
-            plane_ids,
-        ) = self.get_radial_profile(
-            ra,
-            dec,
-            active_planes,
-            max_radius,
-            n_radial_bins,
-            model_to_subtract,
-            subtract_model_from_model,
+        (radii, excess_model, excess_data, excess_error, plane_ids, weights) = (
+            self._get_radial_profile(
+                ra,
+                dec,
+                active_planes,
+                max_radius,
+                n_radial_bins,
+                model_to_subtract,
+                subtract_model_from_model,
+            )
         )
 
         # add a dataframe for easy retrieval for calculations of surface
         # brighntess, if necessary.
-        df = pd.DataFrame(columns=["Excess", "Bkg", "Model"], index=radii)
+        df = pd.DataFrame(columns=["Excess", "Error", "Model"], index=radii)
         df.index.name = "Radii"
         df["Excess"] = excess_data
-        df["Bkg"] = excess_error
+        df["Error"] = excess_error
         df["Model"] = excess_model
 
         fig, ax = plt.subplots(figsize=(10, 8))
@@ -673,7 +675,9 @@ def plot_radial_profile(
 
         plt.plot(radii, excess_model, color="red", label="Model")
 
-        plt.legend(bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16)
+        plt.legend(
+            bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16
+        )
         plt.axhline(0, color="deepskyblue", linestyle="--")
 
         x_limits = [0, max_radius]
@@ -780,7 +784,9 @@ def display_spectrum(self):
 
         yerr = [yerr_high, yerr_low]
 
-        return self._plot_spectrum(net_counts, yerr, model_only, residuals, residuals_err)
+        return self._plot_spectrum(
+            net_counts, yerr, model_only, residuals, residuals_err
+        )
 
     def _plot_spectrum(self, net_counts, yerr, model_only, residuals, residuals_err):
         fig, subs = plt.subplots(
@@ -848,7 +854,9 @@ def get_log_like(self, individual_bins=False, return_null=False):
             assert (
                 n_point_sources == self._convolved_point_sources.n_sources_in_cache
                 and n_ext_sources == self._convolved_ext_sources.n_sources_in_cache
-            ), "The number of sources has changed. Please re-assign the model to the plugin."
+            ), (
+                "The number of sources has changed. Please re-assign the model to the plugin."
+            )
 
         # This will hold the total log-likelihood
 
@@ -1053,7 +1061,9 @@ def _get_expectation(
             )
 
             # Now multiply by the pixel area of the new map to go back to flux
-            this_model_map_hpx *= hp.nside2pixarea(data_analysis_bin.nside, degrees=True)
+            this_model_map_hpx *= hp.nside2pixarea(
+                data_analysis_bin.nside, degrees=True
+            )
 
         else:
             # No sources
@@ -1178,7 +1188,9 @@ def display_fit(self, smoothing_kernel_sigma=0.1, display_colorbar=False):
             subs[i][0].set_title("model, bin {}".format(data_analysis_bin.name))
 
             # Plot data map
-            images[1] = subs[i][1].imshow(proj_data, origin="lower", vmin=vmin, vmax=vmax)
+            images[1] = subs[i][1].imshow(
+                proj_data, origin="lower", vmin=vmin, vmax=vmax
+            )
             subs[i][1].set_title("excess, bin {}".format(data_analysis_bin.name))
 
             # Plot background map.
@@ -1289,7 +1301,7 @@ def get_number_of_data_points(self):
 
         return n_points
 
-    def _get_model_map(self, plane_id, n_pt_src, n_ext_src):
+    def _get_model_map(self, plane_id, n_pt_src, n_ext_src) -> SparseHealpix:
         """
         This function returns a model map for a particular bin
         """
@@ -1298,7 +1310,9 @@ def _get_model_map(self, plane_id, n_pt_src, n_ext_src):
             raise ValueError(f"{plane_id} not a plane in the current model")
 
         model_map = SparseHealpix(
-            self._get_expectation(self._maptree[plane_id], plane_id, n_pt_src, n_ext_src),
+            self._get_expectation(
+                self._maptree[plane_id], plane_id, n_pt_src, n_ext_src
+            ),
             self._active_pixels[plane_id],
             self._maptree[plane_id].observation_map.nside,
         )
@@ -1371,13 +1385,17 @@ def _write_a_map(self, file_name, which, fluctuate=False, return_map=False):
         if return_map:
             return new_map_tree
 
-    def write_model_map(self, file_name, poisson_fluctuate=False, test_return_map=False):
+    def write_model_map(
+        self, file_name, poisson_fluctuate=False, test_return_map=False
+    ):
         """
         This function writes the model map to a file.
         The interface is based off of HAWCLike for consistency
         """
         if test_return_map:
-            log.warning("test_return_map=True should only be used for testing purposes!")
+            log.warning(
+                "test_return_map=True should only be used for testing purposes!"
+            )
         return self._write_a_map(file_name, "model", poisson_fluctuate, test_return_map)
 
     def write_residual_map(self, file_name, test_return_map=False):
@@ -1386,5 +1404,7 @@ def write_residual_map(self, file_name, test_return_map=False):
         The interface is based off of HAWCLike for consistency
         """
         if test_return_map:
-            log.warning("test_return_map=True should only be used for testing purposes!")
+            log.warning(
+                "test_return_map=True should only be used for testing purposes!"
+            )
         return self._write_a_map(file_name, "residual", False, test_return_map)
diff --git a/pyproject.toml b/pyproject.toml
index 97491fb..c55f705 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -14,7 +14,7 @@ description = "HAWC Accelerated Likelihood; Read and handle HAWC data"
 license = "BSD-3-Clause"
 packages = [{ include = "hawc_hal" }, { include = "scripts" }]
 
-[tool.poetry.urls]
+[tool.project.urls]
 homepage = "https://threeml.readthedocs.io/en/stable/index.html"
 repository = "https://github.com/threeML/hawc_hal"
 documentation = "https://threeml.readthedocs.io/en/stable/notebooks/hal_example.html"
@@ -84,10 +84,11 @@ exclude = [
 ]
 line-length = 88
 indent-width = 4
-docstring-code-format = false
 
 
 [tool.ruff.format]
+docstring-code-format = false
+
 # Like Black, use double quotes for strings.
 quote-style = "double"
 

From 9ee22c8ba30508ab79a326d36573b17f26b59c81 Mon Sep 17 00:00:00 2001
From: Ramiro Torres-Escobedo <rtorrese@proton.me>
Date: Tue, 21 Jan 2025 22:47:34 +0800
Subject: [PATCH 3/7] acquiring changes from alternative radial profile branch

---
 hawc_hal/HAL.py                        | 81 +++++++++++---------------
 hawc_hal/tests/test_radial_profiles.py | 11 ++--
 2 files changed, 41 insertions(+), 51 deletions(-)

diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py
index d12861e..1063f50 100644
--- a/hawc_hal/HAL.py
+++ b/hawc_hal/HAL.py
@@ -413,12 +413,15 @@ def _get_excess_background(
         the `get_radial_profile` method
         """
 
+        from math import radians
+
         if self._active_planes is None:
             raise ValueError(
                 "No active planes have been set. Please use set_active_measurements()"
             )
 
-        radius_radians = np.deg2rad(radius)
+        # healpy likes to work with radians
+        radius_radians: float = radians(radius)
 
         total_counts = np.zeros(len(self._active_planes), dtype=float)
         background = np.zeros_like(total_counts)
@@ -441,14 +444,14 @@ def _get_excess_background(
         )
 
         # select the pixels that are only within the radial bin
-        pixels_within_rad_bin = np.isin(
+        pixels_within_rad_bin = np.in1d(
             self._active_pixels[self._active_planes[0]], radial_bin_pixels
         )
 
         # calculate the areas per bin by the product
         # of pixel area by the number of pixels at each radial bin
         this_area = hp.nside2pixarea(this_nside) * radial_bin_pixels.shape[0]
-        area = np.full(this_area, len(self._active_planes))
+        area = np.full(len(self._active_planes), this_area)
 
         for i, energy_id in enumerate(self._active_planes):
             data_analysis_bin: DataAnalysisBin = self._maptree[energy_id]
@@ -506,8 +509,8 @@ def _get_radial_profile(
         good_planes = [plane_id in active_planes for plane_id in self._active_planes]
         plane_ids = set(active_planes) & set(self._active_planes)
 
-        offset = 0.5
-        delta_r = 1.0 * max_radius / n_radial_bins
+        offset = 0.50
+        delta_r = (1.0 * max_radius) / n_radial_bins
         radii = np.array([delta_r * (r + offset) for r in range(n_radial_bins)])
 
         # Get area of all pixels in a given circle
@@ -593,9 +596,8 @@ def _get_radial_profile(
         total_bkg = self._get_excess_background(ra, dec, max_radius)[2][good_planes]
         total_model = self._get_excess_background(ra, dec, max_radius)[3][good_planes]
 
-        # w = np.divide(total_model, total_bkg)
-        w = total_model / total_bkg
-        weight = np.array([w / np.sum(w) for _ in radii])
+        w = np.divide(total_model, total_bkg)
+        weight = np.array([w / w.sum() for _ in radii])
 
         # restrict profiles to the user-specified analysis bins
         area = area[:, good_planes]
@@ -615,7 +617,7 @@ def _get_radial_profile(
             excess_data,
             excess_error,
             sorted(plane_ids),  # pyright:ignore
-            w[0, :],
+            weight[0, :],
         )
 
     def plot_radial_profile(
@@ -627,7 +629,7 @@ def plot_radial_profile(
         n_radial_bins: int = 30,
         model_to_subtract: astromodels.Model | None = None,
         subtract_model_from_model: bool = False,
-    ) -> tuple[Figure, pd.DataFrame]:
+    ) -> tuple[Figure, pd.DataFrame, pd.DataFrame]:
         """Plot radial prfiles for a source in units of excess counts per steradian
 
         :param ra: RA of origin of radial profile (J2000)
@@ -655,15 +657,19 @@ def plot_radial_profile(
 
         # add a dataframe for easy retrieval for calculations of surface
         # brighntess, if necessary.
-        df = pd.DataFrame(columns=["Excess", "Error", "Model"], index=radii)
+        df = pd.DataFrame(
+            {"Excess": excess_data, "Error": excess_error, "Model": excess_model},
+            index=radii,
+            dtype=float,
+        )
         df.index.name = "Radii"
-        df["Excess"] = excess_data
-        df["Error"] = excess_error
-        df["Model"] = excess_model
 
-        fig, ax = plt.subplots(figsize=(10, 8))
+        dfw = pd.DataFrame({"Bins": plane_ids, "Weights": weights})
 
-        plt.errorbar(
+        fig = plt.figure(figsize=(10, 8))
+        ax = fig.add_subplot()
+
+        ax.errorbar(
             radii,
             excess_data,
             yerr=excess_error,
@@ -673,51 +679,34 @@ def plot_radial_profile(
             fmt=".",
         )
 
-        plt.plot(radii, excess_model, color="red", label="Model")
+        ax.plot(radii, excess_model, color="red", label="Model")
 
-        plt.legend(
+        ax.legend(
             bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16
         )
-        plt.axhline(0, color="deepskyblue", linestyle="--")
-
-        x_limits = [0, max_radius]
-        plt.xlim(x_limits)
-        plt.xticks(fontsize=18)
-        plt.yticks(fontsize=18)
+        ax.axhline(0, color="deepskyblue", linestyle="--")
 
-        plt.ylabel(r"Apparent Radial Excess [sr$^{-1}$]", fontsize=18)
-        plt.xlabel(
-            f"Distance from source at ({ra:0.2f} $^{{\circ}}$, {dec:0.2f} $^{{\circ}}$)",
-            fontsize=18,
-        )
+        ax.set_xlim(left=0, right=max_radius)
+        ax.tick_params(axis="both", labelsize=18)
 
         if len(plane_ids) == 1:
             title = f"Radial Profile, bin {plane_ids[0]}"
 
         else:
             title = "Radial Profile"
-            # tmptitle = f"Radial Profile, bins \n{plane_ids}"
-            # width = 80
-            # title = "\n".join(
-            # tmptitle[i : i + width] for i in range(0, len(tmptitle), width)
-            # )
-            # title = tmptitle
-
-        plt.title(title)
 
+        plt.ylabel(r"Apparent Radial Excess [sr$^{-1}$]", fontsize=18)
+        plt.xlabel(
+            f"Distance from source at ({ra:0.2f} $^{{\circ}}$, {dec:0.2f} $^{{\circ}}$)",
+            fontsize=18,
+        )
+        ax.set_title(title)
         ax.grid(True)
 
         with contextlib.suppress(Exception):
             plt.tight_layout()
-        # try:
-        #
-        # plt.tight_layout()
-        #
-        # except Exception:
-        #
-        # pass
-
-        return fig, df
+
+        return fig, df, dfw
 
     def display_spectrum(self):
         """
diff --git a/hawc_hal/tests/test_radial_profiles.py b/hawc_hal/tests/test_radial_profiles.py
index 73518da..82f0299 100644
--- a/hawc_hal/tests/test_radial_profiles.py
+++ b/hawc_hal/tests/test_radial_profiles.py
@@ -42,28 +42,29 @@ def test_simulation(test_fit):
 
 def test_plots(test_fit):
 
-    jl, hawc, pts_model, param_df, like_df, data = test_fit
+    _, hawc, pts_model, *_ = test_fit
 
     ra = pts_model.pts.position.ra.value
     dec = pts_model.pts.position.dec.value
-    radius = 1.5
+    radius = 3.0
 
     bins = hawc._active_planes
 
-    fig, table = hawc.plot_radial_profile(
+    fig, table, _ = hawc.plot_radial_profile(
         ra,
         dec,
         bins,
         radius,
-        n_radial_bins=25,
+        n_radial_bins=30,
     )
+
     fig.savefig("hal_src_radial_profile.png")
     table.to_hdf("hal_src_radial_table.hd5", key="radial")
 
     prog_bar = tqdm(total=len(hawc._active_planes), desc="Smoothing planes")
     for bin in hawc._active_planes:
 
-        fig, table = hawc.plot_radial_profile(ra, dec, f"{bin}", radius)
+        fig, table, _ = hawc.plot_radial_profile(ra, dec, f"{bin}", radius)
         fig.savefig(f"hal_src_radial_profile_bin{bin}.png")
 
         table.to_hdf(f"hal_src_radial_table_{bin}.hd5", key="radial")

From 69d2186b9f07c2e990c609f5cf0b0253029ec4ed Mon Sep 17 00:00:00 2001
From: Ramiro Torres-Escobedo <rtorrese@proton.me>
Date: Wed, 5 Feb 2025 21:33:22 +0800
Subject: [PATCH 4/7] improving the performance of computing the radial
 profiles even more

---
 hawc_hal/HAL.py | 132 +++++++++++++++++++++++++-----------------------
 1 file changed, 69 insertions(+), 63 deletions(-)

diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py
index 1063f50..6a6b404 100644
--- a/hawc_hal/HAL.py
+++ b/hawc_hal/HAL.py
@@ -12,7 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-from astromodels import Parameter
+from astromodels import Model, Parameter
 from astropy.convolution import Gaussian2DKernel
 from astropy.convolution import convolve_fft as convolve
 from matplotlib.figure import Figure
@@ -122,7 +122,7 @@ def __init__(
         # python3 new way of doing things
         super().__init__(name, self._nuisance_parameters)
 
-        self._likelihood_model = None
+        self._likelihood_model: Model | None = None
 
         # These lists will contain the maps for the point sources
         self._convolved_point_sources = ConvolvedSourcesContainer()
@@ -209,10 +209,9 @@ def psf_integration_method(self):
 
     @psf_integration_method.setter
     def psf_integration_method(self, mode):
-        assert mode.lower() in [
-            "exact",
-            "fast",
-        ], "PSF integration method must be either 'exact' or 'fast'"
+        assert mode.lower() in ["exact", "fast"], (
+            "PSF integration method must be either 'exact' or 'fast'"
+        )
 
         self._psf_integration_method = mode.lower()
 
@@ -290,9 +289,7 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non
             for this_bin in range(bin_id_min, bin_id_max + 1):
                 this_bin = str(this_bin)
                 if this_bin not in self._all_planes:
-                    raise ValueError(
-                        f"Bin {this_bin} is not contained in this maptree."
-                    )
+                    raise ValueError(f"Bin {this_bin} is not contained in this maptree.")
 
                 self._active_planes.append(this_bin)
 
@@ -308,9 +305,7 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non
             for this_bin in bin_list:
                 # if not this_bin in self._all_planes:
                 if this_bin not in self._all_planes:
-                    raise ValueError(
-                        f"Bin {this_bin} is not contained in this maptree."
-                    )
+                    raise ValueError(f"Bin {this_bin} is not contained in this maptree.")
 
                 self._active_planes.append(this_bin)
 
@@ -398,6 +393,40 @@ def set_model(self, likelihood_model_instance):
 
                 self._convolved_ext_sources.append(this_convolved_ext_source)
 
+    def _prepare_data_for_radial_profile(self) -> None:
+        """Organize the data for radial profiles. This requires getting the observation,
+        background, and model maps for the the current region of interest. We use this in the
+        `get_excess_background` method to retrieve the values for each map that coincide
+        with every radial bin.
+
+        :raises ValueError: Check that the likelihood model and active planes have been set.
+        """
+        self._data_for_radial_profile: collections.OrderedDict[
+            str, collections.OrderedDict[str, ndarray]
+        ] = collections.OrderedDict()
+
+        if self._likelihood_model is None:
+            raise ValueError("No likelihood model has been set.")
+
+        if self._active_planes is None:
+            raise ValueError("Active planes have not been set")
+
+        n_point_sources = self._likelihood_model.get_number_of_point_sources()
+        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
+        for i, energy_id in enumerate(self._active_planes):
+            data_analysis_bin: DataAnalysisBin = self._maptree[energy_id]
+
+            # obtain the excess, background, and expected excess at
+            # each radial bin
+            data: ndarray = data_analysis_bin.observation_map.as_partial()
+            bkg: ndarray = data_analysis_bin.background_map.as_partial()
+            mdl: ndarray = self._get_model_map(
+                energy_id, n_point_sources, n_ext_sources
+            ).as_partial()
+            self._data_for_radial_profile[energy_id] = collections.OrderedDict(
+                {"data": data, "bkg": bkg, "mdl": mdl}
+            )
+
     def _get_excess_background(
         self, ra: float, dec: float, radius: float
     ) -> tuple[ndarray, ...]:
@@ -454,22 +483,24 @@ def _get_excess_background(
         area = np.full(len(self._active_planes), this_area)
 
         for i, energy_id in enumerate(self._active_planes):
-            data_analysis_bin: DataAnalysisBin = self._maptree[energy_id]
-
+            # data_analysis_bin: DataAnalysisBin = self._maptree[energy_id]
             # obtain the excess, background, and expected excess at
             # each radial bin
-            data: ndarray = data_analysis_bin.observation_map.as_partial()
-            bkg: ndarray = data_analysis_bin.background_map.as_partial()
-            mdl: ndarray = self._get_model_map(
-                energy_id, n_point_sources, n_ext_sources
-            ).as_partial()
+            # data: ndarray = data_analysis_bin.observation_map.as_partial()
+            # bkg: ndarray = data_analysis_bin.background_map.as_partial()
+            # mdl: ndarray = self._get_model_map(
+            #     energy_id, n_point_sources, n_ext_sources
+            # ).as_partial()
+            data: ndarray = self._data_for_radial_profile[energy_id]["data"]
+            bkg: ndarray = self._data_for_radial_profile[energy_id]["bkg"]
+            mdl: ndarray = self._data_for_radial_profile[energy_id]["mdl"]
 
             # select the information only from the pixels that are within the radial
             # bin from origin of radial profile
 
-            this_data_tot = data[pixels_within_rad_bin].sum()
-            this_bkg_tot = bkg[pixels_within_rad_bin].sum()
-            this_model_tot = mdl[pixels_within_rad_bin].sum()
+            this_data_tot: float = data[pixels_within_rad_bin].sum()
+            this_bkg_tot: float = bkg[pixels_within_rad_bin].sum()
+            this_model_tot: float = mdl[pixels_within_rad_bin].sum()
 
             background[i] = this_bkg_tot
             observation[i] = this_data_tot
@@ -505,22 +536,22 @@ def _get_radial_profile(
         if active_planes is None:
             active_planes = self._active_planes
 
+        # call this here to get the data in order
+        self._prepare_data_for_radial_profile()
+
         # Make sure we use bins with data
         good_planes = [plane_id in active_planes for plane_id in self._active_planes]
         plane_ids = set(active_planes) & set(self._active_planes)
 
         offset = 0.50
         delta_r = (1.0 * max_radius) / n_radial_bins
-        radii = np.array([delta_r * (r + offset) for r in range(n_radial_bins)])
+        radii: ndarray = np.array([delta_r * (r + offset) for r in range(n_radial_bins)])
 
         # Get area of all pixels in a given circle
         # The area of each ring is then given by the difference between two
         # subsequent circe areas.
         area = np.array(
-            [
-                self._get_excess_background(ra, dec, r + offset * delta_r)[0]
-                for r in radii
-            ]
+            [self._get_excess_background(ra, dec, r + offset * delta_r)[0] for r in radii]
         )
 
         temp = area[1:] - area[:-1]
@@ -528,10 +559,7 @@ def _get_radial_profile(
 
         # signals
         signal = np.array(
-            [
-                self._get_excess_background(ra, dec, r + offset * delta_r)[1]
-                for r in radii
-            ]
+            [self._get_excess_background(ra, dec, r + offset * delta_r)[1] for r in radii]
         )
 
         temp = signal[1:] - signal[:-1]
@@ -539,10 +567,7 @@ def _get_radial_profile(
 
         # backgrounds
         bkg = np.array(
-            [
-                self._get_excess_background(ra, dec, r + offset * delta_r)[2]
-                for r in radii
-            ]
+            [self._get_excess_background(ra, dec, r + offset * delta_r)[2] for r in radii]
         )
 
         temp = bkg[1:] - bkg[:-1]
@@ -553,10 +578,7 @@ def _get_radial_profile(
         # model
         # convert 'top hat' excess into 'ring' excesses.
         model = np.array(
-            [
-                self._get_excess_background(ra, dec, r + offset * delta_r)[3]
-                for r in radii
-            ]
+            [self._get_excess_background(ra, dec, r + offset * delta_r)[3] for r in radii]
         )
 
         temp = model[1:] - model[:-1]
@@ -681,9 +703,7 @@ def plot_radial_profile(
 
         ax.plot(radii, excess_model, color="red", label="Model")
 
-        ax.legend(
-            bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16
-        )
+        ax.legend(bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16)
         ax.axhline(0, color="deepskyblue", linestyle="--")
 
         ax.set_xlim(left=0, right=max_radius)
@@ -773,9 +793,7 @@ def display_spectrum(self):
 
         yerr = [yerr_high, yerr_low]
 
-        return self._plot_spectrum(
-            net_counts, yerr, model_only, residuals, residuals_err
-        )
+        return self._plot_spectrum(net_counts, yerr, model_only, residuals, residuals_err)
 
     def _plot_spectrum(self, net_counts, yerr, model_only, residuals, residuals_err):
         fig, subs = plt.subplots(
@@ -1050,9 +1068,7 @@ def _get_expectation(
             )
 
             # Now multiply by the pixel area of the new map to go back to flux
-            this_model_map_hpx *= hp.nside2pixarea(
-                data_analysis_bin.nside, degrees=True
-            )
+            this_model_map_hpx *= hp.nside2pixarea(data_analysis_bin.nside, degrees=True)
 
         else:
             # No sources
@@ -1177,9 +1193,7 @@ def display_fit(self, smoothing_kernel_sigma=0.1, display_colorbar=False):
             subs[i][0].set_title("model, bin {}".format(data_analysis_bin.name))
 
             # Plot data map
-            images[1] = subs[i][1].imshow(
-                proj_data, origin="lower", vmin=vmin, vmax=vmax
-            )
+            images[1] = subs[i][1].imshow(proj_data, origin="lower", vmin=vmin, vmax=vmax)
             subs[i][1].set_title("excess, bin {}".format(data_analysis_bin.name))
 
             # Plot background map.
@@ -1299,9 +1313,7 @@ def _get_model_map(self, plane_id, n_pt_src, n_ext_src) -> SparseHealpix:
             raise ValueError(f"{plane_id} not a plane in the current model")
 
         model_map = SparseHealpix(
-            self._get_expectation(
-                self._maptree[plane_id], plane_id, n_pt_src, n_ext_src
-            ),
+            self._get_expectation(self._maptree[plane_id], plane_id, n_pt_src, n_ext_src),
             self._active_pixels[plane_id],
             self._maptree[plane_id].observation_map.nside,
         )
@@ -1374,17 +1386,13 @@ def _write_a_map(self, file_name, which, fluctuate=False, return_map=False):
         if return_map:
             return new_map_tree
 
-    def write_model_map(
-        self, file_name, poisson_fluctuate=False, test_return_map=False
-    ):
+    def write_model_map(self, file_name, poisson_fluctuate=False, test_return_map=False):
         """
         This function writes the model map to a file.
         The interface is based off of HAWCLike for consistency
         """
         if test_return_map:
-            log.warning(
-                "test_return_map=True should only be used for testing purposes!"
-            )
+            log.warning("test_return_map=True should only be used for testing purposes!")
         return self._write_a_map(file_name, "model", poisson_fluctuate, test_return_map)
 
     def write_residual_map(self, file_name, test_return_map=False):
@@ -1393,7 +1401,5 @@ def write_residual_map(self, file_name, test_return_map=False):
         The interface is based off of HAWCLike for consistency
         """
         if test_return_map:
-            log.warning(
-                "test_return_map=True should only be used for testing purposes!"
-            )
+            log.warning("test_return_map=True should only be used for testing purposes!")
         return self._write_a_map(file_name, "residual", False, test_return_map)

From 55e2452df0ddb41a113f8115c6d2e9ccf3c2bb14 Mon Sep 17 00:00:00 2001
From: Ramiro Torres-Escobedo <rtorrese@proton.me>
Date: Wed, 5 Feb 2025 22:18:24 +0800
Subject: [PATCH 5/7] cleaning the code a bit

---
 hawc_hal/HAL.py | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py
index 6a6b404..b2e3887 100644
--- a/hawc_hal/HAL.py
+++ b/hawc_hal/HAL.py
@@ -483,21 +483,12 @@ def _get_excess_background(
         area = np.full(len(self._active_planes), this_area)
 
         for i, energy_id in enumerate(self._active_planes):
-            # data_analysis_bin: DataAnalysisBin = self._maptree[energy_id]
-            # obtain the excess, background, and expected excess at
-            # each radial bin
-            # data: ndarray = data_analysis_bin.observation_map.as_partial()
-            # bkg: ndarray = data_analysis_bin.background_map.as_partial()
-            # mdl: ndarray = self._get_model_map(
-            #     energy_id, n_point_sources, n_ext_sources
-            # ).as_partial()
             data: ndarray = self._data_for_radial_profile[energy_id]["data"]
             bkg: ndarray = self._data_for_radial_profile[energy_id]["bkg"]
             mdl: ndarray = self._data_for_radial_profile[energy_id]["mdl"]
 
             # select the information only from the pixels that are within the radial
             # bin from origin of radial profile
-
             this_data_tot: float = data[pixels_within_rad_bin].sum()
             this_bkg_tot: float = bkg[pixels_within_rad_bin].sum()
             this_model_tot: float = mdl[pixels_within_rad_bin].sum()

From 01e072afea7da1ba11ba742a3b4a7cfaa04b6433 Mon Sep 17 00:00:00 2001
From: Ramiro Torres-Escobedo <rtorrese@proton.me>
Date: Wed, 5 Feb 2025 22:18:35 +0800
Subject: [PATCH 6/7] adding more type hinting information

---
 hawc_hal/HAL.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py
index b2e3887..0fcbda6 100644
--- a/hawc_hal/HAL.py
+++ b/hawc_hal/HAL.py
@@ -620,9 +620,11 @@ def _get_radial_profile(
         bkg = bkg[:, good_planes]
 
         # average over the analysis bins
-        excess_data = np.average(signal / area, weights=weight, axis=1)
-        excess_error = np.sqrt(np.sum(counts * weight * weight / (area * area), axis=1))
-        excess_model = np.average(model / area, weights=weight, axis=1)
+        excess_data: ndarray = np.average(signal / area, weights=weight, axis=1)
+        excess_error: ndarray = np.sqrt(
+            np.sum(counts * weight * weight / (area * area), axis=1)
+        )
+        excess_model: ndarray = np.average(model / area, weights=weight, axis=1)
 
         return (
             radii,

From 761f487f6f0472ad393d52bc7417febe55c22a97 Mon Sep 17 00:00:00 2001
From: Ramiro Torres-Escobedo <rtorrese@proton.me>
Date: Wed, 5 Feb 2025 22:50:13 +0800
Subject: [PATCH 7/7] more cleanup

---
 hawc_hal/HAL.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py
index 0fcbda6..e86b587 100644
--- a/hawc_hal/HAL.py
+++ b/hawc_hal/HAL.py
@@ -459,9 +459,6 @@ def _get_excess_background(
         signal = np.zeros_like(total_counts)
         area = np.zeros_like(total_counts)
 
-        n_point_sources = self._likelihood_model.get_number_of_point_sources()  # type: ignore
-        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()  # type: ignore
-
         longitude = ra_to_longitude(ra)
         latitude = dec
         center = hp.ang2vec(longitude, latitude, lonlat=True)