# [Scipy-svn] r4659 - branches/interpolate/interpSNd

scipy-svn@scip... scipy-svn@scip...
Wed Aug 20 13:54:39 CDT 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 @@

-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

```