working on the guarantee simplex routine. I have a test case to try it out against, and need to look at visiting neighbors of faces adjacent to points.

bin/test.py

@@ -3,7 +3,7 @@
 import sys
 from grid.DDD import random_grid
 from baker import get_phis_3D, run_baker_3D
-from baker.tools import exact_func_3D, smberror, evaluate_answer 
+from baker.tools import exact_func_3D, smberror, improved_answer 
 try:
   total_points = int(sys.argv[1])
 except:
@@ -32,13 +32,13 @@ if [i for i in phis if i < 0.0]:
 
 try:
   r = run_baker_3D(X, R, S)
-  evaluate_answer(r, exact)
+  improved_answer(r, exact)
 except smberror as e:
   print e
 
 
 try:
   answer =  g.run_baker(X)
-  evaluate_answer(answer, exact)
+  improved_answer(answer, exact)
 except smberror as e:
   print e

bin/test2d-connectivity.py

@@ -0,0 +1,42 @@
+#!/usr/bin/python2.6
+
+import sys
+from grid.DD import random_grid
+from baker import get_phis, run_baker
+from baker.tools import exact_func, smberror, improved_answer
+from glob import glob
+from os import remove
+
+FILE_PREFIX='/tmp/qhull-'
+
+for g in glob('%s*' % FILE_PREFIX):
+  remove(g)
+
+try:
+  total_points = int(sys.argv[1])
+  total_tries  = int(sys.argv[2])
+except:
+  print "usage: app.py total_points total_tries"
+  sys.exit(1)
+
+
+print "total points", total_points
+
+
+
+d = {True: 0, False: 0}
+for i in xrange(total_tries):
+  g = random_grid(total_points)
+
+  open('%s%0.3d.txt' % (FILE_PREFIX, i), 'w').write(g.for_qhull())
+
+  X = [0.3, 0.3]
+  exact = exact_func(X)
+
+  try:
+    answer =  g.run_baker(X)
+    d[improved_answer(answer, exact)] += 1
+  except smberror as e:
+    print e
+
+print d

bin/test2d-qhull-specific.py

@@ -0,0 +1,36 @@
+#!/usr/bin/python2.6
+
+import sys
+from grid.DD import grid
+from grid.qhull import parse_qhull_file
+from baker import get_phis, run_baker
+from baker.tools import exact_func, smberror, improved_answer
+from glob import glob
+from os import remove
+
+FILE_PREFIX='/tmp/qhull-'
+
+for g in glob('%s*' % FILE_PREFIX):
+  remove(g)
+
+try:
+  verts = parse_qhull_file(sys.argv[1], verbose=True)
+except:
+  print "usage: app.py qhullinputfile"
+  sys.exit(1)
+
+
+q = [exact_func(i) for i in verts]
+
+g = grid(verts,q)
+
+X = [0.3, 0.3]
+
+g.get_containing_simplex(X)
+exact = exact_func(X)
+
+try:
+  answer =  g.run_baker(X)
+  d[improved_answer(answer, exact)] += 1
+except smberror as e:
+  print e

bin/test_baker.py

@@ -0,0 +1,49 @@
+#!/usr/bin/python
+
+import math
+from baker        import run_baker
+
+from grid.DD      import grid
+from grid.simplex import contains
+
+def exact_func(X):
+  x = X[0]
+  y = X[0]
+  return 1 - math.sin((x-0.5)**2 + (y-0.5)**2)
+
+if __name__ == '__main__':
+  points = [
+             [ 0.25, 0.40], # 0
+             [ 0.60, 0.80], # 1
+             [ 0.65, 0.28], # 2
+             [ 0.28, 0.65], # 3
+             [ 1.00, 0.75], # 4
+             [ 0.30, 0.95], # 5
+             [ 0.80, 0.50], # 6
+             [ 0.35, 0.15], # 7
+           ]
+  q = [exact_func(p) for p in points]
+
+  X = [0.55, 0.45]
+
+  g = grid(points, q)
+  g.construct_connectivity()
+  R = g.create_mesh(range(3))
+
+  print contains(X, R.points)
+
+  extra = 3
+  S = g.create_mesh(range(3, 3 + extra))
+
+  answer          = run_baker(X, R, S, 3)
+  answer['exact'] = exact_func(X)
+  t = ('qlin', 'error', 'exact')
+
+  for i in t:
+    print "%s %s" % (
+                      i.ljust(10),
+                      ("%0.4f" % answer[i]),
+                    )
+  lin_err   = abs(answer['exact'] - answer['qlin'])
+  final_err = abs(answer['exact'] - answer['final'])
+  print lin_err >= final_err

lib/baker/__init__.py

@@ -198,14 +198,6 @@ def run_baker(X, R, S, order=2):
   # calculate values only for the simplex triangle
   phi, qlin = qlinear(X, R)
 
-  if [i for i in phi if i <= 0.0]:
-    s = "this is not a containing simplex:\n"
-    s += "  X: %s\n" % X
-    s += "  R: %s\n" % R
-    s += "  phi: %s, sum(%0.4e)\n" % (phi, sum(phi))
-    print >> sys.stderr, s
-    raise smberror("simplex does not contain point")
-
   if len(S.points) == 0:
     answer = {
                'a': None,

lib/baker/tools.py

@@ -36,12 +36,15 @@ def exact_func_3D(X):
   z = X[2]
   return np.power((np.sin(x * np.pi / 2.0) * np.sin(y * np.pi / 2.0) * np.sin(z * np.pi / 2.0)), 2)
 
-def evaluate_answer(answer, exact):
+def improved_answer(answer, exact, verbose=False):
+  if verbose:
     print 'qlin' , answer['qlin']
     print 'error', answer['error']
     print 'final', answer['final']
 
   if abs(answer['final'] - exact) <= abs(answer['qlin'] - exact):
-    print ":) improved result"
+    if verbose: print ":) improved result"
+    return True
   else:
-    print ":( damaged result"
+    if verbose: print ":( damaged result"
+    return False

lib/grid/__init__.py

@@ -7,7 +7,7 @@ import scipy.spatial
 
 from baker        import run_baker
 from baker.tools  import exact_func, smberror
-from simplex      import face
+from simplex      import face, contains
 from smcqdelaunay import *
 
 
@@ -49,6 +49,9 @@ class grid(object):
     q = [self.q[i]      for i in indicies]
     return grid(p, q)
 
+  def get_containing_simplex(self, X):
+    pass
+
 
   def get_simplex_and_nearest_points(self, X, extra_points = 3, simplex_size = 3):
     """
@@ -115,6 +118,8 @@ class grid(object):
 
     try:
       (R, S) = self.get_simplex_and_nearest_points(X)
+      if not contains(X, R.points):
+        raise smberror("run_baker with get_simplex_and_nearest_points returned non-containing simplex")
       answer = run_baker(X, R, S)
     except smberror, e:
       print >> sys.stderr, "caught error: %s, trying with connectivity-based mesh" % e

lib/grid/qhull.py

@@ -0,0 +1,15 @@
+
+def parse_qhull_file(filename, verbose=False):
+  f = open(filename, 'r')
+
+  if verbose:
+    print 'filename: ', filename
+    degree = int(f.readline().strip())
+    print "degree:", degree
+    print "number of points", f.readline().strip()
+
+  verts = []
+  for p in f:
+    v = [float(i) for i in p.strip().split()]
+    verts.append(v)
+  return verts

tools/blender/plot_cloud.py

@@ -1,7 +1,7 @@
 from Blender import *
 import bpy
 
-from grid.test3d import rect_grid
+from grid.DDD import rect_grid
 
 g = rect_grid(10, 10, 20)
 

tools/blender/qhull-plt.py

@@ -0,0 +1,31 @@
+"""
+Name: 'qhull mesh display'
+Blender: 249
+Group: 'Add'
+Tooltip:'this lets you display a qhull mesh'
+"""
+from Blender import *
+import bpy
+
+from grid.DDD import rect_grid
+
+def draw_file(filename):
+  f = open(filename, 'r')
+  print 'filename: ', filename
+  print "degree:", f.readline().strip()
+  print "number of points", f.readline().strip()
+
+  verts = []
+  for p in f:
+    v = [float(i) for i in p.strip().split()]
+    if len(v) == 2:
+      v.append(0.0)
+    verts.append(v)
+
+  me = bpy.data.meshes.new('points')
+  me.verts.extend(verts)
+  # me.faces.extend([i.verts for i in p.faces.itervalues()])
+  scn = bpy.data.scenes.active
+  ob  = scn.objects.new(me, 'points_obj')
+
+Window.FileSelector(draw_file, 'input file')