diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index d83fb15..e86b587 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 @@ -11,9 +12,11 @@ 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 +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): @@ -113,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() @@ -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 @@ -200,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() @@ -269,9 +277,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) @@ -286,9 +294,9 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non 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 @@ -385,28 +393,64 @@ 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 _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, ...]: + """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 """ - radius_radians = np.deg2rad(radius) + from math import radians + + if self._active_planes is None: + raise ValueError( + "No active planes have been set. Please use set_active_measurements()" + ) + + # 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) @@ -415,47 +459,36 @@ 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() - 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 - - radial_bin_pixels = hp.query_disc( - this_nside, center, radius_radians, inclusive=False - ) + # 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 + ) - # 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] + # select the pixels that are only within the radial bin + pixels_within_rad_bin = np.in1d( + self._active_pixels[self._active_planes[0]], radial_bin_pixels + ) - # 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] + # 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(len(self._active_planes), this_area) - # 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( - 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]) + for i, energy_id in enumerate(self._active_planes): + 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"] - this_data_tot = np.sum(bin_data) - this_bkg_tot = np.sum(bin_bkg) - this_model_tot = np.sum(bin_model) + # 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() background[i] = this_bkg_tot observation[i] = this_data_tot @@ -464,55 +497,49 @@ 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: 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.5 - delta_r = 1.0 * max_radius / n_radial_bins - radii = np.array([delta_r * (r + offset) for r in range(n_radial_bins)]) + offset = 0.50 + delta_r = (1.0 * max_radius) / 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] @@ -520,7 +547,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] @@ -528,7 +555,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] @@ -539,7 +566,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] @@ -551,7 +578,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,18 +600,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) - weight = np.array([w / np.sum(w) for _ in radii]) + weight = np.array([w / w.sum() for _ in radii]) # restrict profiles to the user-specified analysis bins area = area[:, good_planes] @@ -594,74 +617,71 @@ 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, excess_model, excess_data, excess_error, sorted(plane_ids) + return ( + radii, + excess_model, + excess_data, + excess_error, + sorted(plane_ids), # pyright:ignore + weight[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, 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( + {"Excess": excess_data, "Error": excess_error, "Model": excess_model}, + index=radii, + dtype=float, + ) df.index.name = "Radii" - df["Excess"] = excess_data - df["Bkg"] = 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, @@ -671,49 +691,32 @@ 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(bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16) - plt.axhline(0, color="deepskyblue", linestyle="--") + ax.legend(bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16) + ax.axhline(0, color="deepskyblue", linestyle="--") - x_limits = [0, max_radius] - plt.xlim(x_limits) - plt.xticks(fontsize=18) - plt.yticks(fontsize=18) - - 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): """ @@ -848,7 +851,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 @@ -1289,7 +1294,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 """ 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") 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"