Changeset 1750

Timestamp:

Aug 24, 2005, 11:45:12 AM (19 years ago)

Author:

ole

Message:

Work towards new set_quantity API

Location:

inundation

Files:

• pyvolution/quantity.py (modified) (8 diffs)
• pyvolution/test_quantity.py (modified) (1 diff)
• pyvolution/util.py (modified) (1 diff)
• wiki/future_directions.txt (modified) (2 diffs)

inundation/pyvolution/quantity.py

@@ -82 +82 @@
 
     #New leaner interface to setting values
-    def set_quantity(self,
-                     numeric = None,  #List, numeric array or constant
-                     quantity = None, #Another quantity 
-                     function = None, #Callable object: f(x,y)  
-                     points = None, values = None, #Input for least squares
-                     filename = None, attribute_name = None, #Input from file
-                     alpha = None,
-                     indices = None,
-                     location = None, #Really needed
-                     verbose = None):
+    def set_values_new(self,
+                   numeric = None,  #List, numeric array or constant
+                   quantity = None, #Another quantity 
+                   function = None, #Callable object: f(x,y)  
+                   points = None, values = None, #Input for least squares
+                   filename = None, attribute_name = None, #Input from file
+                   alpha = None,
+                   location = None,                   
+                   indices = None,
+                   verbose = None):
         
         """Set values for quantity based on different sources.
-
+        
         numeric:
         
@@ -109 +109 @@
           Any callable object that takes two 1d arrays x and y
           each of length N and returns an array also of length N.
+          The function will be evaluated at points determined by
+          location and indices.
 
         points:
@@ -125 +127 @@
           Smoothing parameter to be used with least squares fits.
           See module least_squares for further details about alpha.
+          Alpha will only be used with points, values or filename.
+          Otherwise it will be ignored.
                
           
@@ -131 +135 @@
                   Default is 'vertices'
 
-        In case of location == 'centroids' the dimension values must
-        be a list of a Numerical array of length N, N being the number
-        of elements. Otherwise it must be of dimension Nx3
-
-        The values will be stored in elements following their
-        internal ordering.
-
-        If values are described a function, it will be evaluated at
-        specified points
-
-        If location is not 'unique vertices' Indices is the set of element ids
-        that the operation applies to.
-        If location is 'unique vertices' Indices is the set of vertex ids
-        that the operation applies to.
-
-        If selected location is vertices, values for centroid and edges
-        will be assigned interpolated values.
-        In any other case, only values for the specified locations
-        will be assigned and the others will be left undefined.
-        
-        """
-        pass
-
-    
-    
-    
-
-
-    #Old interface to setting values
-    def set_values(self, X,
-                   location='vertices',
-                   indexes = None):
-        """Set values for quantity
-
-        X: Compatible list, Numeric array (see below), constant or function
-        location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroids
-                  Default is "vertices"
-
-        In case of location == 'centroids' the dimension values must
-        be a list of a Numerical array of length N, N being the number
-        of elements. Otherwise it must be of dimension Nx3
-
-        The values will be stored in elements following their
-        internal ordering.
-
-        If values are described a function, it will be evaluated at
-        specified points
-
-        If indexex is not 'unique vertices' Indexes is the set of element ids
-        that the operation applies to.
-        If indexex is 'unique vertices' Indexes is the set of vertex ids
-        that the operation applies to.
-
-        If selected location is 'vertices', values for centroid and edges
-        will be assigned interpolated values.
-        In any other case, only values for the specified locations
-        will be assigned and the others will be left undefined.
-
-        """
+                  In case of location == 'centroids' the dimension values must
+                  be a list of a Numerical array of length N,
+                  N being the number of elements.
+                  Otherwise it must be of dimension Nx3
+
+
+                  The values will be stored in elements following their
+                  internal ordering.
+
+                  If location is not 'unique vertices' Indices is the
+                  set of element ids that the operation applies to.
+                  If location is 'unique vertices' Indices is the set
+                  of vertex ids that the operation applies to.
+
+                  If selected location is vertices, values for
+                  centroid and edges will be assigned interpolated
+                  values.  In any other case, only values for the
+                  specified locations will be assigned and the others
+                  will be left undefined.
+
+
+        Exactly one of the arguments
+          numeric, quantity, function, points, filename
+        must be present.
+        """
+
+        from types import FloatType, IntType, LongType, ListType, NoneType
+        from Numeric import ArrayType
+
+        #General input checks
+        msg = 'Exactly one of the arguments '+\
+              'numeric, quantity, function, points, or filename '+\
+              'must be present.'
+        L = [numeric, quantity, function, points, filename]
+        assert L.count(None) == len(L)-1, msg
+
 
         if location not in ['vertices', 'centroids', 'edges',
@@ -198 +178 @@
 
 
-        if X is None:
-            msg = 'Given values are None'
-            raise msg
-
-        import types, Numeric
-        assert type(indexes) in [types.ListType, types.NoneType,
-                                 Numeric.ArrayType],\
-                                 'Indices must be a list or None'
-
-
-        if callable(X):
-            #Use function specific method
-            self.set_function_values(X, location, indexes = indexes)
-        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
-            if location == 'centroids':
-                if (indexes ==  None):
-                    self.centroid_values[:] = X
-                else:
-                    #Brute force
-                    for i in indexes:
-                        self.centroid_values[i,:] = X
-
-            elif location == 'edges':
-                if (indexes ==  None):
-                    self.edge_values[:] = X
-                else:
-                    #Brute force
-                    for i in indexes:
-                        self.edge_values[i,:] = X
-
-            elif location == 'unique vertices':
-                if (indexes ==  None):
-                    self.edge_values[:] = X
-                else:
-
-                    #Go through list of unique vertices
-                    for unique_vert_id in indexes:
-                        triangles = self.domain.vertexlist[unique_vert_id]
-
-                        #In case there are unused points
-                        if triangles is None: continue
-
-                        #Go through all triangle, vertex pairs
-                        #and set corresponding vertex value
-                        for triangle_id, vertex_id in triangles:
-                            self.vertex_values[triangle_id, vertex_id] = X
-
-                        #Intialise centroid and edge_values
-                        self.interpolate()
+        msg = 'Indices must be a list or None'
+        assert type(indices) in [ListType, NoneType, ArrayType], msg
+               
+
+        
+        #Determine which 'set_values_from_...' to use
+
+        if numeric is not None:
+            if type(numeric) in [FloatType, IntType, LongType]:
+                self.set_values_from_constant(numeric,
+                                              location, indices, verbose)
+            elif type(numeric) in [ArrayType, ListType]:
+                self.set_values_from_array(numeric,
+                                           location, indices, verbose)
+            elif callable(numeric):
+                self.set_values_from_function(numeric,
+                                              location, indices, verbose)
             else:
-                if (indexes ==  None):
-                    self.vertex_values[:] = X
-                else:
-                    #Brute force
-                    for i_vertex in indexes:
-                        self.vertex_values[i_vertex,:] = X
-
-        elif type(X) in [Numeric.ArrayType, types.ListType]:
-            #Use array specific method
-            self.set_array_values(X, location, indexes = indexes)
-        elif type(X) == types.StringType:
-            #Assume X is a filename
-            self.set_values_from_file(X) #FIXME: More parameters
-
-
-        if location == 'vertices' or location == 'unique vertices':
-            #Intialise centroid and edge_values
-            self.interpolate()
-
-        if location == 'centroids':
-            #Extrapolate 1st order - to capture notion of area being specified
-            self.extrapolate_first_order()
-
-    def get_values(self, location='vertices', indexes = None):
-        """get values for quantity
-
-        return X, Compatible list, Numeric array (see below)
-        location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
-                  Default is "vertices"
-
-        In case of location == 'centroids' the dimension values must
-        be a list of a Numerical array of length N, N being the number
-        of elements. Otherwise it must be of dimension Nx3
-
-        The returned values with be a list the length of indexes
-        (N if indexes = None).  Each value will be a list of the three
-        vertex values for this quantity.
-
-        Indexes is the set of element ids that the operation applies to.
-
-        """
-        from Numeric import take
-
-        if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
-            msg = 'Invalid location: %s' %location
-            raise msg
-
-        import types, Numeric
-        assert type(indexes) in [types.ListType, types.NoneType,
-                                 Numeric.ArrayType],\
-                                 'Indices must be a list or None'
-
-        if location == 'centroids':
-            if (indexes ==  None):
-                indexes = range(len(self))
-            return take(self.centroid_values,indexes)
-        elif location == 'edges':
-            if (indexes ==  None):
-                indexes = range(len(self))
-            return take(self.edge_values,indexes)
-        elif location == 'unique vertices':
-            if (indexes ==  None):
-                indexes=range(self.domain.coordinates.shape[0])
-            vert_values = []
-            #Go through list of unique vertices
-            for unique_vert_id in indexes:
-                triangles = self.domain.vertexlist[unique_vert_id]
-
-                #In case there are unused points
-                if triangles is None:
-                    msg = 'Unique vertex not associated with triangles'
-                    raise msg
-
-                # Go through all triangle, vertex pairs
-                # Average the values
-                sum = 0
-                for triangle_id, vertex_id in triangles:
-                    sum += self.vertex_values[triangle_id, vertex_id]
-                vert_values.append(sum/len(triangles))
-            return Numeric.array(vert_values)
-        else:
-            if (indexes ==  None):
-                indexes = range(len(self))
-            return take(self.vertex_values,indexes)
-
-
-    def set_function_values(self, f, location='vertices', indexes = None):
+                msg = 'Illegal type for argument numeric: %s' %str(numeric)
+                raise msg
+        elif quantity is not None:
+            self.set_values_from_quantity(quantity,
+                                          location, indices, verbose)
+        elif function is not None:
+            msg = 'Argument function must be callable'
+            assert callable(function), msg            
+            self.set_values_from_function(function,
+                                          location, indices, verbose)
+        elif points is not None:
+            msg = 'When points are specified, associated values must also be.'
+            assert values is not None, msg
+            self.set_values_from_points(points, values, alpha,
+                                        location, indices, verbose)
+        elif filename is not None:
+            self.set_values_from_file(filename, attribute_name, alpha,
+                                      location, indices, verbose)
+        else:
+            raise 'This can\'t happen :-)'
+                
+
+    
+    #Specific functions for setting values
+    def set_values_from_function(self, function,
+                                 location, indices, verbose):
         """Set values for quantity using specified function
 
         f: x, y -> z Function where x, y and z are arrays
         location: Where values are to be stored.
-                  Permissible options are: vertices, centroid
+                  Permissible options are: vertices, centroid, edges,
+                  unique vertices
                   Default is "vertices"
         """
 
         #FIXME: Should check that function returns something sensible and
-        #raise a meaningfol exception if it returns None for example
+        #raise a meaningfull exception if it returns None for example
 
         from Numeric import take
 
-        if (indexes ==  None):
+        if (indexes is None):
             indexes = range(len(self))
             is_subset = False
@@ -473 +360 @@
 
 
-    def set_values_from_file(self, filename,
-                             attribute_name = None,
-                             verbose = False):
-        """Set quantity based on arbitrary points in .pts file using least_squares
-
-        attribute_name selects name of attribute present in file.
+    def set_values_from_points(self, points, values, alpha,
+                               location, indices, verbose):
+        """
+        """
+
+        
+        #FIXME: Needs unit test         
+        from util import ensure_numeric
+        from least_squares import fit_to_mesh
+        
+        points = ensure_numeric(points, Float)
+        values = ensure_numeric(values, Float)
+
+        if location != 'vertices':
+            msg = 'set_values_from_points is only defined for'+\
+                  'location=\'vertices\''
+            raise msg
+
+        coordinates = self.domain.coordinates
+        triangles = self.domain.triangles
+
+        #FIXME Pass and use caching here
+        vertex_attributes = fit_to_mesh(coordinates,
+                                        triangles,
+                                        points,
+                                        values,
+                                        alpha = alpha,
+                                        verbose = verbose)
+        
+        self.set_values_from_array(vertex_attributes,
+                                   location, indices, verbose)        
+        
+
+        
+        
+    
+    def set_values_from_file(self, filename, attribute_name, alpha,
+                             location, indices, verbose):
+        """Set quantity based on arbitrary points in .pts file
+        using least_squares attribute_name selects name of attribute
+        present in file. 
         If not specified try to use whatever is available in file.
         """
 
-        #FIXME location, indices
-
-
-        import util, least_squares
-
+
+        #FIXME: Needs unit test 
+        from types import StringType
+        msg = 'Filename must be a text string'
+        assert type(filename) == StringType, msg
+
+
+        #Read from (NetCDF) file
+        import util
         points, attributes = util.read_xya(filename)
-
+        
         if attribute_name is None:
             names = attributes.keys()
             attribute_name = names[0]
 
+        msg = 'Attribute_name must be a text string'
+        assert type(attribute_name) == StringType, msg            
+
 
         if verbose:
             print 'Using attribute %s from file %s' %(attribute_name, filename)
             print 'Available attributes: %s' %(names)
-
 
         try:
@@ -506 +434 @@
             raise msg
 
-
-        #Fit attributes to mesh
-        X = least_squares.fit_to_mesh(self.domain.coordinates,
-                                      self.domain.triangles,
-                                      points,
-                                      z,
-                                      verbose=verbose)
-
-
-
-        #Recursively set value using fitted array
-        print 'shape', X.shape
-        self.set_values(X, location='vertices')  #FIXME, indexes = None):
+            
+        #Call least squares method    
+        self.set_values_from_points(points, z, alpha,
+                                    location, indices, verbose)
+
+
+
+
+
+    #Old interface to setting values
+    def set_values(self, X,
+                   location='vertices',
+                   indexes = None):
+        """Set values for quantity
+
+        X: Compatible list, Numeric array (see below), constant or function
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroids
+                  Default is "vertices"
+
+        In case of location == 'centroids' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements. Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If values are described a function, it will be evaluated at
+        specified points
+
+        If indexex is not 'unique vertices' Indexes is the set of element ids
+        that the operation applies to.
+        If indexex is 'unique vertices' Indexes is the set of vertex ids
+        that the operation applies to.
+
+        If selected location is 'vertices', values for centroid and edges
+        will be assigned interpolated values.
+        In any other case, only values for the specified locations
+        will be assigned and the others will be left undefined.
+
+        """
+
+        if location not in ['vertices', 'centroids', 'edges',
+                            'unique vertices']:
+            msg = 'Invalid location: %s' %location
+            raise msg
+
+
+        if X is None:
+            msg = 'Given values are None'
+            raise msg
+
+        import types, Numeric
+        assert type(indexes) in [types.ListType, types.NoneType,
+                                 Numeric.ArrayType],\
+                                 'Indices must be a list or None'
+
+
+        if callable(X):
+            #Use function specific method
+            self.set_function_values(X, location, indexes = indexes)
+        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
+            if location == 'centroids':
+                if (indexes ==  None):
+                    self.centroid_values[:] = X
+                else:
+                    #Brute force
+                    for i in indexes:
+                        self.centroid_values[i,:] = X
+
+            elif location == 'edges':
+                if (indexes ==  None):
+                    self.edge_values[:] = X
+                else:
+                    #Brute force
+                    for i in indexes:
+                        self.edge_values[i,:] = X
+
+            elif location == 'unique vertices':
+                if (indexes ==  None):
+                    self.edge_values[:] = X
+                else:
+
+                    #Go through list of unique vertices
+                    for unique_vert_id in indexes:
+                        triangles = self.domain.vertexlist[unique_vert_id]
+
+                        #In case there are unused points
+                        if triangles is None: continue
+
+                        #Go through all triangle, vertex pairs
+                        #and set corresponding vertex value
+                        for triangle_id, vertex_id in triangles:
+                            self.vertex_values[triangle_id, vertex_id] = X
+
+                        #Intialise centroid and edge_values
+                        self.interpolate()
+            else:
+                if (indexes ==  None):
+                    self.vertex_values[:] = X
+                else:
+                    #Brute force
+                    for i_vertex in indexes:
+                        self.vertex_values[i_vertex,:] = X
+
+        elif type(X) in [Numeric.ArrayType, types.ListType]:
+            #Use array specific method
+            self.set_array_values(X, location, indexes = indexes)
+        elif type(X) == types.StringType:
+            #Assume X is a filename
+            self.set_values_from_file(X) #FIXME: More parameters
+
+
+        if location == 'vertices' or location == 'unique vertices':
+            #Intialise centroid and edge_values
+            self.interpolate()
+
+        if location == 'centroids':
+            #Extrapolate 1st order - to capture notion of area being specified
+            self.extrapolate_first_order()
+
+
+
+    def get_values(self, location='vertices', indexes = None):
+        """get values for quantity
+
+        return X, Compatible list, Numeric array (see below)
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroid
+                  Default is "vertices"
+
+        In case of location == 'centroids' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements. Otherwise it must be of dimension Nx3
+
+        The returned values with be a list the length of indexes
+        (N if indexes = None).  Each value will be a list of the three
+        vertex values for this quantity.
+
+        Indexes is the set of element ids that the operation applies to.
+
+        """
+        from Numeric import take
+
+        if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
+            msg = 'Invalid location: %s' %location
+            raise msg
+
+        import types, Numeric
+        assert type(indexes) in [types.ListType, types.NoneType,
+                                 Numeric.ArrayType],\
+                                 'Indices must be a list or None'
+
+        if location == 'centroids':
+            if (indexes ==  None):
+                indexes = range(len(self))
+            return take(self.centroid_values,indexes)
+        elif location == 'edges':
+            if (indexes ==  None):
+                indexes = range(len(self))
+            return take(self.edge_values,indexes)
+        elif location == 'unique vertices':
+            if (indexes ==  None):
+                indexes=range(self.domain.coordinates.shape[0])
+            vert_values = []
+            #Go through list of unique vertices
+            for unique_vert_id in indexes:
+                triangles = self.domain.vertexlist[unique_vert_id]
+
+                #In case there are unused points
+                if triangles is None:
+                    msg = 'Unique vertex not associated with triangles'
+                    raise msg
+
+                # Go through all triangle, vertex pairs
+                # Average the values
+                sum = 0
+                for triangle_id, vertex_id in triangles:
+                    sum += self.vertex_values[triangle_id, vertex_id]
+                vert_values.append(sum/len(triangles))
+            return Numeric.array(vert_values)
+        else:
+            if (indexes ==  None):
+                indexes = range(len(self))
+            return take(self.vertex_values,indexes)
+
+
+    def set_function_values(self, f, location='vertices', indexes = None):
+        """Set values for quantity using specified function
+
+        f: x, y -> z Function where x, y and z are arrays
+        location: Where values are to be stored.
+                  Permissible options are: vertices, centroid
+                  Default is "vertices"
+        """
+
+        #FIXME: Should check that function returns something sensible and
+        #raise a meaningfol exception if it returns None for example
+
+        from Numeric import take
+
+        if (indexes ==  None):
+            indexes = range(len(self))
+            is_subset = False
+        else:
+            is_subset = True
+        if location == 'centroids':
+            P = take(self.domain.centroid_coordinates,indexes)
+            if is_subset:
+                self.set_values(f(P[:,0], P[:,1]), location, indexes = indexes)
+            else:
+                self.set_values(f(P[:,0], P[:,1]), location)
+        elif location == 'vertices':
+            P = self.domain.vertex_coordinates
+            if is_subset:
+                #Brute force
+                for e in indexes:
+                    for i in range(3):
+                        self.vertex_values[e,i] = f(P[e,2*i], P[e,2*i+1])
+            else:
+                for i in range(3):
+                    self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
+        else:
+            raise 'Not implemented: %s' %location
+
+
+    def set_array_values(self, values, location='vertices', indexes = None):
+        """Set values for quantity
+
+        values: Numeric array
+        location: Where values are to be stored.
+        Permissible options are: vertices, edges, centroid, unique vertices
+        Default is "vertices"
+
+        indexes - if this action is carried out on a subset of
+        elements or unique vertices
+        The element/unique vertex indexes are specified here.
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements.
+
+        Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If selected location is vertices, values for centroid and edges
+        will be assigned interpolated values.
+        In any other case, only values for the specified locations
+        will be assigned and the others will be left undefined.
+        """
+
+        from Numeric import array, Float, Int, allclose
+
+        values = array(values).astype(Float)
+
+        if (indexes <> None):
+            indexes = array(indexes).astype(Int)
+            msg = 'Number of values must match number of indexes'
+            assert values.shape[0] == indexes.shape[0], msg
+
+        N = self.centroid_values.shape[0]
+
+        if location == 'centroids':
+            assert len(values.shape) == 1, 'Values array must be 1d'
+
+            if indexes == None:
+                msg = 'Number of values must match number of elements'
+                assert values.shape[0] == N, msg
+
+                self.centroid_values = values
+            else:
+                msg = 'Number of values must match number of indexes'
+                assert values.shape[0] == indexes.shape[0], msg
+
+                #Brute force
+                for i in range(len(indexes)):
+                    self.centroid_values[indexes[i]] = values[i]
+
+        elif location == 'edges':
+            assert len(values.shape) == 2, 'Values array must be 2d'
+
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
+
+            msg = 'Array must be N x 3'
+            assert values.shape[1] == 3, msg
+
+            self.edge_values = values
+
+        elif location == 'unique vertices':
+            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
+                   'Values array must be 1d'
+
+            self.set_vertex_values(values.flat, indexes = indexes)
+        else:
+            if len(values.shape) == 1:
+                self.set_vertex_values(values, indexes = indexes)
+                #if indexes == None:
+                    #Values are being specified once for each unique vertex
+                #    msg = 'Number of values must match number of vertices'
+                #    assert values.shape[0] == self.domain.coordinates.shape[0], msg
+                 #   self.set_vertex_values(values)
+                #else:
+                #    for element_index, value in map(None, indexes, values):
+                #        self.vertex_values[element_index, :] = value
+
+            elif len(values.shape) == 2:
+                #Vertex values are given as a triplet for each triangle
+
+                msg = 'Array must be N x 3'
+                assert values.shape[1] == 3, msg
+
+                if indexes == None:
+                    self.vertex_values = values
+                else:
+                    for element_index, value in map(None, indexes, values):
+                        self.vertex_values[element_index] = value
+            else:
+                msg = 'Values array must be 1d or 2d'
+                raise msg
+
+
+
 
 
@@ -729 +969 @@
 
 
+    
+
 class Conserved_quantity(Quantity):
     """Class conserved quantity adds to Quantity:

inundation/pyvolution/test_quantity.py

@@ -979 +979 @@
         indexes = [1]
         quantity.set_values(value,
-                                  location = 'centroids',
-                                  indexes = indexes)
+                            location = 'centroids',
+                            indexes = indexes)
         #print "quantity.centroid_values",quantity.centroid_values
         #print "quantity.get_values(location = 'centroids') ",\

inundation/pyvolution/util.py

@@ -125 +125 @@
         if type(A) == ArrayType:
             if A.typecode == typecode:
-                return array(A)
+                return array(A)  #FIXME: Shouldn't this be just return A?
             else:
                 return A.astype(typecode)

inundation/wiki/future_directions.txt

@@ -156 +156 @@
 
 
+Introduce create_quantity in domain.py:
+  It will make a new named instance and populate it by calling set_quantity 
+  if desired.
+  Create_quantity would be called automatically by shallow_water.         
+  
+
 Make stage appear as any other quantity:
 Either 
@@ -171 +177 @@
 
 
+
 Boundary (Dirichlet): Think about specifying values by quantity name (and having defaults)
      Dirichlet(stage = 1.0, ymomentum = 0.2)