gstlal-inspiral: port to lal swig time and frequency series

gstlal-inspiral/python/cbc_template_fir.py

@@ -59,8 +59,6 @@ from lal import LIGOTimeGPS
 import lalsimulation as lalsim
 
 
-from pylal import datatypes as laltypes
-from pylal import lalfft
 from pylal import spawaveform
 
 
@@ -135,14 +133,16 @@ def generate_template(template_bank_row, approximant, sample_rate, duration, f_l
 	else:
 		raise ValueError("Unsupported approximant given %s" % approximant)
 
-	return laltypes.COMPLEX16FrequencySeries(
+	fseries = lal.CreateCOMPLEX16FrequencySeries(
 		name = "template",
 		epoch = LIGOTimeGPS(hplus.epoch.gpsSeconds, hplus.epoch.gpsNanoSeconds),
 		f0 = 0.0,
 		deltaF = 1.0 / duration,
-		sampleUnits = laltypes.LALUnit("strain"),
-		data = z
+		sampleUnits = lal.Unit("strain"),
+		length = len(z)
 	)
+	fseries.data.data = z
+	return fseries
 
 def condition_imr_template(approximant, data, epoch_time, sample_rate_max, max_ringtime):
 	assert -len(data) / sample_rate_max <= epoch_time < 0.0, "Epoch returned follows a different convention"
@@ -288,17 +288,17 @@ def generate_templates(template_table, approximant, psd, f_low, time_slices, aut
 	if psd is not None:
 		psd = condition_psd(psd, 1.0 / working_duration, minfs = (working_f_low, f_low), maxfs = (sample_rate_max / 2.0 * 0.90, sample_rate_max / 2.0))
 
-	revplan = lalfft.XLALCreateReverseCOMPLEX16FFTPlan(working_length, 1)
-	fwdplan = lalfft.XLALCreateForwardREAL8FFTPlan(working_length, 1)
-	tseries = laltypes.COMPLEX16TimeSeries(
-		data = numpy.zeros((working_length,), dtype = "cdouble")
+	revplan = lal.CreateReverseCOMPLEX16FFTPlan(working_length, 1)
+	fwdplan = lal.CreateForwardREAL8FFTPlan(working_length, 1)
+	tseries = lal.CreateCOMPLEX16TimeSeries(
+		length = working_length
 	)
-	fworkspace = laltypes.COMPLEX16FrequencySeries(
+	fworkspace = lal.CreateCOMPLEX16FrequencySeries(
 		name = "template",
 		epoch = LIGOTimeGPS(0),
 		f0 = 0.0,
 		deltaF = 1.0 / working_duration,
-		data = numpy.zeros((working_length//2 + 1,), dtype = "cdouble")
+		length = working_length // 2 + 1
 	)
 
 	# Check parity of autocorrelation length
@@ -337,7 +337,7 @@ def generate_templates(template_table, approximant, psd, f_low, time_slices, aut
 		#
 
 		if psd is not None:
-			lalfft.XLALWhitenCOMPLEX16FrequencySeries(fseries, psd)
+			lal.WhitenCOMPLEX16FrequencySeries(fseries, psd)
 		fseries = templates.QuadradurePhase.add_quadrature_phase(fseries, working_length)
 
 		#
@@ -352,7 +352,7 @@ def generate_templates(template_table, approximant, psd, f_low, time_slices, aut
 		# transform template to time domain
 		#
 
-		lalfft.XLALCOMPLEX16FreqTimeFFT(tseries, fseries, revplan)
+		lal.COMPLEX16FreqTimeFFT(tseries, fseries, revplan)
 
 		data = tseries.data
 		epoch_time = fseries.epoch.seconds + fseries.epoch.nanoseconds*1.e-9

gstlal-inspiral/python/inspiral.py

@@ -100,8 +100,6 @@ from glue.ligolw.utils import time_slide as ligolw_time_slide
 import lal
 from lal import LIGOTimeGPS
 from pylal import rate
-from pylal.datatypes import LALUnit
-from pylal.datatypes import REAL8FrequencySeries
 
 from gstlal import bottle
 from gstlal import reference_psd
@@ -869,14 +867,16 @@ class Data(object):
 				psddict = {}
 				for instrument in self.seglistdicts["triggersegments"]:
 					elem = self.pipeline.get_by_name("lal_whiten_%s" % instrument)
-					psddict[instrument] = REAL8FrequencySeries(
+					data = numpy.array(elem.get_property("mean-psd"))
+					psddict[instrument] = lal.CreateREAL8FrequencySeries(
 						name = "PSD",
 						epoch = LIGOTimeGPS(lal.UTCToGPS(time.gmtime()), 0),
 						f0 = 0.0,
 						deltaF = elem.get_property("delta-f"),
-						sampleUnits = LALUnit("s strain^2"),	# FIXME:  don't hard-code this
-						data = numpy.array(elem.get_property("mean-psd"))
+						sampleUnits = lal.Unit("s strain^2"),	# FIXME:  don't hard-code this
+						length = len(data)
 					)
+					psddict[instrument].data.data = data
 				fobj = StringIO.StringIO()
 				reference_psd.write_psd_fileobj(fobj, psddict, gz = True, trap_signals = None)
 				common_messages.append(("strain spectral densities", "psd.xml.gz", "psd", fobj.getvalue()))

gstlal-inspiral/python/templates.py

@@ -48,8 +48,6 @@ import lal
 import lalsimulation as lalsim
 
 
-from pylal import datatypes as laltypes
-from pylal import lalfft
 from pylal import spawaveform
 from gstlal import chirptime
 
@@ -109,21 +107,21 @@ class QuadraturePhase(object):
 	Example:
 
 	>>> import numpy
-	>>> from pylal.datatypes import REAL8TimeSeries
+	>>> import lal
 	>>> q = QuadraturePhase(128) # initialize for 128-sample templates
-	>>> inseries = REAL8TimeSeries(deltaT = 1.0 / 128, data = numpy.cos(numpy.arange(128, dtype = "double") * 2 * numpy.pi / 128)) # one cycle of cos(t)
+	>>> inseries = lal.CreateREAL8TimeSeries(deltaT = 1.0 / 128, length = 128)
+	>>> inseries.data.data = numpy.cos(numpy.arange(128, dtype = "double") * 2 * numpy.pi / 128) # one cycle of cos(t)
 	>>> outseries = q(inseries) # output has cos(t) in real part, sin(t) in imaginary part
 	"""
-
 	def __init__(self, n):
 		"""
 		Initialize.  n is the size, in samples, of the templates to
 		be processed.  This is used to pre-allocate work space.
 		"""
 		self.n = n
-		self.fwdplan = lalfft.XLALCreateForwardREAL8FFTPlan(n, 1)
-		self.revplan = lalfft.XLALCreateReverseCOMPLEX16FFTPlan(n, 1)
-		self.in_fseries = lalfft.prepare_fseries_for_real8tseries(laltypes.REAL8TimeSeries(deltaT = 1.0, data = numpy.zeros((n,), dtype = "double")))
+		self.fwdplan = lal.CreateForwardREAL8FFTPlan(n, 1)
+		self.revplan = lal.CreateReverseCOMPLEX16FFTPlan(n, 1)
+		self.in_fseries = lalfft.prepare_fseries_for_real8tseries(lal.CreateREAL8TimeSeries(deltaT = 1.0, length = n))
 
 
 	@staticmethod
@@ -137,18 +135,6 @@ class QuadraturePhase(object):
 		COMPLEX16FrequencySeries and n is the number of samples in
 		the original time series.
 		"""
-		#
-		# prepare output frequency series
-		#
-
-		out_fseries = laltypes.COMPLEX16FrequencySeries(
-			name = fseries.name,
-			epoch = fseries.epoch,
-			f0 = fseries.f0,	# caution: only 0 is supported
-			deltaF = fseries.deltaF,
-			sampleUnits = fseries.sampleUnits
-		)
-
 		#
 		# positive frequencies include Nyquist if n is even
 		#
@@ -159,14 +145,26 @@ class QuadraturePhase(object):
 		# shuffle frequency bins
 		#
 
-		positive_frequencies = fseries.data
+		positive_frequencies = numpy.array(fseries.data.data) # work with copy
 		positive_frequencies[0] = 0	# set DC to zero
 		zeros = numpy.zeros((len(positive_frequencies),), dtype = "cdouble")
 		if have_nyquist:
 			# complex transform never includes positive Nyquist
 			positive_frequencies = positive_frequencies[:-1]
 
-		out_fseries.data = numpy.concatenate((zeros, 2 * positive_frequencies[1:]))
+		#
+		# prepare output frequency series
+		#
+
+		out_fseries = lal.CreateCOMPLEX16FrequencySeries(
+			name = fseries.name,
+			epoch = fseries.epoch,
+			f0 = fseries.f0,	# caution: only 0 is supported
+			deltaF = fseries.deltaF,
+			sampleUnits = fseries.sampleUnits,
+			length = len(zeros) + len(positive_frequencies) - 1
+		)
+		out_fseries.data.data = numpy.concatenate((zeros, 2 * positive_frequencies[1:]))
 
 		return out_fseries
 
@@ -183,14 +181,14 @@ class QuadraturePhase(object):
 		# transform to frequency series
 		#
 
-		lalfft.XLALREAL8TimeFreqFFT(self.in_fseries, tseries, self.fwdplan)
+		lal.REAL8TimeFreqFFT(self.in_fseries, tseries, self.fwdplan)
 
 		#
 		# transform to complex time series
 		#
 
-		tseries = laltypes.COMPLEX16TimeSeries(data = numpy.zeros((self.n,), dtype = "cdouble"))
-		lalfft.XLALCOMPLEX16FreqTimeFFT(tseries, self.add_quadrature_phase(self.in_fseries, self.n), self.revplan)
+		tseries = lal.CreateCOMPLEX16TimeSeries(length = self.n)
+		lal.COMPLEX16FreqTimeFFT(tseries, self.add_quadrature_phase(self.in_fseries, self.n), self.revplan)
 
 		#
 		# done
@@ -200,21 +198,23 @@ class QuadraturePhase(object):
 
 
 def normalized_autocorrelation(fseries, revplan):
-	data = fseries.data
-	fseries = laltypes.COMPLEX16FrequencySeries(
+	data = fseries.data.data
+	fseries = lal.CreateCOMPLEX16FrequencySeries(
 		name = fseries.name,
 		epoch = fseries.epoch,
 		f0 = fseries.f0,
 		deltaF = fseries.deltaF,
 		sampleUnits = fseries.sampleUnits,
-		data = data * numpy.conj(data)
+		length = len(data)
 	)
-	tseries = laltypes.COMPLEX16TimeSeries(
-		data = numpy.empty((len(data),), dtype = "cdouble")
+	fseries.data.data = data * numpy.conj(data)
+	tseries = lal.CreateCOMPLEX16TimeSeries(
+		length = len(data)
 	)
-	lalfft.XLALCOMPLEX16FreqTimeFFT(tseries, fseries, revplan)
-	data = tseries.data
-	tseries.data = data / data[0]
+	tseries.data.data = numpy.empty((len(data),), dtype = "cdouble")
+	lal.COMPLEX16FreqTimeFFT(tseries, fseries, revplan)
+	data = tseries.data.data
+	tseries.data.data = data / data[0]
 	return tseries