add metric overlap code

gstlal-ugly/bin/Makefile.am

@@ -11,8 +11,11 @@ dist_bin_SCRIPTS = \
 	gstlal_harmonic_mean_psd \
 	gstlal_mock_data_server \
 	gstlal_ilwdify \
+	gstlal_inspiral_add_metric_overlaps \
 	gstlal_inspiral_best_coinc_file \
 	gstlal_inspiral_check_livetimes \
+	gstlal_inspiral_metric_overlap \
+	gstlal_inspiral_metric_overlap_dag \
 	gstlal_inspiral_treebank \
 	gstlal_inspiral_treebank_dag \
 	gstlal_inspiral_bankviz \

gstlal-ugly/bin/gstlal_inspiral_add_metric_overlaps

+#!/usr/bin/python
+import sys
+import numpy
+import argparse
+import h5py
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--out-h5-file", required = True, help = "provide the output hdf5 file name")
+parser.add_argument("infile", nargs="+", help = "provide the output hdf5 file name")
+
+args = parser.parse_args()
+
+output = []
+for f in args.infile:
+	h5f = h5py.File(f, 'r')
+	output += list(numpy.array(h5f["overlaps"])) 
+output = numpy.array(output)
+h5f = h5py.File(args.out_h5_file, 'w')
+h5f.create_dataset('overlaps', data = output)
+h5f.close()

gstlal-ugly/bin/gstlal_inspiral_metric_overlap

+#!/usr/bin/python
+import sys
+from gstlal import metric as metric_module
+from ligo.lw import ligolw
+from ligo.lw import utils as ligolw_utils
+from ligo.lw import lsctables
+import numpy
+import argparse
+import h5py
+
+class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
+	pass
+lsctables.use_in(LIGOLWContentHandler)
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--psd-xml-file", help = "provide a psd xml file")
+parser.add_argument("--bank-file", help = "provide the bank file for which overlaps will be calculated")
+parser.add_argument("--out-h5-file", required = True, help = "provide the output hdf5 file name")
+parser.add_argument("--start-row", type=int, help = "The starting row to calculate overlaps")
+parser.add_argument("--num-rows", type=int, help = "The number of rows to calculate overlaps")
+
+args = parser.parse_args()
+
+g_ij = metric_module.Metric(
+        args.psd_xml_file,
+        coord_func = metric_module.m1_m2_s1z_s2z_func,
+        duration = 1.0, # FIXME!!!!!
+        flow = 10,
+        fhigh = 1024,
+        approximant = "IMRPhenomD"
+)
+
+xmldoc = ligolw_utils.load_filename(args.bank_file, verbose = True, contenthandler = LIGOLWContentHandler)
+sngl_inspiral_table = lsctables.SnglInspiralTable.get_table(xmldoc)
+
+vec1s = numpy.array([[row.template_id, row.mass1, row.mass2, row.spin1z, row.spin2z] for row in sngl_inspiral_table[args.start_row:args.start_row+args.num_rows]])
+vec2s = numpy.array([[row.template_id, row.mass1, row.mass2, row.spin1z, row.spin2z] for row in sngl_inspiral_table])
+
+output = []
+
+for n, vec1 in enumerate(vec1s):
+	g, det = g_ij(vec1[1:])
+	def match(vec2, vec1 = vec1, g = g):
+		return (vec1[0], vec2[0], g_ij.metric_match(g, vec1[1:], vec2[1:]))
+	thisoutput = [row for row in map(match, vec2s)	if row[2] > 0.9]
+	print n, len(thisoutput)
+	output += thisoutput
+
+output = numpy.array(output)
+h5f = h5py.File(args.out_h5_file, 'w')
+h5f.create_dataset('overlaps', data = output)
+h5f.close()

gstlal-ugly/bin/gstlal_inspiral_metric_overlap_dag

+#!/usr/bin/python
+import sys
+from ligo.lw import ligolw
+from ligo.lw import utils as ligolw_utils
+from ligo.lw import lsctables
+import numpy
+import argparse
+import h5py
+from gstlal import pipeline
+from gstlal import dagparts
+
+class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
+	pass
+lsctables.use_in(LIGOLWContentHandler)
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--psd-xml-file", help = "provide a psd xml file")
+parser.add_argument("--bank-file", help = "provide the bank file for which overlaps will be calculated")
+args = parser.parse_args()
+
+xmldoc = ligolw_utils.load_filename(args.bank_file, verbose = True, contenthandler = LIGOLWContentHandler)
+sngl_inspiral_table = lsctables.SnglInspiralTable.get_table(xmldoc)
+
+try:
+	os.mkdir("logs")
+except:
+	pass
+dag = dagparts.DAG("metric_overlap")
+
+overlapJob = dagparts.DAGJob("gstlal_inspiral_metric_overlap", condor_commands = {"want_graceful_removal":"True", "kill_sig":"15", "accounting_group":"ligo.prod.o3.cbc.uber.gstlaloffline"})
+addJob = dagparts.DAGJob("gstlal_inspiral_add_metric_overlaps", condor_commands = {"want_graceful_removal":"True", "kill_sig":"15", "accounting_group":"ligo.prod.o3.cbc.uber.gstlaloffline"})
+
+num = 1000
+overlapnodes = []
+# FIXME dont hardcode 3345408, it comes from number of tiles in TimePhaseSNR
+for start in range(0, len(sngl_inspiral_table), num):
+        stop = start + num
+        overlapnodes.append(dagparts.DAGNode(overlapJob, dag, parent_nodes = [], opts = {"start-row":str(start), "num-rows":str(num)}, output_files = {"out-h5-file":"%s/metric_overlaps_%d_%d.h5" % (overlapJob.output_path, start, num)}, input_files = {"psd-xml-file": args.psd_xml_file, "bank-file": args.bank_file}))
+
+addnode = dagparts.DAGNode(addJob, dag, parent_nodes = overlapnodes, input_files = {"": [n.output_files["out-h5-file"] for n in overlapnodes]})
+
+dag.write_sub_files()
+dag.write_dag()
+dag.write_script()