Merge branch 'master' of git+ssh://ligo-vcs.phys.uwm.edu/usr/local/git/gstlal

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

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:])