Changeset 5907

Timestamp:

Nov 6, 2008, 4:15:04 PM (16 years ago)

Author:

rwilson

Message:

NumPy? conversion.

Location:

anuga_core/source_numpy_conversion/anuga/alpha_shape

Files:

• alpha_shape.py (modified) (6 diffs)
• test_alpha_shape.py (modified) (11 diffs)

anuga_core/source_numpy_conversion/anuga/alpha_shape/alpha_shape.py

@@ -20 +20 @@
 
 import exceptions
-from Numeric import array, Float, divide_safe, sqrt, product
+import numpy
 import random
 
@@ -87 +87 @@
         
         #Convert input to Numeric arrays
-        self.points = array(points).astype(Float)
+        self.points = numpy.array(points).astype(numpy.float)
 
     
@@ -240 +240 @@
         ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]
 
-        x1 = array([ x[j] for j in ind1 ])
-        y1 = array([ y[j] for j in ind1 ])
-        x2 = array([ x[j] for j in ind2 ])
-        y2 = array([ y[j] for j in ind2 ])
-        x3 = array([ x[j] for j in ind3 ])
-        y3 = array([ y[j] for j in ind3 ])
+        x1 = numpy.array([ x[j] for j in ind1 ])
+        y1 = numpy.array([ y[j] for j in ind1 ])
+        x2 = numpy.array([ x[j] for j in ind2 ])
+        y2 = numpy.array([ y[j] for j in ind2 ])
+        x3 = numpy.array([ x[j] for j in ind3 ])
+        y3 = numpy.array([ y[j] for j in ind3 ])
 
         x21 = x2-x1
@@ -286 +286 @@
             zeroind = [k for k in range(len(denom)) if \
                        (denom[k]< EPSILON and  denom[k] > -EPSILON)]
-        try:
-            dx = divide_safe(y31*dist21 - y21*dist31,denom)
-            dy = divide_safe(x21*dist31 - x31*dist21,denom)
-        except ZeroDivisionError:
+##NumPy        try:
+##NumPy            dx = divide_safe(y31*dist21 - y21*dist31,denom)
+##NumPy            dy = divide_safe(x21*dist31 - x31*dist21,denom)
+##NumPy        except ZeroDivisionError:
+##NumPy            raise  AlphaError
+# In NumPy there is no divide_safe(), so above code becomes:
+        dx = numpy.divide(y31*dist21 - y21*dist31,denom)
+        dy = numpy.divide(x21*dist31 - x31*dist21,denom)
+        if not numpy.alltrue(numpy.isfinite(dx)) or \
+           not numpy.alltrue(numpy.isfinite(dy)):
             raise  AlphaError
             
-        self.triradius = 0.5*sqrt(dx*dx + dy*dy)
+        self.triradius = 0.5*numpy.sqrt(dx*dx + dy*dy)
         #print "triangle radii", self.triradius
 
@@ -313 +319 @@
             tri = self.deltri[t]
             trinbr = self.deltrinbr[t]
-            dx = array([self.points[tri[(i+1)%3],0] -
+            dx = numpy.array([self.points[tri[(i+1)%3],0] -
                         self.points[tri[(i+2)%3],0] for i in [0,1,2]])
-            dy = array([self.points[tri[(i+1)%3],1] -
+            dy = numpy.array([self.points[tri[(i+1)%3],1] -
                         self.points[tri[(i+2)%3],1] for i in [0,1,2]])
-            elen = sqrt(dx*dx+dy*dy)
+            elen = numpy.sqrt(dx*dx+dy*dy)
             # really only need sign - not angle value:
-            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-
+            anglesign = numpy.array([(-dx[(i+1)%3]*dx[(i+2)%3]-
                                 dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
             
@@ -594 +600 @@
             tri = self.deltri[tind[0]]
             
-            dx = array([self.points[tri[(i+1)%3],0] - \
+            dx = numpy.array([self.points[tri[(i+1)%3],0] - \
                         self.points[tri[(i+2)%3],0] for i in [0,1,2]])
-            dy = array([self.points[tri[(i+1)%3],1] - \
+            dy = numpy.array([self.points[tri[(i+1)%3],1] - \
                         self.points[tri[(i+2)%3],1] for i in [0,1,2]])
-            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
+            anglesign = numpy.array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
                                 dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
             # record any triangle that has an acute angle spanned by

anuga_core/source_numpy_conversion/anuga/alpha_shape/test_alpha_shape.py

@@ -4 +4 @@
 import sys
 import unittest
-from Numeric import allclose
+import numpy
 
 try:
@@ -32 +32 @@
         result = alpha.get_delaunay()
         answer = [(0, 1, 5), (5, 1, 4), (4, 2, 3), (2, 4, 1)]
-        assert allclose(answer, result) 
+        assert numpy.allclose(answer, result) 
 
     def test_3_points_on_line(self):
@@ -63 +63 @@
         #print "result",result
         answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
-        assert allclose(answer, result) 
+        assert numpy.allclose(answer, result) 
 
 
@@ -80 +80 @@
         #print "result",result
         answer = [(0, 3), (3, 6), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
-        assert allclose(answer, result) 
+        assert numpy.allclose(answer, result) 
 
 
@@ -97 +97 @@
         #print "result",result
         answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
-        assert allclose(answer, result) 
+        assert numpy.allclose(answer, result) 
 
 
@@ -114 +114 @@
         #print "result",result
         answer = [(0, 3), (3, 6), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
-        assert allclose(answer, result) 
+        assert numpy.allclose(answer, result) 
    
 
@@ -132 +132 @@
         #print "result",result
         answer = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 0)]
-        assert allclose(answer, result) 
+        assert numpy.allclose(answer, result) 
 
 
@@ -148 +148 @@
         result = alpha.get_boundary()
         answer = [(1, 0), (0, 5), (3, 2), (4, 3), (2, 6), (6, 1), (5, 4)]
-        assert allclose(answer, result)
+        assert numpy.allclose(answer, result)
 
 
@@ -222 +222 @@
         #print "result",result
         answer = [(1, 0), (4, 3), (0, 4), (3, 1)]
-        assert allclose(answer, result)
+        assert numpy.allclose(answer, result)
   
     def test_sharp_indents(self):
@@ -239 +239 @@
         #print "result",result
         answer = [(3, 4), (2, 3), (0, 1), (1, 2), (4, 0)]
-        assert allclose(answer, result)
+        assert numpy.allclose(answer, result)
 
         
@@ -392 +392 @@
                   (61, 62), (62, 63), (59, 60), (58, 59), (60, 61), \
                   (46, 45)]
-        assert allclose(answer, result)    
+        assert numpy.allclose(answer, result)    
 #-------------------------------------------------------------
 if __name__ == "__main__":