Skip to content
Snippets Groups Projects
Commit 8e70cb87 authored by ChiWai Chan's avatar ChiWai Chan
Browse files

gstlal_inspiral_calc_snr & svd_bank_snr.py 

	-- allow analysis for multiple instrument
	-- add save analysis file option
parent 017e45e1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
gstlal-inspiral/bin/gstlal_inspiral_calc_snr
+ 107
− 101
View file @ 8e70cb87
......@@ -67,6 +67,7 @@ from optparse import OptionParser, OptionGroup, IndentedHelpFormatter
import os
from gstlal import datasource
from gstlal import inspiral
from gstlal import pipeparts
from gstlal import reference_psd
from gstlal import svd_bank
......@@ -119,7 +120,7 @@ def parse_command_line():
parser.add_option_group(group)
group = OptionGroup(parser, "Template Options", "Choose a template from a SVD bank file / a single SnglInspiral Table")
group.add_option("--svd-bank", metavar = "filename", help = "A LIGO light-weight xml / xml.gz file containing svd bank information (require)." )
group.add_option("--svd-bank", metavar = "filename", help = "A LIGO light-weight xml / xml.gz file containing svd bank information. These can be given as a comma separated list such as H1:file1,H2:file2,L1:file3 to analyze multiple instruments (require)." )
group.add_option("--sub-bank-id", type = "int", help = "Bank id is of the form <int>ID_<int>N where N is the sub bank id. (require).")
group.add_option("--row-number", type = "int", help = "The row number of the template (optional). All the SNRs will be outputed if it is not given.")
group.add_option("--table", metavar = "filename", help = "A LIGO light-weight xml.gz file containing SnglInspiral Table. Expecting one template for each instrument only.")
......@@ -138,19 +139,20 @@ def parse_command_line():
group = OptionGroup(parser, "GraceDb Event Options", "Produce SNR time series for gstlal gracedb event.")
group.add_option("--gid", metavar = "gracedb event id", type = "str", help = "The gracedb event id.")
group.add_option("--observatory", metavar = "OBS", type = "str", help = "See gwdatafind.")
group.add_option("--type", metavar = "frame type", type = "str", help = "See gwdatafind.")
group.add_option("--observatory", metavar = "OBS", type = "str", action = "append", help = "See gwdatafind.")
group.add_option("--type", metavar = "frame type", type = "str", action = "append", help = "See gwdatafind.")
group.add_option("--time-span", metavar = "seconds", type = "int", default = 1000, help = "The time span around the event's trigger time (default = 1000).")
parser.add_option_group(group)
group = OptionGroup(parser, "Output Control Options", "Control SNR output")
group.add_option("--outdir", metavar = "directory", type = "str", help = "Output directory for SNR(s) (require).")
group.add_option("--outdir", metavar = "directory", default = ".", type = "str", help = "Output directory for SNR(s) (default = .).")
group.add_option("--save", action = "store_true", default = False, help = "Save frame cache / svd bank / psd (default = True).")
group.add_option("--mode", metavar = "method", type = "int", default = 0, help = "The method (0 = LLOID / 1 = FIR) that is used to calculate SNR (default = 0).")
group.add_option("--complex", action = "store_true", help = "Choose whether to output the complex snr or not.")
group.add_option("--start", metavar = "seconds", type = "float", help = "Start SNR time series at GPS time '--start' (require).")
group.add_option("--end", metavar = "seconds", type = "float", help = "End SNR time series at GPS time '--end' (require).")
group.add_option("--output-width", metavar = "bits", type = "int", default = 32, help = "The size of the output data, can only be 32 or 64 bits (default = 32 bits).")
group.add_option("--instrument", metavar = "name", help = "The detector from which the --reference-psd and --frame-cache are loaded (require).")
group.add_option("--instrument", metavar = "name", type = "str", action = "append", help = "The detector from which the --reference-psd and --frame-cache are loaded (require).")
parser.add_option_group(group)
parser.add_option("--verbose", action = "store_true", help = "Be verbsoe.")
......@@ -160,42 +162,52 @@ def parse_command_line():
# initialize gw_data_source_info
gw_data_source_info = None
if options.outdir is None:
missing_required_options.append("--outdir")
# Setting up GW data
if options.gid is not None:
# Setting up files for GraceDb event
# Setting up inputs for GraceDb event, this will overwrite svd-bank/template/psd related inputs.
if options.observatory is None or options.type is None:
raise ValueError("When using --gid, --observatory and --type must be provided.")
else:
trigger_time, gps_start_time, gps_end_time, channel_name = svd_bank_snr.gwdata_from_event(options.gid, options.observatory, options.type, time_span = options.time_span, outdir = options.outdir, verbose = options.verbose)
gwdata_metavars = svd_bank_snr.framecache_from_event(options.gid, options.observatory, options.type, time_span = options.time_span, outdir = options.outdir, verbose = options.verbose)
# We can hardcoded here, since we know all the information from gracedb.
# This assume everything on the gracedb are correct and complete which could go wrong in future.
# This assume everything on the gracedb are correct and complete which could go wrong in the future.
# (e.g files are deleted)
options.data_source = "frames"
options.frame_cache = os.path.join(options.outdir, "frame.cache")
options.gps_start_time = gps_start_time
options.gps_end_time = gps_end_time
options.channel_name = [channel_name]
options.gps_start_time = gwdata_metavars["gps_start_time"][0]
options.gps_end_time = gwdata_metavars["gps_end_time"][0]
options.channel_name = gwdata_metavars["channels_name"]
options.instrument = gwdata_metavars["instruments"]
gw_data_source_info = datasource.GWDataSourceInfo(options)
# FIXME: This could change in the future
obs2ifo = {"H": "H1", "L": "L1", "V": "V1"}
options.instrument = obs2ifo[options.observatory]
# FIXME: Adjustable parameters, hardcoded here for simplicity.
trigger_time = min(gwdata_metavars["trigger_times"])
options.start = trigger_time - 5
options.end = trigger_time + 5
psd = svd_bank_snr.psd_from_event(options.gid)
# It must be using LLOID
options.mode = 0
banks_dict, options.sub_bank_id, options.row_number = svd_bank_snr.svd_banks_from_event(options.gid, options.verbose)
bank = banks_dict[options.instrument][options.sub_bank_id]
psds_dict = svd_bank_snr.psd_from_event(options.gid, save = options.save, verbose = options.verbose)
banks_dict, options.sub_bank_id, options.row_number = svd_bank_snr.svd_banks_from_event(options.gid, save = options.save, verbose = options.verbose)
return options, gw_data_source_info, bank, psd
return options, gw_data_source_info, banks_dict, psds_dict
else:
# Setting up GW data for non-GraceDb event
gw_data_source_info = datasource.GWDataSourceInfo(options)
if options.instrument not in gw_data_source_info.channel_dict.keys():
raise ValueError("No such instrument: %s in GWDataSourceInfo: (%s)"% (options.instrument, ", ".join(gw_data_source_info.channel_dict.keys())))
# Setting up how many instruments to analyze
if options.instrument is not None:
# Avoid repeated --instrument
options.instrument = list(set(options.instrument))
else:
raise ValueError("Must provide at least one --instrument.")
# Check if --frame-cache contains the data from --instrument
for instrument in options.instrument:
if instrument not in gw_data_source_info.channel_dict.keys():
raise ValueError("No such instrument: %s in GWDataSourceInfo: (%s)"% (instrument, ", ".join(gw_data_source_info.channel_dict.keys())))
# Check SNRs series output
if options.start is None or options.end is None:
......@@ -210,8 +222,10 @@ def parse_command_line():
# Setting up PSD
if options.reference_psd:
psd = lal.series.read_psd_xmldoc(ligolw_utils.load_url(options.reference_psd, contenthandler = lal.series.PSDContentHandler))
if options.instrument not in set(psd):
raise ValueError("No such instrument: %s in psd: (%s)"% (options.instrument, ", ".join(set(psd))))
# Check if --reference-psd contains the psd from --instrument
for instrument in options.instrument:
if instrument not in set(psd):
raise ValueError("No such instrument: %s in psd: (%s)"% (instrument, ", ".join(set(psd))))
else:
options.track_psd = True
......@@ -223,27 +237,26 @@ def parse_command_line():
missing_required_options.append("--svd-bank")
if options.sub_bank_id is None:
missing_required_options.append("--sub-bank-id")
if options.outdir is None:
missing_required_options.append("--outdir")
if options.instrument is None:
missing_required_options.append("--instrument")
# Raise VauleError is required option(s) is/are missing
if missing_required_options:
raise ValueError("Missing required option(s) %s" % ", ".join(sorted(missing_required_options)))
# Setting up SVD bank
banks = svd_bank.read_banks(options.svd_bank, svd_bank.DefaultContentHandler)
if banks is None:
raise ValueError("svd bank is empty: Invaild --svd-bank %s" % options.svd_bank)
else:
if 0 <= options.sub_bank_id < len(banks) :
bank = banks[options.sub_bank_id]
else:
raise ValueError("Invaild --sub-bank-id %d. Possible id [0-%d)\n" % (options.sub_bank_id, len(banks)))
if not( 0 <= options.row_number < len(bank.sngl_inspiral_table) ):
raise ValueError("No such template: Invaild --row-number %d. Possible range [0-%d)\n" % (options.row_number, len(bank.sngl_inspiral_table)))
return options, gw_data_source_info, bank, psd
svdbanks_dict = inspiral.parse_svdbank_string(options.svd_bank)
# Check if --svd-bank contains svd banks from --instrument
for instrument in options.instrument:
if instrument not in set(svdbanks_dict):
raise ValueError("No SVD Banks for --instrument=%s." % instrument)
banks = dict([(ifo, svd_bank.read_banks(svdbanks, svd_bank.DefaultContentHandler)) for ifo, svdbanks in svdbanks_dict.items()])
# Check if --sub-bank-id and --row-number is valid
for bank in banks.values():
if not (0 <= options.sub_bank_id < len(bank)) :
raise ValueError("Invaild --sub-bank-id %d. Possible id [0-%d)\n" % (options.sub_bank_id, len(bank)))
if options.row_number is not None and not (0 <= options.row_number < len(bank[options.sub_bank_id].sngl_inspiral_table)):
raise ValueError("No such template: Invaild --row-number %d. Possible range [0-%d)\n" % (options.row_number, len(bank[options.sub_bank_id].sngl_inspiral_table)))
return options, gw_data_source_info, banks, psd
# Use Finite Impulse Response
elif options.mode == 1:
......@@ -253,10 +266,6 @@ def parse_command_line():
missing_required_options.append("--table")
if options.approximant is None:
missing_required_options.append("--approximant")
if options.outdir is None:
missing_required_options.append("--outdir")
if options.instrument is None:
missing_required_options.append("--instrument")
if options.template_psd is None:
missing_required_options.append("--template-psd")
# Raise VauleError is required option(s) is/are missing
......@@ -272,7 +281,7 @@ def parse_command_line():
else:
raise ValueError("Invalid mode: %d" % options.mode)
options, gw_data_source_info, template, psd = parse_command_line()
options, gw_data_source_info, template, psds_dict = parse_command_line()
if options.veto_segments_file is not None:
veto_segments = ligolw_segments.segmenttable_get_by_name(ligolw_utils.load_filename(options.veto_segments_file, verbose = options.verbose, contenthandler = LIGOLWContentHandler), options.veto_segments_name).coalesce()
......@@ -286,64 +295,61 @@ else:
#====================================================================================================
if options.mode == 0:
bank = template
num_of_row = bank.bank_fragments[0].mix_matrix.shape[1] / 2
lloid_snr = svd_bank_snr.LLOID_SNR(
gw_data_source_info,
bank,
options.instrument,
options.row_number,
options.start,
options.end,
psd = psd,
psd_fft_length = options.psd_fft_length,
ht_gate_threshold = options.ht_gate_threshold,
veto_segments = veto_segments,
track_psd = options.track_psd,
width = options.output_width,
verbose = options.verbose
)
if options.row_number is None:
for index, snr in enumerate(lloid_snr(COMPLEX = options.complex)):
snr.epoch += bank.sngl_inspiral_table[index].end
snrdict = {options.instrument: [snr]}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR_%d-%d-%d.xml.gz" % (options.instrument, index, int(snr.epoch), int(snr.data.length * snr.deltaT))), verbose = options.verbose)
else:
lloidsnr = lloid_snr(COMPLEX = options.complex)
lloidsnr[0].epoch += bank.sngl_inspiral_table[options.row_number].end
snrdict = {options.instrument: lloidsnr}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR_%d-%d-%d.xml.gz" % (options.instrument, options.row_number, int(lloidsnr[0].epoch), int(lloidsnr[0].data.length * lloidsnr[0].deltaT))), verbose = options.verbose)
#
# uncomment to save all snrs in one single XML
#
#snrdict = {options.instrument : lloid_snr(COMPLEX = options.complex)}
#svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR-%d-%d.xml.gz" % (options.instrument, int(snrdict.[options.instrument][0].epoch), int(snrdict[options.instrument][0].data.length * snrdict[options.instrument][0].deltaT))), verbose = options.verbose)
banks_dict = template
for instrument in options.instrument:
lloid_snr = svd_bank_snr.LLOID_SNR(
gw_data_source_info,
banks_dict[instrument][options.sub_bank_id],
instrument,
options.row_number,
options.start,
options.end,
psd = psds_dict,
psd_fft_length = options.psd_fft_length,
ht_gate_threshold = options.ht_gate_threshold,
veto_segments = veto_segments,
track_psd = options.track_psd,
width = options.output_width,
verbose = options.verbose
)
snrdict = {}
if options.row_number is None:
snrdict.setdefault(instrument, [])
for index, snr in enumerate(lloid_snr(COMPLEX = options.complex)):
snr.epoch += banks_dict[instrument][options.sub_bank_id].sngl_inspiral_table[index].end
snrdict[instrument].append(snr)
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SVD_BANK_SNR_%d-%d-%d.xml.gz" % (instrument, options.sub_bank_id, int(snr.epoch), int(snr.data.length * snr.deltaT))), verbose = options.verbose)
else:
snrdict = {instrument : lloid_snr(COMPLEX = options.complex)}
snrdict[instrument][0].epoch += banks_dict[instrument][options.sub_bank_id].sngl_inspiral_table[options.row_number].end
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SVD_BANK_SNR_%d_%d-%d-%d.xml.gz" % (instrument, options.sub_bank_id, options.row_number, int(snrdict[instrument][0].epoch), int(snrdict[instrument][0].data.length * snrdict[instrument][0].deltaT))), verbose = options.verbose)
elif options.mode == 1:
template_table = template
template, time_offset = svd_bank_snr.FIR_SNR.make_template(template_table, options.template_psd, options.sample_rate, options.approximant, options.instrument, options.f_low, f_high = options.f_high, verbose = options.verbose)
#FIXME: proper handle for latency
fir_snr = svd_bank_snr.FIR_SNR(
gw_data_source_info,
template,
options.instrument,
options.sample_rate,
0,
options.start,
options.end,
psd = psd,
psd_fft_length = options.psd_fft_length,
ht_gate_threshold = options.ht_gate_threshold,
veto_segments = veto_segments,
width = options.output_width,
track_psd = options.track_psd,
verbose = options.verbose
)
firsnr = fir_snr(COMPLEX = options.complex)
firsnr[0].epoch += time_offset
snrdict = {options.instrument : firsnr}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict),os.path.join(options.outdir, "%s-SNR-%d-%d.xml.gz" % (options.instrument, int(snrdict[options.instrument][0].epoch), int(snrdict[options.instrument][0].data.length * snrdict[options.instrument][0].deltaT))), verbose = options.verbose)
for instrument in options.instrument:
template, time_offset = svd_bank_snr.FIR_SNR.make_template(template_table, options.template_psd, options.sample_rate, options.approximant, instrument, options.f_low, f_high = options.f_high, verbose = options.verbose)
#FIXME: proper handle for latency
fir_snr = svd_bank_snr.FIR_SNR(
gw_data_source_info,
template,
instrument,
options.sample_rate,
0,
options.start,
options.end,
psd = psds_dict,
psd_fft_length = options.psd_fft_length,
ht_gate_threshold = options.ht_gate_threshold,
veto_segments = veto_segments,
width = options.output_width,
track_psd = options.track_psd,
verbose = options.verbose
)
firsnr = fir_snr(COMPLEX = options.complex)
snrdict = {instrument : firsnr}
snrdict[instrument][0].epoch += time_offset
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict),os.path.join(options.outdir, "%s-FIR_SNR-%d-%d.xml.gz" % (instrument, int(snrdict[instrument][0].epoch), int(snrdict[instrument][0].data.length * snrdict[instrument][0].deltaT))), verbose = options.verbose)
gstlal-inspiral/python/svd_bank_snr.py
+ 42
− 11
View file @ 8e70cb87
......@@ -5,6 +5,7 @@ A gstlal-based direct matched filter in time domain is also implemented.
import os
import sys
import shutil
import numpy
import threading
......@@ -273,9 +274,9 @@ class FIR_SNR(SNR_Pipeline):
@staticmethod
def make_template(template_table, template_psd, sample_rate, approximant, instrument, f_low, f_high = None, autocorrelation_length = None, verbose = False):
row = [row for row in template_table if row.ifo == instrument]
row = [row for row in template_table]
if len(row) != 1 :
raise ValueError("Expecting only one template for instrument=%s or cannot find template for instrument=%s" % (instrument, instrument))
raise ValueError("Expecting only one template or cannot find any template.")
template_psd = lal.series.read_psd_xmldoc(ligolw_utils.load_url(template_psd, contenthandler = lal.series.PSDContentHandler))
if instrument not in set(template_psd):
......@@ -363,7 +364,7 @@ def read_url(filename, contenthandler = SNRContentHandler, verbose = False):
#
#=============================================================================================
def svd_banks_from_event(gid, verbose = False):
def svd_banks_from_event(gid, outdir = ".", save = True, verbose = False):
gracedb_client = gracedb.GraceDb()
coinc_xmldoc = lvalert_helper.get_coinc_xmldoc(gracedb_client, gid)
eventid_trigger_dict = dict((row.event_id, row) for row in lsctables.SnglInspiralTable.get_table(coinc_xmldoc))
......@@ -377,6 +378,17 @@ def svd_banks_from_event(gid, verbose = False):
sys.stderr.write("Files Not Found! Make sure you are on the LIGO-Caltech Computing Cluster or check if file exist.\nAbortting...\n")
sys.exit()
if save:
try:
for bank_url in bank_urls.values():
outname =os.path.join(outdir, os.path.basename(bank_url))
if verbose:
sys.stderr.write("saving SVD bank file to %s ...\n" % outname)
shutil.copyfile(bank_url, outname)
# FIXME: in python > 2.7, OSError will be raised if destination is not writable.
except IOError as e:
raise e
# Just get one of the template bank from any instrument,
# the templates should all have the same template_id because they are exact-matched.
sub_bank_id = None
......@@ -391,21 +403,35 @@ def svd_banks_from_event(gid, verbose = False):
return banks_dict, sub_bank_id, row_number
def gwdata_from_event(gid, observatory, frame_type, time_span = 1000, outdir = ".", verbose = False):
def framecache_from_event(gid, observatories, frame_types, time_span = 1000, outdir = ".", filename = "frame.cache", verbose = False):
assert time_span >= 1000., "Please use time_span larger or equal to 1000."
obs2ifo = {"H": "H1", "L": "L1", "V": "V1"}
observatories = set(observatories)
frame_types = set(frame_types)
if len(observatories) != len(frame_types):
raise ValueError("Must have as many frame_types as observatories.")
# FIXME: This is not reliable, have a better way to map frame_type to observatory?
obs_type_dict = dict([(obs, frame_type) for obs in observatories for frame_type in frame_types if obs == frame_type[0]])
gracedb_client = gracedb.GraceDb()
coinc_xmldoc = lvalert_helper.get_coinc_xmldoc(gracedb_client, gid)
eventid_trigger_dict = dict((row.ifo, row) for row in lsctables.SnglInspiralTable.get_table(coinc_xmldoc))
channel_names_dict = dict([(row.value.split("=")[0], row.value) for row in lsctables.ProcessParamsTable.get_table(coinc_xmldoc) if row.param == "--channel-name"])
trigger_time = eventid_trigger_dict[obs2ifo[observatory]].end
gps_start_time = int(trigger_time - time_span)
gps_end_time = int(trigger_time + time_span)
gwdata_metavar_headers = ["instruments", "trigger_times", "gps_start_time", "gps_end_time", "channels_name"]
gwdata_metavar_values = []
urls = []
for observatory, frame_type in obs_type_dict.items():
trigger_time = eventid_trigger_dict[obs2ifo[observatory]].end
gps_start_time = int(trigger_time - time_span)
gps_end_time = int(trigger_time + time_span)
gwdata_metavar_values.append((obs2ifo[observatory], trigger_time, gps_start_time, gps_end_time, channel_names_dict[obs2ifo[observatory]]))
urls += gwdatafind.find_urls(observatory, frame_type, gps_start_time, gps_end_time)
urls = gwdatafind.find_urls(observatory, frame_type, gps_start_time, gps_end_time)
with open(os.path.join(outdir, "frame.cache"), "w") as cache:
for url in urls:
filename = str(CacheEntry.from_T050017(url))
......@@ -415,9 +441,14 @@ def gwdata_from_event(gid, observatory, frame_type, time_span = 1000, outdir = "
if verbose:
sys.stderr.write("Done.\n")
return trigger_time, gps_start_time, gps_end_time, channel_names_dict[obs2ifo[observatory]]
return dict(zip(gwdata_metavar_headers, zip(*gwdata_metavar_values)))
def psd_from_event(gid, filename = "psd.xml.gz"):
def psd_from_event(gid, outdir = ".", save = True, filename = "psd.xml.gz", verbose = False):
gracedb_client = gracedb.GraceDb()
psd_fileobj = lvalert_helper.get_filename(gracedb_client, gid, filename)
return lal.series.read_psd_xmldoc(ligolw_utils.load_fileobj(psd_fileobj, contenthandler = lal.series.PSDContentHandler))
xmldoc = ligolw_utils.load_fileobj(psd_fileobj, contenthandler = lal.series.PSDContentHandler)
if save:
if verbose:
sys.stderr.write("saving psd file to %s ...\n" % os.path.join(outdir, filename))
ligolw_utils.write_filename(xmldoc, filename, gz = filename.endswith("gz"))
return lal.series.read_psd_xmldoc(xmldoc)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment