From efaad0b157e42482768b78e4b5b0bea1b47351b8 Mon Sep 17 00:00:00 2001
From: "chad.hanna" <crh184@psu.edu>
Date: Sun, 13 Oct 2019 21:03:25 -0700
Subject: [PATCH] gstlal_inspiral_activation_counts: improve initial version
---
gstlal-ugly/bin/Makefile.am | 1 +
.../bin/gstlal_inspiral_activation_counts | 50 +++++++++----------
.../bin/gstlal_inspiral_activation_counts_dag | 41 +++++++++++++++
3 files changed, 67 insertions(+), 25 deletions(-)
create mode 100755 gstlal-ugly/bin/gstlal_inspiral_activation_counts_dag
diff --git a/gstlal-ugly/bin/Makefile.am b/gstlal-ugly/bin/Makefile.am
index b285edb312..f591355479 100644
--- a/gstlal-ugly/bin/Makefile.am
+++ b/gstlal-ugly/bin/Makefile.am
@@ -12,6 +12,7 @@ dist_bin_SCRIPTS = \
gstlal_mock_data_server \
gstlal_ilwdify \
gstlal_inspiral_activation_counts \
+ gstlal_inspiral_activation_counts_dag \
gstlal_inspiral_add_metric_overlaps \
gstlal_inspiral_best_coinc_file \
gstlal_inspiral_calc_likelihood_by_bin \
diff --git a/gstlal-ugly/bin/gstlal_inspiral_activation_counts b/gstlal-ugly/bin/gstlal_inspiral_activation_counts
index 9338b323bf..b86bd37198 100755
--- a/gstlal-ugly/bin/gstlal_inspiral_activation_counts
+++ b/gstlal-ugly/bin/gstlal_inspiral_activation_counts
@@ -33,8 +33,9 @@ parser = argparse.ArgumentParser()
parser.add_argument("--output", help = "provide the output file")
parser.add_argument("--psd-xml-file", help = "provide a psd xml file")
parser.add_argument("--svd-file", help = "provide the bank file for which overlaps will be calculated")
-parser.add_argument("--m1", action = "append", help = "provide the start:stop:num parameters for mass 1, e.g., 1:3:100")
-parser.add_argument("--m2", action = "append", help = "provide the start:stop:num parameters for mass 2, e.g., 1:3:100")
+parser.add_argument("--m1", action = "append", help = "provide the num:start:stop parameters for mass 1, e.g., 100:1:3")
+parser.add_argument("--m2", action = "append", help = "provide the num:start:stop parameters for mass 2, e.g., 100:1:3")
+parser.add_argument("--s1", action = "append", help = "provide the num:start:stop parameters for spin 1, e.g., 11:-1:1")
args = parser.parse_args()
@@ -66,46 +67,45 @@ def x_y_z_zn_from_row(row):
print >> outfile, "m1 start, m1 stop, m2 start, m2 stop, activation counts"
# compute template match functions:
-
-# put the templates in a useful form
-tmps = numpy.array([x_y_z_zn_from_row(row) for row in sngl_inspiral_table])
-
print "computing metrics"
-metrics = [(tmp, g_ij(tmp)[0]) for tmp in tmps]
-
+match_funcs = []
+for tmp in (x_y_z_zn_from_row(row) for row in sngl_inspiral_table):
+ def match(vec2, vec1 = tmp, g = g_ij(tmp)[0]):
+ return g_ij.pseudo_match(g, vec1, vec2)
+ match_funcs.append(match)
-for m1, m2 in zip(args.m1, args.m2):
+snrs = numpy.linspace(8., 20., 13)
+snrnorm = sum(snrs**-4)
- print m1, m2
+for m1, m2, s1 in zip(args.m1, args.m2, args.s1):
- m1s,m1e,m1n = [float(x) for x in m1.split(":")]
- m2s,m2e,m2n = [float(x) for x in m2.split(":")]
+ print m1, m2
+ m1n,m1s,m1e = [float(x) for x in m1.split(":")]
+ m2n,m2s,m2e = [float(x) for x in m2.split(":")]
+ s1n,s1s,s1e = [float(x) for x in s1.split(":")]
# uniform population in m1 / m2 space
signals = []
weights = []
len_signals = 0
- print "calculating signal pop"
- for m1,m2 in itertools.product(numpy.logspace(numpy.log10(m1s), numpy.log10(m1e), m1n), numpy.logspace(numpy.log10(m2s), numpy.log10(m2e), m2n)):
+ for m1,m2,s1 in itertools.product(numpy.logspace(numpy.log10(m1s), numpy.log10(m1e), m1n), numpy.logspace(numpy.log10(m2s), numpy.log10(m2e), m2n), numpy.linspace(s1s, s1e, s1n)):
# assume zero spin
len_signals += 1
if mchirp(m1, m2) <= max_mchirp and mchirp(m1, m2) >= min_mchirp:
- signals.append(numpy.array([metric_module.x_from_m1_m2_s1_s2(m1, m2, 0, 0), metric_module.y_from_m1_m2_s1_s2(m1, m2, 0, 0), 0, 0]))
+ signals.append(numpy.array([metric_module.x_from_m1_m2_s1_s2(m1, m2, s1, 0), metric_module.y_from_m1_m2_s1_s2(m1, m2, s1, 0), metric_module.z_from_m1_m2_s1_s2(m1, m2, s1, 0), 0]))
# Correct for the logarithmic distribution
weights.append(m1*m2)
+ weights = numpy.array(weights)
+ signals = numpy.array(signals)
count = 0
print "computing counts for %d signals" % len(signals)
-
- for n, (tmp, g) in enumerate(metrics):
- print "tmp %d/%d" % (n, len(metrics))
- def match(vec2, vec1 = tmp, g = g):
- return g_ij.pseudo_match(g, vec1, vec2)
- # Assume an SNR 10 signal...What should this be?? Some function like this?
- matches = map(match, signals)
- snrs = numpy.linspace(8., 20., 13)
- for rho in snrs:
- count += rho**-4 * sum(numpy.exp(-rho**2*(1.0 - x)**2/2.) * w for x, w in zip(matches, weights) if x <= 1.0) / len_signals / sum(snrs**-4)
+
+ if len(signals) > 0:
+ for n, match_func in enumerate(match_funcs):
+ matches = numpy.apply_along_axis(match_func, 1, signals)
+ for rho in snrs:
+ count += rho**-4 / snrnorm / len_signals * (numpy.exp(-rho**2*(1.0 - matches)**2/2.) * weights).sum()
print >> outfile, "%f, %f, %f, %f, %e" % (m1s, m1e, m2s, m2e, count)
diff --git a/gstlal-ugly/bin/gstlal_inspiral_activation_counts_dag b/gstlal-ugly/bin/gstlal_inspiral_activation_counts_dag
new file mode 100755
index 0000000000..c1dde62e27
--- /dev/null
+++ b/gstlal-ugly/bin/gstlal_inspiral_activation_counts_dag
@@ -0,0 +1,41 @@
+#!/usr/bin/python
+import os
+import argparse
+from gstlal import dagparts
+from lal.utils import CacheEntry
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--output", help = "provide the output file")
+parser.add_argument("--psd-xml-file", help = "provide a psd xml file")
+parser.add_argument("--svd-cache", help = "provide the svd bank cache.")
+parser.add_argument("--m1", action = "append", help = "provide the num:start:stop parameters for mass 1, e.g., 100:1:3")
+parser.add_argument("--m2", action = "append", help = "provide the num:start:stop parameters for mass 2, e.g., 100:1:3")
+parser.add_argument("--s1", action = "append", help = "provide the num:start:stop parameters for spin 1, e.g., 11:-1:1")
+parser.add_argument("--accounting-group", help = "set accounting group")
+parser.add_argument("--accounting-group-user", help = "set accounting group user")
+
+args = parser.parse_args()
+
+try:
+ os.mkdir("logs")
+except:
+ pass
+
+dag = dagparts.DAG("ac_counts")
+
+acJob = dagparts.DAGJob("gstlal_inspiral_activation_counts", condor_commands = {"accounting_group": args.accounting_group, "accounting_group_user": args.accounting_group_user})
+
+for bin, line in enumerate(open(args.svd_cache)):
+ dagparts.DAGNode(acJob,
+ dag,
+ parent_nodes = [],
+ opts = {"m1": args.m1, "m2": args.m2, "s1":args.s1},
+ input_files = {"psd-xml-file": args.psd_xml_file, "svd-file": CacheEntry(line).path},
+ output_files = {"output": "%04d-COUNTS.txt" % bin}
+ )
+
+dag.write_sub_files()
+dag.write_dag()
+dag.write_script()
+dag.write_cache()
+
--
GitLab