Skip to content
Snippets Groups Projects
Commit 209c5ed1 authored by Chad Hanna's avatar Chad Hanna
Browse files

WIP

parent c0d37aed
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
gstlal-ugly/bin/gstlal_inspiral_treebank
+ 16
− 10
View file @ 209c5ed1
......@@ -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:
......
gstlal-ugly/python/metric.py
+ 17
− 34
View file @ 209c5ed1
......@@ -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
......
gstlal-ugly/python/tree.py
+ 42
− 80
View file @ 209c5ed1
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment