diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index d83fb15..5bcc2d1 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -354,26 +354,30 @@ def set_model(self, likelihood_model_instance): self._convolved_point_sources.reset() self._convolved_ext_sources.reset() + self._setup_psf_convolutors() + # For each point source in the model, build the convolution class + pnt_sources = self._likelihood_model.point_sources.values() - for source in list(self._likelihood_model.point_sources.values()): - this_convolved_point_source = ConvolvedPointSource( - source, self._response, self._flat_sky_projection - ) + if pnt_sources: + + for source in list(pnt_sources): + this_convolved_point_source = ConvolvedPointSource( + source, self._response, self._flat_sky_projection, self._active_planes + ) + + self._convolved_point_sources.append(this_convolved_point_source) - self._convolved_point_sources.append(this_convolved_point_source) - # Samewise for extended sources + # Likewise for extended sources ext_sources = list(self._likelihood_model.extended_sources.values()) # NOTE: ext_sources evaluate to False if empty if ext_sources: - # We will need to convolve - - self._setup_psf_convolutors() for source in ext_sources: if source.spatial_shape.n_dim == 2: + this_convolved_ext_source = ConvolvedExtendedSource2D( source, self._response, self._flat_sky_projection ) @@ -861,6 +865,7 @@ def get_log_like(self, individual_bins=False, return_null=False): this_model_map_hpx = self._get_expectation( data_analysis_bin, bin_id, n_point_sources, n_ext_sources ) + # Now compare with observation bkg_renorm = list(self._nuisance_parameters.values())[0].value @@ -921,20 +926,13 @@ def get_simulated_dataset(self, name): data_analysis_bin = self._maptree[bin_id] if bin_id not in self._active_planes: expectations[bin_id] = None - else: - expectations[bin_id] = ( - self._get_expectation( - data_analysis_bin, bin_id, n_point_sources, n_ext_sources - ) - + data_analysis_bin.background_map.as_partial() - ) + expectations[bin_id] = ( self._get_expectation( data_analysis_bin, bin_id, n_point_sources, n_ext_sources ) + + data_analysis_bin.background_map.as_partial() ) if parallel_client.is_parallel_computation_active(): # Do not clone, as the parallel environment already makes clones - clone = self - else: clone = copy.deepcopy(self) @@ -946,13 +944,10 @@ def get_simulated_dataset(self, name): if bin_id not in self._active_planes: continue - else: # Active plane. Generate new data expectation = self._clone[1][bin_id] - new_data = np.random.poisson( - expectation, size=(1, expectation.shape[0]) - ).flatten() + new_data = np.random.poisson( expectation, size=(1, expectation.shape[0]) ).flatten() # Substitute data data_analysis_bin.observation_map.set_new_values(new_data) @@ -962,95 +957,80 @@ def get_simulated_dataset(self, name): # Adjust the name of the nuisance parameter old_name = list(self._clone[0]._nuisance_parameters.keys())[0] new_name = old_name.replace(self.name, name) - self._clone[0]._nuisance_parameters[new_name] = self._clone[ - 0 - ]._nuisance_parameters.pop(old_name) + self._clone[0]._nuisance_parameters[new_name] = self._clone[0]._nuisance_parameters.pop(old_name) # Recompute biases self._clone[0]._compute_likelihood_biases() return self._clone[0] - def _get_expectation( - self, data_analysis_bin, energy_bin_id, n_point_sources, n_ext_sources - ): + def _get_expectation( self, data_analysis_bin, energy_bin_id, n_point_sources, n_ext_sources ): # Compute the expectation from the model - this_model_map = None - for pts_id in range(n_point_sources): - this_conv_pnt_src: ConvolvedPointSource = self._convolved_point_sources[ - pts_id - ] + #If there are sources in the map + if n_point_sources > 0 or n_ext_sources > 0: - expectation_per_transit = this_conv_pnt_src.get_source_map( - energy_bin_id, - tag=None, - psf_integration_method=self._psf_integration_method, - ) + this_unconv_model_map = None + this_conv_ps_model_map = None + this_conv_ex_model_map = None - expectation_from_this_source = ( - expectation_per_transit * data_analysis_bin.n_transits - ) + #Process point sources + if n_point_sources > 0: - if this_model_map is None: - # First addition + for src_id in range(n_point_sources): - this_model_map = expectation_from_this_source + this_conv_src: ConvolvedPointSource = self._convolved_point_sources[ src_id ] - else: - this_model_map += expectation_from_this_source + expectation_per_transit = this_conv_src.get_source_map( energy_bin_id ) - # Now process extended sources - if n_ext_sources > 0: - this_ext_model_map = None + if this_conv_ps_model_map is None: + # First addition + this_conv_ps_model_map = expectation_per_transit + else: + this_conv_ps_model_map += expectation_per_transit - for ext_id in range(n_ext_sources): - this_conv_src: Union[ - ConvolvedExtendedSource2D, ConvolvedExtendedSource3D - ] = self._convolved_ext_sources[ext_id] + #Process extended sources + if n_ext_sources > 0: - expectation_per_transit = this_conv_src.get_source_map(energy_bin_id) + for ext_id in range(n_ext_sources): + this_conv_src: Union[ + ConvolvedExtendedSource2D, ConvolvedExtendedSource3D + ] = self._convolved_ext_sources[ ext_id ] - if this_ext_model_map is None: - # First addition + expectation_per_transit = this_conv_src.get_source_map( energy_bin_id ) - this_ext_model_map = expectation_per_transit + if this_unconv_model_map is None: + # First addition + this_unconv_model_map = expectation_per_transit + else: + this_unconv_model_map += expectation_per_transit - else: - this_ext_model_map += expectation_per_transit + # Now convolve with the PSF + this_conv_ex_model_map = self._psf_convolutors[ energy_bin_id ].extended_source_image( this_unconv_model_map ) - # Now convolve with the PSF - if this_model_map is None: - # Only extended sources - this_model_map = ( - self._psf_convolutors[energy_bin_id].extended_source_image( - this_ext_model_map - ) - * data_analysis_bin.n_transits - ) - else: - this_model_map += ( - self._psf_convolutors[energy_bin_id].extended_source_image( - this_ext_model_map - ) - * data_analysis_bin.n_transits - ) + if n_point_sources > 0 and n_ext_sources > 0: + + this_model_map = ( this_conv_ps_model_map + this_conv_ex_model_map ) * data_analysis_bin.n_transits + + elif n_point_sources > 0: + + this_model_map = this_conv_ps_model_map * data_analysis_bin.n_transits + + elif n_ext_sources > 0: + + this_model_map = this_conv_ex_model_map * data_analysis_bin.n_transits + # Now transform from the flat sky projection to HEALPiX if this_model_map is not None: # First divide for the pixel area because we need to interpolate brightness - # this_model_map = old_div(this_model_map, self._flat_sky_projection.project_plane_pixel_area) - this_model_map = ( - this_model_map / self._flat_sky_projection.project_plane_pixel_area - ) + this_model_map = ( this_model_map / self._flat_sky_projection.project_plane_pixel_area ) - this_model_map_hpx = self._flat_sky_to_healpix_transform[energy_bin_id]( - this_model_map, fill_value=0.0 - ) + this_model_map_hpx = self._flat_sky_to_healpix_transform[energy_bin_id]( this_model_map, fill_value=0.0 ) # 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) @@ -1061,14 +1041,10 @@ def _get_expectation( this_model_map_hpx = 0.0 return this_model_map_hpx - @staticmethod - def _represent_healpix_map( - fig, hpx_map, longitude, latitude, xsize, resolution, smoothing_kernel_sigma - ): - proj = get_gnomonic_projection( - fig, hpx_map, rot=(longitude, latitude, 0.0), xsize=xsize, reso=resolution - ) + def _represent_healpix_map( fig, hpx_map, longitude, latitude, xsize, resolution, smoothing_kernel_sigma ): + + proj = get_gnomonic_projection( fig, hpx_map, rot=(longitude, latitude, 0.0), xsize=xsize, reso=resolution ) if smoothing_kernel_sigma is not None: # Get the sigma in pixels diff --git a/hawc_hal/convolved_source/convolved_point_source.py b/hawc_hal/convolved_source/convolved_point_source.py index 219a640..4936d53 100644 --- a/hawc_hal/convolved_source/convolved_point_source.py +++ b/hawc_hal/convolved_source/convolved_point_source.py @@ -1,135 +1,348 @@ from __future__ import division from builtins import zip from builtins import object -from past.utils import old_div import os import collections import numpy as np -from astromodels import PointSource -from ..psf_fast import PSFInterpolator -from ..interpolation.log_log_interpolation import LogLogInterpolator +from hawc_hal.psf_fast import PSFInterpolator +from hawc_hal import HAL +from astromodels import use_astromodels_memoization +from astromodels.utils.angular_distance import angular_distance from threeML.io.logging import setup_logger log = setup_logger(__name__) log.propagate = False +from scipy.ndimage import shift + class ConvolvedPointSource(object): - def __init__(self, source, response, flat_sky_projection): + def __init__(self, source, response, flat_sky_projection, _active_planes): - assert isinstance(source, PointSource) + self._response = response + self._flat_sky_projection = flat_sky_projection + self._active_planes = _active_planes + # Get name + self._name = source.name self._source = source - # Get name - self._name = self._source.name + # Get the RA and Dec of our source + lon = self._source.position.ra.value + lat = self._source.position.dec.value + + # Get the defined dec bins lower edges + self._lower_edges = np.array([x[0] for x in response.dec_bins]) + self._upper_edges = np.array([x[-1] for x in response.dec_bins]) + self._centers = np.array([x[1] for x in response.dec_bins]) + + #Now find the dec bins we need to consider and pull out the psf's + if (self._source.position.dec.free): + dec_min = self._source.position.dec.min_value + dec_max = self._source.position.dec.max_value + else: + #if we have a fixed source + dec_min = dec_max = lat + + # Find the dec bins within which the source can move + self._dec_bins_to_consider_idx = np.flatnonzero( (self._upper_edges >= dec_min) & (self._lower_edges <= dec_max) ) + # Add one dec bin to cover the last part + if (self._centers[self._dec_bins_to_consider_idx[-1]] < dec_max): + # If not in very last bin + if ( self._dec_bins_to_consider_idx[-1] != len(self._centers) - 1 ): + self._dec_bins_to_consider_idx = np.append(self._dec_bins_to_consider_idx, [self._dec_bins_to_consider_idx[-1] + 1]) + # Add one dec bin to cover the first part + if (self._centers[self._dec_bins_to_consider_idx[0]] > dec_min): + # If not in very first bin + if ( self._dec_bins_to_consider_idx[0] != 0 ): + self._dec_bins_to_consider_idx = np.insert(self._dec_bins_to_consider_idx, 0, [self._dec_bins_to_consider_idx[0] - 1]) + + self._dec_bins_to_consider = [self._response.response_bins[self._centers[x]] for x in self._dec_bins_to_consider_idx] + + log.info("Considering %i dec bins for point source %s" % (len(self._dec_bins_to_consider), self._name)) + log.info("Dec bins are %s" % (self._dec_bins_to_consider_idx) ) + + self._decbin_to_idx = {} + for i, decbin in enumerate(self._dec_bins_to_consider_idx): + self._decbin_to_idx[ decbin ] = i + + # Prepare somewhere to save convoluters and convolved images for all dec bins and energy bins + # Convoluters + self._psf_convolutors = [] + for response_bin in self._dec_bins_to_consider: + psf_convolutor = collections.OrderedDict() + for energy_bin_id in self._active_planes: + psf_convolutor[ energy_bin_id ] = PSFInterpolator( response_bin[ energy_bin_id ].psf, self._flat_sky_projection ) + self._psf_convolutors.append( psf_convolutor ) - self._response = response + # This stores the convoluted sources + self._this_model_image_norm_conv = [] + for response_bin in self._dec_bins_to_consider: + this_model_image_norm_conv = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + self._this_model_image_norm_conv.append( this_model_image_norm_conv ) - self._flat_sky_projection = flat_sky_projection + # Store weights and Dec bin indices for current position + self._weights = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + self._dec_bins_idx = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + + # This stores the convoluted sources + self._this_model_image_norm_conv_shifted = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + + # Get the initial skymask and energies + self._energy_centers_keV = self._update_energy_centers( lat ) + + # Prepare array for fluxes + self._all_fluxes = np.zeros( self._energy_centers_keV.shape[0] ) + + # This stores the position that was processed during the last step, store a dummy value here for the time being + self._last_processed_position = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + for key, value in self._last_processed_position.items(): + self._last_processed_position[ key ] = (0., 0.) + + self._last_processed_position_idx = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + for key, value in self._last_processed_position_idx.items(): + self._last_processed_position_idx[ key ] = [0, 0] + + # We implement a caching system so that the source flux is evaluated only when strictly needed, + # because it is the most computationally intense part otherwise. + self._recompute = True + + # Register callback to keep track of the parameter changes + self._setup_callbacks(self._parameter_change_callback) + + # This is a simple counter that counts the number of iterations in each energy bin + self.new_bool = collections.OrderedDict( dict.fromkeys( self._active_planes ) ) + for key, value in self.new_bool.items(): + self.new_bool[ key ] = True + + + def _get_relevant_bins( self, dec_current ): + + dec_bins_idx = np.flatnonzero((self._upper_edges >= dec_current) & (self._lower_edges <= dec_current)) + + # Add one dec bin to cover the last part + if (self._centers[dec_bins_idx[-1]] < dec_current): + # If not in very last bin + if ( dec_bins_idx[-1] != len(self._centers) - 1 ): + dec_bins_idx = np.append(dec_bins_idx, [dec_bins_idx[-1] + 1]) + else: + dec_bins_idx = np.append(dec_bins_idx, [dec_bins_idx[-1]]) + + # Add one dec bin to cover the first part + if (self._centers[dec_bins_idx[0]] > dec_current): + # If not in very first bin + if ( dec_bins_idx[0] != 0 ): + dec_bins_idx = np.insert(dec_bins_idx, 0, [dec_bins_idx[0] - 1]) + else: + dec_bins_idx = np.insert(dec_bins_idx, 0, [dec_bins_idx[0]]) + + assert dec_bins_idx[0] <= dec_bins_idx[1], "Dec bins are in the wrong order!" + + return dec_bins_idx + + + def _get_weights( self, dec_bins_idx, energy_bin_id ): + + dec_bin1, dec_bin2 = self._dec_bins_to_consider[ + self._decbin_to_idx[ dec_bins_idx[0] ] ], self._dec_bins_to_consider[ + self._decbin_to_idx[ dec_bins_idx[1] ] ] + + # Get the two response bins to consider + this_response_bin1, this_response_bin2 = dec_bin1[ energy_bin_id ], dec_bin2[ energy_bin_id ] + + # Figure out which pixels are between the centers of the dec bins we are considering + c1, c2 = this_response_bin1.declination_center, this_response_bin2.declination_center + + # Compute the interpolation weights for the two responses + w1 = (self._pix_ctr_coords[1] - c2) / (c1 - c2) + w2 = (self._pix_ctr_coords[1] - c1) / (c2 - c1) + + return np.array( [w1, w2 ] ) + + + def _update_energy_centers( self, dec_current): + + # Find central bin for the PSF + self._central_response_bins = self._response.get_response_dec_bin( dec_current, interpolate=False ) + + # Get the energies needed for the computation of the flux + energy_centers_keV = self._central_response_bins[list(self._central_response_bins.keys())[0]].sim_energy_bin_centers * 1e9 + + return energy_centers_keV + + + def _parameter_change_callback(self, this_parameter): + + # A parameter has changed, we need to recompute the function. + # NOTE: we do not recompute it here because if more than one parameter changes at the time (like for example + # during sampling) we do not want to recompute the function multiple time before we get to the convolution + # stage. Therefore we will compute the function in the get_source_map method + + #This is a general callback - any parameter changes this gets set to True + + # print("%s has changed" % this_parameter.name) + self._recompute = True + + + def _setup_callbacks(self, callback): + + # Register call back with all free parameters and all linked parameters. If a parameter is linked to another + # one, the other one might be in a different source, so we add a callback to that one so that we recompute + # this source when needed. + for parameter in list(self._source.parameters.values()): + + if parameter.free: + parameter.add_callback(callback) + + if parameter.has_auxiliary_variable: + # Add a callback to the auxiliary variable so that when that is changed, we need + # to recompute the model + aux_variable, _ = parameter.auxiliary_variable + + aux_variable.add_callback(callback) + + + def get_source_map(self, energy_bin_id, tag=None ): + + # We do not use the memoization in astromodels because we are doing caching by ourselves, + # so the astromodels memoization would turn into 100% cache miss and use a lot of RAM for nothing, + # given that we are evaluating the function on many points and many energies + with use_astromodels_memoization(False): + + ra_current, dec_current = self._source.position.ra.value, self._source.position.dec.value + + if self.new_bool[ energy_bin_id ]: + + # Only need to do this once + self.new_bool[ energy_bin_id ] = False + + # Using one of the ordered dicts to check whether this is the first energy bin in the collection as only need to + # calc some things once for each iteration across all bins + if ( energy_bin_id == self._active_planes[0] ): + + # Find the index in RA and Dec of the pixel containing the source and the pixel centre coordinates + self._idx = self._flat_sky_projection.wcs.world_to_array_index_values( [[ ra_current, dec_current ]] )[0] + self._pix_ctr_coords = self._flat_sky_projection.wcs.array_index_to_world_values( [self._idx] )[0] + + self._deltaidx = ( self._flat_sky_projection.wcs.world_to_pixel_values( [[ ra_current, dec_current ]] )[0] - + self._flat_sky_projection.wcs.world_to_pixel_values( [self._pix_ctr_coords] )[0] ) - # This will store the position of the source - # right now, let's use a fake value so that at the first iteration the source maps will be filled - # (see get_expected_signal_per_transit) - self._last_processed_position = (-999, -999) + # Calculate the PSF in each dec bin for the energy bin, we only need to do this once + for i in range( len( self._dec_bins_to_consider ) ): - self._response_energy_bins = None - self._psf_interpolators = None + self._this_model_image_norm_conv[ i ][ energy_bin_id ] = self._psf_convolutors[ i ][ + energy_bin_id ].point_source_image( self._pix_ctr_coords[0], self._pix_ctr_coords[1], "fast") - @property - def name(self): + # The convoluter returns negative pixel values in some cases, this needs to be fixed really, we just botch it here + # We don't just truncate at 0. but shift the image, just in case there is gradient information in there that the + # minimizer might need. +# if ( np.any(self._this_model_image_norm_conv[ i ][ energy_bin_id ] < 0.) ): + self._this_model_image_norm_conv[ i ][ energy_bin_id ] = ( + self._this_model_image_norm_conv[ i ][ energy_bin_id ] + + np.abs( np.min( self._this_model_image_norm_conv[ i ][ energy_bin_id ] ) ) ) + self._this_model_image_norm_conv[ i ][ energy_bin_id ] = ( + self._this_model_image_norm_conv[ i ][ energy_bin_id ].clip(min=0.) / + np.sum( self._this_model_image_norm_conv[ i ][ energy_bin_id ].clip(min=0.) ) ) - return self._name - def _update_dec_bins(self, dec_src): + # Get the two decbins that relevant for this points source convolution + self._dec_bins_idx[ energy_bin_id ] = self._get_relevant_bins( dec_current ) - if abs(dec_src - self._last_processed_position[1]) > 0.1: + # Calculate the weight between the two bins + self._weights[ energy_bin_id ] = self._get_weights( self._dec_bins_idx[ energy_bin_id ], energy_bin_id ) - # Source moved by more than 0.1 deg, let's recompute the response and the PSF - self._response_energy_bins = self._response.get_response_dec_bin(dec_src, interpolate=True) + # Calculate the weighted sum of the convoluted images + this_model_image_norm_conv_weighted = ( + self._weights[ energy_bin_id ][0] * self._this_model_image_norm_conv[ + self._decbin_to_idx[ self._dec_bins_idx[ energy_bin_id ][0] ] ][ energy_bin_id ] + + self._weights[ energy_bin_id ][1] * self._this_model_image_norm_conv[ + self._decbin_to_idx[ self._dec_bins_idx[ energy_bin_id ][1] ] ][ energy_bin_id ] ) - # Setup the PSF interpolators - self._psf_interpolators = collections.OrderedDict() - for bin_id in self._response_energy_bins: - self._psf_interpolators[bin_id] = PSFInterpolator(self._response_energy_bins[bin_id].psf, - self._flat_sky_projection) - def get_source_map(self, response_bin_id, tag=None, integrate=False, psf_integration_method='fast'): + self._this_model_image_norm_conv_shifted[ energy_bin_id ] = shift( this_model_image_norm_conv_weighted, + self._deltaidx[::-1], order=2, mode='grid-constant', cval=0.0, prefilter=True ) - # Get current point source position - # NOTE: this might change if the point source position is free during the fit, - # that's why it is here + this_model_image_norm_conv = self._this_model_image_norm_conv_shifted[ energy_bin_id ] - ra_src, dec_src = self._source.position.ra.value, self._source.position.dec.value + # Store the last processed position + self._last_processed_position[ energy_bin_id ] = ( ra_current, dec_current ) - if (ra_src, dec_src) != self._last_processed_position: + # Store the indices of the position in the RA/Dec grid + self._last_processed_position_idx[ energy_bin_id ] = self._idx - # Position changed (or first iteration), let's update the dec bins - self._update_dec_bins(dec_src) + # same position so don't need to update the convolution + elif ( ( ra_current, dec_current ) == self._last_processed_position[ energy_bin_id ] ): - self._last_processed_position = (ra_src, dec_src) + this_model_image_norm_conv = self._this_model_image_norm_conv_shifted[ energy_bin_id ] - # Get the current response bin - response_energy_bin = self._response_energy_bins[response_bin_id] - psf_interpolator = self._psf_interpolators[response_bin_id] + # Position has changed, so we just shift the convolved source in the image. We can do this at the moment because we are + # using only a single psf in the ROI, to improve this, for large ROIs we should really update the PSF and run the first + # item in the structure + else: - # Get the PSF image - # This is cached inside the PSF class, so that if the position doesn't change this line - # is very fast - this_map = psf_interpolator.point_source_image(ra_src, dec_src, psf_integration_method) + if ( energy_bin_id == self._active_planes[0] ): - # Check that the point source image is entirely contained in the ROI, if not print a warning - map_sum = this_map.sum() + # Find the new pixel and its centre coordinates + self._idx = self._flat_sky_projection.wcs.world_to_array_index_values( [[ ra_current, dec_current ]] )[0] + self._pix_ctr_coords = self._flat_sky_projection.wcs.array_index_to_world_values( [self._idx] )[0] - if not np.isclose(map_sum, 1.0, rtol=1e-2): + # Determine the difference in position and save this for the other energy bins + self._deltaidx = ( self._flat_sky_projection.wcs.world_to_pixel_values( [[ ra_current, dec_current ]] )[0] - + self._last_processed_position_idx[ energy_bin_id ] ) - log.warning("PSF for source %s is not entirely contained " - "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)) + # get the two Dec bins between which to interpolate the PSF, get the weights and calculate the interpolates image + self._dec_bins_idx[ energy_bin_id ] = self._get_relevant_bins( dec_current ) - # 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 + self._weights[ energy_bin_id ] = self._get_weights( self._dec_bins_idx[ energy_bin_id ], energy_bin_id ) - # This call needs to be here because the parameters of the model might change, - # for example during a fit + this_model_image_norm_conv_weighted = ( + self._weights[ energy_bin_id ][0] * self._this_model_image_norm_conv[ + self._decbin_to_idx[ self._dec_bins_idx[ energy_bin_id ][0] ] ][ energy_bin_id ] + + self._weights[ energy_bin_id ][1] * self._this_model_image_norm_conv[ + self._decbin_to_idx[ self._dec_bins_idx[ energy_bin_id ][1] ] ][ energy_bin_id ] ) - source_diff_spectrum = self._source(energy_centers_keV, tag=tag) + # Shift the image with sub-pixel precision + this_model_image_norm_conv = shift( this_model_image_norm_conv_weighted, + self._deltaidx[::-1], order=2, mode='grid-constant', cval=0.0, prefilter=True ) - # 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 + # If we need to recompute the flux, let's do it + if self._recompute: - integrate = True + #Check that dec is still between decbin centres, if not need to calculate dec bins again + if ( dec_current != self._last_processed_position[ energy_bin_id ][1] ): - if integrate: # pragma: no cover + self._energy_centers_keV = self._update_energy_centers( dec_current ) - # 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) + # Recompute the fluxes for the pixels that are covered by this source, we do this no matter what + self._all_fluxes = self._source( self._energy_centers_keV, tag=tag ) + # 1 / (keV cm^2 s rad^2) - interp_spectrum = LogLogInterpolator(response_energy_bin.sim_energy_bin_centers, - source_diff_spectrum * 1e9, - k=2) + # We don't need to recompute the function anymore until a parameter changes so reset the callback to False + self._recompute = False - 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)] + # Loop over the Dec bins that cover this source and compute the expected flux, interpolating between + # two dec bins for each point + dec_bin1, dec_bin2 = self._dec_bins_to_consider[ + self._decbin_to_idx[ self._dec_bins_idx[ energy_bin_id ][0] ] ], self._dec_bins_to_consider[ + self._decbin_to_idx[ self._dec_bins_idx[ energy_bin_id ][1] ] ] - 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)] + # Get the two response bins to consider + this_response_bin1, this_response_bin2 = dec_bin1[energy_bin_id], dec_bin2[energy_bin_id] - scale = old_div(np.array(src_spectrum), np.array(sim_spectrum)) + # Reweight the spectrum separately for the two bins + # NOTE: the scale is the same because the sim_differential_photon_fluxes are the same (the simulation + # used to make the response used the same spectrum for each bin). What changes between the two bins + # is the observed signal per bin (the .sim_signal_events_per_bin member) + scale = self._all_fluxes / this_response_bin1.sim_differential_photon_fluxes - else: + # This scales the unit PSF image to the actual flux + factor = (self._weights[ energy_bin_id ][0] * np.sum(scale * this_response_bin1.sim_signal_events_per_bin) + + self._weights[ energy_bin_id ][1] * np.sum(scale * this_response_bin2.sim_signal_events_per_bin)) * 1e9 - # Transform from keV^-1 cm^-2 s^-1 to TeV^-1 cm^-2 s^-1 and re-weight the detected counts - scale = old_div(source_diff_spectrum, response_energy_bin.sim_differential_photon_fluxes) * 1e9 + this_model_image_conv = this_model_image_norm_conv * factor + + return this_model_image_conv - # Now return the map multiplied by the scale factor - return np.sum(scale * response_energy_bin.sim_signal_events_per_bin) * this_map