Changeset 485 for inundation/ga


Ignore:
Timestamp:
Nov 5, 2004, 10:51:54 AM (20 years ago)
Author:
ole
Message:

Cleaned up is sparse and least squares.
Verified that quad trees work.
Implemented sparse matrix x matrix mult and simplified
interpolate in least_squares.
Added more unit testing.

Location:
inundation/ga/storm_surge/pyvolution
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r484 r485  
    1919
    2020
     21#FIXME (Ole): Currently datapoints outside the triangular mesh are ignored.
     22#             Is there a clean way of inlcuding them?
     23
     24
    2125import exceptions
    2226class ShapeError(exceptions.Exception): pass
    … …  
    3034try:
    3135    from util import gradient
    32 
    3336except ImportError, e:
    3437    #FIXME reduce the dependency of modules in pyvolution
    … …  
    5053DEFAULT_ALPHA = 0.001
    5154   
    52 
    5355def fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha=DEFAULT_ALPHA):
    5456    """
    … …  
    6769    old_title_list = mesh_dict['generatedpointattributetitlelist']
    6870
    69     
     71   
    7072    # load in the .xya file
    71    
    7273    point_dict = load_xya_file(point_file)
    7374    point_coordinates = point_dict['pointlist']
    7475    point_attributes = point_dict['pointattributelist']
    7576    title_string = point_dict['title']
    76     title_list = title_string.split(',') #iffy! Hard coding title delimiter 
     77    title_list = title_string.split(',') #FIXME iffy! Hard coding title delimiter 
    7778    for i in range(len(title_list)):
    7879        title_list[i] = title_list[i].strip()
    … …  
    218219        self.AtA = Sparse(m,m)
    219220
    220         #Build quad tree of vertices
    221         #print 'Building quad tree'
     221        #Build quad tree of vertices (FIXME: Is this the right spot for that?)
    222222        root = build_quadtree(self.mesh)
    223         #print 'Quad tree built'
    224         #root.show()
    225223
    226224        #Compute matrix elements
    … …  
    235233            candidate_vertices = root.search(x[0], x[1])
    236234
    237             #Find triangle containing x
     235            #Find triangle containing x:
    238236            element_found = False           
     237           
     238            #For all vertices in same cell as point x
    239239            for v in candidate_vertices:
    240                 for triangle_id, vertex_id in self.mesh.vertexlist[v]:
    241 
    242                     k = triangle_id
     240           
     241                #for each triangle id (k) which has v as a vertex
     242                for k, _ in self.mesh.vertexlist[v]:
    243243                   
    244244                    #Get the three vertex_points of candidate triangle
    245245                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
    246246                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
    247                     xi2 = self.mesh.get_vertex_coordinate(k, 2)                
     247                    xi2 = self.mesh.get_vertex_coordinate(k, 2)     
    248248
    249249                    #Get the three normals
    … …  
    269269                    #Don't look for any other triangle
    270270                    break
    271 
     271                   
     272
     273            #Update interpolation matrix A if necessary     
    272274            if element_found is True:       
    273275                #Assign values to matrix A
    … …  
    291293            else:
    292294                pass
    293                 #Ok if there is no triangle for datapoint (as in brute force version)
     295                #Ok if there is no triangle for datapoint
     296                #(as in brute force version)
    294297                #raise 'Could not find triangle for point', x
    295298
    … …  
    301304        m is the number of basis functions phi_k (one per vertex)
    302305
    303 
    304         This is the brute force which is too slow for large problems, but good for testing
     306        This is the brute force which is too slow for large problems,
     307        but could be used for testing
    305308        """
    306309
    … …  
    413416        #"element stiffness matrices:
    414417
    415        
    416418        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
    417419
    418         #self.D = zeros((m,m), Float)
    419420        self.D = Sparse(m,m)
    420421
    … …  
    463464            self.D[v0,v2] += e20             
    464465
    465         #self.D = (self.D).tocsc()
    466466           
    467467    def fit(self, z):
    … …  
    488488        Atz = self.A.trans_mult(z)
    489489
    490 
    491         #print 'fit: Atz',Atz
    492490       
    493491        #Check sanity
    … …  
    500498
    501499
    502         #Solve and return
    503         #return solve_linear_equations(self.get_B(), Atz)
    504 
    505         # caused errors
    506         #return sparse.solve(self.B,Atz)
    507 
    508500        return conjugate_gradient(self.B, Atz, Atz,imax=2*len(Atz) )
    509501        #FIXME: Should we store the result here for later use? (ON)       
    510    
     502
     503           
    511504    def fit_points(self, z):
    512         """
    513         Like fit, but more robust when each point has two or more attributes
    514         """
     505        """Like fit, but more robust when each point has two or more attributes
     506        FIXME(Ole): The name fit_points doesn't carry any meaning for me.
     507        How about something like fit_multiple or fit_columns?
     508        """
     509       
    515510        try:
    516511            return self.fit(z)
    … …  
    525520            n = z.shape[1]               #Number of data points         
    526521
    527             f = zeros((m,n), Float)
    528             #f = sparse.dok_matrix() # even though it wont be sparse?
     522            f = zeros((m,n), Float) #Resulting columns
    529523           
    530524            for i in range(z.shape[1]):
    531525                f[:,i] = self.fit(z[:,i])
     526               
    532527            return f
    533528           
    … …  
    548543        """
    549544
    550         try:
    551             return self.A * f
    552         except ValueError:
    553             # We are here, so it probalby means that f is 2 dimensional
    554             # and so will do each column separately due to problems in
    555             # sparse matrix, 2d array multiplication
    556             (N , M) = self.A.shape
    557             f = array(f).astype(Float)
    558             (m , n) = f.shape
    559             #print "m",m
    560             #print "M",M
    561             if m != M :
    562                 #print "!!self.A",self.A
    563                 #print "!!self.A.todense() ",self.A.todense()
    564                 #print "self.get_A()",self.get_A()
    565                 #return dot(self.get_A(),f)
    566                 raise VectorShapeError, 'Mismatch between A and f dimensions'
    567            
    568             y = zeros( (N,n), Float)
    569 
    570             for i in range(n):
    571                 y[:,i] = self.A * f[:,i]
    572 
    573             #print "!!self.A.todense() ",self.A.todense()
    574             return y
     545        return self.A * f
     546       
     547           
    575548     
    576549    #FIXME: We will need a method 'evaluate(self):' that will interpolate
    … …  
    597570    """
    598571    import os, sys
    599     usage = "usage: %s mesh_input.tsh point.xya  mesh_output.tsh alpha" % os.path.basename(sys.argv[0])
     572    usage = "usage: %s mesh_input.tsh point.xya mesh_output.tsh alpha"\
     573            %os.path.basename(sys.argv[0])
    600574
    601575    if len(sys.argv) < 4:

  • inundation/ga/storm_surge/pyvolution/sparse.py

    r474 r485  
    116116            print 'FIXME: Only Numeric types implemented so far'
    117117
    118 
    119118        #Assume numeric types from now on
    120         R = zeros((self.M,), Float) #Result
    121 
     119       
    122120        if len(B.shape) == 0:
    123121            #Scalar - use __rmul__ method
    124122            R = B*self
     123           
    125124        elif len(B.shape) == 1:
    126125            #Vector
    127126            assert B.shape[0] == self.N, 'Mismatching dimensions'
    128127
     128            R = zeros(self.M, Float) #Result
     129           
    129130            #Multiply nonzero elements
    130131            for key in self.A.keys():
    … …  
    133134                R[i] += self.A[key]*B[j]
    134135        elif len(B.shape) == 2:
    135             raise ValueError, 'Numeric matrix not yet implemented'           
     136       
     137           
     138            R = zeros((self.M, B.shape[1]), Float) #Result matrix
     139
     140            #Multiply nonzero elements
     141            for col in range(R.shape[1]):
     142                #For each column
     143               
     144                for key in self.A.keys():
     145                    i, j = key
     146
     147                    R[i, col] += self.A[key]*B[j, col]
     148           
     149           
    136150        else:
    137151            raise ValueError, 'Dimension too high: d=%d' %len(B.shape)

  • inundation/ga/storm_surge/pyvolution/test_least_squares.py

    r484 r485  
    168168        #print "answer\n",answer
    169169       
    170         assert allclose(f, answer,atol=1e-7)
     170        assert allclose(f, answer, atol=1e-7)
    171171       
    172172       
    … …  
    200200
    201201        vertices = [a, b, c, d,e,f]
    202         triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
     202        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
    203203
    204204        point_coords = [[-2.0, 2.0],
    … …  
    238238        d1 = [1.0, 1.0]
    239239        d2 = [1.0, 3.0]
    240         d3 = [3.0,1.0]
    241         z1 = [2,4]
    242         z2 = [4,8]
    243         z3 = [4,8]
    244         data_coords = [ d1, d2, d3]
     240        d3 = [3.0, 1.0]
     241        z1 = [2, 4]
     242        z2 = [4, 8]
     243        z3 = [4, 8]
     244        data_coords = [d1, d2, d3]
    245245       
    246246        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
    … …  
    260260        d0 = [1.0, 1.0]
    261261        d1 = [1.0, 2.0]
    262         d2 = [3.0,1.0]
     262        d2 = [3.0, 1.0]
    263263        point_coords = [ d0, d1, d2]
    264264       
    … …  
    275275        a = [-1.0, 0.0]
    276276        b = [3.0, 4.0]
    277         c = [4.0,1.0]
     277        c = [4.0, 1.0]
    278278        d = [-3.0, 2.0] #3
    279         e = [-1.0,-2.0]
     279        e = [-1.0, -2.0]
    280280        f = [1.0, -2.0] #5
    281281
    282282        vertices = [a, b, c, d,e,f]
    283         triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
     283        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
    284284
    285285         
    286286        point_coords = [[-2.0, 2.0],
    287287                        [-1.0, 1.0],
    288                         [0.0,2.0],
     288                        [0.0, 2.0],
    289289                        [1.0, 1.0],
    290290                        [2.0, 1.0],
    291                         [0.0,0.0],
     291                        [0.0, 0.0],
    292292                        [1.0, 0.0],
    293293                        [0.0, -1.0],
    294                         [-0.2,-0.5],
     294                        [-0.2, -0.5],
    295295                        [-0.9, -1.5],
    296296                        [0.5, -1.9],
    297                         [3.0,1.0]]
     297                        [3.0, 1.0]]
    298298       
    299299        interp = Interpolation(vertices, triangles, point_coords)
    … …  
    335335        a = [-1.0, 0.0]
    336336        b = [3.0, 4.0]
    337         c = [4.0,1.0]
     337        c = [4.0, 1.0]
    338338        d = [-3.0, 2.0] #3
    339         e = [-1.0,-2.0]
     339        e = [-1.0, -2.0]
    340340        f = [1.0, -2.0] #5
    341341
    342342        vertices = [a, b, c, d,e,f]
    343         triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
     343        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
    344344
    345345         
    346346        point_coords = [[-2.0, 2.0],
    347347                        [-1.0, 1.0],
    348                         [0.0,2.0],
     348                        [0.0, 2.0],
    349349                        [1.0, 1.0],
    350350                        [2.0, 1.0],
    351                         [0.0,0.0],
     351                        [0.0, 0.0],
    352352                        [1.0, 0.0],
    353353                        [0.0, -1.0],
    354                         [-0.2,-0.5],
     354                        [-0.2, -0.5],
    355355                        [-0.9, -1.5],
    356356                        [0.5, -1.9],
    357                         [3.0,1.0]]
     357                        [3.0, 1.0]]
    358358       
    359359        interp = Interpolation(vertices, triangles, point_coords)
    … …  
    384384        d1 = [1.0, 1.0]
    385385        d2 = [1.0, 3.0]
    386         d3 = [3.0,1.0]
    387         z1 = [2,4]
    388         z2 = [4,8]
    389         z3 = [4,8]
    390         data_coords = [ d1, d2, d3]
     386        d3 = [3.0, 1.0]
     387        z1 = [2, 4]
     388        z2 = [4, 8]
     389        z3 = [4, 8]
     390        data_coords = [d1, d2, d3]
    391391        z = [z1, z2, z3]
    392392       

  • inundation/ga/storm_surge/pyvolution/test_sparse.py

    r479 r485  
    33import unittest
    44from math import sqrt
    5 
    65
    76from sparse import *
    … …  
    1514    def tearDown(self):
    1615        pass
    17 
    1816
    1917    def test_init1(Self):
    … …  
    7270        assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])       
    7371
    74     def test_sparse_multiplication(self):
     72    def test_sparse_multiplication_vector(self):
    7573        A = Sparse(3,3)
    7674       
    … …  
    8381        v = [2,3,4]
    8482
    85        
    8683        u = A*v
    8784        assert allclose(u, [6,14,4])
    88 
    89 
    9085
    9186        #Right hand side column
    … …  
    9792        u = A*v[:,1]
    9893        assert allclose(u, [12,16,4])       
     94
     95
     96    def test_sparse_multiplication_matrix(self):
     97        A = Sparse(3,3)
     98       
     99        A[0,0] = 3
     100        A[1,1] = 2
     101        A[1,2] = 2
     102        A[2,2] = 1
     103
     104        #Right hand side matrix
     105        v = array([[2,4],[3,4],[4,4]])
     106
     107        u = A*v
     108        assert allclose(u, [[6,12], [14,16], [4,4]]) 
     109
    99110
    100111
Note: See TracChangeset for help on using the changeset viewer.