From 3dcc10ea0e9707c471b7de9896619e00013363ba Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Sat, 24 Apr 2010 18:26:44 -0600 Subject: [PATCH] working on the containing simplex problem. the get containing simplex function works, if the point is in the domain. I am trying to force the domain to contain the point (placing points along perimiter of the mesh). I am not doing that right, but have to run --HG-- rename : bin/qhull-029.txt => data/qhull-029.txt rename : bin/qhull-029.txt.gv => data/qhull-029.txt.gv --- bin/test.py | 3 +- bin/test2d-connectivity.py | 32 ++- bin/test2d-qhull-specific.py | 3 +- data/qhull-000.txt | 502 +++++++++++++++++++++++++++++++++ {bin => data}/qhull-029.txt | 0 {bin => data}/qhull-029.txt.gv | 0 lib/baker/__init__.py | 49 ++-- lib/baker/tools.py | 10 + lib/grid/DD.py | 9 + lib/grid/__init__.py | 52 +++- 10 files changed, 604 insertions(+), 56 deletions(-) create mode 100644 data/qhull-000.txt rename {bin => data}/qhull-029.txt (100%) rename {bin => data}/qhull-029.txt.gv (100%) diff --git a/bin/test.py b/bin/test.py index 66c7851..e05aa47 100755 --- a/bin/test.py +++ b/bin/test.py @@ -3,6 +3,7 @@ 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 try: total_points = int(sys.argv[1]) @@ -10,7 +11,7 @@ except: total_points = 20 -print total_points +log.info(total_points) g = random_grid(total_points) open('/tmp/for_qhull.txt', 'w').write(g.for_qhull()) diff --git a/bin/test2d-connectivity.py b/bin/test2d-connectivity.py index 121c74d..3ffbfb8 100755 --- a/bin/test2d-connectivity.py +++ b/bin/test2d-connectivity.py @@ -7,16 +7,19 @@ from baker.tools import exact_func, smberror, improved_answer from glob import glob from os import remove +import numpy as np + FILE_PREFIX='/tmp/qhull-' for g in glob('%s*' % FILE_PREFIX): remove(g) try: - total_points = int(sys.argv[1]) - total_tries = int(sys.argv[2]) + total_points = int(sys.argv[1]) + random_points = int(sys.argv[2]) + total_tries = int(sys.argv[3]) except: - print "usage: app.py total_points total_tries" + print "usage: app.py [total points in random grid] [number of random points to attempt] [total attempts]" sys.exit(1) @@ -25,18 +28,19 @@ print "total points", total_points d = {True: 0, False: 0} -for i in xrange(total_tries): +for cur_try in xrange(total_tries): g = random_grid(total_points) + open('%s%0.3d.txt' % (FILE_PREFIX, cur_try), 'w').write(g.for_qhull()) - open('%s%0.3d.txt' % (FILE_PREFIX, i), 'w').write(g.for_qhull()) + for i in xrange(random_points): + X = [np.random.rand(), np.random.rand()] + exact = exact_func(X) - X = [0.3, 0.3] - exact = exact_func(X) + try: + answer = g.run_baker(X) + d[improved_answer(answer, exact)] += 1 + except smberror as e: + print e + print d - try: - answer = g.run_baker(X) - d[improved_answer(answer, exact)] += 1 - except smberror as e: - print e - -print d + print d diff --git a/bin/test2d-qhull-specific.py b/bin/test2d-qhull-specific.py index fe7a958..c62ac6b 100755 --- a/bin/test2d-qhull-specific.py +++ b/bin/test2d-qhull-specific.py @@ -24,9 +24,8 @@ q = [exact_func(i) for i in verts] g = grid(verts,q) -X = [0.3, 0.3] +X = [0.13206472640187261, 0.98644929172054607] -g.get_containing_simplex(X) exact = exact_func(X) try: diff --git a/data/qhull-000.txt b/data/qhull-000.txt new file mode 100644 index 0000000..8e77127 --- /dev/null +++ b/data/qhull-000.txt @@ -0,0 +1,502 @@ +2 +500 +0.501822 0.099459 +0.398871 0.145078 +0.717727 0.819884 +0.218163 0.456135 +0.692657 0.874376 +0.698942 0.726672 +0.514508 0.713024 +0.635646 0.506584 +0.838798 0.506751 +0.090156 0.534329 +0.763002 0.330018 +0.325797 0.177421 +0.929666 0.349459 +0.903110 0.758269 +0.788727 0.796974 +0.827882 0.283563 +0.932292 0.840258 +0.730115 0.024117 +0.050927 0.682638 +0.340828 0.344778 +0.265962 0.860156 +0.254063 0.747634 +0.066011 0.086948 +0.124991 0.722423 +0.774609 0.825563 +0.743862 0.624939 +0.587882 0.527069 +0.003707 0.400760 +0.606581 0.508782 +0.733931 0.823112 +0.875321 0.387433 +0.421334 0.179495 +0.487737 0.244955 +0.986045 0.549714 +0.803826 0.122495 +0.513748 0.295921 +0.862030 0.670105 +0.855385 0.848034 +0.254363 0.034489 +0.244426 0.550130 +0.039767 0.946977 +0.324643 0.651447 +0.325729 0.827348 +0.209335 0.544249 +0.579131 0.622810 +0.174603 0.157842 +0.872530 0.154341 +0.588218 0.573515 +0.259025 0.733801 +0.573148 0.073529 +0.781172 0.728942 +0.995370 0.779436 +0.002421 0.341314 +0.288521 0.965963 +0.947036 0.016416 +0.548763 0.242006 +0.138231 0.932469 +0.695652 0.832567 +0.553380 0.794949 +0.495854 0.870735 +0.336602 0.132771 +0.682331 0.430044 +0.475190 0.954946 +0.294329 0.497394 +0.180446 0.687540 +0.627226 0.322490 +0.107537 0.650209 +0.325204 0.625028 +0.883858 0.984149 +0.472804 0.533670 +0.906653 0.295008 +0.463240 0.246410 +0.313726 0.493177 +0.613654 0.141180 +0.567792 0.034890 +0.294997 0.296258 +0.432614 0.445397 +0.656278 0.341600 +0.358028 0.681119 +0.322254 0.325132 +0.954686 0.952094 +0.219197 0.059499 +0.973052 0.038154 +0.079889 0.703179 +0.953068 0.369875 +0.781529 0.466225 +0.156998 0.618783 +0.504962 0.588302 +0.036981 0.711112 +0.225776 0.494560 +0.595356 0.032457 +0.287397 0.647716 +0.555859 0.800712 +0.103839 0.396361 +0.856306 0.169453 +0.899026 0.579813 +0.253680 0.732919 +0.934093 0.780840 +0.052220 0.267688 +0.088222 0.827587 +0.024444 0.578090 +0.771479 0.741172 +0.472561 0.162611 +0.830357 0.115662 +0.832320 0.364346 +0.853613 0.047028 +0.261261 0.744922 +0.635028 0.241084 +0.975904 0.320519 +0.247161 0.979842 +0.227195 0.664762 +0.414410 0.678235 +0.849322 0.057320 +0.889820 0.773301 +0.852028 0.219397 +0.961608 0.956863 +0.780950 0.161184 +0.984412 0.948117 +0.109250 0.003821 +0.828057 0.102844 +0.882761 0.392760 +0.482812 0.284440 +0.897279 0.752738 +0.783414 0.125775 +0.137902 0.965423 +0.443151 0.894608 +0.602202 0.120583 +0.902763 0.527426 +0.052316 0.440764 +0.181373 0.671437 +0.340350 0.361405 +0.761933 0.780775 +0.337599 0.242759 +0.672853 0.540019 +0.606246 0.113563 +0.666268 0.990161 +0.519951 0.759898 +0.464147 0.660421 +0.907478 0.106734 +0.959816 0.896417 +0.350590 0.283244 +0.844312 0.342034 +0.780929 0.692151 +0.937083 0.711286 +0.042409 0.901348 +0.830176 0.200075 +0.532453 0.391126 +0.737095 0.328319 +0.386647 0.824122 +0.770762 0.234036 +0.090658 0.794761 +0.542149 0.652688 +0.447155 0.403339 +0.393642 0.302263 +0.419523 0.642328 +0.164542 0.013951 +0.480443 0.505623 +0.824197 0.518905 +0.462652 0.864183 +0.442192 0.980663 +0.249635 0.157728 +0.236748 0.711805 +0.806916 0.851598 +0.968949 0.659096 +0.061502 0.498107 +0.197414 0.645504 +0.435848 0.921071 +0.758112 0.038039 +0.723839 0.350983 +0.745719 0.558800 +0.391832 0.771099 +0.015951 0.379046 +0.076760 0.583620 +0.779949 0.715885 +0.378045 0.961058 +0.208069 0.476748 +0.525158 0.050896 +0.429921 0.754940 +0.324172 0.825312 +0.338008 0.128663 +0.515633 0.738803 +0.287414 0.276497 +0.525973 0.061000 +0.975127 0.625654 +0.140107 0.734905 +0.671401 0.789625 +0.695956 0.098697 +0.930942 0.698302 +0.334685 0.738623 +0.603622 0.716263 +0.879529 0.306514 +0.115903 0.967280 +0.746742 0.843177 +0.713765 0.886153 +0.115388 0.333187 +0.626863 0.175279 +0.648830 0.067674 +0.322128 0.191630 +0.843209 0.929941 +0.113607 0.429130 +0.764912 0.949248 +0.853386 0.604863 +0.802790 0.303906 +0.077434 0.340228 +0.958669 0.742775 +0.659572 0.005397 +0.981930 0.649579 +0.638143 0.652816 +0.294420 0.963091 +0.692114 0.081922 +0.102524 0.420897 +0.852326 0.646664 +0.569550 0.608715 +0.109439 0.207973 +0.921732 0.071777 +0.480585 0.712810 +0.087866 0.027396 +0.482599 0.864175 +0.119039 0.355474 +0.560596 0.391998 +0.765300 0.861846 +0.755712 0.392778 +0.415602 0.117571 +0.111080 0.486659 +0.673324 0.986224 +0.707238 0.020572 +0.347888 0.484752 +0.770161 0.468600 +0.667128 0.574563 +0.296673 0.036848 +0.463483 0.049365 +0.157612 0.089058 +0.853492 0.621741 +0.354271 0.059681 +0.074932 0.665053 +0.576566 0.385005 +0.677528 0.653185 +0.975808 0.765341 +0.182598 0.836122 +0.934129 0.413901 +0.715217 0.526190 +0.524517 0.970114 +0.462338 0.496160 +0.172289 0.911521 +0.052674 0.448721 +0.130097 0.741582 +0.220472 0.096307 +0.170737 0.859070 +0.988636 0.191693 +0.303820 0.151587 +0.639676 0.628261 +0.422134 0.222226 +0.499111 0.458441 +0.832883 0.971754 +0.996332 0.280903 +0.404805 0.347382 +0.219804 0.815241 +0.872310 0.039327 +0.085328 0.117189 +0.034204 0.162739 +0.122386 0.387208 +0.221209 0.315134 +0.490247 0.692898 +0.244592 0.911100 +0.274404 0.544079 +0.897967 0.654600 +0.723768 0.098982 +0.215147 0.464956 +0.295471 0.199515 +0.785139 0.102539 +0.915978 0.678969 +0.485031 0.737332 +0.396784 0.937990 +0.295549 0.911359 +0.817506 0.697086 +0.558902 0.825836 +0.990089 0.854087 +0.680328 0.660951 +0.190960 0.255092 +0.227429 0.906142 +0.840200 0.232338 +0.412978 0.003388 +0.281234 0.702691 +0.606999 0.359344 +0.142862 0.809069 +0.095922 0.259706 +0.678389 0.998062 +0.085304 0.556810 +0.778608 0.391243 +0.376031 0.954876 +0.585315 0.829717 +0.632199 0.555069 +0.992517 0.210612 +0.851749 0.815774 +0.929535 0.214854 +0.737507 0.910058 +0.593620 0.223300 +0.483844 0.063215 +0.797046 0.674707 +0.002823 0.276132 +0.212396 0.760201 +0.226260 0.110551 +0.331908 0.002332 +0.683207 0.797607 +0.218008 0.399797 +0.811012 0.770548 +0.513188 0.444371 +0.547434 0.449839 +0.421245 0.608438 +0.294447 0.275174 +0.804497 0.425950 +0.475786 0.468890 +0.758350 0.751438 +0.947073 0.186788 +0.680574 0.273894 +0.816370 0.779451 +0.808707 0.992098 +0.597573 0.280983 +0.992554 0.351677 +0.296248 0.313174 +0.289828 0.714345 +0.904771 0.832176 +0.968310 0.490132 +0.870036 0.021450 +0.376132 0.170474 +0.377137 0.223095 +0.055540 0.426383 +0.147409 0.741300 +0.300452 0.108651 +0.032884 0.313330 +0.867730 0.276854 +0.484683 0.581226 +0.271740 0.187032 +0.802124 0.885425 +0.467229 0.098781 +0.925030 0.484433 +0.929984 0.772211 +0.349099 0.123479 +0.678556 0.241998 +0.584291 0.283861 +0.621535 0.769794 +0.712339 0.869903 +0.121415 0.362363 +0.139690 0.576022 +0.623300 0.284221 +0.874903 0.308943 +0.050968 0.315062 +0.415732 0.861086 +0.269109 0.632107 +0.915507 0.504078 +0.263057 0.955175 +0.080084 0.384505 +0.825620 0.145982 +0.696184 0.599655 +0.041665 0.753284 +0.935032 0.058134 +0.919558 0.320587 +0.996719 0.935324 +0.034179 0.175121 +0.853912 0.697859 +0.237573 0.737667 +0.267024 0.466481 +0.752191 0.412488 +0.573358 0.776042 +0.117738 0.008809 +0.391351 0.799168 +0.351743 0.448648 +0.662856 0.038363 +0.155716 0.234058 +0.138554 0.485890 +0.172528 0.261335 +0.667910 0.841149 +0.704714 0.702788 +0.326822 0.296431 +0.165563 0.696065 +0.535877 0.037066 +0.790105 0.399294 +0.863059 0.179309 +0.946292 0.349932 +0.190070 0.030322 +0.548990 0.765724 +0.189105 0.326289 +0.593151 0.625304 +0.080545 0.071429 +0.609873 0.315857 +0.065005 0.557161 +0.172772 0.854238 +0.236891 0.340104 +0.243522 0.490912 +0.521110 0.131008 +0.265345 0.106485 +0.556598 0.480571 +0.583463 0.318790 +0.422577 0.524559 +0.047355 0.832263 +0.015258 0.359675 +0.087080 0.749055 +0.586842 0.711002 +0.500708 0.511499 +0.906278 0.753350 +0.191065 0.742964 +0.850002 0.747699 +0.810984 0.958665 +0.745445 0.101771 +0.812433 0.356154 +0.349673 0.482697 +0.309605 0.448379 +0.906286 0.305023 +0.871195 0.236735 +0.880085 0.619723 +0.252215 0.218488 +0.501502 0.938435 +0.571883 0.839314 +0.491132 0.508797 +0.613232 0.418147 +0.838755 0.849814 +0.319751 0.694594 +0.102121 0.280414 +0.493146 0.520976 +0.080032 0.935737 +0.978841 0.426193 +0.572612 0.929880 +0.562613 0.788561 +0.162098 0.497810 +0.187087 0.439069 +0.080479 0.414313 +0.900542 0.203189 +0.746217 0.777058 +0.733827 0.872460 +0.735957 0.728743 +0.129533 0.625628 +0.579763 0.036779 +0.481897 0.045606 +0.909063 0.879236 +0.415542 0.769329 +0.550507 0.500253 +0.726148 0.975612 +0.000193 0.603004 +0.487489 0.695720 +0.954665 0.714453 +0.362283 0.292757 +0.948384 0.300804 +0.395100 0.822704 +0.135503 0.304895 +0.592481 0.260597 +0.021018 0.757291 +0.592562 0.782326 +0.138078 0.735031 +0.638480 0.457732 +0.587009 0.255635 +0.355952 0.747673 +0.721680 0.346664 +0.302100 0.822532 +0.451586 0.424316 +0.795754 0.368634 +0.074531 0.085938 +0.171169 0.610738 +0.169327 0.627517 +0.170210 0.139028 +0.716423 0.333147 +0.757550 0.343675 +0.547498 0.134499 +0.797433 0.285402 +0.056388 0.965889 +0.876258 0.041214 +0.940046 0.599479 +0.414035 0.393592 +0.240744 0.965824 +0.848756 0.776561 +0.699504 0.236643 +0.378730 0.003658 +0.659677 0.101798 +0.855023 0.541043 +0.206409 0.835201 +0.864257 0.742651 +0.507091 0.737711 +0.278080 0.176298 +0.439871 0.586918 +0.325005 0.662985 +0.898581 0.664883 +0.973856 0.807850 +0.767325 0.718386 +0.616289 0.183717 +0.002939 0.230201 +0.925781 0.241950 +0.500550 0.574635 +0.622247 0.624035 +0.260977 0.344729 +0.644911 0.765267 +0.253137 0.819781 +0.241926 0.709515 +0.359244 0.740632 +0.109408 0.105881 +0.306583 0.231504 +0.197410 0.586292 +0.066206 0.518310 +0.375724 0.151438 +0.459500 0.107776 +0.088232 0.245200 +0.558558 0.076231 diff --git a/bin/qhull-029.txt b/data/qhull-029.txt similarity index 100% rename from bin/qhull-029.txt rename to data/qhull-029.txt diff --git a/bin/qhull-029.txt.gv b/data/qhull-029.txt.gv similarity index 100% rename from bin/qhull-029.txt.gv rename to data/qhull-029.txt.gv diff --git a/lib/baker/__init__.py b/lib/baker/__init__.py index e59c2a0..48b092b 100644 --- a/lib/baker/__init__.py +++ b/lib/baker/__init__.py @@ -1,4 +1,5 @@ from baker import * +from baker.tools import logging as log import numpy as np import sys @@ -30,7 +31,7 @@ def get_phis(X, R): 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 - # TODO: log this -- > print >> sys.stderr, msg + log.error(msg) raise smberror(msg) phi = np.dot(np.linalg.pinv(A), b) @@ -136,7 +137,7 @@ def get_error_quadratic(phi, R, S): + b * phi[1] * phi[2]\ + c * phi[2] * phi[0] - return error_term, a, b, c + return error_term def get_error_cubic(phi, R, S): B = [] # baker eq 9 @@ -181,7 +182,7 @@ def get_error_cubic(phi, R, S): + f * phi[2] * phi[1] * phi[1]\ + g * phi[0] * phi[1] * phi[2]\ - return error_term, a, b, c + return error_term def run_baker(X, R, S, order=2): @@ -195,37 +196,29 @@ def run_baker(X, R, S, order=2): S = extra points """ + answer = { + 'qlin': None, + 'error': None, + 'final': None, + } # calculate values only for the simplex triangle phi, qlin = qlinear(X, R) - if len(S.points) == 0: - answer = { - 'a': None, - 'b': None, - 'c': None, - 'qlin': qlin, - 'error': None, - 'final': None, - } + if order == 1: + answer['qlin'] = qlin return answer - - if order == 2: - error_term, a, b, c = get_error_quadratic(phi, R, S) + elif order == 2: + error_term = get_error_quadratic(phi, R, S) elif order == 3: - error_term, a, b, c = get_error_cubic(phi, R, S) + error_term = get_error_cubic(phi, R, S) else: raise smberror('unacceptable order for baker method') q_final = qlin + error_term - answer = { - 'a': a, - 'b': b, - 'c': c, - 'qlin': qlin, - 'error': error_term, - 'final': q_final, - } + answer['qlin' ] = qlin + answer['error'] = error_term + answer['final'] = q_final return answer @@ -243,14 +236,6 @@ def run_baker_3D(X, R, S): # calculate values only for the triangle phi, qlin = qlinear_3D(X, R) - if [i for i in phi if i <= 0.0]: - s = "this is not a containing simplex:\n" - s += " X: %s\n" % X - s += " R: %s\n" % R - s += " phi: %s, sum(%0.4e)\n" % (phi, sum(phi)) - print >> sys.stderr, s - raise smberror("not containing simplex") - if len(S.points) == 0: answer = { 'a': None, diff --git a/lib/baker/tools.py b/lib/baker/tools.py index 172c3f6..0a8b7fb 100644 --- a/lib/baker/tools.py +++ b/lib/baker/tools.py @@ -1,5 +1,15 @@ +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 numpy as np + class smberror(Exception): """ this is a silly little exception subclass diff --git a/lib/grid/DD.py b/lib/grid/DD.py index 0959be2..fe73f74 100644 --- a/lib/grid/DD.py +++ b/lib/grid/DD.py @@ -61,6 +61,15 @@ class random_grid(rect_grid): q = [] r = np.random + appx_side_res = int(np.sqrt(num_points)) + 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 + new_point = [cur_x, cur_y] + points.append(new_point) + q.append(exact_func(new_point)) for i in xrange(num_points): cur_x = r.rand() diff --git a/lib/grid/__init__.py b/lib/grid/__init__.py index 037dede..1400858 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 +from baker.tools import exact_func, smberror, logging from simplex import face, contains from smcqdelaunay import * @@ -42,7 +42,7 @@ 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] @@ -50,7 +50,45 @@ class grid(object): return grid(p, q) def get_containing_simplex(self, X): - pass + if not self.faces: + logging.debug('get_containing_simplex: 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]) + + simplex = None + checked_facets = [] + facets_to_check = self.facets_for_point[closest_point] + + attempts = 0 + while not simplex: + logging.debug('attempt: %d' % attempts) + attempts += 1 + if attempts > 10: + 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]) + + if cur_facet.contains(X, self): + simplex = cur_facet + + if not simplex: + raise AssertionError('no containing simplex found') + + R = self.create_mesh(simplex.verts) + + logging.debug('this must be R: %s' % R) + return R def get_simplex_and_nearest_points(self, X, extra_points = 3, simplex_size = 3): @@ -64,9 +102,8 @@ class grid(object): """ (dist, indicies) = self.tree.query(X, simplex_size + extra_points) - - # get the containing simplex - r_mesh = self.create_mesh(indicies[:simplex_size]) + # r_mesh = self.create_mesh(indicies[:simplex_size]) + r_mesh = self.get_containing_simplex(X) # and some extra points s_mesh = self.create_mesh(indicies[simplex_size:]) @@ -103,7 +140,7 @@ class grid(object): if not simplex: raise AssertionError('no containing simplex found') - R = self.create_mesh(simplex.verts) + R = get_containing_simplex(X)# self.create_mesh(simplex.verts) s = [] @@ -136,6 +173,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') qdelaunay_string = get_qdelaunay_dump_str(self) facet_to_facets = [] for matcher in grid.facet_re.finditer(qdelaunay_string):