__init__.py 17.2 KB
Newer Older
1
import os
Min-A Cho's avatar
Min-A Cho committed
2
import shutil
3
import tempfile
Leo Pound Singer's avatar
Leo Pound Singer committed
4
import urllib
5 6
import webbrowser

7
import astropy.coordinates as coord
8
import astropy.time
9 10
import astropy.units as u
import healpy as hp
11
import lxml.etree
12
import numpy as np
Tri V Nguyen's avatar
Tri V Nguyen committed
13
from astropy.io.fits import getheader
14 15 16 17
from ligo.gracedb import rest
from ligo.skymap.io.fits import read_sky_map
from ligo.skymap.postprocess.ellipse import find_ellipse
from ligo.skymap.postprocess.find_injection import find_injection_moc
18 19

from .jinja import env
Leo Pound Singer's avatar
Leo Pound Singer committed
20 21 22 23
from ._version import get_versions

__version__ = get_versions()['version']
del get_versions
24 25


Leo Pound Singer's avatar
Leo Pound Singer committed
26
def 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 31 32 33 34 35 36 37 38 39


def guess_skyloc_pipeline(comments):
    comments = comments.upper()
    skyloc_pipelines = ['cWB', 'BAYESTAR', 'LIB', 'LALInference', 'UNKNOWN']
    for skyloc_pipeline in skyloc_pipelines:
        if skyloc_pipeline.upper() in comments:
            break
    return skyloc_pipeline


40 41
def main_dict(gracedb_id, client):
    """Create general dictionary to pass to compose circular"""
42

Min-A Cho's avatar
Min-A Cho committed
43 44 45
    event = client.superevent(gracedb_id).json()
    preferred_event_id = event['preferred_event']
    preferred_event = client.event(preferred_event_id).json()
46 47 48
    preferred_pipeline = preferred_event['pipeline']
    gw_events = event['gw_events']
    other_pipelines = []
49

50 51
    for gw_event in gw_events:
        pipeline = client.event(gw_event).json()['pipeline']
52
        if pipeline not in other_pipelines and pipeline != preferred_pipeline:
53
            other_pipelines.append(pipeline)
54 55
    voevents = client.voevents(gracedb_id).json()['voevents']
    log = client.logs(gracedb_id).json()['log']
56

Min-A Cho's avatar
Min-A Cho committed
57
    gpstime = float(preferred_event['gpstime'])
58 59
    event_time = astropy.time.Time(gpstime, format='gps').utc

60
    skymaps = {}
61
    for voevent in voevents:
62 63
        voevent_text = client.files(gracedb_id, voevent['filename']).read()
        root = lxml.etree.fromstring(voevent_text)
64 65
        alert_type = root.find(
            './What/Param[@name="AlertType"]').attrib['value'].lower()
Min-A Cho's avatar
Min-A Cho committed
66
        url = root.find('./What/Group/Param[@name="skymap_fits"]')
67 68 69 70
        if url is None:
            continue
        url = url.attrib['value']
        _, filename = os.path.split(url)
71
        skyloc_pipeline = guess_skyloc_pipeline(filename)
72
        issued_time = astropy.time.Time(root.find('./Who/Date').text).datetime
73
        if filename not in skymaps:
Min-A Cho's avatar
Min-A Cho committed
74 75 76
            for message in log:
                if filename == message['filename']:
                    tag_names = message['tag_names']
77
                    if 'sky_loc' in tag_names and 'public' in tag_names:
Min-A Cho's avatar
Min-A Cho committed
78 79 80 81 82 83
                        skymaps[filename] = dict(
                            alert_type=alert_type,
                            pipeline=skyloc_pipeline,
                            filename=filename,
                            latency=issued_time-event_time.datetime)
                    break
84
    skymaps = list(skymaps.values())
85
    preferred_event_files = client.files(preferred_event_id).json()
86
    em_brightfile = 'em_bright.json'
87
    if em_brightfile in preferred_event_files:
88
        source_classification = client.files(
89
            preferred_event_id, em_brightfile).json()
90 91 92
        source_classification = {
            key: 100 * value for key, value in source_classification.items()}

93 94 95
    else:
        source_classification = {}

96
    # adding the p_atro informations if available
97
    p_astro_file = 'p_astro.json'
98
    if p_astro_file in preferred_event_files:
99
        classifications = client.files(preferred_event_id, p_astro_file).json()
100

101 102 103
        # Convert to percent for consistency with em_bright
        classifications = {
            key: 100 * value for key, value in classifications.items()}
104
    else:
105
        classifications = {}
106

Min-A Cho's avatar
Min-A Cho committed
107 108
    o = urllib.parse.urlparse(client.service_url)

109
    kwargs = dict(
110
        subject='Identification',
111
        gracedb_id=gracedb_id,
Min-A Cho's avatar
Min-A Cho committed
112
        gracedb_service_url=urllib.parse.urlunsplit(
Min-A Cho's avatar
Min-A Cho committed
113 114
            (o.scheme, o.netloc, '/superevents/', '', '')),
        group=preferred_event['group'],
115 116
        pipeline=preferred_pipeline,
        other_pipelines=other_pipelines,
Min-A Cho's avatar
Min-A Cho committed
117
        all_pipelines=[preferred_pipeline] + other_pipelines,
Min-A Cho's avatar
Min-A Cho committed
118 119 120
        gpstime='{0:.03f}'.format(round(float(preferred_event['gpstime']), 3)),
        search=preferred_event.get('search', ''),
        far=preferred_event['far'],
121
        utctime=event_time.iso,
Min-A Cho's avatar
Min-A Cho committed
122
        instruments=preferred_event['instruments'].split(','),
123
        skymaps=skymaps,
124 125
        prob_has_ns=source_classification.get('HasNS'),
        prob_has_remnant=source_classification.get('HasRemnant'),
126
        include_ellipse=None,
127
        classifications=classifications)
128 129 130 131

    if skymaps:
        preferred_skymap = skymaps[-1]['filename']
        cl = 90
132 133
        include_ellipse, ra, dec, a, b, pa, area, greedy_area = \
            uncertainty_ellipse(gracedb_id, preferred_skymap, client, cl=cl)
134 135 136 137 138 139 140 141 142
        kwargs.update(
            preferred_skymap=preferred_skymap,
            cl=cl,
            include_ellipse=include_ellipse,
            ra=coord.Longitude(ra*u.deg),
            dec=coord.Latitude(dec*u.deg),
            a=coord.Angle(a*u.deg),
            b=coord.Angle(b*u.deg),
            pa=coord.Angle(pa*u.deg),
143 144
            ellipse_area=area,
            greedy_area=greedy_area)
145 146 147 148 149 150 151 152 153
        try:
            distmu, distsig = get_distances_skymap_gracedb(gracedb_id,
                                                           preferred_skymap,
                                                           client)
            kwargs.update(
                distmu=distmu,
                distsig=distsig)
        except TypeError:
            pass
154

155 156 157 158 159 160 161 162 163 164 165
    return kwargs


def compose(gracedb_id, authors=(), mailto=False,
            service=rest.DEFAULT_SERVICE_URL, client=None):
    """Compose GCN Circular draft"""

    if client is None:
        client = rest.GraceDb(service)

    kwargs = main_dict(gracedb_id, client=client)
166 167
    kwargs.update(authors=authors)
    kwargs.update(change_significance_statement=False)
168

169
    subject = env.get_template('subject.jinja2').render(**kwargs).strip()
170
    body = env.get_template('initial_circular.jinja2').render(**kwargs).strip()
171

172
    if mailto:
173 174 175 176 177 178
        pattern = 'mailto:emfollow@ligo.org,dac@ligo.org?subject={0}&body={1}'
        url = pattern.format(
            urllib.parse.quote(subject),
            urllib.parse.quote(body))
        webbrowser.open(url)
    else:
179
        return '{0}\n{1}'.format(subject, body)
180 181


182
def compose_raven(gracedb_id, authors=(),
183 184 185 186 187
                  service=rest.DEFAULT_SERVICE_URL, client=None):
    """Compose EM_COINC RAVEN GCN Circular draft"""

    if client is None:
        client = rest.GraceDb(service)
188

189
    event = client.superevent(gracedb_id).json()
190

191
    if 'EM_COINC' not in event['labels']:
192
        return
193

194
    preferred_event_id = event['preferred_event']
195 196
    preferred_event = client.event(preferred_event_id).json()
    group = preferred_event['group']
197 198
    superevent_far = preferred_event['far']
    gpstime = float(preferred_event['gpstime'])
199

200
    for em_event_id in event['em_events']:
201
        em_event = client.event(em_event_id).json()
202
        em_event_gpstime = float(em_event['gpstime'])
203
        external_pipeline = em_event['pipeline']
204
        # FIXME in GraceDb: Even SNEWS triggers have an extra attribute GRB.
Min-A Cho's avatar
Min-A Cho committed
205
        external_trigger_id = em_event['extra_attributes']['GRB']['trigger_id']
206
        snews = (em_event['search'] == 'Supernova')
207 208
        grb = em_event['search'] in ['GRB', 'SubGRB']
        subthreshold = (em_event['search'] == 'SubGRB')
209

210
    o = urllib.parse.urlparse(client.service_url)
211 212 213
    kwargs = dict(
        gracedb_service_url=urllib.parse.urlunsplit(
            (o.scheme, o.netloc, '/superevents/', '', '')),
214 215 216
        gracedb_id=gracedb_id,
        group=group,
        superevent_far=preferred_event['far'],
217
        external_pipeline=external_pipeline,
Min-A Cho's avatar
Min-A Cho committed
218
        external_trigger=external_trigger_id,
219 220
        snews=snews,
        grb=grb,
221
        subthreshold=subthreshold,
222 223
        latency=abs(round(em_event_gpstime-gpstime, 1)),
        beforeafter='before' if gpstime > em_event_gpstime else 'after')
224

225
    if grb:
226 227 228 229 230
        # Grab GRB coincidence FARs
        files = client.files(gracedb_id).json()
        coinc_far_file = 'coincidence_far.json'
        if coinc_far_file in files:
            coincidence_far = client.files(gracedb_id, coinc_far_file).json()
231 232
            time_coinc_far = coincidence_far.get('temporal_coinc_far')
            space_time_coinc_far = coincidence_far.get(
233 234
                'spatiotemporal_coinc_far')
            kwargs.update(
235 236
                time_coinc_far=time_coinc_far,
                space_time_coinc_far=space_time_coinc_far)
237

238
        # Find combined sky map for GRB
239 240 241 242 243 244
        logs = client.logs(gracedb_id).json()['log']
        for message in reversed(logs):
            comment = message['comment']
            filename = message['filename']
            if (filename.endswith('-gbm.fits.gz')
                    or 'LVC-Fermi sky map' in comment):
245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
                cl = 90
                include_ellipse, ra, dec, a, b, pa, area, greedy_area = \
                    uncertainty_ellipse(gracedb_id, filename, client, cl=cl)
                kwargs.update(
                    combined_skymap=filename,
                    cl=cl,
                    combined_skymap_include_ellipse=include_ellipse,
                    combined_skymap_ra=coord.Longitude(ra*u.deg),
                    combined_skymap_dec=coord.Latitude(dec*u.deg),
                    combined_skymap_a=coord.Angle(a*u.deg),
                    combined_skymap_b=coord.Angle(b*u.deg),
                    combined_skymap_pa=coord.Angle(pa*u.deg),
                    combined_skymap_ellipse_area=area,
                    combined_skymap_greedy_area=greedy_area)

    kwargs.update(main_dict(gracedb_id, client=client))
261 262
    if (group.lower() == 'cbc' and superevent_far < 1 / (60 * 86400)) or \
            (group.lower() == 'burst' and superevent_far < 1 / (365 * 86400)):
263
        kwargs.update(change_significance_statement=False)
264
    else:
265
        kwargs.update(change_significance_statement=True)
266 267 268 269 270 271 272 273 274 275

    subject = (
        env.get_template('RAVEN_subject.jinja2').render(**kwargs)
        .strip())
    body = (
        env.get_template('RAVEN_circular.jinja2').render(**kwargs)
        .strip())
    return '{0}\n{1}'.format(subject, body)


276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291
def compose_grb_medium_latency(
        gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL, client=None):
    """Compose GRB Medium Latency GCN Circular draft.
    Here, gracedb_id will be a GRB external trigger ID in GraceDb."""

    if client is None:
        client = rest.GraceDb(service)

    event = client.event(gracedb_id).json()
    search = event['search']

    if search != 'GRB':
        return

    group = event['group']
    pipeline = event['pipeline']
Min-A Cho's avatar
Min-A Cho committed
292
    external_trigger = event['extra_attributes']['GRB']['trigger_id']
293 294 295 296 297 298 299

    o = urllib.parse.urlparse(client.service_url)
    kwargs = dict(
        gracedb_service_url=urllib.parse.urlunsplit(
            (o.scheme, o.netloc, '/events/', '', '')),
        gracedb_id=gracedb_id,
        group=group,
300
        grb=True,
301
        pipeline=pipeline,
Min-A Cho's avatar
Min-A Cho committed
302
        external_trigger=external_trigger,
303 304 305 306 307 308 309 310 311 312 313 314 315 316
        exclusions=[],
        detections=[])

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

    xpipeline_fap_file = 'false_alarm_probability_x.json'
    if xpipeline_fap_file in files:
        xpipeline_fap = client.files(gracedb_id, xpipeline_fap_file).json()
        online_xpipeline_fap = xpipeline_fap.get('Online Xpipeline')
        if online_xpipeline_fap < 0.001:
            kwargs['detections'].append('Burst')
            kwargs.update(online_xpipeline_fap=online_xpipeline_fap)
        else:
            kwargs['exclusions'].append('Burst')
317 318 319 320 321
            xpipeline_distances_file = 'distances_x.json'
            xpipeline_distances = client.files(gracedb_id,
                                               xpipeline_distances_file).json()
            burst_exclusion = xpipeline_distances.get('Burst Exclusion')
            kwargs.update(burst_exclusion=burst_exclusion)
322 323 324 325 326 327 328 329 330 331

    pygrb_fap_file = 'false_alarm_probability_pygrb.json'
    if pygrb_fap_file in files:
        pygrb_fap = client.files(gracedb_id, pygrb_fap_file).json()
        online_pygrb_fap = pygrb_fap.get('Online PyGRB')
        if online_pygrb_fap < 0.01:
            kwargs['detections'].append('CBC')
            kwargs.update(online_pygrb_fap=online_pygrb_fap)
        else:
            kwargs['exclusions'].append('CBC')
332 333 334 335 336 337 338
            pygrb_distances_file = 'distances_pygrb.json'
            pygrb_distances = client.files(gracedb_id,
                                           pygrb_distances_file).json()
            nsns_exclusion = pygrb_distances.get('NSNS Exclusion')
            nsbh_exclusion = pygrb_distances.get('NSBH Exclusion')
            kwargs.update(nsbh_exclusion=nsbh_exclusion,
                          nsns_exclusion=nsns_exclusion)
339 340 341 342 343 344 345 346 347 348

    subject = (
        env.get_template('medium_latency_grb_subject.jinja2').render(**kwargs)
        .strip())
    body = (
        env.get_template('medium_latency_grb_circular.jinja2').render(**kwargs)
        .strip())
    return '{0}\n{1}'.format(subject, body)


349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368
def compose_update(gracedb_id, authors=(),
                   service=rest.DEFAULT_SERVICE_URL,
                   update_types=['sky_localization', 'p_astro', 'em_bright'],
                   client=None):
    """Compose GCN update circular"""
    if client is None:
        client = rest.GraceDb(service)

    kwargs = main_dict(gracedb_id, client=client)
    kwargs.update(authors=authors)
    kwargs.update(change_significance_statement=False)
    kwargs.update(subject='Update')
    kwargs.update(update_types=update_types)

    subject = env.get_template('subject.jinja2').render(**kwargs).strip()
    body = env.get_template(
               'update_circular.jinja2').render(**kwargs).strip()
    return '{0}\n{1}'.format(subject, body)


369 370 371 372 373
def compose_retraction(gracedb_id, authors=(),
                       service=rest.DEFAULT_SERVICE_URL, client=None):
    """Compose GCN retraction circular"""
    if client is None:
        client = rest.GraceDb(service)
374 375
    event = client.superevent(gracedb_id).json()
    preferred_event = client.event(event['preferred_event']).json()
376 377 378 379 380 381 382 383 384 385 386 387 388

    kwargs = dict(
                 subject='Retraction',
                 gracedb_id=gracedb_id,
                 group=preferred_event['group'],
                 authors=authors
             )

    subject = env.get_template('subject.jinja2').render(**kwargs).strip()
    body = env.get_template('retraction.jinja2').render(**kwargs).strip()
    return '{0}\n{1}'.format(subject, body)


Min-A Cho's avatar
Min-A Cho committed
389
def read_map_gracedb(graceid, filename, client):
390
    with tempfile.NamedTemporaryFile(mode='w+b') as localfile:
Min-A Cho's avatar
Min-A Cho committed
391 392
        remotefile = client.files(graceid, filename, raw=True)
        shutil.copyfileobj(remotefile, localfile)
393
        localfile.flush()
Min-A Cho's avatar
Min-A Cho committed
394
        return read_sky_map(localfile.name, moc=True)
395 396


397 398
def get_distances_skymap_gracedb(graceid, filename, client):
    with tempfile.NamedTemporaryFile(mode='w+b') as localfile:
Tri V Nguyen's avatar
Tri V Nguyen committed
399
        remotefile = client.files(graceid, filename, raw=True)
400 401
        shutil.copyfileobj(remotefile, localfile)
        localfile.flush()
Tri V Nguyen's avatar
Tri V Nguyen committed
402
        header = getheader(localfile.name, 1)
403 404 405 406 407 408
        try:
            return header['distmean'], header['diststd']
        except KeyError:
            pass


409
def read_map_from_path(path, client):
410
    return read_map_gracedb(*path.split('/'), client)[0]
411 412


413 414 415 416 417 418 419 420 421 422
def mask_cl(p, level=90):
    pflat = p.ravel()
    i = np.flipud(np.argsort(p))
    cs = np.cumsum(pflat[i])
    cls = np.empty_like(pflat)
    cls[i] = cs
    cls = cls.reshape(p.shape)
    return cls <= 1e-2 * level


423
def compare_skymaps(paths, service=rest.DEFAULT_SERVICE_URL, client=None):
424
    """Produce table of sky map overlaps"""
425 426
    if client is None:
        client = rest.GraceDb(service)
427
    filenames = [path.split('/')[1] for path in paths]
428
    pipelines = [guess_skyloc_pipeline(filename) for filename in filenames]
429
    probs = [read_map_from_path(path, client) for path in paths]
430 431 432 433 434 435 436 437 438 439
    npix = max(len(prob) for prob in probs)
    nside = hp.npix2nside(npix)
    deg2perpix = hp.nside2pixarea(nside, degrees=True)
    probs = [hp.ud_grade(prob, nside, power=-2) for prob in probs]
    masks = [mask_cl(prob) for prob in probs]
    areas = [mask.sum() * deg2perpix for mask in masks]
    joint_areas = [(mask & masks[-1]).sum() * deg2perpix for mask in masks]

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

440
    return env.get_template('compare_skymaps.jinja2').render(**kwargs)
441 442 443


def uncertainty_ellipse(graceid, filename, client, cl=90):
444
    """Compute uncertainty ellipses for a given sky map"""
Tri V Nguyen's avatar
Tri V Nguyen committed
445
    if filename.endswith('.gz'):
446
        # Try using the multi-res sky map if it exists
447
        try:
Tri V Nguyen's avatar
Tri V Nguyen committed
448 449
            new_filename = filename.replace('.fits.gz', '.multiorder.fits')
            skymap = read_map_gracedb(graceid, new_filename, client)
450
        except (IOError, rest.HTTPError):
Min-A Cho's avatar
Min-A Cho committed
451
            skymap = read_map_gracedb(graceid, filename, client)
452
    else:
Min-A Cho's avatar
Min-A Cho committed
453 454 455
        skymap = read_map_gracedb(graceid, filename, client)

    result = find_injection_moc(skymap, contours=[cl / 100])
456
    greedy_area = result.contour_areas[0]
Min-A Cho's avatar
Min-A Cho committed
457
    ra, dec, a, b, pa, ellipse_area = find_ellipse(skymap, cl=cl)
458 459
    return ellipse_area <= 1.35*greedy_area, ra, dec, a, b, pa, \
        ellipse_area, greedy_area