From 83bfb64e60c0b5efc99b350aaa28a91fbd93e975 Mon Sep 17 00:00:00 2001
From: Chad Hanna <crh184@psu.edu>
Date: Thu, 18 May 2017 12:34:15 -0400
Subject: [PATCH] metric.py: put in proper t spacing
---
gstlal-ugly/python/metric.py | 67 +++++++++++++++++++++++++++---------
1 file changed, 50 insertions(+), 17 deletions(-)
diff --git a/gstlal-ugly/python/metric.py b/gstlal-ugly/python/metric.py
index eb8937e86f..8b79b2c0f0 100644
--- a/gstlal-ugly/python/metric.py
+++ b/gstlal-ugly/python/metric.py
@@ -25,6 +25,18 @@ import itertools
import scipy
from lal import LIGOTimeGPS
import sys
+import math
+
+# Round a number up to the nearest power of 2
+def ceil_pow_2(x):
+ x = int(math.ceil(x))
+ x -= 1
+ n = 1
+ while n and (x & (x + 1)):
+ x |= x >> n
+ n *= 2
+ return x + 1
+
def add_quadrature_phase(fseries, n):
"""
@@ -110,10 +122,15 @@ class Metric(object):
self.metric_tensor = None
self.metric_is_valid = False
self.revplan = lal.CreateReverseCOMPLEX16FFTPlan(self.working_length, 1)
- # 3 microsecond offset
- self.delta_t = 3e-6
- self.t_factor = numpy.exp(-2j * numpy.pi * (numpy.arange(self.working_length) * self.df - self.fhigh) * self.delta_t)
- self.neg_t_factor = numpy.exp(-2j * numpy.pi * (numpy.arange(self.working_length) * self.df - self.fhigh) * (-self.delta_t))
+ #FIXME
+ self.delta_t = {}
+ self.t_factor = {}
+ self.neg_t_factor = {}
+ for t in numpy.array([1.,2.,4.,8.,16.,32.,64.,128.,256.,512.,1024]):
+ self.delta_t[t] = 5e-6 * 1 * t # 1 M time spacing
+ print self.delta_t[t]
+ self.t_factor[t] = numpy.exp(-2j * numpy.pi * (numpy.arange(self.working_length) * self.df - self.fhigh) * self.delta_t[t])
+ self.neg_t_factor[t] = numpy.exp(-2j * numpy.pi * (numpy.arange(self.working_length) * self.df - self.fhigh) * (-self.delta_t[t]))
self.tseries = lal.CreateCOMPLEX16TimeSeries(
name = "workspace",
epoch = 0,
@@ -179,18 +196,16 @@ class Metric(object):
return None
return fseries
- def match(self, w1, w2, dt = None):
+ def match(self, w1, w2, t_factor = None):
def norm(w):
n = numpy.real((numpy.conj(w) * w).sum())**.5
return n
try:
- if dt is None:
+ if t_factor is None:
self.w1w2.data.data[:] = numpy.conj(w1.data.data) * w2.data.data
- elif dt == "neg":
- self.w1w2.data.data[:] = numpy.conj(w1.data.data) * w2.data.data * self.neg_t_factor
- elif dt == "pos":
- self.w1w2.data.data[:] = numpy.conj(w1.data.data) * w2.data.data * self.t_factor
+ else:
+ self.w1w2.data.data[:] = numpy.conj(w1.data.data) * w2.data.data * t_factor
m = numpy.real(numpy.abs(self.w1w2.data.data[:].sum())) / norm(w1.data.data) / norm(w2.data.data)
if m > 1.0000001:
raise ValueError("Match is greater than 1 : %f" % m)
@@ -231,10 +246,13 @@ class Metric(object):
def __set_tt_metric_tensor_component(self, center, w1):
- minus_match = self.match(w1, w1, dt = "neg")
- plus_match = self.match(w1, w1, dt = "pos")
- d2mbydx2 = (plus_match + minus_match - 2.0) / self.delta_t**2
- return -0.5 * d2mbydx2, self.delta_t
+ # FIXME FIXME assumes m1 and m2 are first coords
+ M = center[0] + center[1]
+ ix = ceil_pow_2(M)
+ minus_match = self.match(w1, w1, t_factor = self.neg_t_factor[ix])
+ plus_match = self.match(w1, w1, t_factor = self.t_factor[ix])
+ d2mbydx2 = (plus_match + minus_match - 2.0) / self.delta_t[ix]**2
+ return -0.5 * d2mbydx2, self.delta_t[ix]
def __set_offdiagonal_metric_tensor_component(self, (i,j), center, deltas, g, w1):
# evaluate symmetrically
@@ -261,6 +279,9 @@ class Metric(object):
return None
def __set_offdiagonal_time_metric_tensor_component(self, j, center, deltas, g, g_tt, delta_t, w1):
+ # FIXME FIXME assumes m1 and m2 are first coords
+ M = center[0] + center[1]
+ ix = ceil_pow_2(M)
# make the vector to solve for the metric by choosing
# either a principle axis or a bisector depending on if
# this is a diagonal component or not
@@ -270,8 +291,8 @@ class Metric(object):
ftt = -2 * g_tt * delta_t**2 + 2.0
# Second order
- plus_match = self.match(w1, self.waveform(center+x), dt="pos")
- minus_match = self.match(w1, self.waveform(center-x), dt="neg")
+ plus_match = self.match(w1, self.waveform(center+x), t_factor = self.t_factor[ix])
+ minus_match = self.match(w1, self.waveform(center-x), t_factor = self.neg_t_factor[ix])
return -0.5 * (plus_match + minus_match - ftt - fjj + 2.0) / 2 / delta_t / deltas[j]
def __call__(self, center, deltas = None, thresh = 1. * numpy.finfo(numpy.float32).eps):
@@ -323,4 +344,16 @@ class Metric(object):
def explicit_match(self, c1, c2):
- return self.match(self.waveform(c1), self.waveform(c2))
+ def fftmatch(w1, w2):
+ def norm(w):
+ n = numpy.real((numpy.conj(w) * w).sum())**.5 / self.duration**.5
+ return n
+
+ self.w1w2.data.data[:] = numpy.conj(w1.data.data) * w2.data.data
+ lal.COMPLEX16FreqTimeFFT(self.tseries, self.w1w2, self.revplan)
+ m = numpy.real(numpy.abs(numpy.array(self.tseries.data.data)).max()) / norm(w1.data.data) / norm(w2.data.data)
+ if m > 1.0000001:
+ raise ValueError("Match is greater than 1 : %f" % m)
+ return m
+
+ return fftmatch(self.waveform(c1), self.waveform(c2))
--
GitLab