From e86269e88091f2b08f5cc86d8729a748ae79d92f Mon Sep 17 00:00:00 2001
From: Chad Hanna <crh184@psu.edu>
Date: Tue, 2 May 2017 10:08:49 -0400
Subject: [PATCH] metric.py: implement central finite differencing for time
coordinate
---
gstlal-ugly/python/metric.py | 39 +++++++++++++++++++++++++-----------
1 file changed, 27 insertions(+), 12 deletions(-)
diff --git a/gstlal-ugly/python/metric.py b/gstlal-ugly/python/metric.py
index eecb94f582..101fd8b235 100644
--- a/gstlal-ugly/python/metric.py
+++ b/gstlal-ugly/python/metric.py
@@ -110,9 +110,10 @@ class Metric(object):
self.metric_tensor = None
self.metric_is_valid = False
self.revplan = lal.CreateReverseCOMPLEX16FFTPlan(self.working_length, 1)
- # 1 microsecond offset
- self.delta_t = 1e-6
+ # 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))
self.tseries = lal.CreateCOMPLEX16TimeSeries(
name = "workspace",
epoch = 0,
@@ -139,8 +140,8 @@ class Metric(object):
mc = (m1*m2)**.6 / (m1+m2)**.2 * 5e-6
return 1./numpy.pi / mc * (5./256. * mc / 1.)**(3./8.)
- flow = self.flow
- #flow = max(min(fmin(p[0], p[1]), self.flow), 10)
+ #flow = self.flow
+ flow = max(min(fmin(p[0], p[1]), self.flow), 10)
try:
parameters = {}
@@ -186,7 +187,9 @@ class Metric(object):
try:
if dt is None:
self.w1w2.data.data[:] = numpy.conj(w1.data.data) * w2.data.data
- else:
+ 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
m = numpy.real(numpy.abs(self.w1w2.data.data[:].sum())) / norm(w1.data.data) / norm(w2.data.data)
if m > 1.0000001:
@@ -209,14 +212,18 @@ class Metric(object):
plus_match = self.match(w1, self.waveform(center+x))
if plus_match is None:
print "\nhere plus\n"
- plus_match = self.match(w1, self.waveform(center+10*x))
+ center += 31.159*x
+ plus_match = self.match(w1, self.waveform(center))
#return self.__set_diagonal_metric_tensor_component(i, center, deltas * 2, g, w1)
+ factor = 4 * numpy.finfo(numpy.float32).eps / (1. - plus_match)
+ x *= factor**.5
+ plus_match = self.match(w1, self.waveform(center+x))
# get the negative side of the central difference
minus_match = self.match(w1, self.waveform(center-x))
if minus_match is None:
print "\nhere negative\n"
- minus_match = self.match(w1, self.waveform(center-10*x))
+ minus_match = self.match(w1, self.waveform(center-31.159*x))
#return self.__set_diagonal_metric_tensor_component(i, center, deltas * 2, g, w1)
# second order
@@ -232,8 +239,10 @@ class Metric(object):
def __set_tt_metric_tensor_component(self, center, w1):
- d2 = 1. - self.match(w1, w1, dt = True)
- return d2 / self.delta_t / self.delta_t, self.delta_t
+ 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
def __set_offdiagonal_metric_tensor_component(self, (i,j), center, deltas, g, w1):
# evaluate symmetrically
@@ -278,10 +287,16 @@ class Metric(object):
# this is a diagonal component or not
x = numpy.zeros(len(deltas))
x[j] = deltas[j]
+ fjj = -2 * g[j,j] * deltas[j]**2 + 2.0
+ 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")
+ return -0.5 * (plus_match + minus_match - ftt - fjj + 2.0) / 2 / delta_t / deltas[j]
- # Check the match
- d2 = 1 - self.match(w1, self.waveform(center+x), dt = True)
- return (d2 - g_tt * delta_t**2 - g[j,j] * deltas[j]**2) / 2 / delta_t / deltas[j]
+ #d2 = 1 - self.match(w1, self.waveform(center+x), dt = "pos")
+ #return (d2 - g_tt * delta_t**2 - g[j,j] * deltas[j]**2) / 2 / delta_t / deltas[j]
#def __call__(self, center, deltas = None, thresh = 1. * numpy.finfo(numpy.float32).eps):
def __call__(self, center, deltas = None, thresh = 1. * numpy.finfo(numpy.float32).eps):
--
GitLab