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

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

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

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

See Kapadia et al (2019), arXiv:1903.06881, for details. 

""" 

import io 

import json 

 

from matplotlib import pyplot as plt 

import numpy as np 

import pkg_resources 

 

from celery.utils.log import get_task_logger 

from ligo import p_astro_computation as pastrocomp 

 

from . import gracedb, lvalert 

from .. import app 

 

log = get_task_logger(__name__) 

 

 

@app.task(shared=False) 

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

""" 

Task to compute `p_astro` by source category. 

 

Parameters 

---------- 

snr : float 

event's SNR 

far : float 

event's cfar 

mass1 : float 

event's mass1 

mass2 : float 

event's mass2 

instruments : set 

set of instruments that detected the event 

 

Returns 

------- 

p_astros : str 

JSON dump of the p_astro by source category 

 

Example 

------- 

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

>>> p_astros 

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

""" 

 

# Read mean values from file 

with open(pkg_resources.resource_filename( 

'ligo.data', 'H1L1V1-mean_counts-1126051217-61603201.json')) as f: 

mean_values_dict = json.load(f) 

 

# Read thresholds on FAR and SNR from file 

with open(pkg_resources.resource_filename( 

'ligo.data', 'H1L1V1-pipeline-far_snr-thresholds.json')) as f: 

thresholds_dict = json.load(f) 

 

# Read livetime corresponding to mean values from file 

with open(pkg_resources.resource_filename( 

'ligo.data', 'p_astro_livetime.json')) as f: 

p_astro_livetime = json.load(f)['p_astro_livetime'] 

 

# Ensure SNR does not increase indefinitely beyond limiting FAR 

# for MBTA and PyCBC events 

snr_choice = pastrocomp.choose_snr(far, 

snr, 

pipeline, 

instruments, 

thresholds_dict) 

 

# Define constants to compute bayesfactors 

snr_star = 8.5 

far_star = 1 / (30 * 86400) 

 

# Compute astrophysical bayesfactor for 

# GraceDB event 

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

bground = far / far_star 

astro_bayesfac = fground / bground 

 

# Update terrestrial count based on far threshold 

lam_0 = far_star * p_astro_livetime 

mean_values_dict["counts_Terrestrial"] = lam_0 

 

# Compute categorical p_astro values 

p_astro_values = \ 

pastrocomp.evaluate_p_astro_from_bayesfac(astro_bayesfac, 

mean_values_dict, 

mass1, 

mass2) 

# Dump mean values in json file 

return json.dumps(p_astro_values) 

 

 

def _format_prob(prob): 

if prob >= 1: 

return '100%' 

elif prob <= 0: 

return '0%' 

elif prob > 0.99: 

return '>99%' 

elif prob < 0.01: 

return '<1%' 

else: 

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

 

 

@app.task(shared=False) 

def plot(contents): 

"""Make a visualization of the source classification. 

 

Parameters 

---------- 

contents : str, bytes 

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

 

Returns 

------- 

png : bytes 

The contents of a PNG file. 

 

Notes 

----- 

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

viewing in GraceDB's image display widget. 

 

Examples 

-------- 

 

.. plot:: 

:include-source: 

 

>>> from gwcelery.tasks import p_astro 

>>> contents = ''' 

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

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

... ''' 

>>> p_astro.plot(contents) 

""" 

classification = json.loads(contents) 

outfile = io.BytesIO() 

 

probs, names = zip( 

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

 

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

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

ax.barh(names, probs) 

for i, prob in enumerate(probs): 

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

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

ax.set_xlim(0, 1) 

ax.set_xticks([]) 

ax.tick_params(left=False) 

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

ax.spines[side].set_visible(False) 

fig.tight_layout() 

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

return outfile.getvalue() 

 

 

@lvalert.handler('superevent', 

'mdc_superevent', 

shared=False) 

def handle(alert): 

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

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

graceid = alert['uid'] 

filename = 'p_astro.json' 

 

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

( 

gracedb.download.s(filename, graceid) 

| 

plot.s() 

| 

gracedb.upload.s( 

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

message=( 

'Source classification visualization from ' 

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

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

graceid=graceid, filename=filename), 

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

) 

).delay()