gstlal_inspiral_compute_dtdphideff_cov_matrix: initial commit

gstlal-inspiral/bin/Makefile.am

@@ -9,6 +9,7 @@ dist_bin_SCRIPTS = \
 	gstlal_inspiral_calc_likelihood \
 	gstlal_inspiral_calc_rank_pdfs \
 	gstlal_inspiral_coinc_extractor \
+	gstlal_inspiral_compute_dtdphideff_cov_matrix \
 	gstlal_inspiral_add_dt_dphi_snr_ratio_pdfs \
 	gstlal_inspiral_calc_snr \
 	gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs \

gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix

+#!/usr/bin/python
+#
+# Copyright (C) 2019 Chad Hanna Amit Reza
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+import argparse
+import itertools
+import numpy, scipy.interpolate
+from lal import series
+from ligo.lw import utils as ligolw_utils
+import sys
+
+
+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')
+args = parser.parse_args()
+
+refpsd = args.psd_xml
+rho = {"L1": args.L_snr, "H1": args.H_snr, "V1": args.V_snr}
+
+
+psd = series.read_psd_xmldoc(ligolw_utils.load_filename(refpsd, verbose = True, contenthandler = series.PSDContentHandler))
+psddict = {}
+for ifo in psd:
+	f = numpy.arange(len(psd[ifo].data.data)) * psd[ifo].deltaF
+	psddict[ifo] = scipy.interpolate.interp1d(f, psd[ifo].data.data)
+
+def fn(fmin, fmax, n, ifo, psddict = psddict):
+	farr = numpy.linspace(fmin, fmax, 100000)
+	df = farr[1] - farr[0]
+	return 4 * numpy.sum(farr**(-7/3.) * farr**n * df / psddict[ifo](farr))
+
+def moment(fmin, fmax, n, ifo, psddict = psddict):
+	return fn(fmin, fmax, n, ifo, psddict) / fn(fmin, fmax, 0, ifo, psddict)
+
+sigsqtt = {}
+sigsqpp = {}
+sigsqtp = {}
+
+for ifo in rho:
+	sigsqf = moment(10, 1024, 2, ifo) - moment(10, 1024, 1, ifo)**2
+	sigsqtt[ifo] = (1. / (2 * 3.14 * rho[ifo] * sigsqf**.5)**2)
+	sigsqpp[ifo] = moment(10, 1024, 1, ifo)**2 / (rho[ifo]**2 * sigsqf)
+	sigsqtp[ifo] = moment(10, 1024, 1, ifo) / (2 * 3.14 * rho[ifo]**2 * sigsqf)
+
+for combo in itertools.combinations(("H1", "L1", "V1"), 2):
+	a,b = combo
+	m11 = sigsqtt[a] + sigsqtt[b]
+	m22 = sigsqpp[a] + sigsqpp[b]
+	m12 = m21 = sigsqtp[a] + sigsqtp[b]
+	mat = numpy.array([[m11, m12], [m21, m22]])
+	print combo
+
+	#print "cov mat"
+	#print mat
+	#print "cov mat inv"
+	#print matinv
+
+	print "cholesky on "
+	matinv = numpy.linalg.inv(mat)
+	cholesky_transpose = numpy.linalg.cholesky(matinv).T
+	print cholesky_transpose
+
+	#val, vec = scipy.linalg.eigh(mat)
+	#print "sanity check"
+	#print mat
+	#print numpy.dot(vec.T, numpy.dot(numpy.diag(val), vec))
+	#answer = numpy.dot(vec, numpy.diag(val)**.5)
+	#print answer
+	#print "sanity check"
+	#print numpy.dot(answer, answer.T)