Changeset 2897 for inundation/fit_interpolate

Timestamp:

May 17, 2006, 4:13:37 PM (19 years ago)

Author:

duncan

Message:

close to finishing fit, which replaces least squares

Location:

inundation/fit_interpolate

Files:

• fit.py (modified) (15 diffs)
• general_fit_interpolate.py (modified) (2 diffs)
• run_long_benchmark.py (modified) (1 diff)
• spike_least_squares.py (modified) (1 diff)
• spike_test_least_squares.py (modified) (1 diff)
• test_fit.py (modified) (5 diffs)

inundation/fit_interpolate/fit.py

@@ -16 +16 @@
    Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
    Geoscience Australia, 2004.
+
+   TO DO
+   * test geo_ref, geo_spatial
 """
 
-from geospatial_data.geospatial_data import Geospatial_data
+from Numeric import zeros, Float
+
+from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
 from fit_interpolate.general_fit_interpolate import FitInterpolate
+from utilities.sparse import Sparse, Sparse_CSR
+from utilities.polygon import in_and_outside_polygon
+from fit_interpolate.search_functions import search_tree_of_vertices
+from utilities.cg_solve import conjugate_gradient
+from utilities.numerical_tools import ensure_numeric, gradient
+
+import exceptions
+class ToFewPointsError(exceptions.Exception): pass
 
 DEFAULT_ALPHA = 0.001
@@ -66 +79 @@
 
         if alpha is None:
+
             self.alpha = DEFAULT_ALPHA
         else:    
             self.alpha = alpha
-        
         FitInterpolate.__init__(self,
                  vertex_coordinates,
@@ -79 +92 @@
         m = self.mesh.coordinates.shape[0] #Nbr of basis functions (vertices)
         
-        #Build Atz and AtA matrix
-        self.AtA = Sparse(m,m)
-        self.Atz = zeros((m,att_num), Float)
-      
-
+        self.AtA = None
+        self.Atz = None
+
+        self.point_count = 0
+        if self.alpha <> 0:
+            if verbose: print 'Building smoothing matrix'
+            self._build_smoothing_matrix_D()
+            
     def _build_coefficient_matrix_B(self,
                                   verbose = False):
-        """Build final coefficient matrix"""
-
+        """
+        Build final coefficient matrix
+
+        Precon
+        If alpha is not zero, matrix D has been built
+        Matrix Ata has been built
+        """
+
+        if self.alpha <> 0:
+            #if verbose: print 'Building smoothing matrix'
+            #self._build_smoothing_matrix_D()
+            self.B = self.AtA + self.alpha*self.D
+        else:
+            self.B = self.AtA
+
+        #Convert self.B matrix to CSR format for faster matrix vector
+        self.B = Sparse_CSR(self.B)
 
     def _build_smoothing_matrix_D(self):
@@ -114 +145 @@
         are obtained by computing the gradient a_k, b_k for basis function k
         """
+        
+        #FIXME: algorithm might be optimised by computing local 9x9
+        #"element stiffness matrices:
+
+        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+
+        self.D = Sparse(m,m)
+
+        #For each triangle compute contributions to D = D1+D2
+        for i in range(len(self.mesh)):
+
+            #Get area
+            area = self.mesh.areas[i]
+
+            #Get global vertex indices
+            v0 = self.mesh.triangles[i,0]
+            v1 = self.mesh.triangles[i,1]
+            v2 = self.mesh.triangles[i,2]
+
+            #Get the three vertex_points
+            xi0 = self.mesh.get_vertex_coordinate(i, 0)
+            xi1 = self.mesh.get_vertex_coordinate(i, 1)
+            xi2 = self.mesh.get_vertex_coordinate(i, 2)
+
+            #Compute gradients for each vertex
+            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
+                              1, 0, 0)
+
+            a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
+                              0, 1, 0)
+
+            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
+                              0, 0, 1)
+
+            #Compute diagonal contributions
+            self.D[v0,v0] += (a0*a0 + b0*b0)*area
+            self.D[v1,v1] += (a1*a1 + b1*b1)*area
+            self.D[v2,v2] += (a2*a2 + b2*b2)*area
+
+            #Compute contributions for basis functions sharing edges
+            e01 = (a0*a1 + b0*b1)*area
+            self.D[v0,v1] += e01
+            self.D[v1,v0] += e01
+
+            e12 = (a1*a2 + b1*b2)*area
+            self.D[v1,v2] += e12
+            self.D[v2,v1] += e12
+
+            e20 = (a2*a0 + b2*b0)*area
+            self.D[v2,v0] += e20
+            self.D[v0,v2] += e20
+
+
+    def get_D(self):
+        return self.D.todense()
+
 
     def _build_matrix_AtA_Atz(self,
@@ -131 +218 @@
         z and points are numeric
         Point_coordindates and mesh vertices have the same origin.
-        """
-        
+
+        The number of attributes of the data points does not change
+        """
+        #Build n x m interpolation matrix
+
+        if self.AtA == None:
+            # AtA and Atz need ot be initialised.
+            m = self.mesh.coordinates.shape[0] #Nbr of vertices
+            if len(z.shape) > 1:
+                att_num = z.shape[1]
+                self.Atz = zeros((m,att_num), Float)
+            else:
+                att_num = 1
+                self.Atz = zeros((m,), Float)
+            assert z.shape[0] == point_coordinates.shape[0] 
+
+            self.AtA = Sparse(m,m)
+            #self.Atz = zeros((m,att_num), Float)
+        self.point_count += point_coordinates.shape[0]
+        #print "_build_matrix_AtA_Atz - self.point_count", self.point_count
         if verbose: print 'Getting indices inside mesh boundary'
-        #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon() 
+        #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon()
         self.inside_poly_indices, self.outside_poly_indices  = \
                      in_and_outside_polygon(point_coordinates,
@@ -143 +248 @@
 
         
+        n = len(self.inside_poly_indices)
         #Compute matrix elements for points inside the mesh
         for i in self.inside_poly_indices:
@@ -162 +268 @@
 
                 for j in js:
-                    self.Atz[j] +=  sigmas[j]*z[i]                  
+                    self.Atz[j] +=  sigmas[j]*z[i]
+                    #print "self.Atz building", self.Atz
+                    #print "self.Atz[j]", self.Atz[j]
+                    #print " sigmas[j]", sigmas[j]
+                    #print "z[i]",z[i]
+                    #print "result", sigmas[j]*z[i]
+                    
                     for k in js:
-                        if interp_only == False:
-                            self.AtA[j,k] += sigmas[j]*sigmas[k]
+                        self.AtA[j,k] += sigmas[j]*sigmas[k]
             else:
                 msg = 'Could not find triangle for point', x 
@@ -171 +282 @@
     
         
-    def fit(self, point_coordinates=point_coordinates, z=z):
+    def fit(self, point_coordinates=None, z=None,
+                              verbose = False,
+                              point_origin = None):
         """Fit a smooth surface to given 1d array of data points z.
 
@@ -184 +297 @@
           
         """
-        # build ata and atz
-        # solve fit
-
-        
-    def build_fit_subset(self, point_coordinates, z):
+        if point_coordinates is None:
+            assert self.AtA <> None
+            assert self.Atz <> None
+            #FIXME (DSG) - do  a message
+        else:
+            point_coordinates = ensure_absolute(point_coordinates,
+                                                geo_reference=point_origin)
+            self.build_fit_subset(point_coordinates, z, verbose)
+       
+        #Check sanity
+        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+        n = self.point_count
+        if n<m and self.alpha == 0.0:
+            msg = 'ERROR (least_squares): Too few data points\n'
+            msg += 'There are only %d data points and alpha == 0. ' %n
+            msg += 'Need at least %d\n' %m
+            msg += 'Alternatively, set smoothing parameter alpha to a small '
+            msg += 'positive value,\ne.g. 1.0e-3.'
+            raise ToFewPointsError(msg)
+
+        self._build_coefficient_matrix_B(verbose)
+        return conjugate_gradient(self.B, self.Atz, self.Atz,
+                                  imax=2*len(self.Atz) )
+
+        
+    def build_fit_subset(self, point_coordinates, z,
+                              verbose = False):
         """Fit a smooth surface to given 1d array of data points z.
 
@@ -202 +337 @@
         """
         #Note: Don't get the z info from Geospatial_data.attributes yet.
-        # That means fit has to handle attribute title info.
+        # If we did fit would have to handle attribute title info.
 
         #FIXME(DSG-DSG): Check that the vert and point coords
         #have the same zone.
         if isinstance(point_coordinates,Geospatial_data):
-            point_coordinates = vertex_coordinates.get_data_points( \
+            point_coordinates = point_coordinates.get_data_points( \
                 absolute = True)
         
@@ -214 +349 @@
         point_coordinates = ensure_numeric(point_coordinates, Float)
 
+        self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
 
 
@@ -225 +361 @@
                 verbose = False,
                 acceptable_overshoot = 1.01,
-                expand_search = False,
+                mesh_origin = None,
                 data_origin = None,
-                mesh_origin = None,
-                precrop = False,
                 use_cache = False):
     """
@@ -236 +370 @@
 
         Inputs:
-
-          vertex_coordinates: List of coordinate pairs [xi, eta] of points
-          constituting mesh (or a an m x 2 Numeric array)
+        vertex_coordinates: List of coordinate pairs [xi, eta] of
+              points constituting a mesh (or an m x 2 Numeric array or
+              a geospatial object)
+              Points may appear multiple times
+              (e.g. if vertices have discontinuities)
 
           triangles: List of 3-tuples (or a Numeric array) of
@@ -252 +388 @@
           as min(z) - acceptable_overshoot*delta z and upper limit
           as max(z) + acceptable_overshoot*delta z
+
+          mesh_origin: A geo_reference object or 3-tuples consisting of
+              UTM zone, easting and northing.
+              If specified vertex coordinates are assumed to be
+              relative to their respective origins.
           
 
-          point_attributes: Vector or array of data at the point_coordinates.
-
-          data_origin and mesh_origin are 3-tuples consisting of
-          UTM zone, easting and northing. If specified
-          point coordinates and vertex coordinates are assumed to be
-          relative to their respective origins.
+          point_attributes: Vector or array of data at the
+                            point_coordinates.
 
     """
-
+    #Since this is a wrapper for fit, lets handle the geo_spatial att's
+    
     if use_cache is True:
         interp = cache(_fit,
@@ -268 +406 @@
                         triangles),
                        {'verbose': verbose,
-                        'mesh_origin': mesh_origin},
+                        'mesh_origin': mesh_origin,
+                        'alpha':alpha},
                        verbose = verbose)        
         
     else:
-        interp = Interpolation(vertex_coordinates,
-                               triangles,
-                               verbose = verbose,
-                               mesh_origin = mesh_origin)
-        
-    vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
-
-        
-#                               point_coordinates,
-#                               data_origin = data_origin,
-
-
-    #Sanity check
-    point_coordinates = ensure_numeric(point_coordinates)
-    vertex_coordinates = ensure_numeric(vertex_coordinates)
-
-    #Data points
-    X = point_coordinates[:,0]
-    Y = point_coordinates[:,1]  
-    Z = ensure_numeric(point_attributes)
-    if len(Z.shape) == 1:
-        Z = Z[:, NewAxis]
-        
-
-    #Data points inside mesh boundary
-    indices = interp.point_indices
-    if indices is not None:    
-        Xc = take(X, indices)
-        Yc = take(Y, indices)   
-        Zc = take(Z, indices)
-    else:
-        Xc = X
-        Yc = Y  
-        Zc = Z        
+        interp = Fit(vertex_coordinates,
+                     triangles,
+                     verbose = verbose,
+                     mesh_origin = mesh_origin,
+                     alpha=alpha)
+        
+    vertex_attributes = interp.fit(point_coordinates,
+                                   point_attributes,
+                                   point_origin = data_origin,
+                                   verbose = verbose)
+
+        
+    # Add the value checking stuff that's in least squares.
+    # Maybe this stuff should get pushed down into Fit.
+    # at least be a method of Fit.
+    # Or intigrate it into the fit method, saving teh max and min's
+    # as att's.
     
-    #Vertex coordinates
-    Xi = vertex_coordinates[:,0]
-    Eta = vertex_coordinates[:,1]       
-    Zeta = ensure_numeric(vertex_attributes)
-    if len(Zeta.shape) == 1:
-        Zeta = Zeta[:, NewAxis]    
-
-    for i in range(Zeta.shape[1]): #For each attribute
-        zeta = Zeta[:,i]
-        z = Z[:,i]                
-        zc = Zc[:,i]
-
-        max_zc = max(zc)
-        min_zc = min(zc)
-        delta_zc = max_zc-min_zc
-        upper_limit = max_zc + delta_zc*acceptable_overshoot
-        lower_limit = min_zc - delta_zc*acceptable_overshoot        
-        
-
-        if max(zeta) > upper_limit or min(zeta) < lower_limit:
-            msg = 'Least sqares produced values outside the allowed '
-            msg += 'range [%f, %f].\n' %(lower_limit, upper_limit)
-            msg += 'z in [%f, %f], zeta in [%f, %f].\n' %(min_zc, max_zc,
-                                                          min(zeta), max(zeta))
-            msg += 'If greater range is needed, increase the value of '
-            msg += 'acceptable_fit_overshoot (currently %.2f).\n' %(acceptable_overshoot)
-
-
-            offending_vertices = (zeta > upper_limit or zeta < lower_limit)
-            Xi_c = compress(offending_vertices, Xi)
-            Eta_c = compress(offending_vertices, Eta)
-            offending_coordinates = concatenate((Xi_c[:, NewAxis],
-                                                 Eta_c[:, NewAxis]),
-                                                axis=1)
-
-            msg += 'Offending locations:\n %s' %(offending_coordinates)
-            
-            raise FittingError, msg
-
-
-    
-        if verbose:
-            print '+------------------------------------------------'
-            print 'Least squares statistics'
-            print '+------------------------------------------------'   
-            print 'points: %d points' %(len(z))
-            print '    x in [%f, %f]'%(min(X), max(X))
-            print '    y in [%f, %f]'%(min(Y), max(Y))
-            print '    z in [%f, %f]'%(min(z), max(z))
-            print
-
-            if indices is not None:
-                print 'Cropped points: %d points' %(len(zc))
-                print '    x in [%f, %f]'%(min(Xc), max(Xc))
-                print '    y in [%f, %f]'%(min(Yc), max(Yc))
-                print '    z in [%f, %f]'%(min(zc), max(zc))
-                print
-            
-
-            print 'Mesh: %d vertices' %(len(zeta))
-            print '    xi in [%f, %f]'%(min(Xi), max(Xi))
-            print '    eta in [%f, %f]'%(min(Eta), max(Eta))
-            print '    zeta in [%f, %f]'%(min(zeta), max(zeta))
-            print '+------------------------------------------------'
-
     return vertex_attributes
 

inundation/fit_interpolate/general_fit_interpolate.py

@@ -32 +32 @@
 from utilities.numerical_tools import ensure_numeric, mean, INF
 from utilities.polygon import in_and_outside_polygon
-from geospatial_data.geospatial_data import Geospatial_data
+from geospatial_data.geospatial_data import Geospatial_data, \
+     ensure_absolute
 from search_functions import search_tree_of_vertices
 
@@ -75 +76 @@
         #Convert input to Numeric arrays
         triangles = ensure_numeric(triangles, Int)
-        
-        if isinstance(vertex_coordinates,Geospatial_data):
-            vertex_coordinates = vertex_coordinates.get_data_points( \
-                absolute = True)
-            msg = "Use a Geospatial_data object or a mesh origin. Not both."
-            assert mesh_origin == None, msg
-            
-        else:
-            vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
-        #Build underlying mesh
-        if verbose: print 'Building mesh'
-        
-        if mesh_origin is None:
-            geo = None #Geo_reference()
-        else:
-            if isinstance(mesh_origin, Geo_reference):
-                geo = mesh_origin
-            else:
-                geo = Geo_reference(mesh_origin[0],
-                                    mesh_origin[1],
-                                    mesh_origin[2])
-            vertex_coordinates = geo.get_absolute(vertex_coordinates)
+        vertex_coordinates = ensure_absolute(vertex_coordinates,
+                                         geo_reference = mesh_origin)
+
         #Don't pass geo_reference to mesh.  It doesn't work.
         self.mesh = Mesh(vertex_coordinates, triangles)

inundation/fit_interpolate/run_long_benchmark.py

@@ -8 +8 @@
 
 use_least_squares_list = [True,False]
-is_fit_list = [False]
+is_fit_list = [True]
 num_of_points_list = [10]
 maxArea_list = [0.1, 0.001]

inundation/fit_interpolate/spike_least_squares.py

@@ -565 +565 @@
         self.AtA = Sparse(m,m)
         self.Atz = zeros((m,att_num), Float)
+        print "self.Atz",self.Atz 
+        self.Atz[0] = 100
+        print "self.Atz",self.Atz 
 
         #Build quad tree of vertices

inundation/fit_interpolate/spike_test_least_squares.py

@@ -12 +12 @@
 from utilities.sparse import Sparse, Sparse_CSR
 
-from coordinate_transforms.geo_reference import Geo_reference
-
-def distance(x, y):
-    return sqrt( sum( (array(x)-array(y))**2 ))
+
 
 def linear_function(point):

inundation/fit_interpolate/test_fit.py

@@ -6 +6 @@
 from math import sqrt
 
-
-from spike_least_squares import *
-from Numeric import allclose, array, transpose
-from Numeric import zeros, take, compress, array, Float, Int, dot, transpose, concatenate, ArrayType
+from Numeric import zeros, take, compress, Float, Int, dot, concatenate, \
+     ArrayType, allclose, array
+
+from fit import *
 from utilities.sparse import Sparse, Sparse_CSR
-
 from coordinate_transforms.geo_reference import Geo_reference
+from utilities.numerical_tools import ensure_numeric
+from geospatial_data.geospatial_data import Geospatial_data
 
 def distance(x, y):
@@ -22 +23 @@
 
 
-class Test_Least_Squares(unittest.TestCase):
+class Test_Fit(unittest.TestCase):
 
     def setUp(self):
@@ -29 +30 @@
     def tearDown(self):
         pass
+
 
     def test_smooth_attributes_to_mesh(self):
@@ -46 +48 @@
 
         z = [z1, z2, z3]
-        interp = Interpolation(points, triangles, z,data_coords, alpha=0)
+        fit = Fit(points, triangles, alpha=0)
         #print "interp.get_A()", interp.get_A()
-        A = interp.A
-        #print "A",A
+        fit._build_matrix_AtA_Atz(ensure_numeric(data_coords),
+                                  ensure_numeric(z))
+        #print "Atz - from fit", fit.Atz
+        #print "AtA - from fit", fit.AtA.todense()
         #print "z",z 
-        Atz = A.trans_mult(z)
-        #print "Atz",Atz
-        
-        f = interp.fit(z)
+
+        assert allclose(fit.Atz, [2.8, 3.6, 3.6], atol=1e-7)
+
+        f = fit.fit()
+        
         answer = [0, 5., 5.]
 
@@ -62 +67 @@
         assert allclose(f, answer, atol=1e-7)
 
-       
+    def test_smooth_att_to_meshII(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 2.0]
+        d3 = [3.0,1.0]
+        data_coords = [d1, d2, d3]
+        z = linear_function(data_coords)
+        #print "z",z 
+
+        interp = Fit(points, triangles, alpha=0.0)
+        f = interp.fit(data_coords, z)
+        answer = linear_function(points)
+        #print "f\n",f
+        #print "answer\n",answer
+
+        assert allclose(f, answer)
+
+    def test_smooth_attributes_to_meshIII(self):
+
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0,1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0,-2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
+
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [0.0,2.0],
+                        [1.0, 1.0],
+                        [2.0, 1.0],
+                        [0.0,0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2,-0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [3.0,1.0]]
+
+        z = linear_function(point_coords)
+        interp = Fit(vertices, triangles,
+                                alpha=0.0)
+
+        #print 'z',z
+        f = interp.fit(point_coords,z)
+        answer = linear_function(vertices)
+        #print "f\n",f
+        #print "answer\n",answer
+        assert allclose(f, answer)
+
+   
+    def test_smooth_attributes_to_meshIV(self):
+        """ Testing 2 attributes smoothed to the mesh
+        """
+
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 3.0]
+        d3 = [3.0, 1.0]
+        z1 = [2, 4]
+        z2 = [4, 8]
+        z3 = [4, 8]
+        data_coords = [d1, d2, d3]
+
+        z = [z1, z2, z3]
+        fit = Fit(points, triangles, alpha=0)
+        
+        f =  fit.fit(data_coords,z)
+        answer = [[0,0], [5., 10.], [5., 10.]]
+        assert allclose(f, answer)
+
+    def test_smooth_attributes_to_mesh_build_fit_subset(self):
+
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0,1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0,-2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
+
+        interp = Fit(vertices, triangles,
+                                alpha=0.0)
+
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [0.0,2.0],
+                        [1.0, 1.0],
+                       ]
+
+        z = linear_function(point_coords)
+
+        f = interp.build_fit_subset(point_coords,z)
+        
+        point_coords = [
+                        [2.0, 1.0],
+                        [0.0,0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2,-0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [3.0,1.0]]
+        
+        z = linear_function(point_coords)
+
+        f = interp.build_fit_subset(point_coords,z)
+        
+        #print 'z',z
+        f = interp.fit()
+        answer = linear_function(vertices)
+        #print "f\n",f
+        #print "answer\n",answer
+        assert allclose(f, answer)
+
+    def test_fit_and_interpolation(self):
+        from mesh import Mesh
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Get (enough) datapoints
+        data_points = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667],
+                       [ 0.0, 1.0],
+                       [ 0.0, 3.0],
+                       [ 1.0, 0.0],
+                       [ 1.0, 1.0],
+                       [ 1.0, 2.0],
+                       [ 1.0, 3.0],
+                       [ 2.0, 1.0],
+                       [ 3.0, 0.0],
+                       [ 3.0, 1.0]]
+
+        z = linear_function(data_points)
+        interp = Fit(points, triangles,
+                                alpha=0.0)
+
+        answer = linear_function(points)
+
+        f = interp.fit(data_points, z)
+
+        #print "f",f
+        #print "answer",answer
+        assert allclose(f, answer)
+
+        
+    def test_smoothing_and_interpolation(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Get (too few!) datapoints
+        data_points = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667]]
+
+        z = linear_function(data_points)
+        answer = linear_function(points)
+
+        #Make interpolator with too few data points and no smoothing
+        interp = Fit(points, triangles, alpha=0.0)
+        #Must raise an exception
+        try:
+            f = interp.fit(data_points,z)
+        except ToFewPointsError:
+            pass
+
+        #Now try with smoothing parameter
+        interp = Fit(points, triangles, alpha=1.0e-13)
+
+        f = interp.fit(data_points,z)
+        #f will be different from answer due to smoothing
+        assert allclose(f, answer,atol=5)
+
+
+    #Tests of smoothing matrix
+    def test_smoothing_matrix_one_triangle(self):
+        from Numeric import dot
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        points = [a, b, c]
+
+        vertices = [ [1,0,2] ]   #bac
+
+        interp = Fit(points, vertices)
+
+        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
+                                   [-0.5, 0.5, 0],
+                                   [-0.5, 0, 0.5]])
+
+        #Define f(x,y) = x
+        f = array([0,0,2]) #Value at global vertex 2
+
+        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
+        #           int 1 dx dy = area = 2
+        assert dot(dot(f, interp.get_D()), f) == 2
+
+        #Define f(x,y) = y
+        f = array([0,2,0])  #Value at global vertex 1
+
+        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
+        #           int 1 dx dy = area = 2
+        assert dot(dot(f, interp.get_D()), f) == 2
+
+        #Define f(x,y) = x+y
+        f = array([0,2,2])  #Values at global vertex 1 and 2
+
+        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
+        #           int 2 dx dy = 2*area = 4
+        assert dot(dot(f, interp.get_D()), f) == 4
+
+
+    def test_smoothing_matrix_more_triangles(self):
+        from Numeric import dot
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        interp = Fit(points, vertices)
+
+
+        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
+        #                           [-0.5, 0.5, 0],
+        #                           [-0.5, 0, 0.5]])
+
+        #Define f(x,y) = x
+        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
+
+        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
+        #           int 1 dx dy = total area = 8
+        assert dot(dot(f, interp.get_D()), f) == 8
+
+        #Define f(x,y) = y
+        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
+
+        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
+        #           int 1 dx dy = area = 8
+        assert dot(dot(f, interp.get_D()), f) == 8
+
+        #Define f(x,y) = x+y
+        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
+
+        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
+        #           int 2 dx dy = 2*area = 16
+        assert dot(dot(f, interp.get_D()), f) == 16
+
+
+
+    def test_fit_and_interpolation_with_different_origins(self):
+        """Fit a surface to one set of points. Then interpolate that surface
+        using another set of points.
+        This test tests situtaion where points and mesh belong to a different
+        coordinate system as defined by origin.
+        """
+        from mesh import Mesh
+
+        #Setup mesh used to represent fitted function
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Datapoints to fit from
+        data_points1 = [[ 0.66666667, 0.66666667],
+                        [ 1.33333333, 1.33333333],
+                        [ 2.66666667, 0.66666667],
+                        [ 0.66666667, 2.66666667],
+                        [ 0.0, 1.0],
+                        [ 0.0, 3.0],
+                        [ 1.0, 0.0],
+                        [ 1.0, 1.0],
+                        [ 1.0, 2.0],
+                        [ 1.0, 3.0],
+                        [ 2.0, 1.0],
+                        [ 3.0, 0.0],
+                        [ 3.0, 1.0]]
+
+
+        #First check that things are OK when using same origin
+        mesh_origin = (56, 290000, 618000) #zone, easting, northing
+        data_origin = (56, 290000, 618000) #zone, easting, northing
+
+
+        #Fit surface to mesh
+        interp = Fit(points, triangles, 
+                               alpha=0.0,
+                               mesh_origin = mesh_origin)
+
+        data_geo_spatial = Geospatial_data(data_points1,
+                         geo_reference = Geo_reference(56, 290000, 618000))
+        z = linear_function(data_points1) #Example z-values
+        f = interp.fit(data_geo_spatial, z)   #Fitted values at vertices
+
+        #Shift datapoints according to new origins
+        for k in range(len(data_points1)):
+            data_points1[k][0] += mesh_origin[1] - data_origin[1]
+            data_points1[k][1] += mesh_origin[2] - data_origin[2]
+
+
+
+        #Fit surface to mesh
+        interp = Fit(points, triangles, 
+                               alpha=0.0) #,
+                               # mesh_origin = mesh_origin)
+
+        f1 = interp.fit(data_points1,z) #Fitted values at vertices (using same z as before)
+
+        assert allclose(f,f1), 'Fit should have been unaltered'
+
+
+    def test_smooth_attributes_to_mesh_function(self):
+        """ Testing 2 attributes smoothed to the mesh
+        """
+
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 3.0]
+        d3 = [3.0, 1.0]
+        z1 = [2, 4]
+        z2 = [4, 8]
+        z3 = [4, 8]
+        data_coords = [d1, d2, d3]
+        z = [z1, z2, z3]
+
+        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
+        answer = [[0, 0], [5., 10.], [5., 10.]]
+
+        assert allclose(f, answer)
+
+    def test_fit_to_mesh_w_georef(self):
+        """Simple check that georef works at the fit_to_mesh level
+        """
+        
+        from coordinate_transforms.geo_reference import Geo_reference
+
+        #Mesh
+        vertex_coordinates = [[0.76, 0.76],
+                              [0.76, 5.76],
+                              [5.76, 0.76]]
+        triangles = [[0,2,1]]
+
+        mesh_geo = Geo_reference(56,-0.76,-0.76)
+        #print "mesh_geo.get_absolute(vertex_coordinates)", \
+         #     mesh_geo.get_absolute(vertex_coordinates)
+
+        #Data                       
+        data_points = [[ 201.0, 401.0],
+                       [ 201.0, 403.0],
+                       [ 203.0, 401.0]]
+
+        z = [2, 4, 4]
+        
+        data_geo = Geo_reference(56,-200,-400)
+
+        #print "data_geo.get_absolute(data_points)", \
+        #      data_geo.get_absolute(data_points)
+        
+        #Fit
+        zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
+                         data_origin = data_geo.get_origin(),
+                         mesh_origin = mesh_geo.get_origin(),
+                         alpha = 0)
+        assert allclose( zz, [0,5,5] )
+
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Least_Squares,'test')
-    #suite = unittest.makeSuite(Test_Least_Squares,'test_smoothing_and_interpolation')
-    #suite = unittest.makeSuite(Test_Least_Squares,'test_smooth_attributes_to_mesh_one_point')
+    suite = unittest.makeSuite(Test_Fit,'test')
+    #suite = unittest.makeSuite(Test_Fit,'test_smoothing_and_interpolation')
+    #suite = unittest.makeSuite(Test_Fit,'test_smooth_attributes_to_mesh_one_point')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)