gstlal_iir_bank

#!/usr/bin/env python
#
# Copyright (C) 2011-2012  Shaun Hooper
#
# 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.

## @file
# Compute an IIR bank for inspiral waveforms
#
# ### Command line interface
#
#   + `--flow` [Hz] (float): Set the template low-frequency cut-off (default = 40.0).
#   + `--sampleRate` [Hz] (float): Set the sample rate of the IIR template bank (optional).
#   + `--padding` [pad] (float): Fractional amount to pad time slices (default = 1.1).
#   + `--epsilon` [pad] (float): Second order correction factor (default = 0.02).
#   + `--reference-psd` [filename]: Load the spectrum from this LIGO light-weight XML file (required).
#   + `--template-bank` [filename]: Set the name of the LIGO light-weight XML file from which to load the template bank (required).
#   + `--instrument` [ifo]: Set the instrument.
#   + `--output` [filename]: Set the filename in which to save the template bank (required).
#   + `--verbose`: Be verbose.
#   + `--downsample`: Choose if you want to downsample IIR bank (recommended).


import sys
import scipy
import numpy
from optparse import OptionParser

parser = OptionParser(description = __doc__)
parser.add_option("--flow", metavar = "Hz", type = "float", default = 40.0, help = "Set the template low-frequency cut-off (default = 40.0).")
parser.add_option("--sampleRate", metavar = "Hz", type = "float", help = "Set the sample rate of the IIR template bank (optional).")
parser.add_option("--padding", metavar = "pad", type = "float", default = 1.1, help = "Fractional amount to pad time slices.")
parser.add_option("--epsilon", metavar = "pad", type = "float", default = 0.02, help = "Second order correction factor.")
# FIXME add this when we know how to (by changing the parameters accoring to a function or optimizing on the fly or something)
#parser.add_option("--iir-tolerance", metavar = "match", type = "float", default = 0.99, help = "Set the SVD reconstruction tolerance (default = 0.99).")
parser.add_option("--reference-psd", metavar = "filename", help = "load the spectrum from this LIGO light-weight XML file (required).")
parser.add_option("--template-bank", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load the template bank (required).")
#FIXME figure out how to encode instrument info
parser.add_option("--instrument", metavar = "ifo", help = "Set the instrument")

parser.add_option("--output", metavar = "filename", help = "Set the filename in which to save the template bank (required).")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose (optional).")
parser.add_option("--downsample", action = "store_true", help = "Choose if you want to downsample IIR bank (recommended)")

options, filenames = parser.parse_args()

required_options = ("template_bank", "output")

missing_options = [option for option in required_options if getattr(options, option) is None]
if missing_options:
    raise ValueError, "missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in sorted(missing_options))


from glue.ligolw import ligolw
from glue.ligolw import array
from glue.ligolw import param
array.use_in(ligolw.LIGOLWContentHandler)
param.use_in(ligolw.LIGOLWContentHandler)
from glue.ligolw import utils
from gstlal import cbc_template_iir
from pylal.series import read_psd_xmldoc
from glue.ligolw import utils, lsctables

# read bank file
bank_xmldoc = utils.load_filename(options.template_bank, gz=options.template_bank.endswith('.gz'))
sngl_inspiral_table = lsctables.table.get_table(bank_xmldoc, lsctables.SnglInspiralTable.tableName)
fFinal = max(sngl_inspiral_table.getColumnByName("f_final"))

# read psd file
if options.reference_psd:
    ALLpsd = read_psd_xmldoc(utils.load_filename(options.reference_psd, verbose=options.verbose, contenthandler=ligolw.LIGOLWContentHandler))
    bank_sngl_table = lsctables.table.get_table( bank_xmldoc,lsctables.SnglInspiralTable.tableName )
    psd = ALLpsd[bank_sngl_table[0].ifo]
    # smooth and create an interp object
    psd = cbc_template_iir.smooth_and_interp(psd)
else:
    psd = None


# get the iir coefficients
A, B, D, snrs = cbc_template_iir.makeiirbank(bank_xmldoc, sampleRate = 4096, psd_interp = psd, verbose=options.verbose, padding=options.padding, flower=options.flow, downsample = options.downsample, output_to_xml = True, epsilon = options.epsilon)

utils.write_filename(bank_xmldoc, options.output, gz=options.output.endswith('.gz'))

#import pylab
#pylab.hist(snrs)
#pylab.show()