[Scipy-svn] r4659 - branches/interpolate/interpSNd
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Aug 20 14:54:39 EDT 2008
Author: fcady
Date: 2008-08-20 13:54:38 -0500 (Wed, 20 Aug 2008)
New Revision: 4659
Modified:
branches/interpolate/interpSNd/dewall.py
branches/interpolate/interpSNd/dewall.pyc
branches/interpolate/interpSNd/interp.py
branches/interpolate/interpSNd/interpolateSNd.py
branches/interpolate/interpSNd/script.py
branches/interpolate/interpSNd/test.py
Log:
fixed a bug in dewall. It hadn't alwayds been calculating the complete convex hull. Still not perfect, but appears to only miss interior points
Modified: branches/interpolate/interpSNd/dewall.py
===================================================================
--- branches/interpolate/interpSNd/dewall.py 2008-08-20 00:35:58 UTC (rev 4658)
+++ branches/interpolate/interpSNd/dewall.py 2008-08-20 18:54:38 UTC (rev 4659)
@@ -10,6 +10,10 @@
# In particular, calculation of the circumcircle is a purely
# mathematical operation that really should be made into C.
+# WARNING
+# This code doesn't always work. Occasionally a point that is
+# interior to the convex hull is missed.
+
import numpy as np
from numpy.linalg import norm
@@ -49,6 +53,7 @@
Sigma= []
alpha = select_alpha(P, AFL)
+ print "\nalpha:\n", alpha
# divide points into two sets separated by alpha
P1, P2 = pointset_partition(P, alpha) # both lists of points
@@ -58,6 +63,7 @@
if len(AFL) == 0: # This is only executed once, at the start of the algorithm
just_starting = True
first_simplex = make_first_simplex(P, alpha)
+ print "\nfirst simplex:\n",first_simplex
AFL = [ (face, get_out_vec(face,first_simplex))\
for face in faces(first_simplex)] # d+1 of them
Sigma.append(first_simplex)
@@ -71,6 +77,7 @@
if is_subset(face, P2):
AFL2.append((face,outvec))
while len(AFL_alpha) != 0:
+ print "\nAFL_alpha start:",[face for face, vec in AFL_alpha]
face, outvec = AFL_alpha.pop()
@@ -78,9 +85,10 @@
outward_points = filter( lambda p: (np.dot(p,outvec)>np.dot(face[0],outvec)),\
P)
else:
- outward_points = []#filter( lambda p: np.all([not point_in_list(p,vertex) for vertex in face]), P)
+ outward_points = filter( lambda p: np.all([not point_in_list(p,vertex) for vertex in face]), P)
t = make_simplex(face, outward_points) # make only over outer half space
+ print "\nnew simplex:\n",t
if t is not None:
Sigma.append(t)
@@ -96,6 +104,7 @@
if is_subset(f0, P2):
AFL2 = update(new_pair, AFL2)
+ print "\nAFL_alpha end:",[face for face, vec in AFL_alpha]
# now Sigma contains all simplices that intersect alpha
# Recursive Triangulation
@@ -161,14 +170,20 @@
# both are lists of arrays
return np.alltrue([ point_in_list(s1, S2) for s1 in S1])
-def update (face_pair, face_pair_list):
+def update (face_pair, face_pair_list): # this func has been problematic
# returns face_list with face_pair added if it wasn't there, else deleted
face, outvec = face_pair
- face_list = [face for face, outvec in face_pair_list]
+ face=face_pair[0]
+ print "face_pair: ", face_pair
+ face_list = [Face for Face, outvec in face_pair_list]
+ print "face: ", face
+ print "face_list: ", face_list
if face_in_list(face, face_list):
+ print "face in list"
f_not_equal_face = lambda face_pair : not np.alltrue([ point_in_list(p, face_pair[0]) for p in face ])
face_pair_list = filter(f_not_equal_face, face_pair_list)
else:
+ print "face not in list"
face_pair_list.append(face_pair)
return face_pair_list
Modified: branches/interpolate/interpSNd/dewall.pyc
===================================================================
(Binary files differ)
Modified: branches/interpolate/interpSNd/interp.py
===================================================================
--- branches/interpolate/interpSNd/interp.py 2008-08-20 00:35:58 UTC (rev 4658)
+++ branches/interpolate/interpSNd/interp.py 2008-08-20 18:54:38 UTC (rev 4659)
@@ -4,7 +4,7 @@
reload(SN)
-if False:
+if True:
points = np.array([[ 1.,0.,1.,0.],[0.,0.,1.,1.]])
z = np.array([1.,0.,2.,1.])
interp = SN.InterpolateSNd(points,z)
Modified: branches/interpolate/interpSNd/interpolateSNd.py
===================================================================
--- branches/interpolate/interpSNd/interpolateSNd.py 2008-08-20 00:35:58 UTC (rev 4658)
+++ branches/interpolate/interpSNd/interpolateSNd.py 2008-08-20 18:54:38 UTC (rev 4659)
@@ -3,6 +3,7 @@
import dewall as dw
import numpy as np
+from numpy.linalg import norm
class InterpolateSNd:
""" Interpolation of scatter data in arbitrarily many dimensions
@@ -135,7 +136,7 @@
assert vals.shape == (self.ndim+1,), \
"vals wrong shape: "+str(vals.shape)+"need: "+str((self.ndim,))
- weights = self.calculate_weights(simplex, point)
+ weights = self.calculate_weights(simplex, point) # None if point is not weighted average
return np.dot(np.ravel(weights), np.ravel(vals))
def calculate_weights(self, simplex, points):
@@ -143,6 +144,10 @@
of columns in simplex. Returns matrix where
jth column gives these weights for the jth column
of points.
+
+ If columns of simplex don't span the space, returns None.
+ FIXME : return correct weight if columns don't span space,
+ so long as point is still a weighted average
"""
assert simplex.shape == (self.ndim, self.ndim+1), "simplex shape: "+str(simplex.shape)
d, V = simplex.shape
@@ -159,30 +164,43 @@
axis = 0
)
- weights_vecs = np.linalg.solve(matrix_to_solve, vec_to_solve)
+ if rank(matrix_to_solve) < d:
+ print "matrix rank: "+str(np.rank(matrix_to_solve))
+ print "needed dim: "+str(d)
+ weights_vec = None
+ else:
+ weights_vec = np.linalg.solve(matrix_to_solve, vec_to_solve)
- return weights_vecs
+ return weights_vec
-def point_is_in_simplex(point, simplex):
- # point = array
- # simplex = matrix
- #assert point.shape == (self.ndim, 1), "wrong shape of point: "+str(point.shape)
- #assert simplex.shape == (self.ndim, self.ndim+1), "wrong shape of simplex: "+str(simplex.shape)
- weights_vec = calculate_weights(simplex, point)
- print "weights:\n", weights_vec
- weight_in_range = (0.<= weights_vec) & \
- (weights_vec <= 1.)
- print "in_range:\n", weight_in_range
- return np.alltrue(weight_in_range, axis=0)
+def rank(matrix):
+ matrix = matrix.copy()
+ n, m = matrix.shape
+ # Graham-Schmidt orthonormalization of input vectors
+ for j in range(m):
+ for k in range(j):
+ matrix[:,j] = matrix[:,j] - np.dot(matrix[:,k],matrix[:,j])*matrix[:,k]
+ if norm(matrix[:,j]) > dw.eps:
+ matrix[:,j]=matrix[:,j]/norm(matrix[:,j])
+
+ norm_of_each_vec = np.sum( matrix*matrix , axis=0)
+
+ return int(np.sum(norm_of_each_vec))
+
def calculate_weights(simplex, points):
""" Each column in points is a weighted average
of columns in simplex. Returns matrix where
jth column gives these weights for the jth column
of points.
+
+ If columns of simplex don't span the space, returns None.
+ FIXME : return correct weight if columns don't span space,
+ so long as point is still a weighted average
"""
- N, V = simplex.shape
- N, P = points.shape
+ #assert simplex.shape == (self.ndim, self.ndim+1), "simplex shape: "+str(simplex.shape)
+ d, V = simplex.shape
+ d, P = points.shape
matrix_to_solve = np.concatenate((simplex, \
np.ones((1,V))
@@ -195,6 +213,11 @@
axis = 0
)
- weights_vecs = np.linalg.solve(matrix_to_solve, vec_to_solve)
+ if np.rank(matrix_to_solve) < d:
+ print "matrix rank: "+str(rank(matrix_to_solve))
+ print "needed dim: "+str(d)
+ weights_vec = None
+ else:
+ weights_vec = np.linalg.solve(matrix_to_solve, vec_to_solve)
- return weights_vecs
\ No newline at end of file
+ return weights_vec
\ No newline at end of file
Modified: branches/interpolate/interpSNd/script.py
===================================================================
--- branches/interpolate/interpSNd/script.py 2008-08-20 00:35:58 UTC (rev 4658)
+++ branches/interpolate/interpSNd/script.py 2008-08-20 18:54:38 UTC (rev 4659)
@@ -23,12 +23,12 @@
y = center[1] + y_offset
pyplot.plot(x,y)
-Pb=[array([.25, -.25]), array([0,.75])]
+Pb=[]#[array([.25, -.25]), array([0,.75])]
P = Pa+Pb
P = [ np.array([np.random.gamma(1), np.random.gamma(1)]) \
- for j in range(6) ]
+ for j in range(10) ]
triangul = dw.dewall(P)
@@ -42,8 +42,8 @@
# plotting the circumcircles
for tri in triangul:
- plot_circle(dw.circumcircle(tri))
+ pass#plot_circle(dw.circumcircle(tri))
pyplot.show()
-print triangul
\ No newline at end of file
+print "triangulation:\n",triangul
\ No newline at end of file
Modified: branches/interpolate/interpSNd/test.py
===================================================================
--- branches/interpolate/interpSNd/test.py 2008-08-20 00:35:58 UTC (rev 4658)
+++ branches/interpolate/interpSNd/test.py 2008-08-20 18:54:38 UTC (rev 4659)
@@ -9,61 +9,59 @@
class Test(unittest.TestCase):
def compare_arrays(self, a, b):
- return np.allclose(a,b)
+ return np.allclose(a,b,rtol=1e-3) | (np.isnan(a)&np.isnan(b))
## test Delaunay triangulation itself
# testing pathological cases with
- def test_square(self):
+ def _test_square(self):
P = [array([0.,1.]), array([0.,0.]), array([1.,0.]), array([1.,1.])]
tri = dw.dewall(P)
self.assert_( len(tri)==2 )
self.assert_( len(tri[0])==3 )
self.assert_( len(tri[1])==3 )
- def test_linear(self):
+ def _test_linear(self):
P = [array([0.,1.]), array([0.,0.]), array([0.,-1.])]
tri = dw.dewall(P)
# testing general case using random data
- def test_2d(self):
- print "TESTING 2D"
- P = [np.random.random_sample(2) for i in range(15)]
+ def _test_2d(self):
+ ndim = 2
+ nold = 15
+ print "TESTING %iD"%ndim
+ corner_points_matrix = corners(ndim)
+ corner_points = [corner_points_matrix[:,i] for i in range(2**ndim)]
+ interior_points = [np.random.random_sample(ndim) for i in range(nold)]
+ P = corner_points+interior_points
+
+ # not checking if its correct, just if it runs.
tri = dw.dewall(P)
- # not checking if its correct, just if it runs.
+
+ # making sure it's correct
- def test_3d(self):
+ def _test_3d(self):
print "TESTING 3D"
P = [np.random.random_sample(3) for i in range(9)]
tri = dw.dewall(P)
# not checking if its correct, just if it runs.
- def test_4d(self):
+ def _test_4d(self):
print "TESTING 4D"
P = [np.random.random_sample(4) for i in range(9)]
tri = dw.dewall(P)
# not checking if its correct, just if it runs.
- def test_5d(self):
+ def _test_5d(self):
P = [np.random.random_sample(5) for i in range(9)]
tri = dw.dewall(P)
# not checking if its correct, just if it runs.
## test interpolation, and thus also triangulation by extension
- def test_2d(self):
- points = np.array([[ 1.,0.,1.,0.],[0.,0.,1.,1.]])
- z = np.array([1.,0.,2.,1.])
- interp = SNd.InterpolateSNd(points,z)
-
- X=np.array([[.4,.1,.55],[.4,.1,.3]])
-
- output = interp(X)
- self.assert_(self.compare_arrays(output, array([.8,.2,.85])))
-
- def test_linear_on_cube(self):
+ def _test_linear_on_cube(self):
x = array([0., 1, 0, 1, 0, 1, 0, 1])
y = array([0., 0, 1, 1, 0, 0, 1, 1])
z = array([0., 0, 0, 0, 1, 1, 1, 1])
@@ -78,18 +76,49 @@
self.assert_(self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)))
- def test_linear_5d(self):
- P = [np.random.random_sample(5) for i in range(9)]
- points = array(P).reshape((5,9))
- fvals = points[:,0]+points[:,1]+points[:,2]+points[:,3]+points[:,4]
+ def test_linear_2d(self):
+ ndim = 2 # num dimensions
+ nold = 5 # num known data points
+ nnew = 5 # num points at which to interpolate
+ print "%iD Interpolation"%ndim
+
+ P = [np.random.random_sample(ndim) for i in range(nold)]
+ # points at corners of hypercube and radnimly scattered in the interior
+ points = np.concatenate((corners(ndim) , array(P).reshape((ndim,nold))), axis=1)
+ fvals = np.zeros((1,points.shape[1]))
+ for i in range(ndim):
+ fvals = fvals+points[i,:]
+ fvals = fvals.reshape((points.shape[1]))
+
+ print "points:\n",points
+ print "fvals:\n",fvals
+
interp = SNd.InterpolateSNd(points, fvals)
- newdata = np.random.random_sample((5,8))
+ print "\ntriang:"
+ for x in interp._triangulation: print x
+
+ newdata = np.random.random_sample((ndim,nnew))
+ print "\nnewdata:\n",newdata
interpvals = interp(newdata)
realvals = np.sum(newdata, axis=0)
- self.assert_(self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)))
+ print "%iD interpvals: "%ndim, interpvals
+ print "%iD realvals: "%ndim, realvals
+ #self.assert_(self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)))
+ assert self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)), "wrong data"
+
+def corners(ndim):
+ # returns matrix indicating corners of unit cube
+ result = np.zeros((ndim,2**ndim))
+ for i in range(ndim):
+ # each block is 2**i spaces wide
+ for j in range(2**ndim):
+ if ( j%(2**(i+1)) )/(2**i) == 1: result[i,j]=1.
+ return result
+
+
if __name__ == "__main__":
unittest.main()
\ No newline at end of file
More information about the Scipy-svn
mailing list