Changeset 5907


Ignore:
Timestamp:
Nov 6, 2008, 4:15:04 PM (15 years ago)
Author:
rwilson
Message:

NumPy? conversion.

Location:
anuga_core/source_numpy_conversion/anuga/alpha_shape
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source_numpy_conversion/anuga/alpha_shape/alpha_shape.py

    r4955 r5907  
    2020
    2121import exceptions
    22 from Numeric import array, Float, divide_safe, sqrt, product
     22import numpy
    2323import random
    2424
     
    8787       
    8888        #Convert input to Numeric arrays
    89         self.points = array(points).astype(Float)
     89        self.points = numpy.array(points).astype(numpy.float)
    9090
    9191   
     
    240240        ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]
    241241
    242         x1 = array([ x[j] for j in ind1 ])
    243         y1 = array([ y[j] for j in ind1 ])
    244         x2 = array([ x[j] for j in ind2 ])
    245         y2 = array([ y[j] for j in ind2 ])
    246         x3 = array([ x[j] for j in ind3 ])
    247         y3 = array([ y[j] for j in ind3 ])
     242        x1 = numpy.array([ x[j] for j in ind1 ])
     243        y1 = numpy.array([ y[j] for j in ind1 ])
     244        x2 = numpy.array([ x[j] for j in ind2 ])
     245        y2 = numpy.array([ y[j] for j in ind2 ])
     246        x3 = numpy.array([ x[j] for j in ind3 ])
     247        y3 = numpy.array([ y[j] for j in ind3 ])
    248248
    249249        x21 = x2-x1
     
    286286            zeroind = [k for k in range(len(denom)) if \
    287287                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
    288         try:
    289             dx = divide_safe(y31*dist21 - y21*dist31,denom)
    290             dy = divide_safe(x21*dist31 - x31*dist21,denom)
    291         except ZeroDivisionError:
     288##NumPy        try:
     289##NumPy            dx = divide_safe(y31*dist21 - y21*dist31,denom)
     290##NumPy            dy = divide_safe(x21*dist31 - x31*dist21,denom)
     291##NumPy        except ZeroDivisionError:
     292##NumPy            raise  AlphaError
     293# In NumPy there is no divide_safe(), so above code becomes:
     294        dx = numpy.divide(y31*dist21 - y21*dist31,denom)
     295        dy = numpy.divide(x21*dist31 - x31*dist21,denom)
     296        if not numpy.alltrue(numpy.isfinite(dx)) or \
     297           not numpy.alltrue(numpy.isfinite(dy)):
    292298            raise  AlphaError
    293299           
    294         self.triradius = 0.5*sqrt(dx*dx + dy*dy)
     300        self.triradius = 0.5*numpy.sqrt(dx*dx + dy*dy)
    295301        #print "triangle radii", self.triradius
    296302
     
    313319            tri = self.deltri[t]
    314320            trinbr = self.deltrinbr[t]
    315             dx = array([self.points[tri[(i+1)%3],0] -
     321            dx = numpy.array([self.points[tri[(i+1)%3],0] -
    316322                        self.points[tri[(i+2)%3],0] for i in [0,1,2]])
    317             dy = array([self.points[tri[(i+1)%3],1] -
     323            dy = numpy.array([self.points[tri[(i+1)%3],1] -
    318324                        self.points[tri[(i+2)%3],1] for i in [0,1,2]])
    319             elen = sqrt(dx*dx+dy*dy)
     325            elen = numpy.sqrt(dx*dx+dy*dy)
    320326            # really only need sign - not angle value:
    321             anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-
     327            anglesign = numpy.array([(-dx[(i+1)%3]*dx[(i+2)%3]-
    322328                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
    323329           
     
    594600            tri = self.deltri[tind[0]]
    595601           
    596             dx = array([self.points[tri[(i+1)%3],0] - \
     602            dx = numpy.array([self.points[tri[(i+1)%3],0] - \
    597603                        self.points[tri[(i+2)%3],0] for i in [0,1,2]])
    598             dy = array([self.points[tri[(i+1)%3],1] - \
     604            dy = numpy.array([self.points[tri[(i+1)%3],1] - \
    599605                        self.points[tri[(i+2)%3],1] for i in [0,1,2]])
    600             anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
     606            anglesign = numpy.array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
    601607                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
    602608            # record any triangle that has an acute angle spanned by
  • anuga_core/source_numpy_conversion/anuga/alpha_shape/test_alpha_shape.py

    r4663 r5907  
    44import sys
    55import unittest
    6 from Numeric import allclose
     6import numpy
    77
    88try:
     
    3232        result = alpha.get_delaunay()
    3333        answer = [(0, 1, 5), (5, 1, 4), (4, 2, 3), (2, 4, 1)]
    34         assert allclose(answer, result)
     34        assert numpy.allclose(answer, result)
    3535
    3636    def test_3_points_on_line(self):
     
    6363        #print "result",result
    6464        answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
    65         assert allclose(answer, result)
     65        assert numpy.allclose(answer, result)
    6666
    6767
     
    8080        #print "result",result
    8181        answer = [(0, 3), (3, 6), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
    82         assert allclose(answer, result)
     82        assert numpy.allclose(answer, result)
    8383
    8484
     
    9797        #print "result",result
    9898        answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
    99         assert allclose(answer, result)
     99        assert numpy.allclose(answer, result)
    100100
    101101
     
    114114        #print "result",result
    115115        answer = [(0, 3), (3, 6), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
    116         assert allclose(answer, result)
     116        assert numpy.allclose(answer, result)
    117117   
    118118
     
    132132        #print "result",result
    133133        answer = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 0)]
    134         assert allclose(answer, result)
     134        assert numpy.allclose(answer, result)
    135135
    136136
     
    148148        result = alpha.get_boundary()
    149149        answer = [(1, 0), (0, 5), (3, 2), (4, 3), (2, 6), (6, 1), (5, 4)]
    150         assert allclose(answer, result)
     150        assert numpy.allclose(answer, result)
    151151
    152152
     
    222222        #print "result",result
    223223        answer = [(1, 0), (4, 3), (0, 4), (3, 1)]
    224         assert allclose(answer, result)
     224        assert numpy.allclose(answer, result)
    225225 
    226226    def test_sharp_indents(self):
     
    239239        #print "result",result
    240240        answer = [(3, 4), (2, 3), (0, 1), (1, 2), (4, 0)]
    241         assert allclose(answer, result)
     241        assert numpy.allclose(answer, result)
    242242
    243243       
     
    392392                  (61, 62), (62, 63), (59, 60), (58, 59), (60, 61), \
    393393                  (46, 45)]
    394         assert allclose(answer, result)   
     394        assert numpy.allclose(answer, result)   
    395395#-------------------------------------------------------------
    396396if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.