From 2cbd92e15b1cfb3218306ad6f9853e9314283459 Mon Sep 17 00:00:00 2001
From: Stephen Mardson McQuay <stephen@mcquay.me>
Date: Thu, 29 Apr 2010 23:29:35 -0600
Subject: [PATCH] wrapping up a night. there isn't enough consistent
 improvement to merit using this method. i must have a bug somewhere.

---
 bin/test.py                | 30 +++++++++++------------
 bin/test2d-connectivity.py |  2 +-
 bin/test_baker.py          | 27 ++++++++++----------
 lib/baker/__init__.py      | 14 +++++------
 lib/baker/tools.py         | 50 ++++++++++++++++++++++++++++----------
 lib/grid/DD.py             | 15 +++++++++---
 lib/grid/__init__.py       | 34 ++++++++++++--------------
 lib/grid/simplex.py        | 11 +++++++--
 test/baker.test.py         |  2 +-
 test/baker3D.test.py       |  2 +-
 test/cubic.test.py         |  2 +-
 test/qhull.test.py         |  2 +-
 test/quad.test.py          |  2 +-
 13 files changed, 115 insertions(+), 78 deletions(-)

diff --git a/bin/test.py b/bin/test.py
index e05aa47..0bb926f 100755
--- a/bin/test.py
+++ b/bin/test.py
@@ -1,45 +1,45 @@
-#!/usr/bin/python2.6
+#!/usr/bin/env python
 
 import sys
 from grid.DDD import random_grid
 from baker import get_phis_3D, run_baker_3D
-from baker.tools import logging as log
-from baker.tools import exact_func_3D, smberror, improved_answer 
+from baker.tools import smblog
+from baker.tools import exact_func_3D, smberror, improved_answer
+from grid.simplex import contains
+
 try:
   total_points = int(sys.argv[1])
 except:
   total_points = 20
 
 
-log.info(total_points)
+smblog.info(total_points)
 g = random_grid(total_points)
 
 open('/tmp/for_qhull.txt', 'w').write(g.for_qhull())
 
-X = [0.3, 0.3, 0.4]
+X = [0.5, 0.3, 0.4]
 
 R, S =  g.get_simplex_and_nearest_points(X, extra_points = 6, simplex_size=4)
-
-print "R\n", R
-print "S\n", S
+smblog.info("Containing Simplex: %s" % contains(X, R.points))
 
 exact = exact_func_3D(X)
-print "exact solution: %0.6f" % exact
+smblog.info("exact solution: %0.6f" % exact)
 
 phis = get_phis_3D(X, R.points)
-print "phi values (should all be positive): ", phis, sum(phis)
+smblog.info("phi values (should all be positive): %s %s" % (phis, sum(phis)))
 if [i for i in phis if i < 0.0]:
-  print "problems"
+  smblog.error("problems")
 
 try:
   r = run_baker_3D(X, R, S)
-  improved_answer(r, exact)
+  smblog.info("run_baker manual: %s" % improved_answer(r, exact))
 except smberror as e:
-  print e
+  smblog.error(e)
 
 
 try:
   answer =  g.run_baker(X)
-  improved_answer(answer, exact)
+  smblog.info("run_baker from grid: %s" % improved_answer(answer, exact))
 except smberror as e:
-  print e
+  smblog.error(e)
diff --git a/bin/test2d-connectivity.py b/bin/test2d-connectivity.py
index 3ffbfb8..5c07ac1 100755
--- a/bin/test2d-connectivity.py
+++ b/bin/test2d-connectivity.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python2.6
+#!/usr/bin/env python
 
 import sys
 from grid.DD import random_grid
diff --git a/bin/test_baker.py b/bin/test_baker.py
index bd59e44..a6d1eb4 100755
--- a/bin/test_baker.py
+++ b/bin/test_baker.py
@@ -1,7 +1,10 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import math
 from baker        import run_baker
+from baker.tools import improved_answer
+
+from baker.tools import smblog
 
 from grid.DD      import grid
 from grid.simplex import contains
@@ -32,18 +35,16 @@ if __name__ == '__main__':
 
   print contains(X, R.points)
 
-  extra = 3
+  try:
+    extra = int(sys.argv[1])
+  except:
+    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 xrange(1,4):
+    answer  = run_baker(X, R, S, i)
+    exact   = exact_func(X)
+    smblog.debug(improved_answer(answer, 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
+  smblog.debug(improved_answer(g.run_baker(X), exact))
diff --git a/lib/baker/__init__.py b/lib/baker/__init__.py
index 48b092b..bcdea62 100644
--- a/lib/baker/__init__.py
+++ b/lib/baker/__init__.py
@@ -1,5 +1,5 @@
 from baker import *
-from baker.tools import logging as log
+from baker.tools import smblog
 import numpy as np
 import sys
 
@@ -30,8 +30,8 @@ def get_phis(X, R):
   try:
     phi = np.linalg.solve(A,b)
   except np.linalg.LinAlgError as e:
-    msg = "warning: get_phis: calculation of phis yielded a linearly dependant system (%s)" % e
-    log.error(msg)
+    msg = "calculation of phis yielded a linearly dependant system (%s)" % e
+    smblog.error(msg)
     raise smberror(msg)
     phi = np.dot(np.linalg.pinv(A), b)
 
@@ -69,7 +69,7 @@ def get_phis_3D(X, R):
   try:
     phi = np.linalg.solve(A,b)
   except np.linalg.LinAlgError as e:
-    print >> sys.stderr, "warning: get_phis_3D: calculation of phis yielded a linearly dependant system", e
+    smblog.error("calculation of phis yielded a linearly dependant system: %s" % e)
     phi = np.dot(np.linalg.pinv(A), b)
 
   return phi
@@ -130,7 +130,7 @@ def get_error_quadratic(phi, R, S):
   try:
     (a, b, c) = np.linalg.solve(A,b)
   except np.linalg.LinAlgError as e:
-    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e
+    smblog.error("linear calculation went bad, resorting to np.linalg.pinv: %s" % e)
     (a, b, c) = np.dot(np.linalg.pinv(A), b)
 
   error_term =  a * phi[0] * phi[1]\
@@ -171,7 +171,7 @@ def get_error_cubic(phi, R, S):
   try:
     (a, b, c, d, e, f, g) = np.linalg.solve(A,b)
   except np.linalg.LinAlgError as e:
-    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e
+    smblog.error("linear calculation went bad, resorting to np.linalg.pinv: %s" % e)
     (a, b, c, d, e, f, g) = np.dot(np.linalg.pinv(A), b)
 
   error_term =  a * phi[0] * phi[1] * phi[1]\
@@ -280,7 +280,7 @@ def run_baker_3D(X, R, S):
   try:
     (a, b, c, d, e, f) = np.linalg.solve(A,b)
   except np.linalg.LinAlgError as e:
-    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e
+    smblog.error("linear calculation went bad, resorting to np.linalg.pinv: %s", e)
     (a, b, c, d, e, f) = np.dot(np.linalg.pinv(A), b)
 
   error_term =  a * phi[0] * phi[1]\
diff --git a/lib/baker/tools.py b/lib/baker/tools.py
index 0a8b7fb..dccd75b 100644
--- a/lib/baker/tools.py
+++ b/lib/baker/tools.py
@@ -1,15 +1,37 @@
 import os
 import logging
-
-logging.basicConfig(
-                     level = logging.DEBUG,
-                     format = '%(asctime)s %(levelname)s %(message)s',
-                     filename = os.path.join(os.sep, 'tmp', 'baker.lol'),
-                   )
-
+import inspect
 import numpy as np
 
 
+class smbLog(object):
+  interpolator = "%s ==> %s"
+  def __init__(self, level = logging.DEBUG):
+    logging.basicConfig(
+                         level = level,
+                         format = '%(asctime)s %(levelname)s %(message)s',
+                         filename = os.path.join(os.sep, 'tmp', 'baker.lol'),
+                       )
+    self.log = logging.getLogger()
+  def debug(self, message = None):
+    msg = smbLog.interpolator % (inspect.stack()[1][3], message)
+    self.log.debug(msg)
+
+  def info(self, message = None):
+    msg = smbLog.interpolator % (inspect.stack()[1][3], message)
+    self.log.info(msg)
+
+  def warn(self, message = None):
+    msg = smbLog.interpolator % (inspect.stack()[1][3], message)
+    self.log.warn(msg)
+
+  def error(self, message = None):
+    msg = smbLog.interpolator % (inspect.stack()[1][3], message)
+    self.log.error(msg)
+
+
+smblog = smbLog(logging.DEBUG)
+
 class smberror(Exception):
   """
     this is a silly little exception subclass
@@ -47,14 +69,16 @@ def exact_func_3D(X):
   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 improved_answer(answer, exact, verbose=False):
-  if verbose:
-    print 'qlin' , answer['qlin']
-    print 'error', answer['error']
-    print 'final', answer['final']
+  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('final: %s' % answer['final'])
 
   if abs(answer['final'] - exact) <= abs(answer['qlin'] - exact):
-    if verbose: print ":) improved result"
+    smblog.debug(":) improved result")
     return True
   else:
-    if verbose: print ":( damaged result"
+    smblog.debug(":( damaged result")
     return False
diff --git a/lib/grid/DD.py b/lib/grid/DD.py
index fe73f74..89f8f7e 100644
--- a/lib/grid/DD.py
+++ b/lib/grid/DD.py
@@ -1,6 +1,6 @@
 from grid import grid as basegrid
 
-from baker.tools import exact_func
+from baker.tools import exact_func, smblog
 
 import numpy as np
 
@@ -57,16 +57,25 @@ class rect_grid(grid):
 
 class random_grid(rect_grid):
   def __init__(self, num_points = 10):
+    smblog.debug("number of points: %d" % num_points)
     points = []
     q = []
 
     r = np.random
     appx_side_res = int(np.sqrt(num_points))
+    smblog.debug("appx_side_res: %d" % appx_side_res)
     delta = 1.0 / float(appx_side_res)
+
     for x in xrange(appx_side_res):
       cur_x = x * delta
-      for y in xrange(appx_side_res):
-        cur_y = y * 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):
+      cur_y = y * delta
+      for cur_x in (0, 1):
         new_point = [cur_x, cur_y]
         points.append(new_point)
         q.append(exact_func(new_point))
diff --git a/lib/grid/__init__.py b/lib/grid/__init__.py
index 1400858..25f5fc6 100644
--- a/lib/grid/__init__.py
+++ b/lib/grid/__init__.py
@@ -6,7 +6,7 @@ import numpy as np
 import scipy.spatial
 
 from baker        import run_baker
-from baker.tools  import exact_func, smberror, logging
+from baker.tools  import exact_func, smberror, smblog
 from simplex      import face, contains
 from smcqdelaunay import *
 
@@ -42,7 +42,6 @@ class grid(object):
     self.tree             = scipy.spatial.KDTree(self.points)
     self.faces            = {}
     self.facets_for_point = defaultdict(list)
-    self.counter = 0
 
   def create_mesh(self, indicies):
     p = [self.points[i] for i in indicies]
@@ -51,19 +50,17 @@ class grid(object):
 
   def get_containing_simplex(self, X):
     if not self.faces:
-      logging.debug('get_containing_simplex: setting up connectivity')
+      smblog.debug('setting up connectivity')
       self.construct_connectivity()
 
     # get closest point
     (dist, indicies) = self.tree.query(X, 2)
     closest_point = indicies[0]
 
-    logging.debug('counter: %d' % self.counter)
-    self.counter += 1
-    logging.debug('X: %s' % X)
-    logging.debug('point index: %d' % closest_point)
-    logging.debug('actual point %s' % self.points[closest_point])
-    logging.debug('distance = %0.4f' % dist[0])
+    smblog.debug('X: %s' % X)
+    smblog.debug('point index: %d' % closest_point)
+    smblog.debug('actual point %s' % self.points[closest_point])
+    smblog.debug('distance = %0.4f' % dist[0])
 
     simplex = None
     checked_facets = []
@@ -71,10 +68,9 @@ class grid(object):
 
     attempts = 0
     while not simplex:
-      logging.debug('attempt: %d' % attempts)
       attempts += 1
-      if attempts > 10:
-        raise smberror("probably recursing to many times")
+      # if attempts > 20:
+      #   raise smberror("probably recursing to many times")
       cur_facet = facets_to_check.pop()
       checked_facets.append(cur_facet)
       facets_to_check.extend([i for i in cur_facet.neighbors if i not in checked_facets])
@@ -87,7 +83,7 @@ class grid(object):
 
     R = self.create_mesh(simplex.verts)
 
-    logging.debug('this must be R: %s' % R)
+    smblog.debug('total attempts before finding simplex: %d' % attempts)
     return R
 
 
@@ -140,7 +136,7 @@ class grid(object):
     if not simplex:
       raise AssertionError('no containing simplex found')
 
-    R = get_containing_simplex(X)# self.create_mesh(simplex.verts)
+    R = self.get_containing_simplex(X)# self.create_mesh(simplex.verts)
 
 
     s = []
@@ -150,18 +146,18 @@ class grid(object):
 
     return R, S
 
-  def run_baker(self, X):
+  def run_baker(self, X, extra_points = 3, order = 2):
     answer = None
 
     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)
+      answer = run_baker(X, R, S, order)
     except smberror, e:
-      print >> sys.stderr, "caught error: %s, trying with connectivity-based mesh" % e
+      smblog.error("caught error: %s, trying with connectivity-based mesh" % e)
       (R, S) = self.get_points_conn(X)
-      answer = run_baker(X, R, S)
+      answer = run_baker(X, R, S, order)
 
     return answer
 
@@ -173,7 +169,7 @@ class grid(object):
 
       this is part of the __init__ for a rect_grid, but can be called from any grid object
     """
-    logging.debug('calling construct_connectivity')
+    smblog.debug()
     qdelaunay_string = get_qdelaunay_dump_str(self)
     facet_to_facets = []
     for matcher in grid.facet_re.finditer(qdelaunay_string):
diff --git a/lib/grid/simplex.py b/lib/grid/simplex.py
index 40f60c4..d91fea8 100644
--- a/lib/grid/simplex.py
+++ b/lib/grid/simplex.py
@@ -1,4 +1,5 @@
-from baker import get_phis
+from baker import get_phis, get_phis_3D
+from baker.tools import smblog
 
 TOL = 1e-8
 
@@ -6,8 +7,14 @@ def contains(X, R):
   """
     tests if X (point) is in R (a simplex,
     represented by a list of n-degree coordinates)
+
+    it now correctly checks for 2/3-D points
   """
-  phis = get_phis(X, R)
+  if len(X) == 2:
+    phis = get_phis(X, R)
+  elif len(X) == 3:
+    phis = get_phis_3D(X, R)
+
   r = True
   if [i for i in phis if i < 0.0 - TOL]:
     r = False
diff --git a/test/baker.test.py b/test/baker.test.py
index 3c1382d..df4b4cb 100755
--- a/test/baker.test.py
+++ b/test/baker.test.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import unittest
 import baker
diff --git a/test/baker3D.test.py b/test/baker3D.test.py
index bf2e6d7..49a29cd 100755
--- a/test/baker3D.test.py
+++ b/test/baker3D.test.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import unittest
 import baker
diff --git a/test/cubic.test.py b/test/cubic.test.py
index 98568b2..b6ae673 100755
--- a/test/cubic.test.py
+++ b/test/cubic.test.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import unittest
 
diff --git a/test/qhull.test.py b/test/qhull.test.py
index fdb02ca..48bd45c 100755
--- a/test/qhull.test.py
+++ b/test/qhull.test.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import unittest
 
diff --git a/test/quad.test.py b/test/quad.test.py
index ea8f67e..8de48f7 100755
--- a/test/quad.test.py
+++ b/test/quad.test.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import unittest