sm/smbinterp
1
0
Fork 0
Code Issues Pull Requests Activity

wrapping up a night. there isn't enough consistent improvement to merit using this method. i must have a bug somewhere.

Browse Source
This commit is contained in:
Stephen McQuay 2010-04-29 23:29:35 -06:00
parent 3dcc10ea0e
commit 2cbd92e15b
13 changed files with 115 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
bin/test.py
View File

@ -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)

2
bin/test2d-connectivity.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/python2.6
#!/usr/bin/env python
import sys
from grid.DD import random_grid

27
bin/test_baker.py
View File

@ -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))

14
lib/baker/__init__.py
View File

@ -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]\

50
lib/baker/tools.py
View File

@ -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

15
lib/grid/DD.py
View File

@ -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))

34
lib/grid/__init__.py
View File

@ -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):

11
lib/grid/simplex.py
View File

@ -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

2
test/baker.test.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/python
#!/usr/bin/env python
import unittest
import baker

2
test/baker3D.test.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/python
#!/usr/bin/env python
import unittest
import baker

2
test/cubic.test.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/python
#!/usr/bin/env python
import unittest

2
test/qhull.test.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/python
#!/usr/bin/env python
import unittest

2
test/quad.test.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/python
#!/usr/bin/env python
import unittest