From 961b3b7b2624a5464da1bde25e1f6349f5c62070 Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Sun, 27 Dec 2009 10:48:27 -0700 Subject: [PATCH] init --- .hgignore | 1 + .setup | 2 + bin/driver-random.py | 43 ++++++++++++++++ bin/driver-structured.py | 45 ++++++++++++++++ bin/generate_qhull.pl | 16 ++++++ data/mesh.shelve | Bin 0 -> 12288 bytes data/qhull_input.txt | 11 ++++ lib/baker.py | 106 ++++++++++++++++++++++++++++++++++++++ lib/grid.py | 77 +++++++++++++++++++++++++++ lib/smcqhull.py | 4 ++ lib/tools.py | 12 +++++ test/baker.test.py | 75 +++++++++++++++++++++++++++ test/qhull.test.py | 33 ++++++++++++ test/utest.py | 29 +++++++++++ tools/mkpoints.py | 26 ++++++++++ tools/multi/forktest.py | 22 ++++++++ tools/multi/pool.py | 12 +++++ tools/multi/smp.py | 24 +++++++++ tools/multi/threadtest.py | 22 ++++++++ 19 files changed, 560 insertions(+) create mode 100644 .hgignore create mode 100644 .setup create mode 100755 bin/driver-random.py create mode 100755 bin/driver-structured.py create mode 100755 bin/generate_qhull.pl create mode 100644 data/mesh.shelve create mode 100644 data/qhull_input.txt create mode 100644 lib/baker.py create mode 100755 lib/grid.py create mode 100644 lib/smcqhull.py create mode 100644 lib/tools.py create mode 100755 test/baker.test.py create mode 100755 test/qhull.test.py create mode 100755 test/utest.py create mode 100755 tools/mkpoints.py create mode 100755 tools/multi/forktest.py create mode 100755 tools/multi/pool.py create mode 100755 tools/multi/smp.py create mode 100755 tools/multi/threadtest.py diff --git a/.hgignore b/.hgignore new file mode 100644 index 0000000..493485d --- /dev/null +++ b/.hgignore @@ -0,0 +1 @@ +\.pyc$ diff --git a/.setup b/.setup new file mode 100644 index 0000000..b25a3fb --- /dev/null +++ b/.setup @@ -0,0 +1,2 @@ +export PATH=~/src/baker/bin:~/src/baker/test:$PATH +export PYTHONPATH=~/src/baker/lib diff --git a/bin/driver-random.py b/bin/driver-random.py new file mode 100755 index 0000000..2658efe --- /dev/null +++ b/bin/driver-random.py @@ -0,0 +1,43 @@ +#!/usr/bin/python + +import sys +from optparse import OptionParser + +import numpy as np +import scipy.spatial + +from grid import simple_random_grid, exact_func +import baker + +def rms(errors): + r = 0.0 + for i in errors: + r += np.power(i, 2) + r = np.sqrt(r / len(errors)) + return r + +if __name__ == '__main__': + parser = OptionParser() + parser.add_option("-o", "--output-file", dest="output", type='str', default = '/tmp/for_qhull.txt', help = "qhull output file") + parser.add_option("-e", "--extra-points", dest="extra", type='int', default = 3, help = "how many extra points") + parser.add_option("-s", "--source-total", dest="source_total", type='int', default = 100, help = "total number of source points") + parser.add_option("-d", "--destination-total",dest="destination_total",type='int', default = 100, help = "total number of destination points") + parser.add_option("-v", "--verbose", action = 'store_true', default = False, help = "verbosity") + + (options, args) = parser.parse_args() + print >> sys.stderr, options + + mesh_source = simple_random_grid(options.source_total) + tree = scipy.spatial.KDTree(mesh_source.points) + mesh_dest = simple_random_grid(options.destination_total) + + open(options.output, 'w').write(mesh_source.for_qhull()) + print >> sys.stderr, "wrote output to %s" % options.output + + errors = [] + for x in mesh_dest.points: + (final, exact) = baker.run_baker(x, mesh_source, tree, options.extra, options.verbose) + cur_error = np.abs(final - exact) + errors.append(cur_error) + + print rms(errors) diff --git a/bin/driver-structured.py b/bin/driver-structured.py new file mode 100755 index 0000000..b547b05 --- /dev/null +++ b/bin/driver-structured.py @@ -0,0 +1,45 @@ +#!/usr/bin/python + +import sys +from optparse import OptionParser + +import numpy as np +import scipy.spatial + +from grid import simple_rect_grid, exact_func +import baker +from tools import rms + + +if __name__ == '__main__': + parser = OptionParser() + parser.add_option("-e", "--extra-points", dest="extra", type='int', default = 3, help = "how many extra points") + parser.add_option("-s", "--source-density", dest="source", type='int', default = 11, help = "resolution of source mesh") + parser.add_option("-d", "--dest-density", dest="dest", type='int', default = 11, help = "resolution of dest mesh") + parser.add_option("-t", "--random-total", dest="random_total", type='int', default = 100, help = "total number of random points") + parser.add_option("-r", "--random-dest", action = 'store_true', default = False, help = "verbosity") + parser.add_option("-v", "--verbose", action = 'store_true', default = False, help = "verbosity") + + (options, args) = parser.parse_args() + print >> sys.stderr, options + + mesh_source = simple_rect_grid(options.source, options.source) + tree = scipy.spatial.KDTree(mesh_source.points) + + mesh_dest = simple_rect_grid(options.dest, options.dest) + + if options.random_dest: + x = [] + for i in xrange(options.random_total): + rx = np.random.rand() * 2 - 1 + ry = np.random.rand() * 2 - 1 + x.append([rx, ry]) + mesh_dest.points = np.array(x) + + errors = [] + for x in mesh_dest.points: + (final, exact) = baker.run_baker(x, mesh_source, tree, options.extra, options.verbose) + cur_error = np.abs(final - exact) + errors.append(cur_error) + + print rms(errors) diff --git a/bin/generate_qhull.pl b/bin/generate_qhull.pl new file mode 100755 index 0000000..437a6a9 --- /dev/null +++ b/bin/generate_qhull.pl @@ -0,0 +1,16 @@ +#!/usr/bin/perl + +use strict; +use warnings; + +my $output_filename = $ARGV[0] or die "no args"; +my $tmp_grid_file = "$output_filename.grid.txt"; + +system("grid.py > $tmp_grid_file"); + + +system("cat $tmp_grid_file | qdelaunay GD2 s > $output_filename.txt\n" ); +system("cat $tmp_grid_file | qdelaunay GD2 s Qt > ${output_filename}_qt.txt\n"); +system("cat $tmp_grid_file | qdelaunay GD2 s QJ > ${output_filename}_qj.txt\n"); + +unlink($tmp_grid_file); diff --git a/data/mesh.shelve b/data/mesh.shelve new file mode 100644 index 0000000000000000000000000000000000000000..bf3a15579b3446a026b54f61f3edfb92d49ed731 GIT binary patch literal 12288 zcmeI%OAdlC5P;zU6CZT}8}Hyo1HL^1Hy|!d*l9wLM6coDJcVa)XG>cnhNa8DNgvZ0 zhUD8;L?qDDjAGK4p*h5QHV2W=#6M$@$I{dJ=O`Z^N}Du2d>>!_xSx()_dd^#00Iag zfB*srAbI|2wGfB*srAbZw5(I5mxan^sVXu(3mq<9nAmXD&aPFa6U_&v hcw<9XyS0&Pv{R<#r*v;aSDGd^G_~`fzMQMu^$meoFcts+ literal 0 HcmV?d00001 diff --git a/data/qhull_input.txt b/data/qhull_input.txt new file mode 100644 index 0000000..4985325 --- /dev/null +++ b/data/qhull_input.txt @@ -0,0 +1,11 @@ +2 +9 +-1.0000 -1.0000 +-1.0000 0.0000 +-1.0000 1.0000 + 0.0000 -1.0000 + 0.0000 0.0000 + 0.0000 1.0000 + 1.0000 -1.0000 + 1.0000 0.0000 + 1.0000 1.0000 diff --git a/lib/baker.py b/lib/baker.py new file mode 100644 index 0000000..5e2e29a --- /dev/null +++ b/lib/baker.py @@ -0,0 +1,106 @@ +from grid import exact_func +import numpy as np +import sys + +def get_phis(X, r): + """ + The get_phis function is used to get barycentric coordonites for a point on a triangle. + + X -- the destination point (2D) + X = [0,0] + r -- the three points that make up the triangle (2D) + r = [[-1, -1], [0, 2], [1, -1]] + + this will return [0.333, 0.333, 0.333] + """ + + # baker: eq 7 + 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]]) + try: + phi = np.linalg.solve(A,b) + except: + print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system" + phi = np.dot(np.linalg.pinv(A), b) + + return phi + +def qlinear(X, r, q): + """ + this calculates the linear portion of q from X to r + + X = destination point + r = simplex points + q = CFD quantities of interest at the simplex points + """ + + phis = get_phis(X, r) + qlin = sum([q_i * phi_i for q_i, phi_i in zip(q[:len(phis)], phis)]) + return qlin + +def run_baker(X, g, tree, extra_points = 3, verbose = False): + """ + This is the main function to call to get an interpolation to X from the tree + + X -- the destination point (2D) + X = [0,0] + + g -- the grid object + + tree -- the kdtree search object (built from the g mesh) + """ + + (dist, indicies) = tree.query(X, 3 + extra_points) + + nn = [g.points[i] for i in indicies] + nq = [g.q[i] for i in indicies] + + phi = get_phis(X, nn[:3]) + qlin = nq[0] * phi[0] + nq[1] * phi[1] + nq[2] * phi[2] + + error_term = 0.0 + + if extra_points != 0: + B = [] # baker eq 9 + w = [] # baker eq 11 + + for index in indicies[3:]: + (phi1,phi2,phi3) = get_phis(g.points[index], nn) + B.append([phi1 * phi2, phi2*phi3, phi3*phi1]) + + w.append(g.q[index] - qlinear(g.points[index], nn, nq)) + + B = np.array(B) + w = np.array(w) + + A = np.dot(B.T, B) + b = np.dot(B.T, w) + + # baker solve eq 10 + try: + (a, b, c) = np.linalg.solve(A,b) + except: + print >> sys.stderr, "warning: linear calculation went bad, resorting to np.linalg.pinv" + (a, b, c) = np.dot(np.linalg.pinv(A), b) + + error_term = a * phi[0] * phi[1]\ + + b * phi[1] * phi[2]\ + + c * phi[2] * phi[0] + + exact = exact_func(X[0], X[1]) + q_final = qlin + error_term + + if verbose: + print "current point : %s" % X + print "exact : %0.4f" % exact + print "qlin : %0.4f" % qlin + print "qlinerr : %0.4f" % np.abs(exact - qlin) + print "q_final : %0.4f" % q_final + print "q_final_err : %0.4f" % np.abs(exact - q_final) + print + + return (q_final, exact) diff --git a/lib/grid.py b/lib/grid.py new file mode 100755 index 0000000..12aa0c5 --- /dev/null +++ b/lib/grid.py @@ -0,0 +1,77 @@ +#!/usr/bin/python + +import sys +import numpy as np +import scipy.spatial + +def exact_func(x, y): + return np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) + return np.sin(x * np.pi) * np.cos(y * np.pi) + + +class grid(object): + def __init__(self, points, q): + self.points = np.array(points) + self.q = np.array(q) + + def __str__(self): + r = '' + assert( len(self.points) == len(self.q) ) + for i in xrange(len(self.points)): + r += "%r: %0.4f\n" % ( self.points[i], self.q[i] ) + return r + +class simple_rect_grid(grid): + def __init__(self, xres = 5, yres = 5): + xmin = -1.0 + xmax = 1.0 + xspan = xmax - xmin + xdel = xspan / float(xres - 1) + + ymin = -1.0 + ymay = 1.0 + yspan = ymay - ymin + ydel = yspan / float(yres - 1) + + + self.points = [] + self.q = [] + for x in xrange(xres): + cur_x = xmin + (x * xdel) + for y in xrange(yres): + cur_y = ymin + (y * ydel) + self.points.append([cur_x, cur_y]) + self.q.append(exact_func(cur_x, cur_y)) + self.points = np.array(self.points) + self.q = np.array(self.q) + + + def for_qhull(self): + r = '2\n' + r += '%d\n' % len(self.points) + for p in self.points: + r += "%f %f\n" % (p[0], p[1]) + return r + + +class simple_random_grid(simple_rect_grid): + def __init__(self, num_points = 10): + self.points = [] + self.q = [] + + r = np.random + + for i in xrange(num_points): + cur_x = r.rand() + cur_y = r.rand() + + self.points.append([cur_x, cur_y]) + self.q.append(exact_func(cur_x, cur_y)) + + self.points = np.array(self.points) + self.q = np.array(self.q) + + +if __name__ == '__main__': + g = simple_random_grid(100) + print g.for_qhull() diff --git a/lib/smcqhull.py b/lib/smcqhull.py new file mode 100644 index 0000000..a50160d --- /dev/null +++ b/lib/smcqhull.py @@ -0,0 +1,4 @@ +#!/usr/bin/python + +if __name__ == '__main__': + print "hello world" diff --git a/lib/tools.py b/lib/tools.py new file mode 100644 index 0000000..65c3b66 --- /dev/null +++ b/lib/tools.py @@ -0,0 +1,12 @@ +import numpy as np + + +def rms(errors): + """ + root mean square calculation + """ + r = 0.0 + for i in errors: + r += np.power(i, 2) + r = np.sqrt(r / len(errors)) + return r diff --git a/test/baker.test.py b/test/baker.test.py new file mode 100755 index 0000000..bd512b5 --- /dev/null +++ b/test/baker.test.py @@ -0,0 +1,75 @@ +#!/usr/bin/python + +import unittest +import baker + +import numpy as np +import scipy.spatial + +class TestSequenceFunctions(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.approx_fmt = "%0.6f" + + def testGetPhis(self): + + X = [0,0] + r = [[-1, -1], [0, 2], [1, -1]] + + result = baker.get_phis(X, r) + result = [self.approx_fmt % i for i in result] + + right_answer = [self.approx_fmt % i for i in [1/3.0, 1/3.0, 1/3.0]] + + for a,b in zip(result, right_answer): + self.assertEqual(a,b) + + def testGetPhis2(self): + + X = [0.5,0.25] + r = [[0, 0], [1, 0], [1, 1]] + + result = baker.get_phis(X, r) + + right_answer = [0.5, 0.25, 0.25] + + 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] + + result = baker.qlinear(X, r, q) + + right_answer = 0.5 + + self.assertEqual(result, right_answer) + + def testRunBaker(self): + X = [0.5, 0.25] + + 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 + ] + q = [1, 0, 0, 0, 0, 0, 0, 0, 0] + + tree = scipy.spatial.KDTree(all_points) + baker.run_baker(X, + result = 3 + right_answer = 5 + + self.assertEqual(result, right_answer) + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions) + unittest.TextTestRunner(verbosity=2).run(suite) diff --git a/test/qhull.test.py b/test/qhull.test.py new file mode 100755 index 0000000..bd0b0a1 --- /dev/null +++ b/test/qhull.test.py @@ -0,0 +1,33 @@ +#!/usr/bin/python + +import delaunay +import subprocess + +QFILE = '/tmp/forQhull.txt' + +l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]] + +qdelauney_file = open(QFILE, 'w') + +qdelauney_file.write("%d\n" % len(l[0])) +qdelauney_file.write("%d\n" % len(l)) +for i in l: + qdelauney_file.write("%0.3f %0.3f\n" % (i[0], i[1])) + +dt = delaunay.Triangulation(l) +print dt.indices + +answer = [ + [4,1,3], + [1,5,0], + [5,1,4], + [7,3,6], + [7,4,3], + [7,5,4], + [5,7,8], + ] + +print answer == dt.indices +print dt.indices + +print "now run /usr/bin/qdelaunay Qt i < %s" % QFILE diff --git a/test/utest.py b/test/utest.py new file mode 100755 index 0000000..5a12ba0 --- /dev/null +++ b/test/utest.py @@ -0,0 +1,29 @@ +#!/usr/bin/python + + +import random +import unittest +import delaunay + +class TestSequenceFunctions(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]] + + + def testQhull(self): + dt = delaunay.Triangulation(self.l) + answer = [ + [4,1,3], + [1,5,0], + [5,1,4], + [7,3,6], + [7,4,3], + [7,5,4], + [5,7,8], + ] + + self.assertEqual(dt.indices, answer) + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions) + unittest.TextTestRunner(verbosity=5).run(suite) diff --git a/tools/mkpoints.py b/tools/mkpoints.py new file mode 100755 index 0000000..f670a8b --- /dev/null +++ b/tools/mkpoints.py @@ -0,0 +1,26 @@ +#!/usr/bin/python + +import sqlite3 +import sys, os +import numpy as np + +output_file = 'points.db' + +if os.path.exists(output_file): + print "found old db, erasing" + os.unlink(output_file) + +create = 'CREATE TABLE points (point_id integer primary key, x float, y float)' + +sqconn = sqlite3.connect(output_file) +sqc = sqconn.cursor() +sqc.execute(create) + +for i in xrange(int(1e6)): + rx = np.random.rand() * 2 - 1 + ry = np.random.rand() * 2 - 1 + sqc.execute('INSERT INTO points VALUES (NULL, ?, ?)', (rx, ry)) + + +sqconn.commit() +sqconn.close() diff --git a/tools/multi/forktest.py b/tools/multi/forktest.py new file mode 100755 index 0000000..9189ba0 --- /dev/null +++ b/tools/multi/forktest.py @@ -0,0 +1,22 @@ +#!/usr/bin/python + +import os +import time + +def work(): + j = 0.0 + for i in xrange(10000000): + j = j + j/2.0 + + + +pid = os.fork() +if pid != 0: + print "i spawned %d and am working" % pid + work() +else: + print "i am spawn, and am working" + work() + os._exit(0) + +print "parent done" diff --git a/tools/multi/pool.py b/tools/multi/pool.py new file mode 100755 index 0000000..e910798 --- /dev/null +++ b/tools/multi/pool.py @@ -0,0 +1,12 @@ +#!/usr/bin/python + +from multiprocessing import Pool + +def f(x): + return x*x + +if __name__ == '__main__': + pool = Pool(processes = 4) + result = pool.apply_async(f, (10,)) + print result.get() + print pool.map(f, range(10)) diff --git a/tools/multi/smp.py b/tools/multi/smp.py new file mode 100755 index 0000000..013126c --- /dev/null +++ b/tools/multi/smp.py @@ -0,0 +1,24 @@ +#!/usr/bin/python + +from multiprocessing import Process, Lock + +def f(l, i, d): + d['a'] = 'shutup' + print 'hello world', i, d + j = 0.0 + for i in xrange(1000000): + j = j + j/2.0 + +if __name__ == '__main__': + lock = Lock() + + d = {'a': 1, 'b': 2} + + ps = [] + for num in range(2): + ps.append(Process(target=f, args=(lock, num, d))) + + for p in ps: p.start() + for p in ps: p.join() + + print d diff --git a/tools/multi/threadtest.py b/tools/multi/threadtest.py new file mode 100755 index 0000000..fb2eab3 --- /dev/null +++ b/tools/multi/threadtest.py @@ -0,0 +1,22 @@ +#!/usr/bin/python + +import thread, time + +def work(): + j = 0.0 + for i in xrange(10000000): + j = j + j/2.0 + +def child(tid): + print "%d start" % tid + work() + print "%d stop" % tid + +def parent(): + i = 0 + while True: + i += 1 + thread.start_new(child, (i,)) + if raw_input() == 'q': break + +parent()