[Scipy-svn] r4656 - branches/interpolate/interpSNd
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Aug 19 19:39:20 EDT 2008
Author: fcady
Date: 2008-08-19 18:39:17 -0500 (Tue, 19 Aug 2008)
New Revision: 4656
Removed:
branches/interpolate/interpSNd/interpolateSNd.pyc
branches/interpolate/interpSNd/output.txt
Modified:
branches/interpolate/interpSNd/dewall.py
branches/interpolate/interpSNd/dewall.pyc
branches/interpolate/interpSNd/interp.py
branches/interpolate/interpSNd/interpolateSNd.py
branches/interpolate/interpSNd/test.py
Log:
fixed a bug in 3D and higher
Modified: branches/interpolate/interpSNd/dewall.py
===================================================================
--- branches/interpolate/interpSNd/dewall.py 2008-08-19 22:26:16 UTC (rev 4655)
+++ branches/interpolate/interpSNd/dewall.py 2008-08-19 23:39:17 UTC (rev 4656)
@@ -23,14 +23,13 @@
def dewall (P, #set of points
AFL = [], # list of faces: (d-1)face list
):
-
+
# checking input
assert isinstance(P, list)
if len(P)>0:
assert isinstance(P[0], np.ndarray)
assert isinstance(AFL, list)
if len(AFL)>0:
- #print "AFL[0]: ", AFL[0]
assert isinstance(AFL[0],tuple)
assert isinstance(AFL[0][0], list)
assert isinstance(AFL[0][0][0], np.ndarray)
@@ -53,10 +52,10 @@
# divide points into two sets separated by alpha
P1, P2 = pointset_partition(P, alpha) # both lists of points
-
+
# Simplex Wall Construction
just_starting = False #source of problem?
- if len(AFL) == 0:
+ 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)
AFL = [ (face, get_out_vec(face,first_simplex))\
@@ -79,9 +78,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
+
if t is not None:
Sigma.append(t)
# for f0 != f in faces(t) , ie for new outward faces
@@ -90,6 +90,7 @@
if is_intersected(f0, alpha):
# continue building wall out in that direction
AFL_alpha = update(new_pair, AFL_alpha)
+ #np.random.shuffle(AFL_alpha)
if is_subset(f0, P1):
AFL1 = update(new_pair, AFL1)
if is_subset(f0, P2):
@@ -109,16 +110,6 @@
# must intersect active face if there is one
# must not intersect any points
- #assert isinstance(P, list)
- #if len(P)>0:
- # assert isinstance(P[0], np.ndarray)
- #assert isinstance(AFL, list)
- #if len(AFL)>0:
- # assert isinstance(AFL[0],list)
- # assert isinstance(AFL[0][0], list)
- # assert isinstance(AFL[0][0][0], np.ndarray)
- # assert isinstance(AFL[0][1], np.ndarray)
-
d = len(P[0])
# plane through avg of cluster. Guarantees separation
@@ -174,8 +165,8 @@
# 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]
- if face_in_list(face, face_pair_list):
- f_not_equal_face = lambda f : not np.alltrue([ point_in_list(p, f[0]) for p in face ])
+ if face_in_list(face, face_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:
face_pair_list.append(face_pair)
Modified: branches/interpolate/interpSNd/dewall.pyc
===================================================================
(Binary files differ)
Modified: branches/interpolate/interpSNd/interp.py
===================================================================
--- branches/interpolate/interpSNd/interp.py 2008-08-19 22:26:16 UTC (rev 4655)
+++ branches/interpolate/interpSNd/interp.py 2008-08-19 23:39:17 UTC (rev 4656)
@@ -1,18 +1,40 @@
import interpolateSNd as SN
import numpy as np
+import dewall as dw
reload(SN)
-points = np.array([[ 1.,0.,1.,0.],[0.,0.,1.,1.]])
-z = np.array([1.,0.,2.,1.])
-interp = SN.InterpolateSNd(points,z)
+if False:
+ points = np.array([[ 1.,0.,1.,0.],[0.,0.,1.,1.]])
+ z = np.array([1.,0.,2.,1.])
+ interp = SN.InterpolateSNd(points,z)
-print "triangulation:\n",interp._triangulation
+ print "triangulation:\n",interp._triangulation
-X=np.array([[1.4,.1,.55],[.4,.1,.3]])
-print "X values:\n",X
-# last component is .1 too much
+ X=np.array([[1.4,.1,.55],[.4,.1,.3]])
+ print "X values:\n",X
-output = interp(X)
+ output = interp(X)
-print "output:\n", output
\ No newline at end of file
+ print "output:\n", output
+
+
+if True:
+ points = np.array([ [0., 0, 0, 1.,.2],
+ [0., 1., 0, 0, .2],
+ [0., 0, 1., 0, .2] ])
+ z = np.array([.1,1.,1.,1.,.6])
+ interp = SN.InterpolateSNd(points,z)
+
+ X = np.array([ [.1,.2,.1,.1],
+ [.1,.1,.2,.1],
+ [.1,.1,.1,.2] ])
+
+ output = interp(X)
+
+ print "output:\n",output
+
+if False:
+ P = [np.random.random_sample(3) for i in range(7)]
+ print "P:",P
+ tri = dw.dewall(P)
\ No newline at end of file
Modified: branches/interpolate/interpSNd/interpolateSNd.py
===================================================================
--- branches/interpolate/interpSNd/interpolateSNd.py 2008-08-19 22:26:16 UTC (rev 4655)
+++ branches/interpolate/interpSNd/interpolateSNd.py 2008-08-19 23:39:17 UTC (rev 4656)
@@ -41,8 +41,8 @@
assert P.ndim == 2, "P must be 2-dimensional"
d, n = P.shape
- assert len(fvals)==n, "fvals must have length n,\
- where n is number of points"
+ assert len(fvals)==n, \
+ "fvals must have length n, where n is number of points"
# remember dimensionality of space
self.ndim = d
Deleted: branches/interpolate/interpSNd/interpolateSNd.pyc
===================================================================
(Binary files differ)
Deleted: branches/interpolate/interpSNd/output.txt
===================================================================
--- branches/interpolate/interpSNd/output.txt 2008-08-19 22:26:16 UTC (rev 4655)
+++ branches/interpolate/interpSNd/output.txt 2008-08-19 23:39:17 UTC (rev 4656)
@@ -1,40 +0,0 @@
-
-face: [array([ 0.78256165, 1.25794206]), array([ 5.5867473 , 0.30561524])]
-points: []
-
-face: [array([ 5.5867473 , 0.30561524]), array([ 0.84738355, 0.4417885 ])]
-points: [array([ 0.49630621, 0.38945595])]
-delaunay distances:
-[(14.481574504701088, array([ 0.49630621, 0.38945595]))]
-
-face: [array([ 5.5867473 , 0.30561524]), array([ 0.49630621, 0.38945595])]
-points: []
-
-face: [array([ 0.78256165, 1.25794206]), array([ 0.84738355, 0.4417885 ])]
-points: [array([ 0.45697308, 1.26267539]), array([ 0.49630621, 0.38945595]), array([ 0.15363154, 0.794239 ])]
-delaunay distances:
-[(0.45650505613020009, array([ 0.45697308, 1.26267539])), (0.45830428089834485, array([ 0.49630621, 0.38945595])), (0.45808830610514073, array([ 0.15363154, 0.794239 ]))]
-
-face: [array([ 0.84738355, 0.4417885 ]), array([ 0.45697308, 1.26267539])]
-points: [array([ 0.49630621, 0.38945595]), array([ 0.15363154, 0.794239 ])]
-delaunay distances:
-[(0.45691821155399881, array([ 0.49630621, 0.38945595])), (0.45699647524242137, array([ 0.15363154, 0.794239 ]))]
-
-face: [array([ 0.45697308, 1.26267539]), array([ 0.49630621, 0.38945595])]
-points: [array([ 0.15363154, 0.794239 ])]
-delaunay distances:
-[(-0.45659643493482976, array([ 0.15363154, 0.794239 ]))]
-
-face: [array([ 0.45697308, 1.26267539]), array([ 0.15363154, 0.794239 ])]
-points: []
-
-face: [array([ 0.84738355, 0.4417885 ]), array([ 0.49630621, 0.38945595])]
-points: []
-
-face: [array([ 0.84738355, 0.4417885 ]), array([ 0.49630621, 0.38945595])]
-points: [array([ 0.15363154, 0.794239 ])]
-delaunay distances:
-[(0.45765190740816952, array([ 0.15363154, 0.794239 ]))]
-
-face: [array([ 0.84738355, 0.4417885 ]), array([ 0.15363154, 0.794239 ])]
-points: []
Modified: branches/interpolate/interpSNd/test.py
===================================================================
--- branches/interpolate/interpSNd/test.py 2008-08-19 22:26:16 UTC (rev 4655)
+++ branches/interpolate/interpSNd/test.py 2008-08-19 23:39:17 UTC (rev 4656)
@@ -8,7 +8,7 @@
class Test(unittest.TestCase):
- def compare_array(self, a, b):
+ def compare_arrays(self, a, b):
return np.allclose(a,b)
## test Delaunay triangulation itself
@@ -21,16 +21,48 @@
self.assert_( len(tri)==2 )
self.assert_( len(tri[0])==3 )
self.assert_( len(tri[1])==3 )
- print "square triangulation:\n", tri
def test_linear(self):
P = [array([0.,1.]), array([0.,0.]), array([0.,-1.])]
tri = dw.dewall(P)
- print "line triang:\n", tri
# testing general case using random data
+
+ def test_2d(self):
+ print "TESTING 2D"
+ P = [np.random.random_sample(2) for i in range(15)]
+ tri = dw.dewall(P)
+ # not checking if its correct, just if it runs.
+
+ 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):
+ 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):
+ 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):
x = array([0., 1, 0, 1, 0, 1, 0, 1])
y = array([0., 0, 1, 1, 0, 0, 1, 1])
@@ -40,14 +72,24 @@
interp = SNd.InterpolateSNd(points, fvals)
- newdata = np.random.random_sample((3,20))
+ newdata = np.random.random_sample((3,8))
interpvals = interp(newdata)
-
realvals = newdata[0,:]+newdata[1,:]-newdata[2,:]
- self.assert_(compare_array(np.ravel(interpvals), np.ravel(realvals)))
+ 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]
+ interp = SNd.InterpolateSNd(points, fvals)
+ newdata = np.random.random_sample((5,8))
+ interpvals = interp(newdata)
+ realvals = np.sum(newdata, axis=0)
+ self.assert_(self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)))
+
if __name__ == "__main__":
unittest.main()
\ No newline at end of file
More information about the Scipy-svn
mailing list