[Scipy-svn] r4656 - branches/interpolate/interpSNd

scipy-svn@scip... scipy-svn@scip...
Tue Aug 19 18:39:20 CDT 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

-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

```