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 } 

41 try: 

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

43 except KeyError: 

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

45 

46 

47def text_width(remove_text_wrap): 

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

49 not.""" 

50 return 9999 if remove_text_wrap else 79 

51 

52 

53def main_dict(gracedb_id, client, raven_coinc=False): 

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

55 

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

57 preferred_event = superevent['preferred_event_data'] 

58 pipeline = preferred_event['pipeline'] 

59 search = preferred_event['search'] 

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

61 early_warning_pipeline_searches = [] 

62 pipeline_searches = [] 

63 gw_events = superevent['gw_events'] 

64 early_warning_alert = False 

65 

66 for gw_event in gw_events: 

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

68 pipeline = gw_event_dict['pipeline'] 

69 search = gw_event_dict['search'] 

70 # Remap MDC to allsky to not double count 

71 if search == 'MDC': 

72 search = 'AllSky' 

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

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

75 if pipeline_search not in pipeline_searches and \ 

76 'EarlyWarning' not in pipeline_search: 

77 pipeline_searches.append(pipeline_search) 

78 if pipeline_search not in early_warning_pipeline_searches and \ 

79 'EarlyWarning' in pipeline_search: 

80 # Remap MBTA early warning to the all sky citation 

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

82 pipeline_search = 'MBTA_AllSky' 

83 early_warning_pipeline_searches.append(pipeline_search) 

84 # Sort to get alphabetical order 

85 pipeline_searches.sort(key=str.lower) 

86 early_warning_pipeline_searches.sort(key=str.lower) 

87 

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

89 if not voevents: 

90 raise ValueError( 

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

92 gracedb_id)) 

93 

94 citation_index = { 

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

96 for pipeline_search in pipeline_searches} 

97 for pipeline_search in early_warning_pipeline_searches: 

98 if pipeline_search.lower() not in citation_index: 

99 citation_index[pipeline_search.lower()] = \ 

100 max(citation_index.values()) + 1 

101 

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

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

104 

105 # Grab latest p_astro and em_bright values from lastest VOEvent 

106 voevent_text_latest = \ 

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

108 root = lxml.etree.fromstring(voevent_text_latest) 

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

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

111 classifications = {} 

112 source_classification = {} 

113 mchirp_bin = {} 

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

115 p_astro_pipeline = None 

116 

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

118 if p_astros or em_brights: 

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

120 

121 if p_astros: 

122 for p_astro in p_astros.getchildren(): 

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

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

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

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

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

128 for message in reversed(logs): 

129 filename = message['filename'] 

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

131 filename != 'p_astro.json': 

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

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

134 for key in classifications.keys()): 

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

136 break 

137 if p_astro_pipeline == 'rapidpe_rift': 

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

139 if em_brights: 

140 for em_bright in em_brights.getchildren(): 

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

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

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

144 if len(source_classification) > 0: 

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

146 

147 # chirp mass estimates included when em-bright is 

148 for message in reversed(logs): 

149 filename = message['filename'] 

150 # use most recent mchirp estimate 

151 if filename and 'mchirp_source' in filename: 

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

153 mchirp_bin_edges = cgmi_json['bin_edges'] 

154 mchirp_probs = cgmi_json['probabilities'] 

155 # find the highest probability bin 

156 max_prob_idx = np.argmax(mchirp_probs) 

157 left_bin = mchirp_bin_edges[max_prob_idx] 

158 right_bin = mchirp_bin_edges[max_prob_idx+1] 

159 mchirp_bin = left_bin, right_bin 

160 break 

161 

162 skymaps = {} 

163 voevents_text = [] 

164 for voevent in voevents: 

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

166 if voevent == voevents[-1]: 

167 voevent_text = voevent_text_latest 

168 # If earlier voevent, load 

169 else: 

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

171 root = lxml.etree.fromstring(voevent_text) 

172 alert_type = root.find( 

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

174 if alert_type == 'earlywarning': 

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

176 early_warning_alert = True 

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

178 if url is None: 

179 continue 

180 url = url.attrib['value'] 

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

182 skyloc_pipeline = guess_skyloc_pipeline(filename) 

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

184 if filename not in skymaps: 

185 skymaps[filename] = dict( 

186 alert_type=alert_type, 

187 pipeline=skyloc_pipeline, 

188 filename=filename, 

189 latency=issued_time-event_time.gps) 

190 if skyloc_pipeline.lower() not in citation_index: 

191 citation_index[skyloc_pipeline.lower()] = \ 

192 max(citation_index.values()) + 1 

193 voevents_text.append(voevent_text) 

194 skymaps = list(skymaps.values()) 

195 

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

197 

198 kwargs = dict( 

199 subject='Identification', 

200 gracedb_id=gracedb_id, 

201 gracedb_service_url=urllib.parse.urlunsplit( 

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

203 group=preferred_event['group'], 

204 pipeline_search=preferred_pipeline_search, 

205 all_pipeline_searches=pipeline_searches, 

206 early_warning_alert=early_warning_alert, 

207 early_warning_pipeline_searches=early_warning_pipeline_searches, 

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

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

210 far=preferred_event['far'], 

211 utctime=event_time.iso, 

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

213 skymaps=skymaps, 

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

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

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

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

218 source_classification=source_classification, 

219 mchirp_bin=mchirp_bin, 

220 include_ellipse=None, 

221 classifications=classifications, 

222 p_astro_pipeline=p_astro_pipeline, 

223 citation_index=citation_index) 

224 

225 if skymaps: 

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

227 cls = [50, 90] 

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

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

230 kwargs.update( 

231 preferred_skymap=preferred_skymap, 

232 cl=cls[-1], 

233 include_ellipse=include_ellipse, 

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

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

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

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

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

239 ellipse_area=area, 

240 greedy_area=greedy_area) 

241 try: 

242 distmu, distsig = get_distances_skymap_gracedb(gracedb_id, 

243 preferred_skymap, 

244 client) 

245 kwargs.update( 

246 distmu=distmu, 

247 distsig=distsig) 

248 except TypeError: 

249 pass 

250 

251 if raven_coinc: 

252 kwargs = _update_raven_parameters(superevent, kwargs, client, 

253 voevents_text) 

254 return kwargs 

255 

256 

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

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

259 """Compose GCN Circular draft""" 

260 

261 if client is None: 

262 client = rest.GraceDb(service) 

263 

264 kwargs = main_dict(gracedb_id, client=client) 

265 kwargs.update(authors=authors) 

266 kwargs.update(change_significance_statement=False) 

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

268 

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

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

271 

272 if mailto: 

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

274 url = pattern.format( 

275 urllib.parse.quote(subject), 

276 urllib.parse.quote(body)) 

277 webbrowser.open(url) 

278 else: 

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

280 

281 

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

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

284 """Compose EM_COINC RAVEN GCN Circular draft""" 

285 

286 if client is None: 

287 client = rest.GraceDb(service) 

288 

289 kwargs = dict() 

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

291 kwargs.update(update_alert=False) 

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

293 # Add RAVEN citation 

294 citation_index = kwargs['citation_index'] 

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

296 kwargs['citation_index'] = citation_index 

297 

298 subject = ( 

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

300 .strip()) 

301 body = ( 

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

303 .strip()) 

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

305 

306 

307def compose_llama( 

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

309 icecube_alert=None, remove_text_wrap=False, 

310 client=None): 

311 """Compose GRB LLAMA GCN Circular draft. 

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

313 

314 if client is None: 

315 client = rest.GraceDb(service) 

316 

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

318 

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

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

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

322 tl_datetime = str(astropy.time.Time( 

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

324 th_datetime = str(astropy.time.Time( 

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

326 

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

328 kwargs = dict( 

329 gracedb_service_url=urllib.parse.urlunsplit( 

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

331 gracedb_id=gracedb_id, 

332 llama=True, 

333 icecube_alert=icecube_alert, 

334 event_time=event_time, 

335 tl_datetime=tl_datetime, 

336 th_datetime=th_datetime, 

337 authors=authors) 

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

339 

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

341 kwargs.update(citation_index=citation_index) 

342 

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

344 

345 llama_stat_filename = 'significance_subthreshold_lvc-i3.json' 

346 if llama_stat_filename in files: 

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

348 llama_fap = llama_stat_file["p_value"] 

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

350 lines = [] 

351 for neutrino in neutrinos: 

352 # Build aligned string 

353 vals = [] 

354 dt = (event_time - 

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

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

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

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

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

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

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

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

363 

364 kwargs.update(llama_fap=llama_fap, 

365 neutrinos=lines) 

366 

367 subject = ( 

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

369 .strip()) 

370 if icecube_alert: 

371 body = ( 

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

373 .strip()) 

374 else: 

375 body = ( 

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

377 .strip()) 

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

379 

380 

381def compose_grb_medium_latency( 

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

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

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

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

386 

387 if client is None: 

388 client = rest.GraceDb(service) 

389 

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

391 search = event['search'] 

392 

393 if search != 'GRB': 

394 return 

395 

396 group = event['group'] 

397 pipeline = event['pipeline'] 

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

399 

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

401 kwargs = dict( 

402 gracedb_service_url=urllib.parse.urlunsplit( 

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

404 gracedb_id=gracedb_id, 

405 group=group, 

406 grb=True, 

407 pipeline=pipeline, 

408 external_trigger=external_trigger, 

409 exclusions=[], 

410 detections=[]) 

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

412 

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

414 

415 citation_index = {} 

416 index = 1 

417 xpipeline_fap_file = 'false_alarm_probability_x.json' 

418 if xpipeline_fap_file in files: 

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

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

421 # Create detection/exclusion circular based on given argument 

422 # Use default cutoff if not provided 

423 xpipeline_detection = (use_detection_template if use_detection_template 

424 is not None else online_xpipeline_fap < 0.001) 

425 if xpipeline_detection: 

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

427 kwargs.update(online_xpipeline_fap=online_xpipeline_fap) 

428 else: 

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

430 xpipeline_distances_file = 'distances_x.json' 

431 xpipeline_distances = client.files(gracedb_id, 

432 xpipeline_distances_file).json() 

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

434 kwargs.update(burst_exclusion=burst_exclusion) 

435 citation_index['xpipeline'] = index 

436 index += 1 

437 

438 pygrb_fap_file = 'false_alarm_probability_pygrb.json' 

439 if pygrb_fap_file in files: 

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

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

442 # Create detection/exclusion circular based on given argument 

443 # Use default cutoff if not provided 

444 pygrb_detection = (use_detection_template if use_detection_template 

445 is not None else online_pygrb_fap < 0.01) 

446 if pygrb_detection: 

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

448 kwargs.update(online_pygrb_fap=online_pygrb_fap) 

449 else: 

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

451 pygrb_distances_file = 'distances_pygrb.json' 

452 pygrb_distances = client.files(gracedb_id, 

453 pygrb_distances_file).json() 

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

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

456 kwargs.update(nsbh_exclusion=nsbh_exclusion, 

457 nsns_exclusion=nsns_exclusion) 

458 citation_index['pygrb'] = index 

459 

460 kwargs.update(citation_index=citation_index) 

461 

462 subject = ( 

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

464 .strip()) 

465 body = ( 

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

467 .strip()) 

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

469 

470 

471def compose_update(gracedb_id, authors=(), 

472 service=rest.DEFAULT_SERVICE_URL, 

473 update_types=['sky_localization', 'p_astro', 

474 'em_bright', 'raven'], 

475 remove_text_wrap=False, 

476 client=None): 

477 """Compose GCN update circular""" 

478 if client is None: 

479 client = rest.GraceDb(service) 

480 

481 kwargs = main_dict(gracedb_id, client=client, 

482 raven_coinc='raven' in update_types) 

483 kwargs.pop('citation_index', None) 

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

485 

486 if isinstance(update_types, str): 

487 update_types = update_types.split(',') 

488 

489 # Adjust files for update type alert 

490 citation_index = {} 

491 skymaps = [] 

492 index = 1 

493 

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

495 kwargs.update(updated_skymap=updated_skymap) 

496 skymaps = [updated_skymap] 

497 if 'sky_localization' in update_types: 

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

499 index += 1 

500 if 'p_astro' in update_types and \ 

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

502 citation_index['rapidpe_rift'] = index 

503 index += 1 

504 if 'em_bright' in update_types: 

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

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

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

508 index += 1 

509 citation_index['em_bright'] = index 

510 index += 1 

511 if 'raven' in update_types: 

512 citation_index['raven'] = index 

513 

514 kwargs.update(skymaps=skymaps, 

515 citation_index=citation_index, 

516 all_pipeline_searches=[], 

517 update_alert=True) 

518 

519 kwargs.update(authors=authors) 

520 kwargs.update(change_significance_statement=False) 

521 kwargs.update(subject='Update') 

522 kwargs.update(update_types=update_types) 

523 

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

525 body = env.get_template( 

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

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

528 

529 

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

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

532 """Compose GCN retraction circular""" 

533 if client is None: 

534 client = rest.GraceDb(service) 

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

536 preferred_event = event['preferred_event_data'] 

537 labels = event['labels'] 

538 earlywarning = \ 

539 ('EARLY_WARNING' in labels and 

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

541 

542 kwargs = dict( 

543 subject='Retraction', 

544 gracedb_id=gracedb_id, 

545 group=preferred_event['group'], 

546 earlywarning=earlywarning, 

547 authors=authors 

548 ) 

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

550 

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

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

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

554 

555 

556def read_map_gracedb(graceid, filename, client): 

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

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

559 shutil.copyfileobj(remotefile, localfile) 

560 localfile.flush() 

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

562 

563 

564def get_distances_skymap_gracedb(graceid, filename, client): 

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

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

567 shutil.copyfileobj(remotefile, localfile) 

568 localfile.flush() 

569 header = getheader(localfile.name, 1) 

570 try: 

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

572 except KeyError: 

573 pass 

574 

575 

576def read_map_from_path(path, client): 

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

578 

579 

580def align_number_string(nums, positions): 

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

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

583 for i, val in enumerate(nums)) 

584 return ''.join(gen) 

585 

586 

587def mask_cl(p, level=90): 

588 pflat = p.ravel() 

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

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

591 cls = np.empty_like(pflat) 

592 cls[i] = cs 

593 cls = cls.reshape(p.shape) 

594 return cls <= 1e-2 * level 

595 

596 

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

598 """Produce table of sky map overlaps""" 

599 if client is None: 

600 client = rest.GraceDb(service) 

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

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

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

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

605 nside = hp.npix2nside(npix) 

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

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

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

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

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

611 

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

613 

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

615 

616 

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

618 ratio_ellipse_cl_areas=1.35): 

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

620 

621 Parameters 

622 ---------- 

623 graceid: str 

624 ID of the trigger used by GraceDB 

625 filename: str 

626 File name of sky map 

627 client: class 

628 REST API client for HTTP connection 

629 cls: array-like 

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

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

632 ratio_ellipse_cl_areas: float 

633 Ratio between ellipse area and minimal credible area from cl 

634 """ 

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

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

637 try: 

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

639 skymap = read_map_gracedb(graceid, new_filename, client) 

640 except (IOError, rest.HTTPError): 

641 skymap = read_map_gracedb(graceid, filename, client) 

642 else: 

643 skymap = read_map_gracedb(graceid, filename, client) 

644 

645 # Convert to an array if necessary 

646 if np.isscalar(cls): 

647 cls = [cls] 

648 cls = np.asarray(cls) 

649 

650 # Pass array of contour inteverals to get areas 

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

652 greedy_areas = np.asarray(result.contour_areas) 

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

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

655 

656 # Only use ellipse if every confidence interval passes 

657 use_ellipse = \ 

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

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

660 greedy_areas[-1]) 

661 

662 

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

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

665 

666 gracedb_id = superevent['superevent_id'] 

667 

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

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

670 gracedb_id)) 

671 

672 preferred_event = superevent['preferred_event_data'] 

673 group = preferred_event['group'] 

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

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

676 em_event_id = superevent['em_type'] 

677 

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

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

680 external_pipeline = em_event['pipeline'] 

681 # Get all other pipelines 

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

683 in superevent['em_events']] 

684 # Remove duplicates and vetoed events 

685 other_ext_pipelines = \ 

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

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

688 # Remove preferred pipeline 

689 other_ext_pipelines.remove(external_pipeline) 

690 # FIXME in GraceDb: Even SNEWS triggers have an extra attribute GRB. 

691 external_trigger_id = em_event['extra_attributes']['GRB']['trigger_id'] 

692 snews = (em_event['pipeline'] == 'SNEWS') 

693 grb = (em_event['search'] in ['GRB', 'SubGRB', 'SubGRBTargeted', 'MDC'] 

694 and not snews) 

695 subthreshold = em_event['search'] in ['SubGRB', 'SubGRBTargeted'] 

696 subthreshold_targeted = em_event['search'] == 'SubGRBTargeted' 

697 far_grb = em_event['far'] 

698 

699 voevent_text_latest = voevents_text[-1] 

700 root = lxml.etree.fromstring(voevent_text_latest) 

701 time_diff = \ 

702 root.find('./What/Group/Param[@name="Time_Difference"]') 

703 time_diff = float(time_diff.attrib['value']) 

704 

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

706 kwargs.update( 

707 gracedb_service_url=urllib.parse.urlunsplit( 

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

709 gracedb_id=gracedb_id, 

710 group=group, 

711 external_pipeline=external_pipeline, 

712 external_trigger=external_trigger_id, 

713 snews=snews, 

714 grb=grb, 

715 subthreshold=subthreshold, 

716 subthreshold_targeted=subthreshold_targeted, 

717 other_ext_pipelines=sorted(other_ext_pipelines), 

718 far_grb=far_grb, 

719 latency=abs(round(time_diff, 1)),