From 837a72b2460990f87924dbb99ebd81daf15fcd9e Mon Sep 17 00:00:00 2001 From: "Stephen M. McQuay" Date: Sat, 17 Sep 2011 15:38:49 -0600 Subject: [PATCH] Major: made scripts pass pep8 and pyflakes --- interp/__init__.py | 39 ---- interp/baker/__init__.py | 350 ++++++++++++++++----------------- interp/bootstrap.py | 7 +- interp/cluster/__init__.py | 38 ++-- interp/config.py | 19 ++ interp/grid/DD.py | 5 +- interp/grid/DDD.py | 8 +- interp/grid/__init__.py | 390 +++++++++++++++++++------------------ interp/grid/gmsh.py | 5 +- interp/tools.py | 111 +++++------ test/all.py | 20 +- test/baker2d.py | 216 ++++++++++---------- test/baker2dorder.py | 141 +++++++------- test/baker3d.py | 47 ++--- test/cubic2d.py | 117 +++++------ test/pattern.py | 68 ++++--- test/quadratic2d.py | 120 ++++++------ 17 files changed, 835 insertions(+), 866 deletions(-) create mode 100644 interp/config.py diff --git a/interp/__init__.py b/interp/__init__.py index 2bf2848..b650ceb 100644 --- a/interp/__init__.py +++ b/interp/__init__.py @@ -1,40 +1 @@ -import os - -import logging -import logging.handlers - -import json - - -LEVELS = {'debug': logging.DEBUG, - 'info': logging.INFO, - 'warning': logging.WARNING, - 'error': logging.ERROR, - 'critical': logging.CRITICAL} - -default_config = { - 'filename': '/tmp/interp.log', - 'level': 'debug', - 'size' : 102400, - 'logbackup': 10, - 'pypath': None, -} - -try: - with open(os.path.expanduser('~/.config/interp.json')) as config_file: - d = json.load(config_file) -except IOError as e: - d = {} - -config = dict(default_config.items() + d.items()) - -logger = logging.getLogger('interp') -logger.setLevel(LEVELS[config['level']]) -my_format = logging.Formatter('%(asctime)s %(levelname)s (%(process)d) %(filename)s %(funcName)s:%(lineno)d %(message)s') -handler = logging.handlers.RotatingFileHandler( - config['filename'], maxBytes = config['size'] * 1024, backupCount = config['logbackup']) -handler.setFormatter(my_format) -logger.addHandler(handler) - - __version__ = '0.2' diff --git a/interp/baker/__init__.py b/interp/baker/__init__.py index 0b9ab68..ddd9e4e 100644 --- a/interp/baker/__init__.py +++ b/interp/baker/__init__.py @@ -1,220 +1,204 @@ -import sys - import numpy as np from functools import wraps import itertools import interp -import logging -log = logging.getLogger('interp') + +AGGRESSIVE_ERROR_SOLVE = True +RAISE_PATHOLOGICAL_EXCEPTION = False + +__version__ = interp.__version__ + def get_phis(X, R): - """ - The get_phis function is used to get barycentric coordonites for a - point on a triangle or tetrahedron. This is equation (*\ref{eq:qlinarea}*) + """ + The get_phis function is used to get barycentric coordonites for a + point on a triangle or tetrahedron (Equation (*\ref{eq:qlinarea}*)) - in 2D: + in 2D: - X - the destination point (2D) - X = [0,0] - R - the three points that make up the 2-D triangular simplex - R = [[-1, -1], [0, 2], [1, -1]] + X - the destination point (2D) + X = [0,0] + R - the three points that make up the 2-D triangular simplex + R = [[-1, -1], [0, 2], [1, -1]] - this will return [0.333, 0.333, 0.333] + this will return [0.333, 0.333, 0.333] - in 3D: + in 3D: - X - the destination point (3D) - X = [0,0,0] - R - the four points that make up the 3-D simplex (tetrahedron) - R = [ - [ 0.0000, 0.0000, 1.0000], - [ 0.9428, 0.0000, -0.3333], - [-0.4714, 0.8165, -0.3333], - [-0.4714, -0.8165, -0.3333], - ] + X - the destination point (3D) + X = [0,0,0] + R - the four points that make up the 3-D simplex (tetrahedron) + R = [ + [ 0.0000, 0.0000, 1.0000], + [ 0.9428, 0.0000, -0.3333], + [-0.4714, 0.8165, -0.3333], + [-0.4714, -0.8165, -0.3333], + ] - this will return [0.25, 0.25, 0.25, 0.25] - """ + this will return [0.25, 0.25, 0.25, 0.25] + """ - # equations (*\ref{eq:lin3d}*) and (*\ref{eq:lin2d}*) - if len(X) == 2: - log.debug("running 2D") - A = np.array([ - [ 1, 1, 1], - [R[0][0], R[1][0], R[2][0]], - [R[0][1], R[1][1], R[2][1]], - ]) - b = np.array([ 1, - X[0], - X[1] - ]) - elif len(X) == 3: - log.debug("running 3D") - A = np.array([ - [ 1, 1, 1, 1 ], - [R[0][0], R[1][0], R[2][0], R[3][0]], - [R[0][1], R[1][1], R[2][1], R[3][1]], - [R[0][2], R[1][2], R[2][2], R[3][2]], - ]) - b = np.array([ 1, - X[0], - X[1], - X[2] - ]) - else: - raise Exception("inapropriate demension on X") - - try: - phi = np.linalg.solve(A,b) - except np.linalg.LinAlgError as e: - msg = "calculation of phis yielded a linearly dependant system (%s)" % e - log.error(msg) - # raise Exception(msg) - phi = np.dot(np.linalg.pinv(A), b) - - log.debug("phi: %s", phi) - - return phi - -def qlinear(X, R): - """ - this calculates the linear portion of q from R to X - - This is equation (*\ref{eq:qlinbasis}*) - - X = destination point - R = a inter.grid object; must have R.points and R.q - """ - - phis = get_phis(X, R.verts) - qlin = np.sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)]) - - log.debug("phis: %s", phis) - log.debug("qlin: %s", qlin) - - return phis, qlin - -def get_error(phi, R, S, order = 2): - """ - Calculate the error approximation terms, returning the unknowns - a,b, and c in equation (*\ref{eq:quadratic2d}*). - """ - B = [] # equation ((*\ref{eq:B2d}*) - w = [] # equation ((*\ref{eq:w}*) - - cur_pattern = pattern(len(phi), order) - log.info("pattern: %s" % cur_pattern) - - for (s,q) in zip(S.verts, S.q): - cur_phi, cur_qlin = qlinear(s, R) - l = [] - for i in cur_pattern: - cur_sum = cur_phi[i[0]] - for j in i[1:]: - cur_sum *= cur_phi[j] - l.append(cur_sum) - - B.append(l) - w.append(q - cur_qlin) - - log.info("B: %s" % B) - log.info("w: %s" % w) + # equations (*\ref{eq:lin3d}*) and (*\ref{eq:lin2d}*) + if len(X) == 2: + A = np.array([ + [1, 1, 1], + [R[0][0], R[1][0], R[2][0]], + [R[0][1], R[1][1], R[2][1]], + ]) + b = np.array([1, X[0], X[1]]) + elif len(X) == 3: + A = np.array([ + [1, 1, 1, 1], + [R[0][0], R[1][0], R[2][0], R[3][0]], + [R[0][1], R[1][1], R[2][1], R[3][1]], + [R[0][2], R[1][2], R[2][2], R[3][2]], + ]) + b = np.array([1, X[0], X[1], X[2]]) + else: + raise Exception("inapropriate demension on X") + phi = np.linalg.solve(A, b) + return phi - B = np.array(B) - w = np.array(w) +def qlinear(X, R, q): + """ + this calculates the linear portion of q from R to X - A = np.dot(B.T, B) - b = np.dot(B.T, w) + This is equation (*\ref{eq:qlinbasis}*) - try: - abc = np.linalg.solve(A,b) - except np.linalg.LinAlgError as e: - log.error("linear calculation went bad, resorting to np.linalg.pinv: %s" % e) - abc = np.dot(np.linalg.pinv(A), b) + X = destination point + R = a inter.grid object; must have R.points and R.q + """ - error_term = 0.0 - for (a, i) in zip(abc, cur_pattern): - cur_sum = a - for j in i: - cur_sum *= phi[j] - error_term += cur_sum + phis = get_phis(X, R) + qlin = np.sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) - log.debug("error_term: %s" % error_term) - return error_term, abc + return phis, qlin -def run_baker(X, R, S, order=2): - """ - This is the main function to call to get an interpolation to X from the - input meshes - X -- the destination point +def get_error(phi, R, R_q, S, S_q, order=2): + """ + Calculate the error approximation terms, returning the unknowns + a,b, and c in equation (*\ref{eq:quadratic2d}*). + """ + B = [] # equation ((*\ref{eq:B2d}*) + w = [] # equation ((*\ref{eq:w}*) - R = Simplex - S = extra points - """ - log.debug("order = %d" % order) - log.debug("extra points = %d" % len(S.verts)) + cur_pattern = pattern(len(phi), order) - answer = { - 'qlin': None, - 'error': None, - 'final': None, - } - # calculate values only for the simplex triangle - phi, qlin = qlinear(X, R) + for (s, cur_q) in zip(S, S_q): + cur_phi, cur_qlin = qlinear(s, R, R_q) + l = [] + for i in cur_pattern: + cur_sum = cur_phi[i[0]] + for j in i[1:]: + cur_sum *= cur_phi[j] + l.append(cur_sum) + + B.append(l) + w.append(cur_q - cur_qlin) + + B = np.array(B) + w = np.array(w) + + A = np.dot(B.T, B) + b = np.dot(B.T, w) + + try: + abc = np.linalg.solve(A, b) + except np.linalg.LinAlgError: + if not AGGRESSIVE_ERROR_SOLVE: + return None, None + abc = np.dot(np.linalg.pinv(A), b) + + error_term = 0.0 + for (a, i) in zip(abc, cur_pattern): + cur_sum = a + for j in i: + cur_sum *= phi[j] + error_term += cur_sum + + return error_term, abc + + +def run_baker(X, R, R_q, S, S_q, order=2): + """ + This is the main function to call to get an interpolation to X from the + input meshes + + X -- the destination point + + R = Simplex + S = extra points + """ + + answer = { + 'qlin': None, + 'error': None, + 'final': None, + } + + # calculate values only for the simplex triangle + phi, qlin = qlinear(X, R, R_q) + + if order == 1: + answer['qlin'] = qlin + answer['final'] = qlin + return answer + elif order in xrange(2, 11): + error_term, abc = get_error(phi, R, R_q, S, S_q, order) + + # if a pathological vertex configuration was encountered and + # AGGRESSIVE_ERROR_SOLVE is False, get_error will return (None, None) + # indicating that only linear interpolation should be performed + if (error_term is None) and (abc is None): + if RAISE_PATHOLOGICAL_EXCEPTION: + raise np.linalg.LinAlgError("Pathological Vertex Config") + answer['qlin'] = qlin + answer['final'] = qlin + return answer + else: + raise Exception('unsupported order "%d" for baker method' % order) + + q_final = qlin + error_term - if order == 1: answer['qlin'] = qlin - answer['final'] = qlin + answer['error'] = error_term + answer['final'] = q_final + answer['abc'] = abc + return answer - elif order in xrange(2,11): - error_term, abc = get_error(phi, R, S, order) - else: - raise Exception('unsupported order "%d" for baker method' % order) - - q_final = qlin + error_term - - answer['qlin' ] = qlin - answer['error'] = error_term - answer['final'] = q_final - answer['abc' ] = abc - - log.debug(answer) - - return answer def memoize(f): - """ - for more information on what I'm doing here, please read: - http://en.wikipedia.org/wiki/Memoize - """ - cache = {} - @wraps(f) - def memf(simplex_size, nu): - x = (simplex_size, nu) - if x not in cache: - log.debug("adding to cache: %s", x) - cache[x] = f(simplex_size, nu) - return cache[x] - return memf + """ + for more information on what I'm doing here, please read: + http://en.wikipedia.org/wiki/Memoize + """ + cache = {} + + @wraps(f) + def memf(simplex_size, nu): + x = (simplex_size, nu) + if x not in cache: + cache[x] = f(simplex_size, nu) + return cache[x] + return memf @memoize def pattern(simplex_size, nu): - """ - This function returns the pattern requisite to compose the error - approximation function, and the matrix B. - """ - log.debug("pattern: simplex: %d, order: %d" % (simplex_size, nu)) + """ + This function returns the pattern requisite to compose the error + approximation function, and the matrix B. + """ - r = [] - for i in itertools.product(xrange(simplex_size), repeat = nu): - if len(set(i)) !=1: - r.append(tuple(sorted(i))) - unique_r = list(set(r)) - return unique_r + r = [] + for i in itertools.product(xrange(simplex_size), repeat=nu): + if len(set(i)) != 1: + r.append(tuple(sorted(i))) + unique_r = list(set(r)) + return unique_r diff --git a/interp/bootstrap.py b/interp/bootstrap.py index 7742642..a9868e6 100644 --- a/interp/bootstrap.py +++ b/interp/bootstrap.py @@ -5,12 +5,13 @@ import rlcompleter historyPath = os.path.expanduser("~/.pyhistory") + def save_history(historyPath=historyPath): - import readline - readline.write_history_file(historyPath) + import readline + readline.write_history_file(historyPath) if os.path.exists(historyPath): - readline.read_history_file(historyPath) + readline.read_history_file(historyPath) atexit.register(save_history) del os, atexit, readline, rlcompleter, save_history, historyPath diff --git a/interp/cluster/__init__.py b/interp/cluster/__init__.py index 5203f3a..9fb1eeb 100644 --- a/interp/cluster/__init__.py +++ b/interp/cluster/__init__.py @@ -1,28 +1,30 @@ from multiprocessing.managers import BaseManager import Queue -tasks_q = Queue.Queue() +tasks_q = Queue.Queue() results_q = Queue.Queue() minions_q = Queue.Queue() -master_q = Queue.Queue() +master_q = Queue.Queue() + class QueueManager(BaseManager): - """ - One QueueManager to rule all network Queues - """ - pass + """ + One QueueManager to rule all network Queues + """ + pass + +QueueManager.register('get_tasks_q', callable=lambda: tasks_q) +QueueManager.register('get_results_q', callable=lambda: results_q) +QueueManager.register('get_minions_q', callable=lambda: minions_q) +QueueManager.register('get_master_q', callable=lambda: master_q) -QueueManager.register('get_tasks_q' , callable=lambda:tasks_q ) -QueueManager.register('get_results_q', callable=lambda:results_q ) -QueueManager.register('get_minions_q', callable=lambda:minions_q ) -QueueManager.register('get_master_q' , callable=lambda:master_q ) def get_qs(qm): - """ - pass in a QueueManager, and this function returns all relevant - queues attached to that QueueManager. - """ - return (qm.get_tasks_q(), - qm.get_results_q(), - qm.get_master_q(), - qm.get_minions_q()) + """ + pass in a QueueManager, and this function returns all relevant + queues attached to that QueueManager. + """ + return (qm.get_tasks_q(), + qm.get_results_q(), + qm.get_master_q(), + qm.get_minions_q()) diff --git a/interp/config.py b/interp/config.py new file mode 100644 index 0000000..85460c5 --- /dev/null +++ b/interp/config.py @@ -0,0 +1,19 @@ +import os + +import json + +default_config = { + 'filename': '/tmp/interp.log', + 'level': 'debug', + 'size': 102400, + 'logbackup': 10, + 'pypath': None, +} + +try: + with open(os.path.expanduser('~/.config/interp.json')) as config_file: + d = json.load(config_file) +except IOError as e: + d = {} + +config = dict(default_config.items() + d.items()) diff --git a/interp/grid/DD.py b/interp/grid/DD.py index 669155b..a107cc4 100644 --- a/interp/grid/DD.py +++ b/interp/grid/DD.py @@ -1,10 +1,9 @@ -from interp.grid.delaunay import dgrid as basegrid -from interp.tools import baker_exact_2D as exact_func - from itertools import product import numpy as np +from interp.grid.delaunay import dgrid as basegrid + class rect_grid(basegrid): def __init__(self, xres = 5, yres = 5): xmin = 0.0 diff --git a/interp/grid/DDD.py b/interp/grid/DDD.py index 4ec579f..394cd6c 100644 --- a/interp/grid/DDD.py +++ b/interp/grid/DDD.py @@ -1,10 +1,9 @@ -from interp.grid.delaunay import dgrid as basegrid -from interp.tools import baker_exact_3D, log - from itertools import product import numpy as np +from interp.grid.delaunay import dgrid as basegrid + class rect_grid(basegrid): def __init__(self, xres = 5, yres = 5, zres = 5): xmin = 0.0 @@ -22,7 +21,6 @@ class rect_grid(basegrid): zspan = zmaz - zmin zdel = zspan / float(zres - 1) - verts = [] q = np.zeros(xres * yres * zres) for x in xrange(xres): @@ -41,8 +39,6 @@ class random_grid(rect_grid): def __init__(self, num_verts = 100): verts = [] - r = np.random - appx_side_res = int(np.power(num_verts, 1/3.0)) delta = 1.0 / float(appx_side_res) diff --git a/interp/grid/__init__.py b/interp/grid/__init__.py index deec79e..3c58fc0 100644 --- a/interp/grid/__init__.py +++ b/interp/grid/__init__.py @@ -1,5 +1,4 @@ -import sys -from collections import defaultdict +from collections import defaultdict import pickle from xml.dom.minidom import Document @@ -9,256 +8,265 @@ from scipy.spatial import KDTree from interp.baker import run_baker from interp.baker import get_phis +import interp import logging log = logging.getLogger("interp") MAX_SEARCH_COUNT = 256 +TOL = 1e-8 + +__version__ = interp.__version__ + class grid(object): - def __init__(self, verts = None, q = None): - """ - verts = array of arrays (if passed in, will convert to numpy.array) - [ - [x0,y0 <, z0>], - [x1,y1 <, z1>], - ... - ] + def __init__(self, verts=None, q=None): + """ + verts = array of arrays (if passed in, will convert to numpy.array) + [ + [x0,y0 <, z0>], + [x1,y1 <, z1>], + ... + ] - q = array (1D) of physical values - """ + q = array (1D) of physical values + """ - if verts != None: - self.verts = np.array(verts) - self.tree = KDTree(self.verts) + if verts != None: + self.verts = np.array(verts) + self.tree = KDTree(self.verts) - if q != None: - self.q = np.array(q) + if q != None: + self.q = np.array(q) - self.cells = {} - self.cells_for_vert = defaultdict(list) + self.cells = {} + self.cells_for_vert = defaultdict(list) - def get_containing_simplex(self, X): - if not self.cells: - raise Exception("cell connectivity is not set up") + def get_containing_simplex(self, X): + if not self.cells: + raise Exception("cell connectivity is not set up") - # get closest point - (dist, indicies) = self.tree.query(X, 2) - closest_point = indicies[0] + # get closest point + (dist, indicies) = self.tree.query(X, 2) + closest_point = indicies[0] - log.debug('X: %s' % X) - log.debug('point index: %d' % closest_point) - log.debug('actual point %s' % self.verts[closest_point]) - log.debug('distance = %0.4f' % dist[0]) + log.debug('X: %s' % X) + log.debug('point index: %d' % closest_point) + log.debug('actual point %s' % self.verts[closest_point]) + log.debug('distance = %0.4f' % dist[0]) - simplex = None - checked_cells = [] - cells_to_check = list(self.cells_for_vert[closest_point]) + simplex = None + checked_cells = [] + cells_to_check = list(self.cells_for_vert[closest_point]) - attempts = 0 - while not simplex and cells_to_check: - attempts += 1 + attempts = 0 + while not simplex and cells_to_check: + attempts += 1 - if attempts > MAX_SEARCH_COUNT: - raise Exception("Is the search becoming exhaustive? (%d attempts)" % attempts) + if attempts > MAX_SEARCH_COUNT: + raise Exception("Is the search becoming exhaustive?'\ + '(%d attempts)" % attempts) - cur_cell = cells_to_check.pop(0) - checked_cells.append(cur_cell) + cur_cell = cells_to_check.pop(0) + checked_cells.append(cur_cell) - if cur_cell.contains(X, self): - simplex = cur_cell - continue + if cur_cell.contains(X, self): + simplex = cur_cell + continue - for neighbor in cur_cell.neighbors: - if (neighbor not in checked_cells) and (neighbor not in cells_to_check): - cells_to_check.append(neighbor) + for neighbor in cur_cell.neighbors: + if (neighbor not in checked_cells) \ + and (neighbor not in cells_to_check): + cells_to_check.append(neighbor) - if not simplex: - raise Exception('no containing simplex found') + if not simplex: + raise Exception('no containing simplex found') - log.debug("simplex vert indicies: %s" % simplex.verts) - R = self.create_mesh(simplex.verts) - log.debug("R:\n%s", R) + log.debug("simplex vert indicies: %s" % simplex.verts) + R = self.create_mesh(simplex.verts) + log.debug("R:\n%s", R) - log.debug('total attempts before finding simplex: %d' % attempts) - return R + log.debug('total attempts before finding simplex: %d' % attempts) + return R - def create_mesh(self, indicies): - """ - this function takes a list of indicies, and then creates and returns a - grid object (collection of verts and q). + def create_mesh(self, indicies): + """ + this function takes a list of indicies, and then creates and + returns a grid object (collection of verts and q). - note: the input is indicies, the grid contains verts - """ + note: the input is indicies, the grid contains verts + """ - return grid(self.verts[indicies], self.q[indicies]) + return grid(self.verts[indicies], self.q[indicies]) - def get_simplex_and_nearest_points(self, X, extra_points = 3): - """ - this returns two grid objects: R and S. + def get_simplex_and_nearest_points(self, X, extra_points=3): + """ + this returns two grid objects: R and S. - R is a grid object that is a containing simplex around point X + R is a grid object that is a containing simplex around point X - S : some verts from all points that are not the simplex - """ - simplex_size = self.dim + 1 - log.debug("extra verts: %d" % extra_points) - log.debug("simplex size: %d" % simplex_size) + S : some verts from all points that are not the simplex + """ + simplex_size = self.dim + 1 + log.debug("extra verts: %d" % extra_points) + log.debug("simplex size: %d" % simplex_size) - r_mesh = self.get_containing_simplex(X) + r_mesh = self.get_containing_simplex(X) - # and some UNIQUE extra verts - (dist, indicies) = self.tree.query(X, simplex_size + extra_points) - log.debug("extra indicies: %s" % indicies) + # and some UNIQUE extra verts + (dist, indicies) = self.tree.query(X, simplex_size + extra_points) + log.debug("extra indicies: %s" % indicies) - unique_indicies = [] - for index in indicies: - close_point_in_R = False - for rvert in r_mesh.verts: - if all(rvert == self.verts[index]): - close_point_in_R = True - break + unique_indicies = [] + for index in indicies: + close_point_in_R = False + for rvert in r_mesh.verts: + if all(rvert == self.verts[index]): + close_point_in_R = True + break - if not close_point_in_R: - unique_indicies.append(index) - else: - log.debug('throwing out %s: %s' % (index, self.verts[index])) + if not close_point_in_R: + unique_indicies.append(index) + else: + log.debug('throwing out %s: %s' % (index, self.verts[index])) - log.debug("indicies: %s" % indicies) - log.debug("unique indicies: %s" % unique_indicies) - s_mesh = self.create_mesh(unique_indicies) + log.debug("indicies: %s" % indicies) + log.debug("unique indicies: %s" % unique_indicies) + s_mesh = self.create_mesh(unique_indicies) - return (r_mesh, s_mesh) + return (r_mesh, s_mesh) - def run_baker(self, X, order = 2, extra_points = 3): - (R, S) = self.get_simplex_and_nearest_points(X, extra_points) - answer = run_baker(X, R, S, order) - return answer + def run_baker(self, X, order=2, extra_points=3): + (R, S) = self.get_simplex_and_nearest_points(X, extra_points) + answer = run_baker(X, R, S, order) + return answer - def for_qhull_generator(self): - """ - this returns a generator that should be fed into qdelaunay - """ + def for_qhull_generator(self): + """ + this returns a generator that should be fed into qdelaunay + """ - yield str(len(self.verts[0])); - yield '%d' % len(self.verts) + yield str(len(self.verts[0])) + yield '%d' % len(self.verts) - for p in self.verts: - yield "%f %f %f" % tuple(p) + for p in self.verts: + yield "%f %f %f" % tuple(p) - def for_qhull(self): - """ - this returns a single string that should be fed into qdelaunay - """ - r = '%d\n' % len(self.verts[0]) - r += '%d\n' % len(self.verts) - for p in self.verts: - # r += "%f %f %f\n" % tuple(p) - r += "%s\n" % " ".join("%f" % i for i in p) - return r + def for_qhull(self): + """ + this returns a single string that should be fed into qdelaunay + """ + r = '%d\n' % len(self.verts[0]) + r += '%d\n' % len(self.verts) + for p in self.verts: + # r += "%f %f %f\n" % tuple(p) + r += "%s\n" % " ".join("%f" % i for i in p) + return r - def __str__(self): - r = '' - assert( len(self.verts) == len(self.q) ) - for c, i in enumerate(zip(self.verts, self.q)): - r += "%d vert(%s): q(%0.4f)" % (c,i[0], i[1]) - cell_str = ", ".join([str(f.name) for f in self.cells_for_vert[c]]) - r += " cells: [%s]" % cell_str - r += "\n" - if self.cells: - for v in self.cells.itervalues(): - r += "%s\n" % v - return r + def __str__(self): + r = '' + assert(len(self.verts) == len(self.q)) + for c, i in enumerate(zip(self.verts, self.q)): + r += "%d vert(%s): q(%0.4f)" % (c, i[0], i[1]) + cell_str = ", ".join([str(f.name) for f in self.cells_for_vert[c]]) + r += " cells: [%s]" % cell_str + r += "\n" + if self.cells: + for v in self.cells.itervalues(): + r += "%s\n" % v + return r - def normalize_q(self, new_max = 0.1): - largest_number = np.max(np.abs(self.q)) - self.q *= new_max/largest_number + def normalize_q(self, new_max=0.1): + largest_number = np.max(np.abs(self.q)) + self.q *= new_max / largest_number + def dump_to_blender_files(self, + pfile='/tmp/points.p', cfile='/tmp/cells.p'): + if len(self.verts[0]) == 2: + pickle.dump([(p[0], p[1], 0.0) for p in self.verts], + open(pfile, 'w')) + else: + pickle.dump([(p[0], p[1], p[2]) for p in self.verts], + open(pfile, 'w')) - def dump_to_blender_files(self, pfile = '/tmp/points.p', cfile = '/tmp/cells.p'): - if len(self.verts[0]) == 2: - pickle.dump([(p[0], p[1], 0.0) for p in self.verts], open(pfile, 'w')) - else: - pickle.dump([(p[0], p[1], p[2]) for p in self.verts], open(pfile, 'w')) + pickle.dump([f.verts for f in self.cells.itervalues()], + open(cfile, 'w')) - pickle.dump([f.verts for f in self.cells.itervalues()], open(cfile, 'w')) + def get_xml(self): + doc = Document() + ps = doc.createElement("points") + doc.appendChild(ps) + for i in zip(self.verts, self.q): + p = doc.createElement("point") - def get_xml(self): - doc = Document() - ps = doc.createElement("points") - doc.appendChild(ps) - for i in zip(self.verts, self.q): - p = doc.createElement("point") + p.setAttribute("x", str(i[0][0])) + p.setAttribute('y', str(i[0][1])) + p.setAttribute('z', str(i[0][2])) + p.setAttribute('q', str(i[1])) + ps.appendChild(p) - p.setAttribute("x", str(i[0][0])) - p.setAttribute('y', str(i[0][1])) - p.setAttribute('z', str(i[0][2])) - p.setAttribute('q', str(i[1] )) - ps.appendChild(p) + return doc - return doc + def toxml(self): + return self.get_xml().toxml() - def toxml(self): - return self.get_xml().toxml() - def toprettyxml(self): - return self.get_xml().toprettyxml() + def toprettyxml(self): + return self.get_xml().toprettyxml() class cell(object): - def __init__(self, name): - self.name = name - self.verts = [] - self.neighbors = [] + def __init__(self, name): + self.name = name + self.verts = [] + self.neighbors = [] - def add_vert(self, v): - """ - v should be an index into grid.verts - """ - self.verts.append(v) + def add_vert(self, v): + """ + v should be an index into grid.verts + """ + self.verts.append(v) - def add_neighbor(self, n): - """ - reference to another cell object - """ - self.neighbors.append(n) + def add_neighbor(self, n): + """ + reference to another cell object + """ + self.neighbors.append(n) - def contains(self, X, G): - """ - X = point of interest - G = corrensponding grid object (G.verts) + def contains(self, X, G): + """ + X = point of interest + G = corrensponding grid object (G.verts) - because of the way i'm storing things, a cell simply stores indicies, - and so one must pass in a reference to the grid object containing real - verts. + because of the way i'm storing things, a cell simply stores + indicies, and so one must pass in a reference to the grid object + containing real verts. - this simply calls grid.simplex.contains - """ - return contains(X, [G.verts[i] for i in self.verts]) + this simply calls grid.simplex.contains + """ + return contains(X, [G.verts[i] for i in self.verts]) - def __str__(self): - # neighbors = [str(i.name) for i in self.neighbors] - return '' %\ - ( - self.name, - self.verts, - len(self.neighbors), - # ", ".join(neighbors) - ) + def __str__(self): + # neighbors = [str(i.name) for i in self.neighbors] + return '' %\ + ( + self.name, + self.verts, + len(self.neighbors), + # ", ".join(neighbors) + ) - __repr__ = __str__ + __repr__ = __str__ -TOL = 1e-8 - def contains(X, R): - """ - tests if X (point) is in R + """ + tests if X (point) is in R - R is a simplex, represented by a list of n-degree coordinates - """ - phis = get_phis(X, R) + R is a simplex, represented by a list of n-degree coordinates + """ + phis = get_phis(X, R) - r = True - if [i for i in phis if i < 0.0 - TOL]: - r = False - return r + r = True + if [i for i in phis if i < 0.0 - TOL]: + r = False + return r diff --git a/interp/grid/gmsh.py b/interp/grid/gmsh.py index 58c8f78..6b8e818 100644 --- a/interp/grid/gmsh.py +++ b/interp/grid/gmsh.py @@ -1,7 +1,4 @@ -import pickle - from itertools import combinations -from collections import defaultdict import numpy as np from scipy.spatial import KDTree @@ -36,7 +33,7 @@ class ggrid(grid): gmsh_file.readline() # $MeshFormat - fmat = gmsh_file.readline() + gmsh_file.readline() gmsh_file.readline() # $EndMeshFormat gmsh_file.readline() # $Nodes diff --git a/interp/tools.py b/interp/tools.py index 523687d..b7fdd19 100644 --- a/interp/tools.py +++ b/interp/tools.py @@ -1,82 +1,75 @@ -import os - import numpy as np -import logging -log = logging.getLogger("interp") def rms(errors): - """ - root mean square calculation - """ + """ + root mean square calculation + """ - # slow pure python way for reference: - # r = 0.0 - # for i in errors: - # r += np.power(i, 2) - # r = np.sqrt(r / len(errors)) - # return r + # slow pure python way for reference: + # r = 0.0 + # for i in errors: + # r += np.power(i, 2) + # r = np.sqrt(r / len(errors)) + # return r + + return np.sqrt((errors ** 2).mean()) - return np.sqrt((errors**2).mean()) def baker_exact_2D(X): - """ - the exact function (2D) used from baker's article (for testing, slightly - modified) - """ - x ,y = X + """ + the exact function (2D) used from baker's article (for testing, + slightly modified) + """ + x, y = X + + answer = np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) + return answer - answer = np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) - log.debug(answer) - return answer def friendly_exact_2D(X): - """ - A friendlier 2D func - """ - x ,y = X - answer = 1.0 + x*x + y*y - log.debug(answer) - return answer + """ + A friendlier 2D func + """ + x, y = X + answer = 1.0 + x * x + y * y + return answer + def baker_exact_3D(X): - """ - the exact function (3D) used from baker's article (for testing) - """ - x = X[0] - y = X[1] - z = X[2] - answer = np.power((np.sin(x * np.pi / 2.0) * np.sin(y * np.pi / 2.0) * np.sin(z * np.pi / 2.0)), 2) - log.debug(answer) - return answer + """ + the exact function (3D) used from baker's article (for testing) + """ + x, y, z = X + answer = np.power((np.sin(x * np.pi / 2.0) * np.sin(y * np.pi / 2.0) * + np.sin(z * np.pi / 2.0)), 2) + return answer + def friendly_exact_3D(X): - x,y,z = X - return 1 + x*x + y*y + z*z + x, y, z = X + return 1 + x * x + y * y + z * z + def scipy_exact_2D(X): - x,y = X - return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2 + x, y = X + return x * (1 - x) * np.cos(4 * np.pi * x) *\ + np.sin(4 * np.pi * y ** 2) ** 2 + def improved_answer(answer, exact): - if not answer['error']: - # was probably just a linear interpolation - return False + if not answer['error']: + # was probably just a linear interpolation + return False - log.debug('qlin: %s' % answer['qlin']) - log.debug('error: %s' % answer['error']) - log.debug('final: %s' % answer['final']) - log.debug('exact: %s' % exact) + if np.abs(answer['final'] - exact) <= np.abs(answer['qlin'] - exact): + return True + else: + return False - if np.abs(answer['final'] - exact) <= np.abs(answer['qlin'] - exact): - log.debug(":) improved result") - return True - else: - log.debug(":( damaged result") - return False def improved(qlin, err, final, exact): - if np.abs(final - exact) <= np.abs(qlin - exact): - return True - else: - return False + if np.abs(final - exact) <= np.abs(qlin - exact): + return True + else: + return False diff --git a/test/all.py b/test/all.py index 88856e8..1888a26 100644 --- a/test/all.py +++ b/test/all.py @@ -11,14 +11,14 @@ import quadratic2d if __name__ == '__main__': - tests = [ - unittest.TestLoader().loadTestsFromTestCase(baker2dorder.Test), - unittest.TestLoader().loadTestsFromTestCase(baker2d.Test), - unittest.TestLoader().loadTestsFromTestCase(baker3d.Test), - unittest.TestLoader().loadTestsFromTestCase(cubic2d.Test), - unittest.TestLoader().loadTestsFromTestCase(pattern.Test), - unittest.TestLoader().loadTestsFromTestCase(quadratic2d.Test), - ] + tests = [ + unittest.TestLoader().loadTestsFromTestCase(baker2dorder.Test), + unittest.TestLoader().loadTestsFromTestCase(baker2d.Test), + unittest.TestLoader().loadTestsFromTestCase(baker3d.Test), + unittest.TestLoader().loadTestsFromTestCase(cubic2d.Test), + unittest.TestLoader().loadTestsFromTestCase(pattern.Test), + unittest.TestLoader().loadTestsFromTestCase(quadratic2d.Test), + ] - for test in tests: - unittest.TextTestRunner(verbosity=3).run(test) + for test in tests: + unittest.TextTestRunner(verbosity=3).run(test) diff --git a/test/baker2d.py b/test/baker2d.py index 5123262..599b6e2 100644 --- a/test/baker2d.py +++ b/test/baker2d.py @@ -2,154 +2,154 @@ import unittest -from interp import baker -from interp import grid +from interp import baker -import numpy as np -import scipy.spatial class Test(unittest.TestCase): - def setUp(self): - self.l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]] - self.all_points = [ - [ 0, 0], # 0 - [ 1, 0], # 1 - [ 1, 1], # 2 - [ 0, 1], # 3 - [ 1,-1], # 4 - [ 0,-1], # 5 - [-1, 1], # 6 - [-1, 0], # 7 - [-1,-1], # 8 - ] - self.q = [1, 0, 0, 0, 0, 0, 0, 0, 0] - self.X = [0.5, 0.25] - self.accuracy = 8 + def setUp(self): + self.l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], + [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]] + self.all_points = [ + [0, 0], # 0 + [1, 0], # 1 + [1, 1], # 2 + [0, 1], # 3 + [1, -1], # 4 + [0, -1], # 5 + [-1, 1], # 6 + [-1, 0], # 7 + [-1, -1], # 8 + ] + self.q = [1, 0, 0, 0, 0, 0, 0, 0, 0] + self.X = [0.5, 0.25] + self.accuracy = 8 - def testImports(self): - import numpy - import scipy - import interp.grid - import interp.baker + def testImports(self): + import numpy + import scipy + import interp.grid as gv + import interp.baker as bv - def testGetPhis(self): + numpy.__version__ + scipy.__version__ - X = [0,0] - r = [[-1, -1], [0, 2], [1, -1]] + gv, bv - result = baker.get_phis(X, r) + def testGetPhis(self): + X = [0, 0] + r = [[-1, -1], [0, 2], [1, -1]] - right_answer = [1/3.0, 1/3.0, 1/3.0] + result = baker.get_phis(X, r) - for a,b in zip(result, right_answer): - self.assertAlmostEqual(a,b) + right_answer = [1 / 3.0, 1 / 3.0, 1 / 3.0] - def testGetPhis2(self): + for a, b in zip(result, right_answer): + self.assertAlmostEqual(a, b) - X = [0.5,0.25] - r = [[0, 0], [1, 0], [1, 1]] + def testGetPhis2(self): + X = [0.5, 0.25] + r = [[0, 0], [1, 0], [1, 1]] - result = baker.get_phis(X, r) + result = baker.get_phis(X, r) - right_answer = [0.5, 0.25, 0.25] + right_answer = [0.5, 0.25, 0.25] - for a,b in zip(result, right_answer): - self.assertEqual(a,b) + for a, b in zip(result, right_answer): + self.assertEqual(a, b) - def testQlinear(self): - X = [0.5, 0.25] - r = [[0, 0], [1, 0], [1, 1]] - q = [1, 0, 0] + def testQlinear(self): + X = [0.5, 0.25] + r = [[0, 0], [1, 0], [1, 1]] + q = [1, 0, 0] - phi, result = baker.qlinear(X, grid.grid(r,q)) + phi, result = baker.qlinear(X, r, q) - right_answer = 0.5 + right_answer = 0.5 - self.assertAlmostEqual(result, right_answer) + self.assertAlmostEqual(result, right_answer) - def testRunBaker_1(self): - size_of_simplex = 3 - extra_points = 3 + def testRunBaker_1(self): + size_of_simplex = 3 + extra_points = 3 - R = grid.grid(self.all_points[:size_of_simplex], - self.q[:size_of_simplex]) + R, R_q = (self.all_points[:size_of_simplex], + self.q[:size_of_simplex]) - S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], - self.q[size_of_simplex:size_of_simplex + extra_points]) + S, S_q = (self.all_points[size_of_simplex:size_of_simplex \ + + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + answer = baker.run_baker(self.X, R, R_q, S, S_q) - answer = baker.run_baker(self.X, R, S) + a = answer['abc'][0] + b = answer['abc'][1] + c = answer['abc'][2] - a = answer['abc'][0] - b = answer['abc'][1] - c = answer['abc'][2] + self.assertEqual(sorted((a, b, c)), sorted((0, 0.0, 1 / 3.))) - self.assertEqual(sorted((a,b,c)), sorted((0,0.0,1/3.))) + def testRunBaker_2(self): + size_of_simplex = 3 + extra_points = 4 - def testRunBaker_2(self): - size_of_simplex = 3 - extra_points = 4 + R, R_q = (self.all_points[:size_of_simplex], self.q[:size_of_simplex]) - R = grid.grid(self.all_points[:size_of_simplex], - self.q[:size_of_simplex]) + S, S_q = (self.all_points[size_of_simplex:size_of_simplex \ + + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) - S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], - self.q[size_of_simplex:size_of_simplex + extra_points]) + answer = baker.run_baker(self.X, R, R_q, S, S_q) - answer = baker.run_baker(self.X, R, S) + a, b, c = sorted(answer['abc']) + aa, bb, cc = sorted((2 / 3.0, 2 / 3.0, 1 / 3.0)) - a, b, c = sorted(answer['abc']) - aa,bb,cc = sorted((2/3.0, 2/3.0, 1/3.0)) + self.assertAlmostEqual(a, aa) + self.assertAlmostEqual(b, bb) + self.assertAlmostEqual(c, cc) - self.assertAlmostEqual(a,aa) - self.assertAlmostEqual(b,bb) - self.assertAlmostEqual(c,cc) + def testRunBaker_3(self): + size_of_simplex = 3 + extra_points = 5 - def testRunBaker_3(self): - size_of_simplex = 3 - extra_points = 5 + R, R_q = (self.all_points[:size_of_simplex], self.q[:size_of_simplex]) - R = grid.grid(self.all_points[:size_of_simplex], - self.q[:size_of_simplex]) + S, S_q = (self.all_points[size_of_simplex:size_of_simplex \ + + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + answer = baker.run_baker(self.X, R, R_q, S, S_q) - S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], - self.q[size_of_simplex:size_of_simplex + extra_points]) + a = answer['abc'][0] + b = answer['abc'][1] + c = answer['abc'][2] - answer = baker.run_baker(self.X, R, S) + a, b, c = sorted((a, b, c)) + aa, bb, cc = sorted((13 / 14., 2 / 7., 15 / 14.)) - a = answer['abc'][0] - b = answer['abc'][1] - c = answer['abc'][2] + self.assertAlmostEqual(a, aa) + self.assertAlmostEqual(b, bb) + self.assertAlmostEqual(c, cc) - a,b,c = sorted((a,b,c)) - aa, bb, cc = sorted((13/14., 2/7., 15/14.)) + def testRunBaker_4(self): + size_of_simplex = 3 + extra_points = 6 - self.assertAlmostEqual(a,aa) - self.assertAlmostEqual(b,bb) - self.assertAlmostEqual(c,cc) + R, R_q = (self.all_points[:size_of_simplex], + self.q[:size_of_simplex]) + S, S_q = (self.all_points[size_of_simplex:size_of_simplex \ + + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + answer = baker.run_baker(self.X, R, R_q, S, S_q) - def testRunBaker_4(self): - size_of_simplex = 3 - extra_points = 6 + a = answer['abc'][0] + b = answer['abc'][1] + c = answer['abc'][2] - R = grid.grid(self.all_points[:size_of_simplex], - self.q[:size_of_simplex]) + a, b, c = sorted((a, b, c)) + aa, bb, cc = sorted((48 / 53.0, 15 / 53.0, 54 / 53.0)) - S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], - self.q[size_of_simplex:size_of_simplex + extra_points]) - - answer = baker.run_baker(self.X, R, S) - a = answer['abc'][0] - b = answer['abc'][1] - c = answer['abc'][2] - - a,b,c = sorted((a,b,c)) - aa,bb,cc = sorted((48/53.0, 15/53.0, 54/53.0)) - - self.assertAlmostEqual(a, aa) - self.assertAlmostEqual(b, bb) - self.assertAlmostEqual(c, cc) + self.assertAlmostEqual(a, aa) + self.assertAlmostEqual(b, bb) + self.assertAlmostEqual(c, cc) if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(Test) - unittest.TextTestRunner(verbosity=3).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase(Test) + unittest.TextTestRunner(verbosity=3).run(suite) diff --git a/test/baker2dorder.py b/test/baker2dorder.py index 974f924..bbc7315 100644 --- a/test/baker2dorder.py +++ b/test/baker2dorder.py @@ -8,97 +8,102 @@ import numpy as np from interp.grid import contains -def exact_func(point): - x = point[0] - y = point[1] - return 0.5 + x*x + y -def calculate_error_term(self, a,b,c,d,e,f): +def exact_func(point): + x = point[0] + y = point[1] + return 0.5 + x * x + y + + +def calculate_error_term(self, a, b, c, d, e, f): B = np.array([ - self.p1[a] * self.p1[b], self.p1[c] * self.p1[d], self.p1[e] * self.p1[f], - self.p2[a] * self.p2[b], self.p2[c] * self.p2[d], self.p2[e] * self.p2[f], - self.p3[a] * self.p3[b], self.p3[c] * self.p3[d], self.p3[e] * self.p3[f], - self.p4[a] * self.p4[b], self.p4[c] * self.p4[d], self.p4[e] * self.p4[f], - ]) - B.shape = (4,3) + self.p1[a] * self.p1[b], self.p1[c] * self.p1[d], self.p1[e] * self.p1[f], + self.p2[a] * self.p2[b], self.p2[c] * self.p2[d], self.p2[e] * self.p2[f], + self.p3[a] * self.p3[b], self.p3[c] * self.p3[d], self.p3[e] * self.p3[f], + self.p4[a] * self.p4[b], self.p4[c] * self.p4[d], self.p4[e] * self.p4[f], + ]) + + B.shape = (4, 3) A = np.dot(B.T, B) rhs = np.dot(B.T, self.w) - abc = np.linalg.solve(A,rhs) + abc = np.linalg.solve(A, rhs) err = \ - abc[0] * self.phis[a] * self.phis[b] + \ - abc[1] * self.phis[c] * self.phis[d] + \ - abc[2] * self.phis[e] * self.phis[f] + abc[0] * self.phis[a] * self.phis[b] + \ + abc[1] * self.phis[c] * self.phis[d] + \ + abc[2] * self.phis[e] * self.phis[f] return err + class Test(unittest.TestCase): - def setUp(self): - self.verts = [ - [ 2, 3], # 0 - [ 7, 4], # 1 - [ 4, 8], # 2 - [ 0, 7], # 3, 1 - [ 5, 0], # 4, 2 - [10, 5], # 5, 3 - [ 8, 9], # 6, 4 - ] + def setUp(self): + self.verts = [ + [2, 3], # 0 + [7, 4], # 1 + [4, 8], # 2 + [0, 7], # 3, 1 + [5, 0], # 4, 2 + [0, 5], # 5, 3 + [8, 9], # 6, 4 + ] + self.q = [exact_func(v) for v in self.verts] - self.q = [exact_func(v) for v in self.verts] + self.g = grid(self.verts, self.q) + self.R, self.R_q = (self.verts[:3], self.q[:3]) + self.S, self.S_q = (self.verts[3:], self.q[3:]) - self.g = grid(self.verts, self.q) - self.R = grid(self.verts[:3], self.q[:3]) - self.S = grid(self.verts[3:], self.q[3:]) + self.p1, self.ql1 = baker.qlinear(self.verts[3], self.R, self.q) + self.p2, self.ql2 = baker.qlinear(self.verts[4], self.R, self.q) + self.p3, self.ql3 = baker.qlinear(self.verts[5], self.R, self.q) + self.p4, self.ql4 = baker.qlinear(self.verts[6], self.R, self.q) - self.p1, self.ql1 = baker.qlinear(self.verts[3], self.R) - self.p2, self.ql2 = baker.qlinear(self.verts[4], self.R) - self.p3, self.ql3 = baker.qlinear(self.verts[5], self.R) - self.p4, self.ql4 = baker.qlinear(self.verts[6], self.R) + self.q1 = exact_func(self.verts[3]) + self.q2 = exact_func(self.verts[4]) + self.q3 = exact_func(self.verts[5]) + self.q4 = exact_func(self.verts[6]) - self.q1 = exact_func(self.verts[3]) - self.q2 = exact_func(self.verts[4]) - self.q3 = exact_func(self.verts[5]) - self.q4 = exact_func(self.verts[6]) + self.w = np.array([ + self.q1 - self.ql1, + self.q2 - self.ql2, + self.q3 - self.ql3, + self.q4 - self.ql4, + ]) + self.X = [4, 5] - self.w = np.array([ - self.q1 - self.ql1, - self.q2 - self.ql2, - self.q3 - self.ql3, - self.q4 - self.ql4, - ]) + self.g = grid(self.verts, self.q) - self.X = [4,5] + self.phis, self.qlin = baker.qlinear(self.X, self.R, self.q) + self.exact = exact_func(self.X) + self.answer = baker.run_baker(self.X, self.R, + self.R_q, self.S, self.S_q) - self.g = grid(self.verts, self.q) + def test_R_contains_X(self): + self.assertTrue(contains(self.X, self.R)) - self.phis, self.qlin = baker.qlinear(self.X, self.R) - self.exact = exact_func(self.X) - self.answer = baker.run_baker(self.X,self.R,self.S) + def test_1(self): + a, b, c, d, e, f = (0, 1, 1, 2, 2, 0) + err = calculate_error_term(self, a, b, c, d, e, f) + self.assertAlmostEqual(err, self.answer['error']) + def test_swap_first_elements(self): + a, b, c, d, e, f = (1, 0, 1, 2, 2, 0) + err = calculate_error_term(self, a, b, c, d, e, f) + self.assertAlmostEqual(err, self.answer['error']) - def test_R_contains_X(self): - self.assertTrue(contains(self.X, self.R.verts)) + def test_swap_two_pairs(self): + a, b, c, d, e, f = (1, 2, 0, 1, 2, 0) + err = calculate_error_term(self, a, b, c, d, e, f) + self.assertAlmostEqual(err, self.answer['error']) - def test_1(self): - a,b,c,d,e,f = (0,1, 1,2, 2,0) - err = calculate_error_term(self, a,b,c,d,e,f) - self.assertAlmostEqual(err, self.answer['error']) - def test_swap_first_elements(self): - a,b,c,d,e,f = (1,0, 1,2, 2,0) - err = calculate_error_term(self, a,b,c,d,e,f) - self.assertAlmostEqual(err, self.answer['error']) - def test_swap_two_pairs(self): - a,b,c,d,e,f = (1,2, 0,1, 2,0) - err = calculate_error_term(self, a,b,c,d,e,f) - self.assertAlmostEqual(err, self.answer['error']) - def test_swap_all_pairs(self): - a,b,c,d,e,f = (0,2, 0,1, 2,1) - err = calculate_error_term(self, a,b,c,d,e,f) - self.assertAlmostEqual(err, self.answer['error']) + def test_swap_all_pairs(self): + a, b, c, d, e, f = (0, 2, 0, 1, 2, 1) + err = calculate_error_term(self, a, b, c, d, e, f) + self.assertAlmostEqual(err, self.answer['error']) if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(Test) - unittest.TextTestRunner(verbosity=3).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase(Test) + unittest.TextTestRunner(verbosity=3).run(suite) diff --git a/test/baker3d.py b/test/baker3d.py index 87fefa7..0e607ef 100644 --- a/test/baker3d.py +++ b/test/baker3d.py @@ -2,38 +2,33 @@ import unittest from interp.baker import get_phis, qlinear -from interp.grid import grid -import numpy as np -import scipy.spatial class Test(unittest.TestCase): - def setUp(self): - self.X = [0.0, 0.0, 0.0] - self.r = [ - [0.0, 0.0, 1.0], - [0.94280904333606508, 0.0, -0.3333333283722672], - [-0.47140452166803232, 0.81649658244673617, -0.3333333283722672], - [-0.47140452166803298, -0.81649658244673584, -0.3333333283722672], - ] - self.q = [0.0, 0.0, 0.0, 4] + def setUp(self): + self.X = [0.0, 0.0, 0.0] + self.r = [ + [0.0, 0.0, 1.0], + [0.94280904333606508, 0.0, -0.3333333283722672], + [-0.47140452166803232, 0.81649658244673617, -0.3333333283722672], + [-0.47140452166803298, -0.81649658244673584, -0.3333333283722672], + ] + self.q = [0.0, 0.0, 0.0, 4] + def testGetPhis(self): + result = get_phis(self.X, self.r) + right_answer = [0.25, 0.25, 0.25, 0.25] - def testGetPhis(self): - result = get_phis(self.X, self.r) - right_answer = [0.25, 0.25, 0.25, 0.25] + for a, b in zip(result, right_answer): + self.assertAlmostEqual(a, b) - for a,b in zip(result, right_answer): - self.assertAlmostEqual(a,b) - - - def testQlinear(self): - phi, result = qlinear(self.X, grid(self.r, self.q)) - result = result - right_answer = 1.0 - self.assertAlmostEqual(result, right_answer) + def testQlinear(self): + phi, result = qlinear(self.X, self.r, self.q) + result = result + right_answer = 1.0 + self.assertAlmostEqual(result, right_answer) if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(Test) - unittest.TextTestRunner(verbosity=2).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase(Test) + unittest.TextTestRunner(verbosity=2).run(suite) diff --git a/test/cubic2d.py b/test/cubic2d.py index 27ffad0..53d64ed 100644 --- a/test/cubic2d.py +++ b/test/cubic2d.py @@ -4,71 +4,76 @@ import unittest from interp.baker import run_baker -from interp.grid import grid -from interp.grid import contains +from interp.grid import contains + def exact_func(X): - x = X[0] - y = X[0] - return 1 + x + y + x = X[0] + y = X[0] + return 1 + x + y + class Test(unittest.TestCase): - def setUp(self): - self.verts = [ - [ 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 - ] - self.q = [exact_func(p) for p in self.verts] + def setUp(self): + self.g = [[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 + ] + self.q = [exact_func(p) for p in self.g] + self.X = [0.55, 0.45] + self.R = self.g[0:3] + self.R_q = self.q[0:3] + self.exact = exact_func(self.X) - self.X = [0.55, 0.45] + def test_R_contains_X(self): + self.assertTrue(contains(self.X, self.R)) - self.g = grid(self.verts, self.q) - # self.g.construct_connectivity() - self.R = self.g.create_mesh(range(3)) + def test_RunBaker_1_extra_point(self, extra=1): + S = self.g[3:3 + extra] + S_q = self.q[3:3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q, order=3) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + # expected failure ... + self.assertTrue(lin_err >= final_err) - self.exact = exact_func(self.X) + def test_RunBaker_2_extra_point(self, extra=2): + S = self.g[3: 3 + extra] + S_q = self.q[3:3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q, order=3) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) + def test_RunBaker_3_extra_point(self, extra=3): + S = self.g[3: 3 + extra] + S_q = self.q[3:3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q, order=3) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) - def test_R_contains_X(self): - self.assertTrue(contains(self.X, self.R.verts)) + def test_RunBaker_4_extra_point(self, extra=4): + S = self.g[3: 3 + extra] + S_q = self.q[3:3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q, order=3) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) - def test_RunBaker_1_extra_point(self, extra=1): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S, order=3) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_2_extra_point(self, extra=2): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S, order=3) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_3_extra_point(self, extra=3): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S, order=3) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_4_extra_point(self, extra=4): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S, order=3) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_5_extra_point(self, extra=5): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S, order=3) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) + def test_RunBaker_5_extra_point(self, extra=5): + S = self.g[3: 3 + extra] + S_q = self.q[3:3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q, order=3) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(Test) - unittest.TextTestRunner(verbosity=3).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase(Test) + unittest.TextTestRunner(verbosity=3).run(suite) diff --git a/test/pattern.py b/test/pattern.py index f55b067..b951833 100644 --- a/test/pattern.py +++ b/test/pattern.py @@ -4,49 +4,47 @@ import unittest from interp.baker import pattern - class Test(unittest.TestCase): - def setUp(self): - pass + def setUp(self): + pass - def testImports(self): - from interp.baker import pattern + def testImports(self): + from interp.baker import pattern as ppp + ppp - def test_baker_eq_8(self): - b = sorted([tuple(sorted(i)) for i in ((0,1),(1,2),(2,0))]) - p = sorted(pattern(3,2)) - self.assertEqual(b,p) + def test_baker_eq_8(self): + b = sorted([tuple(sorted(i)) for i in ((0, 1), (1, 2), (2, 0))]) + p = sorted(pattern(3, 2)) + self.assertEqual(b, p) - def test_baker_eq_17(self): - b = sorted([tuple(sorted(i)) for i in ((0,1,1), (0,2,2), (1,0,0), (1,2,2), (2,0,0), (2,1,1), (0,1,2))]) - p = sorted(pattern(3,3)) - self.assertEqual(b,p) + def test_baker_eq_17(self): + b = sorted([tuple(sorted(i)) for i in ((0, 1, 1), (0, 2, 2), (1, 0, 0), + (1, 2, 2), (2, 0, 0), (2, 1, 1), (0, 1, 2))]) + p = sorted(pattern(3, 3)) + self.assertEqual(b, p) - def test_baker_eq_15(self): - b = sorted([tuple(sorted(i)) for i in ( - (0,1), (0,2), (0,3), - (1,2), (1,3), (2,3))]) + def test_baker_eq_15(self): + b = sorted([tuple(sorted(i)) for i in ( + (0, 1), (0, 2), (0, 3), + (1, 2), (1, 3), (2, 3))]) - p = sorted(pattern(4,2)) - - self.assertEqual(b,p) - - def test_smcquay_(self): - b = sorted([tuple(sorted(i)) for i in ( - (0,1,2), (1,2,3), (0,1,3), (0,2,3), - (0,0,1), (0,1,1), - (1,2,2), (1,1,2), - (0,2,2), (0,0,2), - (1,3,3), (1,1,3), - (2,2,3), (2,3,3), - (0,3,3), (0,0,3))]) - - p = sorted(pattern(4,3)) - self.assertEqual(b,p) + p = sorted(pattern(4, 2)) + self.assertEqual(b, p) + def test_smcquay_(self): + b = sorted([tuple(sorted(i)) for i in ( + (0, 1, 2), (1, 2, 3), (0, 1, 3), (0, 2, 3), + (0, 0, 1), (0, 1, 1), + (1, 2, 2), (1, 1, 2), + (0, 2, 2), (0, 0, 2), + (1, 3, 3), (1, 1, 3), + (2, 2, 3), (2, 3, 3), + (0, 3, 3), (0, 0, 3))]) + p = sorted(pattern(4, 3)) + self.assertEqual(b, p) if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(Test) - unittest.TextTestRunner(verbosity=3).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase(Test) + unittest.TextTestRunner(verbosity=3).run(suite) diff --git a/test/quadratic2d.py b/test/quadratic2d.py index fe45cc7..52fbb5e 100644 --- a/test/quadratic2d.py +++ b/test/quadratic2d.py @@ -4,76 +4,82 @@ import unittest from interp.baker import run_baker -from interp.grid import grid -from interp.grid import contains +from interp.grid import grid +from interp.grid import contains + def exact_func(X): - x = X[0] - y = X[0] - return 1 - x*x + y*y + x = X[0] + y = X[0] + return 1 - x * x + y * y + class Test(unittest.TestCase): - def setUp(self): - self.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 + def setUp(self): + self.g = [ + [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 ] - self.q = [exact_func(p) for p in self.points] + self.q = [exact_func(p) for p in self.g] - self.X = [0.25, 0.4001] - self.X = [0.55, 0.45] + self.X = [0.25, 0.4001] + self.X = [0.55, 0.45] - self.g = grid(self.points, self.q) - self.R = self.g.create_mesh(range(3)) + self.R = self.g[0:3] + self.R_q = self.q[0:3] + self.exact = exact_func(self.X) - self.exact = exact_func(self.X) + def test_R_contains_X(self): + self.assertTrue(contains(self.X, self.R)) + def test_RunBaker_1_extra_point(self, extra=1): + S = self.g[3: 3 + extra] + S_q = self.q[3: 3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) - self.accuracy = 8 + #XXX: not sure about this one: + self.assertEqual(lin_err, final_err) - def test_R_contains_X(self): - self.assertTrue(contains(self.X, self.R.verts)) + def test_RunBaker_2_extra_point(self, extra=2): + S = self.g[3: 3 + extra] + S_q = self.q[3: 3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) - def test_RunBaker_1_extra_point(self, extra=1): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) + def test_RunBaker_3_extra_point(self, extra=3): + S = self.g[3: 3 + extra] + S_q = self.q[3: 3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) - # I expect this one to be bad: - # self.assertTrue(lin_err >= final_err) + def test_RunBaker_4_extra_point(self, extra=4): + S = self.g[3: 3 + extra] + S_q = self.q[3: 3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) - def test_RunBaker_2_extra_point(self, extra=2): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_3_extra_point(self, extra=3): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_4_extra_point(self, extra=4): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) - def test_RunBaker_5_extra_point(self, extra=5): - S = self.g.create_mesh(range(3, 3 + extra)) - answer = run_baker(self.X, self.R, S) - lin_err = abs(self.exact - answer['qlin']) - final_err = abs(self.exact - answer['final']) - self.assertTrue(lin_err >= final_err) + def test_RunBaker_5_extra_point(self, extra=5): + S = self.g[3: 3 + extra] + S_q = self.q[3: 3 + extra] + answer = run_baker(self.X, self.R, self.R_q, S, S_q) + lin_err = abs(self.exact - answer['qlin']) + final_err = abs(self.exact - answer['final']) + self.assertTrue(lin_err >= final_err) if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(Test) - unittest.TextTestRunner(verbosity=3).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase(Test) + unittest.TextTestRunner(verbosity=3).run(suite)