From b089d4ec9b59ab22c24f5e0728acc36dffb7e2dd Mon Sep 17 00:00:00 2001
From: Chad Hanna <chad.hanna@ligo.org>
Date: Fri, 24 May 2019 12:36:10 -0400
Subject: [PATCH] gstlal_inspiral_compute_dtdphideff_cov_matrix: initial commit
---
gstlal-inspiral/bin/Makefile.am | 1 +
...lal_inspiral_compute_dtdphideff_cov_matrix | 86 +++++++++++++++++++
2 files changed, 87 insertions(+)
create mode 100755 gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
diff --git a/gstlal-inspiral/bin/Makefile.am b/gstlal-inspiral/bin/Makefile.am
index 6a3daf771b..7734a2dc41 100644
--- a/gstlal-inspiral/bin/Makefile.am
+++ b/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 \
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
new file mode 100755
index 0000000000..9ae89a9f8b
--- /dev/null
+++ b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
@@ -0,0 +1,86 @@
+#!/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)
--
GitLab