Changeset 5891

Timestamp:

Nov 5, 2008, 4:19:05 PM (16 years ago)

Author:

rwilson

Message:

Converted to full numpy.

Location:

anuga_core/source/anuga/utilities

Files:

• cg_solve.py (modified) (5 diffs)
• compile.py (modified) (2 diffs)
• interp.py (modified) (5 diffs)
• numerical_tools.py (modified) (9 diffs)
• polygon.py (modified) (25 diffs)
• polygon_ext.c (modified) (4 diffs)
• quad.py (modified) (3 diffs)
• sparse.py (modified) (15 diffs)
• sparse_ext.c (modified) (1 diff)
• test_cg_solve.py (modified) (12 diffs)
• test_data_audit.py (modified) (1 diff)
• test_numerical_tools.py (modified) (8 diffs)
• test_polygon.py (modified) (41 diffs)
• test_quad.py (modified) (3 diffs)
• test_sparse.py (modified) (11 diffs)
• test_system_tools.py (modified) (4 diffs)
• test_xml_tools.py (modified) (1 diff)
• util_ext.c (modified) (1 diff)
• util_ext.h (modified) (1 diff)

anuga_core/source/anuga/utilities/cg_solve.py

@@ +1 @@
+
 import exceptions
 class VectorShapeError(exceptions.Exception): pass
 class ConvergenceError(exceptions.Exception): pass
 
-from Numeric import dot, array, Float, zeros
+import numpy
    
 import logging, logging.config
@@ -24 +25 @@
     
     if x0 is None:
-        x0 = zeros(b.shape, typecode=Float)
+        x0 = numpy.zeros(b.shape, dtype=numpy.float)
     else:
-        x0 = array(x0, typecode=Float)
+        x0 = numpy.array(x0, dtype=numpy.float)
 
-    b  = array(b, typecode=Float)
+    b  = numpy.array(b, dtype=numpy.float)
     if len(b.shape) != 1 :
        
@@ -57 +58 @@
 
 
-   b  = array(b, typecode=Float)
+   b  = numpy.array(b, dtype=numpy.float)
    if len(b.shape) != 1 :
       raise VectorShapeError, 'input vector should consist of only one column'
 
    if x0 is None:
-      x0 = zeros(b.shape, typecode=Float)
+      x0 = numpy.zeros(b.shape, dtype=numpy.float)
    else:
-      x0 = array(x0, typecode=Float)
+      x0 = numpy.array(x0, dtype=numpy.float)
 
 
@@ -75 +76 @@
    r = b - A*x
    d = r
-   rTr = dot(r,r)
+   rTr = numpy.dot(r,r)
    rTr0 = rTr
 
    while (i<imax and rTr>tol**2*rTr0):
        q = A*d
-       alpha = rTr/dot(d,q)
+       alpha = rTr/numpy.dot(d,q)
        x = x + alpha*d
        if i%50 :
@@ -87 +88 @@
            r = r - alpha*q
        rTrOld = rTr
-       rTr = dot(r,r)
+       rTr = numpy.dot(r,r)
        bt = rTr/rTrOld
 

anuga_core/source/anuga/utilities/compile.py

@@ -10 +10 @@
    Ole Nielsen, Duncan Gray Oct 2001      
 """     
+
+#NumPy ------------------------------------
+# Something like these lines recommended in "Converting from NUMARRAY to NUMPY"
+import numpy
+numpyincludedirs = numpy.get_include()
+I_dirs = '-I"%s" ' % numpyincludedirs
+#NumPy ------------------------------------
 
 # FIXME (Ole): Although this script says it works with a range of compilers,
@@ -255 +262 @@
     else:
       if FN == "triangle.c" or FN == "mesh_engine_c_layer.c":
-        s = '%s -c %s -I"%s" -I"%s" -o "%s.o" -O3 -DTRILIBRARY=1 -DNO_TIMER=1'\
-            %(compiler, FN, python_include, utilities_include_dir, root)
+#NumPy        s = '%s -c %s -I"%s" -I"%s" -o "%s.o" -O3 -DTRILIBRARY=1 -DNO_TIMER=1'\
+        s = '%s -c %s %s -I"%s" -I"%s" -o "%s.o" -O3 -DTRILIBRARY=1 -DNO_TIMER=1'\
+            %(compiler, FN, I_dirs, python_include, utilities_include_dir, root)
       else:
-        s = '%s -c %s -I"%s" -I"%s" -o "%s.o" -Wall -O3'\
-            %(compiler, FN, python_include, utilities_include_dir, root)
+#NumPy        s = '%s -c %s -I"%s" -I"%s" -o "%s.o" -Wall -O3'\
+        s = '%s -c %s %s -I"%s" -I"%s" -o "%s.o" -Wall -O3'\
+            %(compiler, FN, I_dirs, python_include, utilities_include_dir, root)
 
     if os.name == 'posix' and os.uname()[4] == 'x86_64':

anuga_core/source/anuga/utilities/interp.py

@@ -111 +111 @@
 
     >>> from interp import interp
-    >>> import Numeric as N
+    >>> import numpy.oldnumeric as N
     >>> x = N.array([1., 2., 3., 4., 5.])
     >>> y = N.array([3., 6., 2.,-5.,-3.])
@@ -138 +138 @@
     """
     import arrayfns
-    import MA
-    import Numeric as N
+    import numpy.ma as MA
+    import numpy as N
     from where_close import where_close
 
@@ -216 +216 @@
 
     >>> from interp import interp
-    >>> import Numeric as N
+    >>> import numpy.oldnumeric as N
     >>> x = N.array([1.,    2., 3.,  4.,  5.,  6.])
     >>> y = N.array([3., 1e+20, 2., -5., -3., -4.])
@@ -284 +284 @@
     >>> ['%.7g' % yint[i] for i in range(len(yint))]
     ['3.4', '2.3']
-    >>> yint.typecode()
+    >>> yint.dtype.char
     'd'
     >>> x = N.arange(6)
@@ -292 +292 @@
     >>> ['%.7g' % yint[i] for i in range(len(yint))]
     ['3', '2']
-    >>> yint.typecode()
+    >>> yint.dtype.char
     'd'
     """}

anuga_core/source/anuga/utilities/numerical_tools.py

@@ -7 +7 @@
 from warnings import warn
 
-#Establish which Numeric package to use
-#(this should move to somewhere central)
-#try:
-#    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt,
-# searchsorted, sort, concatenate, Float, arange    
-#except:
-#    #print 'Could not find scipy - using Numeric'
-#    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt,
-#searchsorted, sort, concatenate, Float, arange
-
-from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt,\
-     searchsorted, sort, concatenate, Float, arange    
+##from numpy import ndarray, array, sum, inner, ravel, sqrt, searchsorted, sort, concatenate, float, arange
+import numpy
 
 # Getting an infinite number to use when using Numeric
 #INF = (array([1])/0.)[0]
 
-NAN = (array([1])/0.)[0]
+NAN = (numpy.array([1])/0.)[0]
 # Note, INF is used instead of NAN (Not a number), since Numeric has no NAN
 # if we use a package that has NAN, this should be updated to use NAN.
@@ -89 +79 @@
         v2 = [1.0, 0.0] # Unit vector along the x-axis
         
-    v1 = ensure_numeric(v1, Float)
-    v2 = ensure_numeric(v2, Float)    
+    v1 = ensure_numeric(v1, numpy.float)
+    v2 = ensure_numeric(v2, numpy.float)    
     
     # Normalise
-    v1 = v1/sqrt(sum(v1**2))
-    v2 = v2/sqrt(sum(v2**2))
+    v1 = v1/numpy.sqrt(numpy.sum(v1**2))
+    v2 = v2/numpy.sqrt(numpy.sum(v2**2))
 
     # Compute angle
-    p = innerproduct(v1, v2)
-    c = innerproduct(v1, normal_vector(v2)) # Projection onto normal
+    p = numpy.inner(v1, v2)
+    c = numpy.inner(v1, normal_vector(v2)) # Projection onto normal
                                             # (negative cross product)
         
@@ -140 +130 @@
     """
     
-    return array([-v[1], v[0]], Float)
+    return numpy.array([-v[1], v[0]], numpy.float)
 
     
@@ -150 +140 @@
     """Mean value of a vector
     """
-    return(float(sum(x))/len(x))
+    return(float(numpy.sum(x))/len(x))
 
 
@@ -171 +161 @@
     cy = y - mean(y)  
 
-    p = innerproduct(cx,cy)/N
+    p = numpy.inner(cx,cy)/N
     return(p)
 
@@ -220 +210 @@
     """
   
-    y = ravel(x)
-    p = sqrt(innerproduct(y,y))
+    y = numpy.ravel(x)
+    p = numpy.sqrt(numpy.inner(y,y))
     return p
     
@@ -262 +252 @@
 
     if typecode is None:
-        if type(A) == ArrayType:
+##NumPy        if isinstance(A, ArrayType):
+        if type(A) == numpy.ndarray:
             return A
         else:
-            return array(A)
+            return numpy.array(A)
     else:
-        if type(A) == ArrayType:
-            if A.typecode == typecode:
-                return array(A)  #FIXME: Shouldn't this just return A?
+##NumPy        if isinstance(A, ArrayType):
+        if type(A) == numpy.ndarray:
+##NumPy            if A.typecode == typecode:
+            if A.dtype == typecode:
+                return numpy.array(A)  #FIXME: Shouldn't this just return A?
             else:
-                return array(A,typecode)
+                return numpy.array(A, typecode)
         else:
-            return array(A,typecode)
-
+            import types                            ##
+            from numpy import str                   ##
+            if isinstance(A, types.StringType):     ##
+                return numpy.array(A, dtype=int)          ##
+            return numpy.array(A, typecode)
 
 
@@ -285 +281 @@
     """
 
-    n = searchsorted(sort(a), bins)
-    n = concatenate( [n, [len(a)]] )
+    n = numpy.searchsorted(numpy.sort(a), bins)
+    n = numpy.concatenate( [n, [len(a)]] )
 
     hist = n[1:]-n[:-1]
 
     if relative is True:
-        hist = hist/float(sum(hist))
+        hist = hist/float(numpy.sum(hist))
         
     return hist 
@@ -301 +297 @@
     """
 
-    mx = max(data)
-    mn = min(data)
+    mx = numpy.max(data)
+    mn = numpy.min(data)
 
     if mx == mn:
-        bins = array([mn])
+        bins = numpy.array([mn])
     else:
         if number_of_bins is None:
             number_of_bins = 10
             
-        bins = arange(mn, mx, (mx-mn)/number_of_bins)
+        bins = numpy.arange(mn, mx, (mx-mn)/number_of_bins)
 
     return bins

anuga_core/source/anuga/utilities/polygon.py

@@ -10 +10 @@
 #    #print 'Could not find scipy - using Numeric'
 
-from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose
+##from numpy import float, int, zeros, ones, array, concatenate, reshape, dot, allclose, newaxis, ascontiguousarray
+import numpy
 
 
@@ -76 +77 @@
     # FIXME (Ole): Write this in C
 
-    line0 = ensure_numeric(line0, Float)
-    line1 = ensure_numeric(line1, Float)    
+    line0 = ensure_numeric(line0, numpy.float)
+    line1 = ensure_numeric(line1, numpy.float)    
 
     x0 = line0[0,0]; y0 = line0[0,1]
@@ -89 +90 @@
     u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
         
-    if allclose(denom, 0.0):
+    if numpy.allclose(denom, 0.0):
         # Lines are parallel - check if they coincide on a shared a segment
 
-        if allclose( [u0, u1], 0.0 ):
+        if numpy.allclose( [u0, u1], 0.0 ):
             # We now know that the lines if continued coincide
             # The remaining check will establish if the finite lines share a segment
@@ -120 +121 @@
             if line0_starts_on_line1 and line0_ends_on_line1:
                 # Shared segment is line0 fully included in line1
-                segment = array([[x0, y0], [x1, y1]])                
+                segment = numpy.array([[x0, y0], [x1, y1]])                
 
             if line1_starts_on_line0 and line1_ends_on_line0:
                 # Shared segment is line1 fully included in line0
-                segment = array([[x2, y2], [x3, y3]])
+                segment = numpy.array([[x2, y2], [x3, y3]])
             
 
@@ -130 +131 @@
             if line0_starts_on_line1 and line1_ends_on_line0:
                 # Shared segment from line0 start to line 1 end
-                segment = array([[x0, y0], [x3, y3]])
+                segment = numpy.array([[x0, y0], [x3, y3]])
 
             if line1_starts_on_line0 and line0_ends_on_line1:
                 # Shared segment from line1 start to line 0 end
-                segment = array([[x2, y2], [x1, y1]])                                
+                segment = numpy.array([[x2, y2], [x1, y1]])                                
 
 
@@ -140 +141 @@
             if line0_starts_on_line1 and line1_starts_on_line0:
                 # Shared segment from line0 start to line 1 end
-                segment = array([[x0, y0], [x2, y2]])
+                segment = numpy.array([[x0, y0], [x2, y2]])
 
             if line0_ends_on_line1 and line1_ends_on_line0:
                 # Shared segment from line0 start to line 1 end
-                segment = array([[x3, y3], [x1, y1]])                
+                segment = numpy.array([[x3, y3], [x1, y1]])                
 
                 
@@ -161 +162 @@
 
         # Sanity check - can be removed to speed up if needed
-        assert allclose(x, x2 + u1*(x3-x2))
-        assert allclose(y, y2 + u1*(y3-y2))        
+        assert numpy.allclose(x, x2 + u1*(x3-x2))
+        assert numpy.allclose(y, y2 + u1*(y3-y2))        
 
         # Check if point found lies within given line segments
@@ -168 +169 @@
             # We have intersection
 
-            return 1, array([x, y])
+            return 1, numpy.array([x, y])
         else:
             # No intersection
@@ -203 +204 @@
 
 
-    line0 = ensure_numeric(line0, Float)
-    line1 = ensure_numeric(line1, Float)    
+    line0 = ensure_numeric(line0, numpy.float)
+    line1 = ensure_numeric(line1, numpy.float)    
 
     status, value = _intersection(line0[0,0], line0[0,1],
@@ -256 +257 @@
         # converted to a numeric array.
         msg = 'Points could not be converted to Numeric array' 
-        raise msg
+        raise TypeError, msg
 
     try:
@@ -265 +266 @@
         # If this fails it is going to be because the points can't be
         # converted to a numeric array.
-        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
-        raise msg
+        msg = 'Polygon %s could not be converted to Numeric array' % (str(polygon))
+        raise TypeError, msg
 
     if len(points.shape) == 1:
         # Only one point was passed in. Convert to array of points
-        points = reshape(points, (1,2))
+        points = numpy.reshape(points, (1,2))
 
     indices, count = separate_points_by_polygon(points, polygon,
@@ -312 +313 @@
     #if verbose: print 'Checking input to outside_polygon'
     try:
-        points = ensure_numeric(points, Float)
+        points = ensure_numeric(points, numpy.float)
     except NameError, e:
         raise NameError, e
@@ -320 +321 @@
 
     try:
-        polygon = ensure_numeric(polygon, Float)
+        polygon = ensure_numeric(polygon, numpy.float)
     except NameError, e:
         raise NameError, e
@@ -330 +331 @@
     if len(points.shape) == 1:
         # Only one point was passed in. Convert to array of points
-        points = reshape(points, (1,2))
+        points = numpy.reshape(points, (1,2))
 
     indices, count = separate_points_by_polygon(points, polygon,
@@ -339 +340 @@
     if count == len(indices):
         # No points are outside
-        return array([])
+        return numpy.array([])
     else:
         return indices[count:][::-1]  #return reversed
@@ -354 +355 @@
     #if verbose: print 'Checking input to outside_polygon'
     try:
-        points = ensure_numeric(points, Float)
+        points = ensure_numeric(points, numpy.float)
     except NameError, e:
         raise NameError, e
@@ -362 +363 @@
 
     try:
-        polygon = ensure_numeric(polygon, Float)
+        polygon = ensure_numeric(polygon, numpy.float)
     except NameError, e:
         raise NameError, e
@@ -371 +372 @@
     if len(points.shape) == 1:
         # Only one point was passed in. Convert to array of points
-        points = reshape(points, (1,2))
+        points = numpy.reshape(points, (1,2))
 
 
@@ -439 +440 @@
     assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
 
-
+##    print 'Before: points=%s, flags=%s' % (type(points), str(points.flags))
     try:
-        points = ensure_numeric(points, Float)
+##        points = numpy.ascontiguousarray(ensure_numeric(points, numpy.float))
+        points = ensure_numeric(points, numpy.float)
     except NameError, e:
         raise NameError, e
     except:
         msg = 'Points could not be converted to Numeric array'
-        raise msg
+        raise TypeError, msg
+##    print 'After: points=%s, flags=%s' % (type(points), str(points.flags))
 
     #if verbose: print 'Checking input to separate_points_by_polygon 2'
     try:
-        polygon = ensure_numeric(polygon, Float)
+##        polygon = numpy.ascontiguousarray(ensure_numeric(polygon, numpy.float))
+        polygon = ensure_numeric(polygon, numpy.float)
     except NameError, e:
         raise NameError, e
     except:
         msg = 'Polygon could not be converted to Numeric array'
-        raise msg
+        raise TypeError, msg
 
     msg = 'Polygon array must be a 2d array of vertices'
@@ -472 +476 @@
         # Only one point was passed in.
         # Convert to array of points
-        points = reshape(points, (1,2))
+        points = numpy.reshape(points, (1,2))
 
     
     msg = 'Point array must have two columns (x,y), '
-    msg += 'I got points.shape[1] == %d' %points.shape[0]
+    msg += 'I got points.shape[1] == %d' % points.shape[1]
     assert points.shape[1] == 2, msg
 
@@ -491 +495 @@
 
 
-    indices = zeros( M, Int )
+    indices = numpy.zeros( M, numpy.int )
 
     count = _separate_points_by_polygon(points, polygon, indices,
@@ -618 +622 @@
 
     try:
-        polygon = ensure_numeric(polygon, Float)
+        polygon = ensure_numeric(polygon, numpy.float)
     except NameError, e:
         raise NameError, e
@@ -627 +631 @@
     x = polygon[:,0]
     y = polygon[:,1]
-    x = concatenate((x, [polygon[0,0]]), axis = 0)
-    y = concatenate((y, [polygon[0,1]]), axis = 0)
+    x = numpy.concatenate((x, [polygon[0,0]]), axis = 0)
+    y = numpy.concatenate((y, [polygon[0,1]]), axis = 0)
     
     return x, y
@@ -680 +684 @@
     """
 
-    def __init__(self, regions, default=0.0, geo_reference=None):
+    def __init__(self, regions, default=0.0, geo_reference=None, verbose=False):
 
         try:
@@ -723 +727 @@
             self.regions.append( (P, value) )
 
+        self.verbose = verbose
+
 
 
 
     def __call__(self, x, y):
-        x = array(x).astype(Float)
-        y = array(y).astype(Float)
-
-        N = len(x)
-        assert len(y) == N
-
-        points = concatenate( (reshape(x, (N, 1)),
-                               reshape(y, (N, 1))), axis=1 )
-
+        x = numpy.array(x).astype(numpy.float)
+        y = numpy.array(y).astype(numpy.float)
+
+        assert len(x.shape) == 1 and len(y.shape) == 1
+        
+        N = x.shape[0]
+        assert y.shape[0] == N
+
+        points = numpy.ascontiguousarray(numpy.concatenate( (x[:,numpy.newaxis], y[:,numpy.newaxis]), axis=1 ))
+        
         if callable(self.default):
-            z = self.default(x,y)
+            z = self.default(x, y)
         else:
-            z = ones(N, Float) * self.default
+            z = numpy.ones(N, numpy.float) * self.default
 
         for polygon, value in self.regions:
-            indices = inside_polygon(points, polygon)
+            indices = inside_polygon(points, polygon, verbose=self.verbose)
 
             # FIXME: This needs to be vectorised
             if callable(value):
                 for i in indices:
-                    xx = array([x[i]])
-                    yy = array([y[i]])
+                    xx = numpy.array([x[i]])
+                    yy = numpy.array([y[i]])
                     z[i] = value(xx, yy)[0]
             else:

anuga_core/source/anuga/utilities/polygon_ext.c

@@ -15 +15 @@
 
 #include "Python.h"
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"
 
@@ -222 +222 @@
   }  
   for (k=0; k<M; k++) {
+    x = points[2*k];
+    y = points[2*k + 1];
+
     if (verbose){
-      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
+      if (k %((M+10)/10)==0) printf("Doing %d of %d, x=%f, y=%f\n", k, M, x, y);
     }
     
-    x = points[2*k];
-    y = points[2*k + 1];
-
     inside = 0;
 
@@ -376 +376 @@
     *polygon,
     *indices;
+//  PyObject *xxxx;
 
   int closed, verbose; //Flags
@@ -393 +394 @@
     return NULL;
   }
+  
+//  points = (PyArrayObject *) PyArray_ContiguousFromObject(xxxx, PyArray_DOUBLE, 1, 1);
 
   M = points -> dimensions[0];   //Number of points

anuga_core/source/anuga/utilities/quad.py

@@ +1 @@
+
 """quad.py - quad tree data structure for fast indexing of points in the plane
 
@@ -222 +223 @@
             triangles = {}
             verts = self.retrieve_vertices()
-            # print "verts", verts
+            print "verts", verts
             for vert in verts:
                 triangle_list = self.mesh.get_triangles_and_vertices_per_node(vert)
+                print 'triangle_list=%s' % str(triangle_list)
                 for k, _ in triangle_list:
                     if not triangles.has_key(k):
                         # print 'k',k
-                        tri = self.mesh.get_vertex_coordinates(k,
-                                                               absolute=True)
+                        tri = self.mesh.get_vertex_coordinates(k, absolute=True)
                         n0 = self.mesh.get_normal(k, 0)
                         n1 = self.mesh.get_normal(k, 1)
@@ -435 +436 @@
     """
 
-    from Numeric import minimum, maximum
-
-
     #Make root cell
     #print mesh.coordinates

anuga_core/source/anuga/utilities/sparse.py

@@ +1 @@
+
 """Proof of concept sparse matrix code
 """
 
+import numpy
 
 class Sparse:
@@ -17 +19 @@
             
         if len(args) == 1:
-            from Numeric import array
             try:
-                A = array(args[0])
+                A = numpy.array(args[0])
             except:
                 raise 'Input must be convertable to a Numeric array'
@@ -92 +93 @@
 
     def todense(self):
-        from Numeric import zeros, Float
-
-        D = zeros( (self.M, self.N), Float)
+        D = numpy.zeros( (self.M, self.N), numpy.float)
         
         for i in range(self.M):
@@ -109 +108 @@
         """
 
-        from Numeric import array, zeros, Float
-        
         try:
-            B = array(other)
+            B = numpy.array(other)
         except:
             msg = 'FIXME: Only Numeric types implemented so far'
@@ -130 +127 @@
             assert B.shape[0] == self.N, msg
 
-            R = zeros(self.M, Float) #Result
+            R = numpy.zeros(self.M, numpy.float) #Result
             
             # Multiply nonzero elements
@@ -140 +137 @@
         
             
-            R = zeros((self.M, B.shape[1]), Float) #Result matrix
+            R = numpy.zeros((self.M, B.shape[1]), numpy.float) #Result matrix
 
             # Multiply nonzero elements
@@ -162 +159 @@
         """
 
-        from Numeric import array, zeros, Float
-        
         new = other.copy()
         for key in self.Data.keys():
@@ -177 +172 @@
         """
 
-        from Numeric import array, zeros, Float
-        
         try:
             other = float(other)
@@ -200 +193 @@
         """
 
-        from Numeric import array, zeros, Float
-        
         try:
-            B = array(other)
+            B = numpy.array(other)
         except:
             print 'FIXME: Only Numeric types implemented so far'
@@ -214 +205 @@
             assert B.shape[0] == self.M, 'Mismatching dimensions'
 
-            R = zeros((self.N,), Float) #Result
+            R = numpy.zeros((self.N,), numpy.float) #Result
 
             #Multiply nonzero elements
@@ -251 +242 @@
         """
 
-        from Numeric import array, Float, Int
-
         if isinstance(A,Sparse):
 
-            from Numeric import zeros
             keys = A.Data.keys()
             keys.sort()
             nnz = len(keys)
-            data    = zeros ( (nnz,), Float)
-            colind  = zeros ( (nnz,), Int)
-            row_ptr = zeros ( (A.M+1,), Int)
+            data    = numpy.zeros ( (nnz,), numpy.float)
+            colind  = numpy.zeros ( (nnz,), numpy.int)
+            row_ptr = numpy.zeros ( (A.M+1,), numpy.int)
             current_row = -1
             k = 0
@@ -299 +287 @@
 
     def todense(self):
-        from Numeric import zeros, Float
-
-        D = zeros( (self.M, self.N), Float)
+
+        D = numpy.zeros( (self.M, self.N), numpy.float)
         
         for i in range(self.M):
@@ -314 +301 @@
         """
 
-        from Numeric import array, zeros, Float
-        
         try:
-            B = array(other)
+            B = numpy.array(other)
         except:
             print 'FIXME: Only Numeric types implemented so far'
@@ -333 +318 @@
 if __name__ == '__main__':
     # A little selftest
-    
-    from Numeric import allclose, array, Float 
-    
     A = Sparse(3,3)
 
@@ -366 +348 @@
     u = A*v
     print u
-    assert allclose(u, [6,14,4])
+    assert numpy.allclose(u, [6,14,4])
 
     u = A.trans_mult(v)
     print u
-    assert allclose(u, [6,6,10])
+    assert numpy.allclose(u, [6,6,10])
 
     #Right hand side column
-    v = array([[2,4],[3,4],[4,4]])
+    v = numpy.array([[2,4],[3,4],[4,4]])
 
     u = A*v[:,0]
-    assert allclose(u, [6,14,4])
+    assert numpy.allclose(u, [6,14,4])
 
     #u = A*v[:,1]

anuga_core/source/anuga/utilities/sparse_ext.c

@@ -11 +11 @@
         
 #include "Python.h"
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"
 #include "stdio.h"

anuga_core/source/anuga/utilities/test_cg_solve.py

@@ -5 +5 @@
 import unittest
 
-
-from Numeric import dot, allclose, array, transpose, arange, ones, Float
+import numpy
 from anuga.utilities.cg_solve import *
 from anuga.utilities.cg_solve import _conjugate_gradient
@@ -30 +29 @@
         x = conjugate_gradient(A,b,x,iprint=0)
 
-        assert allclose(x,xe)
+        assert numpy.allclose(x,xe)
 
     def test_max_iter(self):
@@ -61 +60 @@
         A = Sparse(n,n)
 
-        for i in arange(0,n):
+        for i in numpy.arange(0,n):
             A[i,i] = 1.0
             if i > 0 :
@@ -68 +67 @@
                 A[i,i+1] = -0.5
 
-        xe = ones( (n,), Float)
+        xe = numpy.ones( (n,), numpy.float)
 
         b  = A*xe
         x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=1)
 
-        assert allclose(x,xe)
+        assert numpy.allclose(x,xe)
 
     def test_solve_large_2d(self):
@@ -83 +82 @@
         A = Sparse(m*n, m*n)
 
-        for i in arange(0,n):
-            for j in arange(0,m):
+        for i in numpy.arange(0,n):
+            for j in numpy.arange(0,m):
                 I = j+m*i
                 A[I,I] = 4.0
@@ -96 +95 @@
                     A[I,I+1] = -1.0
 
-        xe = ones( (n*m,), Float)
+        xe = numpy.ones( (n*m,), numpy.float)
 
         b  = A*xe
         x = conjugate_gradient(A,b,b,iprint=0)
 
-        assert allclose(x,xe)
+        assert numpy.allclose(x,xe)
 
     def test_solve_large_2d_csr_matrix(self):
@@ -112 +111 @@
         A = Sparse(m*n, m*n)
 
-        for i in arange(0,n):
-            for j in arange(0,m):
+        for i in numpy.arange(0,n):
+            for j in numpy.arange(0,m):
                 I = j+m*i
                 A[I,I] = 4.0
@@ -125 +124 @@
                     A[I,I+1] = -1.0
 
-        xe = ones( (n*m,), Float)
+        xe = numpy.ones( (n*m,), numpy.float)
 
         # Convert to csr format
@@ -134 +133 @@
         x = conjugate_gradient(A,b,b,iprint=20)
 
-        assert allclose(x,xe)
+        assert numpy.allclose(x,xe)
 
 
@@ -145 +144 @@
         A = Sparse(m*n, m*n)
 
-        for i in arange(0,n):
-            for j in arange(0,m):
+        for i in numpy.arange(0,n):
+            for j in numpy.arange(0,m):
                 I = j+m*i
                 A[I,I] = 4.0
@@ -158 +157 @@
                     A[I,I+1] = -1.0
 
-        xe = ones( (n*m,), Float)
+        xe = numpy.ones( (n*m,), numpy.float)
 
         b  = A*xe
         x = conjugate_gradient(A,b)
 
-        assert allclose(x,xe)
+        assert numpy.allclose(x,xe)
 
 
@@ -202 +201 @@
         x = conjugate_gradient(A,b,x,iprint=0)
 
-        assert allclose(x,xe)
+        assert numpy.allclose(x,xe)
 
 #-------------------------------------------------------------

anuga_core/source/anuga/utilities/test_data_audit.py

@@ -3 +3 @@
 
 import unittest
-from Numeric import zeros, array, allclose, Float
 from tempfile import mkstemp
 import os

anuga_core/source/anuga/utilities/test_numerical_tools.py

@@ -3 +3 @@
 
 import unittest
-from Numeric import zeros, array, allclose
-from Numeric import ArrayType, Float, Int, array, alltrue
+import numpy
 
 from math import sqrt, pi
@@ -23 +22 @@
         """Test angles between one vector and the x-axis
         """
-        assert allclose(angle([1.0, 0.0])/pi*180, 0.0)      
-        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
-        assert allclose(angle([0.0, 1.0])/pi*180, 90.0)         
-        assert allclose(angle([-1.0, 1.0])/pi*180, 135.0)               
-        assert allclose(angle([-1.0, 0.0])/pi*180, 180.0)
-        assert allclose(angle([-1.0, -1.0])/pi*180, 225.0)
-        assert allclose(angle([0.0, -1.0])/pi*180, 270.0)
-        assert allclose(angle([1.0, -1.0])/pi*180, 315.0)
+        assert numpy.allclose(angle([1.0, 0.0])/pi*180, 0.0)        
+        assert numpy.allclose(angle([1.0, 1.0])/pi*180, 45.0)
+        assert numpy.allclose(angle([0.0, 1.0])/pi*180, 90.0)           
+        assert numpy.allclose(angle([-1.0, 1.0])/pi*180, 135.0)         
+        assert numpy.allclose(angle([-1.0, 0.0])/pi*180, 180.0)
+        assert numpy.allclose(angle([-1.0, -1.0])/pi*180, 225.0)
+        assert numpy.allclose(angle([0.0, -1.0])/pi*180, 270.0)
+        assert numpy.allclose(angle([1.0, -1.0])/pi*180, 315.0)
                 
                                                           
@@ -37 +36 @@
         """    
         
-        assert allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
-        assert allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
+        assert numpy.allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
+        assert numpy.allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
                 
-        assert allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)    
-        assert allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0)  
+        assert numpy.allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)      
+        assert numpy.allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0)    
         
-        assert allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
-        assert allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)   
+        assert numpy.allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
+        assert numpy.allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)     
         
-        assert allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)   
-        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)
+        assert numpy.allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)     
+        assert numpy.allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)
 
         #From test_get_boundary_polygon_V
         v_prev = [-0.5, -0.5]
         vc = [ 0.0,  -0.5]
-        assert allclose(angle(vc, v_prev)/pi*180, 45.0)
+        assert numpy.allclose(angle(vc, v_prev)/pi*180, 45.0)
 
         vc = [ 0.5,  0.0]
-        assert allclose(angle(vc, v_prev)/pi*180, 135.0)
+        assert numpy.allclose(angle(vc, v_prev)/pi*180, 135.0)
 
         vc = [ -0.5,  0.5]
-        assert allclose(angle(vc, v_prev)/pi*180, 270.0)                
+        assert numpy.allclose(angle(vc, v_prev)/pi*180, 270.0)                
 
 
     def test_anglediff(self):
-        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
+        assert numpy.allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
 
         
@@ -68 +67 @@
         A = [1,2,3,4]
         B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'l'
+        assert isinstance(B, numpy.ndarray)
+        assert B.dtype.char == 'l'
         assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
 
         A = [1,2,3.14,4]
         B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
+        assert isinstance(B, numpy.ndarray)
+        assert B.dtype.char == 'd'
         assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
 
         A = [1,2,3,4]
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
+        B = ensure_numeric(A, numpy.float)
+        assert isinstance(B, numpy.ndarray)
+        assert B.dtype.char == 'd'
         assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
 
         A = [1,2,3,4]
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
+        B = ensure_numeric(A, numpy.float)
+        assert isinstance(B, numpy.ndarray)
+        assert B.dtype.char == 'd'
         assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
 
-        A = array([1,2,3,4])
+        A = numpy.array([1,2,3,4])
         B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'l'        
-        assert alltrue(A == B)    
+        assert isinstance(B, numpy.ndarray)
+        assert B.dtype.char == 'l'        
+        assert numpy.alltrue(A == B)    
         assert A is B   #Same object
 
-        A = array([1,2,3,4])
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'        
-        assert alltrue(A == B)    
+        A = numpy.array([1,2,3,4])
+        B = ensure_numeric(A, numpy.float)
+        assert isinstance(B, numpy.ndarray)
+        assert B.dtype.char == 'd'        
+        assert numpy.alltrue(A == B)    
         assert A is not B   #Not the same object
 
         # Check scalars
         A = 1
-        B = ensure_numeric(A, Float)
+        B = ensure_numeric(A, numpy.float)
         #print A, B[0], len(B), type(B) 
         #print B.shape
-        assert alltrue(A == B)
-
-        B = ensure_numeric(A, Int)        
+        assert numpy.alltrue(A == B)
+
+        B = ensure_numeric(A, numpy.int)        
         #print A, B
         #print B.shape
-        assert alltrue(A == B)
-
+        assert numpy.alltrue(A == B)
+
+    def NO_test_ensure_numeric_char(self):
         # Error situation
-
-        B = ensure_numeric('hello', Int)                
-        assert allclose(B, [104, 101, 108, 108, 111])
+        B = ensure_numeric('hello', numpy.int)
+        print 'B=%s %s' % (str(B), type(B))
+        assert numpy.allclose(B, [104, 101, 108, 108, 111])
 
 
@@ -160 +160 @@
         zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
         a, b = gradient2(x0, y0, x1, y1, z0, z1)
-
+        
         assert zx == a
         assert zy == b
@@ -208 +208 @@
         #There are four elements greater than or equal to 3
         bins = [3]
-        assert allclose(histogram(a, bins), [4])
+        assert numpy.allclose(histogram(a, bins), [4])
 
         bins = [ min(a) ]
-        assert allclose(histogram(a, bins), [len(a)])
+        assert numpy.allclose(histogram(a, bins), [len(a)])
 
         bins = [ max(a)+0.00001 ]
-        assert allclose(histogram(a, bins), [0])        
+        assert numpy.allclose(histogram(a, bins), [0])        
         
         bins = [1,2,3,4]
-        assert allclose(histogram(a, bins), [8,3,3,1])
+        assert numpy.allclose(histogram(a, bins), [8,3,3,1])
 
         bins = [1.1,2,3.1,4]
         #print histogram(a, bins)
-        assert allclose(histogram(a, bins), [0,6,0,1])
+        assert numpy.allclose(histogram(a, bins), [0,6,0,1])
 
         bins = [0,1.5,2,3]
-        assert allclose(histogram(a, bins), [8,0,3,4])
-        assert allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
+        assert numpy.allclose(histogram(a, bins), [8,0,3,4])
+        assert numpy.allclose(histogram(a, [0,3]), histogram(a, [-0.5,3]))
 
         # Check situation with #bins >= #datapoints
         a = [1.7]
         bins = [0,1.5,2,3]
-        assert allclose(histogram(a, bins), [0,1,0,0])
+        assert numpy.allclose(histogram(a, bins), [0,1,0,0])
 
         a = [1.7]
         bins = [0]
-        assert allclose(histogram(a, bins), [1])
+        assert numpy.allclose(histogram(a, bins), [1])
 
         a = [-1.7]
         bins = [0]
-        assert allclose(histogram(a, bins), [0])
+        assert numpy.allclose(histogram(a, bins), [0])
 
         a = [-1.7]
         bins = [-1.7]
-        assert allclose(histogram(a, bins), [1])        
+        assert numpy.allclose(histogram(a, bins), [1])        
         
 
@@ -281 +281 @@
         from util_ext import gradient as gradient_c
 
-        from RandomArray import uniform, seed
-        seed(17, 53)
+        from numpy.random import uniform, seed
+        seed(53)    # numeric was seed(17, 53)
 
         x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
@@ -312 +312 @@
         # sqrt(sum_over_x&y((xi - yi)^2))
         err__1 = err(x,y,2,False)
+        
         assert err__1 == sqrt(20)
         #print "err_", err_

anuga_core/source/anuga/utilities/test_polygon.py

@@ -3 +3 @@
 
 import unittest
-from Numeric import zeros, array, allclose
+import numpy
 from math import sqrt, pi
 from anuga.utilities.numerical_tools import ensure_numeric
@@ -43 +43 @@
         p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
 
-        f = Polygon_function( [(p1, 1.0)] )
+        f = Polygon_function( [(p1, 1.0)], verbose=False )
         z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
-        assert allclose(z, [1,1,0,0])
+        assert numpy.allclose(z, [1,1,0,0])
 
 
         f = Polygon_function( [(p2, 2.0)] )
         z = f([5, 5, 27, 35], [5, 9, 8, -5]) # First and last inside p2
-        assert allclose(z, [2,0,0,2])
+        assert numpy.allclose(z, [2,0,0,2])
 
 
@@ -56 +56 @@
         f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
         z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,1,0,2])
+        assert numpy.allclose(z, [2,1,0,2])
+
+
+    def test_polygon_function_constants_tuple(self):
+        polygon = [[0,0], [10,0], [10,10], [0,10]]
+        points = ((5, 5), (5, 9), (27, 8), (35, -5))
+
+        (indices, count) = separate_points_by_polygon(points, polygon, verbose=False)
 
     def test_polygon_function_csvfile(self):
@@ -72 +79 @@
         z = f([430000,480000], [490000,7720000]) # first outside, second inside
         
-        assert allclose(z, [0,10])
+        assert numpy.allclose(z, [0,10])
 
     def test_polygon_function_georef(self):
@@ -90 +97 @@
         z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
 
-        assert allclose(z, [1,1,0,0])
+        assert numpy.allclose(z, [1,1,0,0])
 
 
         f = Polygon_function( [(p2, 2.0)], geo_reference=geo)
         z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
-        assert allclose(z, [2,0,0,2])
+        assert numpy.allclose(z, [2,0,0,2])
 
 
@@ -101 +108 @@
         f = Polygon_function( [(p1, 1.0), (p2, 2.0)], geo_reference=geo)
         z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,1,0,2])
+        assert numpy.allclose(z, [2,1,0,2])
 
 
@@ -107 +114 @@
         f = Polygon_function( [(p1, 1.0), (p2, 2.0)])
         z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert not allclose(z, [2,1,0,2])        
+        assert not numpy.allclose(z, [2,1,0,2])        
 
 
@@ -120 +127 @@
         f = Polygon_function( [(p1, test_function)] )
         z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
-        assert allclose(z, [10,14,0,0])
+        assert numpy.allclose(z, [10,14,0,0])
 
         # Combined
         f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
         z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,14,0,2])
+        assert numpy.allclose(z, [2,14,0,2])
 
 
@@ -131 +138 @@
         f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
         z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,14,3.14,2])
+        assert numpy.allclose(z, [2,14,3.14,2])
 
 
@@ -138 +145 @@
                               default = test_function)
         z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,14,35,2])
+        assert numpy.allclose(z, [2,14,35,2])
 
 
@@ -210 +217 @@
         res = inside_polygon(points, polygon)
         assert len(res) == 2
-        assert allclose(res, [0,1])
+        assert numpy.allclose(res, [0,1])     ## alltrue?
 
 
@@ -305 +312 @@
         res = inside_polygon( points, polygon, verbose=False )
 
-        assert allclose( res, [0,1,2] )
+        assert numpy.allclose( res, [0,1,2] )
 
     def test_outside_polygon(self):
@@ -317 +324 @@
         
         indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-        assert allclose( indices, [1] )
+        assert numpy.allclose( indices, [1] )
         
         # One more test of vector formulation returning indices
@@ -324 +331 @@
         res = outside_polygon( points, polygon )
 
-        assert allclose( res, [3, 4] )
+        assert numpy.allclose( res, [3, 4] )
 
 
@@ -332 +339 @@
         res = outside_polygon( points, polygon )
 
-        assert allclose( res, [0, 4, 5] )        
+        assert numpy.allclose( res, [0, 4, 5] )        
      
     def test_outside_polygon2(self):
@@ -355 +362 @@
         #print indices, count
         assert count == 0 #None inside
-        assert allclose(indices, [3,2,1,0])
+        assert numpy.allclose(indices, [3,2,1,0])
 
         indices = outside_polygon(points, U, closed = True)
-        assert allclose(indices, [0,1,2,3])
+        assert numpy.allclose(indices, [0,1,2,3])
 
         indices = inside_polygon(points, U, closed = True)
-        assert allclose(indices, [])                
+        assert numpy.allclose(indices, [])                
 
 
@@ -375 +382 @@
         indices, count = separate_points_by_polygon(points, U)
         assert count == 3 #All inside
-        assert allclose(indices, [0,1,2])
+        assert numpy.allclose(indices, [0,1,2])
 
         indices = outside_polygon(points, U, closed = True)
-        assert allclose(indices, [])
+        assert numpy.allclose(indices, [])
 
         indices = inside_polygon(points, U, closed = True)
-        assert allclose(indices, [0,1,2])
+        assert numpy.allclose(indices, [0,1,2])
         
 
@@ -388 +395 @@
 
         indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-        assert allclose( indices, [0,2,1] )
+        assert numpy.allclose( indices, [0,2,1] )
         assert count == 2
         
@@ -396 +403 @@
         res, count = separate_points_by_polygon( points, polygon )
 
-        assert allclose( res, [0,1,2,4,3] )
+        assert numpy.allclose( res, [0,1,2,4,3] )
         assert count == 3
 
@@ -404 +411 @@
         res, count = separate_points_by_polygon( points, polygon )
 
-        assert allclose( res, [1,2,3,5,4,0] )        
+        assert numpy.allclose( res, [1,2,3,5,4,0] )        
         assert count == 3
         
@@ -588 +595 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [0.0, 0.0])
+        assert numpy.allclose(value, [0.0, 0.0])
 
     def test_intersection2(self):
@@ -596 +603 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [12.0, 6.0])
+        assert numpy.allclose(value, [12.0, 6.0])
 
         # Swap direction of one line
@@ -603 +610 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [12.0, 6.0])
+        assert numpy.allclose(value, [12.0, 6.0])
 
         # Swap order of lines
         status, value = intersection(line1, line0)
         assert status == 1
-        assert allclose(value, [12.0, 6.0])        
+        assert numpy.allclose(value, [12.0, 6.0])        
         
     def test_intersection3(self):
@@ -616 +623 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [14.068965517, 7.0344827586])
+        assert numpy.allclose(value, [14.068965517, 7.0344827586])
 
         # Swap direction of one line
@@ -623 +630 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [14.068965517, 7.0344827586])        
+        assert numpy.allclose(value, [14.068965517, 7.0344827586])        
 
         # Swap order of lines
         status, value = intersection(line1, line0)
         assert status == 1        
-        assert allclose(value, [14.068965517, 7.0344827586])        
+        assert numpy.allclose(value, [14.068965517, 7.0344827586])        
 
 
@@ -641 +648 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [1.0, 1.0])
+        assert numpy.allclose(value, [1.0, 1.0])
 
 
@@ -649 +656 @@
         status, value = intersection(line0, line1)
         assert status == 1
-        assert allclose(value, [1.0, 1.0])        
+        assert numpy.allclose(value, [1.0, 1.0])        
         
 
@@ -679 +686 @@
                                                           common_end_point[1])
             msg += ' gave %s, \nbut when reversed we got %s' %(p1, p2)
-            assert allclose(p1, p2), msg
+            assert numpy.allclose(p1, p2), msg
 
             # Swap order of lines
@@ -685 +692 @@
             assert status == 1                        
             msg = 'Order of lines gave different results'
-            assert allclose(p1, p3), msg
+            assert numpy.allclose(p1, p3), msg
             
 
@@ -725 +732 @@
         status, value = intersection(line0, line1)
         assert status == 2
-        assert allclose(value, [[0,0], [3,0]])
+        assert numpy.allclose(value, [[0,0], [3,0]])
 
         # Overlap 2
@@ -733 +740 @@
         status, value = intersection(line0, line1)
         assert status == 2
-        assert allclose(value, [[-3, 0], [5,0]])        
+        assert numpy.allclose(value, [[-3, 0], [5,0]])        
 
         # Inclusion 1
@@ -741 +748 @@
         status, value = intersection(line0, line1)
         assert status == 2        
-        assert allclose(value, line1)
+        assert numpy.allclose(value, line1)
 
         # Inclusion 2
@@ -749 +756 @@
         status, value = intersection(line0, line1)
         assert status == 2        
-        assert allclose(value, line0)                                        
+        assert numpy.allclose(value, line0)                                        
 
 
@@ -768 +775 @@
         status, value = intersection(line0, line1)
         assert status == 2                
-        assert allclose(value, [[1, 7], [7, 19]])
+        assert numpy.allclose(value, [[1, 7], [7, 19]])
 
         status, value = intersection(line1, line0)
         assert status == 2
-        assert allclose(value, [[1, 7], [7, 19]])
+        assert numpy.allclose(value, [[1, 7], [7, 19]])
 
         # Swap direction
@@ -779 +786 @@
         status, value = intersection(line0, line1)
         assert status == 2
-        assert allclose(value, [[7, 19], [1, 7]])
+        assert numpy.allclose(value, [[7, 19], [1, 7]])
 
         line0 = [[0,5], [7,19]]
@@ -785 +792 @@
         status, value = intersection(line0, line1)
         assert status == 2
-        assert allclose(value, [[1, 7], [7, 19]])        
+        assert numpy.allclose(value, [[1, 7], [7, 19]])        
         
 
@@ -793 +800 @@
         status, value = intersection(line0, line1)
         assert status == 2                        
-        assert allclose(value, [[1,7], [7, 19]])                
+        assert numpy.allclose(value, [[1,7], [7, 19]])                
 
         line0 = [[0,5], [10,25]]
@@ -799 +806 @@
         status, value = intersection(line0, line1)
         assert status == 2                        
-        assert allclose(value, [[1,7], [7, 19]])
+        assert numpy.allclose(value, [[1,7], [7, 19]])
 
 
@@ -806 +813 @@
         status, value = intersection(line0, line1)
         assert status == 2                        
-        assert allclose(value, [[7, 19], [1, 7]])                       
+        assert numpy.allclose(value, [[7, 19], [1, 7]])                       
         
         
@@ -843 +850 @@
         res = inside_polygon(points, polygon)
         assert len(res) == 2
-        assert allclose(res, [0,1])
+        assert numpy.allclose(res, [0,1])
 
     def test_polygon_area(self):
@@ -990 +997 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Polygon,'test')
-    #suite = unittest.makeSuite(Test_Polygon,'test_inside_polygon_geo_ref')
+##    suite = unittest.makeSuite(Test_Polygon,'test_polygon_function_constantsX')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/utilities/test_quad.py

@@ -1 +1 @@
 import unittest
-from Numeric import array, allclose
+##import numpy
 
 from quad import Cell, build_quadtree
@@ -236 +236 @@
                      [ 0.,  -1.]]
                      
-        # assert allclose(array(results),[[[ 2.,  1.],
+        # assert numpy.allclose(numpy.array(results),[[[ 2.,  1.],
         #[ 4.,  1.], [ 4.,  4.]], [[ 4.,  1.],[ 5.,  4.],[ 4.,  4.]]] )
         results = Q.search(5,4.)
@@ -253 +253 @@
                      [ 5.,  4.],
                      [ 4.,  4.]]
-        #assert allclose(array(results),[[[ 2.,  1.],[ 4.,  1.], [ 4.,  4.]]
+        #assert numpy.allclose(numpy.array(results),[[[ 2.,  1.],[ 4.,  1.], [ 4.,  4.]]
          #                               ,[[ 2.,  1.],[ 4.,  4.], [ 2.,  4.]],
         #[[ 4.,  1.],  [ 5.,  4.], [ 4.,  4.]],

anuga_core/source/anuga/utilities/test_sparse.py

@@ -5 +5 @@
 
 from sparse import *
-from Numeric import allclose, array, transpose, Float
+import numpy
 
 class Test_Sparse(unittest.TestCase):
@@ -43 +43 @@
         C = Sparse(B)
 
-        assert allclose(C.todense(), B)
+        assert numpy.allclose(C.todense(), B)
 
 
@@ -50 +50 @@
         A[1,1] = 4
 
-        assert allclose(A.todense(), [[0,0,0], [0,4,0], [0,0,0], [0,0,0]])
+        assert numpy.allclose(A.todense(), [[0,0,0], [0,4,0], [0,0,0], [0,0,0]])
 
 
@@ -63 +63 @@
 
         assert len(A) == 0
-        assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
+        assert numpy.allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
 
         #Set an existing zero element to zero
         A[1,2] = 0
         assert len(A) == 0
-        assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
+        assert numpy.allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
 
     def test_sparse_multiplication_vector(self):
@@ -82 +82 @@
 
         u = A*v
-        assert allclose(u, [6,14,4])
+        assert numpy.allclose(u, [6,14,4])
 
         #Right hand side column
-        v = array([[2,4],[3,4],[4,4]])
+        v = numpy.array([[2,4],[3,4],[4,4]])
 
         u = A*v[:,0]
-        assert allclose(u, [6,14,4])
+        assert numpy.allclose(u, [6,14,4])
 
         u = A*v[:,1]
-        assert allclose(u, [12,16,4])
+        assert numpy.allclose(u, [12,16,4])
 
 
@@ -103 +103 @@
 
         #Right hand side matrix
-        v = array([[2,4],[3,4],[4,4]])
+        v = numpy.array([[2,4],[3,4],[4,4]])
 
         u = A*v
-        assert allclose(u, [[6,12], [14,16], [4,4]])
+        assert numpy.allclose(u, [[6,12], [14,16], [4,4]])
 
 
@@ -122 +122 @@
 
         u = A.trans_mult(v)
-        assert allclose(u, [6,6,10])
+        assert numpy.allclose(u, [6,6,10])
 
 
@@ -137 +137 @@
 
         B = 3*A
-        assert allclose(B.todense(), 3*A.todense())
+        assert numpy.allclose(B.todense(), 3*A.todense())
 
         B = A*3
-        assert allclose(B.todense(), 3*A.todense())
+        assert numpy.allclose(B.todense(), 3*A.todense())
 
         try:
@@ -166 +166 @@
         C = A+B
 
-        assert allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]])
+        assert numpy.allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]])
 
     def test_sparse_tocsr(self):
@@ -190 +190 @@
         C = [1, 2, 3]
 
-        assert allclose(B*C, [15.0, 10.0 ,8.0, 0.0])
+        assert numpy.allclose(B*C, [15.0, 10.0 ,8.0, 0.0])
 
         C2 = [[1,2],[2,4],[3,6]]
@@ -196 +196 @@
         #print B*C2
 
-        assert allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]])
+        assert numpy.allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]])
 
 

anuga_core/source/anuga/utilities/test_system_tools.py

@@ -1 +1 @@
 #!/usr/bin/env python
 
-
 import unittest
-from Numeric import zeros, array, allclose, Float
+import numpy
 import zlib
 from os.path import join, split, sep
@@ -81 +80 @@
             pass
         else:
-            test_array = array([[7.0, 3.14], [-31.333, 0.0]])
+            test_array = numpy.array([[7.0, 3.14], [-31.333, 0.0]])
 
             # First file
@@ -87 +86 @@
             fid = NetCDFFile(filename1, 'w')
             fid.createDimension('two', 2)
-            fid.createVariable('test_array', Float,
-                               ('two', 'two'))
+            fid.createVariable('test_array', numpy.dtype(numpy.float).char, ('two', 'two'))
             fid.variables['test_array'][:] = test_array
             fid.close()
@@ -96 +94 @@
             fid = NetCDFFile(filename2, 'w')
             fid.createDimension('two', 2)
-            fid.createVariable('test_array', Float,
-                               ('two', 'two'))
+            fid.createVariable('test_array', numpy.dtype(numpy.float).char, ('two', 'two'))
             fid.variables['test_array'][:] = test_array
             fid.close()

anuga_core/source/anuga/utilities/test_xml_tools.py

@@ -3 +3 @@
 
 import unittest
-from Numeric import zeros, array, allclose, Float
 from tempfile import mkstemp, mktemp
 

anuga_core/source/anuga/utilities/util_ext.c

@@ -16 +16 @@
 
 #include "Python.h"
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"
 

anuga_core/source/anuga/utilities/util_ext.h

@@ -11 +11 @@
         
 #include "Python.h"     
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"