gstlal_inspiral_create_prior_diststats: add seed likelihood option

gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats

@@ -32,6 +32,7 @@ from optparse import OptionParser
 import numpy
 import scipy
 from collections import defaultdict
+import warnings
 
 from ligo.lw import ligolw
 from ligo.lw import lsctables
@@ -76,6 +77,7 @@ def parse_command_line():
 	# FIXME:  default must be identical to gstlal_inspiral's default
 	parser.add_option("--min-instruments", metavar = "count", type = "int", default = 2, help = "Set the minimum number of instruments that must contribute triggers to form a candidate (default = 2).")
 	parser.add_option("--write-likelihood", metavar = "filename", help = "Write merged raw likelihood data to this file.")
+	parser.add_option("--seed-likelihood", metavar = "filename", help = "Start with a likelihood file and only update certain components. This is incompatible with --coincidence-threshold, --min-instruments, --instrument, and --background-prior so these options will be ignored")
 	parser.add_option("--instrument", action = "append", help = "Append to a list of instruments to create dist stats for.  List 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")
 	parser.add_option("--df", metavar = "N", default = 40, help = "set the degrees of freedom for the background chisq prior: default 40. You can also use template bandwidth to set this by setting it to 'bandwidth' or the values analytically derived from auto-correlations by setting it to 'analytic'.")
@@ -91,15 +93,6 @@ def parse_command_line():
 
 	process_params = dict(options.__dict__)
 
-	if not options.instrument:
-		raise ValueError("must specify at least one --instrument")
-	options.instrument = set(options.instrument)
-
-	if options.min_instruments < 1:
-		raise ValueError("--min-instruments must be >= 1")
-	if options.min_instruments > len(options.instrument):
-		raise ValueError("--min-instruments is greater than the number of unique --instrument's")
-
 	if filenames:
 		raise ValueError("unrecognized arguments after options: %s" % " ".join(filenames))
 
@@ -158,6 +151,23 @@ def parse_command_line():
 		# don't let the bandwidth get too small
 		options.df = max(int(min(bandwidths)) + 1, 10) * 3
 
+	if options.seed_likelihood:
+		warnings.warn("--seed-likelihood given, the following options  will be ignored: --coincidence-threshold, --min-instruments, --instrument, and --background-prior")
+		options.coincidence_threshold = None
+		options.min_instruments = None
+		options.instrument = None
+		options.background_prior = None
+	else:
+		if not options.instrument:
+			raise ValueError("must specify at least one --instrument")
+		options.instrument = set(options.instrument)
+		if options.min_instruments < 1:
+			raise ValueError("--min-instruments must be >= 1")
+		if options.min_instruments > len(options.instrument):
+			raise ValueError("--min-instruments is greater than the number of unique --instrument's")
+
+
+
 	return options, process_params, filenames, template_ids, horizon_factors, lambda_eta_sums_median
 
 
@@ -179,26 +189,44 @@ options, process_params, filenames, template_ids, horizon_factors, lambda_eta_su
 
 
 #
-# initialize output document (records process start time)
+# Either make a new file or use a seed which contains background data from a real analysis
 #
 
 
-xmldoc = ligolw.Document()
-xmldoc.appendChild(ligolw.LIGO_LW())
-process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_inspiral_create_prior_diststats", instruments = options.instrument, paramdict = process_params)
+if not options.seed_likelihood:
+	#
+	# initialize output document (records process start time),
+	# create parameter distribution priors, and
+	# add the denominator if starting from scratch
+	#
+	xmldoc = ligolw.Document()
+	xmldoc.appendChild(ligolw.LIGO_LW())
+	process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_inspiral_create_prior_diststats", instruments = options.instrument, paramdict = process_params)
+	rankingstat = far.RankingStat(template_ids = template_ids, instruments = options.instrument, delta_t = options.coincidence_threshold, min_instruments = options.min_instruments, population_model_file = options.mass_model_file, dtdphi_file = options.dtdphi_file, horizon_factors = horizon_factors, idq_file = options.idq_file)
+	if options.background_prior > 0:
+		rankingstat.denominator.add_noise_model(number_of_events = options.background_prior)
+else:
+	#
+	# Otherwise we have a seed and we will only overide signal options for e.g., a rerank
+	#
+	rankingstat, _ = far.parse_likelihood_control_doc(far.ligolw_utils.load_filename(options.seed_likelihood, verbose = options.verbose, contenthandler = far.RankingStat.LIGOLWContentHandler))
+	rankingstat.numerator = inspiral_lr.LnSignalDensity(
+		template_ids = rankingstat.template_ids,
+		instruments = rankingstat.instruments,
+		delta_t = rankingstat.delta_t,
+		population_model_file = rankingstat.population_model_file if options.mass_model_file is None else options.mass_model_file,
+		dtdphi_file = rankingstat.dtdphi_file if options.dtdphi_file is None else options.dtdphi_file,
+		min_instruments = rankingstat.min_instruments,
+		horizon_factors = rankingstat.horizon_factors,
+		idq_file = rankingstat.idq_file if options.idq_file is None else options.idq_file
+		)
 
 
 #
-# create parameter distribution priors
+# Decide which signal model to add
 #
 
 
-rankingstat = far.RankingStat(template_ids = template_ids, instruments = options.instrument, delta_t = options.coincidence_threshold, min_instruments = options.min_instruments, population_model_file = options.mass_model_file, dtdphi_file = options.dtdphi_file, horizon_factors = horizon_factors, idq_file = options.idq_file)
-
-if options.background_prior > 0:
-	rankingstat.denominator.add_noise_model(number_of_events = options.background_prior)
-
-# Add the numerator
 if options.df == "analytic":
 	mismatch_range = (options.mismatch_min, options.mismatch_max)
 	rankingstat.numerator.add_signal_model_analytic(lambda_sum=lambda_eta_sums_median["lambda_sum"],