Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Computation of ``p_astro`` by source category and utilities 

2related to ``p_astro.json`` source classification files. 

3See Kapadia et al (2019), :doi:`10.1088/1361-6382/ab5f2d`, for details. 

4""" 

5import io 

6import json 

7 

8from celery.utils.log import get_task_logger 

9from ligo import p_astro_computation as pastrocomp 

10from matplotlib import pyplot as plt 

11import numpy as np 

12 

13from . import gracedb, lvalert 

14from .. import app 

15from ..util import PromiseProxy, read_json 

16 

17MEAN_VALUES_DICT = PromiseProxy( 

18 read_json, ('ligo.data', 'H1L1V1-mean_counts-1126051217-61603201.json')) 

19 

20THRESHOLDS_DICT = PromiseProxy( 

21 read_json, ('ligo.data', 'H1L1V1-pipeline-far_snr-thresholds.json')) 

22 

23P_ASTRO_LIVETIME = PromiseProxy( 

24 read_json, ('ligo.data', 'p_astro_livetime.json')) 

25 

26 

27log = get_task_logger(__name__) 

28 

29 

30@app.task(shared=False) 

31def compute_p_astro(snr, far, mass1, mass2, pipeline, instruments): 

32 """ 

33 Task to compute `p_astro` by source category. 

34 

35 Parameters 

36 ---------- 

37 snr : float 

38 event's SNR 

39 far : float 

40 event's cfar 

41 mass1 : float 

42 event's mass1 

43 mass2 : float 

44 event's mass2 

45 instruments : set 

46 set of instruments that detected the event 

47 

48 Returns 

49 ------- 

50 p_astros : str 

51 JSON dump of the p_astro by source category 

52 

53 Example 

54 ------- 

55 >>> p_astros = json.loads(compute_p_astro(files)) 

56 >>> p_astros 

57 {'BNS': 0.999, 'BBH': 0.0, 'NSBH': 0.0, 'Terrestrial': 0.001} 

58 

59 """ 

60 # Ensure SNR does not increase indefinitely beyond limiting FAR 

61 # for MBTA and PyCBC events 

62 snr_choice = pastrocomp.choose_snr(far, 

63 snr, 

64 pipeline, 

65 instruments, 

66 THRESHOLDS_DICT) 

67 

68 # Define constants to compute bayesfactors 

69 snr_star = 8.5 

70 far_star = 1 / (30 * 86400) 

71 

72 # Compute astrophysical bayesfactor for 

73 # GraceDB event 

74 fground = 3 * snr_star**3 / (snr_choice**4) 

75 bground = far / far_star 

76 astro_bayesfac = fground / bground 

77 

78 # Update terrestrial count based on far threshold 

79 lam_0 = far_star * P_ASTRO_LIVETIME['p_astro_livetime'] 

80 mean_values_dict = dict(MEAN_VALUES_DICT) 

81 mean_values_dict["counts_Terrestrial"] = lam_0 

82 

83 # Compute categorical p_astro values 

84 p_astro_values = \ 

85 pastrocomp.evaluate_p_astro_from_bayesfac(astro_bayesfac, 

86 mean_values_dict, 

87 mass1, 

88 mass2) 

89 # Dump mean values in json file 

90 return json.dumps(p_astro_values) 

91 

92 

93def _format_prob(prob): 

94 if prob >= 1: 

95 return '100%' 

96 elif prob <= 0: 

97 return '0%' 

98 elif prob > 0.99: 

99 return '>99%' 

100 elif prob < 0.01: 

101 return '<1%' 

102 else: 

103 return '{}%'.format(int(np.round(100 * prob))) 

104 

105 

106@app.task(shared=False) 

107def plot(contents): 

108 """Make a visualization of the source classification. 

109 

110 Parameters 

111 ---------- 

112 contents : str, bytes 

113 The contents of the ``p_astro.json`` file. 

114 

115 Returns 

116 ------- 

117 png : bytes 

118 The contents of a PNG file. 

119 

120 Notes 

121 ----- 

122 The unusually small size of the plot (2.5 x 2 inches) is optimized for 

123 viewing in GraceDB's image display widget. 

124 

125 Examples 

126 -------- 

127 .. plot:: 

128 :include-source: 

129 

130 >>> from gwcelery.tasks import p_astro 

131 >>> contents = ''' 

132 ... {"Terrestrial": 0.001, "BNS": 0.65, "NSBH": 0.20, 

133 ... "MassGap": 0.10, "BBH": 0.059} 

134 ... ''' 

135 >>> p_astro.plot(contents) 

136 

137 """ 

138 classification = json.loads(contents) 

139 outfile = io.BytesIO() 

140 

141 probs, names = zip( 

142 *sorted(zip(classification.values(), classification.keys()))) 

143 

144 with plt.style.context('seaborn-white'): 

145 fig, ax = plt.subplots(figsize=(2.5, 2)) 

146 ax.barh(names, probs) 

147 for i, prob in enumerate(probs): 

148 ax.annotate(_format_prob(prob), (0, i), (4, 0), 

149 textcoords='offset points', ha='left', va='center') 

150 ax.set_xlim(0, 1) 

151 ax.set_xticks([]) 

152 ax.tick_params(left=False) 

153 for side in ['top', 'bottom', 'right']: 

154 ax.spines[side].set_visible(False) 

155 fig.tight_layout() 

156 fig.savefig(outfile, format='png') 

157 return outfile.getvalue() 

158 

159 

160@lvalert.handler('superevent', 

161 'mdc_superevent', 

162 shared=False) 

163def handle(alert): 

164 """LVAlert handler to plot and upload a visualization of every 

165 ``p_astro.json`` that is added to a superevent. 

166 """ 

167 graceid = alert['uid'] 

168 filename = 'p_astro.json' 

169 

170 if alert['alert_type'] == 'log' and alert['data']['filename'] == filename: 

171 ( 

172 gracedb.download.s(filename, graceid) 

173 | 

174 plot.s() 

175 | 

176 gracedb.upload.s( 

177 filename.replace('.json', '.png'), graceid, 

178 message=( 

179 'Source classification visualization from ' 

180 '<a href="/api/superevents/{graceid}/files/{filename}">' 

181 '{filename}</a>').format( 

182 graceid=graceid, filename=filename), 

183 tags=['em_follow', 'p_astro', 'public'] 

184 ) 

185 ).delay()