Changeset 441


Ignore:
Timestamp:
Oct 24, 2004, 10:08:40 AM (20 years ago)
Author:
steve
Message:

Dealing with 2d data arrays in fit and interpolate

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

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r436 r441  
    2020#FIXME: Current implementation uses full matrices and a general solver.
    2121#Later on we may consider using sparse techniques
     22
     23import exceptions
     24class ShapeError(exceptions.Exception): pass
    2225
    2326from general_mesh import General_Mesh
     
    178181       
    179182    def build_coefficient_matrix_B(self, point_coordinates=None):
    180         """Build final coefficient matrix
    181         """
    182         #self.alpha = 0
     183        """Build final coefficient matrix"""
     184       
     185
    183186        if self.alpha <> 0:
    184187            self.build_smoothing_matrix_D()
    185188       
    186        
    187189        if point_coordinates:
     190
    188191            self.build_interpolation_matrix_A(point_coordinates)
    189             #self.At = self.At.tocsc()
    190             #print "self.At ",self.At
    191             #print "self.A",self.A
    192            
    193             #AtA = self.At.todense() * self.A.todense()
    194             #print "AtA dense",AtA
    195             #AtA = sparse.dok_matrix(AtA)
    196             AtA = self.At * self.A
    197             #print "AtA",AtA
     192
    198193            if self.alpha <> 0:
    199                 self.B = AtA + self.alpha*self.D
     194                self.B = self.AtA + self.alpha*self.D
    200195            else:
    201                 self.B = AtA
    202         #print "end of building b"
     196                self.B = self.AtA
     197
    203198       
    204199    def build_interpolation_matrix_A(self, point_coordinates):
    205200        """Build n x m interpolation matrix, where
    206201        n is the number of data points and
    207         m is the number of basis functions phi_k (one per vertex)
    208         """
     202        m is the number of basis functions phi_k (one per vertex)"""
     203       
    209204        #Convert input to Numeric arrays
    210205        point_coordinates = array(point_coordinates).astype(Float)
     
    216211        #self.A = zeros((n,m), Float)
    217212        self.A = sparse.dok_matrix()
     213        self.AtA = sparse.dok_matrix()
    218214
    219215        #Compute matrix elements
     
    249245
    250246                    #Assign values to matrix A
    251                    
     247
     248
    252249                    j0 = self.mesh.triangles[k,0] #Global vertex id
    253                     self.A[i, j0] = sigma0
     250                    #self.A[i, j0] = sigma0
    254251
    255252                    j1 = self.mesh.triangles[k,1] #Global vertex id
    256                     self.A[i, j1] = sigma1
     253                    #self.A[i, j1] = sigma1
    257254
    258255                    j2 = self.mesh.triangles[k,2] #Global vertex id
    259                     self.A[i, j2] = sigma2
    260 
    261                    
    262         #self.A = self.A.todense()       
    263         #Precompute
    264         self.At = self.A.transp()
     256                    #self.A[i, j2] = sigma2
     257
     258                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
     259                    js     = [j0,j1,j2]
     260
     261                    for j in js:
     262                        self.A[i,j] = sigmas[j]
     263                        for k in js:
     264                            self.AtA[j,k] += sigmas[j]*sigmas[k]
     265
     266       
     267        self.A   = (self.A).tocsc()
     268        self.AtA = (self.AtA).tocsc()
     269        self.At  =  self.A.transp()
    265270
    266271    def get_A(self):
     
    354359            self.D[v2,v0] += e20
    355360            self.D[v0,v2] += e20             
    356            
     361
     362        self.D = (self.D).tocsc()
    357363           
    358364    def fit(self, z):
    359         """Fit a smooth surface to given data points z.
     365        """Fit a smooth surface to given 1d array of data points z.
    360366
    361367        The smooth surface is computed at each vertex in the underlying
     
    366372         
    367373        Inputs:
    368           z: Vector or array of data at the point_coordinates.
    369           If z is an array, smoothing will be done for each column
     374          z: Single 1d vector or array of data at the point_coordinates.
    370375        """
    371376
    372377        #Convert input to Numeric arrays
    373378        z = array(z).astype(Float)
     379
     380
     381        if len(z.shape) > 1 :
     382            raise VectorShapeError, 'Can only deal with 1d data vector'
    374383       
    375384        #Compute right hand side based on data
     
    392401
    393402        return conjugate_gradient(self.B, Atz, Atz)
    394         #FIXME: Should we store the result here for later use? (ON)
     403        #FIXME: Should we store the result here for later use? (ON)       
    395404   
    396405    def fit_points(self, z):
     
    432441          If f is an array, interpolation will be done for each column
    433442        """
    434        
    435         return self.A * f
    436 
     443
     444        try:
     445            return self.A * f
     446        except ValueError:
     447            # We are here, so it probalby means that f is 2 dimensional
     448            # and so will do each column separately due to problems in
     449            # sparse matrix, 2d array multiplication
     450           
     451            (N , M) = self.A.shape
     452            f = array(f).astype(Float)
     453            (m , n) = f.shape
     454            if m != M :
     455                raise VectorShapeError, 'Mismatch between A and f dimensions'
     456           
     457            y = zeros( (N,n), Float)
     458
     459            for i in range(n):
     460                y[:,i] = self.A * f[:,i]
     461
     462            return y
    437463     
    438464    #FIXME: We will need a method 'evaluate(self):' that will interpolate
  • inundation/ga/storm_surge/pyvolution/test_least_squares.py

    r436 r441  
    1717class TestCase(unittest.TestCase):
    1818
    19     def setUp(self):
    20         pass
    21        
    22     def tearDown(self):
    23         pass
     19##     def setUp(self):
     20##         pass
     21       
     22##     def tearDown(self):
     23##         pass
    2424
    2525
    2626    def test_datapoint_at_centroid(self):
     27       
    2728        a = [0.0, 0.0]
    2829        b = [0.0, 2.0]
     
    3940
    4041    def test_datapoints_at_vertices(self):
    41         """Test that data points coinciding with vertices yield a diagonal matrix
    42         """
    43         print "aa"
     42
     43
    4444        a = [0.0, 0.0]
    4545        b = [0.0, 2.0]
     
    5858
    5959    def test_datapoints_on_edge_midpoints(self):
    60         """Try datapoints midway on edges -
    61         each point should affect two matrix entries equally
    62         """
    63         print "bb"
    64        
     60
    6561        a = [0.0, 0.0]
    6662        b = [0.0, 2.0]
     
    7975
    8076    def test_datapoints_on_edges(self):
    81         """Try datapoints on edges -
    82         each point should affect two matrix entries in proportion
    83         """
    84         print "cc"
    8577       
    8678        a = [0.0, 0.0]
     
    9991
    10092    def test_arbitrary_datapoints(self):
    101         """Try arbitrary datapoints
    102         """
     93
    10394
    10495        from Numeric import sum
     
    117108       
    118109
    119     # this causes a memory error in scipy.sparce       
    120     def bad_test_more_triangles(self):
    121         print "dd"
     110    def test_more_triangles(self):
     111
    122112        a = [-1.0, 0.0]
    123113        b = [3.0, 4.0]
     
    131121
    132122        #Data points
    133         data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
     123        data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3], [-1.0,-1.5 ], [1.0,-1.0]]
    134124        interp = Interpolation(points, triangles, data_points)
    135125       
     
    147137
    148138    def test_smooth_attributes_to_mesh(self):
    149         print "ee"
     139       
    150140        a = [0.0, 0.0]
    151141        b = [0.0, 5.0]
     
    170160       
    171161    def test_smooth_att_to_meshII(self):
    172         print "ff"
     162       
    173163        a = [0.0, 0.0]
    174164        b = [0.0, 5.0]
     
    189179   
    190180    def test_smooth_attributes_to_meshIII(self):
    191         print "gg"
     181
    192182        a = [-1.0, 0.0]
    193183        b = [3.0, 4.0]
     
    221211
    222212    def test_smooth_attributes_to_meshIV(self):
    223         print "hh"
    224         """ Testing 2 attributes smoothed to the mesh
    225         """
    226         print "test_smooth_attributes_to_meshIV"
     213
    227214        a = [0.0, 0.0]
    228215        b = [0.0, 5.0]
     
    246233       
    247234    def test_interpolate_attributes_to_points(self):
     235       
    248236        v0 = [0.0, 0.0]
    249237        v1 = [0.0, 5.0]
     
    268256       
    269257    def test_interpolate_attributes_to_pointsII(self):
     258
    270259        a = [-1.0, 0.0]
    271260        b = [3.0, 4.0]
     
    301290   
    302291    def test_interpolate_attributes_to_pointsIII(self):
     292       
    303293        v0 = [0.0, 0.0]
    304294        v1 = [0.0, 5.0]
     
    328318   
    329319    def test_interpolate_attributes_to_pointsIV(self):
     320
    330321        a = [-1.0, 0.0]
    331322        b = [3.0, 4.0]
     
    366357       
    367358    def test_smooth_attributes_to_mesh_function(self):
    368         """ Testing 2 attributes smoothed to the mesh
    369         """
    370         """Test multiple attributes
    371         """
    372359       
    373360        a = [0.0, 0.0]
     
    433420
    434421    def test_smoothing_matrix_more_triangles(self):
     422
     423
    435424        from Numeric import dot
    436425       
     
    476465
    477466    def test_fit_and_interpolation(self):
     467
    478468        from mesh import Mesh
    479469       
     
    663653if __name__ == "__main__":
    664654    suite = unittest.makeSuite(TestCase,'test')
    665     runner = unittest.TextTestRunner()
     655    runner = unittest.TextTestRunner(verbosity=2)
    666656    runner.run(suite)
    667657
Note: See TracChangeset for help on using the changeset viewer.