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

1import math as mth 

2import os 

3import shutil 

4import tempfile 

5import urllib 

6import webbrowser 

7 

8import astropy.coordinates as coord 

9import astropy.time 

10import astropy.units as u 

11import healpy as hp 

12import lxml.etree 

13import numpy as np 

14from astropy.io.fits import getheader 

15from ligo.gracedb import rest 

16from ligo.skymap.io.fits import read_sky_map 

17from ligo.skymap.postprocess.ellipse import find_ellipse 

18from ligo.skymap.postprocess.crossmatch import crossmatch 

19 

20from .jinja import env 

21import importlib.metadata 

22__version__ = importlib.metadata.version(__name__) 

23 

24 

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

26 """Write GCN Circular author list""" 

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

28 

29 

30def guess_skyloc_pipeline(filename): 

31 skyloc_pipelines_dict = { 

32 'cwb': 'cWB_AllSky', 

33 'bayestar': 'BAYESTAR', 

34 'bilby': 'Bilby', 

35 'lib': 'LIB', 

36 'lalinference': 'LALInference', 

37 'olib': 'oLIB_AllSky', 

38 'mly': 'MLy_AllSky', 

39 'rapidpe_rift': 'RapidPE-RIFT', 

40 'amplfi': 'AMPLFI' 

41 } 

42 try: 

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

44 except KeyError: 

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

46 

47 

48def text_width(remove_text_wrap): 

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

50 not.""" 

51 return 9999 if remove_text_wrap else 79 

52 

53 

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

55 cgmi_filename=None): 

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

57 

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

59 preferred_event = superevent['preferred_event_data'] 

60 pipeline = preferred_event['pipeline'] 

61 search = preferred_event['search'] 

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

63 early_warning_pipeline_searches = [] 

64 pipeline_searches = [] 

65 gw_events = superevent['gw_events'] 

66 early_warning_alert = False 

67 

68 for gw_event in gw_events: 

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

70 pipeline = gw_event_dict['pipeline'] 

71 search = gw_event_dict['search'] 

72 # Remap MDC to allsky to not double count 

73 if search == 'MDC': 

74 search = 'AllSky' 

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

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

77 if pipeline_search not in pipeline_searches and \ 

78 'EarlyWarning' not in pipeline_search: 

79 pipeline_searches.append(pipeline_search) 

80 if 'EarlyWarning' in pipeline_search: 

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

82 # checking if there is a previous entry 

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

84 pipeline_search = 'MBTA_AllSky' 

85 if pipeline_search not in early_warning_pipeline_searches: 

86 early_warning_pipeline_searches.append(pipeline_search) 

87 

88 if not pipeline_searches: 

89 raise ValueError( 

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

91 gracedb_id)) 

92 

93 # Sort to get alphabetical order 

94 pipeline_searches.sort(key=str.lower) 

95 early_warning_pipeline_searches.sort(key=str.lower) 

96 

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

98 if not voevents: 

99 raise ValueError( 

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

101 gracedb_id)) 

102 

103 citation_index = { 

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

105 for pipeline_search in pipeline_searches} 

106 for pipeline_search in early_warning_pipeline_searches: 

107 if pipeline_search.lower() not in citation_index: 

108 citation_index[pipeline_search.lower()] = \ 

109 max(citation_index.values()) + 1 

110 

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

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

113 

114 # Grab latest p_astro and em_bright values from lastest VOEvent 

115 voevent_text_latest = \ 

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

117 root = lxml.etree.fromstring(voevent_text_latest) 

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

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

120 classifications = {} 

121 source_classification = {} 

122 mchirp_bin = {} 

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

124 p_astro_pipeline = None 

125 search_for_p_astro = False 

126 search_for_cgmi_file = False 

127 logs = [] 

128 

129 # Grab em_bright values if present 

130 for em_bright in em_brights.getchildren(): 

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

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

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

134 

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

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

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

138 update_alert and cgmi_filename is None 

139 

140 # Grab p_astro values if present 

141 for p_astro in p_astros.getchildren(): 

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

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

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

145 search_for_p_astro = True 

146 

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

148 if search_for_p_astro or search_for_cgmi_file: 

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

150 

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

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

153 for message in reversed(logs): 

154 filename = message['filename'] 

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

156 filename != 'p_astro.json': 

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

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

159 for key in classifications.keys()): 

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

161 break 

162 

163 # Adjust citations if needed 

164 if p_astro_pipeline == 'rapidpe_rift': 

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

166 if len(source_classification) > 0: 

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

168 

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

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

171 if search_for_cgmi_file or cgmi_filename is not None: 

172 cgmi_json = {} 

173 # chirp mass estimates included when em-bright is 

174 if cgmi_filename is not None: 

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

176 else: 

177 for message in reversed(logs): 

178 filename = message['filename'] 

179 # use most recent mchirp estimate 

180 if filename and 'mchirp_source' in filename and \ 

181 '.json' in filename: 

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

183 break 

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

185 if cgmi_json: 

186 mchirp_bin_edges = cgmi_json['bin_edges'] 

187 mchirp_probs = cgmi_json['probabilities'] 

188 # find the highest probability bin 

189 max_prob_idx = np.argmax(mchirp_probs) 

190 left_bin = mchirp_bin_edges[max_prob_idx] 

191 right_bin = mchirp_bin_edges[max_prob_idx+1] 

192 mchirp_bin = left_bin, right_bin 

193 

194 skymaps = {} 

195 voevents_text = [] 

196 for voevent in voevents: 

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

198 if voevent == voevents[-1]: 

199 voevent_text = voevent_text_latest 

200 # If earlier voevent, load 

201 else: 

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

203 root = lxml.etree.fromstring(voevent_text) 

204 alert_type = root.find( 

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

206 if alert_type == 'earlywarning': 

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

208 early_warning_alert = True 

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

210 if url is None: 

211 continue 

212 url = url.attrib['value'] 

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

214 skyloc_pipeline = guess_skyloc_pipeline(filename) 

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

216 if filename not in skymaps: 

217 skymaps[filename] = dict( 

218 alert_type=alert_type, 

219 pipeline=skyloc_pipeline, 

220 filename=filename, 

221 latency=issued_time-event_time.gps) 

222 if skyloc_pipeline.lower() not in citation_index: 

223 citation_index[skyloc_pipeline.lower()] = \ 

224 max(citation_index.values()) + 1 

225 voevents_text.append(voevent_text) 

226 skymaps = list(skymaps.values()) 

227 

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

229 

230 kwargs = dict( 

231 subject='Identification', 

232 gracedb_id=gracedb_id, 

233 gracedb_service_url=urllib.parse.urlunsplit( 

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

235 update_alert=update_alert, 

236 cgmi_filename=cgmi_filename, 

237 group=preferred_event['group'], 

238 pipeline_search=preferred_pipeline_search, 

239 post_merger_pipeline_searches=pipeline_searches, 

240 early_warning_alert=early_warning_alert, 

241 early_warning_pipeline_searches=early_warning_pipeline_searches, 

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

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

244 far=preferred_event['far'], 

245 utctime=event_time.iso, 

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

247 skymaps=skymaps, 

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

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

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

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

252 source_classification=source_classification, 

253 mchirp_bin=mchirp_bin, 

254 include_ellipse=None, 

255 classifications=classifications, 

256 p_astro_pipeline=p_astro_pipeline, 

257 citation_index=citation_index) 

258 

259 if skymaps: 

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

261 cls = [50, 90] 

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

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

264 kwargs.update( 

265 preferred_skymap=preferred_skymap, 

266 cl=cls[-1], 

267 include_ellipse=include_ellipse, 

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

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

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

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

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

273 ellipse_area=area, 

274 greedy_area=greedy_area) 

275 try: 

276 distmu, distsig = get_distances_skymap_gracedb(gracedb_id, 

277 preferred_skymap, 

278 client) 

279 kwargs.update( 

280 distmu=distmu, 

281 distsig=distsig) 

282 except TypeError: 

283 pass 

284 

285 if raven_coinc: 

286 kwargs = _update_raven_parameters(superevent, kwargs, client, 

287 voevents_text) 

288 return kwargs 

289 

290 

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

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

293 """Compose GCN Circular draft""" 

294 

295 if client is None: 

296 client = rest.GraceDb(service) 

297 

298 kwargs = main_dict(gracedb_id, client=client) 

299 kwargs.update(authors=authors) 

300 kwargs.update(gw_is_subthreshold=False) 

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

302 

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

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

305 

306 if mailto: 

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

308 url = pattern.format( 

309 urllib.parse.quote(subject), 

310 urllib.parse.quote(body)) 

311 webbrowser.open(url) 

312 else: 

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

314 

315 

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

317 service=rest.DEFAULT_SERVICE_URL, client=None, 

318 gw_is_subthreshold=False): 

319 """Compose EM_COINC RAVEN GCN Circular draft""" 

320 

321 if client is None: 

322 client = rest.GraceDb(service) 

323 

324 kwargs = dict() 

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

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

327 kwargs.update(gw_is_subthreshold=gw_is_subthreshold) 

328 # Add RAVEN citation 

329 citation_index = kwargs['citation_index'] 

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

331 kwargs['citation_index'] = citation_index 

332 

333 subject = ( 

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

335 .strip()) 

336 body = ( 

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

338 .strip()) 

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

340 

341 

342def compose_llama( 

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

344 icecube_alert=None, remove_text_wrap=False, 

345 client=None): 

346 """Compose GRB LLAMA GCN Circular draft. 

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

348 

349 if client is None: 

350 client = rest.GraceDb(service) 

351 

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

353 

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

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

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

357 tl_datetime = str(astropy.time.Time( 

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

359 th_datetime = str(astropy.time.Time( 

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

361 

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

363 kwargs = dict( 

364 gracedb_service_url=urllib.parse.urlunsplit( 

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

366 gracedb_id=gracedb_id, 

367 llama=True, 

368 icecube_alert=icecube_alert, 

369 event_time=event_time, 

370 tl_datetime=tl_datetime, 

371 th_datetime=th_datetime, 

372 authors=authors) 

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

374 

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

376 kwargs.update(citation_index=citation_index) 

377 

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

379 

380 llama_stat_filename = 'significance_subthreshold_lvc-i3.json' 

381 if llama_stat_filename in files: 

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

383 llama_fap = llama_stat_file["p_value"] 

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

385 lines = [] 

386 for neutrino in neutrinos: 

387 # Build aligned string 

388 vals = [] 

389 dt = (event_time - 

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

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

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

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

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

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

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

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

398 

399 kwargs.update(llama_fap=llama_fap, 

400 neutrinos=lines) 

401 

402 subject = ( 

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

404 .strip()) 

405 if icecube_alert: 

406 body = ( 

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

408 .strip()) 

409 else: 

410 body = ( 

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

412 .strip()) 

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

414 

415 

416def compose_grb_medium_latency( 

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

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

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

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

421 

422 if client is None: 

423 client = rest.GraceDb(service) 

424 

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

426 search = event['search'] 

427 

428 if search != 'GRB': 

429 return 

430 

431 group = event['group'] 

432 pipeline = event['pipeline'] 

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

434 

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

436 kwargs = dict( 

437 gracedb_service_url=urllib.parse.urlunsplit( 

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

439 gracedb_id=gracedb_id, 

440 group=group, 

441 grb=True, 

442 pipeline=pipeline, 

443 external_trigger=external_trigger, 

444 exclusions=[], 

445 detections=[]) 

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

447 

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

449 

450 citation_index = {} 

451 index = 1 

452 xpipeline_fap_file = 'false_alarm_probability_x.json' 

453 if xpipeline_fap_file in files: 

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

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

456 # Create detection/exclusion circular based on given argument 

457 # Use default cutoff if not provided 

458 xpipeline_detection = (use_detection_template if use_detection_template 

459 is not None else online_xpipeline_fap < 0.001) 

460 if xpipeline_detection: 

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

462 kwargs.update(online_xpipeline_fap=online_xpipeline_fap) 

463 else: 

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

465 xpipeline_distances_file = 'distances_x.json' 

466 xpipeline_distances = client.files(gracedb_id, 

467 xpipeline_distances_file).json() 

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

469 kwargs.update(burst_exclusion=burst_exclusion) 

470 citation_index['xpipeline'] = index 

471 index += 1 

472 

473 pygrb_fap_file = 'false_alarm_probability_pygrb.json' 

474 if pygrb_fap_file in files: 

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

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

477 # Create detection/exclusion circular based on given argument 

478 # Use default cutoff if not provided 

479 pygrb_detection = (use_detection_template if use_detection_template 

480 is not None else online_pygrb_fap < 0.01) 

481 if pygrb_detection: 

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

483 kwargs.update(online_pygrb_fap=online_pygrb_fap) 

484 else: 

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

486 pygrb_distances_file = 'distances_pygrb.json' 

487 pygrb_distances = client.files(gracedb_id, 

488 pygrb_distances_file).json() 

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

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

491 kwargs.update(nsbh_exclusion=nsbh_exclusion, 

492 nsns_exclusion=nsns_exclusion) 

493 citation_index['pygrb'] = index 

494 

495 kwargs.update(citation_index=citation_index) 

496 

497 subject = ( 

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

499 .strip()) 

500 body = ( 

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

502 .strip()) 

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

504 

505 

506def compose_update(gracedb_id, authors=(), 

507 service=rest.DEFAULT_SERVICE_URL, 

508 update_types=['sky_localization', 'p_astro', 

509 'em_bright', 'raven'], 

510 remove_text_wrap=False, 

511 client=None, 

512 cgmi_filename=None): 

513 """Compose GCN update circular""" 

514 if client is None: 

515 client = rest.GraceDb(service) 

516 

517 kwargs = main_dict(gracedb_id, client=client, 

518 raven_coinc='raven' in update_types, 

519 update_alert=True, 

520 cgmi_filename=cgmi_filename) 

521 kwargs.pop('citation_index', None) 

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

523 

524 if isinstance(update_types, str): 

525 update_types = update_types.split(',') 

526 

527 # Adjust files for update type alert 

528 citation_index = {} 

529 skymaps = [] 

530 index = 1 

531 

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

533 kwargs.update(updated_skymap=updated_skymap) 

534 skymaps = [updated_skymap] 

535 if cgmi_filename: 

536 update_types.append('cgmi') 

537 if 'sky_localization' in update_types: 

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

539 index += 1 

540 if 'p_astro' in update_types and \ 

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

542 citation_index['rapidpe_rift'] = index 

543 index += 1 

544 if 'em_bright' in update_types: 

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

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

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

548 index += 1 

549 citation_index['em_bright'] = index 

550 index += 1 

551 if 'raven' in update_types: 

552 citation_index['raven'] = index 

553 

554 kwargs.update(skymaps=skymaps, 

555 citation_index=citation_index, 

556 post_merger_pipeline_searches=[], 

557 update_alert=True) 

558 

559 kwargs.update(authors=authors) 

560 kwargs.update(gw_is_subthreshold=False) 

561 kwargs.update(subject='Update') 

562 kwargs.update(update_types=update_types) 

563 

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

565 body = env.get_template( 

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

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

568 

569 

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

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

572 """Compose GCN retraction circular""" 

573 if client is None: 

574 client = rest.GraceDb(service) 

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

576 preferred_event = event['preferred_event_data'] 

577 labels = event['labels'] 

578 earlywarning = \ 

579 ('EARLY_WARNING' in labels and 

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

581 

582 kwargs = dict( 

583 subject='Retraction', 

584 gracedb_id=gracedb_id, 

585 group=preferred_event['group'], 

586 earlywarning=earlywarning, 

587 authors=authors 

588 ) 

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

590 

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

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

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

594 

595 

596def read_map_gracedb(graceid, filename, client): 

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

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

599 shutil.copyfileobj(remotefile, localfile) 

600 localfile.flush() 

601 return read_sky_map(localfile.name, moc=True) 

602 

603 

604def get_distances_skymap_gracedb(graceid, filename, client): 

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

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

607 shutil.copyfileobj(remotefile, localfile) 

608 localfile.flush() 

609 header = getheader(localfile.name, 1) 

610 try: 

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

612 except KeyError: 

613 pass 

614 

615 

616def read_map_from_path(path, client): 

617 return read_map_gracedb(*path.split('/'), client)[0] 

618 

619 

620def align_number_string(nums, positions): 

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

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

623 for i, val in enumerate(nums)) 

624 return ''.join(gen) 

625 

626 

627def mask_cl(p, level=90): 

628 pflat = p.ravel() 

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

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

631 cls = np.empty_like(pflat) 

632 cls[i] = cs 

633 cls = cls.reshape(p.shape) 

634 return cls <= 1e-2 * level 

635 

636 

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

638 """Produce table of sky map overlaps""" 

639 if client is None: 

640 client = rest.GraceDb(service) 

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

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

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

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

645 nside = hp.npix2nside(npix) 

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

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

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

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

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

651 

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

653 

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

655 

656 

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

658 ratio_ellipse_cl_areas=1.35): 

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

660 

661 Parameters 

662 ---------- 

663 graceid: str 

664 ID of the trigger used by GraceDB 

665 filename: str 

666 File name of sky map 

667 client: class 

668 REST API client for HTTP connection 

669 cls: array-like 

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

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

672 ratio_ellipse_cl_areas: float 

673 Ratio between ellipse area and minimal credible area from cl 

674 """ 

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

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

677 try: 

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

679 skymap = read_map_gracedb(graceid, new_filename, client) 

680 except (IOError, rest.HTTPError): 

681 skymap = read_map_gracedb(graceid, filename, client) 

682 else: 

683 skymap = read_map_gracedb(graceid, filename, client) 

684 

685 # Convert to an array if necessary 

686 if np.isscalar(cls): 

687 cls = [cls] 

688 cls = np.asarray(cls) 

689 

690 # Pass array of contour inteverals to get areas 

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

692 greedy_areas = np.asarray(result.contour_areas) 

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

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

695 

696 # Only use ellipse if every confidence interval passes 

697 use_ellipse = \ 

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

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

700 greedy_areas[-1]) 

701 

702 

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

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

705 

706 gracedb_id = superevent['superevent_id'] 

707 

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

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

710 gracedb_id)) 

711 

712 preferred_event = superevent['preferred_event_data'] 

713 group = preferred_event['group'] 

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

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

716 em_event_id = superevent['em_type'] 

717 

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

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

720 external_pipeline = em_event['pipeline'] 

721 # Get all other pipelines 

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

723 in superevent['em_events']] 

724 # Remove duplicates and vetoed events 

725 other_ext_pipelines = \ 

726 [*set(event['pipeline'] for event in ext_events 

727 if 'NOT_GRB' not in event['labels'])] 

728 # Remove preferred pipeline if present 

729 # This is to cover a corner case where NOT_GRB gets added to a preferred 

730 # external event after RAVEN_ALERT is applied 

731 if external_pipeline in other_ext_pipelines: 

732 other_ext_pipelines.remove(external_pipeline)