Changeset 4108

Timestamp:

Dec 20, 2006, 3:58:01 PM (16 years ago)

Author:

duncan

Message:

mainly comments - Getting ready for blocking..

Location:

anuga_core/source/anuga/fit_interpolate

Files:

• benchmark_least_squares.py (modified) (1 diff)
• fit.py (modified) (8 diffs)

anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py

@@ -93 +93 @@
                 print "Interpolate!"
                 calc = interp.interpolate(mesh_dict['vertex_attributes'])
-        else:
-            # need to change to fit_interpolate code
-            interp = Fit(mesh_dict['vertices'],
+        else: 
+            if is_fit is True:
+                interp = Fit(mesh_dict['vertices'],
                                  mesh_dict['triangles'], 
-                                 max_vertices_per_cell = max_points_per_cell) 
-            if is_fit is True:
+                                 max_vertices_per_cell = max_points_per_cell)
                 print "Fit in Fit"
                 calc = interp.fit(points_dict['points'],
                                   points_dict['point_attributes'])
-                pass
             else:
                 # run an interploate problem.

anuga_core/source/anuga/fit_interpolate/fit.py

@@ -78 +78 @@
           Note: Don't supply a vertex coords as a geospatial object and
               a mesh origin, since geospatial has its own mesh origin.
+
+
+        Usage,
+        To use this in a blocking way, call  build_fit_subset, with z info,
+        and then fit, with no point coord, z info.
+        
         """
         # Initialise variabels
@@ -232 +238 @@
 
         if self.AtA == None:
-            # AtA and Atz need ot be initialised.
+            # AtA and Atz need to be initialised.
             m = self.mesh.number_of_nodes
             if len(z.shape) > 1:
@@ -243 +249 @@
 
             self.AtA = Sparse(m,m)
+            # The memory damage has been done by now.
+            
         self.point_count += point_coordinates.shape[0]
         #print "_build_matrix_AtA_Atz - self.point_count", self.point_count
@@ -250 +258 @@
         #      self.mesh.get_boundary_polygon()
 
+        # Why are these global?
         self.inside_poly_indices, self.outside_poly_indices  = \
                      in_and_outside_polygon(point_coordinates,
@@ -311 +320 @@
             #FIXME (DSG) - do  a message
         else:
+            # This is where build fit subset can be looped over,
+            # if a file name is passed in.
             point_coordinates = ensure_absolute(point_coordinates,
                                                 geo_reference=point_origin)
+            #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)
 
@@ -345 +358 @@
 
         
-    def build_fit_subset(self, point_coordinates, z,
-                              verbose = False):
+    def build_fit_subset(self, point_coordinates, z=None, attribute_name=None,
+                              verbose=False):
         """Fit a smooth surface to given 1d array of data points z.
 
@@ -356 +369 @@
               List of coordinate pairs [x, y] of
               data points or an nx2 Numeric array or a Geospatial_data object
-          z: Single 1d vector or array of data at the point_coordinates.
-
-        """
-        #Note: Don't get the z info from anuga.geospatial_data.attributes yet.
-        # If we did fit would have to handle attribute title info.
+        z: Single 1d vector or array of data at the point_coordinates.
+        attribute_name: Used to get the z values from the
+              geospatial object if no attribute_name is specified,
+              it's a bit of a lucky dip as to what attributes you get.
+              If there is only one attribute it will be that one.
+
+        """
 
         #FIXME(DSG-DSG): Check that the vert and point coords
@@ -369 +384 @@
         
         #Convert input to Numeric arrays
-        z = ensure_numeric(z, Float)
+        if z is not None:
+            z = ensure_numeric(z, Float)
+        else:
+            msg = 'z not specified'
+            assert isinstance(point_coordinates,Geospatial_data), msg
+            z = point_coordinates.get_attributes(attribute_name)
+            
         point_coordinates = ensure_numeric(point_coordinates, Float)