Changeset 4341

Timestamp:

Mar 29, 2007, 5:54:23 PM (18 years ago)

Author:

ole

Message:

Refactored file_function so that instances of domain is no longer pass into cached internal function _file_function. This should fix the issues with caching for good (I hope).

All tests pass.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/test_util.py (modified) (1 diff)
• abstract_2d_finite_volumes/util.py (modified) (6 diffs)
• utilities/numerical_tools.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -1293 +1293 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Util,'test')
-    suite = unittest.makeSuite(Test_Util,'test_get_data_from_file')
+    suite = unittest.makeSuite(Test_Util,'test')
+    #suite = unittest.makeSuite(Test_Util,'test_get_data_from_file')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -50 +50 @@
     quantities - the name of the quantity to be interpolated or a
                  list of quantity names. The resulting function will return
-                 a tuple of values - one for each quantity  
+                 a tuple of values - one for each quantity
+                 If quantities are None, domain's conserved quantities
+                 are used.
 
     interpolation_points - list of absolute UTM coordinates for points (N x 2)
@@ -63 +65 @@
 
 
+    #FIXME (OLE): Should check origin of domain against that of file
+    #In fact, this is where origin should be converted to that of domain
+    #Also, check that file covers domain fully.
+
+    #Take into account:
+    #- domain's georef
+    #- sww file's georef
+    #- interpolation points as absolute UTM coordinates
+
+
+
+
+    # Use domain's conserved_quantity names as defaults
+    if domain is not None:    
+        if quantities is None: 
+            quantities = domain.conserved_quantities
+            
+        domain_starttime = domain.starttime
+    else:
+        domain_starttime = None
+
+
     # Build arguments and keyword arguments for use with caching or apply.
     args = (filename,)
-    
-    kwargs = {'domain': domain,
-              'quantities': quantities,
+
+
+    # FIXME (Ole): Caching this function will not work well
+    # if domain is passed in as instances change hash code.
+    # Instead we pass in those attributes that are needed (and return them
+    # if modified)
+    kwargs = {'quantities': quantities,
               'interpolation_points': interpolation_points,
+              'domain_starttime': domain_starttime,
               'time_thinning': time_thinning,                   
               'verbose': verbose}
@@ -82 +111 @@
             raise msg
 
-        f = cache(_file_function,
-                  args, kwargs,
-                  dependencies=[filename],
-                  compression=False,                  
-                  verbose=verbose)
+        f, starttime = cache(_file_function,
+                             args, kwargs,
+                             dependencies=[filename],
+                             compression=False,                  
+                             verbose=verbose)
 
     else:
-        f = apply(_file_function,
-                  args, kwargs)
+        f, starttime = apply(_file_function,
+                             args, kwargs)
 
 
     #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
     #structure
-        
-
-    return f
-
-
-
-def _file_function(filename,
-                   domain=None,
-                   quantities=None,
-                   interpolation_points=None,
-                   time_thinning=1,                                                
-                   verbose=False):
-    """Internal function
-    
-    See file_function for documentatiton
-    """
-    
-
-    #FIXME (OLE): Should check origin of domain against that of file
-    #In fact, this is where origin should be converted to that of domain
-    #Also, check that file covers domain fully.
-
-    #Take into account:
-    #- domain's georef
-    #- sww file's georef
-    #- interpolation points as absolute UTM coordinates
-
-
-    assert type(filename) == type(''),\
-               'First argument to File_function must be a string'
-
-    try:
-        fid = open(filename)
-    except Exception, e:
-        msg = 'File "%s" could not be opened: Error="%s"'\
-                  %(filename, e)
-        raise msg
-
-    line = fid.readline()
-    fid.close()
-
-    if quantities is None: 
-        if domain is not None:
-            quantities = domain.conserved_quantities
-
-
-
-    if line[:3] == 'CDF':
-        return get_netcdf_file_function(filename, domain, quantities,
-                                        interpolation_points,
-                                        time_thinning=time_thinning,
-                                        verbose=verbose)
-    else:
-        raise 'Must be a NetCDF File'
-
-
-
-def get_netcdf_file_function(filename,
-                             domain=None,
-                             quantity_names=None,
-                             interpolation_points=None,
-                             time_thinning=1,                             
-                             verbose=False):
-    """Read time history of spatial data from NetCDF sww file and
-    return a callable object f(t,x,y)
-    which will return interpolated values based on the input file.
-
-    If domain is specified, model time (domain.starttime)
-    will be checked and possibly modified
-    
-    All times are assumed to be in UTC
-
-    See Interpolation function for further documetation
-    
-    """
-    
-    
-    #FIXME: Check that model origin is the same as file's origin
-    #(both in UTM coordinates)
-    #If not - modify those from file to match domain
-    #Take this code from e.g. dem2pts in data_manager.py
-    #FIXME: Use geo_reference to read and write xllcorner...
-        
-
-    #FIXME: Maybe move caching out to this level rather than at the
-    #Interpolation_function level (below)
-
-    import time, calendar, types
-    from anuga.config import time_format
-    from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
-
-    #Open NetCDF file
-    if verbose: print 'Reading', filename
-    fid = NetCDFFile(filename, 'r')
-
-    if type(quantity_names) == types.StringType:
-        quantity_names = [quantity_names]        
-    
-    if quantity_names is None or len(quantity_names) < 1:
-        #If no quantities are specified get quantities from file
-        #x, y, time are assumed as the independent variables so
-        #they are excluded as potentiol quantities
-        quantity_names = []
-        for name in fid.variables.keys():
-            if name not in ['x', 'y', 'time']:
-                quantity_names.append(name)
-
-    if len(quantity_names) < 1:                
-        msg = 'ERROR: At least one independent value must be specified'
-        raise msg
-
-
-    if interpolation_points is not None:
-        interpolation_points = ensure_absolute(interpolation_points)
-        msg = 'Points must by N x 2. I got %d' %interpolation_points.shape[1]
-        assert interpolation_points.shape[1] == 2, msg
-
-
-    #Now assert that requested quantitites (and the independent ones)
-    #are present in file 
-    missing = []
-    for quantity in ['time'] + quantity_names:
-        if not fid.variables.has_key(quantity):
-            missing.append(quantity)
-
-    if len(missing) > 0:
-        msg = 'Quantities %s could not be found in file %s'\
-              %(str(missing), filename)
-        fid.close()
-        raise Exception, msg
-
-    #Decide whether this data has a spatial dimension
-    spatial = True
-    for quantity in ['x', 'y']:
-        if not fid.variables.has_key(quantity):
-            spatial = False
-
-    if filename[-3:] == 'tms' and spatial is True:
-        msg = 'Files of type tms must not contain spatial information'
-        raise msg
-
-    if filename[-3:] == 'sww' and spatial is False:
-        msg = 'Files of type sww must contain spatial information'        
-        raise msg        
-
-    #Get first timestep
-    try:
-        starttime = fid.starttime[0]
-    except ValueError:
-        msg = 'Could not read starttime from file %s' %filename
-        raise msg
-
-    #Get variables
-    #if verbose: print 'Get variables'    
-    time = fid.variables['time'][:]    
-
-    # Get time independent stuff
-    if spatial:
-        #Get origin
-        xllcorner = fid.xllcorner[0]
-        yllcorner = fid.yllcorner[0]
-        zone = fid.zone[0]        
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        triangles = fid.variables['volumes'][:]
-
-        x = reshape(x, (len(x),1))
-        y = reshape(y, (len(y),1))
-        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
-
-        if interpolation_points is not None:
-            #Adjust for georef
-            interpolation_points[:,0] -= xllcorner
-            interpolation_points[:,1] -= yllcorner        
-        
-
 
 
     if domain is not None:
-        #Update domain.startime if it is *earlier* than starttime
+        #Update domain.startime if it is *earlier* than starttime from file
         if starttime > domain.starttime:
             msg = 'WARNING: Start time as specified in domain (%f)'\
@@ -289 +140 @@
                %domain.starttime
 
-
-        #If domain.startime is *later* than starttime,
+    return f
+
+
+
+def _file_function(filename,
+                   quantities=None,
+                   interpolation_points=None,
+                   domain_starttime=None,
+                   time_thinning=1,
+                   verbose=False):
+    """Internal function
+    
+    See file_function for documentatiton
+    """
+    
+
+    assert type(filename) == type(''),\
+               'First argument to File_function must be a string'
+
+    try:
+        fid = open(filename)
+    except Exception, e:
+        msg = 'File "%s" could not be opened: Error="%s"'\
+                  %(filename, e)
+        raise msg
+
+    line = fid.readline()
+    fid.close()
+
+
+    if line[:3] == 'CDF':
+        return get_netcdf_file_function(filename,
+                                        quantities,
+                                        interpolation_points,
+                                        domain_starttime,
+                                        time_thinning=time_thinning,
+                                        verbose=verbose)
+    else:
+        raise 'Must be a NetCDF File'
+
+
+
+def get_netcdf_file_function(filename,
+                             quantity_names=None,
+                             interpolation_points=None,
+                             domain_starttime=None,                            
+                             time_thinning=1,                             
+                             verbose=False):
+    """Read time history of spatial data from NetCDF sww file and
+    return a callable object f(t,x,y)
+    which will return interpolated values based on the input file.
+
+    Model time (domain_starttime)
+    will be checked, possibly modified and returned
+    
+    All times are assumed to be in UTC
+
+    See Interpolation function for further documetation
+    
+    """
+    
+    
+    #FIXME: Check that model origin is the same as file's origin
+    #(both in UTM coordinates)
+    #If not - modify those from file to match domain
+    #(origin should be passed in)
+    #Take this code from e.g. dem2pts in data_manager.py
+    #FIXME: Use geo_reference to read and write xllcorner...
+        
+
+    import time, calendar, types
+    from anuga.config import time_format
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
+
+    #Open NetCDF file
+    if verbose: print 'Reading', filename
+    fid = NetCDFFile(filename, 'r')
+
+    if type(quantity_names) == types.StringType:
+        quantity_names = [quantity_names]        
+    
+    if quantity_names is None or len(quantity_names) < 1:
+        #If no quantities are specified get quantities from file
+        #x, y, time are assumed as the independent variables so
+        #they are excluded as potentiol quantities
+        quantity_names = []
+        for name in fid.variables.keys():
+            if name not in ['x', 'y', 'time']:
+                quantity_names.append(name)
+
+    if len(quantity_names) < 1:                
+        msg = 'ERROR: At least one independent value must be specified'
+        raise msg
+
+
+    if interpolation_points is not None:
+        interpolation_points = ensure_absolute(interpolation_points)
+        msg = 'Points must by N x 2. I got %d' %interpolation_points.shape[1]
+        assert interpolation_points.shape[1] == 2, msg
+
+
+    #Now assert that requested quantitites (and the independent ones)
+    #are present in file 
+    missing = []
+    for quantity in ['time'] + quantity_names:
+        if not fid.variables.has_key(quantity):
+            missing.append(quantity)
+
+    if len(missing) > 0:
+        msg = 'Quantities %s could not be found in file %s'\
+              %(str(missing), filename)
+        fid.close()
+        raise Exception, msg
+
+    #Decide whether this data has a spatial dimension
+    spatial = True
+    for quantity in ['x', 'y']:
+        if not fid.variables.has_key(quantity):
+            spatial = False
+
+    if filename[-3:] == 'tms' and spatial is True:
+        msg = 'Files of type tms must not contain spatial information'
+        raise msg
+
+    if filename[-3:] == 'sww' and spatial is False:
+        msg = 'Files of type sww must contain spatial information'        
+        raise msg        
+
+    #Get first timestep
+    try:
+        starttime = fid.starttime[0]
+    except ValueError:
+        msg = 'Could not read starttime from file %s' %filename
+        raise msg
+
+    #Get variables
+    #if verbose: print 'Get variables'    
+    time = fid.variables['time'][:]    
+
+    # Get time independent stuff
+    if spatial:
+        #Get origin
+        xllcorner = fid.xllcorner[0]
+        yllcorner = fid.yllcorner[0]
+        zone = fid.zone[0]        
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        triangles = fid.variables['volumes'][:]
+
+        x = reshape(x, (len(x),1))
+        y = reshape(y, (len(y),1))
+        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+
+        if interpolation_points is not None:
+            #Adjust for georef
+            interpolation_points[:,0] -= xllcorner
+            interpolation_points[:,1] -= yllcorner        
+        
+
+
+
+    if domain_starttime is not None:
+
+        #If domain_startime is *later* than starttime,
         #move time back - relative to domain's time
-        if domain.starttime > starttime:
-            time = time - domain.starttime + starttime
+        if domain_starttime > starttime:
+            time = time - domain_starttime + starttime
 
         #FIXME Use method in geo to reconcile
@@ -325 +340 @@
         vertex_coordinates = triangles = interpolation_points = None         
 
+
+    # Return Interpolation_function instance as well as
+    # starttime for use to possible modify that of domain
     return Interpolation_function(time,
                                   quantities,
@@ -332 +350 @@
                                   interpolation_points,
                                   time_thinning=time_thinning,
-                                  verbose=verbose)
+                                  verbose=verbose), starttime
+
+    # NOTE (Ole): Caching Interpolation function is too slow as
+    # the very long parameters need to be hashed.
+
+
+
 
 

anuga_core/source/anuga/utilities/numerical_tools.py

@@ -163 +163 @@
 
     x = ensure_numeric(x)
-    y = ensure_numeric(y)       
-    assert(len(x)==len(y))
+    y = ensure_numeric(y)
+    msg = 'Lengths must be equal: len(x) == %d, len(y) == %d' %(len(x), len(y))
+    assert(len(x)==len(y)), msg
 
     N = len(x)