Revision 891994c5246e5e41805b7aaae6539fb79595c86d authored by Alex Nitz on 11 July 2017, 08:22:10 UTC, committed by GitHub on 11 July 2017, 08:22:10 UTC
* Revert "Add option in pycbc_ringinj to taper time-domain ringdown (#1766)" This reverts commit 31db6f8f015dd4c4605172803cc3266b67bcad95. * Revert "Add build for Debian virtual environment (#1756)" This reverts commit d4eb7ee7b69940a54b1a505bead656e98618f34c.
1 parent 31db6f8
live.py
import logging
import os
import pycbc
import numpy
import lal
from pycbc_glue.ligolw import ligolw
from pycbc_glue.ligolw import lsctables
from pycbc_glue.ligolw import utils as ligolw_utils
from pycbc_glue.ligolw.utils import process as ligolw_process
from pycbc_glue.ligolw import param as ligolw_param
from pycbc import version as pycbc_version
from pycbc import pnutils
from pycbc.tmpltbank import return_empty_sngl
from pycbc.filter import compute_followup_snr_series
#FIXME Legacy build PSD xml helpers, delete me when we move away entirely from
# xml formats
def _build_series(series, dim_names, comment, delta_name, delta_unit):
from pycbc_glue.ligolw import array as ligolw_array
Attributes = ligolw.sax.xmlreader.AttributesImpl
elem = ligolw.LIGO_LW(Attributes({u"Name": unicode(series.__class__.__name__)}))
if comment is not None:
elem.appendChild(ligolw.Comment()).pcdata = comment
elem.appendChild(ligolw.Time.from_gps(series.epoch, u"epoch"))
elem.appendChild(ligolw_param.from_pyvalue(u"f0", series.f0, unit=u"s^-1"))
delta = getattr(series, delta_name)
if numpy.iscomplexobj(series.data.data):
data = numpy.row_stack((numpy.arange(len(series.data.data)) * delta,
series.data.data.real, series.data.data.imag))
else:
data = numpy.row_stack((numpy.arange(len(series.data.data)) * delta, series.data.data))
a = ligolw_array.from_array(series.name, data, dim_names=dim_names)
a.Unit = str(series.sampleUnits)
dim0 = a.getElementsByTagName(ligolw.Dim.tagName)[0]
dim0.Unit = delta_unit
dim0.Start = series.f0
dim0.Scale = delta
elem.appendChild(a)
return elem
def snr_series_to_xml(snr_series, document, sngl_inspiral_id):
"""Save an SNR time series into an XML document, in a format compatible
with BAYESTAR.
"""
snr_lal = snr_series.lal()
snr_lal.name = 'snr'
snr_lal.sampleUnits = ''
snr_xml = _build_series(snr_lal, (u'Time', u'Time,Real,Imaginary'), None,
'deltaT', 's')
snr_node = document.childNodes[-1].appendChild(snr_xml)
eid_param = ligolw_param.new_param(u'event_id', u'ilwd:char',
unicode(sngl_inspiral_id))
snr_node.appendChild(eid_param)
def make_psd_xmldoc(psddict):
Attributes = ligolw.sax.xmlreader.AttributesImpl
xmldoc = ligolw.Document()
root_name = u"psd"
lw = xmldoc.appendChild(ligolw.LIGO_LW(Attributes({u"Name": root_name})))
for instrument, psd in psddict.items():
xmlseries = _build_series(psd, (u"Frequency,Real", u"Frequency"),
None, 'deltaF', 's^-1')
fs = lw.appendChild(xmlseries)
fs.appendChild(ligolw_param.from_pyvalue(u"instrument", instrument))
return xmldoc
class SingleCoincForGraceDB(object):
"""Create xml files and submit them to gracedb from PyCBC Live"""
def __init__(self, ifos, coinc_results, **kwargs):
"""Initialize a ligolw xml representation of a zerolag trigger
for upload from pycbc live to gracedb.
Parameters
----------
ifos: list of strs
A list of the ifos pariticipating in this trigger
coinc_results: dict of values
A dictionary of values. The format is defined in
pycbc/events/coinc.py and matches the on disk representation
in the hdf file for this time.
"""
followup_ifos = kwargs.get('followup_ifos') or []
self.template_id = coinc_results['foreground/%s/template_id' % ifos[0]]
# remember if this should be marked as HWINJ
self.is_hardware_injection = ('HWINJ' in coinc_results)
# remember if we want to use a non-standard gracedb server
self.gracedb_server = kwargs.get('gracedb_server')
# compute SNR time series if needed, and figure out which of
# the followup detectors are usable
subthreshold_sngl_time = numpy.mean(
[coinc_results['foreground/%s/end_time' % ifo]
for ifo in ifos])
self.upload_snr_series = kwargs.get('upload_snr_series')
usable_ifos = []
if self.upload_snr_series:
self.snr_series = {}
self.snr_series_psd = {}
htilde = kwargs['bank'][self.template_id]
for ifo in ifos + followup_ifos:
if ifo in ifos:
trig_time = coinc_results['foreground/%s/end_time' % ifo]
else:
trig_time = subthreshold_sngl_time
# NOTE we only check the state/DQ of followup IFOs here.
# IFOs producing the coincidence are assumed to also
# produce valid SNR series.
snr_series, snr_series_psd = compute_followup_snr_series(
kwargs['data_readers'][ifo], htilde, trig_time,
check_state=(ifo in followup_ifos))
if snr_series is not None:
self.snr_series[ifo] = snr_series
self.snr_series_psd[ifo] = snr_series_psd
usable_ifos.append(ifo)
else:
usable_ifos = ifos
# Set up the bare structure of the xml document
outdoc = ligolw.Document()
outdoc.appendChild(ligolw.LIGO_LW())
proc_id = ligolw_process.register_to_xmldoc(
outdoc, 'pycbc', {}, ifos=usable_ifos, comment='',
version=pycbc_version.git_hash,
cvs_repository='pycbc/'+pycbc_version.git_branch,
cvs_entry_time=pycbc_version.date).process_id
# Set up coinc_definer table
coinc_def_table = lsctables.New(lsctables.CoincDefTable)
coinc_def_id = lsctables.CoincDefID(0)
coinc_def_row = lsctables.CoincDef()
coinc_def_row.search = "inspiral"
coinc_def_row.description = "sngl_inspiral<-->sngl_inspiral coincs"
coinc_def_row.coinc_def_id = coinc_def_id
coinc_def_row.search_coinc_type = 0
coinc_def_table.append(coinc_def_row)
outdoc.childNodes[0].appendChild(coinc_def_table)
# Set up coinc inspiral and coinc event tables
coinc_id = lsctables.CoincID(0)
coinc_event_table = lsctables.New(lsctables.CoincTable)
coinc_event_row = lsctables.Coinc()
coinc_event_row.coinc_def_id = coinc_def_id
coinc_event_row.nevents = len(usable_ifos)
coinc_event_row.instruments = ','.join(usable_ifos)
coinc_event_row.time_slide_id = lsctables.TimeSlideID(0)
coinc_event_row.process_id = proc_id
coinc_event_row.coinc_event_id = coinc_id
coinc_event_row.likelihood = 0.
coinc_event_table.append(coinc_event_row)
outdoc.childNodes[0].appendChild(coinc_event_table)
# Set up sngls
sngl_inspiral_table = lsctables.New(lsctables.SnglInspiralTable)
coinc_event_map_table = lsctables.New(lsctables.CoincMapTable)
sngl_populated = None
for sngl_id, ifo in enumerate(usable_ifos):
sngl = return_empty_sngl(nones=True)
sngl.event_id = lsctables.SnglInspiralID(sngl_id)
sngl.process_id = proc_id
sngl.ifo = ifo
names = [n.split('/')[-1] for n in coinc_results
if 'foreground/%s' % ifo in n]
for name in names:
val = coinc_results['foreground/%s/%s' % (ifo, name)]
if name == 'end_time':
sngl.set_end(lal.LIGOTimeGPS(val))
else:
try:
setattr(sngl, name, val)
except AttributeError:
pass
if sngl.mass1 and sngl.mass2:
sngl.mtotal, sngl.eta = pnutils.mass1_mass2_to_mtotal_eta(
sngl.mass1, sngl.mass2)
sngl.mchirp, _ = pnutils.mass1_mass2_to_mchirp_eta(
sngl.mass1, sngl.mass2)
sngl_populated = sngl
if sngl.snr:
sngl.eff_distance = (sngl.sigmasq)**0.5 / sngl.snr
sngl_inspiral_table.append(sngl)
# Set up coinc_map entry
coinc_map_row = lsctables.CoincMap()
coinc_map_row.table_name = 'sngl_inspiral'
coinc_map_row.coinc_event_id = coinc_id
coinc_map_row.event_id = sngl.event_id
coinc_event_map_table.append(coinc_map_row)
if self.upload_snr_series:
snr_series_to_xml(self.snr_series[ifo], outdoc, sngl.event_id)
# for subthreshold detectors, respect BAYESTAR's assumptions and checks
bayestar_check_fields = ('mass1 mass2 mtotal mchirp eta spin1x '
'spin1y spin1z spin2x spin2y spin2z').split()
for sngl in sngl_inspiral_table:
if sngl.ifo in followup_ifos:
for bcf in bayestar_check_fields:
setattr(sngl, bcf, getattr(sngl_populated, bcf))
sngl.set_end(lal.LIGOTimeGPS(subthreshold_sngl_time))
outdoc.childNodes[0].appendChild(coinc_event_map_table)
outdoc.childNodes[0].appendChild(sngl_inspiral_table)
# Set up the coinc inspiral table
coinc_inspiral_table = lsctables.New(lsctables.CoincInspiralTable)
coinc_inspiral_row = lsctables.CoincInspiral()
# This seems to be used as FAP, which should not be in gracedb
coinc_inspiral_row.false_alarm_rate = 0
coinc_inspiral_row.minimum_duration = 0.
coinc_inspiral_row.set_ifos(usable_ifos)
coinc_inspiral_row.coinc_event_id = coinc_id
coinc_inspiral_row.mchirp = sngl_populated.mchirp
coinc_inspiral_row.mass = sngl_populated.mtotal
coinc_inspiral_row.end_time = sngl_populated.end_time
coinc_inspiral_row.end_time_ns = sngl_populated.end_time_ns
coinc_inspiral_row.snr = coinc_results['foreground/stat']
far = 1.0 / (lal.YRJUL_SI * coinc_results['foreground/ifar'])
coinc_inspiral_row.combined_far = far
coinc_inspiral_table.append(coinc_inspiral_row)
outdoc.childNodes[0].appendChild(coinc_inspiral_table)
self.outdoc = outdoc
self.time = sngl_populated.get_end()
def save(self, filename):
"""Write this trigger to gracedb compatible xml format
Parameters
----------
filename: str
Name of file to write to disk.
"""
ligolw_utils.write_filename(self.outdoc, filename)
def upload(self, fname, psds, low_frequency_cutoff,
testing=True,
extra_strings=None,
):
"""Upload this trigger to gracedb
Parameters
----------
fname: str
The name to give the xml file associated with this trigger
pds: dict of pybc.types.FrequencySeries
A ifo keyed dictionary of psds to be uploaded in association
with this trigger.
low_frequency_cutoff: float
The low frequency cutoff of the psds.
testing: bool
Switch to determine if the upload should be sent to gracedb as a
test trigger (True) or a production trigger (False)
"""
from ligo.gracedb.rest import GraceDb
# first of all, make sure the event and PSDs are saved on disk
# as GraceDB operations can fail later
self.save(fname)
psds_lal = {}
for ifo in psds:
psd = psds[ifo]
kmin = int(low_frequency_cutoff / psd.delta_f)
fseries = lal.CreateREAL8FrequencySeries(
"psd", psd.epoch, low_frequency_cutoff, psd.delta_f,
lal.StrainUnit**2 / lal.HertzUnit, len(psd) - kmin)
fseries.data.data = psd.numpy()[kmin:] / pycbc.DYN_RANGE_FAC ** 2.0
psds_lal[ifo] = fseries
psd_xmldoc = make_psd_xmldoc(psds_lal)
psd_xml_path = os.path.splitext(fname)[0] + '-psd.xml.gz'
ligolw_utils.write_filename(psd_xmldoc, psd_xml_path, gz=True)
if self.upload_snr_series:
snr_series_fname = os.path.splitext(fname)[0] + '.hdf'
for ifo in self.snr_series:
self.snr_series[ifo].save(snr_series_fname,
group='%s/snr' % ifo)
self.snr_series_psd[ifo].save(snr_series_fname,
group='%s/psd' % ifo)
# try connecting to GraceDB
try:
gracedb = GraceDb(self.gracedb_server) \
if self.gracedb_server is not None else GraceDb()
except Exception as exc:
logging.error('Cannot connect to GraceDB')
logging.error(str(exc))
logging.error('Carrying on, but event %s will NOT be uploaded!', fname)
return None
# create GraceDB event
group = 'Test' if testing else 'CBC'
try:
r = gracedb.createEvent(group, "pycbc", fname, "AllSky").json()
except Exception as exc:
logging.error('Cannot create GraceDB event')
logging.error(str(exc))
logging.error('Carrying on, but event %s will NOT be uploaded!', fname)
return None
logging.info("Uploaded event %s", r["graceid"])
if self.is_hardware_injection:
try:
gracedb.writeLabel(r['graceid'], 'INJ')
except Exception as exc:
logging.error("Cannot tag event %s as an injection", r["graceid"])
logging.error(str(exc))
logging.info("Tagging event %s as an injection", r["graceid"])
# upload PSDs
try:
gracedb.writeLog(r["graceid"],
"PyCBC PSD estimate from the time of event",
"psd.xml.gz", open(psd_xml_path, "rb").read(),
"psd").json()
except Exception as exc:
logging.error("Cannot upload PSDs for event %s", r["graceid"])
logging.error(str(exc))
logging.info("Uploaded PSDs for event %s", r["graceid"])
# add other tags and comments
try:
gracedb.writeLog(r["graceid"],
"Using PyCBC code hash %s" % pycbc_version.git_hash).json()
extra_strings = [] if extra_strings is None else extra_strings
for text in extra_strings:
gracedb.writeLog(r["graceid"], text).json()
except Exception as exc:
logging.error("Cannot write comments for event %s", r["graceid"])
logging.error(str(exc))
# upload SNR series in HDF format
if self.upload_snr_series:
try:
gracedb.writeFile(r['graceid'], snr_series_fname)
except Exception as exc:
logging.error("Cannot upload HDF SNR series for event %s",
r["graceid"])
logging.error(str(exc))
return r['graceid']
class SingleForGraceDB(SingleCoincForGraceDB):
"""Create xml files and submit them to gracedb from PyCBC Live"""
def __init__(self, ifo, sngls_dict, hardware_injection=False, **kwds):
"""Initialize a ligolw xml representation of this single trigger for
upload to gracedb
Parameters
----------
ifo: str
The IFO that the trigger came from.
sngls_dict: dict
Dictionary of singles parameters. Must include template parameters
and both 'ifar' and 'stat' values.
"""
fake_coinc = {}
fake_coinc['foreground/stat'] = sngls_dict.pop('stat')
fake_coinc['foreground/ifar'] = sngls_dict.pop('ifar')
for key in sngls_dict:
fake_coinc['foreground/%s/%s' % (ifo, key)] = sngls_dict[key]
if hardware_injection:
fake_coinc['HWINJ'] = True
SingleCoincForGraceDB.__init__(self, [ifo], fake_coinc, **kwds)
Computing file changes ...