diff --git a/src/bin/gstlal_inspiral b/src/bin/gstlal_inspiral
index d49339b07407558d0d38fd7bdc43a0a89fdfc007..f90d37d04698fecf84c3ed49461a7506d87e898b 100755
--- a/src/bin/gstlal_inspiral
+++ b/src/bin/gstlal_inspiral
@@ -10,7 +10,7 @@
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+# MERCHANTABLITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
@@ -57,7 +57,7 @@ from gstlal import pipeio
from gstlal import pipeparts
from gstlal import cbc_template_fir
from gstlal import misc as gstlalmisc
-
+from gstlal import templates
#
# =============================================================================
@@ -152,6 +152,12 @@ def add_cbc_metadata(xmldoc, process, seg):
return search_summary
+
+
+
+
+
+
#
# =============================================================================
#
@@ -195,7 +201,7 @@ class Bank(object):
self.filter_length = max(fragment.end for fragment in bank_fragments)
self.snr_threshold = snr_threshold
self.logname = logname
-
+
sample_rate = max(bank_fragment.rate for bank_fragment in bank_fragments)
template_bank, self.autocorrelation_bank = cbc_template_fir.generate_templates(lsctables.table.get_table(utils.load_filename(template_bank_filename, gz = (template_bank_filename or "stdin").endswith(".gz"), verbose = verbose), lsctables.SnglInspiralTable.tableName), psd, flow, sample_rate, self.filter_length, autocorrelation_length = autocorrelation_length, verbose = verbose)
@@ -350,6 +356,7 @@ def mkLLOIDsrc(pipeline, instrument, detector, rates, psd = None, psd_fft_length
128: 1.03 / 0.93,
256: 1.03 / 0.93,
512: 1.03 / 0.93,
+ 1024: 1.03 / 1.03,
2048: 1.03 / 1.03
}
@@ -804,12 +811,9 @@ banks = [
Bank(
template_bank_filename,
psd,
- [
- Bank.BankFragment(2048, 0.0, 1.0),
- Bank.BankFragment(512, 1.0, 5.0),
- Bank.BankFragment(256, 5.0, 13.0),
- Bank.BankFragment(128, 13.0, 29.0)
- ],
+ [ Bank.BankFragment(rate,begin,end)
+ for rate,begin,end in
+ templates.time_frequency_boundaries(template_bank_filename,flow=options.flow)],
gate_fap = options.ortho_gate_fap,
snr_threshold = options.snr_threshold,
tolerance = options.svd_tolerance,
@@ -850,6 +854,7 @@ sngl_inspiral_table.sync_next_id()
#
+
pipeline = gst.Pipeline("gstlal_inspiral")
mainloop = gobject.MainLoop()
@@ -881,7 +886,6 @@ handler = LLOIDHandler(mainloop, pipeline, verbose = options.verbose)
# process requested segment
#
-
pipeline.set_state(gst.STATE_PAUSED)
pipeline.seek(1.0, gst.Format(gst.FORMAT_TIME), gst.SEEK_FLAG_FLUSH, gst.SEEK_TYPE_SET, options.seg[0].ns(), gst.SEEK_TYPE_SET, options.seg[1].ns())
pipeline.set_state(gst.STATE_PLAYING)
@@ -893,6 +897,7 @@ mainloop.run()
#
+
sngl_inspiral_table.sort(lambda a, b: cmp(a.end_time, b.end_time) or cmp(a.end_time_ns, b.end_time_ns) or cmp(a.ifo, b.ifo))
for row in sngl_inspiral_table:
row.process_id = process.process_id
diff --git a/src/python/templates.py b/src/python/templates.py
index f83f72bb8f20913379159bef2cf7134ed6c6b453..f3b880b92f227ac28ce684d80f4a4d186638edd3 100644
--- a/src/python/templates.py
+++ b/src/python/templates.py
@@ -29,7 +29,9 @@ import numpy
from pylal import datatypes as laltypes
from pylal import lalfft
-
+from pylal import spawaveform
+from glue.ligolw import lsctables
+from glue.ligolw import utils
__author__ = "Kipp Cannon <kipp.cannon@ligo.org>, Chad Hanna <chad.hanna@ligo.org>, Drew Keppel <drew.keppel@ligo.org>"
__version__ = "FIXME"
@@ -160,3 +162,111 @@ def normalized_autocorrelation(fseries, revplan):
data = tseries.data
tseries.data = data / data[0]
return tseries
+
+
+def time_frequency_boundaries( template_bank_filename,
+ segment_samples_max = 2048,
+ flow = 64,
+ sample_rate_max = 2048,
+ padding = 0.9 ):
+ """
+ The function time_frequency_boundaries splits a template bank up by times
+ for which different sampling rates are appropriate. The function returns
+ a list of 3-tuples of the form (rate,begin,end) where rate is the sampling
+ rate and begin/end mark the boundaries during which the given rate is
+ guaranteed to be appropriate (no template exceeds a frequency of padding*Nyquist
+ during these times and no lower sampling rate would work). For computational
+ purposes, no time interval exceeds max_samples_per_segment. The same rate may
+ therefore apply to more than one segment.
+ """
+ # Round a number up to the nearest power of 2
+ def ceil_pow_2( number ):
+ return 2**(numpy.ceil(numpy.log2( number )))
+
+ # We only allow sample rates that are powers of two.
+ #
+ # h(t) is sampled at 16384Hz, which sets the upper limit
+ # and advligo will likely not reach below 10Hz, which
+ # sets the lower limit (32Hz = ceil_pow_2(2*10)Hz )
+ #
+ allowed_rates = [32,64,128,256,512,1024,2048,4096,8192,16384]
+
+ # Load template bank mass params
+ template_bank_table = (
+ lsctables.table.get_table(
+ utils.load_filename(
+ template_bank_filename,
+ gz=template_bank_filename.endswith("gz") ),
+ "sngl_inspiral") )
+ mass1 = template_bank_table.get_column('mass1')
+ mass2 = template_bank_table.get_column('mass2')
+
+ #
+ # Adjust the allowed_rates to fit with the template bank
+ #
+ ffinal_max = max(spawaveform.ffinal(m1,m2,'schwarz_isco') for m1,m2 in zip(mass1,mass2) )
+ ffinal_max = min(allowed_rates[-1],ffinal_max)
+
+ # FIXME: output sample_rate_max may be bigger than input
+ sample_rate_max = min(ceil_pow_2( 2*(1./padding)* ffinal_max ),ceil_pow_2(sample_rate_max))
+ sample_rate_min = ceil_pow_2( 2*(1./padding)* flow )
+
+
+ allowed_rates = allowed_rates[allowed_rates.index(sample_rate_min):allowed_rates.index(sample_rate_max)] # excludes sample_rate_max
+
+ #
+ # Split up templates by time
+ #
+
+ # FIXME: what happens if padding*sample_rate_min/2 == flow?
+ longest_chirp = max(spawaveform.chirptime(m1,m2,7,flow,sample_rate_max/2) for m1,m2 in zip(mass1,mass2))
+ time_partition = [ longest_chirp ]
+ for rate in allowed_rates:
+ longest_chirp = max(spawaveform.chirptime(m1,m2,7,padding*rate/2,sample_rate_max/2) for m1,m2 in zip(mass1,mass2))
+ time_partition.append( longest_chirp )
+
+ time_partition.append(0)
+ allowed_rates.append(sample_rate_max)
+
+ time_partition.reverse()
+ allowed_rates.reverse()
+
+ #
+ # SVD needs more sample points than templates.
+ # Reject some sampling_rates if the templates don't
+ # spend enough time there
+ # e.g.
+ # allowed_rates = [64,128,256,512,1024,2048]
+ # time_partition = [100,25,10,4,1,0.01,0.001]
+
+
+ num_templates = len(mass1)
+ # FIXME: there are undoubtedly some dangerous corner cases here
+ if not 2*num_templates < segment_samples_max:
+ # FIXME this is bad
+ return (None,None,None)
+
+ for rate,begin,end in zip(allowed_rates,time_partition[:-1],time_partition[1:]):
+ segment_samples = (end-begin)*rate
+ if segment_samples < 2*num_templates:
+ allowed_rates.remove(rate/2) #FIXME: danger at right edge of list
+ time_partition.remove(end)
+
+ # Check that no time segment is excessively long.
+ # If it is, chop it up into the minimal number of
+ # equal-sized bits that satisfy the smallness
+ # requirement
+ for rate,begin,end in zip(allowed_rates,time_partition[:-1],time_partition[1:]):
+ segment_samples = (end-begin)*rate
+ if segment_samples > segment_samples_max:
+ num_segments = numpy.ceil(float(segment_samples)/segment_samples_max)
+ increment = float(end-begin)/num_segments
+ start_index = time_partition.index(begin)
+ index = 1
+ while index < num_segments:
+ allowed_rates.insert(start_index,rate)
+ time_partition.insert(start_index+index,begin+index*increment)
+ index += 1
+
+
+ return zip(allowed_rates,time_partition[:-1],time_partition[1:])