diff --git a/bin/grid_random.py b/bin/grid_random.py deleted file mode 100755 index bdafe64..0000000 --- a/bin/grid_random.py +++ /dev/null @@ -1,28 +0,0 @@ -#!/usr/bin/python - -import sys -import pickle - -from grid import simple_rect_grid, simple_random_grid - -pfile = '/tmp/grid_random.p' -qfile = '/tmp/grid_random.txt' -MAX_POINTS = 200 - -if __name__ == '__main__': - sys.setrecursionlimit(4096) - print sys.getrecursionlimit() - try: - total_points = int(sys.argv[1]) - except: - total_points = 10 - - if total_points > MAX_POINTS: - print "too many points" - sys.exit(1) - - g = simple_random_grid(total_points) - g.construct_connectivity() - print g - open(qfile, 'w').write(g.for_qhull()) - pickle.dump(g, open(pfile, 'w')) diff --git a/bin/grid_regular.py b/bin/grid_regular.py index 719afaa..22772ef 100755 --- a/bin/grid_regular.py +++ b/bin/grid_regular.py @@ -5,7 +5,6 @@ import pickle from grid import simple_rect_grid, simple_random_grid -pfile = '/tmp/grid_regular.p' qfile = '/tmp/grid_regular.txt' if __name__ == '__main__': @@ -17,6 +16,6 @@ if __name__ == '__main__': resolution = 3 g = simple_rect_grid(resolution, resolution) - g.construct_connectivity() + print g + print g.get_points_conn([0.4, 0]) open(qfile, 'w').write(g.for_qhull()) - pickle.dump(g, open(pfile, 'w')) diff --git a/lib/baker.py b/lib/baker.py index 7595f99..4adc1a2 100644 --- a/lib/baker.py +++ b/lib/baker.py @@ -67,7 +67,7 @@ def get_phis_3D(X, r): def qlinear(X, R, q): """ this calculates the linear portion of q from X to r - + also, this is baker eq 3 X = destination point diff --git a/lib/grid.py b/lib/grid.py index fbc4361..fd4a192 100755 --- a/lib/grid.py +++ b/lib/grid.py @@ -7,7 +7,7 @@ from collections import defaultdict import numpy as np import scipy.spatial -from baker import run_baker +from baker import run_baker, get_phis from tools import exact_func from smcqdelaunay import * @@ -32,6 +32,28 @@ class face(object): ", ".join(neighbors) ) + + +class simplex(object): + def __init__(self, points): + self.R = points + + def contains(self, X): + phis = get_phis(X, self.R) + + r = True + if [i for i in phis if i < 0.0]: + r = False + return r + + def __str__(self): + print type(self.R) + return "simplex: %d: [%s]" % (len(self.R), ", ".join([str(i) for i in self.R]),) + + + + + class grid(object): facet_re = re.compile(r''' -\s+(?Pf\d+).*? @@ -89,6 +111,20 @@ class grid(object): (R, S) = self.get_simplex_and_nearest_points(X) return run_baker(X, R, S) + def get_points_conn(self, X): + if not self.faces: + self.construct_connectivity() + + (dist, indicies) = self.tree.query(X, 3) + print dist, indicies, [(self.points[i][0], self.points[i][1]) for i in indicies] + print [i.name for i in self.facets_for_point[indicies[0]]] + print self.facets_for_point[indicies[0]] + + smplx = simplex([self.points[i] for i in indicies]) + + if not smplx.contains(X): + raise AssertionError('simplex should contain point') + def construct_connectivity(self): """ a call to this method prepares the internal connectivity structure.