[Scipy-svn] r4660 - in branches/interpolate: docs interpSNd

scipy-svn@scip... scipy-svn@scip...
Wed Aug 20 19:01:13 CDT 2008


Author: fcady
Date: 2008-08-20 19:01:08 -0500 (Wed, 20 Aug 2008)
New Revision: 4660

Modified:
   branches/interpolate/docs/tutorial.rst
   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:
more work, but there are still MAJOR bugs in dewall.py

Modified: branches/interpolate/docs/tutorial.rst
===================================================================
--- branches/interpolate/docs/tutorial.rst	2008-08-20 18:54:38 UTC (rev 4659)
+++ branches/interpolate/docs/tutorial.rst	2008-08-21 00:01:08 UTC (rev 4660)
@@ -737,8 +737,14 @@
 ND Scattered Interpolation
 ================================================
  
- Still in development.
+ InterpolateND has the significant disadvantage of requiring uniformly spaced
+ data, which is often not available in applications.  To remedy this problem, 
+ there is a callable class in development that does linear interpolation over the
+ convex hull of the known data points.
  
- Ideally the range of interpolation would be the convex hull of the known
- data points, and a Delaunay tesselation would be determined and stored
- at instantiation.  Then again, that would be very expensive.
\ No newline at end of file
+ The class is InterpolateSNd, for "scattered n-dimensional" data.  However, the
+ code is still extremely slow, and more importantly there are significant bugs remaining.
+ Depending on chance (the algorithm uses randomness) and the data used, the convex
+ hull can be improperly calculated, and there is even the possibility of infinite recursion.
+ That said, the algorithm works in other cases, and is expected to improve greatly in
+ the future.
\ No newline at end of file

Modified: branches/interpolate/interpSNd/dewall.py
===================================================================
--- branches/interpolate/interpSNd/dewall.py	2008-08-20 18:54:38 UTC (rev 4659)
+++ branches/interpolate/interpSNd/dewall.py	2008-08-21 00:01:08 UTC (rev 4660)
@@ -11,8 +11,8 @@
 # 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.
+# This code doesn't always work.  There are still significant
+# bugs to be fixed.
 
 import numpy as np
 from numpy.linalg import norm
@@ -28,6 +28,8 @@
                 AFL = [], # list of faces: (d-1)face list
                 ):
     
+    print "DEWALL: ",len(AFL)
+    
     # checking input
     assert isinstance(P, list)
     if len(P)>0:
@@ -42,6 +44,7 @@
     # base case
     if len(P) == 0: return []
     if len(P) <= len(P[0]): return [] # better to return [] whenever points are co-hyper-planar
+    if len(P) == len(P[0])+1: return P if linearly_independent([ point-P[0] for point in P[1:]]) else [] # simplex iff not on same hyperplane
     
     # lists of active faces
     # elem = ( [list of d points] , outvec )
@@ -53,7 +56,7 @@
     Sigma= []
     
     alpha = select_alpha(P, AFL)
-    print "\nalpha:\n", alpha
+    print "alpha:\n", alpha
                             
     # divide points into two sets separated by alpha
     P1, P2 = pointset_partition(P, alpha) # both lists of points
@@ -61,40 +64,59 @@
     # Simplex Wall Construction
     just_starting = False #source of problem?
     if len(AFL) == 0: # This is only executed once, at the start of the algorithm
-        just_starting = True
+        #just_starting = True
         first_simplex = make_first_simplex(P, alpha)
-        print "\nfirst simplex:\n",first_simplex
+        print "first_simplex:"
+        display(first_simplex)
         AFL = [ (face, get_out_vec(face,first_simplex))\
                 for face in faces(first_simplex)] # d+1 of them
         Sigma.append(first_simplex)
     for face, outvec in AFL:
         if is_intersected(face, alpha): # not counting as intersected
             AFL_alpha.append((face, \
-                            get_out_vec(face,first_simplex) if just_starting \
-                            else outvec))
-        if is_subset(face, P1):
+                            outvec))#get_out_vec(face,first_simplex) if just_starting \
+                            #else outvec))
+        elif is_subset(face, P1):
             AFL1.append((face,outvec))
-        if is_subset(face, P2):
+        elif is_subset(face, P2):
             AFL2.append((face,outvec))
     while len(AFL_alpha) != 0:
-        print "\nAFL_alpha start:",[face for face, vec in AFL_alpha]
+        print "\nLENGTH: ",len(AFL_alpha)
+        print "AFL_alpha"
+        for FACE, VEC in AFL_alpha: 
+            print "face"
+            display(FACE)
+        print "AFL1"
+        for FACE, VEC in AFL1: 
+            print "face"
+            display(FACE)
+        print "AFL2"
+        for FACE, VEC in AFL2: 
+            print "face"
+            display(FACE)
         
         face, outvec = AFL_alpha.pop()
         
+        print '\nprocessing face:'
+        display(face)
+        
         if outvec is not None:
             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)
+        #outward_points = filter( lambda p: (np.dot(p,outvec)>np.dot(face[0],outvec)),\
+            #                                P)
                 
-        t = make_simplex(face, outward_points) # make only over outer half space
-        print "\nnew simplex:\n",t
+        t = make_simplex(face, outward_points)
+        print "\nnew simplex:"
+        display(t)
         
         if t is not None:
             Sigma.append(t)
             # for f0 != f in faces(t) , ie for new outward faces
             for f0 in filter(lambda f: not face_in_list(f,[face]), faces(t)):
-                new_pair = (f0, get_out_vec(f0, t))
+                new_pair = (f0, get_out_vec(f0,t) )
                 if is_intersected(f0, alpha):
                     # continue building wall out in that direction
                     AFL_alpha = update(new_pair, AFL_alpha)
@@ -104,21 +126,34 @@
                 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
+    if len(AFL1)!=0:
+        print "RECURSING P1"
+        print "P1:"
+        display(P1)
+        print "AFL1:"
+        for face, vec in AFL1: display(face)
+        Sigma = Sigma + dewall(P1,AFL1)
     if len(AFL2)!=0:
+        print "RECURSING P2"
+        print "P2:"
+        display(P2)
+        print "AFL2:"
+        for face, vec in AFL2: display(face)
         Sigma = Sigma + dewall(P2,AFL2)
-    if len(AFL1)!=0:
-        Sigma = Sigma + dewall(P1,AFL1) # not doing this branch of recursion
     return Sigma
-
+def display(FACE):
+        if FACE==None: print 'NONE'
+        else:
+            for arr in FACE: print arr
 def select_alpha(P, AFL):
     # dividing plane.  This must divide points into 2 non-empty parts
-    # must intersect active face if there is one
-    # must not intersect any points
     
+    # creates random plane which divides into 2 non-empty parts
+    # and, if AFL not empty, passes through one of its faces
+    
     d = len(P[0])
     
     # plane through avg of cluster.  Guarantees separation
@@ -127,11 +162,14 @@
         mid = sum(AFL[0][0])/(d)
     else:
         mid = sum(P)/len(P)
+    direction = np.random.random_sample((d))
+    direction = direction/norm(direction)
     if norm(mid)==0:
-        direction = np.random.random_sample((d))
-        alpha = ( direction/norm(direction), 0)
+        alpha = ( direction , 0)
     else: 
-        alpha =(mid/norm(mid), norm(mid))
+        alpha =(direction, np.dot(direction,mid))
+    if alpha[1] < 0:
+        alpha = (-1*alpha[0], -1*alpha[1])
         
     return alpha
 
@@ -170,22 +208,22 @@
     # both are lists of arrays
     return np.alltrue([ point_in_list(s1, S2) for s1 in S1])
 
-def update (face_pair, face_pair_list): # this func has been problematic
+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=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)
+    face = face_pair[0] #list of arrays
+    Face_list = [Face_pair[0] for Face_pair in Face_pair_list] #list of lists(or arrays)
+    if face_in_list(face, Face_list):
+        #print "face in list"
+        def Face_pair_not_equal_face(Pair):
+            Face = Pair[0]
+            vertex_of_face_in_Face = [point_in_list(vert,Face) for vert in face]
+            return not np.all(vertex_of_face_in_Face)
+        #Face_pair_not_equal_face = lambda Pair :  not np.alltrue([ point_in_list(p, Pair[0]) for p in face ])
+        Face_pair_list = filter(Face_pair_not_equal_face, Face_pair_list)
     else:
-        print "face not in list"
-        face_pair_list.append(face_pair)
-    return face_pair_list
+        #print "not in list"
+        Face_pair_list.append(face_pair)
+    return Face_pair_list
         
 def pointset_partition(P, alpha): #WORKS
     P1 = [p for p in P if np.dot(p,alpha[0])<  alpha[1]]
@@ -211,7 +249,7 @@
 
 def make_first_simplex(P, alpha): #WORKS
     # alpha = unit_normal_vec, distance
-    # assume no points on plane
+    # ENSURE!! no points on same hyperplane
 
     d = len(P[0])
     unit_normal, distance = alpha
@@ -222,21 +260,33 @@
     first_point = sorted(points_and_dist_from_alpha,
                                 cmp = compare_first_elem \
                                 )[0][1] # closest to alpha
-                                                    
+    print "first points:\n"
+    display(first_point)
     possible_second_pts = [ (circumcircle([first_point, point])[1], point) \
-                                        for point in P if (point != first_point).any() and \
+                                        for point in P if \
+                                        (point != first_point).any() and \
                                         (np.dot(unit_normal, first_point)-distance)*(np.dot(unit_normal, point)-distance)<=0 \
                                     ]
     second_point = sorted(possible_second_pts, \
                                     cmp = compare_first_elem \
                                     )[0][1]
+    print "second point:\n"
+    display(second_point)
     simplex = [first_point, second_point]
     for i in range(d-1):
         radii_of_circumcircles = [(circumcircle(  copy_list(simplex)+[point.copy()]  )[1], point) \
-                                        for point in P if not point_in_list(point, simplex) ]
-        new_point = sorted(radii_of_circumcircles, \
+                                        for point in P if \
+                                            not point_in_list(point, simplex) and \
+                                            linearly_independent([(point-first_point)]+[ (p-first_point) for p in simplex[1:] ]) \
+                                        ]
+        #new_point = sorted(radii_of_circumcircles, \
+          #                          cmp = compare_first_elem \
+            #                        )[0][1]
+        new_point_pair = sorted(radii_of_circumcircles, \
                                     cmp = compare_first_elem \
-                                    )[0][1]
+                                    )[0]
+        new_point = new_point_pair[1]
+        print "new point:\n", new_point
         simplex.append(new_point)
 
     return simplex
@@ -249,14 +299,14 @@
     
 
         
-def make_simplex(f,P): #WORKS
+def make_simplex(f,outward_points): #WORKS
     # must be only in the outer halfspace
     
     # returns the simlex
     
     # clean up by making sure only points not in hyperplane(f) are in P
         
-    valid_points = [p for p in P if not point_in_list(p,f)]
+    valid_points = outward_points#[p for p in P if not point_in_list(p,f)]
     if len(valid_points) == 0: 
         return None
     delaunay_distances = [(delaunay_dist(f,p) ,p) for p in valid_points]
@@ -391,9 +441,24 @@
     d=len(P[0])
     if len(P)>d: return False
     matrix = np.array(P).reshape((d,len(P)))
-    return np.rank(matrix)==len(P)
+    return rank(matrix)==len(P)
+
+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]) > 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))
     
+    
+    
+    
     
\ No newline at end of file

Modified: branches/interpolate/interpSNd/dewall.pyc
===================================================================
(Binary files differ)

Modified: branches/interpolate/interpSNd/interp.py
===================================================================
--- branches/interpolate/interpSNd/interp.py	2008-08-20 18:54:38 UTC (rev 4659)
+++ branches/interpolate/interpSNd/interp.py	2008-08-21 00:01:08 UTC (rev 4660)
@@ -4,7 +4,7 @@
 
 reload(SN)
 
-if True:
+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)
@@ -20,12 +20,14 @@
     
     
 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])
+    points = np.array([  [0., 0, 0, 1., 1., 1., 0., 1., .2],
+                                [0., 1., 0, 0, 1., 0., 1., 1., .2],
+                                [0., 0, 1., 0, 0., 1., 1., 1., .2] ])
+    z = np.sum(points,axis=0).reshape(points.shape[1])
     interp = SN.InterpolateSNd(points,z)
     
+    print "*"*100+'\nMADE IT'
+    
     X = np.array([ [.1,.2,.1,.1],
                         [.1,.1,.2,.1],
                         [.1,.1,.1,.2] ])

Modified: branches/interpolate/interpSNd/interpolateSNd.py
===================================================================
--- branches/interpolate/interpSNd/interpolateSNd.py	2008-08-20 18:54:38 UTC (rev 4659)
+++ branches/interpolate/interpSNd/interpolateSNd.py	2008-08-21 00:01:08 UTC (rev 4660)
@@ -165,7 +165,7 @@
                                                         )
         
         if rank(matrix_to_solve) < d:
-            print "matrix rank: "+str(np.rank(matrix_to_solve))
+            print "matrix rank: "+str(rank(matrix_to_solve))
             print "needed dim:  "+str(d)
             weights_vec = None
         else:

Modified: branches/interpolate/interpSNd/script.py
===================================================================
--- branches/interpolate/interpSNd/script.py	2008-08-20 18:54:38 UTC (rev 4659)
+++ branches/interpolate/interpSNd/script.py	2008-08-21 00:01:08 UTC (rev 4660)
@@ -31,12 +31,15 @@
         for j in range(10) ]
 
 triangul = dw.dewall(P)
+print "triangulation:\n",triangul
 
+
 # plotting the known data points
 for p in P: pyplot.scatter([p[0]],[p[1]])
 
 # plotting the triangulation
 for tri in triangul:
+    print "triangle:\n", tri
     for seg in segments(tri):
         pyplot.plot(seg[0],seg[1])
         
@@ -46,4 +49,3 @@
         
 pyplot.show()
 
-print "triangulation:\n",triangul
\ No newline at end of file

Modified: branches/interpolate/interpSNd/test.py
===================================================================
--- branches/interpolate/interpSNd/test.py	2008-08-20 18:54:38 UTC (rev 4659)
+++ branches/interpolate/interpSNd/test.py	2008-08-21 00:01:08 UTC (rev 4660)
@@ -9,7 +9,7 @@
 class Test(unittest.TestCase):
     
     def compare_arrays(self, a, b):
-        return np.allclose(a,b,rtol=1e-3) | (np.isnan(a)&np.isnan(b))
+        return np.allclose(a,b,rtol=1e-3) or (np.isnan(a)&np.isnan(b)).all()
     
     ## test Delaunay triangulation itself
     
@@ -27,7 +27,6 @@
         tri = dw.dewall(P)
         
     # testing general case using random data
-    
     def _test_2d(self):
         ndim = 2
         nold = 15
@@ -60,7 +59,6 @@
         # not checking if its correct, just if it runs.
         
     ## test interpolation, and thus also triangulation by extension
-    
     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])
@@ -110,6 +108,40 @@
         #self.assert_(self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)))
         assert self.compare_arrays(np.ravel(interpvals), np.ravel(realvals)), "wrong data"
 
+    def test_linear_3d(self):
+        ndim = 3 # num dimensions
+        nold = 1 # 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)
+        
+        #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)
+        
+        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))



More information about the Scipy-svn mailing list