gstlal_inspiral_compute_dtdphideff_cov_matrix: output covariance matrix to .h5 file

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