Changeset 2897


Ignore:
Timestamp:
May 17, 2006, 4:13:37 PM (18 years ago)
Author:
duncan
Message:

close to finishing fit, which replaces least squares

Location:
inundation/fit_interpolate
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • inundation/fit_interpolate/fit.py

    r2802 r2897  
    1616   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
    1717   Geoscience Australia, 2004.
     18
     19   TO DO
     20   * test geo_ref, geo_spatial
    1821"""
    1922
    20 from geospatial_data.geospatial_data import Geospatial_data
     23from Numeric import zeros, Float
     24
     25from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
    2126from fit_interpolate.general_fit_interpolate import FitInterpolate
     27from utilities.sparse import Sparse, Sparse_CSR
     28from utilities.polygon import in_and_outside_polygon
     29from fit_interpolate.search_functions import search_tree_of_vertices
     30from utilities.cg_solve import conjugate_gradient
     31from utilities.numerical_tools import ensure_numeric, gradient
     32
     33import exceptions
     34class ToFewPointsError(exceptions.Exception): pass
    2235
    2336DEFAULT_ALPHA = 0.001
     
    6679
    6780        if alpha is None:
     81
    6882            self.alpha = DEFAULT_ALPHA
    6983        else:   
    7084            self.alpha = alpha
    71        
    7285        FitInterpolate.__init__(self,
    7386                 vertex_coordinates,
     
    7992        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (vertices)
    8093       
    81         #Build Atz and AtA matrix
    82         self.AtA = Sparse(m,m)
    83         self.Atz = zeros((m,att_num), Float)
    84      
    85 
     94        self.AtA = None
     95        self.Atz = None
     96
     97        self.point_count = 0
     98        if self.alpha <> 0:
     99            if verbose: print 'Building smoothing matrix'
     100            self._build_smoothing_matrix_D()
     101           
    86102    def _build_coefficient_matrix_B(self,
    87103                                  verbose = False):
    88         """Build final coefficient matrix"""
    89 
     104        """
     105        Build final coefficient matrix
     106
     107        Precon
     108        If alpha is not zero, matrix D has been built
     109        Matrix Ata has been built
     110        """
     111
     112        if self.alpha <> 0:
     113            #if verbose: print 'Building smoothing matrix'
     114            #self._build_smoothing_matrix_D()
     115            self.B = self.AtA + self.alpha*self.D
     116        else:
     117            self.B = self.AtA
     118
     119        #Convert self.B matrix to CSR format for faster matrix vector
     120        self.B = Sparse_CSR(self.B)
    90121
    91122    def _build_smoothing_matrix_D(self):
     
    114145        are obtained by computing the gradient a_k, b_k for basis function k
    115146        """
     147       
     148        #FIXME: algorithm might be optimised by computing local 9x9
     149        #"element stiffness matrices:
     150
     151        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
     152
     153        self.D = Sparse(m,m)
     154
     155        #For each triangle compute contributions to D = D1+D2
     156        for i in range(len(self.mesh)):
     157
     158            #Get area
     159            area = self.mesh.areas[i]
     160
     161            #Get global vertex indices
     162            v0 = self.mesh.triangles[i,0]
     163            v1 = self.mesh.triangles[i,1]
     164            v2 = self.mesh.triangles[i,2]
     165
     166            #Get the three vertex_points
     167            xi0 = self.mesh.get_vertex_coordinate(i, 0)
     168            xi1 = self.mesh.get_vertex_coordinate(i, 1)
     169            xi2 = self.mesh.get_vertex_coordinate(i, 2)
     170
     171            #Compute gradients for each vertex
     172            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
     173                              1, 0, 0)
     174
     175            a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
     176                              0, 1, 0)
     177
     178            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
     179                              0, 0, 1)
     180
     181            #Compute diagonal contributions
     182            self.D[v0,v0] += (a0*a0 + b0*b0)*area
     183            self.D[v1,v1] += (a1*a1 + b1*b1)*area
     184            self.D[v2,v2] += (a2*a2 + b2*b2)*area
     185
     186            #Compute contributions for basis functions sharing edges
     187            e01 = (a0*a1 + b0*b1)*area
     188            self.D[v0,v1] += e01
     189            self.D[v1,v0] += e01
     190
     191            e12 = (a1*a2 + b1*b2)*area
     192            self.D[v1,v2] += e12
     193            self.D[v2,v1] += e12
     194
     195            e20 = (a2*a0 + b2*b0)*area
     196            self.D[v2,v0] += e20
     197            self.D[v0,v2] += e20
     198
     199
     200    def get_D(self):
     201        return self.D.todense()
     202
    116203
    117204    def _build_matrix_AtA_Atz(self,
     
    131218        z and points are numeric
    132219        Point_coordindates and mesh vertices have the same origin.
    133         """
    134        
     220
     221        The number of attributes of the data points does not change
     222        """
     223        #Build n x m interpolation matrix
     224
     225        if self.AtA == None:
     226            # AtA and Atz need ot be initialised.
     227            m = self.mesh.coordinates.shape[0] #Nbr of vertices
     228            if len(z.shape) > 1:
     229                att_num = z.shape[1]
     230                self.Atz = zeros((m,att_num), Float)
     231            else:
     232                att_num = 1
     233                self.Atz = zeros((m,), Float)
     234            assert z.shape[0] == point_coordinates.shape[0]
     235
     236            self.AtA = Sparse(m,m)
     237            #self.Atz = zeros((m,att_num), Float)
     238        self.point_count += point_coordinates.shape[0]
     239        #print "_build_matrix_AtA_Atz - self.point_count", self.point_count
    135240        if verbose: print 'Getting indices inside mesh boundary'
    136         #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon() 
     241        #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon()
    137242        self.inside_poly_indices, self.outside_poly_indices  = \
    138243                     in_and_outside_polygon(point_coordinates,
     
    143248
    144249       
     250        n = len(self.inside_poly_indices)
    145251        #Compute matrix elements for points inside the mesh
    146252        for i in self.inside_poly_indices:
     
    162268
    163269                for j in js:
    164                     self.Atz[j] +=  sigmas[j]*z[i]                 
     270                    self.Atz[j] +=  sigmas[j]*z[i]
     271                    #print "self.Atz building", self.Atz
     272                    #print "self.Atz[j]", self.Atz[j]
     273                    #print " sigmas[j]", sigmas[j]
     274                    #print "z[i]",z[i]
     275                    #print "result", sigmas[j]*z[i]
     276                   
    165277                    for k in js:
    166                         if interp_only == False:
    167                             self.AtA[j,k] += sigmas[j]*sigmas[k]
     278                        self.AtA[j,k] += sigmas[j]*sigmas[k]
    168279            else:
    169280                msg = 'Could not find triangle for point', x
     
    171282   
    172283       
    173     def fit(self, point_coordinates=point_coordinates, z=z):
     284    def fit(self, point_coordinates=None, z=None,
     285                              verbose = False,
     286                              point_origin = None):
    174287        """Fit a smooth surface to given 1d array of data points z.
    175288
     
    184297         
    185298        """
    186         # build ata and atz
    187         # solve fit
    188 
    189        
    190     def build_fit_subset(self, point_coordinates, z):
     299        if point_coordinates is None:
     300            assert self.AtA <> None
     301            assert self.Atz <> None
     302            #FIXME (DSG) - do  a message
     303        else:
     304            point_coordinates = ensure_absolute(point_coordinates,
     305                                                geo_reference=point_origin)
     306            self.build_fit_subset(point_coordinates, z, verbose)
     307       
     308        #Check sanity
     309        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
     310        n = self.point_count
     311        if n<m and self.alpha == 0.0:
     312            msg = 'ERROR (least_squares): Too few data points\n'
     313            msg += 'There are only %d data points and alpha == 0. ' %n
     314            msg += 'Need at least %d\n' %m
     315            msg += 'Alternatively, set smoothing parameter alpha to a small '
     316            msg += 'positive value,\ne.g. 1.0e-3.'
     317            raise ToFewPointsError(msg)
     318
     319        self._build_coefficient_matrix_B(verbose)
     320        return conjugate_gradient(self.B, self.Atz, self.Atz,
     321                                  imax=2*len(self.Atz) )
     322
     323       
     324    def build_fit_subset(self, point_coordinates, z,
     325                              verbose = False):
    191326        """Fit a smooth surface to given 1d array of data points z.
    192327
     
    202337        """
    203338        #Note: Don't get the z info from Geospatial_data.attributes yet.
    204         # That means fit has to handle attribute title info.
     339        # If we did fit would have to handle attribute title info.
    205340
    206341        #FIXME(DSG-DSG): Check that the vert and point coords
    207342        #have the same zone.
    208343        if isinstance(point_coordinates,Geospatial_data):
    209             point_coordinates = vertex_coordinates.get_data_points( \
     344            point_coordinates = point_coordinates.get_data_points( \
    210345                absolute = True)
    211346       
     
    214349        point_coordinates = ensure_numeric(point_coordinates, Float)
    215350
     351        self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
    216352
    217353
     
    225361                verbose = False,
    226362                acceptable_overshoot = 1.01,
    227                 expand_search = False,
     363                mesh_origin = None,
    228364                data_origin = None,
    229                 mesh_origin = None,
    230                 precrop = False,
    231365                use_cache = False):
    232366    """
     
    236370
    237371        Inputs:
    238 
    239           vertex_coordinates: List of coordinate pairs [xi, eta] of points
    240           constituting mesh (or a an m x 2 Numeric array)
     372        vertex_coordinates: List of coordinate pairs [xi, eta] of
     373              points constituting a mesh (or an m x 2 Numeric array or
     374              a geospatial object)
     375              Points may appear multiple times
     376              (e.g. if vertices have discontinuities)
    241377
    242378          triangles: List of 3-tuples (or a Numeric array) of
     
    252388          as min(z) - acceptable_overshoot*delta z and upper limit
    253389          as max(z) + acceptable_overshoot*delta z
     390
     391          mesh_origin: A geo_reference object or 3-tuples consisting of
     392              UTM zone, easting and northing.
     393              If specified vertex coordinates are assumed to be
     394              relative to their respective origins.
    254395         
    255396
    256           point_attributes: Vector or array of data at the point_coordinates.
    257 
    258           data_origin and mesh_origin are 3-tuples consisting of
    259           UTM zone, easting and northing. If specified
    260           point coordinates and vertex coordinates are assumed to be
    261           relative to their respective origins.
     397          point_attributes: Vector or array of data at the
     398                            point_coordinates.
    262399
    263400    """
    264 
     401    #Since this is a wrapper for fit, lets handle the geo_spatial att's
     402   
    265403    if use_cache is True:
    266404        interp = cache(_fit,
     
    268406                        triangles),
    269407                       {'verbose': verbose,
    270                         'mesh_origin': mesh_origin},
     408                        'mesh_origin': mesh_origin,
     409                        'alpha':alpha},
    271410                       verbose = verbose)       
    272411       
    273412    else:
    274         interp = Interpolation(vertex_coordinates,
    275                                triangles,
    276                                verbose = verbose,
    277                                mesh_origin = mesh_origin)
    278        
    279     vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
    280 
    281        
    282 #                               point_coordinates,
    283 #                               data_origin = data_origin,
    284 
    285 
    286     #Sanity check
    287     point_coordinates = ensure_numeric(point_coordinates)
    288     vertex_coordinates = ensure_numeric(vertex_coordinates)
    289 
    290     #Data points
    291     X = point_coordinates[:,0]
    292     Y = point_coordinates[:,1] 
    293     Z = ensure_numeric(point_attributes)
    294     if len(Z.shape) == 1:
    295         Z = Z[:, NewAxis]
    296        
    297 
    298     #Data points inside mesh boundary
    299     indices = interp.point_indices
    300     if indices is not None:   
    301         Xc = take(X, indices)
    302         Yc = take(Y, indices)   
    303         Zc = take(Z, indices)
    304     else:
    305         Xc = X
    306         Yc = Y 
    307         Zc = Z       
     413        interp = Fit(vertex_coordinates,
     414                     triangles,
     415                     verbose = verbose,
     416                     mesh_origin = mesh_origin,
     417                     alpha=alpha)
     418       
     419    vertex_attributes = interp.fit(point_coordinates,
     420                                   point_attributes,
     421                                   point_origin = data_origin,
     422                                   verbose = verbose)
     423
     424       
     425    # Add the value checking stuff that's in least squares.
     426    # Maybe this stuff should get pushed down into Fit.
     427    # at least be a method of Fit.
     428    # Or intigrate it into the fit method, saving teh max and min's
     429    # as att's.
    308430   
    309     #Vertex coordinates
    310     Xi = vertex_coordinates[:,0]
    311     Eta = vertex_coordinates[:,1]       
    312     Zeta = ensure_numeric(vertex_attributes)
    313     if len(Zeta.shape) == 1:
    314         Zeta = Zeta[:, NewAxis]   
    315 
    316     for i in range(Zeta.shape[1]): #For each attribute
    317         zeta = Zeta[:,i]
    318         z = Z[:,i]               
    319         zc = Zc[:,i]
    320 
    321         max_zc = max(zc)
    322         min_zc = min(zc)
    323         delta_zc = max_zc-min_zc
    324         upper_limit = max_zc + delta_zc*acceptable_overshoot
    325         lower_limit = min_zc - delta_zc*acceptable_overshoot       
    326        
    327 
    328         if max(zeta) > upper_limit or min(zeta) < lower_limit:
    329             msg = 'Least sqares produced values outside the allowed '
    330             msg += 'range [%f, %f].\n' %(lower_limit, upper_limit)
    331             msg += 'z in [%f, %f], zeta in [%f, %f].\n' %(min_zc, max_zc,
    332                                                           min(zeta), max(zeta))
    333             msg += 'If greater range is needed, increase the value of '
    334             msg += 'acceptable_fit_overshoot (currently %.2f).\n' %(acceptable_overshoot)
    335 
    336 
    337             offending_vertices = (zeta > upper_limit or zeta < lower_limit)
    338             Xi_c = compress(offending_vertices, Xi)
    339             Eta_c = compress(offending_vertices, Eta)
    340             offending_coordinates = concatenate((Xi_c[:, NewAxis],
    341                                                  Eta_c[:, NewAxis]),
    342                                                 axis=1)
    343 
    344             msg += 'Offending locations:\n %s' %(offending_coordinates)
    345            
    346             raise FittingError, msg
    347 
    348 
    349    
    350         if verbose:
    351             print '+------------------------------------------------'
    352             print 'Least squares statistics'
    353             print '+------------------------------------------------'   
    354             print 'points: %d points' %(len(z))
    355             print '    x in [%f, %f]'%(min(X), max(X))
    356             print '    y in [%f, %f]'%(min(Y), max(Y))
    357             print '    z in [%f, %f]'%(min(z), max(z))
    358             print
    359 
    360             if indices is not None:
    361                 print 'Cropped points: %d points' %(len(zc))
    362                 print '    x in [%f, %f]'%(min(Xc), max(Xc))
    363                 print '    y in [%f, %f]'%(min(Yc), max(Yc))
    364                 print '    z in [%f, %f]'%(min(zc), max(zc))
    365                 print
    366            
    367 
    368             print 'Mesh: %d vertices' %(len(zeta))
    369             print '    xi in [%f, %f]'%(min(Xi), max(Xi))
    370             print '    eta in [%f, %f]'%(min(Eta), max(Eta))
    371             print '    zeta in [%f, %f]'%(min(zeta), max(zeta))
    372             print '+------------------------------------------------'
    373 
    374431    return vertex_attributes
    375432
  • inundation/fit_interpolate/general_fit_interpolate.py

    r2802 r2897  
    3232from utilities.numerical_tools import ensure_numeric, mean, INF
    3333from utilities.polygon import in_and_outside_polygon
    34 from geospatial_data.geospatial_data import Geospatial_data
     34from geospatial_data.geospatial_data import Geospatial_data, \
     35     ensure_absolute
    3536from search_functions import search_tree_of_vertices
    3637
     
    7576        #Convert input to Numeric arrays
    7677        triangles = ensure_numeric(triangles, Int)
    77        
    78         if isinstance(vertex_coordinates,Geospatial_data):
    79             vertex_coordinates = vertex_coordinates.get_data_points( \
    80                 absolute = True)
    81             msg = "Use a Geospatial_data object or a mesh origin. Not both."
    82             assert mesh_origin == None, msg
    83            
    84         else:
    85             vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
    86         #Build underlying mesh
    87         if verbose: print 'Building mesh'
    88        
    89         if mesh_origin is None:
    90             geo = None #Geo_reference()
    91         else:
    92             if isinstance(mesh_origin, Geo_reference):
    93                 geo = mesh_origin
    94             else:
    95                 geo = Geo_reference(mesh_origin[0],
    96                                     mesh_origin[1],
    97                                     mesh_origin[2])
    98             vertex_coordinates = geo.get_absolute(vertex_coordinates)
     78        vertex_coordinates = ensure_absolute(vertex_coordinates,
     79                                         geo_reference = mesh_origin)
     80
    9981        #Don't pass geo_reference to mesh.  It doesn't work.
    10082        self.mesh = Mesh(vertex_coordinates, triangles)
  • inundation/fit_interpolate/run_long_benchmark.py

    r2563 r2897  
    88
    99use_least_squares_list = [True,False]
    10 is_fit_list = [False]
     10is_fit_list = [True]
    1111num_of_points_list = [10]
    1212maxArea_list = [0.1, 0.001]
  • inundation/fit_interpolate/spike_least_squares.py

    r2884 r2897  
    565565        self.AtA = Sparse(m,m)
    566566        self.Atz = zeros((m,att_num), Float)
     567        print "self.Atz",self.Atz
     568        self.Atz[0] = 100
     569        print "self.Atz",self.Atz
    567570
    568571        #Build quad tree of vertices
  • inundation/fit_interpolate/spike_test_least_squares.py

    r2781 r2897  
    1212from utilities.sparse import Sparse, Sparse_CSR
    1313
    14 from coordinate_transforms.geo_reference import Geo_reference
    15 
    16 def distance(x, y):
    17     return sqrt( sum( (array(x)-array(y))**2 ))
     14
    1815
    1916def linear_function(point):
  • inundation/fit_interpolate/test_fit.py

    r2802 r2897  
    66from math import sqrt
    77
    8 
    9 from spike_least_squares import *
    10 from Numeric import allclose, array, transpose
    11 from Numeric import zeros, take, compress, array, Float, Int, dot, transpose, concatenate, ArrayType
     8from Numeric import zeros, take, compress, Float, Int, dot, concatenate, \
     9     ArrayType, allclose, array
     10
     11from fit import *
    1212from utilities.sparse import Sparse, Sparse_CSR
    13 
    1413from coordinate_transforms.geo_reference import Geo_reference
     14from utilities.numerical_tools import ensure_numeric
     15from geospatial_data.geospatial_data import Geospatial_data
    1516
    1617def distance(x, y):
     
    2223
    2324
    24 class Test_Least_Squares(unittest.TestCase):
     25class Test_Fit(unittest.TestCase):
    2526
    2627    def setUp(self):
     
    2930    def tearDown(self):
    3031        pass
     32
    3133
    3234    def test_smooth_attributes_to_mesh(self):
     
    4648
    4749        z = [z1, z2, z3]
    48         interp = Interpolation(points, triangles, z,data_coords, alpha=0)
     50        fit = Fit(points, triangles, alpha=0)
    4951        #print "interp.get_A()", interp.get_A()
    50         A = interp.A
    51         #print "A",A
     52        fit._build_matrix_AtA_Atz(ensure_numeric(data_coords),
     53                                  ensure_numeric(z))
     54        #print "Atz - from fit", fit.Atz
     55        #print "AtA - from fit", fit.AtA.todense()
    5256        #print "z",z
    53         Atz = A.trans_mult(z)
    54         #print "Atz",Atz
    55        
    56         f = interp.fit(z)
     57
     58        assert allclose(fit.Atz, [2.8, 3.6, 3.6], atol=1e-7)
     59
     60        f = fit.fit()
     61       
    5762        answer = [0, 5., 5.]
    5863
     
    6267        assert allclose(f, answer, atol=1e-7)
    6368
    64        
     69    def test_smooth_att_to_meshII(self):
     70
     71        a = [0.0, 0.0]
     72        b = [0.0, 5.0]
     73        c = [5.0, 0.0]
     74        points = [a, b, c]
     75        triangles = [ [1,0,2] ]   #bac
     76
     77        d1 = [1.0, 1.0]
     78        d2 = [1.0, 2.0]
     79        d3 = [3.0,1.0]
     80        data_coords = [d1, d2, d3]
     81        z = linear_function(data_coords)
     82        #print "z",z
     83
     84        interp = Fit(points, triangles, alpha=0.0)
     85        f = interp.fit(data_coords, z)
     86        answer = linear_function(points)
     87        #print "f\n",f
     88        #print "answer\n",answer
     89
     90        assert allclose(f, answer)
     91
     92    def test_smooth_attributes_to_meshIII(self):
     93
     94        a = [-1.0, 0.0]
     95        b = [3.0, 4.0]
     96        c = [4.0,1.0]
     97        d = [-3.0, 2.0] #3
     98        e = [-1.0,-2.0]
     99        f = [1.0, -2.0] #5
     100
     101        vertices = [a, b, c, d,e,f]
     102        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
     103
     104        point_coords = [[-2.0, 2.0],
     105                        [-1.0, 1.0],
     106                        [0.0,2.0],
     107                        [1.0, 1.0],
     108                        [2.0, 1.0],
     109                        [0.0,0.0],
     110                        [1.0, 0.0],
     111                        [0.0, -1.0],
     112                        [-0.2,-0.5],
     113                        [-0.9, -1.5],
     114                        [0.5, -1.9],
     115                        [3.0,1.0]]
     116
     117        z = linear_function(point_coords)
     118        interp = Fit(vertices, triangles,
     119                                alpha=0.0)
     120
     121        #print 'z',z
     122        f = interp.fit(point_coords,z)
     123        answer = linear_function(vertices)
     124        #print "f\n",f
     125        #print "answer\n",answer
     126        assert allclose(f, answer)
     127
     128   
     129    def test_smooth_attributes_to_meshIV(self):
     130        """ Testing 2 attributes smoothed to the mesh
     131        """
     132
     133        a = [0.0, 0.0]
     134        b = [0.0, 5.0]
     135        c = [5.0, 0.0]
     136        points = [a, b, c]
     137        triangles = [ [1,0,2] ]   #bac
     138
     139        d1 = [1.0, 1.0]
     140        d2 = [1.0, 3.0]
     141        d3 = [3.0, 1.0]
     142        z1 = [2, 4]
     143        z2 = [4, 8]
     144        z3 = [4, 8]
     145        data_coords = [d1, d2, d3]
     146
     147        z = [z1, z2, z3]
     148        fit = Fit(points, triangles, alpha=0)
     149       
     150        f =  fit.fit(data_coords,z)
     151        answer = [[0,0], [5., 10.], [5., 10.]]
     152        assert allclose(f, answer)
     153
     154    def test_smooth_attributes_to_mesh_build_fit_subset(self):
     155
     156        a = [-1.0, 0.0]
     157        b = [3.0, 4.0]
     158        c = [4.0,1.0]
     159        d = [-3.0, 2.0] #3
     160        e = [-1.0,-2.0]
     161        f = [1.0, -2.0] #5
     162
     163        vertices = [a, b, c, d,e,f]
     164        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
     165
     166        interp = Fit(vertices, triangles,
     167                                alpha=0.0)
     168
     169        point_coords = [[-2.0, 2.0],
     170                        [-1.0, 1.0],
     171                        [0.0,2.0],
     172                        [1.0, 1.0],
     173                       ]
     174
     175        z = linear_function(point_coords)
     176
     177        f = interp.build_fit_subset(point_coords,z)
     178       
     179        point_coords = [
     180                        [2.0, 1.0],
     181                        [0.0,0.0],
     182                        [1.0, 0.0],
     183                        [0.0, -1.0],
     184                        [-0.2,-0.5],
     185                        [-0.9, -1.5],
     186                        [0.5, -1.9],
     187                        [3.0,1.0]]
     188       
     189        z = linear_function(point_coords)
     190
     191        f = interp.build_fit_subset(point_coords,z)
     192       
     193        #print 'z',z
     194        f = interp.fit()
     195        answer = linear_function(vertices)
     196        #print "f\n",f
     197        #print "answer\n",answer
     198        assert allclose(f, answer)
     199
     200    def test_fit_and_interpolation(self):
     201        from mesh import Mesh
     202
     203        a = [0.0, 0.0]
     204        b = [0.0, 2.0]
     205        c = [2.0, 0.0]
     206        d = [0.0, 4.0]
     207        e = [2.0, 2.0]
     208        f = [4.0, 0.0]
     209
     210        points = [a, b, c, d, e, f]
     211        #bac, bce, ecf, dbe, daf, dae
     212        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     213
     214        #Get (enough) datapoints
     215        data_points = [[ 0.66666667, 0.66666667],
     216                       [ 1.33333333, 1.33333333],
     217                       [ 2.66666667, 0.66666667],
     218                       [ 0.66666667, 2.66666667],
     219                       [ 0.0, 1.0],
     220                       [ 0.0, 3.0],
     221                       [ 1.0, 0.0],
     222                       [ 1.0, 1.0],
     223                       [ 1.0, 2.0],
     224                       [ 1.0, 3.0],
     225                       [ 2.0, 1.0],
     226                       [ 3.0, 0.0],
     227                       [ 3.0, 1.0]]
     228
     229        z = linear_function(data_points)
     230        interp = Fit(points, triangles,
     231                                alpha=0.0)
     232
     233        answer = linear_function(points)
     234
     235        f = interp.fit(data_points, z)
     236
     237        #print "f",f
     238        #print "answer",answer
     239        assert allclose(f, answer)
     240
     241       
     242    def test_smoothing_and_interpolation(self):
     243
     244        a = [0.0, 0.0]
     245        b = [0.0, 2.0]
     246        c = [2.0, 0.0]
     247        d = [0.0, 4.0]
     248        e = [2.0, 2.0]
     249        f = [4.0, 0.0]
     250
     251        points = [a, b, c, d, e, f]
     252        #bac, bce, ecf, dbe, daf, dae
     253        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     254
     255        #Get (too few!) datapoints
     256        data_points = [[ 0.66666667, 0.66666667],
     257                       [ 1.33333333, 1.33333333],
     258                       [ 2.66666667, 0.66666667],
     259                       [ 0.66666667, 2.66666667]]
     260
     261        z = linear_function(data_points)
     262        answer = linear_function(points)
     263
     264        #Make interpolator with too few data points and no smoothing
     265        interp = Fit(points, triangles, alpha=0.0)
     266        #Must raise an exception
     267        try:
     268            f = interp.fit(data_points,z)
     269        except ToFewPointsError:
     270            pass
     271
     272        #Now try with smoothing parameter
     273        interp = Fit(points, triangles, alpha=1.0e-13)
     274
     275        f = interp.fit(data_points,z)
     276        #f will be different from answer due to smoothing
     277        assert allclose(f, answer,atol=5)
     278
     279
     280    #Tests of smoothing matrix
     281    def test_smoothing_matrix_one_triangle(self):
     282        from Numeric import dot
     283        a = [0.0, 0.0]
     284        b = [0.0, 2.0]
     285        c = [2.0,0.0]
     286        points = [a, b, c]
     287
     288        vertices = [ [1,0,2] ]   #bac
     289
     290        interp = Fit(points, vertices)
     291
     292        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
     293                                   [-0.5, 0.5, 0],
     294                                   [-0.5, 0, 0.5]])
     295
     296        #Define f(x,y) = x
     297        f = array([0,0,2]) #Value at global vertex 2
     298
     299        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
     300        #           int 1 dx dy = area = 2
     301        assert dot(dot(f, interp.get_D()), f) == 2
     302
     303        #Define f(x,y) = y
     304        f = array([0,2,0])  #Value at global vertex 1
     305
     306        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
     307        #           int 1 dx dy = area = 2
     308        assert dot(dot(f, interp.get_D()), f) == 2
     309
     310        #Define f(x,y) = x+y
     311        f = array([0,2,2])  #Values at global vertex 1 and 2
     312
     313        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
     314        #           int 2 dx dy = 2*area = 4
     315        assert dot(dot(f, interp.get_D()), f) == 4
     316
     317
     318    def test_smoothing_matrix_more_triangles(self):
     319        from Numeric import dot
     320
     321        a = [0.0, 0.0]
     322        b = [0.0, 2.0]
     323        c = [2.0,0.0]
     324        d = [0.0, 4.0]
     325        e = [2.0, 2.0]
     326        f = [4.0,0.0]
     327
     328        points = [a, b, c, d, e, f]
     329        #bac, bce, ecf, dbe, daf, dae
     330        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     331
     332        interp = Fit(points, vertices)
     333
     334
     335        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
     336        #                           [-0.5, 0.5, 0],
     337        #                           [-0.5, 0, 0.5]])
     338
     339        #Define f(x,y) = x
     340        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
     341
     342        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
     343        #           int 1 dx dy = total area = 8
     344        assert dot(dot(f, interp.get_D()), f) == 8
     345
     346        #Define f(x,y) = y
     347        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
     348
     349        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
     350        #           int 1 dx dy = area = 8
     351        assert dot(dot(f, interp.get_D()), f) == 8
     352
     353        #Define f(x,y) = x+y
     354        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
     355
     356        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
     357        #           int 2 dx dy = 2*area = 16
     358        assert dot(dot(f, interp.get_D()), f) == 16
     359
     360
     361
     362    def test_fit_and_interpolation_with_different_origins(self):
     363        """Fit a surface to one set of points. Then interpolate that surface
     364        using another set of points.
     365        This test tests situtaion where points and mesh belong to a different
     366        coordinate system as defined by origin.
     367        """
     368        from mesh import Mesh
     369
     370        #Setup mesh used to represent fitted function
     371        a = [0.0, 0.0]
     372        b = [0.0, 2.0]
     373        c = [2.0, 0.0]
     374        d = [0.0, 4.0]
     375        e = [2.0, 2.0]
     376        f = [4.0, 0.0]
     377
     378        points = [a, b, c, d, e, f]
     379        #bac, bce, ecf, dbe, daf, dae
     380        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     381
     382        #Datapoints to fit from
     383        data_points1 = [[ 0.66666667, 0.66666667],
     384                        [ 1.33333333, 1.33333333],
     385                        [ 2.66666667, 0.66666667],
     386                        [ 0.66666667, 2.66666667],
     387                        [ 0.0, 1.0],
     388                        [ 0.0, 3.0],
     389                        [ 1.0, 0.0],
     390                        [ 1.0, 1.0],
     391                        [ 1.0, 2.0],
     392                        [ 1.0, 3.0],
     393                        [ 2.0, 1.0],
     394                        [ 3.0, 0.0],
     395                        [ 3.0, 1.0]]
     396
     397
     398        #First check that things are OK when using same origin
     399        mesh_origin = (56, 290000, 618000) #zone, easting, northing
     400        data_origin = (56, 290000, 618000) #zone, easting, northing
     401
     402
     403        #Fit surface to mesh
     404        interp = Fit(points, triangles,
     405                               alpha=0.0,
     406                               mesh_origin = mesh_origin)
     407
     408        data_geo_spatial = Geospatial_data(data_points1,
     409                         geo_reference = Geo_reference(56, 290000, 618000))
     410        z = linear_function(data_points1) #Example z-values
     411        f = interp.fit(data_geo_spatial, z)   #Fitted values at vertices
     412
     413        #Shift datapoints according to new origins
     414        for k in range(len(data_points1)):
     415            data_points1[k][0] += mesh_origin[1] - data_origin[1]
     416            data_points1[k][1] += mesh_origin[2] - data_origin[2]
     417
     418
     419
     420        #Fit surface to mesh
     421        interp = Fit(points, triangles,
     422                               alpha=0.0) #,
     423                               # mesh_origin = mesh_origin)
     424
     425        f1 = interp.fit(data_points1,z) #Fitted values at vertices (using same z as before)
     426
     427        assert allclose(f,f1), 'Fit should have been unaltered'
     428
     429
     430    def test_smooth_attributes_to_mesh_function(self):
     431        """ Testing 2 attributes smoothed to the mesh
     432        """
     433
     434        a = [0.0, 0.0]
     435        b = [0.0, 5.0]
     436        c = [5.0, 0.0]
     437        points = [a, b, c]
     438        triangles = [ [1,0,2] ]   #bac
     439
     440        d1 = [1.0, 1.0]
     441        d2 = [1.0, 3.0]
     442        d3 = [3.0, 1.0]
     443        z1 = [2, 4]
     444        z2 = [4, 8]
     445        z3 = [4, 8]
     446        data_coords = [d1, d2, d3]
     447        z = [z1, z2, z3]
     448
     449        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
     450        answer = [[0, 0], [5., 10.], [5., 10.]]
     451
     452        assert allclose(f, answer)
     453
     454    def test_fit_to_mesh_w_georef(self):
     455        """Simple check that georef works at the fit_to_mesh level
     456        """
     457       
     458        from coordinate_transforms.geo_reference import Geo_reference
     459
     460        #Mesh
     461        vertex_coordinates = [[0.76, 0.76],
     462                              [0.76, 5.76],
     463                              [5.76, 0.76]]
     464        triangles = [[0,2,1]]
     465
     466        mesh_geo = Geo_reference(56,-0.76,-0.76)
     467        #print "mesh_geo.get_absolute(vertex_coordinates)", \
     468         #     mesh_geo.get_absolute(vertex_coordinates)
     469
     470        #Data                       
     471        data_points = [[ 201.0, 401.0],
     472                       [ 201.0, 403.0],
     473                       [ 203.0, 401.0]]
     474
     475        z = [2, 4, 4]
     476       
     477        data_geo = Geo_reference(56,-200,-400)
     478
     479        #print "data_geo.get_absolute(data_points)", \
     480        #      data_geo.get_absolute(data_points)
     481       
     482        #Fit
     483        zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
     484                         data_origin = data_geo.get_origin(),
     485                         mesh_origin = mesh_geo.get_origin(),
     486                         alpha = 0)
     487        assert allclose( zz, [0,5,5] )
     488
     489
     490
    65491#-------------------------------------------------------------
    66492if __name__ == "__main__":
    67     suite = unittest.makeSuite(Test_Least_Squares,'test')
    68     #suite = unittest.makeSuite(Test_Least_Squares,'test_smoothing_and_interpolation')
    69     #suite = unittest.makeSuite(Test_Least_Squares,'test_smooth_attributes_to_mesh_one_point')
     493    suite = unittest.makeSuite(Test_Fit,'test')
     494    #suite = unittest.makeSuite(Test_Fit,'test_smoothing_and_interpolation')
     495    #suite = unittest.makeSuite(Test_Fit,'test_smooth_attributes_to_mesh_one_point')
    70496    runner = unittest.TextTestRunner(verbosity=1)
    71497    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.