1# (c) Stefan Countryman, 2017 

2 

3""" 

4Make a MAT file that can be read in by LoadCWbData.m, in the style of the 

5FormatToCWbData.m function written by Rainer. Adapted from code and 

6information provided by Rainer Corley and Imre Bartos. 

7""" 

8 

9import logging 

10from llama.filehandler import FileHandler 

11from llama.files.lvc_skymap import LvcSkymapFits 

12from llama.files.skymap_info import SkymapInfo 

13 

14LOGGER = logging.getLogger(__name__) 

15 

16 

17def argsort_end_of_cr(pixel_values, total_probability=0.9): 

18 """Take input values, sort them in descending order, and perform a binary 

19 search to find how many values must be included such that the sum of those 

20 top values gives the desired total_probability. Return a list of indices 

21 which, when used to get a slice from the original list of pixel values, 

22 will give the cut off list. This can be used to sort multiple array fields 

23 (while cutting unneeded values) as if they were rows in a spreadsheet.""" 

24 sorted_rows = pixel_values.argsort()[::-1] 

25 sorted_pixel_values = pixel_values[sorted_rows] 

26 mini = 0 

27 maxi = len(sorted_pixel_values) - 1 

28 while mini + 1 != maxi: 

29 midpoint = (mini + maxi) // 2 

30 if sorted_pixel_values[:midpoint].sum() > total_probability: 

31 maxi = midpoint 

32 else: 

33 mini = midpoint 

34 return sorted_rows[:maxi+1] 

35 

36 

37@FileHandler.set_class_attributes 

38class LvcSkymapMat(FileHandler): 

39 """ 

40 A list of sky positions and probability densities in text form generated 

41 from the original Healpix-encoded LVC skymap. 

42 """ 

43 

44 FILENAME = "lvc_skymap.mat" 

45 DEPENDENCIES = (LvcSkymapFits, SkymapInfo) 

46 

47 def _generate(self): # pylint: disable=arguments-differ 

48 # get pixel coordinates and pixel values; we will only keep the pixels 

49 # that together form the 90% confidence region. 

50 import numpy as np 

51 import healpy as hp 

52 pixel_values = hp.read_map(LvcSkymapFits(self).fullpath, 

53 verbose=False) 

54 npix = len(pixel_values) 

55 # used to calculate area of pixels 

56 stepsize = np.sqrt(4 * np.pi / npix) * 180 / np.pi 

57 pixel_indices = range(npix) 

58 nside = hp.pixelfunc.npix2nside(npix) 

59 # coordinates are in radians 

60 theta_rad, phi_rad = hp.pixelfunc.pix2ang(nside, pixel_indices) 

61 # sort by probability in descending order and cut off all pixels 

62 # falling bellow the 90% confidence threshold 

63 # also, reshape them into column vectors 

64 reduced_rows = argsort_end_of_cr(pixel_values) 

65 reduced_pixel_values = np.array([pixel_values[reduced_rows]]) 

66 # .transpose() 

67 reduced_theta_rad = np.array([theta_rad[reduced_rows]]) # .transpose() 

68 reduced_phi_rad = np.array([phi_rad[reduced_rows]]) # .transpose() 

69 # cut 

70 # start making the skymap variable 

71 skymap = { 

72 'RA': reduced_phi_rad * 180 / np.pi, 

73 'DEC': 90 - reduced_theta_rad * 180 / np.pi, 

74 'probability': reduced_pixel_values, 

75 # every pixel has the same average step distance, so just use that. 

76 'step': reduced_phi_rad * 0 + stepsize, 

77 'step2': reduced_phi_rad * 0 + stepsize, 

78 # in the future, the below 2 lines might need to be generalized. 

79 'theta': 0, 

80 'phi': 0, 

81 'cumulative': np.cumsum(reduced_pixel_values) 

82 } 

83 # start making the CWbData variable. 

84 cwbdata = { 

85 # get the time and false alarm rate of the event from the initial 

86 # notice 

87 'time': SkymapInfo(self).event_time_gps, 

88 'rho': 1. / SkymapInfo(self).far, 

89 'map_length': len(reduced_rows), 

90 'skymap': skymap, 

91 'lag': [0.0, 0.0, 0.0], 

92 'inj_hrss': 0, 

93 'phi_RA': 0, 

94 'theta_DEC': 0, 

95 'frequency': 0, 

96 'filename': np.array('filename', dtype=object)} 

97 # write files to a .mat file, old format (version 6) 

98 # pylint doesn't know what it's talking about; savemat is in scipy.io 

99 # pylint: disable=no-member 

100 import scipy.io 

101 scipy.io.savemat(self.fullpath, {'CWbData': cwbdata}, 

102 do_compression=True)