From 209c5ed1e079f9386e7a167f776bbe4a7d1da02b Mon Sep 17 00:00:00 2001
From: Chad Hanna <chad.hanna@ligo.org>
Date: Sun, 24 Mar 2019 13:17:02 -0400
Subject: [PATCH] WIP
---
gstlal-ugly/bin/gstlal_inspiral_treebank | 26 +++--
gstlal-ugly/python/metric.py | 51 ++++------
gstlal-ugly/python/tree.py | 122 ++++++++---------------
3 files changed, 75 insertions(+), 124 deletions(-)
diff --git a/gstlal-ugly/bin/gstlal_inspiral_treebank b/gstlal-ugly/bin/gstlal_inspiral_treebank
index 3a5f775a5c..0c3d368326 100755
--- a/gstlal-ugly/bin/gstlal_inspiral_treebank
+++ b/gstlal-ugly/bin/gstlal_inspiral_treebank
@@ -19,12 +19,18 @@
import itertools
import numpy
-from ligo.lw import ligolw
-from ligo.lw import lsctables
-from ligo.lw import table
-from ligo.lw import utils
-from ligo.lw import ilwd
-from ligo.lw.utils import process as ligolw_process
+#from ligo.lw import ligolw
+#from ligo.lw import lsctables
+#from ligo.lw import table
+#from ligo.lw import utils
+#from ligo.lw import ilwd
+#from ligo.lw.utils import process as ligolw_process
+from glue.ligolw import ligolw
+from glue.ligolw import lsctables
+from glue.ligolw import table
+from glue.ligolw import utils
+from glue.ligolw import ilwd
+from glue.ligolw.utils import process as ligolw_process
from gstlal import metric as metric_module
from gstlal import tree
# FIXME dont do this
@@ -165,7 +171,7 @@ def coord_func(coords, positions = positions):
out[pi] = coords[i]
return out
-if True:#args.mcm2:
+if False:#args.mcm2:
coord_func = metric_module.M_q_func
else:
coord_func = coord_func
@@ -187,17 +193,17 @@ def trans_M_q(coord_limits):
m2limis = [coord_limits[1][0] * 0.9, coord_limits[1][1] * 1.256789]
return [Mlimits] + [m2limis] + coord_limits[2:]
-if True:#args.mcm2:
+if False:#args.mcm2:
coord_limits_2 = trans_M_q(coord_limits)
if args.min_mc:
coord_limits_2[0][0] = max(coord_limits_2[0][0], args.min_mc)
if args.max_mc:
coord_limits_2[0][1] = min(coord_limits_2[0][1], args.max_mc)
print coord_limits_2
- bank = Node(HyperCube(numpy.array(coord_limits_2), mismatch, constraint_func = lambda x: tree.mass_sym_constraint_mc(x, args.max_q, args.max_mtotal, coord_limits[0], coord_limits[1], coord_limits_2[0]), metric = g_ij, singularity = tree.mc_m2_singularity))
+ bank = Node(HyperCube(numpy.array(coord_limits_2), mismatch, constraint_func = lambda x: tree.mass_sym_constraint_mc(x, args.max_q, args.max_mtotal, coord_limits[0], coord_limits[1], coord_limits_2[0]), metric = g_ij))
else:
coord_limits_2 = coord_limits
- bank = Node(HyperCube(numpy.array(coord_limits_2), mismatch, constraint_func = lambda x: tree.mass_sym_constraint(x, args.max_q, args.max_mtotal), metric = g_ij, singularity = tree.m1_m2_singularity))
+ bank = Node(HyperCube(numpy.array(coord_limits_2), mismatch, constraint_func = lambda x: tree.mass_sym_constraint(x, args.max_q, args.max_mtotal), metric = g_ij))
if args.verbose:
diff --git a/gstlal-ugly/python/metric.py b/gstlal-ugly/python/metric.py
index bf40750f6e..63a9584e81 100644
--- a/gstlal-ugly/python/metric.py
+++ b/gstlal-ugly/python/metric.py
@@ -21,14 +21,15 @@ from lal import series
from scipy import integrate
import numpy
from gstlal import reference_psd
-from ligo.lw import utils as ligolw_utils
+#from ligo.lw import utils as ligolw_utils
+from glue.ligolw import utils as ligolw_utils
import itertools
import scipy
from lal import LIGOTimeGPS
import sys
import math
-DELTA = 2.5e-7#2.0e-7
+DELTA = 2e-7#2.0e-7
EIGEN_DELTA_DET = DELTA
# Round a number up to the nearest power of 2
@@ -156,7 +157,6 @@ class Metric(object):
self.flow = flow
self.fhigh = fhigh
self.working_length = int(round(self.duration * 2 * self.fhigh)) + 1
- print self.working_length
self.psd = reference_psd.interpolate_psd(series.read_psd_xmldoc(ligolw_utils.load_filename(psd_xml, verbose = True, contenthandler = series.PSDContentHandler)).values()[0], self.df)
self.metric_tensor = None
self.metric_is_valid = False
@@ -301,12 +301,13 @@ class Metric(object):
minus_match_minus_1 = self.match_minus_1(w1, wm[j,j], t_factor = self.neg_t_factor)
return -0.5 * (plus_match_minus_1 + minus_match_minus_1 - ftt - fjj) / 2 / delta_t / deltas[j]
- def __call__(self, center, deltas = None):
-
+ def __call__(self, center):#, deltas = None):
+ center = numpy.array(center)
g = numpy.zeros((len(center), len(center)), dtype=numpy.double)
w1 = self.waveform(center)
wp = {}
wm = {}
+ deltas = abs(center)**.5 * DELTA + DELTA
for i, x in enumerate(deltas):
for j, y in enumerate(deltas):
dx = numpy.zeros(len(deltas))
@@ -335,41 +336,21 @@ class Metric(object):
for i, j in itertools.product(range(len(deltas)), range(len(deltas))):
g[i,j] = g[i,j] - g_tj[i] * g_tj[j] / g_tt
- #print "before ", g
- #g = numpy.array(nearPD(g, nit=3))
- try:
- U, S, V = numpy.linalg.svd(g)
- except numpy.linalg.linalg.LinAlgError:
- newc = center.copy()
- newc[0:2] *=0.99
- return self.__call__(newc, deltas)
+ U, S, V = numpy.linalg.svd(g)
+ #try:
+ # U, S, V = numpy.linalg.svd(g)
+ #except numpy.linalg.linalg.LinAlgError:
+ # print "svd error"
+ # newc = center.copy()
+ # newc[0:2] *=0.99
+ # return self.__call__(newc)#, deltas)
condition = S < max(S) * EIGEN_DELTA_DET
- S[condition] = 0.0
+ S[condition] = EIGEN_DELTA_DET
g = numpy.dot(U, numpy.dot(numpy.diag(S), V))
det = numpy.product(S[S>0])
eff_dimension = len(S[S>0])
w = S
- #try:
- # w, v = numpy.linalg.eigh(g)
- #except numpy.linalg.linalg.LinAlgError:
- # newc = center.copy()
- # newc[0:2] -= deltas[0:2]
- # return self.__call__(newc, deltas)
- #mxw = numpy.max(w)
- #assert numpy.all(w > 0)
- #if numpy.any(w < 0):
- # print w
- # return self.__call__(center - deltas, deltas)
self.metric_is_valid = True
- #w[w<EIGEN_DELTA_DET * mxw] = EIGEN_DELTA_DET * mxw
- #det = numpy.product(w)
-
- #eff_dimension = len(w[w >= EIGEN_DELTA_DET * mxw])
- #w[w<EIGEN_DELTA_METRIC * mxw] = EIGEN_DELTA_METRIC * mxw
- #g = numpy.dot(numpy.dot(v, numpy.abs(numpy.diag(w))), v.T)
-
- #return g, eff_dimension, numpy.product(S[S>0])
- #print "after ", g
return g, eff_dimension, det, self.metric_is_valid, w
@@ -386,6 +367,8 @@ class Metric(object):
delta = x-y
return (dot(delta, delta, metric_tensor))**.5
+ def volume_element(self, metric_tensor):
+ return abs(numpy.linalg.det(metric_tensor))**.5
def metric_match(self, metric_tensor, c1, c2):
d2 = self.distance(metric_tensor, c1, c2)**2
diff --git a/gstlal-ugly/python/tree.py b/gstlal-ugly/python/tree.py
index a0e0adab80..daf185f31b 100644
--- a/gstlal-ugly/python/tree.py
+++ b/gstlal-ugly/python/tree.py
@@ -98,29 +98,9 @@ def packing_density(n):
if n==8:
return prefactor * numpy.pi**4 / 384
-def mc_m2_singularity(c):
- center = c.copy()
- m1, m2 = metric_module.m1_from_mc_m2(center[0], center[1]), center[1]
- if .80 < m1 / m2 <= 1.0:
- m1 = 0.80 * m2
- elif 1.0 < m1 / m2 <= 1.25:
- m2 = 0.80 * m1
- mc = metric_module.mc_from_m1_m2(m1, m2)
- center[0] = mc
- center[1] = m2
- return center
-
-def m1_m2_singularity(c):
- center = c.copy()
- return center
- F = 1.
- if F*.95 < center[0] / center[1] <= F * 1.05:
- center[1] = 0.95 * center[0]
- return center
-
class HyperCube(object):
- def __init__(self, boundaries, mismatch, constraint_func = mass_sym_constraint, metric = None, metric_tensor = None, effective_dimension = None, det = None, singularity = None, metric_is_valid = False, eigv = None):
+ def __init__(self, boundaries, mismatch, constraint_func = mass_sym_constraint, metric = None, metric_tensor = None, effective_dimension = None, det = None, metric_is_valid = False, eigv = None):
"""
Define a hypercube with boundaries given by boundaries, e.g.,
@@ -141,24 +121,8 @@ class HyperCube(object):
self.center = numpy.array([c[0] + (c[1] - c[0]) / 2. for c in boundaries])
self.deltas = numpy.array([c[1] - c[0] for c in boundaries])
self.metric = metric
- # FIXME don't assume m1 m2 and the spin coords are the coordinates we have here.
- deltas = DELTA * numpy.ones(len(self.center))
- deltas[0:2] *= self.center[0:2]**.5
- self.singularity = singularity
-
if self.metric is not None and metric_tensor is None:
- #try:
- if self.singularity is not None:
- center = self.singularity(self.center)
- else:
- center = self.center
- try:
- self.metric_tensor, self.effective_dimension, self.det, self.metric_is_valid, self.eigv = self.metric(center, deltas)
- except ValueError:
- center *= 0.99
- self.metric_tensor, self.effective_dimension, self.det, self.metric_is_valid, self.eigv = self.metric(center, deltas)
- # print "metric @", self.center, " failed, trying, ", self.center - self.deltas / 2.
- # self.metric_tensor, self.effective_dimension, self.det = self.metric(self.center - self.deltas / 2., deltas)
+ self.metric_tensor, self.effective_dimension, self.det, self.metric_is_valid, self.eigv = self.metric(self.center)
else:
self.metric_tensor = metric_tensor
self.effective_dimension = effective_dimension
@@ -212,9 +176,9 @@ class HyperCube(object):
leftbound[dim,1] = self.center[dim]
rightbound[dim,0] = self.center[dim]
if reuse_metric:
- return HyperCube(leftbound, self.__mismatch, self.constraint_func, metric = self.metric, metric_tensor = self.metric_tensor, effective_dimension = self.effective_dimension, det = self.det, singularity = self.singularity, metric_is_valid = self.metric_is_valid, eigv = self.eigv), HyperCube(rightbound, self.__mismatch, self.constraint_func, metric = self.metric, metric_tensor = self.metric_tensor, effective_dimension = self.effective_dimension, det = self.det, singularity = self.singularity, metric_is_valid = self.metric_is_valid, eigv = self.eigv)
+ return HyperCube(leftbound, self.__mismatch, self.constraint_func, metric = self.metric, metric_tensor = self.metric_tensor, effective_dimension = self.effective_dimension, det = self.det, metric_is_valid = self.metric_is_valid, eigv = self.eigv), HyperCube(rightbound, self.__mismatch, self.constraint_func, metric = self.metric, metric_tensor = self.metric_tensor, effective_dimension = self.effective_dimension, det = self.det, metric_is_valid = self.metric_is_valid, eigv = self.eigv)
else:
- return HyperCube(leftbound, self.__mismatch, self.constraint_func, metric = self.metric, singularity = self.singularity), HyperCube(rightbound, self.__mismatch, self.constraint_func, metric = self.metric, singularity = self.singularity)
+ return HyperCube(leftbound, self.__mismatch, self.constraint_func, metric = self.metric), HyperCube(rightbound, self.__mismatch, self.constraint_func, metric = self.metric)
def tile(self, mismatch, stochastic = False):
self.tiles.append(self.center)
@@ -241,8 +205,6 @@ class HyperCube(object):
def volume(self, metric_tensor = None):
if metric_tensor is None:
metric_tensor = self.metric_tensor
- #return numpy.product(self.deltas) * numpy.linalg.det(metric_tensor)**.5
- #print "volume ", numpy.product(self.deltas) * self.det**.5
return numpy.product(self.deltas) * self.det**.5
def coord_volume(self):
@@ -284,54 +246,54 @@ class Node(object):
if vertex not in self.boundary:
self.on_boundary = True
print "\n\non boundary!!\n\n"
- #self.on_boundary = numpy.any((self.cube.center + self.cube.deltas) == (self.boundary.center + self.boundary.deltas)) or numpy.any((self.cube.center - self.cube.deltas) == (self.boundary.center - self.boundary.deltas))
-
- def split(self, split_num_templates, mismatch, bifurcation = 0, verbose = True, metric_tol = 0.01, max_coord_vol = float(10)):
- size = self.cube.num_tmps_per_side(mismatch)
- #if max(self.cube.deltas) / min(self.cube.deltas) > 100.0:
- # self.splitdim = numpy.argmax(self.cube.deltas)
- #else:
- self.splitdim = numpy.argmax(size)
-
- if not self.parent:
- numtmps = float("inf")
- par_numtmps = float("inf")
- volume_split_condition = False
- metric_diff = 1.0
- metric_cond = True
- else:
- # Get the number of parent templates
- par_numtmps = self.parent.cube.num_templates(mismatch)
- # get the number of sibling templates
- sib_numtmps = self.sibling.cube.num_templates(mismatch)
+ def split(self, split_num_templates, mismatch, bifurcation = 0, verbose = True, metric_tol = 0.03, max_coord_vol = float(10)):
+ if self.cube.constraint_func(self.cube.vertices + [self.cube.center]):
+ size = self.cube.num_tmps_per_side(mismatch)
+ self.splitdim = numpy.argmax(size)
+
+ if not self.parent:
+ numtmps = float("inf")
+ par_numtmps = float("inf")
+ volume_split_condition = False
+ metric_diff = 1.0
+ metric_cond = True
+ else:
+ # Get the number of parent templates
+ par_numtmps = self.parent.cube.num_templates(mismatch)
- # get our number of templates
- numtmps = self.cube.num_templates(mismatch)
+ # get the number of sibling templates
+ sib_numtmps = self.sibling.cube.num_templates(mismatch)
+ # get our number of templates
+ numtmps = self.cube.num_templates(mismatch)
- metric_diff = max(abs(self.cube.deltas / self.cube.size - self.sibling.cube.deltas / self.sibling.cube.size) / (self.cube.deltas / self.cube.size))
- #metric_diff = self.cube.metric_tensor - self.sibling.cube.metric_tensor
- #metric_diff = numpy.linalg.norm(metric_diff) / numpy.linalg.norm(self.cube.metric_tensor)**.5 / numpy.linalg.norm(self.sibling.cube.metric_tensor)**.5
+ #print "self metric:\n", self.cube.metric_tensor
+ #print "sibling metric:\n", self.sibling.cube.metric_tensor
+ #print "parent metric:\n", self.parent.cube.metric_tensor
- metric_diff2 = max(abs(self.cube.deltas / self.cube.size - self.parent.cube.deltas / self.parent.cube.size) / self.cube.deltas / self.cube.size)
- #metric_diff2 = self.cube.metric_tensor - self.parent.cube.metric_tensor
- #metric_diff2 = numpy.linalg.norm(metric_diff2) / numpy.linalg.norm(self.cube.metric_tensor)**.5 / numpy.linalg.norm(self.parent.cube.metric_tensor)**.5
+ #print self.cube.metric_tensor, self.sibling.cube.metric_tensor, self.parent.cube.metric_tensor
+ metric_diff = self.cube.metric_tensor - self.sibling.cube.metric_tensor
+ metric_diff = numpy.linalg.norm(metric_diff) / numpy.linalg.norm(self.cube.metric_tensor)**.5 / numpy.linalg.norm(self.sibling.cube.metric_tensor)**.5
+ metric_diff2 = self.cube.metric_tensor - self.parent.cube.metric_tensor
+ metric_diff2 = numpy.linalg.norm(metric_diff2) / numpy.linalg.norm(self.cube.metric_tensor)**.5 / numpy.linalg.norm(self.parent.cube.metric_tensor)**.5
- metric_diff = max(metric_diff, metric_diff2)
- # take the bigger of self, sibling and parent
- numtmps = max(max(numtmps, par_numtmps/2.0), sib_numtmps)
+ #print metric_diff, metric_diff2
- #mts = [(numtmps, self.cube), (sib_numtmps, self.sibling.cube), (par_numtmps / 2., self.parent.cube)]
- reuse_metric = self.cube #max(mts)[1]
+ metric_diff = max(metric_diff, metric_diff2)
+ #print "metric diff: ", metric_diff, "\n"
+ # take the bigger of self, sibling and parent
+ numtmps = max(max(numtmps, par_numtmps/2.0), sib_numtmps)
- #if self.cube.constraint_func(self.cube.vertices + [self.cube.center]) and ((numtmps >= split_num_templates) or (metric_diff > metric_tol and numtmps > split_num_templates)) or bifurcation < 2:
- tmps_per_side = self.cube.num_tmps_per_side(mismatch)
- if self.cube.constraint_func(self.cube.vertices + [self.cube.center]):
+ reuse_metric = self.cube #max(mts)[1]
+
+ tmps_per_side = self.cube.num_tmps_per_side(mismatch)
# NOTE FIXME work out correct max templates per side
- if ((numtmps >= split_num_templates) or max(tmps_per_side) > 2 * len(size)**.5 ) or bifurcation < 2:
+ if (numtmps >= split_num_templates) or bifurcation < 2 or metric_diff > metric_tol:
+ #print "splitting"
bifurcation += 1
- if metric_diff <= metric_tol and (numtmps < 5**(len(size))) and self.cube.coord_volume() < max_coord_vol:
+ if metric_diff <= metric_tol:
+ print "reuse metric"
self.cube.metric_tensor = reuse_metric.metric_tensor
self.cube.effective_dimension = reuse_metric.effective_dimension
self.cube.det = reuse_metric.det
--
GitLab