Fixing energy fitting

Author: giacomov

hawc_hal/HAL.py

@@ -32,6 +32,10 @@ def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_sizes=0.17
         # Store ROI
         self._roi = roi
 
+        # Set up the flat-sky projection
+
+        self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_sizes)
+
         # Read map tree (data)
 
         self._maptree = map_tree_factory(maptree, roi=roi)
@@ -65,12 +69,10 @@ def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_sizes=0.17
         self._convolved_ext_sources = ConvolvedSourcesContainer()
 
         # By default all energy/nHit bins are used
-        self._all_planes = range(len(self._maptree))
-        self._active_planes = range(len(self._maptree))
-
-        # Set up the flat-sky projection
+        self._all_planes = list(self._maptree.analysis_bins_labels)
 
-        self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_sizes)
+        # The active planes list always contains the list of *indexes* of the active planes
+        self._active_planes_idx = range(len(self._maptree))
 
         # 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
@@ -132,7 +134,7 @@ def _compute_likelihood_biases(self):
             obs = data_analysis_bin.observation_map.as_partial()
             bkg = data_analysis_bin.background_map.as_partial()
 
-            sat_model = np.maximum(obs - bkg, 1e-30).astype(np.float64)
+            sat_model = np.clip(obs - bkg, 1e-50, None).astype(np.float64)
 
             self._saturated_model_like_per_maptree[i] = log_likelihood(obs, bkg, sat_model) - this_log_factorial
 
@@ -144,11 +146,45 @@ def get_saturated_model_likelihood(self):
         """
         return np.sum(self._saturated_model_like_per_maptree)
 
-    def set_active_measurements(self, bin_id_min, bin_id_max):
+    def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=None):
+
+        # 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"
+
+            # Make sure they are strings
+            bin_id_min = str(bin_id_min)
+            bin_id_max = str(bin_id_max)
+
+            assert bin_id_min in self._all_planes and bin_id_max in self._all_planes, "Illegal bin_names"
+
+            assert bin_list is None, "You can either provide a minimum and maximum bin, or a bin list, but not both"
 
-        assert bin_id_min in self._all_planes and bin_id_max in self._all_planes, "Illegal bin_name numbers"
+            idx1 = self._all_planes.index(bin_id_min)
+            idx2 = self._all_planes.index(bin_id_max)
 
-        self._active_planes = range(bin_id_min, bin_id_max + 1)
+            self._active_planes_idx = range(idx1, idx2 + 1)
+
+        else:
+
+            assert bin_id_max is None, "If you provide a maximum bin, you also need to provide a minimum bin"
+
+            assert bin_list is not None
+
+            self._active_planes_idx = []
+
+            for this_bin in bin_list:
+
+                try:
+
+                    this_idx = self._all_planes.index(str(this_bin))
+
+                except ValueError:
+
+                    raise ValueError("Bin %s it not contained in this response" % this_bin)
+
+                self._active_planes_idx.append(this_idx)
 
     def display(self):
 
@@ -179,7 +215,7 @@ def display(self):
         print("")
         print("Active energy/nHit planes: ")
         print("---------------------------\n")
-        print(self._active_planes)
+        print(self._active_planes_idx)
 
     def set_model(self, likelihood_model_instance):
         """
@@ -231,13 +267,13 @@ def display_spectrum(self):
         n_point_sources = self._likelihood_model.get_number_of_point_sources()
         n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
 
-        total_counts = np.zeros(len(self._active_planes), dtype=float)
+        total_counts = np.zeros(len(self._active_planes_idx), dtype=float)
         total_model = np.zeros_like(total_counts)
         model_only = np.zeros_like(total_counts)
         residuals = np.zeros_like(total_counts)
         net_counts = np.zeros_like(total_counts)
 
-        for i, energy_id in enumerate(self._active_planes):
+        for i, energy_id in enumerate(self._active_planes_idx):
 
             data_analysis_bin = self._maptree[energy_id]
 
@@ -270,11 +306,11 @@ def display_spectrum(self):
 
         fig, subs = plt.subplots(2, 1, gridspec_kw={'height_ratios': [2, 1], 'hspace': 0})
 
-        subs[0].errorbar(self._active_planes, net_counts, yerr=np.sqrt(total_counts),
+        subs[0].errorbar(self._active_planes_idx, net_counts, yerr=np.sqrt(total_counts),
                          capsize=0,
                          color='black', label='Net counts', fmt='.')
 
-        subs[0].plot(self._active_planes, model_only, label='Convolved model')
+        subs[0].plot(self._active_planes_idx, model_only, label='Convolved model')
 
         subs[0].legend(bbox_to_anchor=(1.0, 1.0), loc="upper right",
                        numpoints=1)
@@ -283,21 +319,21 @@ def display_spectrum(self):
         subs[1].axhline(0, linestyle='--')
 
         subs[1].errorbar(
-            self._active_planes, residuals,
+            self._active_planes_idx, residuals,
             yerr=np.ones(residuals.shape),
             capsize=0, fmt='.'
         )
 
-        x_limits = [min(self._active_planes) - 0.5, max(self._active_planes) + 0.5]
+        x_limits = [min(self._active_planes_idx) - 0.5, max(self._active_planes_idx) + 0.5]
 
         subs[0].set_yscale("log", nonposy='clip')
         subs[0].set_ylabel("Counts per bin")
         subs[0].set_xticks([])
 
         subs[1].set_xlabel("Analysis bin")
         subs[1].set_ylabel(r"$\frac{{cts - mod - bkg}}{\sqrt{mod + bkg}}$")
-        subs[1].set_xticks(self._active_planes)
-        subs[1].set_xticklabels(self._active_planes)
+        subs[1].set_xticks(self._active_planes_idx)
+        subs[1].set_xticklabels(self._active_planes_idx)
 
         subs[0].set_xlim(x_limits)
         subs[1].set_xlim(x_limits)
@@ -323,16 +359,19 @@ def get_log_like(self):
 
         for i, data_analysis_bin in enumerate(self._maptree):
 
-            if i not in self._active_planes:
+            if i not in self._active_planes_idx:
                 continue
 
             this_model_map_hpx = self._get_expectation(data_analysis_bin, i, n_point_sources, n_ext_sources)
 
             # Now compare with observation
             bkg_renorm = self._nuisance_parameters.values()[0].value
 
-            this_pseudo_log_like = log_likelihood(data_analysis_bin.observation_map.as_partial(),
-                                                  data_analysis_bin.background_map.as_partial() * bkg_renorm,
+            obs = data_analysis_bin.observation_map.as_partial()  # type: np.array
+            bkg = data_analysis_bin.background_map.as_partial() * bkg_renorm  # type: np.array
+
+            this_pseudo_log_like = log_likelihood(obs,
+                                                  bkg,
                                                   this_model_map_hpx)
 
             total_log_like += this_pseudo_log_like - self._log_factorials[i] - self._saturated_model_like_per_maptree[i]
@@ -364,7 +403,7 @@ def get_simulated_dataset(self, name):
 
             for i, data_analysis_bin in enumerate(self._maptree):
 
-                if i not in self._active_planes:
+                if i not in self._active_planes_idx:
 
                     expectations.append(None)
 
@@ -389,7 +428,7 @@ def get_simulated_dataset(self, name):
         # Substitute the observation and background for each data analysis bin
         for i, (data_analysis_bin, orig_data_analysis_bin) in enumerate(zip(self._clone[0]._maptree, self._maptree)):
 
-            if i not in self._active_planes:
+            if i not in self._active_planes_idx:
 
                 continue
 
@@ -524,13 +563,13 @@ def display_fit(self, smoothing_kernel_sigma=0.1):
         # The image is going to cover the diameter plus 20% padding
         xsize = self._get_optimal_xsize(resolution)
 
-        n_active_planes = len(self._active_planes)
+        n_active_planes = len(self._active_planes_idx)
 
         fig, subs = plt.subplots(n_active_planes, 3, figsize=(8, n_active_planes * 2))
 
-        with progress_bar(len(self._active_planes), title='Smoothing maps') as prog_bar:
+        with progress_bar(len(self._active_planes_idx), title='Smoothing maps') as prog_bar:
 
-            for i, plane_id in enumerate(self._active_planes):
+            for i, plane_id in enumerate(self._active_planes_idx):
 
                 data_analysis_bin = self._maptree[plane_id]
 
@@ -604,7 +643,7 @@ def display_stacked_image(self, smoothing_kernel_sigma=0.5):
         # The image is going to cover the diameter plus 20% padding
         xsize = self._get_optimal_xsize(resolution)
 
-        active_planes_bins = map(lambda x: self._maptree[x], self._active_planes)
+        active_planes_bins = map(lambda x: self._maptree[x], self._active_planes_idx)
 
         # Get the center of the projection for this plane
         this_ra, this_dec = self._roi.ra_dec_center

hawc_hal/convenience_functions/__init__.py

@@ -0,0 +1,72 @@
+from ..HAL import HAL
+
+from threeML import DataList, JointLikelihood
+
+
+def fit_point_source(roi,
+                     maptree,
+                     response,
+                     point_source_model,
+                     bin_list,
+                     confidence_intervals=False,
+                     liff=False,
+                     pixel_size=0.17,
+                     verbose=False):
+
+    data_radius = roi.data_radius.to("deg").value
+
+    if not liff:
+
+        # This is a 3ML plugin
+        hawc = HAL("HAWC",
+                   maptree,
+                   response,
+                   roi,
+                   flat_sky_pixels_sizes=pixel_size)
+
+        hawc.set_active_measurements(bin_list=bin_list)
+
+    else:
+
+        from threeML import HAWCLike
+
+        hawc = HAWCLike("HAWC",
+                        maptree,
+                        response,
+                        fullsky=True)
+
+        hawc.set_bin_list(bin_list)
+
+        ra_roi, dec_roi = roi.ra_dec_center
+
+        hawc.set_ROI(ra_roi, dec_roi, data_radius)
+
+    if not liff:
+
+        hawc.display()
+
+    data = DataList(hawc)
+
+    jl = JointLikelihood(point_source_model, data, verbose=verbose)
+
+    point_source_model.display(complete=True)
+
+    try:
+
+        jl.set_minimizer("minuit")
+
+    except:
+
+        jl.set_minimizer("minuit")
+
+    param_df, like_df = jl.fit()
+
+    if confidence_intervals:
+
+        ci = jl.get_errors()
+
+    else:
+
+        ci = None
+
+    return param_df, like_df, ci, jl.results

hawc_hal/convenience_functions/fit_point_source.py

@@ -8,7 +8,7 @@ def _select_with_wrap_around(arr, start, stop, wrap=(360, 0)):
     if start > stop:
 
         # This can happen if for instance lon_start = 280 and lon_stop = 10, which means we
-        # should keep between 280 and 360 and then between 360 to 10
+        # should keep between 280 and 360 and then between 0 to 10
         idx = ((arr >= stop) & (arr <= wrap[0])) | ((arr >= wrap[1]) & (arr <= start))
 
     else:

hawc_hal/convolved_source/convolved_extended_source.py

@@ -1,8 +1,9 @@
 import numpy as np
+import os
 
 from astromodels import PointSource
 from ..psf_fast import PSFInterpolator
-
+from ..interpolation.log_log_interpolation import LogLogInterpolator
 
 from threeML.exceptions.custom_exceptions import custom_warnings
 
@@ -47,11 +48,7 @@ def _update_dec_bins(self, dec_src):
                                                                                self._flat_sky_projection),
                                           self._response_energy_bins)
 
-        # for i, psf_i in enumerate(self._psf_interpolators):
-        #
-        #     print("PSF %i: %.3f, %.3f" % (i, psf_i._psf.truncation_radius, psf_i._psf.kernel_radius))
-
-    def get_source_map(self, response_bin_id, tag=None):
+    def get_source_map(self, response_bin_id, tag=None, integrate=False):
 
         # Get current point source position
         # NOTE: this might change if the point source position is free during the fit,
@@ -75,14 +72,16 @@ def get_source_map(self, response_bin_id, tag=None):
         # is very fast
         this_map = psf_interpolator.point_source_image(ra_src, dec_src)
 
-        # Check that the point source is contained in the ROI, if not print a warning
+        # Check that the point source image is entirely contained in the ROI, if not print a warning
         map_sum = this_map.sum()
+
         if not np.isclose(map_sum, 1.0, rtol=1e-2):
 
             custom_warnings.warn("PSF for source %s is not entirely contained "
-                                 "in ROI for response bin %s. Fraction is %.2f instead of 1.0" % (self._name,
-                                                                                                  response_bin_id,
-                                                                                                  map_sum))
+                                 "in ROI for response bin %s. Fraction is %.2f instead of 1.0. "
+                                 "Consider enlarging your model ROI." % (self._name,
+                                                                         response_bin_id,
+                                                                         map_sum))
 
         # Compute the fluxes from the spectral function at the same energies as the simulated function
         energy_centers_keV = response_energy_bin.sim_energy_bin_centers * 1e9  # astromodels expects energies in keV
@@ -92,11 +91,37 @@ def get_source_map(self, response_bin_id, tag=None):
 
         source_diff_spectrum = self._source(energy_centers_keV, tag=tag)
 
-        # Transform from keV^-1 cm^-2 s^-1 to TeV^-1 cm^-2 s^-1
-        source_diff_spectrum *= 1e9
+        # This provide a way to force the use of integration, for testing
+        if os.environ.get("HAL_INTEGRATE_POINT_SOURCE", '').lower() == 'yes':  # pragma: no cover
+
+            integrate = True
+
+        if integrate:  # pragma: no cover
+
+            # Slower approach, where we integrate the spectra of both the simulation and the source
+            # It is off by default because it is slower and it doesn't provide any improvement in accuracy
+            # according to my simulations (GV)
+
+            interp_spectrum = LogLogInterpolator(response_energy_bin.sim_energy_bin_centers,
+                                                 source_diff_spectrum * 1e9,
+                                                 k=2)
+
+            interp_sim_spectrum = LogLogInterpolator(response_energy_bin.sim_energy_bin_centers,
+                                                     response_energy_bin.sim_differential_photon_fluxes,
+                                                     k=2)
+
+            src_spectrum = [interp_spectrum.integral(a, b) for (a, b) in zip(response_energy_bin.sim_energy_bin_low,
+                                                                             response_energy_bin.sim_energy_bin_hi)]
+
+            sim_spectrum = [interp_sim_spectrum.integral(a, b) for (a, b) in zip(response_energy_bin.sim_energy_bin_low,
+                                                                                 response_energy_bin.sim_energy_bin_hi)]
+
+            scale = np.array(src_spectrum) / np.array(sim_spectrum)
+
+        else:
 
-        # Re-weight the detected counts
-        scale = source_diff_spectrum / response_energy_bin.sim_differential_photon_fluxes
+            # Transform from keV^-1 cm^-2 s^-1 to TeV^-1 cm^-2 s^-1 and re-weight the detected counts
+            scale = source_diff_spectrum / response_energy_bin.sim_differential_photon_fluxes * 1e9
 
         # Now return the map multiplied by the scale factor
         return np.sum(scale * response_energy_bin.sim_signal_events_per_bin) * this_map