From 11d629c66b41ad429c2ff0128990601687efeb18 Mon Sep 17 00:00:00 2001
From: Chad Hanna <chad.hanna@comp-hd-002.gwave.ics.psu.edu>
Date: Mon, 10 Jun 2019 15:45:18 -0400
Subject: [PATCH] gstlal_inspiral_compute_dtdphideff_cov_matrix: output
covariance matrix to .h5 file
---
...lal_inspiral_compute_dtdphideff_cov_matrix | 32 +++++++++++++------
1 file changed, 23 insertions(+), 9 deletions(-)
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
index 66731d651d..27bf499713 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
+++ b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
@@ -19,15 +19,17 @@ import argparse
import itertools
import numpy, scipy.interpolate
from lal import series
+from gstlal.stats.inspiral_extrinsics import TimePhaseSNR
from ligo.lw import utils as ligolw_utils
import sys
-
+import h5py
parser = argparse.ArgumentParser(description = 'generate a dt dphi covariance matrix')
parser.add_argument('--psd-xml', help = 'XML containing HLV psd')
parser.add_argument('--H-snr', type = float, help = 'H characteristic SNR')
parser.add_argument('--L-snr', type = float, help = 'L characteristic SNR')
parser.add_argument('--V-snr', type = float, help = 'V characteristic SNR')
+parser.add_argument("--output", help = "set the output h5 file, e.g., covmat.h5")
args = parser.parse_args()
refpsd = args.psd_xml
@@ -73,23 +75,35 @@ for combo in itertools.combinations(("H1", "L1", "V1"), 2):
m22 = sigsqpp[a] + sigsqpp[b]
m12 = m21 = sigsqtp[a] + sigsqtp[b]
mat = numpy.array([[m11, m12], [m21, m22]])
- #print combo
- #print "cholesky on "
matinv = numpy.linalg.inv(mat)
cholesky_transpose = numpy.linalg.cholesky(matinv).T
- #print cholesky_transpose
transtt[frozenset(combo)] = cholesky_transpose[0,0]
transtp[frozenset(combo)] = cholesky_transpose[0,1]
transpt[frozenset(combo)] = cholesky_transpose[1,0]
transpp[frozenset(combo)] = cholesky_transpose[1,1]
+f = h5py.File(args.output, "w")
+h5_transtt = f.create_group("transtt")
+h5_transtp = f.create_group("transtp")
+h5_transpt = f.create_group("transpt")
+h5_transpp = f.create_group("transpp")
+h5_transdd = f.create_group("transdd")
+for group, mat in zip((h5_transtt, h5_transtp, h5_transpt, h5_transpp, h5_transdd), (transtt, transtp, transpt, transpp, transdd)):
+ for k,v in mat.items():
+ group.create_dataset(",".join(sorted(k)), data = float(v))
+
+combos = TimePhaseSNR.instrument_combos(("H1","L1","V1")) + (("H1",),("L1",),("V1",))
+norm = dict((frozenset(k), 0.) for k in combos)
+norm = {}
+h5_norm = f.create_group("norm")
+for k in combos:
+ h5_norm.create_dataset(",".join(sorted(k)), data = 1.0)
+ norm[k] = 1.0
+f.close()
+
print "transtt =", transtt
print "transtp =", transtp
print "transpt =", transpt
print "transpp =", transpp
print "transdd =", transdd
-
-#transtt = {frozenset(("H1", "L1")): 5.19313416e+04, frozenset(("H1", "V1")): 1.39663777e+04, frozenset(("L1", "V1")): 3.98320074e+04}
-#transpp = {frozenset(("H1", "L1")): 3.50152221e+00, frozenset(("H1", "V1")): 1.56080114e+00, frozenset(("L1", "V1")): 1.63248925e+00}
-#transtp = {frozenset(("H1", "L1")): -6.51410977e+01, frozenset(("H1", "V1")): -2.07581116e+01, frozenset(("L1", "V1")): -6.05964763e+01}
-#transpt = {frozenset(("H1", "L1")): 0.0, frozenset(("H1", "V1")): 0.0, frozenset(("L1", "V1")): 0.0}
+print "norm = ", norm
--
GitLab