WIP: Create initial circular from VOEvent, not GraceDb

This makes it easy to generate sample circulars from the notices
in the User Guide.

Still to do:

* Port remaining circular types.
* Add support for multiple VOEvents so that our update circulars
  can describe the history of the event (including, but not limited
  to, changes to the localization).

ligo/followup_advocate/__init__.py

@@ -17,6 +17,7 @@ from ligo.skymap.postprocess.find_injection import find_injection_moc
 
 from .jinja import env
 from .version import __version__  # noqa
+from .voevent import load
 
 
 def authors(authors, service=rest.DEFAULT_SERVICE_URL):
@@ -147,16 +148,18 @@ def main_dict(gracedb_id, client):
     return kwargs
 
 
-def compose(gracedb_id, authors=(), mailto=False,
+def compose(voevent, 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)
+    # Load Jinja variables from VOEvent params.
+    kwargs = load(voevent)
+
+    # Additional keywrod arguments.
     kwargs.update(authors=authors)
-    kwargs.update(change_significance_statement=False)
 
     subject = env.get_template('subject.jinja2').render(**kwargs).strip()
     body = env.get_template('initial_circular.jinja2').render(**kwargs).strip()

ligo/followup_advocate/templates/initial_body.jinja2

@@ -4,10 +4,10 @@
                             citations, renamegroup %}
 
 {% filter rewrap %}
-The {{pipeline}} {{renamegroup(group)}} analysis {% if pipeline|lower in citations %}
-({{citations[pipeline|lower]}}){% endif %} identified candidate
-{{gracedb_id}} during real-time processing of data from
-{{naturalinstruments(instruments)}} at {{utctime}} UTC (GPS time: {{gpstime}}).
+The {{Pipeline}} {{renamegroup(Group)}} analysis {% if Pipeline|lower in citations %}
+({{citations[Pipeline|lower]}}){% endif %} identified candidate
+{{GraceID}} during real-time processing of data from
+{{naturalinstruments(Instruments.split(','))}} at {{DateObs.iso}} UTC (GPS time: {{'%.3f'|format(DateObs.gps|round(3))}}).
 {% if other_pipelines|length > 0 %}
 This candidate was also found by the {{naturalotherpipelines(other_pipelines)}}
 analysis pipeline{% if other_pipelines|length > 1 %}s{% endif %}.
@@ -17,57 +17,62 @@ analysis pipeline{% if other_pipelines|length > 1 %}s{% endif %}.
 This gravitational wave candidate is not significant enough on its own to produce a public alert
 but its coincidence with the {% if SNEWS %}SNEWS neutrino{% endif %}{% if GRB %}GRB{% endif %} trigger increases its significance.
 {% else %}
-{{gracedb_id}} is an event of interest because its false alarm rate, as
-determined by the online analysis, is {{naturalfar(far)}}.
+{{GraceID}} is an event of interest because its false alarm rate, as
+determined by the online analysis, is {{naturalfar(FAR)}}.
 {% endif %}
 The event's properties can be found at this URL:
 
-{{ gracedb_service_url }}{{ gracedb_id }}
+{{EventPage}}
 
-{% if classifications|length > 0 %}
+{% if Classification|length > 0 %}
 The classification of the signal, in order of descending probability,
-is {{naturalclassifications(classifications)}}.
+is {{naturalclassifications(Classification)}}.
 {% endif %}
 
-{% if prob_has_ns is not none %}
+{% if Properties.HasNS is not none %}
 If the candidate is astrophysical in origin, there
-is {{evidence_for(prob_has_ns)}} the lighter compact object having a mass
-<2.83 solar masses (ProbHasNS: {{probability(prob_has_ns)}}). This 2.83 solar
+is {{evidence_for(Properties.HasNS)}} the lighter compact object having a mass
+<2.83 solar masses (ProbHasNS: {{probability(Properties.HasNS)}}). This 2.83 solar
 mass cutoff corresponds to the maximum neutron star mass assuming a very stiff
 (2H) equation of state. Assuming this neutron star equation of state and the
 masses and spins inferred from the signal, there is
-{{evidence_for(prob_has_remnant)}} matter outside the final compact object
-(ProbHasRemnant: {{probability(prob_has_remnant)}}).
+{{evidence_for(Properties.HasRemnant)}} matter outside the final compact object
+(ProbHasRemnant: {{probability(Properties.HasRemnant)}}).
 {% endif %}
 
-{% if skymaps|length == 0 %}
-No{% else %}{{ skymaps|length|apnumber|capitalize }}{% endif %} skymap{% if skymaps|length == 1 %} is{% else %}s are{% endif %} available at this time{% if skymaps|length > 0 %} and can be retrieved from the GraceDB event page:{% else %}.{% endif %}
+{% if GW_SKYMAP %}
+{# No{% else %}{{ skymaps|length|apnumber|capitalize }}{% endif %} skymap{% if skymaps|length == 1 %} is{% else %}s are{% endif %} available at this time{% if skymaps|length > 0 %} and can be retrieved from the GraceDB event page:{% else %}.{% endif %} #}
+One sky map is available at this time and can be retrieved from the GraceDB event page:
 
-{% for skymap in skymaps %}
- * {{skymap.filename}}, an {% if skymap.alert_type == 'initial' %}initial{% else %}updated{% endif %} localization generated by {{skymap.pipeline}},
-   distributed via GCN notice about {{skymap.latency|naturaldelta}} after the event
-{% endfor %}
+{# {% for skymap in skymaps %} #}
+ * {{GW_SKYMAP.skymap_fits.split('/')[-1]}}, an {% if AlertType == 'Initial' %}initial{% else %}updated{% endif %} localization generated by {{GW_SKYMAP.name}},
+   distributed via GCN notice about {{(Date - DateObs).to('s').value|naturaldelta}} after the event
+{# {% endfor %} #}
 
-{% if skymaps|length != 0 %}
-The preferred skymap at this time is {{preferred_skymap}}.
-For the {{preferred_skymap}} skymap, the {{cl}}% credible region is
-{% if not include_ellipse %}
-{{greedy_area|round|int}} deg2.
+{# {% if skymaps|length != 0 %} #}
+The preferred skymap at this time is {{GW_SKYMAP.skymap_fits.split('/')[-1]}}.
+For the {{GW_SKYMAP.name}} skymap, the {{probability(GW_SKYMAP.credible_level)}} credible region is
+{% if GW_SKYMAP.ellipse.area > 1.35*GW_SKYMAP.area %}
+{{GW_SKYMAP.area|round|int}} deg2.
 {% else %}
-well fit by an ellipse with an area of {{ellipse_area|round|int}} deg2 described
+well fit by an ellipse with an area of {{GW_SKYMAP.ellipse.area|round|int}} deg2 described
 by the following DS9 region (right ascension, declination, semi-major axis,
 semi-minor axis, position angle of the semi-minor axis):
    icrs; ellipse(
-      {{- ra.to('15 arcsec').round().to_string(unit='hourangle', pad=True,
-                                               alwayssign=False) }},
-      {{- ' ' }}{{ dec.to('arcsec').round().to_string(unit='deg', pad=True,
-                                                      alwayssign=True) }},
-      {{- ' ' }}{{ a.round().to_string(fields=1, unit='deg') }},
-      {{- ' ' }}{{ b.round().to_string(fields=1, unit='deg') }},
-      {{- ' ' }}{{ pa.round().to_string(fields=1, unit='deg') }})
+      {{- GW_SKYMAP.ellipse.ra.to('15 arcsec').round().to_string(
+          unit='hourangle', pad=True, alwayssign=False) }},
+      {{- ' ' }}{{ GW_SKYMAP.ellipse.dec.to('arcsec').round().to_string(
+          unit='deg', pad=True, alwayssign=True) }},
+      {{- ' ' }}{{ GW_SKYMAP.ellipse.a.round().to_string(
+          fields=1, unit='deg') }},
+      {{- ' ' }}{{ GW_SKYMAP.ellipse.b.round().to_string(
+          fields=1, unit='deg') }},
+      {{- ' ' }}{{ GW_SKYMAP.ellipse.pa.round().to_string(
+          fields=1, unit='deg') }})
 {% endif %}
-{% if distmu %}
-Marginalized over the whole sky, the a posteriori luminosity distance estimate is {{distmu|round|int}} +/- {{distsig|round|int}} Mpc.
+{% if GW_SKYMAP.distmean is defined %}
+Marginalized over the whole sky, the a posteriori luminosity distance estimate
+is {{GW_SKYMAP.distmean|round|int}} +/- {{GW_SKYMAP.diststd|round|int}} Mpc.
 {% endif %}
 {% endif %}
 {% endfilter %}

ligo/followup_advocate/templates/initial_circular.jinja2

@@ -3,7 +3,6 @@
                             naturalclassifications,
                             citations %}
 {% include 'authors.jinja2' %}
-
 {% filter rewrap %}
 {% include 'initial_body.jinja2' %}
 

ligo/followup_advocate/templates/macros.jinja2

@@ -16,13 +16,13 @@ calculating the odds ratio as K=P/(1-P), and using thresholds based on the
 table of Kass and Raftery (1995, https://doi.org/10.2307%2F2291091).
 
 {% macro evidence_for(value) %}
-{% if 95 <= value %}
+{% if 0.95 <= value %}
 strong evidence for
-{%- elif 75 <= value < 95 %}
+{%- elif 0.75 <= value < 0.95 %}
 evidence for
-{%- elif 25 <= value < 75 %}
+{%- elif 0.25 <= value < 0.75 %}
 indeterminate evidence for
-{%- elif 5 <= value < 25 %}
+{%- elif 0.5 <= value < 0.25 %}
 evidence against
 {%- else %}
 strong evidence against
@@ -33,12 +33,12 @@ strong evidence against
 Macro to render a percentage between 0 and 100, capping at <1% and >99%.
 
 {% macro probability(value) %}
-{% if value > 99 %}
+{% if value > 0.99 %}
 >99
-{%- elif value < 1 %}
+{%- elif value < 0.01 %}
 <1
 {%- else %}
-{{value|round|int}}
+{{(100*value)|round|int}}
 {%- endif %}
 %
 {%- endmacro %}

ligo/followup_advocate/templates/subject.jinja2

 {% from 'macros.jinja2' import renamegroup %}
 
-SUBJECT: LIGO/Virgo {{gracedb_id}}: {{subject}} of a GW {{renamegroup(group)}} candidate
+SUBJECT: LIGO/Virgo {{GraceID}}: Identification of a GW {{renamegroup(Group)}} candidate

ligo/followup_advocate/tool.py

@@ -29,7 +29,9 @@ def main(args=None):
     cmd.add_argument(
         '-m', '--mailto', action='store_true',
         help='Open new message in default e-mail client [default: false]')
-    cmd.add_argument('gracedb_id', metavar='S123456', help='GraceDB ID')
+    cmd.add_argument(
+        'voevent', metavar='VOEVENT.xml', type=argparse.FileType('rb'),
+        default='-', help='VOEvent XML file')
 
     cmd = add_command(followup_advocate.compose_RAVEN, parents=[authors])
     cmd.add_argument('gracedb_id', metavar='S123456', help='GraceDB ID')

ligo/followup_advocate/voevent.py

+from collections import defaultdict
+from functools import partial
+from io import BytesIO
+import os
+from shutil import copyfileobj
+from tempfile import NamedTemporaryFile
+from urllib.error import HTTPError, URLError
+
+from astropy.coordinates import Angle, Latitude, Longitude
+from astropy.time import Time
+from astropy import units as u
+from ligo.gracedb import rest
+from ligo.skymap.io.fits import read_sky_map
+from ligo.skymap.postprocess import find_injection_moc, find_ellipse
+from lxml.etree import iterparse, QName
+
+
+def read_sky_map_from_gracedb(url, *args, **kwargs):
+    """Read a sky map from GraceDb given its URL."""
+    service, apipath, path = url.partition('/api/')
+    service += apipath
+    try:
+        remote = rest.GraceDb(service).get(raw=True)
+        with NamedTemporaryFile(mode='w+b') as local:
+            copyfileobj(remote, local)
+            local.flush()
+            local.seek(0)
+            return read_sky_map(local, *args, **kwargs)
+    except rest.HTTPError as e:
+        # Translate GraceDb errors to Python standard library errors
+        raise HTTPError(None, e.status, e.message, None, None)
+
+
+def read_sky_map_maybe_from_gracedb(url, *args, **kwargs):
+    """Read a sky map from GraceDb or any arbitrary site given its URL.
+    Perform authentication if necessary, and try to download a MOC sky map
+    (e.g. a sky map without the .gz extension) if available."""
+    # Determine file extension.
+    part, ext = os.path.splitext(url)
+
+    # A list of partial function calls to try.
+    attempts = []
+    for func in read_sky_map, read_sky_map_from_gracedb:
+        if ext == '.gz':
+            attempts.append(partial(func, part))
+        attempts.append(partial(func, url))
+
+    # Loop over the attempts. If all of them fail, then raise the last
+    # exception.
+    for attempt in attempts:
+        try:
+            return attempt(*args, **kwargs)
+        except URLError as e:
+            exc = e  # save exception
+    raise exc
+
+
+def load_skymap(params, cl=0.9):
+    # Read the FITS file. If it fails as an ordinary download, then retry
+    # with GraceDb in case authentication is required.
+    skymap = read_sky_map_maybe_from_gracedb(params['skymap_fits'], moc=True)
+
+    # Store FITS header contents.
+    params.update(skymap.meta)
+
+    # Store statics (area, volume).
+    stats = find_injection_moc(skymap, contours=[cl])
+    params['credible_level'] = cl
+    params['area'], = stats.contour_areas
+    params['volume'], = stats.contour_vols
+
+    # Store ellipse.
+    ra, dec, a, b, pa, area = find_ellipse(skymap, 100 * cl)
+    params['ellipse'] = {
+        'ra': Longitude(ra * u.deg),
+        'dec': Latitude(dec * u.deg),
+        'a': Angle(a * u.deg),
+        'b': Angle(b * u.deg),
+        'pa': Angle(pa * u.deg),
+        'area': area
+    }
+
+
+def load(fp, skymaps=True):
+    """
+    Parse a VOEvent file into a simplified hierarchy of Python dictionaries.
+
+    Parameters
+    ----------
+    fp : file-like object
+        The input VOEvent file.
+    skymaps : bool
+        Whether to load and process sky maps too (default: True).
+
+    Returns
+    -------
+    result : dict
+        A dictionary representation of the contents of the VOEvent.
+    """
+
+    # Mapping from Param dataType attribute values to Python types.
+    typemap = defaultdict(
+        lambda: str, {'int': int, 'string': str, 'float': float})
+
+    # Result dictionary and stack to store our current location in it.
+    current = result = {}
+    stack = []
+
+    for action, elem in iterparse(fp, events=('start', 'end')):
+
+        # Ignore XML namespaces.
+        tag = QName(elem.tag).localname
+
+        if tag == 'Group':
+            # Upon encountering a Group element, create a nested dictionary to
+            # store values from enclosed elements.
+            if action == 'start':
+                stack.append(current)
+                current = current.setdefault(elem.attrib['type'], {})
+                try:
+                    current['name'] = elem.attrib['name']
+                except KeyError:
+                    pass
+            elif action == 'end':
+                current = stack.pop()
+        elif action == 'start':
+            # Upon encountering an element with any of the following tags,
+            # store them in the current level in the result dictionary.
+            if tag == 'ISOTime':
+                current['DateObs'] = Time(elem.text)
+            elif tag == 'Param':
+                conversion = typemap[elem.attrib.get('dataType')]
+                name = elem.attrib['name']
+                value = conversion(elem.attrib['value'])
+                current[name] = value
+                if skymaps and name == 'skymap_fits':
+                    load_skymap(current)
+            elif tag == 'VOEvent':
+                current['role'] = elem.attrib['role']
+            elif tag == 'Date':
+                current['Date'] = Time(elem.text)
+
+    # Done!
+    return result
+
+
+def loads(s):
+    """
+    Parse a VOEvent file into a simplified hierarchy of Python dictionaries.
+
+    Parameters
+    ----------
+    s : bytes
+        The input VOEvent contents.
+
+    Returns
+    -------
+    result : dict
+        A dictionary representation of the contents of the VOEvent.
+    """
+    return load(BytesIO(s))