kNNpy.HelperFunctions_2DA
1#################################################################################################### 2 3#------------------- These libraries are required for evaluating the functions ------------------- 4 5import numpy as np 6import scipy 7import healpy as hp 8import time 9import copy 10 11#################################################################################################### 12 13#-------------------------------------- Function Definitions ------------------------------------- 14 15#################################################################################################### 16 17def create_query_2DA(NSIDE_query, mask, tolerance, Verbose=False): 18 r''' 19 Computes the usable query points for the given mask (ie., query points at least a user-defined 20 threshold distance away from the mask edge) and returns the same, along with a HEALPix 21 `query mask` that has the following values: 22 23 0: pixels outside the observational footprint 24 1: pixels inside the footprint but too close to the mask edge (not usable) 25 2: usable pixels 26 27 Parameters 28 ---------- 29 NSIDE_query : int 30 the HEALPix NSIDE of the query grid (needs to be the same as that of the continuous field and the mask). Must be a power of 2 (eg. 128, 256, 512, etc.). 31 mask : numpy float array of shape ``(12*NSIDE_query**2, )`` 32 array encoding the observational footprint associated with the data. The value of the mask should be ``1.0`` for HEALPixels inside the observational footprint and ``healpy.UNSEEN`` for HEALPixels outside the observational footprint. ``healpy.UNSEEN = -1.6375e+30`` is a special value for masked pixels used by the ``healpy`` package. If there is no observational footprint (for example, data such as gravitational wave catalogs that are all-sky, or simulated datasets), please enter an array with all values equal to ``1.0``. 33 tolerance : float 34 the minimum angular distance (in radians) a query point needs to be away from the mask edge 35 to be considered usable. 36 Verbose : bool, optional 37 if set to ``True``, the time taken to complete each step of the calculation will be printed, by default ``False``. 38 39 Returns 40 ------- 41 query_mask : numpy int array of shape ``mask.shape`` 42 the HEALPix query mask, i.e., an array with 0, 1 and 2 indicating that the corresponding HEALPixel is outside the observational footprint, too close to mask boundary and sufficiently inside the observational footprint (far from the boundary), respectively. 43 44 QueryPositions : numpy float array of shape ``(N_usable_pix, 2)`` 45 array of usable query point positions, where 'N_usable_pix' is the number of pixels that are sufficiently far away from the mask edge, as determined by this method. For each query point in the array, the first (second) coordinate is the declination (right ascension) in radians. 46 47 Raises 48 ------ 49 ValueError 50 if `tolerance` is not in `[0, 2*np.pi]` 51 ValueError 52 if `NSIDE_query` is not a power of 2 53 ValueError 54 if `NSIDE_query` is not the same as the NSIDE of the continuous field and the mask 55 56 See Also 57 -------- 58 kNNpy.HelperFunctions.create_query_3D : generates query points in 3D. 59 60 Notes 61 ----- 62 Please refer to Gupta & Banerjee (2024)[^1] for a detailed discussion on creation of query point in presence of observational footprints that do not cover the full sky. The algorithm currently supports only query grids of the same size as the HEALpix grid on which the continuous overdensity field skymap is defined. 63 64 References 65 ---------- 66 [^1]: Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, [Monthly Notices of the Royal Astronomical Society](https://doi.org/10.1093/mnras/stae1424), Volume 531, Issue 4, July 2024, Pages 4619–4639 67 ''' 68 69 #----------------------------------------------------------------------------------------------- 70 71 #Check all inputs are consistent with the function requirement 72 73 if Verbose: 74 total_start = time.perf_counter() 75 print('\nValidating inputs...') 76 77 #Check if 'tolerance' is a valid angle 78 if tolerance<0 or tolerance>2*np.pi: 79 raise ValueError('Invalid threshold distance: please ensure 0 <= tolerance <= 2*pi.') 80 81 #Check if 'NSIDE_query' is a power of 2 82 if not (NSIDE_query > 0 and (NSIDE_query & (NSIDE_query - 1)) == 0): 83 raise ValueError('Invalid NSIDE for the query grid: please ensure NSIDE is a power of 2') 84 85 #Getting the NSIDE for the original mask 86 NSIDE_mask = hp.get_nside(mask) 87 88 #Check if the NSIDEs match 89 if NSIDE_mask!=NSIDE_query: 90 raise ValueError(f'NSIDE of the query grid ({NSIDE_query}) does not match NSIDE of \ 91 the mask ({NSIDE_mask}).') 92 93 if Verbose: print('\tdone.') 94 95 #----------------------------------------------------------------------------------------------- 96 97 #Query points are defined on a Healpix grid, points in or close to the mask are removed 98 99 if Verbose: print('\nDefining a query mask...') 100 101 #Getting number of pixels from the NSIDE for the query points 102 NPIX = hp.nside2npix(NSIDE_query) 103 104 #----------------------------------------------------------------------------------------------- 105 106 #Getting the query mask 107 108 #2 means that the query point is usable 109 query_mask = 2*np.ones(NPIX) 110 111 #Check if non-trivial observational mask is present 112 if np.any(mask == hp.UNSEEN): 113 114 if Verbose: 115 print('\tNon-trivial observational mask detected. Now removing the query points outside the observational footprint and the query points too close to the mask boundary...') 116 117 pixels = np.arange(NPIX) 118 qpos_arr = np.array(hp.pix2vec(NSIDE_query, pixels)) 119 #Looping over the initial query pixels 120 for pix_ind, qpos in enumerate(qpos_arr.T): 121 perc_progress = int((pix_ind+1)*100/len(qpos_arr.T)) 122 perc_progress_2 = int((pix_ind+2)*100/len(qpos_arr.T)) 123 if Verbose and perc_progress%10==0 and perc_progress!=perc_progress_2: print('\t{}% done'.format(perc_progress)) 124 sel_pix_close = hp.query_disc(NSIDE_query, qpos, tolerance) 125 #1 means that the query point is close to the mask edge 126 if np.any(mask[sel_pix_close]==hp.UNSEEN): query_mask[pix_ind] = 1 127 128 pix_inside_arr = hp.vec2pix(NSIDE_query, qpos_arr[0], qpos_arr[1], qpos_arr[2]) 129 #0 means that the query point is outside the footprint 130 query_mask[mask[pix_inside_arr]==hp.UNSEEN] = 0 131 132 #----------------------------------------------------------------------------------------------- 133 134 #Getting the query positions 135 136 if Verbose: print('\nGetting the query positions...') 137 138 query_pixels = np.arange(NPIX)[query_mask==2] 139 140 #Getting the latitudes and longitudes in degrees 141 QueryPositions_Deg = np.transpose(hp.pixelfunc.pix2ang(NSIDE_query, query_pixels, lonlat=True)) 142 QueryPositions_Deg[:, [0, 1]] = QueryPositions_Deg[:, [1, 0]] 143 144 #converting to radians 145 QueryPositions = np.deg2rad(QueryPositions_Deg) 146 147 if Verbose: 148 print('\tdone.') 149 print('\nTotal time taken: {:.2e} s'.format(time.perf_counter()-total_start)) 150 151 return query_mask, QueryPositions 152 153#################################################################################################### 154 155def bl_th(l, ss): 156 r''' 157 Computes Legendre expansion coefficients for the top-hat window function in angular coordinates. 158 159 160 Parameters 161 ---------- 162 l : numpy int array 163 array of multipole numbers. 164 ss : float 165 angular scale (in radians) at which the field is to be smoothed. 166 167 Returns 168 ------- 169 numpy float array of shape ``l.shape`` 170 array of Legendre expansion coefficients at each input multipole number. 171 ''' 172 173 plm1 = scipy.special.eval_legendre(l-1, np.cos(ss)) 174 plp1 = scipy.special.eval_legendre(l+1, np.cos(ss)) 175 num = (plm1-plp1) 176 den = 4*np.pi*(1-np.cos(ss)) 177 return num/den 178 179#################################################################################################### 180 181def top_hat_smoothing_2DA(skymap, scale, Verbose=False): 182 r''' 183 Smooths the given map at the given scale using the top hat window function in harmonic space. 184 185 Parameters 186 ---------- 187 skymap : numpy float array 188 the healpy map of the continuous field that needs to be smoothed. The values of the masked pixels, if any, should be set to `healpy.UNSEEN`. ``healpy.UNSEEN = -1.6375e+30`` is a special value for masked pixels used by the ``healpy`` package. 189 scale : float 190 angular scale (in radians) at which the field is to be smoothed. Please ensure `scale` is between `0` and `2*np.pi`. 191 Verbose : bool, optional 192 if set to `True`, the time taken to complete each step of the calculation will be printed, by default `False`. 193 194 Returns 195 ------- 196 smoothed_map_masked : numpy float array of shape ``skymap.shape`` 197 the smoothed healpy map, keeping the masked pixels of the original map masked. 198 199 Raises 200 ------ 201 ValueError 202 if `scale` is not in `[0, 2*np.pi]` 203 204 See Also 205 -------- 206 kNNpy.HelperFunctions.smoothing_3d : performs smoothing operations in 3D. 207 208 Notes 209 ----- 210 The following expression is used to compute the the spherical harmonic expansion coefficients $\alpha^{\theta}_{\ell m}$ of the field smoothed at angular scale $\theta$ using a top hat window function (See Devaraju (2015)[^1] and Gupta & Banerjee (2024)[^2] for derivations and a detailed discussion) 211 $$\alpha^{\theta}_{\ell m} = 4\pi\frac{b_{\ell}} {2\ell+1}\alpha_{\ell m},$$ 212 where $b_{\ell}$ are the the Legedre expansion coefficients of the top hat function, given by 213 $$b_{\ell} = \frac{1}{4\pi(1-\cos\theta)}\left[P_{\ell-1}(\cos\theta)-P_{\ell+1}(\cos\theta)\right].$$ 214 The smoothed field is reconstructed from $\alpha^{\theta}_{\ell m}$ using healpy's `alm2map` method. 215 216 217 References 218 ---------- 219 [^1]: Devaraju B., 2015, [doctoralThesis](http://elib.uni-stuttgart.de/handle/11682/4002), doi:10.18419/opus-3985. 220 [^2]: Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, [Monthly Notices of the Royal Astronomical Society](https://doi.org/10.1093/mnras/stae1424), Volume 531, Issue 4, July 2024, Pages 4619–4639 221 ''' 222 223 #----------------------------------------------------------------------------------------------- 224 225 if Verbose: total_start_time = time.perf_counter() 226 227 #----------------------------------------------------------------------------------------------- 228 229 #Check all inputs are consistent with the function requirement 230 231 if Verbose: print('Checking inputs ...') 232 233 if scale<0 or scale>2*np.pi: 234 raise ValueError('Invalid angular smoothing scale: please ensure 0 <= scale <= 2*pi.') 235 236 if Verbose: print('\tdone.') 237 238 #----------------------------------------------------------------------------------------------- 239 240 #Getting the NSIDE of the map 241 NSIDE = hp.get_nside(skymap) 242 243 #Using the default l_max for the harmonic expansion 244 l_max = 3*NSIDE-1 245 246 #----------------------------------------------------------------------------------------------- 247 248 #Computing the spherical harmonic expansion for the map 249 if Verbose: 250 print('\nComputing the spherical harmonic expansion for the map ...') 251 l_arr = np.array(range(l_max+1)) 252 map_alm = hp.map2alm(skymap, use_pixel_weights=False, lmax=l_max) 253 if Verbose: print('\tdone.') 254 255 #----------------------------------------------------------------------------------------------- 256 257 #Getting the legendre coefficients of the top hat window function 258 if Verbose: 259 print('\nGetting the legendre coefficients of the top hat window function ...') 260 bl = bl_th(l_arr, scale) 261 if Verbose: print('\tdone.') 262 263 #----------------------------------------------------------------------------------------------- 264 265 #Getting the spherical harmonic expansion of the smoothed map 266 if Verbose: 267 print('\nGetting the spherical harmonic expansion of the smoothed map ...') 268 smoothed_alm = np.zeros(map_alm.shape[0], dtype = 'complex') 269 for l in range(l_max+1): 270 for m in range(l+1): 271 ind_lm = hp.Alm.getidx(l_max, l, m) 272 smoothed_alm[ind_lm] = map_alm[ind_lm]*bl[l]*4*np.pi/(2*l+1) 273 if Verbose: print('\tdone.') 274 275 #----------------------------------------------------------------------------------------------- 276 277 #Generating the smoothed map from the harmonic coeffients 278 if Verbose: 279 print('\nGenerating the smoothed map from the harmonic coeffients ...') 280 smoothed_map = hp.alm2map(smoothed_alm, nside=NSIDE) 281 if Verbose: print('\tdone.') 282 283 #----------------------------------------------------------------------------------------------- 284 285 #If a mask is present 286 smoothed_map_masked = copy.deepcopy(smoothed_map) 287 smoothed_map_masked[skymap==hp.UNSEEN] = hp.UNSEEN 288 289 #----------------------------------------------------------------------------------------------- 290 291 if Verbose: print('\ntotal time taken: {:.2e} s.'.format(time.perf_counter()-total_start_time)) 292 293 return smoothed_map_masked 294 295#################################################################################################### 296 297def create_smoothed_field_dict_2DA(skymap, bins, query_mask, Verbose=False): 298 r''' 299 Creates a dictionary containing the continuous field smoothed at various angular distance scales. 300 301 Parameters 302 ---------- 303 skymap : numpy float array 304 the healpy map of the continuous field that needs to be smoothed. The values of the masked pixels, if any, should be set to `healpy.UNSEEN`. ``healpy.UNSEEN = -1.6375e+30`` is a special value for masked pixels used by the ``healpy`` package. 305 bins : list of numpy float array 306 list of distances for each nearest neighbour. The $i^{th}$ element of the list should contain a numpy array of the desired distance scales for the $i^{th}$ nearest neighbour. 307 query_mask : numpy float array of shape ``skymap.shape`` 308 the HEALPix query mask. 309 Verbose : bool, optional 310 if set to `True`, the time taken to complete each step of the calculation will be printed, by default `False`. 311 312 Returns 313 ------- 314 SmoothedFieldDict : dict 315 dictionary containing the continuous field masked within the observational footprint and smoothed at various angular distance scales. For example, `SmoothedFieldDict['0.215']` represents the continuous map smoothed at a scale of 0.215 radians. 316 317 Notes 318 ----- 319 `query_mask` is a numpy int array with 0, 1 and 2 indicating that the corresponding HEALPixel is outside the mask, too close to mask boundary and sufficiently far away from the boundary, respectively. Please Refer to the helper function method `create_query_2DA()` for creating the query mask. See also Gupta and Banerjee (2024)[^1] for a discussion. 320 321 References 322 ---------- 323 [^1]: Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, [Monthly Notices of the Royal Astronomical Society](https://doi.org/10.1093/mnras/stae1424), Volume 531, Issue 4, July 2024, Pages 4619–4639 324 ''' 325 326 #----------------------------------------------------------------------------------------------- 327 328 if Verbose: 329 total_start = time.perf_counter() 330 print(f'\nSmoothing the density field over the given angular distance scales...\n') 331 332 #----------------------------------------------------------------------------------------------- 333 334 #Initializing the dictionary 335 SmoothedFieldDict = {} 336 337 #----------------------------------------------------------------------------------------------- 338 339 #Looping over the nearest neighbour indices as inferred from the length of 'bins' 340 for i in range(len(bins)): 341 342 #------------------------------------------------------------------------------------------- 343 344 if Verbose: start = time.perf_counter() 345 346 #------------------------------------------------------------------------------------------- 347 348 for j, ss in enumerate(bins[i]): 349 #Square bracket selects only those pixels that are not close to the mask boundaries 350 #Turning verbose off to avoid too much text output 351 SmoothedFieldDict[str(ss)] = \ 352 top_hat_smoothing_2DA(skymap, ss, Verbose=False)[query_mask==2] 353 354 #------------------------------------------------------------------------------------------- 355 356 if Verbose: 357 print('\tdistance scales for {}NN done; time taken: {:.2e} s.'.format(i+1, time.perf_counter()-start)) 358 359 #----------------------------------------------------------------------------------------------- 360 361 if Verbose: print('\ntotal time taken: {:.2e} s.'.format(time.perf_counter()-total_start)) 362 363 return SmoothedFieldDict 364 365#################################################################################################### 366 367#---------------------------------------- END OF PROGRAM! ---------------------------------------- 368 369####################################################################################################
def
create_query_2DA(NSIDE_query, mask, tolerance, Verbose=False):
18def create_query_2DA(NSIDE_query, mask, tolerance, Verbose=False): 19 r''' 20 Computes the usable query points for the given mask (ie., query points at least a user-defined 21 threshold distance away from the mask edge) and returns the same, along with a HEALPix 22 `query mask` that has the following values: 23 24 0: pixels outside the observational footprint 25 1: pixels inside the footprint but too close to the mask edge (not usable) 26 2: usable pixels 27 28 Parameters 29 ---------- 30 NSIDE_query : int 31 the HEALPix NSIDE of the query grid (needs to be the same as that of the continuous field and the mask). Must be a power of 2 (eg. 128, 256, 512, etc.). 32 mask : numpy float array of shape ``(12*NSIDE_query**2, )`` 33 array encoding the observational footprint associated with the data. The value of the mask should be ``1.0`` for HEALPixels inside the observational footprint and ``healpy.UNSEEN`` for HEALPixels outside the observational footprint. ``healpy.UNSEEN = -1.6375e+30`` is a special value for masked pixels used by the ``healpy`` package. If there is no observational footprint (for example, data such as gravitational wave catalogs that are all-sky, or simulated datasets), please enter an array with all values equal to ``1.0``. 34 tolerance : float 35 the minimum angular distance (in radians) a query point needs to be away from the mask edge 36 to be considered usable. 37 Verbose : bool, optional 38 if set to ``True``, the time taken to complete each step of the calculation will be printed, by default ``False``. 39 40 Returns 41 ------- 42 query_mask : numpy int array of shape ``mask.shape`` 43 the HEALPix query mask, i.e., an array with 0, 1 and 2 indicating that the corresponding HEALPixel is outside the observational footprint, too close to mask boundary and sufficiently inside the observational footprint (far from the boundary), respectively. 44 45 QueryPositions : numpy float array of shape ``(N_usable_pix, 2)`` 46 array of usable query point positions, where 'N_usable_pix' is the number of pixels that are sufficiently far away from the mask edge, as determined by this method. For each query point in the array, the first (second) coordinate is the declination (right ascension) in radians. 47 48 Raises 49 ------ 50 ValueError 51 if `tolerance` is not in `[0, 2*np.pi]` 52 ValueError 53 if `NSIDE_query` is not a power of 2 54 ValueError 55 if `NSIDE_query` is not the same as the NSIDE of the continuous field and the mask 56 57 See Also 58 -------- 59 kNNpy.HelperFunctions.create_query_3D : generates query points in 3D. 60 61 Notes 62 ----- 63 Please refer to Gupta & Banerjee (2024)[^1] for a detailed discussion on creation of query point in presence of observational footprints that do not cover the full sky. The algorithm currently supports only query grids of the same size as the HEALpix grid on which the continuous overdensity field skymap is defined. 64 65 References 66 ---------- 67 [^1]: Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, [Monthly Notices of the Royal Astronomical Society](https://doi.org/10.1093/mnras/stae1424), Volume 531, Issue 4, July 2024, Pages 4619–4639 68 ''' 69 70 #----------------------------------------------------------------------------------------------- 71 72 #Check all inputs are consistent with the function requirement 73 74 if Verbose: 75 total_start = time.perf_counter() 76 print('\nValidating inputs...') 77 78 #Check if 'tolerance' is a valid angle 79 if tolerance<0 or tolerance>2*np.pi: 80 raise ValueError('Invalid threshold distance: please ensure 0 <= tolerance <= 2*pi.') 81 82 #Check if 'NSIDE_query' is a power of 2 83 if not (NSIDE_query > 0 and (NSIDE_query & (NSIDE_query - 1)) == 0): 84 raise ValueError('Invalid NSIDE for the query grid: please ensure NSIDE is a power of 2') 85 86 #Getting the NSIDE for the original mask 87 NSIDE_mask = hp.get_nside(mask) 88 89 #Check if the NSIDEs match 90 if NSIDE_mask!=NSIDE_query: 91 raise ValueError(f'NSIDE of the query grid ({NSIDE_query}) does not match NSIDE of \ 92 the mask ({NSIDE_mask}).') 93 94 if Verbose: print('\tdone.') 95 96 #----------------------------------------------------------------------------------------------- 97 98 #Query points are defined on a Healpix grid, points in or close to the mask are removed 99 100 if Verbose: print('\nDefining a query mask...') 101 102 #Getting number of pixels from the NSIDE for the query points 103 NPIX = hp.nside2npix(NSIDE_query) 104 105 #----------------------------------------------------------------------------------------------- 106 107 #Getting the query mask 108 109 #2 means that the query point is usable 110 query_mask = 2*np.ones(NPIX) 111 112 #Check if non-trivial observational mask is present 113 if np.any(mask == hp.UNSEEN): 114 115 if Verbose: 116 print('\tNon-trivial observational mask detected. Now removing the query points outside the observational footprint and the query points too close to the mask boundary...') 117 118 pixels = np.arange(NPIX) 119 qpos_arr = np.array(hp.pix2vec(NSIDE_query, pixels)) 120 #Looping over the initial query pixels 121 for pix_ind, qpos in enumerate(qpos_arr.T): 122 perc_progress = int((pix_ind+1)*100/len(qpos_arr.T)) 123 perc_progress_2 = int((pix_ind+2)*100/len(qpos_arr.T)) 124 if Verbose and perc_progress%10==0 and perc_progress!=perc_progress_2: print('\t{}% done'.format(perc_progress)) 125 sel_pix_close = hp.query_disc(NSIDE_query, qpos, tolerance) 126 #1 means that the query point is close to the mask edge 127 if np.any(mask[sel_pix_close]==hp.UNSEEN): query_mask[pix_ind] = 1 128 129 pix_inside_arr = hp.vec2pix(NSIDE_query, qpos_arr[0], qpos_arr[1], qpos_arr[2]) 130 #0 means that the query point is outside the footprint 131 query_mask[mask[pix_inside_arr]==hp.UNSEEN] = 0 132 133 #----------------------------------------------------------------------------------------------- 134 135 #Getting the query positions 136 137 if Verbose: print('\nGetting the query positions...') 138 139 query_pixels = np.arange(NPIX)[query_mask==2] 140 141 #Getting the latitudes and longitudes in degrees 142 QueryPositions_Deg = np.transpose(hp.pixelfunc.pix2ang(NSIDE_query, query_pixels, lonlat=True)) 143 QueryPositions_Deg[:, [0, 1]] = QueryPositions_Deg[:, [1, 0]] 144 145 #converting to radians 146 QueryPositions = np.deg2rad(QueryPositions_Deg) 147 148 if Verbose: 149 print('\tdone.') 150 print('\nTotal time taken: {:.2e} s'.format(time.perf_counter()-total_start)) 151 152 return query_mask, QueryPositions
Computes the usable query points for the given mask (ie., query points at least a user-defined
threshold distance away from the mask edge) and returns the same, along with a HEALPix
query mask that has the following values:
0: pixels outside the observational footprint
1: pixels inside the footprint but too close to the mask edge (not usable)
2: usable pixels
Parameters
- NSIDE_query (int): the HEALPix NSIDE of the query grid (needs to be the same as that of the continuous field and the mask). Must be a power of 2 (eg. 128, 256, 512, etc.).
- mask (numpy float array of shape
(12*NSIDE_query**2, )): array encoding the observational footprint associated with the data. The value of the mask should be1.0for HEALPixels inside the observational footprint andhealpy.UNSEENfor HEALPixels outside the observational footprint.healpy.UNSEEN = -1.6375e+30is a special value for masked pixels used by thehealpypackage. If there is no observational footprint (for example, data such as gravitational wave catalogs that are all-sky, or simulated datasets), please enter an array with all values equal to1.0. - tolerance (float): the minimum angular distance (in radians) a query point needs to be away from the mask edge to be considered usable.
- Verbose (bool, optional):
if set to
True, the time taken to complete each step of the calculation will be printed, by defaultFalse.
Returns
- query_mask (numpy int array of shape
mask.shape): the HEALPix query mask, i.e., an array with 0, 1 and 2 indicating that the corresponding HEALPixel is outside the observational footprint, too close to mask boundary and sufficiently inside the observational footprint (far from the boundary), respectively. - QueryPositions (numpy float array of shape
(N_usable_pix, 2)): array of usable query point positions, where 'N_usable_pix' is the number of pixels that are sufficiently far away from the mask edge, as determined by this method. For each query point in the array, the first (second) coordinate is the declination (right ascension) in radians.
Raises
- ValueError: if
toleranceis not in[0, 2*np.pi] - ValueError: if
NSIDE_queryis not a power of 2 - ValueError: if
NSIDE_queryis not the same as the NSIDE of the continuous field and the mask
See Also
kNNpy.HelperFunctions.create_query_3D: generates query points in 3D.
Notes
Please refer to Gupta & Banerjee (2024)1 for a detailed discussion on creation of query point in presence of observational footprints that do not cover the full sky. The algorithm currently supports only query grids of the same size as the HEALpix grid on which the continuous overdensity field skymap is defined.
References
-
Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, Monthly Notices of the Royal Astronomical Society, Volume 531, Issue 4, July 2024, Pages 4619–4639 ↩
def
bl_th(l, ss):
156def bl_th(l, ss): 157 r''' 158 Computes Legendre expansion coefficients for the top-hat window function in angular coordinates. 159 160 161 Parameters 162 ---------- 163 l : numpy int array 164 array of multipole numbers. 165 ss : float 166 angular scale (in radians) at which the field is to be smoothed. 167 168 Returns 169 ------- 170 numpy float array of shape ``l.shape`` 171 array of Legendre expansion coefficients at each input multipole number. 172 ''' 173 174 plm1 = scipy.special.eval_legendre(l-1, np.cos(ss)) 175 plp1 = scipy.special.eval_legendre(l+1, np.cos(ss)) 176 num = (plm1-plp1) 177 den = 4*np.pi*(1-np.cos(ss)) 178 return num/den
Computes Legendre expansion coefficients for the top-hat window function in angular coordinates.
Parameters
- l (numpy int array): array of multipole numbers.
- ss (float): angular scale (in radians) at which the field is to be smoothed.
Returns
- numpy float array of shape
l.shape: array of Legendre expansion coefficients at each input multipole number.
def
top_hat_smoothing_2DA(skymap, scale, Verbose=False):
182def top_hat_smoothing_2DA(skymap, scale, Verbose=False): 183 r''' 184 Smooths the given map at the given scale using the top hat window function in harmonic space. 185 186 Parameters 187 ---------- 188 skymap : numpy float array 189 the healpy map of the continuous field that needs to be smoothed. The values of the masked pixels, if any, should be set to `healpy.UNSEEN`. ``healpy.UNSEEN = -1.6375e+30`` is a special value for masked pixels used by the ``healpy`` package. 190 scale : float 191 angular scale (in radians) at which the field is to be smoothed. Please ensure `scale` is between `0` and `2*np.pi`. 192 Verbose : bool, optional 193 if set to `True`, the time taken to complete each step of the calculation will be printed, by default `False`. 194 195 Returns 196 ------- 197 smoothed_map_masked : numpy float array of shape ``skymap.shape`` 198 the smoothed healpy map, keeping the masked pixels of the original map masked. 199 200 Raises 201 ------ 202 ValueError 203 if `scale` is not in `[0, 2*np.pi]` 204 205 See Also 206 -------- 207 kNNpy.HelperFunctions.smoothing_3d : performs smoothing operations in 3D. 208 209 Notes 210 ----- 211 The following expression is used to compute the the spherical harmonic expansion coefficients $\alpha^{\theta}_{\ell m}$ of the field smoothed at angular scale $\theta$ using a top hat window function (See Devaraju (2015)[^1] and Gupta & Banerjee (2024)[^2] for derivations and a detailed discussion) 212 $$\alpha^{\theta}_{\ell m} = 4\pi\frac{b_{\ell}} {2\ell+1}\alpha_{\ell m},$$ 213 where $b_{\ell}$ are the the Legedre expansion coefficients of the top hat function, given by 214 $$b_{\ell} = \frac{1}{4\pi(1-\cos\theta)}\left[P_{\ell-1}(\cos\theta)-P_{\ell+1}(\cos\theta)\right].$$ 215 The smoothed field is reconstructed from $\alpha^{\theta}_{\ell m}$ using healpy's `alm2map` method. 216 217 218 References 219 ---------- 220 [^1]: Devaraju B., 2015, [doctoralThesis](http://elib.uni-stuttgart.de/handle/11682/4002), doi:10.18419/opus-3985. 221 [^2]: Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, [Monthly Notices of the Royal Astronomical Society](https://doi.org/10.1093/mnras/stae1424), Volume 531, Issue 4, July 2024, Pages 4619–4639 222 ''' 223 224 #----------------------------------------------------------------------------------------------- 225 226 if Verbose: total_start_time = time.perf_counter() 227 228 #----------------------------------------------------------------------------------------------- 229 230 #Check all inputs are consistent with the function requirement 231 232 if Verbose: print('Checking inputs ...') 233 234 if scale<0 or scale>2*np.pi: 235 raise ValueError('Invalid angular smoothing scale: please ensure 0 <= scale <= 2*pi.') 236 237 if Verbose: print('\tdone.') 238 239 #----------------------------------------------------------------------------------------------- 240 241 #Getting the NSIDE of the map 242 NSIDE = hp.get_nside(skymap) 243 244 #Using the default l_max for the harmonic expansion 245 l_max = 3*NSIDE-1 246 247 #----------------------------------------------------------------------------------------------- 248 249 #Computing the spherical harmonic expansion for the map 250 if Verbose: 251 print('\nComputing the spherical harmonic expansion for the map ...') 252 l_arr = np.array(range(l_max+1)) 253 map_alm = hp.map2alm(skymap, use_pixel_weights=False, lmax=l_max) 254 if Verbose: print('\tdone.') 255 256 #----------------------------------------------------------------------------------------------- 257 258 #Getting the legendre coefficients of the top hat window function 259 if Verbose: 260 print('\nGetting the legendre coefficients of the top hat window function ...') 261 bl = bl_th(l_arr, scale) 262 if Verbose: print('\tdone.') 263 264 #----------------------------------------------------------------------------------------------- 265 266 #Getting the spherical harmonic expansion of the smoothed map 267 if Verbose: 268 print('\nGetting the spherical harmonic expansion of the smoothed map ...') 269 smoothed_alm = np.zeros(map_alm.shape[0], dtype = 'complex') 270 for l in range(l_max+1): 271 for m in range(l+1): 272 ind_lm = hp.Alm.getidx(l_max, l, m) 273 smoothed_alm[ind_lm] = map_alm[ind_lm]*bl[l]*4*np.pi/(2*l+1) 274 if Verbose: print('\tdone.') 275 276 #----------------------------------------------------------------------------------------------- 277 278 #Generating the smoothed map from the harmonic coeffients 279 if Verbose: 280 print('\nGenerating the smoothed map from the harmonic coeffients ...') 281 smoothed_map = hp.alm2map(smoothed_alm, nside=NSIDE) 282 if Verbose: print('\tdone.') 283 284 #----------------------------------------------------------------------------------------------- 285 286 #If a mask is present 287 smoothed_map_masked = copy.deepcopy(smoothed_map) 288 smoothed_map_masked[skymap==hp.UNSEEN] = hp.UNSEEN 289 290 #----------------------------------------------------------------------------------------------- 291 292 if Verbose: print('\ntotal time taken: {:.2e} s.'.format(time.perf_counter()-total_start_time)) 293 294 return smoothed_map_masked
Smooths the given map at the given scale using the top hat window function in harmonic space.
Parameters
- skymap (numpy float array):
the healpy map of the continuous field that needs to be smoothed. The values of the masked pixels, if any, should be set to
healpy.UNSEEN.healpy.UNSEEN = -1.6375e+30is a special value for masked pixels used by thehealpypackage. - scale (float):
angular scale (in radians) at which the field is to be smoothed. Please ensure
scaleis between0and2*np.pi. - Verbose (bool, optional):
if set to
True, the time taken to complete each step of the calculation will be printed, by defaultFalse.
Returns
- smoothed_map_masked (numpy float array of shape
skymap.shape): the smoothed healpy map, keeping the masked pixels of the original map masked.
Raises
- ValueError: if
scaleis not in[0, 2*np.pi]
See Also
kNNpy.HelperFunctions.smoothing_3d: performs smoothing operations in 3D.
Notes
The following expression is used to compute the the spherical harmonic expansion coefficients $\alpha^{\theta}_{\ell m}$ of the field smoothed at angular scale $\theta$ using a top hat window function (See Devaraju (2015)1 and Gupta & Banerjee (2024)2 for derivations and a detailed discussion)
$$\alpha^{\theta}_{\ell m} = 4\pi\frac{b_{\ell}} {2\ell+1}\alpha_{\ell m},$$
where $b_{\ell}$ are the the Legedre expansion coefficients of the top hat function, given by
$$b_{\ell} = \frac{1}{4\pi(1-\cos\theta)}\left[P_{\ell-1}(\cos\theta)-P_{\ell+1}(\cos\theta)\right].$$
The smoothed field is reconstructed from $\alpha^{\theta}_{\ell m}$ using healpy's alm2map method.
References
-
Devaraju B., 2015, doctoralThesis, doi:10.18419/opus-3985. ↩
-
Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, Monthly Notices of the Royal Astronomical Society, Volume 531, Issue 4, July 2024, Pages 4619–4639 ↩
def
create_smoothed_field_dict_2DA(skymap, bins, query_mask, Verbose=False):
298def create_smoothed_field_dict_2DA(skymap, bins, query_mask, Verbose=False): 299 r''' 300 Creates a dictionary containing the continuous field smoothed at various angular distance scales. 301 302 Parameters 303 ---------- 304 skymap : numpy float array 305 the healpy map of the continuous field that needs to be smoothed. The values of the masked pixels, if any, should be set to `healpy.UNSEEN`. ``healpy.UNSEEN = -1.6375e+30`` is a special value for masked pixels used by the ``healpy`` package. 306 bins : list of numpy float array 307 list of distances for each nearest neighbour. The $i^{th}$ element of the list should contain a numpy array of the desired distance scales for the $i^{th}$ nearest neighbour. 308 query_mask : numpy float array of shape ``skymap.shape`` 309 the HEALPix query mask. 310 Verbose : bool, optional 311 if set to `True`, the time taken to complete each step of the calculation will be printed, by default `False`. 312 313 Returns 314 ------- 315 SmoothedFieldDict : dict 316 dictionary containing the continuous field masked within the observational footprint and smoothed at various angular distance scales. For example, `SmoothedFieldDict['0.215']` represents the continuous map smoothed at a scale of 0.215 radians. 317 318 Notes 319 ----- 320 `query_mask` is a numpy int array with 0, 1 and 2 indicating that the corresponding HEALPixel is outside the mask, too close to mask boundary and sufficiently far away from the boundary, respectively. Please Refer to the helper function method `create_query_2DA()` for creating the query mask. See also Gupta and Banerjee (2024)[^1] for a discussion. 321 322 References 323 ---------- 324 [^1]: Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, [Monthly Notices of the Royal Astronomical Society](https://doi.org/10.1093/mnras/stae1424), Volume 531, Issue 4, July 2024, Pages 4619–4639 325 ''' 326 327 #----------------------------------------------------------------------------------------------- 328 329 if Verbose: 330 total_start = time.perf_counter() 331 print(f'\nSmoothing the density field over the given angular distance scales...\n') 332 333 #----------------------------------------------------------------------------------------------- 334 335 #Initializing the dictionary 336 SmoothedFieldDict = {} 337 338 #----------------------------------------------------------------------------------------------- 339 340 #Looping over the nearest neighbour indices as inferred from the length of 'bins' 341 for i in range(len(bins)): 342 343 #------------------------------------------------------------------------------------------- 344 345 if Verbose: start = time.perf_counter() 346 347 #------------------------------------------------------------------------------------------- 348 349 for j, ss in enumerate(bins[i]): 350 #Square bracket selects only those pixels that are not close to the mask boundaries 351 #Turning verbose off to avoid too much text output 352 SmoothedFieldDict[str(ss)] = \ 353 top_hat_smoothing_2DA(skymap, ss, Verbose=False)[query_mask==2] 354 355 #------------------------------------------------------------------------------------------- 356 357 if Verbose: 358 print('\tdistance scales for {}NN done; time taken: {:.2e} s.'.format(i+1, time.perf_counter()-start)) 359 360 #----------------------------------------------------------------------------------------------- 361 362 if Verbose: print('\ntotal time taken: {:.2e} s.'.format(time.perf_counter()-total_start)) 363 364 return SmoothedFieldDict
Creates a dictionary containing the continuous field smoothed at various angular distance scales.
Parameters
- skymap (numpy float array):
the healpy map of the continuous field that needs to be smoothed. The values of the masked pixels, if any, should be set to
healpy.UNSEEN.healpy.UNSEEN = -1.6375e+30is a special value for masked pixels used by thehealpypackage. - bins (list of numpy float array): list of distances for each nearest neighbour. The $i^{th}$ element of the list should contain a numpy array of the desired distance scales for the $i^{th}$ nearest neighbour.
- query_mask (numpy float array of shape
skymap.shape): the HEALPix query mask. - Verbose (bool, optional):
if set to
True, the time taken to complete each step of the calculation will be printed, by defaultFalse.
Returns
- SmoothedFieldDict (dict):
dictionary containing the continuous field masked within the observational footprint and smoothed at various angular distance scales. For example,
SmoothedFieldDict['0.215']represents the continuous map smoothed at a scale of 0.215 radians.
Notes
query_mask is a numpy int array with 0, 1 and 2 indicating that the corresponding HEALPixel is outside the mask, too close to mask boundary and sufficiently far away from the boundary, respectively. Please Refer to the helper function method create_query_2DA() for creating the query mask. See also Gupta and Banerjee (2024)1 for a discussion.
References
-
Kaustubh Rajesh Gupta, Arka Banerjee, Spatial clustering of gravitational wave sources with k-nearest neighbour distributions, Monthly Notices of the Royal Astronomical Society, Volume 531, Issue 4, July 2024, Pages 4619–4639 ↩