gstlal_fake_frames 15.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
#!/usr/bin/env python
#
# Copyright (C) 2011  Kipp Cannon, Chad Hanna, Drew Keppel
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# 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
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

Kipp Cannon's avatar
Kipp Cannon committed
19
import numpy
20 21 22 23
from optparse import OptionParser
import os
import sys

Chad Hanna's avatar
Chad Hanna committed
24 25 26 27 28
import gi
gi.require_version('Gst', '1.0')
from gi.repository import GObject, Gst
GObject.threads_init()
Gst.init(None)
29

Kipp Cannon's avatar
Kipp Cannon committed
30 31
import lal

32 33 34 35 36
from gstlal import pipeparts
from gstlal import reference_psd
from gstlal import simplehandler
from gstlal import datasource
from gstlal import multirate_datasource
37
from ligo.lw import utils as ligolw_utils
38

Chad Hanna's avatar
Chad Hanna committed
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
### This program will make fake data in a variety of ways.  Its input is anything
### supported by datasource.py.  You can additionally whiten the data or apply a
### frequency domain coloring filter.  It writes its output to frame files.
###
### Overview graph of the pipeline
### ------------------------------
###
### Gray boxes are optional and depend on the command line given
###
### .. graphviz::
###
###    digraph G {
###
###        // graph properties
###
###        rankdir=LR;
###        compound=true;
###        node [shape=record fontsize=10 fontname="Verdana"];
###        edge [fontsize=8 fontname="Verdana"];
###
###        // nodes
###
###        "mkbasicsrc()" [URL="\ref datasource.mkbasicsrc()"];
###        "whitened_multirate_src()" [label="whitened_multirate_src()", URL="\ref multirate_datasource.mkwhitened_multirate_src()", style=filled, color=lightgrey];
###        capsfilter [style=filled, color=lightgrey, URL="\ref pipeparts.mkcapsfilter()"];
###        resample [style=filled, color=lightgrey, URL="\ref pipeparts.mkresample()"];
###        audioamplify [style=filled, color=lightgrey, URL="\ref pipeparts.mkaudioamplify()", label="audioamplify \n if --data-source=white \n in order to have unit variance \n after resampling"];
###        lal_shift[style=filled, color=lightgrey, URL="\ref pipeparts.mkshift()", label="lal_shift \n [iff --shift provided]"];
###        progressreport1 [URL="\ref pipeparts.mkprogressreport()", label="progressreport \n [iff --verbose provided]"];
###        progressreport2 [style=filled, color=lightgrey, URL="\ref pipeparts.mkprogressreport()", label="progressreport \n [iff --verbose provided]"];
###        lal_firbank [style=filled, color=lightgrey, URL="\ref pipeparts.mkfirbank()", label="lal_firbank \n [iff --color-psd provided]"];
###        taginject [URL="\ref pipeparts.mktaginject()"];
###        lal_simulation [style=filled, color=lightgrey, URL="\ref pipeparts.mkinjections()", label="lal_simulation \n [iff --injections provided]"];
###        framecppchannelmux [URL="\ref pipeparts.mkframecppchannelmux()"];
###        framecppfilesink [URL="\ref pipeparts.mkframecppfilesink()"];
###
###        // connect it up
###
###        "mkbasicsrc()" -> progressreport1-> lal_shift  -> progressreport2;
###        progressreport2 -> "whitened_multirate_src()" [label="if --whiten given"];
###        "whitened_multirate_src()" -> lal_firbank;
###        progressreport2 -> resample [label="if --whiten not given"];
###        resample -> capsfilter;
###        capsfilter -> audioamplify
###        audioamplify -> lal_firbank;
###        lal_firbank -> taginject -> lal_simulation -> framecppchannelmux -> framecppfilesink;
###    }
###
###
### Usage cases
### -----------
###
### 1. Making initial LIGO colored noise, Advanced LIGO noise, etc for different
### instruments.   See datasource.append_options() for other choices::
###
###    $ gstlal_fake_frames --data-source=LIGO --channel-name=H1=FAKE-STRAIN \
###    --frame-type=H1_FAKE --gps-start-time=900000000 --gps-end-time=900005000 \
###    --output-path=testframes  --verbose
###
###    $ gstlal_fake_frames --data-source=AdvLIGO --channel-name=L1=FAKE-STRAIN \
###    --frame-type=L1_FAKE --gps-start-time=900000000 --gps-end-time=900005000 \
###    --output-path=testframes  --verbose
###
###
### 2. Custom colored noise.
###
### Obviously it is not practical to code up every possible noise curve to
### simulate as a custom data source.  However, it is possible to provide your
### own custom noise curve as an ASCII file with frequency in one column and
### strain/Hz in the second. e.g., early advanced ligo noise curve
### https://git.ligo.org/lscsoft/gstlal/raw/master/gstlal/share/early_aligo_asd.txt
### You will need to convert ASD text files to PSD xml files using
### gstlal_psd_xml_from_asd_txt first.::
###
###	$ gstlal_fake_frames --data-source=white --color-psd v1psd.xml.gz \
###	--channel-name=V1=FAKE-STRAIN --frame-type=V1_FAKE --gps-start-time=900000000 \
###	--gps-end-time=900005000 --output-path=testframes  --verbose
###
###
### 3. Recoloring existing frame data
###
### This procedure is very similar to the above except that instead of
### using white noise to drive the fake frame generation, we start with
### real frame data and whiten it.  Recoloring can be thought of as first
### whitening data and then applying a coloring filter.  You must first
### make a reference PSD of the data with e.g. gstlal_reference_psd.  You
### will need to make a frame cache of the frames to recolor.::
###
###	$ gstlal_fake_frames --whiten-reference-psd reference_psd.xml.gz \
###	--color-psd recolor_psd.xml.gz --data-source frames \
###	--output-path /path/to/output --output-channel-name h_16384Hz \
###	--gps-start-time 966384031 --frame-type T1300121_V1_EARLY_RECOLORED_V2 \
###	--gps-end-time 966389031 --frame-duration 16 --frames-per-file 256 \
###	--frame-cache frame.cache --channel-name=V1=h_16384Hz --verbose
###
###
### 4. Creating injections into silence (assumes you have an injection xml file from e.g. lalapps_inspinj)::
###
###	$ gstlal_fake_frames --data-source silence --output-path /path/to/output \
###	--gps-start-time 966384031 --frame-type V1_INJECTIONS \
###	--gps-end-time 966389031 --frame-duration 16 --frames-per-file 256 \
###	--verbose --channel-name=V1=INJECTIONS --injections injections.xml
###
###
### 5. Other things are certainly possible, please add some!
###
### Debug
### -----
###
### It is possible to check the pipeline graph by interupting the running process
### with ctrl+c if you have the GST_DEBUG_DUMP_DOT_DIR evironment set.  A dot
### graph will be written to gstlal_fake_frames.  Here is an example::
###
###	$ GST_DEBUG_DUMP_DOT_DIR="." gstlal_fake_frames --data-source=LIGO \
###	--channel-name=H1=FAKE-STRAIN --frame-type=H1_FAKE \
###	--gps-start-time=900000000 --gps-end-time=900005000 \
###	--output-path=testframes  --verbose``
###
### You should see::
###
###	Wrote pipeline to ./gstlal_fake_frames_PLAYING.dot
###
### After Ctrl+c you should see::
###
###	^C*** SIG 2 attempting graceful shutdown (this might take several minutes) ... ***
###		Wrote pipeline to ./gstlal_fake_frames.dot
###
### You can then turn the pipeline graph into an image with dot, e.g.,::
###
###	$ dot -Tpng gstlal_fake_frames_PLAYING.dot > test.png
###
### Review
### ------
###
173

174

175 176 177
##
# Extract and validate the command line options
#
178 179
def parse_command_line():
	parser = OptionParser(description = __doc__)
180

181 182 183
	#
	# Append data source options
	#
184

185
	datasource.append_options(parser)
186

187 188 189
	#
	# Append program specific options
	#
190

191
	parser.add_option("--shift", metavar = "ns", help = "Number of nanoseconds to delay (negative) or advance (positive) the time stream", type = "int")
192
	parser.add_option("--sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate at which to generate the data, should be less than or equal to the sample rate of the measured psds provided, default = 16384 Hz, max 16384 Hz")
193
	parser.add_option("--whiten-reference-psd", metavar = "name", help = "Set the name of psd xml file to whiten the data with")
194
	parser.add_option("--whiten-track-psd", action = "store_true", help = "Calculate PSD from input data and track with time.  Always enabled if --whiten-reference-psd is not given.")
195 196 197 198 199 200 201 202 203 204 205
	parser.add_option("--color-psd", metavar = "name", help = "Set the name of psd xml file to color the data with")
	parser.add_option("--output-path", metavar = "name", default = ".", help = "Path to output frame files (default = \".\").")
	parser.add_option("--output-channel-name", metavar = "name", help = "The name of the channel in the output frames. The default is the same as the channel name")
	parser.add_option("--frame-type", metavar = "name", help = "Frame type, required")
	parser.add_option("--frame-duration", metavar = "s", default = 16, type = "int", help = "Set the duration of the output frames.  The duration of the frame file will be multiplied by --frames-per-file.  Default: 16s")
	parser.add_option("--frames-per-file", metavar = "n", default = 256, type = "int", help = "Set the number of frames per file.  Default: 256")
	parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose (optional).")

	#
	# Parse options
	#
206

207 208
	options, filenames = parser.parse_args()

209 210 211
	if options.sample_rate > 16384:
		raise ValueError("--sample-rate must be <= 16384")

212 213 214
	if options.frame_type is None:
		raise ValueError("--frame-type is required")

215 216 217
	if options.whiten_reference_psd is None:
		options.whiten_track_psd = True

218
	return options, filenames
219 220 221


#
222
# Main
223 224
#

225
options, filenames = parse_command_line()
226

227
## Parse the command line options into a python.datasource.GWDataSourceInfo class instance
228 229
gw_data_source = datasource.GWDataSourceInfo(options)

230
## Assume instrument is the first and only key of the python.datasource.GWDataSourceInfo.channel_dict
231 232 233 234 235 236 237 238 239 240
instrument = gw_data_source.channel_dict.keys()[0]

# disable progress reports if not verbose
if not options.verbose:
	pipeparts.mkprogressreport = lambda pipeline, src, *args: src

# set default output channel if not set by user
if options.output_channel_name is None:
	options.output_channel_name = gw_data_source.channel_dict[instrument]

241
# do not do injections in datasource.mkbasicsrc(), we will do them later
242
injections, gw_data_source.injection_filename = options.injections, None
243

244
## Setup the pipeline
Chad Hanna's avatar
Chad Hanna committed
245
pipeline = Gst.Pipeline(name=os.path.split(sys.argv[0])[1])
246 247

## Main loop 
Chad Hanna's avatar
Chad Hanna committed
248
mainloop = GObject.MainLoop()
249

250 251
## An instance of the python.simplehandler.Handler class
handler = simplehandler.Handler(mainloop, pipeline)
252

253
## 1) Set the pipeline head to basic input from datasource.mkbasicsrc()
254 255
# FIXME: fake source causes problems when making large buffers, so block_size needs to be overwritten
gw_data_source.block_size = 8 * options.sample_rate
256
head, _, _ = datasource.mkbasicsrc(pipeline, gw_data_source, instrument, verbose = options.verbose)
257

258
## 2) Set the pipeline head to be verbose with pipeparts.mkprogressreport()
259
head = pipeparts.mkprogressreport(pipeline, head, "frames")
260 261

if options.shift is not None:
262
	## 3) Set the pipeline head to apply a time shift if requested with a pipeparts.mkshift() element
263 264
	head = pipeparts.mkshift(pipeline, head, shift = options.shift)

265
	## 4) Set the pipeline head to be verbose with a pipeparts.mkprogressreport() element
266
	head = pipeparts.mkprogressreport(pipeline, head, "frames_shifted")
267

268
if options.whiten_reference_psd:
269
	## If whitening read the PSD
Kipp Cannon's avatar
Kipp Cannon committed
270
	wpsd = lal.series.read_psd_xmldoc(ligolw_utils.load_filename(options.whiten_reference_psd, verbose = options.verbose, contenthandler = lal.series.PSDContentHandler))[instrument]
271
else:
272
	## else set wpsd to None
273 274 275
	wpsd = None

if options.whiten_reference_psd or options.whiten_track_psd:
276
	## 5) Set the pipeline head to a whitened data stream if requested using a multirate_datasource.mkwhitened_multirate_src()
277
	head = multirate_datasource.mkwhitened_multirate_src(pipeline, head, [options.sample_rate], instrument, psd = wpsd, track_psd = options.whiten_track_psd)[options.sample_rate]
278
else:
279
	## 6) Otherwise simply add a pipeparts.mkcapsfilter() and pipeparts.mkresample()
Chad Hanna's avatar
Chad Hanna committed
280
	head = pipeparts.mkcapsfilter(pipeline, pipeparts.mkresample(pipeline, head, quality = 9), "audio/x-raw, rate=%d" % options.sample_rate)
281 282 283 284
	# FIXME this is a bit hacky, should datasource.mkbasicsrc be patched to change the sample rate?
	# FIXME don't hardcode sample rate (though this is what datasource generates for all fake data, period 
	# Renormalize if datasource is "white"
	if options.data_source == "white":
285
		head = pipeparts.mkaudioamplify(pipeline, head, 1.0 / (pipeparts.audioresample_variance_gain(9, 16384, options.sample_rate))**.5)
286 287 288 289

# Apply a coloring filter
if options.color_psd:

290
	## read coloring psd file and convert to an FIR filter
Kipp Cannon's avatar
Kipp Cannon committed
291
	rpsd = lal.series.read_psd_xmldoc(ligolw_utils.load_filename(options.color_psd, verbose = options.verbose, contenthandler = lal.series.PSDContentHandler))[instrument]
292 293
	
	## Calculate the maximum sample rate
294 295
	max_sample = int(round(1.0 / rpsd.deltaF * options.sample_rate / 2.0)) + 1

296
	# Truncate to requested output sample rate, if it is higher than the psd provides an assert will fail later
297 298 299
	rpsd_data = numpy.array(rpsd.data.data[:max_sample])
	rpsd = lal.CreateCOMPLEX16FrequencySeries(name = rpsd.name, epoch = rpsd.epoch, f0 = rpsd.f0, deltaF = rpsd.deltaF, length = max_sample, sampleUnits = rpsd.sampleUnits)
	rpsd.data.data = rpsd_data
300
	
301
	# create the coloring FIR kernel from reference_psd.psd_to_fir_kernel()
Chad Hanna's avatar
Chad Hanna committed
302 303
	FIRKernel = reference_psd.PSDFirKernel()
	fir_matrix, latency, measured_sample_rate = FIRKernel.psd_to_linear_phase_whitening_fir_kernel(rpsd, invert = False)
304

305
	## 7) Set the pipeline head to a pipeparts.mkfirbank() recoloring filter
306
	head = pipeparts.mkfirbank(pipeline, head, latency = latency, fir_matrix = [fir_matrix], block_stride = 1 * options.sample_rate)
307 308


309
## Set the tags for the output data
310
tagstr = "units=strain,channel-name=%s,instrument=%s" % (options.output_channel_name, instrument)
311 312

## 8) Put the units back to strain before writing to frames.  Additionally, override the output channel name if provided from the command line
313 314 315 316
head = pipeparts.mktaginject(pipeline, head, tagstr)


if injections is not None:
317
	## 9) Do injections into the final fake data
318 319
	head = pipeparts.mkinjections(pipeline, head, injections)

320 321
if not os.path.isdir(options.output_path):
	try:
Chad Hanna's avatar
Chad Hanna committed
322
		os.makedirs(options.output_path)
323
	except:
Chad Hanna's avatar
Chad Hanna committed
324
		print >> sys.stderr, "Unable to make directory ", options.output_path
325 326 327
		pass
else:
	print "Target output directory already exists."
328

329
## 10) create frames
330
head = pipeparts.mkframecppchannelmux(pipeline, {"%s:%s" % (instrument, options.output_channel_name): head}, frame_duration = options.frame_duration, frames_per_file = options.frames_per_file)
331 332

## 11) Write the frames to disk
333
head = pipeparts.mkframecppfilesink(pipeline, head, frame_type = options.frame_type)
334

335 336
# Put O(100000 s) frames in each directory
head.connect("notify::timestamp", pipeparts.framecpp_filesink_ldas_path_handler, (options.output_path, 5))
337 338

# Run it
339 340 341
if pipeline.set_state(Gst.State.READY) == Gst.StateChangeReturn.FAILURE:
	raise RuntimeError("pipeline failed to enter READY state")
datasource.pipeline_seek_for_gps(pipeline, gw_data_source.seg[0], gw_data_source.seg[1])
Chad Hanna's avatar
Chad Hanna committed
342
if pipeline.set_state(Gst.State.PLAYING) == Gst.StateChangeReturn.FAILURE:
343
	raise RuntimeError("pipeline failed to enter PLAYING state")
344 345 346 347 348 349 350 351 352 353 354 355

## Debugging output
if "GST_DEBUG_DUMP_DOT_DIR" in os.environ:
	pipeparts.write_dump_dot(pipeline, "%s_PLAYING" % pipeline.get_name(), verbose = True)

	## Setup a signal handler to intercept SIGINT in order to write the pipeline graph at ctrl+C before cleanly shutting down
	class SigHandler(simplehandler.OneTimeSignalHandler):
		def do_on_call(self, signum, frame):
			pipeparts.write_dump_dot(pipeline, "%s_SIGINT" % pipeline.get_name(), verbose = True)

	sighandler = SigHandler(pipeline)

356
mainloop.run()