diff --git a/lib/baker/tools.py b/lib/baker/tools.py index dccd75b..2bbae0e 100644 --- a/lib/baker/tools.py +++ b/lib/baker/tools.py @@ -71,12 +71,12 @@ def exact_func_3D(X): def improved_answer(answer, exact, verbose=False): if not answer['error']: return True - smblog.debug('exact: %s' % exact) - smblog.debug('qlin: %s' % answer['qlin']) smblog.debug('error: %s' % answer['error']) + smblog.debug('qlin: %s' % answer['qlin']) smblog.debug('final: %s' % answer['final']) + smblog.debug('exact: %s' % exact) - if abs(answer['final'] - exact) <= abs(answer['qlin'] - exact): + if np.abs(answer['final'] - exact) <= np.abs(answer['qlin'] - exact): smblog.debug(":) improved result") return True else: diff --git a/lib/grid/DD.py b/lib/grid/DD.py index 89f8f7e..9fa0961 100644 --- a/lib/grid/DD.py +++ b/lib/grid/DD.py @@ -66,14 +66,14 @@ class random_grid(rect_grid): smblog.debug("appx_side_res: %d" % appx_side_res) delta = 1.0 / float(appx_side_res) - for x in xrange(appx_side_res): + for x in xrange(appx_side_res + 1): cur_x = x * delta for cur_y in (0, 1): new_point = [cur_x, cur_y] points.append(new_point) q.append(exact_func(new_point)) - for y in xrange(appx_side_res): + for y in xrange(appx_side_res + 1): cur_y = y * delta for cur_x in (0, 1): new_point = [cur_x, cur_y] diff --git a/lib/grid/__init__.py b/lib/grid/__init__.py index 25f5fc6..91d26d0 100644 --- a/lib/grid/__init__.py +++ b/lib/grid/__init__.py @@ -1,6 +1,7 @@ import sys import re from collections import defaultdict +import inspect import numpy as np import scipy.spatial @@ -44,13 +45,18 @@ class grid(object): self.facets_for_point = defaultdict(list) def create_mesh(self, indicies): + """ + this function takes a list of indicies, and then creates + and returns a grid object (collection of points and q). + + note: the input is indicies, the grid contains points + """ p = [self.points[i] for i in indicies] q = [self.q[i] for i in indicies] return grid(p, q) def get_containing_simplex(self, X): if not self.faces: - smblog.debug('setting up connectivity') self.construct_connectivity() # get closest point @@ -67,16 +73,21 @@ class grid(object): facets_to_check = self.facets_for_point[closest_point] attempts = 0 - while not simplex: + while not simplex and facets_to_check: attempts += 1 # if attempts > 20: # raise smberror("probably recursing to many times") - cur_facet = facets_to_check.pop() + cur_facet = facets_to_check.pop(0) checked_facets.append(cur_facet) - facets_to_check.extend([i for i in cur_facet.neighbors if i not in checked_facets]) if cur_facet.contains(X, self): simplex = cur_facet + continue + + new_facest = [] + for neighbor in cur_facet.neighbors: + if (neighbor not in checked_facets) and (neighbor not in facets_to_check): + facets_to_check.append(neighbor) if not simplex: raise AssertionError('no containing simplex found') @@ -96,12 +107,31 @@ class grid(object): S is S_j from baker's paper : some points from all point that are not the simplex """ + smblog.debug(inspect.stack()[1][3]) + smblog.debug("extra points: %d" % extra_points) + smblog.debug("simplex size: %d" % simplex_size) + + r_mesh = self.get_containing_simplex(X) + # smblog.debug("R:\n%s" % r_mesh) + + # and some UNIQUE extra points (dist, indicies) = self.tree.query(X, simplex_size + extra_points) - # r_mesh = self.create_mesh(indicies[:simplex_size]) - r_mesh = self.get_containing_simplex(X) - # and some extra points - s_mesh = self.create_mesh(indicies[simplex_size:]) + unique_indicies = [] + for index in indicies: + if self.points[index] not in r_mesh.points: + unique_indicies.append(index) + + smblog.debug("indicies: %s" % ",".join([str(i) for i in indicies])) + smblog.debug("indicies: %s" % ",".join([str(i) for i in unique_indicies])) + s_mesh = self.create_mesh(unique_indicies)# indicies[simplex_size:]) + + # TODO: eventually remove this test: + for point in s_mesh.points: + if point in r_mesh.points: + smblog.error("ERROR") + smblog.error("\n%s\nin\n%s" % (point, r_mesh)) + raise smberror("repeating point S and R") return (r_mesh, s_mesh) @@ -136,8 +166,8 @@ class grid(object): if not simplex: raise AssertionError('no containing simplex found') - R = self.get_containing_simplex(X)# self.create_mesh(simplex.verts) - + # self.create_mesh(simplex.verts) + R = self.get_containing_simplex(X) s = [] for c,i in enumerate(simplex.neighbors): @@ -169,7 +199,7 @@ class grid(object): this is part of the __init__ for a rect_grid, but can be called from any grid object """ - smblog.debug() + smblog.debug('start') qdelaunay_string = get_qdelaunay_dump_str(self) facet_to_facets = [] for matcher in grid.facet_re.finditer(qdelaunay_string): @@ -194,13 +224,8 @@ class grid(object): if rel[1] in self.faces: self.faces[rel[0]].add_neighbor(self.faces[rel[1]]) - # for matcher in grid.point_re.finditer(qdelaunay_string): - # d = matcher.groupdict() + smblog.debug('end') - # point = d['point'] - # neighboring_facets = d['neigh'] - - # self.facets_for_point[int(point[1:])] = [i for i in neighboring_facets.split() if i in self.faces] def __str__(self): r = '' diff --git a/setup.py b/setup.py index e815ff2..86a87d6 100644 --- a/setup.py +++ b/setup.py @@ -12,4 +12,6 @@ setup( 'numpy', 'scipy', ], + + entry_points = {'console_scripts': ['booze = baker.booze:main',]}, ) diff --git a/tools/mpl_plot.py b/tools/mpl_plot.py new file mode 100644 index 0000000..e90ba45 --- /dev/null +++ b/tools/mpl_plot.py @@ -0,0 +1,14 @@ +#!/usr/bin/env python + +good = [82, 164, 247, 328, 411, 491, 565, 649, 729, 808] +bad = [18, 36, 53, 72, 89, 109, 135, 151, 171, 192] +import pylab + +pylab.xlabel('total runs') +pylab.ylabel('count') +pylab.grid(True) +pylab.plot(bad) +pylab.plot(good) +pylab.legend(('bad', 'good')) + +pylab.show()