gstlal_feature_extractor_template_overlap + idq_utils.py: fixed bug in...

gstlal_feature_extractor_template_overlap + idq_utils.py: fixed bug in template generation causing missing templates, changed way templates are produced as a result, propagated changes to template validation code

gstlal-ugly/bin/gstlal_feature_extractor_template_overlap

@@ -233,7 +233,6 @@ if __name__ == '__main__':
 	aggregator.makedir(options.output_dir)
 
     # common parameters we will use throughout
-	#max_samp_rate = 2048
 	max_samp_rate = 8192
 	min_samp_rate = 32
 	n_rates = int(numpy.log2(max_samp_rate/min_samp_rate) + 1)
@@ -243,22 +242,20 @@ if __name__ == '__main__':
 
 	# generate templates for each rate considered
 	rates = [min_samp_rate*2**i for i in range(n_rates)]
-	waveforms = {}
-	templates = {}
-	basis_params = {}
-	for rate in rates:
-		flow = rate/4.*0.8
-		fhigh = rate/2.*0.8
-		qhigh = options.qhigh
-		qlow = 3.3166
-		if options.waveform == 'sine_gaussian':
-			waveforms[rate] = idq_utils.SineGaussianGenerator((flow, fhigh), (qlow, qhigh), rate, mismatch = options.mismatch)
-		elif options.waveform == 'half_sine_gaussian':
-			waveforms[rate] = idq_utils.HalfSineGaussianGenerator((flow, fhigh), (qlow, qhigh), rate, mismatch = options.mismatch)
-		else:
-			raise NotImplementedError
-		basis_params[rate] = waveforms[rate].parameter_grid
-		templates[rate] = [waveform for waveform in waveforms[rate].generate_templates(quadrature=False)]
+	downsample_factor = 0.8
+	qhigh = options.qhigh
+	qlow = 3.3166
+	fhigh = max_samp_rate / 2.
+	flow = min_samp_rate / 4.
+
+	if options.waveform == 'sine_gaussian':
+		waveforms = idq_utils.SineGaussianGenerator((flow, fhigh), (qlow, qhigh), rates, mismatch = options.mismatch, downsample_factor=downsample_factor)
+	elif options.waveform == 'half_sine_gaussian':
+		waveforms = idq_utils.HalfSineGaussianGenerator((flow, fhigh), (qlow, qhigh), rates, mismatch = options.mismatch, downsample_factor=downsample_factor)
+	else:
+		raise NotImplementedError
+	basis_params = waveforms.parameter_grid
+	templates = {rate: [waveform for waveform in waveforms.generate_templates(rate, quadrature=False)] for rate in rates}
 
 	# get all templates and params
 	all_templates = list(itertools.chain(*templates.values()))
@@ -267,11 +264,16 @@ if __name__ == '__main__':
 	if options.verbose:
 		print >>sys.stderr, "Creating template overlaps..."
 
+	# zero pad templates to make them the same length
+	max_sample_pts = max(len(template) for template in all_templates)
+	all_templates = [numpy.pad(template, ((max_sample_pts - len(template)) // 2, (max_sample_pts - len(template)) // 2), 'constant') for template in all_templates]
+
 	# calculate overlap for each template and find maximum
 	overlaps = []
 	for this_template in all_templates:
 		overlaps.append(max([numpy.dot(this_template,template) for template in all_templates if not numpy.array_equal(template, this_template)]))
 
+	print >>sys.stderr, "total number of templates: %d" % len(all_templates)
 	print >>sys.stderr, "min overlap specified: %f" % (1 - options.mismatch)
 	print >>sys.stderr, "max template overlap: %f" % max(overlaps)
 	print >>sys.stderr, "min template overlap: %f" % min(overlaps)
@@ -291,13 +293,15 @@ if __name__ == '__main__':
 		print >>sys.stderr, "Creating waveform plots..."
 
 	for rate in rates:
-		for template_id, template in enumerate(random.sample(templates[rate], num_samples)):
+		#for template_id, template in enumerate(random.sample(templates[rate], num_samples)):
+		for template_id in random.sample(numpy.arange(len(templates[rate])), num_samples):
 			waveform_params = ["%s: %.3f" % (name, param) for param, name in zip(basis_params[rate][template_id], param_names)]
+			template = templates[rate][template_id]
 
 			if options.verbose:
 				print >>sys.stderr, "\tCreating waveform plot with parameters: %s" % repr(waveform_params)
 
-			series_fig = plot_waveform(waveforms[rate].times, template, waveform_type, waveform_params)
+			series_fig = plot_waveform(waveforms.times[rate], template, waveform_type, waveform_params)
 			fname = 'plot_%s_%s-timeseries.png' % (str(rate).zfill(4), str(template_id).zfill(4))
 			plot_paths.append((template_id*int(rate),fname))
 			series_fig.savefig(os.path.join(options.output_dir, fname))

gstlal-ugly/python/idq_utils.py

@@ -25,6 +25,8 @@
 #################### 
 
 
+import bisect
+from collections import defaultdict
 import glob
 import logging
 import os
@@ -296,32 +298,59 @@ class HalfSineGaussianGenerator(object):
 	"""
 	Generates half sine gaussian templates based on a f, Q range and a sampling frequency.
 	"""
-	def __init__(self, f_range, q_range, f_samp, mismatch=0.2, tolerance=5e-3):
+	def __init__(self, f_range, q_range, rates, mismatch=0.2, tolerance=5e-3, downsample_factor=0.8):
+
 		### define parameter range
-		self.f_low, self.f_high = f_range
+		self.f_low, self.f_high = (downsample_factor * f_range[0], downsample_factor * f_range[1])
 		self.q_low, self.q_high = q_range
-		self.f_samp = f_samp
+
+		### define grid spacing and edge rolloff for templates
 		self.mismatch = mismatch
 		self.tol = tolerance
 
+		### determine how to scale down templates considered
+		### in frequency band to deal with low pass rolloff
+		self.downsample_factor = downsample_factor
+
+		### determine frequency bands considered
+		self.rates = sorted(set(rates))
+
+		### determine lower frequency limits for rate conversion based on downsampling factor
+		self.freq_breakpoints = [self.downsample_factor * rate / 2. for rate in self.rates[:-1]]
+
 		### define grid and template duration
-		self.parameter_grid = [(f, q, self.duration(f, q)) for f, q in self.generate_f_q_grid(self.f_low, self.f_high, self.q_low, self.q_high)]
-		self.max_duration = max(duration for f, q, duration in self.parameter_grid)
+		self.parameter_grid = defaultdict(list)
+		for rate, f, q in self.generate_f_q_grid(self.f_low, self.f_high, self.q_low, self.q_high):
+			self.parameter_grid[rate].append((f, q, self.duration(f, q)))
 
+		self.max_duration = {rate: max(duration for f, q, duration in parameters) for rate, parameters in self.parameter_grid.items()}
+
+		### determine rest of waveform properties
 		self.phases = [0., numpy.pi/2.]
-		self.times = numpy.linspace(-self.max_duration, 0, int(numpy.ceil(self.max_duration*self.f_samp)), endpoint=True)
-		self.latency = 0
+		self.times = {rate: numpy.linspace(-self.max_duration[rate], 0, self.round_to_next_odd(self.max_duration[rate]*rate), endpoint=True) for rate in self.rates}
+		self.latency = {rate: 0 for rate in self.rates}
+
+	def frequency_to_rate(self, frequency):
+		"""
+		Maps a frequency to the correct sampling rate considered.
+		"""
+		idx = bisect.bisect(self.freq_breakpoints, frequency)
+		return self.rates[idx]
+
+	def round_to_next_odd(self, n):
+			return int(numpy.ceil(n) // 2 * 2 + 1)
 
-	def generate_templates(self, quadrature = True):
+	def generate_templates(self, rate, quadrature = True):
 		"""
 		generate all half sine gaussian templates corresponding to a parameter range and template duration
+		for a given sampling rate.
 		"""
-		for f, q in self.generate_f_q_grid(self.f_low, self.f_high, self.q_low, self.q_high):
+		for f, q, _ in self.parameter_grid[rate]:
 			if quadrature:
 				for phase in self.phases:
 					yield self.waveform(f, q, phase)
 			else:
-				yield self.waveform(f, q, self.phases[0])
+				yield self.waveform(f, q, self.phases[0], rate)
 
 	def duration(self, f, q):
 		"""
@@ -329,18 +358,16 @@ class HalfSineGaussianGenerator(object):
 		"""
 		return 0.5 * (q/(2.*numpy.pi*f)) * numpy.log(1./self.tol)
 	
-	def waveform(self, f, q, phase):
+	def waveform(self, f, q, phase, rate):
 		"""
 		construct half sine gaussian waveforms that taper to tolerance at edges of window
 		f is the central frequency of the waveform
 		"""
-		dt = self.times[1] - self.times[0]
-		rate = 1./dt
 		assert f < rate/2. 
 	
 		# phi is the central frequency of the sine gaussian
 		tau = q/(2.*numpy.pi*f)
-		template = numpy.cos(2.*numpy.pi*f*self.times + phase)*numpy.exp(-1.*self.times**2./tau**2.)
+		template = numpy.cos(2.*numpy.pi*f*self.times[rate] + phase)*numpy.exp(-1.*self.times[rate]**2./tau**2.)
 	
 		# normalize sine gaussians to have unit length in their vector space
 		inner_product = numpy.sum(template*template)
@@ -369,23 +396,26 @@ class HalfSineGaussianGenerator(object):
 	
 	def generate_f_q_grid(self, f_min, f_max, q_min, q_max):
 		"""
-		Generates (f, Q) pairs based on f, Q ranges
+		Generates (f_samp, f, Q) pairs based on f, Q ranges
 		"""
 		for q in self.generate_q_values(q_min, q_max):
 			num_f = self.num_f_templates(f_min, f_max, q)
 			for l in range(num_f):
 				f = f_min * (f_max/f_min)**( (0.5+l) /num_f)
-				if f < f_max / (1 + (numpy.sqrt(11)/q)):
-					yield (f, q)
+				rate = self.frequency_to_rate(f)
+				if f < (rate/2) / (1 + (numpy.sqrt(11)/q)):
+					yield rate, f, q
+				elif rate != max(self.rates):
+					yield (2 * rate), f, q
 
 class SineGaussianGenerator(HalfSineGaussianGenerator):
 	"""
 	Generates sine gaussian templates based on a f, Q range and a sampling frequency.
 	"""
-	def __init__(self, f_range, q_range, f_samp, mismatch=0.2, tolerance=5e-3):
-		super(SineGaussianGenerator, self).__init__(f_range, q_range, f_samp, mismatch=mismatch, tolerance=tolerance)
-		self.times = numpy.linspace(-self.max_duration/2., self.max_duration/2., int(numpy.ceil(self.max_duration*self.f_samp)), endpoint=True)
-		self.latency = self.max_duration / 2.
+	def __init__(self, f_range, q_range, rates, mismatch=0.2, tolerance=5e-3, downsample_factor=0.8):
+		super(SineGaussianGenerator, self).__init__(f_range, q_range, rates, mismatch=mismatch, tolerance=tolerance, downsample_factor=0.8)
+		self.times = {rate: numpy.linspace(-self.max_duration[rate]/2., self.max_duration[rate]/2., self.round_to_next_odd(self.max_duration[rate]*rate), endpoint=True) for rate in self.rates}
+		self.latency = {rate: self.max_duration[rate] / 2. for rate in self.rates}
 
 	def duration(self, f, q):
 		"""