Coverage for ligo / followup_advocate / __init__.py: 91%

1from contextlib import contextmanager 

2import math as mth 

3import os 

4import shutil 

5import tempfile 

6import urllib 

7import webbrowser 

8 

9import astropy.coordinates as coord 

10import astropy.time 

11import astropy.units as u 

12import healpy as hp 

13import lxml.etree 

14import numpy as np 

15from astropy.io.fits import getheader 

16from ligo.gracedb import rest 

17from ligo.skymap.io.fits import read_sky_map 

18from ligo.skymap.postprocess.ellipse import find_ellipse 

19from ligo.skymap.postprocess.crossmatch import crossmatch 

20 

21from .jinja import env 

22import importlib.metadata 

23__version__ = importlib.metadata.version(__name__) 

24 

25 

26def authors(authors, service=rest.DEFAULT_SERVICE_URL): 

27 """Write GCN Circular author list""" 

28 return env.get_template('authors.jinja2').render(authors=authors).strip() 

29 

30 

31def guess_skyloc_pipeline(filename): 

32 skyloc_pipelines_dict = { 

33 'cwb': 'cWB_AllSky', 

34 'bayestar': 'BAYESTAR', 

35 'bilby': 'Bilby', 

36 'lib': 'LIB', 

37 'lalinference': 'LALInference', 

38 'olib': 'oLIB_AllSky', 

39 'mly': 'MLy_AllSky', 

40 'rapidpe_rift': 'RapidPE-RIFT', 

41 'amplfi': 'AMPLFI' 

42 } 

43 try: 

44 return skyloc_pipelines_dict[filename.split('.')[0].lower()] 

45 except KeyError: 

46 return filename.split('.')[0] 

47 

48 

49def text_width(remove_text_wrap): 

50 """Return width of text wrap based on whether we wish to wrap the lines or 

51 not.""" 

52 return 9999 if remove_text_wrap else 79 

53 

54 

55def main_dict(gracedb_id, client, raven_coinc=False, update_alert=False, 

56 cgmi_filename=None): 

57 """Create general dictionary to pass to compose circular""" 

58 

59 superevent = client.superevent(gracedb_id).json() 

60 preferred_event = superevent['preferred_event_data'] 

61 pipeline = preferred_event['pipeline'] 

62 search = preferred_event['search'] 

63 preferred_pipeline_search = f'{pipeline}_{search}' 

64 early_warning_pipeline_searches = [] 

65 pipeline_searches = [] 

66 gw_events = superevent['gw_events'] 

67 early_warning_alert = False 

68 

69 for gw_event in gw_events: 

70 gw_event_dict = client.event(gw_event).json() 

71 pipeline = gw_event_dict['pipeline'] 

72 search = gw_event_dict['search'] 

73 # Remap MDC to allsky to not double count 

74 if search == 'MDC': 

75 search = 'AllSky' 

76 pipeline_search = f'{pipeline}_{search}' 

77 # Create separates lists of post merger and early warning pipelines 

78 if pipeline_search not in pipeline_searches and \ 

79 'EarlyWarning' not in pipeline_search: 

80 pipeline_searches.append(pipeline_search) 

81 if 'EarlyWarning' in pipeline_search: 

82 # Remap MBTA early warning to the all sky citation prior to 

83 # checking if there is a previous entry 

84 if 'mbta' in pipeline_search.lower(): 

85 pipeline_search = 'MBTA_AllSky' 

86 if pipeline_search not in early_warning_pipeline_searches: 

87 early_warning_pipeline_searches.append(pipeline_search) 

88 

89 if not pipeline_searches: 

90 raise ValueError( 

91 "{} has no post-merger events to generate circulars from.".format( 

92 gracedb_id)) 

93 

94 # Sort to get alphabetical order 

95 pipeline_searches.sort(key=str.lower) 

96 early_warning_pipeline_searches.sort(key=str.lower) 

97 

98 voevents = client.voevents(gracedb_id).json()['voevents'] 

99 if not voevents: 

100 raise ValueError( 

101 "{} has no VOEvent to generate circulars from.".format( 

102 gracedb_id)) 

103 

104 citation_index = { 

105 pipeline_search.lower(): pipeline_searches.index(pipeline_search) + 1 

106 for pipeline_search in pipeline_searches} 

107 for pipeline_search in early_warning_pipeline_searches: 

108 if pipeline_search.lower() not in citation_index: 

109 citation_index[pipeline_search.lower()] = \ 

110 max(citation_index.values()) + 1 

111 

112 gpstime = float(preferred_event['gpstime']) 

113 event_time = astropy.time.Time(gpstime, format='gps').utc 

114 

115 # Grab latest p_astro and em_bright values from lastest VOEvent 

116 voevent_text_latest = \ 

117 client.files(gracedb_id, voevents[-1]['filename']).read() 

118 root = lxml.etree.fromstring(voevent_text_latest) 

119 p_astros = root.find('./What/Group[@name="Classification"]') 

120 em_brights = root.find('./What/Group[@name="Properties"]') 

121 classifications = {} 

122 source_classification = {} 

123 mchirp_bin = {} 

124 # Only try to load if present to prevent errors with .getchildren() 

125 p_astro_pipeline = None 

126 search_for_p_astro = False 

127 search_for_cgmi_file = False 

128 logs = [] 

129 

130 # Grab em_bright values if present 

131 for em_bright in em_brights.getchildren(): 

132 if em_bright.attrib.get('value'): 

133 source_classification[em_bright.attrib['name']] = \ 

134 float(em_bright.attrib['value']) * 100 

135 

136 # Catch empty filenames or any other values that evaluate to false 

137 if not cgmi_filename: 

138 cgmi_filename = None 

139 

140 # Try to look for CGMI file if there are actual em bright values for a 

141 # preliminary/initial alert, and there is no manual filename given 

142 search_for_cgmi_file = len(source_classification) > 0 and not \ 

143 update_alert and cgmi_filename is None 

144 

145 # Grab p_astro values if present 

146 for p_astro in p_astros.getchildren(): 

147 if p_astro.attrib.get('value'): 

148 classifications[p_astro.attrib['name']] = \ 

149 float(p_astro.attrib['value']) * 100 

150 search_for_p_astro = True 

151 

152 # Grab logs for either RapidPE checks or chirp mass bins 

153 if search_for_p_astro or search_for_cgmi_file: 

154 logs = client.logs(gracedb_id).json()['log'] 

155 

156 # Figure out which pipeline uploaded p_astro, usually the first one 

157 # FIXME: Replace with more efficient method in the future 

158 for message in reversed(logs): 

159 filename = message['filename'] 

160 if filename and '.p_astro.json' in filename and \ 

161 filename != 'p_astro.json': 

162 p_astro = client.files(gracedb_id, filename).json() 

163 if all(mth.isclose(p_astro[key] * 100, classifications[key]) 

164 for key in classifications.keys()): 

165 p_astro_pipeline = filename.split('.')[0].lower() 

166 break 

167 

168 # Adjust citations if needed 

169 if p_astro_pipeline == 'rapidpe_rift': 

170 citation_index['rapidpe_rift'] = max(citation_index.values()) + 1 

171 if len(source_classification) > 0: 

172 citation_index['em_bright'] = max(citation_index.values()) + 1 

173 

174 # Look for CGMI file if a preliminary/inital alert with embright info or 

175 # if given a filename, usually for an update alert 

176 if search_for_cgmi_file or cgmi_filename is not None: 

177 cgmi_json = {} 

178 # chirp mass estimates included when em-bright is 

179 if cgmi_filename is not None: 

180 cgmi_json = client.files(gracedb_id, cgmi_filename).json() 

181 else: 

182 for message in reversed(logs): 

183 filename = message['filename'] 

184 # use most recent mchirp estimate 

185 if filename and 'mchirp_source' in filename and \ 

186 '.json' in filename: 

187 cgmi_json = client.files(gracedb_id, filename).json() 

188 break 

189 # if CGMI file found either way, include in results 

190 if cgmi_json: 

191 mchirp_bin_edges = cgmi_json['bin_edges'] 

192 mchirp_probs = cgmi_json['probabilities'] 

193 # find the highest probability bin 

194 max_prob_idx = np.argmax(mchirp_probs) 

195 left_bin = mchirp_bin_edges[max_prob_idx] 

196 right_bin = mchirp_bin_edges[max_prob_idx+1] 

197 mchirp_bin = left_bin, right_bin 

198 

199 skymaps = {} 

200 voevents_text = [] 

201 for voevent in voevents: 

202 # Don't load latest voevent since already loaded from before 

203 if voevent == voevents[-1]: 

204 voevent_text = voevent_text_latest 

205 # If earlier voevent, load 

206 else: 

207 voevent_text = client.files(gracedb_id, voevent['filename']).read() 

208 root = lxml.etree.fromstring(voevent_text) 

209 alert_type = root.find( 

210 './What/Param[@name="AlertType"]').attrib['value'].lower() 

211 if alert_type == 'earlywarning': 

212 # Add text for early warning detection if one early warning alert 

213 early_warning_alert = True 

214 url = root.find('./What/Group/Param[@name="skymap_fits"]') 

215 if url is None: 

216 continue 

217 url = url.attrib['value'] 

218 _, filename = os.path.split(url) 

219 skyloc_pipeline = guess_skyloc_pipeline(filename) 

220 issued_time = astropy.time.Time(root.find('./Who/Date').text).gps 

221 if filename not in skymaps: 

222 skymaps[filename] = dict( 

223 alert_type=alert_type, 

224 pipeline=skyloc_pipeline, 

225 filename=filename, 

226 latency=issued_time-event_time.gps) 

227 if skyloc_pipeline.lower() not in citation_index: 

228 citation_index[skyloc_pipeline.lower()] = \ 

229 max(citation_index.values()) + 1 

230 voevents_text.append(voevent_text) 

231 skymaps = list(skymaps.values()) 

232 

233 o = urllib.parse.urlparse(client.service_url) 

234 

235 kwargs = dict( 

236 subject='Identification', 

237 gracedb_id=gracedb_id, 

238 gracedb_service_url=urllib.parse.urlunsplit( 

239 (o.scheme, o.netloc, '/superevents/', '', '')), 

240 update_alert=update_alert, 

241 cgmi_filename=cgmi_filename, 

242 group=preferred_event['group'], 

243 pipeline_search=preferred_pipeline_search, 

244 post_merger_pipeline_searches=pipeline_searches, 

245 early_warning_alert=early_warning_alert, 

246 early_warning_pipeline_searches=early_warning_pipeline_searches, 

247 gpstime='{0:.03f}'.format(round(float(preferred_event['gpstime']), 3)), 

248 search=preferred_event.get('search', ''), 

249 far=preferred_event['far'], 

250 utctime=event_time.iso, 

251 instruments=preferred_event['instruments'].split(','), 

252 skymaps=skymaps, 

253 prob_has_ns=source_classification.get('HasNS'), 

254 prob_has_remnant=source_classification.get('HasRemnant'), 

255 prob_has_massgap=source_classification.get('HasMassGap'), 

256 prob_has_ssm=source_classification.get('HasSSM'), 

257 source_classification=source_classification, 

258 mchirp_bin=mchirp_bin, 

259 include_ellipse=None, 

260 classifications=classifications, 

261 p_astro_pipeline=p_astro_pipeline, 

262 citation_index=citation_index) 

263 

264 if skymaps: 

265 preferred_skymap = skymaps[-1]['filename'] 

266 cls = [50, 90] 

267 include_ellipse, ra, dec, a, b, pa, area, greedy_area = \ 

268 uncertainty_ellipse(gracedb_id, preferred_skymap, client, cls=cls) 

269 kwargs.update( 

270 preferred_skymap=preferred_skymap, 

271 cl=cls[-1], 

272 include_ellipse=include_ellipse, 

273 ra=coord.Longitude(ra*u.deg), 

274 dec=coord.Latitude(dec*u.deg), 

275 a=coord.Angle(a*u.deg), 

276 b=coord.Angle(b*u.deg), 

277 pa=coord.Angle(pa*u.deg), 

278 ellipse_area=area, 

279 greedy_area=greedy_area) 

280 try: 

281 distmu, distsig = get_distances_skymap_gracedb(gracedb_id, 

282 preferred_skymap, 

283 client) 

284 kwargs.update( 

285 distmu=distmu, 

286 distsig=distsig) 

287 except TypeError: 

288 pass 

289 

290 if raven_coinc: 

291 kwargs = _update_raven_parameters(superevent, kwargs, client, 

292 voevents_text) 

293 return kwargs 

294 

295 

296def compose(gracedb_id, authors=(), mailto=False, remove_text_wrap=False, 

297 service=rest.DEFAULT_SERVICE_URL, client=None): 

298 """Compose GCN Circular draft""" 

299 

300 if client is None: 

301 client = rest.GraceDb(service) 

302 

303 kwargs = main_dict(gracedb_id, client=client) 

304 kwargs.update(authors=authors) 

305 kwargs.update(gw_is_subthreshold=False) 

306 kwargs.update(text_width=text_width(remove_text_wrap)) 

307 

308 subject = env.get_template('subject.jinja2').render(**kwargs).strip() 

309 body = env.get_template('initial_circular.jinja2').render(**kwargs).strip() 

310 

311 if mailto: 

312 pattern = 'mailto:emfollow@ligo.org,dac@ligo.org?subject={0}&body={1}' 

313 url = pattern.format( 

314 urllib.parse.quote(subject), 

315 urllib.parse.quote(body)) 

316 webbrowser.open(url) 

317 else: 

318 return '{0}\n\n{1}'.format(subject, body) 

319 

320 

321def compose_raven(gracedb_id, authors=(), remove_text_wrap=False, 

322 service=rest.DEFAULT_SERVICE_URL, client=None, 

323 gw_is_subthreshold=False): 

324 """Compose EM_COINC RAVEN GCN Circular draft""" 

325 

326 if client is None: 

327 client = rest.GraceDb(service) 

328 

329 kwargs = dict() 

330 kwargs.update(main_dict(gracedb_id, client=client, raven_coinc=True)) 

331 kwargs.update(text_width=text_width(remove_text_wrap)) 

332 kwargs.update(gw_is_subthreshold=gw_is_subthreshold) 

333 # Add RAVEN citation 

334 citation_index = kwargs['citation_index'] 

335 citation_index['raven'] = max(citation_index.values()) + 1 

336 kwargs['citation_index'] = citation_index 

337 

338 subject = ( 

339 env.get_template('RAVEN_subject.jinja2').render(**kwargs) 

340 .strip()) 

341 body = ( 

342 env.get_template('RAVEN_circular.jinja2').render(**kwargs) 

343 .strip()) 

344 return '{0}\n\n{1}'.format(subject, body) 

345 

346 

347def compose_llama( 

348 gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL, 

349 icecube_alert=None, remove_text_wrap=False, 

350 client=None): 

351 """Compose GRB LLAMA GCN Circular draft. 

352 Here, gracedb_id will be a GRB superevent ID in GraceDb.""" 

353 

354 if client is None: 

355 client = rest.GraceDb(service) 

356 

357 superevent = client.superevent(gracedb_id).json() 

358 

359 gpstime = float(superevent['t_0']) 

360 tl, th = gpstime - 500, gpstime + 500 

361 event_time = astropy.time.Time(gpstime, format='gps').utc 

362 tl_datetime = str(astropy.time.Time( 

363 tl, format='gps').isot).replace('T', ' ') 

364 th_datetime = str(astropy.time.Time( 

365 th, format='gps').isot).replace('T', ' ') 

366 

367 o = urllib.parse.urlparse(client.service_url) 

368 kwargs = dict( 

369 gracedb_service_url=urllib.parse.urlunsplit( 

370 (o.scheme, o.netloc, '/superevents/', '', '')), 

371 gracedb_id=gracedb_id, 

372 llama=True, 

373 icecube_alert=icecube_alert, 

374 event_time=event_time, 

375 tl_datetime=tl_datetime, 

376 th_datetime=th_datetime, 

377 authors=authors) 

378 kwargs.update(text_width=text_width(remove_text_wrap)) 

379 

380 citation_index = {'llama': 1, 'llama_stat': 2} 

381 kwargs.update(citation_index=citation_index) 

382 

383 files = client.files(gracedb_id).json() 

384 

385 llama_stat_filename = 'significance_subthreshold_lvc-i3.json' 

386 if llama_stat_filename in files: 

387 llama_stat_file = client.files(gracedb_id, llama_stat_filename).json() 

388 llama_fap = llama_stat_file["p_value"] 

389 neutrinos = llama_stat_file["inputs"]["neutrino_info"] 

390 lines = [] 

391 for neutrino in neutrinos: 

392 # Build aligned string 

393 vals = [] 

394 dt = (event_time - 

395 astropy.time.Time(neutrino['mjd'], 

396 format='mjd')).to(u.s).value 

397 vals.append('{:.2f}'.format(dt)) 

398 vals.append('{:.2f}'.format(neutrino['ra'])) 

399 vals.append('{:.2f}'.format(neutrino['dec'])) 

400 vals.append('{:.2f}'.format(neutrino['sigma'])) 

401 vals.append('{:.4f}'.format(llama_fap)) 

402 lines.append(align_number_string(vals, [0, 11, 21, 40, 59])) 

403 

404 kwargs.update(llama_fap=llama_fap, 

405 neutrinos=lines) 

406 

407 subject = ( 

408 env.get_template('llama_subject.jinja2').render(**kwargs) 

409 .strip()) 

410 if icecube_alert: 

411 body = ( 

412 env.get_template('llama_alert_circular.jinja2').render(**kwargs) 

413 .strip()) 

414 else: 

415 body = ( 

416 env.get_template('llama_track_circular.jinja2').render(**kwargs) 

417 .strip()) 

418 return '{0}\n\n{1}'.format(subject, body) 

419 

420 

421def compose_grb_medium_latency( 

422 gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL, 

423 use_detection_template=None, remove_text_wrap=False, client=None): 

424 """Compose GRB Medium Latency GCN Circular draft. 

425 Here, gracedb_id will be a GRB external trigger ID in GraceDb.""" 

426 

427 if client is None: 

428 client = rest.GraceDb(service) 

429 

430 event = client.event(gracedb_id).json() 

431 search = event['search'] 

432 

433 if search != 'GRB': 

434 return 

435 

436 group = event['group'] 

437 pipeline = event['pipeline'] 

438 external_trigger = event['extra_attributes']['GRB']['trigger_id'] 

439 

440 o = urllib.parse.urlparse(client.service_url) 

441 kwargs = dict( 

442 gracedb_service_url=urllib.parse.urlunsplit( 

443 (o.scheme, o.netloc, '/events/', '', '')), 

444 gracedb_id=gracedb_id, 

445 group=group, 

446 grb=True, 

447 pipeline=pipeline, 

448 external_trigger=external_trigger, 

449 exclusions=[], 

450 detections=[]) 

451 kwargs.update(text_width=text_width(remove_text_wrap)) 

452 

453 files = client.files(gracedb_id).json() 

454 

455 citation_index = {} 

456 index = 1 

457 xpipeline_fap_file = 'false_alarm_probability_x.json' 

458 if xpipeline_fap_file in files: 

459 xpipeline_fap = client.files(gracedb_id, xpipeline_fap_file).json() 

460 online_xpipeline_fap = xpipeline_fap.get('Online Xpipeline') 

461 # Create detection/exclusion circular based on given argument 

462 # Use default cutoff if not provided 

463 xpipeline_detection = (use_detection_template if use_detection_template 

464 is not None else online_xpipeline_fap < 0.001) 

465 if xpipeline_detection: 

466 kwargs['detections'].append('xpipeline') 

467 kwargs.update(online_xpipeline_fap=online_xpipeline_fap) 

468 else: 

469 kwargs['exclusions'].append('xpipeline') 

470 xpipeline_distances_file = 'distances_x.json' 

471 xpipeline_distances = client.files(gracedb_id, 

472 xpipeline_distances_file).json() 

473 burst_exclusion = xpipeline_distances.get('Burst Exclusion') 

474 kwargs.update(burst_exclusion=burst_exclusion) 

475 citation_index['xpipeline'] = index 

476 index += 1 

477 

478 pygrb_fap_file = 'false_alarm_probability_pygrb.json' 

479 if pygrb_fap_file in files: 

480 pygrb_fap = client.files(gracedb_id, pygrb_fap_file).json() 

481 online_pygrb_fap = pygrb_fap.get('Online PyGRB') 

482 # Create detection/exclusion circular based on given argument 

483 # Use default cutoff if not provided 

484 pygrb_detection = (use_detection_template if use_detection_template 

485 is not None else online_pygrb_fap < 0.01) 

486 if pygrb_detection: 

487 kwargs['detections'].append('pygrb') 

488 kwargs.update(online_pygrb_fap=online_pygrb_fap) 

489 else: 

490 kwargs['exclusions'].append('pygrb') 

491 pygrb_distances_file = 'distances_pygrb.json' 

492 pygrb_distances = client.files(gracedb_id, 

493 pygrb_distances_file).json() 

494 nsns_exclusion = pygrb_distances.get('NSNS Exclusion') 

495 nsbh_exclusion = pygrb_distances.get('NSBH Exclusion') 

496 kwargs.update(nsbh_exclusion=nsbh_exclusion, 

497 nsns_exclusion=nsns_exclusion) 

498 citation_index['pygrb'] = index 

499 

500 kwargs.update(citation_index=citation_index) 

501 

502 subject = ( 

503 env.get_template('medium_latency_grb_subject.jinja2').render(**kwargs) 

504 .strip()) 

505 body = ( 

506 env.get_template('medium_latency_grb_circular.jinja2').render(**kwargs) 

507 .strip()) 

508 return '{0}\n\n{1}'.format(subject, body) 

509 

510 

511def compose_update(gracedb_id, authors=(), 

512 service=rest.DEFAULT_SERVICE_URL, 

513 update_types=['sky_localization', 'p_astro', 

514 'em_bright', 'raven'], 

515 remove_text_wrap=False, 

516 client=None, 

517 cgmi_filename=None): 

518 """Compose GCN update circular""" 

519 if client is None: 

520 client = rest.GraceDb(service) 

521 

522 kwargs = main_dict(gracedb_id, client=client, 

523 raven_coinc='raven' in update_types, 

524 update_alert=True, 

525 cgmi_filename=cgmi_filename) 

526 kwargs.pop('citation_index', None) 

527 kwargs.update(text_width=text_width(remove_text_wrap)) 

528 

529 if isinstance(update_types, str): 

530 update_types = update_types.split(',') 

531 

532 # Adjust files for update type alert 

533 citation_index = {} 

534 skymaps = [] 

535 index = 1 

536 

537 updated_skymap = kwargs.get('skymaps')[-1] 

538 kwargs.update(updated_skymap=updated_skymap) 

539 skymaps = [updated_skymap] 

540 if cgmi_filename: 

541 update_types.append('cgmi') 

542 if 'sky_localization' in update_types: 

543 citation_index[updated_skymap['pipeline'].lower()] = index 

544 index += 1 

545 if 'p_astro' in update_types and \ 

546 kwargs.get('p_astro_pipeline') == 'rapidpe_rift': 

547 citation_index['rapidpe_rift'] = index 

548 index += 1 

549 if 'em_bright' in update_types: 

550 # If not already cited, cite sky map pipeline for em_bright 

551 if updated_skymap['pipeline'].lower() not in citation_index.keys(): 

552 citation_index[updated_skymap['pipeline'].lower()] = index 

553 index += 1 

554 citation_index['em_bright'] = index 

555 index += 1 

556 if 'raven' in update_types: 

557 citation_index['raven'] = index 

558 

559 kwargs.update(skymaps=skymaps, 

560 citation_index=citation_index, 

561 post_merger_pipeline_searches=[], 

562 update_alert=True) 

563 

564 kwargs.update(authors=authors) 

565 kwargs.update(gw_is_subthreshold=False) 

566 kwargs.update(subject='Update') 

567 kwargs.update(update_types=update_types) 

568 

569 subject = env.get_template('subject.jinja2').render(**kwargs).strip() 

570 body = env.get_template( 

571 'update_circular.jinja2').render(**kwargs).strip() 

572 return '{0}\n\n{1}'.format(subject, body) 

573 

574 

575def compose_retraction(gracedb_id, authors=(), remove_text_wrap=False, 

576 service=rest.DEFAULT_SERVICE_URL, client=None): 

577 """Compose GCN retraction circular""" 

578 if client is None: 

579 client = rest.GraceDb(service) 

580 event = client.superevent(gracedb_id).json() 

581 preferred_event = event['preferred_event_data'] 

582 labels = event['labels'] 

583 earlywarning = \ 

584 ('EARLY_WARNING' in labels and 

585 {'EM_SelectedConfident', 'SIGNIF_LOCKED'}.isdisjoint(labels)) 

586 

587 kwargs = dict( 

588 subject='Retraction', 

589 gracedb_id=gracedb_id, 

590 group=preferred_event['group'], 

591 earlywarning=earlywarning, 

592 authors=authors 

593 ) 

594 kwargs.update(text_width=text_width(remove_text_wrap)) 

595 

596 subject = env.get_template('subject.jinja2').render(**kwargs).strip() 

597 body = env.get_template('retraction.jinja2').render(**kwargs).strip() 

598 return '{0}\n\n{1}'.format(subject, body) 

599 

600 

601@contextmanager 

602def open_map_gracedb(graceid, filename, client): 

603 with tempfile.NamedTemporaryFile(mode='w+b') as localfile: 

604 if filename.endswith('.gz'): 

605 # Try using the multi-res sky map if it exists 

606 new_filename = filename.replace('.fits.gz', '.multiorder.fits') 

607 try: 

608 remotefile = client.files(graceid, new_filename, raw=True) 

609 except (IOError, rest.HTTPError): 

610 remotefile = client.files(graceid, filename, raw=True) 

611 else: 

612 remotefile = client.files(graceid, filename, raw=True) 

613 try: 

614 shutil.copyfileobj(remotefile, localfile) 

615 finally: 

616 remotefile.close() 

617 localfile.flush() 

618 localfile.seek(0) 

619 yield localfile 

620 

621 

622def get_distances_skymap_gracedb(graceid, filename, client): 

623 with open_map_gracedb(graceid, filename, client) as f: 

624 header = getheader(f.name, 1) 

625 try: 

626 return header['distmean'], header['diststd'] 

627 except KeyError: 

628 pass 

629 

630 

631def read_map_from_path(path, client): 

632 with open_map_gracedb(*path.split('/'), client) as f: 

633 return read_sky_map(f.name)[0] 

634 

635 

636def align_number_string(nums, positions): 

637 positions.append(len(nums[-1])) 

638 gen = (val + ' ' * (positions[i+1]-positions[i]-len(val)) 

639 for i, val in enumerate(nums)) 

640 return ''.join(gen) 

641 

642 

643def mask_cl(p, level=90): 

644 pflat = p.ravel() 

645 i = np.flipud(np.argsort(p)) 

646 cs = np.cumsum(pflat[i]) 

647 cls = np.empty_like(pflat) 

648 cls[i] = cs 

649 cls = cls.reshape(p.shape) 

650 return cls <= 1e-2 * level 

651 

652 

653def compare_skymaps(paths, service=rest.DEFAULT_SERVICE_URL, client=None): 

654 """Produce table of sky map overlaps""" 

655 if client is None: 

656 client = rest.GraceDb(service) 

657 filenames = [path.split('/')[1] for path in paths] 

658 pipelines = [guess_skyloc_pipeline(filename) for filename in filenames] 

659 probs = [read_map_from_path(path, client) for path in paths] 

660 npix = max(len(prob) for prob in probs) 

661 nside = hp.npix2nside(npix) 

662 deg2perpix = hp.nside2pixarea(nside, degrees=True) 

663 probs = [hp.ud_grade(prob, nside, power=-2) for prob in probs] 

664 masks = [mask_cl(prob) for prob in probs] 

665 areas = [mask.sum() * deg2perpix for mask in masks] 

666 joint_areas = [(mask & masks[-1]).sum() * deg2perpix for mask in masks] 

667 

668 kwargs = dict(params=zip(filenames, pipelines, areas, joint_areas)) 

669 

670 return env.get_template('compare_skymaps.jinja2').render(**kwargs) 

671 

672 

673def uncertainty_ellipse(graceid, filename, client, cls=[50, 90], 

674 ratio_ellipse_cl_areas=1.35): 

675 """Compute uncertainty ellipses for a given sky map 

676 

677 Parameters 

678 ---------- 

679 graceid: str 

680 ID of the trigger used by GraceDB 

681 filename: str 

682 File name of sky map 

683 client: class 

684 REST API client for HTTP connection 

685 cls: array-like 

686 List of percentage of minimal credible area used to check whether the 

687 areas are close to an ellipse, returning the values of the final item 

688 ratio_ellipse_cl_areas: float 

689 Ratio between ellipse area and minimal credible area from cl 

690 """ 

691 with open_map_gracedb(graceid, filename, client) as f: 

692 skymap = read_sky_map(f.name, moc=True) 

693 

694 # Discard any distance ionformation, to prevent `crossmatch()` from doing 

695 # CPU- and memory-intensive credible volume calculations that we don't use. 

696 skymap = skymap['UNIQ', 'PROBDENSITY'] 

697 

698 # Convert to an array if necessary 

699 if np.isscalar(cls): 

700 cls = [cls] 

701 cls = np.asarray(cls) 

702 

703 # Pass array of contour inteverals to get areas 

704 result = crossmatch(skymap, contours=cls / 100) 

705 greedy_areas = np.asarray(result.contour_areas) 

706 ra, dec, a, b, pa, ellipse_areas = find_ellipse(skymap, cl=cls) 

707 a, b = np.asarray(a), np.asarray(b) 

708 

709 # Only use ellipse if every confidence interval passes 

710 use_ellipse = \ 

711 np.all(ellipse_areas <= ratio_ellipse_cl_areas * greedy_areas) 

712 return (use_ellipse, ra, dec, a[-1], b[-1], pa, ellipse_areas[-1], 

713 greedy_areas[-1]) 

714 

715 

716def _update_raven_parameters(superevent, kwargs, client, voevents_text): 

717 """Update kwargs with parameters for RAVEN coincidence""" 

718 

719 gracedb_id = superevent['superevent_id'] 

720 

721 if 'EM_COINC' not in superevent['labels']: 

722 raise ValueError("No EM_COINC label for {}".format( 

723 gracedb_id)) 

724 

725 preferred_event = superevent['preferred_event_data'] 

726 group = preferred_event['group'] 

727 gpstime = float(preferred_event['gpstime']) 

728 event_time = astropy.time.Time(gpstime, format='gps').utc 

729 em_event_id = superevent['em_type'] 

730 

731 # FIXME: Grab more info from the latest VOEvent if deemed suitable 

732 em_event = client.event(em_event_id).json() 

733 external_pipeline = em_event['pipeline'] 

734 # Get all other pipelines 

735 ext_events = [client.event(id).json() for id