Skip to content

new PS convolution #99

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
+358 −169
Open

new PS convolution #99

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
f7d4314
new PS convolution
astromatt42 Apr 3, 2024
66cc02c
increased oder of scipy shit from 1 to 2
astromatt42 Apr 8, 2024
cd3d3db
Final touches 1
astromatt42 Apr 25, 2024
6df31e9
Merge branch 'master' into PS_fix
astromatt42 Jun 7, 2024
919c0e0
fixing the issue with slow reading of root files by specifying the
torresramiro350 Mar 22, 2024
22b91e8
removing the version freeze for uproot, the speed issues should be fixed
torresramiro350 Mar 22, 2024
fa6c38c
updating some information in the README file
torresramiro350 Mar 22, 2024
f4d69fe
removing a duplicate of healpy and ensuring that pytest is installed
torresramiro350 May 23, 2024
625f2cb
Merge branch 'PS_fix' of github.com:astromatt42/hawc_hal into PS_fix
astromatt42 Jun 7, 2024
e7164f7
changed scipy shift to order 2 and some code cleanup
astromatt42 Jun 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
154 changes: 65 additions & 89 deletions hawc_hal/HAL.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down
Loading