Changeset 8578 for trunk/anuga_core/source/anuga/fit_interpolate/fit.py

Timestamp:

Sep 14, 2012, 9:56:39 PM (13 years ago)

Author:

steve

Message:

Added extra unit tests for set_stage and set_elevation operators

File:

• trunk/anuga_core/source/anuga/fit_interpolate/fit.py (modified) (21 diffs)

trunk/anuga_core/source/anuga/fit_interpolate/fit.py

@@ -46 +46 @@
 
 import numpy as num
+import sys
 
 
@@ -100 +101 @@
                  triangles,
                  mesh,
-                 mesh_origin,
-                 verbose)
+                 mesh_origin=mesh_origin,
+                 verbose=verbose)
         
         m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
@@ -110 +111 @@
         self.point_count = 0
         if self.alpha <> 0:
-            if verbose: log.critical('Building smoothing matrix')
-            self._build_smoothing_matrix_D()
-            
+            if verbose: log.critical('Fit: Building smoothing matrix')
+            self._build_smoothing_matrix_D(verbose=verbose)
+
+        if verbose: log.critical('Fit: Get Boundary Polygon')
         bd_poly = self.mesh.get_boundary_polygon()    
         self.mesh_boundary_polygon = ensure_numeric(bd_poly)
+
+        if verbose: log.critical('Fit: Done')
+
             
     def _build_coefficient_matrix_B(self,
@@ -126 +131 @@
         """
 
+        if verbose:
+            print 'Fit: Build Coefficient Matrix B'
+
+
         if self.alpha <> 0:
             #if verbose: log.critical('Building smoothing matrix')
             #self._build_smoothing_matrix_D()
+            # FIXME SR: This is quite time consuming.
+            # As AtA and D have same structure it should be possible
+            # to speed this up.
             self.B = self.AtA + self.alpha*self.D
         else:
             self.B = self.AtA
 
+
+        if verbose:
+            print 'Fit: Convert Coefficient Matrix B to Sparse_CSR'
+            
         # Convert self.B matrix to CSR format for faster matrix vector
         self.B = Sparse_CSR(self.B)
 
-    def _build_smoothing_matrix_D(self):
+    def _build_smoothing_matrix_D(self, verbose=False):
         """Build m x m smoothing matrix, where
         m is the number of basis functions phi_k (one per vertex)
@@ -169 +185 @@
         self.D = Sparse(m,m)
 
+        if verbose :
+            print '['+60*' '+']',
+            sys.stdout.flush()
+
+        N = len(self.mesh)
+        M = N/60
+
         # For each triangle compute contributions to D = D1+D2
-        for i in range(len(self.mesh)):
+        for i in xrange(N):
+
+            if verbose and i%M == 0 :
+                #restart_line()
+                print '\r['+(i/M)*'.'+(60-(i/M))*' ' +']',
+                sys.stdout.flush()
 
             # Get area
@@ -212 +240 @@
             self.D[v2,v0] += e20
             self.D[v0,v2] += e20
+
+        if verbose:
+            print ''
 
     def get_D(self):
@@ -242 +273 @@
         The number of attributes of the data points does not change
         """
-        
+
+
+        if verbose:
+            print 'Fit: Build Matrix AtA Atz'
+
         # Build n x m interpolation matrix
         if self.AtA == None:
@@ -271 +306 @@
         # Compute matrix elements for points inside the mesh
         triangles = self.mesh.triangles # Shorthand
-        for d, i in enumerate(inside_indices):
+
+
+        if verbose :
+            print '['+60*' '+']',
+            sys.stdout.flush()
+
+        m = n/60
+
+
+        for i in xrange(n):
+            d = inside_indices[i]
             # For each data_coordinate point
             # if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d'
                                                             # %(d, n))
-            x = point_coordinates[i]
+
+
+            if verbose and i%m == 0 :
+                print '\r['+(i/m)*'.'+(60-(i/m))*' '+']',
+                sys.stdout.flush()
+
+
+            x = point_coordinates[d]
             
             element_found, sigma0, sigma1, sigma2, k = \
@@ -302 +354 @@
                 
                 # data point has fallen within a hole - so ignore it.
-                
+
+
+        if verbose:
+            print ''
+            
         self.AtA = AtA
 
@@ -324 +380 @@
           
         """
-        
+
+
+        if verbose:
+            print 'Fit.fit: Initializing'
+
         # Use blocking to load in the point info
         if isinstance(point_coordinates_or_filename, basestring):
@@ -357 +417 @@
                 points = geo_block.get_data_points(absolute=True)
                 z = geo_block.get_attributes(attribute_name=attribute_name)
-                self.build_fit_subset(points, z, verbose=verbose)
+                self.build_fit_subset(points, z, attribute_name, verbose)
 
                 # FIXME(Ole): I thought this test would make sense here
@@ -368 +428 @@
         else:
             point_coordinates =  point_coordinates_or_filename
-            
+
+
+
         if point_coordinates is None:
-            if verbose: log.critical('Warning: no data points in fit')
+            if verbose: log.critical('Fit.fit: Warning: no data points in fit')
             msg = 'No interpolation matrix.'
             assert self.AtA is not None, msg
@@ -381 +443 @@
             # if isinstance(point_coordinates,Geospatial_data) and z is None:
             # z will come from the geo-ref
-            self.build_fit_subset(point_coordinates, z, verbose)
+            self.build_fit_subset(point_coordinates, z, verbose=verbose)
+
+
+
+
 
         # Check sanity
@@ -393 +459 @@
             msg += 'positive value,\ne.g. 1.0e-3.'
             raise TooFewPointsError(msg)
+
 
         self._build_coefficient_matrix_B(verbose)
@@ -408 +475 @@
             #raise VertsWithNoTrianglesError(msg)
         
-        
+        if verbose:
+            print 'Fit.fit: Solve Fitting Equation'
+
         return conjugate_gradient(self.B, self.Atz, self.Atz,
                                   imax=2*len(self.Atz) )
@@ -575 +644 @@
 
         if verbose: log.critical('FitInterpolate: Building mesh')
-        mesh = Mesh(vertex_coordinates, triangles)
-        mesh.check_integrity()
+        mesh = Mesh(vertex_coordinates, triangles, verbose=verbose)
+        #mesh.check_integrity() # expensive
     
     
@@ -648 +717 @@
         old_title_list = mesh_dict['vertex_attribute_titles']
 
-    if verbose: log.critical('tsh file %s loaded' % mesh_file)
+    if verbose: log.critical('Fit_to_Mesh_File: tsh file %s loaded' % mesh_file)
 
     # load in the points file
@@ -668 +737 @@
         mesh_origin = None
 
-    if verbose: log.critical("points file loaded")
-    if verbose: log.critical("fitting to mesh")
+    if verbose: log.critical("Fit_to_Mesh_File: points file loaded")
+    if verbose: log.critical("Fit_to_Mesh_File: fitting to mesh")
     f = fit_to_mesh(point_coordinates,
                     vertex_coordinates,
@@ -679 +748 @@
                     data_origin = None,
                     mesh_origin = mesh_origin)
-    if verbose: log.critical("finished fitting to mesh")
+    if verbose: log.critical("Fit_to_Mesh_File: finished fitting to mesh")
 
     # convert array to list of lists
@@ -688 +757 @@
         old_title_list.extend(title_list)
         #FIXME can this be done a faster way? - DSG
-        for i in range(len(old_point_attributes)):
+        for i in xrange(len(old_point_attributes)):
             old_point_attributes[i].extend(new_point_attributes[i])
         mesh_dict['vertex_attributes'] = old_point_attributes
@@ -696 +765 @@
         mesh_dict['vertex_attribute_titles'] = title_list
 
-    if verbose: log.critical("exporting to file %s" % mesh_output_file)
+    if verbose: log.critical("Fit_to_Mesh_File: exporting to file %s" % mesh_output_file)
 
     try: