[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
 
 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