Commit dd8616ed authored by Thomas Almeida's avatar Thomas Almeida
Browse files

Merge branch 'format-pep8' into 'master'

Format python files with yapf (pep8)

See merge request lscsoft/spiir!6
parents e2b4d78a 7a86dca6
......@@ -22,7 +22,6 @@
### ------
###
###
"""
This program calculates characteristic SNRs given an injection xml and a cache of reference PSDs.
"""
......@@ -34,20 +33,43 @@ import numpy
from optparse import OptionParser
from gstlal.pipemodules import stats
def parse_command_line():
parser = OptionParser(description = __doc__)
parser.add_option("--output", metavar = "filename", help = "Set the output filename.")
parser.add_option("--ifos", metavar = "name", help = "set the ifos to plot pdf and cdf")
parser.add_option("--ncx2-dof", default = 2, type = "int", help = "set degree of freedom for non-central chi-square distribution.")
parser.add_option("--ncx2-mean-factor", default = 0.045, type = "float", help = "set mean factor that will be multiplied to the given coherent SNR value for non-central chi-square distribution.")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose (optional).")
options, tmp = parser.parse_args()
parser = OptionParser(description=__doc__)
parser.add_option("--output",
metavar="filename",
help="Set the output filename.")
parser.add_option("--ifos",
metavar="name",
help="set the ifos to plot pdf and cdf")
parser.add_option(
"--ncx2-dof",
default=2,
type="int",
help="set degree of freedom for non-central chi-square distribution.")
parser.add_option(
"--ncx2-mean-factor",
default=0.045,
type="float",
help=
"set mean factor that will be multiplied to the given coherent SNR value for non-central chi-square distribution."
)
parser.add_option("-v",
"--verbose",
action="store_true",
help="Be verbose (optional).")
options, tmp = parser.parse_args()
if options.output is None or options.ifos is None:
raise ValueError("Must specify --output and --ifos")
if options.output is None or options.ifos is None:
raise ValueError("Must specify --output and --ifos")
return options,
return options,
options, = parse_command_line()
stats.signal_stats_to_xml(options.output, options.ifos, ncx2_dof = options.ncx2_dof, ncx2_mean_factor = options.ncx2_mean_factor, verbose = options.verbose)
stats.signal_stats_to_xml(options.output,
options.ifos,
ncx2_dof=options.ncx2_dof,
ncx2_mean_factor=options.ncx2_mean_factor,
verbose=options.verbose)
......@@ -34,7 +34,6 @@
# + `--downsample`: Choose if you want to downsample IIR bank to multirate(recommended).
# + `--verbose`: Be verbose.
import sys
import scipy
import numpy
......@@ -52,39 +51,102 @@ from gstlal.spiirbank.cbc_template_iir import Bank
import json
parser = OptionParser(description = __doc__)
parser.add_option("--flow", metavar = "Hz", type = "float", default = 15.0, help = "Set the template low-frequency cut-off (default = 15.0).")
parser.add_option("--snr-cut", type = "float", default = 0.998, help = "Set the SNR contribution kept for filter placement (default = 0.998).")
parser.add_option("--sampleRate", metavar = "Hz", type = "float", default = 4096.0, help = "Set the sample rate of the IIR template bank (optional).")
parser.add_option("--negative-latency", metavar = "seconds", type = "int", default = 0, help = "Set the number of seconds of negative latency.")
parser.add_option("--padding", metavar = "pad", type = "float", default = 1.3, help = "Fractional amount to pad time slices.")
parser.add_option("--autocorrelation-length", metavar = "len", type = "int", default = 201, help = "Autocorrelation length for chisq.")
parser.add_option("--reference-psd", metavar = "filename", help = "load the spectrum from this LIGO light-weight XML file (required).")
parser.add_option("--template-bank", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load the template bank (required).")
parser = OptionParser(description=__doc__)
parser.add_option(
"--flow",
metavar="Hz",
type="float",
default=15.0,
help="Set the template low-frequency cut-off (default = 15.0).")
parser.add_option(
"--snr-cut",
type="float",
default=0.998,
help="Set the SNR contribution kept for filter placement (default = 0.998)."
)
parser.add_option(
"--sampleRate",
metavar="Hz",
type="float",
default=4096.0,
help="Set the sample rate of the IIR template bank (optional).")
parser.add_option("--negative-latency",
metavar="seconds",
type="int",
default=0,
help="Set the number of seconds of negative latency.")
parser.add_option("--padding",
metavar="pad",
type="float",
default=1.3,
help="Fractional amount to pad time slices.")
parser.add_option("--autocorrelation-length",
metavar="len",
type="int",
default=201,
help="Autocorrelation length for chisq.")
parser.add_option(
"--reference-psd",
metavar="filename",
help="load the spectrum from this LIGO light-weight XML file (required).")
parser.add_option(
"--template-bank",
metavar="filename",
help=
"Set the name of the LIGO light-weight XML file from which to load the template bank (required)."
)
#FIXME figure out how to encode instrument info
parser.add_option("--instrument", metavar = "ifo", help = "Set the instrument")
parser.add_option("--waveform-domain", default = "TD", help = "Set the domain chosen method, TD or FD (optional).")
parser.add_option("--approximant", default = "SpinTaylorT4", help = "Set the approximant (default = SpinTaylorT4).")
parser.add_option("--output", metavar = "filename", help = "Set the filename in which to save the template bank (required).")
parser.add_option("--downsample", action = "store_true", help = "Choose if you want to downsample IIR bank (recommended)")
parser.add_option("--optimizer-options", default="{}", help = "Pass dictionary of options to the optimizer (parsed using json.loads).")
parser.add_option("--epsilon-options", default="{}", help = "Pass dictionary of options for epsilon selection (parsed using json.loads).")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose (optional).")
parser.add_option("--instrument", metavar="ifo", help="Set the instrument")
parser.add_option("--waveform-domain",
default="TD",
help="Set the domain chosen method, TD or FD (optional).")
parser.add_option("--approximant",
default="SpinTaylorT4",
help="Set the approximant (default = SpinTaylorT4).")
parser.add_option(
"--output",
metavar="filename",
help="Set the filename in which to save the template bank (required).")
parser.add_option(
"--downsample",
action="store_true",
help="Choose if you want to downsample IIR bank (recommended)")
parser.add_option(
"--optimizer-options",
default="{}",
help=
"Pass dictionary of options to the optimizer (parsed using json.loads).")
parser.add_option(
"--epsilon-options",
default="{}",
help=
"Pass dictionary of options for epsilon selection (parsed using json.loads)."
)
parser.add_option("-v",
"--verbose",
action="store_true",
help="Be verbose (optional).")
options, filenames = parser.parse_args()
required_options = ("template_bank", "output")
missing_options = [option for option in required_options if getattr(options, option) is None]
missing_options = [
option for option in required_options if getattr(options, option) is None
]
if missing_options:
raise ValueError, "missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in sorted(missing_options))
raise ValueError, "missing required option(s) %s" % ", ".join(
"--%s" % option.replace("_", "-")
for option in sorted(missing_options))
# read psd file
if options.reference_psd:
ALLpsd = read_psd_xmldoc(utils.load_filename(options.reference_psd, verbose=options.verbose, contenthandler=ligolw.LIGOLWContentHandler))
ALLpsd = read_psd_xmldoc(
utils.load_filename(options.reference_psd,
verbose=options.verbose,
contenthandler=ligolw.LIGOLWContentHandler))
else:
ALLpsd = None
ALLpsd = None
# generate the iir coefficients and write to xml
bank = Bank()
......@@ -94,20 +156,18 @@ print >> sys.stderr, epsilon_options
print >> sys.stderr, optimizer_options
bank.build_from_tmpltbank(
options.template_bank,
sampleRate = options.sampleRate,
negative_latency = options.negative_latency,
all_psd = ALLpsd,
verbose=options.verbose,
padding=options.padding,
flower=options.flow,
snr_cut=options.snr_cut,
downsample = options.downsample,
optimizer_options = optimizer_options,
waveform_domain = options.waveform_domain,
approximant = options.approximant,
autocorrelation_length = options.autocorrelation_length,
**epsilon_options)
options.template_bank,
sampleRate=options.sampleRate,
negative_latency=options.negative_latency,
all_psd=ALLpsd,
verbose=options.verbose,
padding=options.padding,
flower=options.flow,
snr_cut=options.snr_cut,
downsample=options.downsample,
optimizer_options=optimizer_options,
waveform_domain=options.waveform_domain,
approximant=options.approximant,
autocorrelation_length=options.autocorrelation_length,
**epsilon_options)
bank.write_to_xml(options.output)
......@@ -26,7 +26,6 @@
# + `--output` [filename]: Set the filename in which to save the inj xml(required).
# + `--verbose`: Be verbose.
import sys
import scipy
import numpy
......@@ -36,30 +35,54 @@ import pdb
from glue.ligolw import ligolw, lsctables, array, param, utils
from glue import iterutils
class DefaultContentHandler(ligolw.LIGOLWContentHandler):
pass
array.use_in(DefaultContentHandler)
param.use_in(DefaultContentHandler)
lsctables.use_in(DefaultContentHandler)
parser = OptionParser(description = __doc__)
parser.add_option("--injxml", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load the inj bank (required).")
parser.add_option("--extract-method", metavar = "index|mchirp", default = "index", help = "Set the extract method")
parser.add_option("--extract-range", metavar = "gpsstart:gpsend", default = None, help = "Set the extract range")
parser.add_option("--output", metavar = "filename", help = "Set the filename in which to save the inj (required).")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose (optional).")
parser = OptionParser(description=__doc__)
parser.add_option(
"--injxml",
metavar="filename",
help=
"Set the name of the LIGO light-weight XML file from which to load the inj bank (required)."
)
parser.add_option("--extract-method",
metavar="index|mchirp",
default="index",
help="Set the extract method")
parser.add_option("--extract-range",
metavar="gpsstart:gpsend",
default=None,
help="Set the extract range")
parser.add_option("--output",
metavar="filename",
help="Set the filename in which to save the inj (required).")
parser.add_option("-v",
"--verbose",
action="store_true",
help="Be verbose (optional).")
options, filenames = parser.parse_args()
required_options = ("injxml", "output")
missing_options = [option for option in required_options if getattr(options, option) is None]
missing_options = [
option for option in required_options if getattr(options, option) is None
]
if missing_options:
raise ValueError, "missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in sorted(missing_options))
raise ValueError, "missing required option(s) %s" % ", ".join(
"--%s" % option.replace("_", "-")
for option in sorted(missing_options))
# load the sngl_inspiral table
inj_xmldoc = utils.load_filename(options.injxml, contenthandler = DefaultContentHandler, verbose = options.verbose)
inj_xmldoc = utils.load_filename(options.injxml,
contenthandler=DefaultContentHandler,
verbose=options.verbose)
# Get sim inspiral table
sim_inspiral_table = lsctables.SimInspiralTable.get_table(inj_xmldoc)
......@@ -74,20 +97,23 @@ new_sim_table = lsctables.New(lsctables.SimInspiralTable)
if options.extract_method == "gpstime":
# select the entries
if options.extract_range is None:
raise ValueError, "do nothing"
raise ValueError, "do nothing"
else:
# sim_inspiral_table index start from 0, the input index start from 1
extract_start = int(options.extract_range.split(":")[0]) - 1
extract_end = int(options.extract_range.split(":")[1]) - 1
if options.verbose:
print "the entries we select is [%d, %d)" % (extract_start, extract_end)
iterutils.inplace_filter(lambda row: row.h_end_time >= extract_start and row.h_end_time < extract_end, sim_inspiral_table)
print "the entries we select is [%d, %d)" % (extract_start,
extract_end)
iterutils.inplace_filter(
lambda row: row.h_end_time >= extract_start and row.h_end_time <
extract_end, sim_inspiral_table)
for row_id in range(0, len(sim_inspiral_table)):
new_sim_table.append(sim_inspiral_table[row_id])
# append the new sim_inspiral table to the new LIGO_LW
lw.appendChild(new_sim_table)
......@@ -95,5 +121,7 @@ lw.appendChild(new_sim_table)
xmldoc.appendChild(lw)
# Write to file
utils.write_filename(xmldoc, options.output, gz = options.output.endswith('.gz'), verbose = options.verbose)
utils.write_filename(xmldoc,
options.output,
gz=options.output.endswith('.gz'),
verbose=options.verbose)
......@@ -25,7 +25,6 @@
### Takes an injection xml and a cache of reference PSDs and calculated characteristic SNRs. The values
### will be written to the alpha1 and alpha2 columns in the sim_inspiral table.
###
"""
This program calculates characteristic SNRs given an injection xml and a cache of reference PSDs.
"""
......@@ -50,112 +49,142 @@ import lalsimulation
from functools import partial
import multiprocessing
def parse_command_line():
parser = OptionParser(description = __doc__)
parser.add_option("--reference-psd-cache", metavar = "filename", help = "Set the reference psd cache file.")
parser.add_option("--injection-file", metavar = "filename", help = "Set the injection xml file.")
parser.add_option("--output", metavar = "filename", help = "Set the output xml file.")
parser.add_option("--flow", metavar = "value", type = "float", help = "Set the low frequency for waveform generation and SNR integral.")
parser.add_option("--npool", metavar = "value", type = "int", default = 2, help = "Set value for pool. Number of parallel sub-processes")
parser.add_option("--premerger-cut", metavar = "value", type = "int", default = 0, help = "Set template truncation cut off. Defaults to full bandwidth.")
options, filenames = parser.parse_args()
if options.reference_psd_cache is None:
raise ValueError("Must specify --reference-psd-cache")
if options.injection_file is None:
raise ValueError("Must specify --injection-file")
if options.flow is None:
raise ValueError("Must specify --flow")
if options.output is None:
options.output = options.injection_file
return options, filenames
def parse_command_line():
parser = OptionParser(description=__doc__)
parser.add_option("--reference-psd-cache",
metavar="filename",
help="Set the reference psd cache file.")
parser.add_option("--injection-file",
metavar="filename",
help="Set the injection xml file.")
parser.add_option("--output",
metavar="filename",
help="Set the output xml file.")
parser.add_option(
"--flow",
metavar="value",
type="float",
help="Set the low frequency for waveform generation and SNR integral.")
parser.add_option(
"--npool",
metavar="value",
type="int",
default=2,
help="Set value for pool. Number of parallel sub-processes")
parser.add_option(
"--premerger-cut",
metavar="value",
type="int",
default=0,
help="Set template truncation cut off. Defaults to full bandwidth.")
options, filenames = parser.parse_args()
if options.reference_psd_cache is None:
raise ValueError("Must specify --reference-psd-cache")
if options.injection_file is None:
raise ValueError("Must specify --injection-file")
if options.flow is None:
raise ValueError("Must specify --flow")
if options.output is None:
options.output = options.injection_file
return options, filenames
def calc_expected_snr(premerger_cut, inj):
# FIXME: don't hard-code detectors
snr = dict.fromkeys(("H1", "L1", "V1"), 0.0)
injtime = inj.time_geocent
# Determine which PSD files have GPS times covering the injection time.
psds = dict((seg, psd) for seg, psd in allPSDs.items() if injtime in seg)
if len(psds) < 1:
# We know no PSD covers the injection.
return snr
elif len(psds) == 1:
# Only one PSD covers the injection time.
seg, chosenPSD = psds.popitem()
else:
# More than one PSD covers the injection time. Find the
# one whose segment is closest to being centred on the
# injection. Compute segment centre using overflow-safe
# arithmetic.
seg, chosenPSD = min(psds.items(), key = lambda (seg, psd): abs(seg[0] + abs(seg) / 2. - injtime))
# FIXME have better scheme for calculating the needed sample_rate
sample_rate = 16384.0
approximant = lalsimulation.GetApproximantFromString(str(inj.waveform))
if approximant == lalsimulation.NR_hdf5:
LALparams = lal.CreateDict()
lalsimulation.SimInspiralWaveformParamsInsertNumRelData(LALparams, str(inj.numrel_data))
f_min = inj.f_lower
else:
LALparams = None
f_min = options.flow
h_plus, h_cross = lalsimulation.SimInspiralTD(
m1 = inj.mass1*lal.MSUN_SI,
m2 = inj.mass2*lal.MSUN_SI,
S1x = inj.spin1x,
S1y = inj.spin1y,
S1z = inj.spin1z,
S2x = inj.spin2x,
S2y = inj.spin2y,
S2z = inj.spin2z,
distance = inj.distance*1e6*lal.PC_SI,
inclination = inj.inclination,
phiRef = inj.coa_phase,
longAscNodes = 0.0,
eccentricity = 0.0,
meanPerAno = 0.0,
deltaT = 1.0 / sample_rate,
f_min = f_min,
f_ref = 0.0,
LALparams = LALparams,
approximant = approximant
)
h_plus.epoch += injtime
h_cross.epoch += injtime
# Compute strain in each detector. If one detector wasn't on, snr will be set to zero.
for instrument in snr:
if instrument not in chosenPSD:
continue
h = lalsimulation.SimDetectorStrainREAL8TimeSeries(h_plus, h_cross, inj.longitude, inj.latitude, inj.polarization, lalsimulation.DetectorPrefixToLALDetector(instrument))
#print 'Samples: ',len(h.data.data)
#print 'Length: ',len(h.data.data)/sample_rate
#print 'Length_struct: ',h.data.length
new_length = len(h.data.data) - premerger_cut*int(sample_rate)
h_n = lal.CreateREAL8TimeSeries(h.name, h.epoch, h.f0, h.deltaT, h.sampleUnits, new_length)
h_n.data.data = h.data.data[:new_length]
snr[instrument] = lalsimulation.MeasureSNR(h_n, chosenPSD[instrument], options.flow, 0.85 * (sample_rate / 2.))
#print premerger_cut
return snr
# FIXME: don't hard-code detectors
snr = dict.fromkeys(("H1", "L1", "V1"), 0.0)
injtime = inj.time_geocent
# Determine which PSD files have GPS times covering the injection time.
psds = dict((seg, psd) for seg, psd in allPSDs.items() if injtime in seg)
if len(psds) < 1:
# We know no PSD covers the injection.
return snr
elif len(psds) == 1:
# Only one PSD covers the injection time.
seg, chosenPSD = psds.popitem()
else:
# More than one PSD covers the injection time. Find the
# one whose segment is closest to being centred on the
# injection. Compute segment centre using overflow-safe
# arithmetic.
seg, chosenPSD = min(psds.items(),
key=lambda
(seg, psd): abs(seg[0] + abs(seg) / 2. - injtime))
# FIXME have better scheme for calculating the needed sample_rate
sample_rate = 16384.0
approximant = lalsimulation.GetApproximantFromString(str(inj.waveform))
if approximant == lalsimulation.NR_hdf5:
LALparams = lal.CreateDict()
lalsimulation.SimInspiralWaveformParamsInsertNumRelData(
LALparams, str(inj.numrel_data))
f_min = inj.f_lower
else:
LALparams = None
f_min = options.flow
h_plus, h_cross = lalsimulation.SimInspiralTD(m1=inj.mass1 * lal.MSUN_SI,
m2=inj.mass2 * lal.MSUN_SI,
S1x=inj.spin1x,
S1y=inj.spin1y,
S1z=inj.spin1z,
S2x=inj.spin2x,
S2y=inj.spin2y,
S2z=inj.spin2z,
distance=inj.distance * 1e6 *
lal.PC_SI,
inclination=inj.inclination,
phiRef=inj.coa_phase,
longAscNodes=0.0,
eccentricity=0.0,
meanPerAno=0.0,
deltaT=1.0 / sample_rate,
f_min=f_min,
f_ref=0.0,
LALparams=LALparams,
approximant=approximant)
h_plus.epoch += injtime
h_cross.epoch += injtime
# Compute strain in each detector. If one detector wasn't on, snr will be set to zero.
for instrument in snr:
if instrument not in chosenPSD:
continue
h = lalsimulation.SimDetectorStrainREAL8TimeSeries(
h_plus, h_cross, inj.longitude, inj.latitude, inj.polarization,
lalsimulation.DetectorPrefixToLALDetector(instrument))
#print 'Samples: ',len(h.data.data)
#print 'Length: ',len(h.data.data)/sample_rate
#print 'Length_struct: ',h.data.length
new_length = len(h.data.data) - premerger_cut * int(sample_rate)
h_n = lal.CreateREAL8TimeSeries(h.name, h.epoch, h.f0, h.deltaT,
h.sampleUnits, new_length)
h_n.data.data = h.data.data[:new_length]
snr[instrument] = lalsimulation.MeasureSNR(h_n, chosenPSD[instrument],
options.flow,
0.85 * (sample_rate / 2.))
#print premerger_cut
return snr
@lsctables.use_in