gstlal_inspiral_create_prior_diststats

#!/usr/bin/env python
#
# Copyright (C) 2010--2013  Kipp Cannon, Chad Hanna
#
# 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 gstlal_inspiral_create_prior_diststats
# A program to create some prior likelihood data to seed an online analysis
#
# ### Command line interface
#	+ `--verbose`: Be verbose.
#	+ `--synthesize-injection-count` [N] (int): Synthesize an injection distribution with N injections (default is 1000000).
#	+ `--write-likelihood` [filename]: Write merged raw likelihood data to this file.
#	+ `--instrument`: Append to a list of instruments to create dist stats for.  Must be whatever instruments you intend to analyze.
#	+ `--trials-far-thresh` (float): Set the far threshold for the thresh parameter in the trials table.
#	+ `--background-prior` [N] (float): Include an exponential background prior with the maximum bin count = N (default is 1).

#
# =============================================================================
#
#                                   Preamble
#
# =============================================================================
#


import sys
from optparse import OptionParser


from glue import segments
from gstlal import far, inspiral
from glue.ligolw import ligolw, lsctables
from glue.ligolw import utils as ligolw_utils
from glue.ligolw.utils import process as ligolw_process
from glue.ligolw.utils import search_summary as ligolw_search_summary
from pylal.datatypes import LIGOTimeGPS


__author__ = "Chad Hanna <chad.hanna@ligo.org>"
__version__ = "git id %s" % ""	# FIXME
__date__ = ""	# FIXME


#
# =============================================================================
#
#                                 Command Line
#
# =============================================================================
#


def parse_command_line():
	parser = OptionParser(
		version = "Name: %%prog\n%s" % "" # FIXME
	)
	parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
	parser.add_option("-s", "--synthesize-injection-count", metavar = "N", default = 10000000, type = "int", help = "Synthesize an injection distribution with N injections. default 1000000")
	parser.add_option("--write-likelihood", metavar = "filename", help = "Write merged raw likelihood data to this file.")
	parser.add_option("--instrument", action = "append", help = "append to a list of instruments to create dist stats for.  Must be whatever instruments you intend to analyze")
	parser.add_option("-p", "--background-prior", metavar = "N", default = 1, type = "float", help = "include an exponential background prior with the maximum bin count = N, default 1")
	options, filenames = parser.parse_args()

	if filenames:
		raise ValueError("unrecognized command line options")

	if not options.instrument or len(options.instrument) < 2:
		raise ValueError("must specify at least two --instrument's")

	return options, filenames



#
# =============================================================================
#
#                                     Main
#
# =============================================================================
#


#
# command line
#


options, filenames = parse_command_line()


#
# create parameter distribution priors
#


coincparamsdistributions = far.ThincaCoincParamsDistributions()

# FIXME:  need segment lists and horizon distances
seglists = segments.segmentlistdict.fromkeys(options.instrument, segments.segmentlist([segments.segment(LIGOTimeGPS(0.), LIGOTimeGPS(1.))]))
horizon_distances = {
	"H1": 120.0,
	"L1": 120.0,
	"V1": 45.0
}

if options.background_prior != 0:
	coincparamsdistributions.add_background_prior(n = options.background_prior, segs = seglists, verbose = options.verbose)
coincparamsdistributions.add_foreground_prior(n = options.synthesize_injection_count, segs = seglists, horizon_distances = horizon_distances, verbose = options.verbose)


#
# done
#


xmldoc = ligolw.Document()
xmldoc.appendChild(ligolw.LIGO_LW())
process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_inspiral_create_prior_diststats", paramdict = {})
search_summary = ligolw_search_summary.append_search_summary(xmldoc, process, ifos = seglists.keys(), inseg = seglists.extent_all(), outseg = seglists.extent_all())
far.gen_likelihood_control_doc(xmldoc, process, coincparamsdistributions, None, seglists, comment = u"background and signal priors (no real data)")
ligolw_process.set_process_end_time(process)


#
# Make a fake horizon distance summ_values table
#


process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_inspiral", {})
# can't use .copy() method, need to make new one, because
# it might be in a database
summ_value_table = lsctables.New(lsctables.SummValueTable)

for instrument, dist in horizon_distances.items():
	row =summ_value_table.RowType()
	row.summ_value_id = summ_value_table.get_next_id()
	row.program = process.program
	row.process_id = process.process_id
	row.frameset_group = None
	row.segment_def_id = None
	# claim the PSD to be valid for a point in time
	row.segment = segments.segment(LIGOTimeGPS(0.), LIGOTimeGPS(1.))
	row.instruments = (instrument,)
	row.name = inspiral.CoincsDocument.summ_value_name_encode(1.4, 1.4, 8.0)
	row.value = dist
	row.error = None
	row.intvalue = None
	row.comment = u"horizon distance (Mpc)"
	summ_value_table.append(row)
summ_value_table._end_of_rows()
xmldoc.childNodes[0].appendChild(summ_value_table)
ligolw_process.set_process_end_time(process)

ligolw_utils.write_filename(xmldoc, options.write_likelihood, gz = options.write_likelihood.endswith(".gz"), verbose = options.verbose)