Changeset 2394

Timestamp:

Feb 14, 2006, 12:56:50 PM (19 years ago)

Author:

duncan

Message:

adding interpolation_function capability

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (5 diffs)
• test_interpolate.py (modified) (3 diffs)

inundation/fit_interpolate/interpolate.py

@@ -21 +21 @@
 
 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
-     ArrayType, allclose, take
+     ArrayType, allclose, take, NewAxis
 
 from pyvolution.mesh import Mesh
@@ -65 +65 @@
 
         """
+
+        # why no alpha in parameter list?
         
         # Initialise variabels
@@ -73 +75 @@
         triangles = ensure_numeric(triangles, Int)
         vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
-
         #Build underlying mesh
         if verbose: print 'Building mesh'
@@ -282 +283 @@
           Interpolated values at inputted points (z).
         """
-        
+        #print "point_coordinates interpolate.interpolate",point_coordinates 
         # Can I interpolate, based on previous point_coordinates? 
         if point_coordinates is None:
@@ -338 +339 @@
     def _get_point_data_z(self, f, verbose=False):
         return self._A * f
-        
+
+
+
+class Interpolation_interface:
+    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
+    which is interpolated from time series defined at vertices of
+    triangular mesh (such as those stored in sww files)
+
+    Let m be the number of vertices, n the number of triangles
+    and p the number of timesteps. 
+
+    Mandatory input
+        time:               px1 array of monotonously increasing times (Float)
+        quantities:         Dictionary of arrays or 1 array (Float) 
+                            The arrays must either have dimensions pxm or mx1.
+                            The resulting function will be time dependent in
+                            the former case while it will be constan with
+                            respect to time in the latter case.
+        
+    Optional input:
+        quantity_names:     List of keys into the quantities dictionary 
+        vertex_coordinates: mx2 array of coordinates (Float)
+        triangles:          nx3 array of indices into vertex_coordinates (Int)
+        interpolation_points: Nx2 array of coordinates to be interpolated to 
+        verbose:            Level of reporting
+    
+    
+    The quantities returned by the callable object are specified by
+    the list quantities which must contain the names of the
+    quantities to be returned and also reflect the order, e.g. for
+    the shallow water wave equation, on would have
+    quantities = ['stage', 'xmomentum', 'ymomentum']
+
+    The parameter interpolation_points decides at which points interpolated
+    quantities are to be computed whenever object is called.
+    If None, return average value
+    """
+
+    
+    
+    def __init__(self,
+                 time,
+                 quantities,
+                 quantity_names = None,  
+                 vertex_coordinates = None,
+                 triangles = None,
+                 interpolation_points = None,
+                 verbose = False):
+        """Initialise object and build spatial interpolation if required
+        """
+
+        from Numeric import array, zeros, Float, alltrue, concatenate,\
+             reshape, ArrayType
+
+
+        from util import mean, ensure_numeric
+        from config import time_format
+        import types
+
+
+        #Check temporal info
+        time = ensure_numeric(time)        
+        msg = 'Time must be a monotonuosly '
+        msg += 'increasing sequence %s' %time
+        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
+
+
+        #Check if quantities is a single array only
+        if type(quantities) != types.DictType:
+            quantities = ensure_numeric(quantities)
+            quantity_names = ['Attribute']
+
+            #Make it a dictionary
+            quantities = {quantity_names[0]: quantities}
+
+
+        #Use keys if no names are specified
+        if quantity_names is None:
+            quantity_names = quantities.keys()
+
+
+        #Check spatial info
+        if vertex_coordinates is None:
+            self.spatial = False
+        else:    
+            vertex_coordinates = ensure_numeric(vertex_coordinates)
+
+            assert triangles is not None, 'Triangles array must be specified'
+            triangles = ensure_numeric(triangles)
+            self.spatial = True            
+
+             
+        #Save for use with statistics
+        self.quantity_names = quantity_names        
+        self.quantities = quantities        
+        self.vertex_coordinates = vertex_coordinates 
+        self.interpolation_points = interpolation_points
+        self.T = time[:]  # Time assumed to be relative to starttime
+        self.index = 0    # Initial time index
+        self.precomputed_values = {}
+            
+        #Precomputed spatial interpolation if requested
+        if interpolation_points is not None:
+            if self.spatial is False:
+                raise 'Triangles and vertex_coordinates must be specified'
+            
+            try:
+                self.interpolation_points = ensure_numeric(interpolation_points)
+            except:
+                msg = 'Interpolation points must be an N x 2 Numeric array '+\
+                      'or a list of points\n'
+                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
+                                      '...')
+                raise msg
+
+
+            m = len(self.interpolation_points)
+            p = len(self.T)
+            
+            for name in quantity_names:
+                self.precomputed_values[name] = zeros((p, m), Float)
+
+            #Build interpolator
+            interpol = Interpolate(vertex_coordinates,
+                                     triangles,
+                                     #point_coordinates = \
+                                     #self.interpolation_points,
+                                     #alpha = 0,
+                                     verbose = verbose)
+
+            if verbose: print 'Interpolate'
+            for i, t in enumerate(self.T):
+                #Interpolate quantities at this timestep
+                if verbose and i%((p+10)/10)==0:
+                    print ' time step %d of %d' %(i, p)
+                    
+                for name in quantity_names:
+                    if len(quantities[name].shape) == 2:
+                        result = interpol.interpolate(quantities[name][i,:],
+                                     point_coordinates = \
+                                     self.interpolation_points)
+                    else:
+                       #Assume no time dependency 
+                       result = interpol.interpolate(quantities[name][:],
+                                     point_coordinates = \
+                                     self.interpolation_points)
+                       
+                    self.precomputed_values[name][i, :] = result
+                   
+            #Report
+            if verbose:
+                print self.statistics()
+                #self.print_statistics()
+            
+        else:
+            #Store quantitites as is
+            for name in quantity_names:
+                self.precomputed_values[name] = quantities[name]
+
+
+    def __repr__(self):
+        #return 'Interpolation function (spatio-temporal)'
+        return self.statistics()
+    
+
+    def __call__(self, t, point_id = None, x = None, y = None):
+        """Evaluate f(t), f(t, point_id) or f(t, x, y)
+
+        Inputs:
+          t: time - Model time. Must lie within existing timesteps
+          point_id: index of one of the preprocessed points.
+          x, y:     Overrides location, point_id ignored
+          
+          If spatial info is present and all of x,y,point_id
+          are None an exception is raised 
+                    
+          If no spatial info is present, point_id and x,y arguments are ignored
+          making f a function of time only.
+
+          
+          FIXME: point_id could also be a slice
+          FIXME: What if x and y are vectors?
+          FIXME: What about f(x,y) without t?
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import zeros, Float
+        from util import mean        
+
+        if self.spatial is True:
+            if point_id is None:
+                if x is None or y is None:
+                    msg = 'Either point_id or x and y must be specified'
+                    raise msg
+            else:
+                if self.interpolation_points is None:
+                    msg = 'Interpolation_function must be instantiated ' +\
+                          'with a list of interpolation points before parameter ' +\
+                          'point_id can be used'
+                    raise msg
+
+
+        msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
+        msg += ' does not match model time: %s\n' %t
+        if t < self.T[0]: raise msg
+        if t > self.T[-1]: raise msg
+
+        oldindex = self.index #Time index
+
+        #Find current time slot
+        while t > self.T[self.index]: self.index += 1
+        while t < self.T[self.index]: self.index -= 1
+
+        if t == self.T[self.index]:
+            #Protect against case where t == T[-1] (last time)
+            # - also works in general when t == T[i]
+            ratio = 0
+        else:
+            #t is now between index and index+1
+            ratio = (t - self.T[self.index])/\
+                    (self.T[self.index+1] - self.T[self.index])
+
+        #Compute interpolated values
+        q = zeros(len(self.quantity_names), Float)
+
+        for i, name in enumerate(self.quantity_names):
+            Q = self.precomputed_values[name]
+
+            if self.spatial is False:
+                #If there is no spatial info                
+                assert len(Q.shape) == 1
+
+                Q0 = Q[self.index]
+                if ratio > 0: Q1 = Q[self.index+1]
+
+            else:
+                if x is not None and y is not None:
+                    #Interpolate to x, y
+                    
+                    raise 'x,y interpolation not yet implemented'
+                else:
+                    #Use precomputed point
+                    Q0 = Q[self.index, point_id]
+                    if ratio > 0: Q1 = Q[self.index+1, point_id]
+
+            #Linear temporal interpolation    
+            if ratio > 0:
+                q[i] = Q0 + ratio*(Q1 - Q0)
+            else:
+                q[i] = Q0
+
+
+        #Return vector of interpolated values
+        #if len(q) == 1:
+        #    return q[0]
+        #else:
+        #    return q
+
+
+        #Return vector of interpolated values
+        #FIXME:
+        if self.spatial is True:
+            return q
+        else:
+            #Replicate q according to x and y
+            #This is e.g used for Wind_stress
+            if x is None or y is None: 
+                return q
+            else:
+                try:
+                    N = len(x)
+                except:
+                    return q
+                else:
+                    from Numeric import ones, Float
+                    #x is a vector - Create one constant column for each value
+                    N = len(x)
+                    assert len(y) == N, 'x and y must have same length'
+                    res = []
+                    for col in q:
+                        res.append(col*ones(N, Float))
+                        
+                return res
+
+
+    def statistics(self):
+        """Output statistics about interpolation_function
+        """
+        
+        vertex_coordinates = self.vertex_coordinates
+        interpolation_points = self.interpolation_points               
+        quantity_names = self.quantity_names
+        quantities = self.quantities
+        precomputed_values = self.precomputed_values                 
+                
+        x = vertex_coordinates[:,0]
+        y = vertex_coordinates[:,1]                
+
+        str =  '------------------------------------------------\n'
+        str += 'Interpolation_function (spatio-temporal) statistics:\n'
+        str += '  Extent:\n'
+        str += '    x in [%f, %f], len(x) == %d\n'\
+               %(min(x), max(x), len(x))
+        str += '    y in [%f, %f], len(y) == %d\n'\
+               %(min(y), max(y), len(y))
+        str += '    t in [%f, %f], len(t) == %d\n'\
+               %(min(self.T), max(self.T), len(self.T))
+        str += '  Quantities:\n'
+        for name in quantity_names:
+            q = quantities[name][:].flat
+            str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
+
+        if interpolation_points is not None:    
+            str += '  Interpolation points (xi, eta):'\
+                   ' number of points == %d\n' %interpolation_points.shape[0]
+            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
+                                            max(interpolation_points[:,0]))
+            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
+                                             max(interpolation_points[:,1]))
+            str += '  Interpolated quantities (over all timesteps):\n'
+        
+            for name in quantity_names:
+                q = precomputed_values[name][:].flat
+                str += '    %s at interpolation points in [%f, %f]\n'\
+                       %(name, min(q), max(q))
+        str += '------------------------------------------------\n'
+
+        return str
+
+def interpolate_sww(sww_file, time, interpolation_points,
+                    quantity_names = None, verbose = False):
+    """
+    obsolete.
+    use file_function in utils
+    """
+    #open sww file
+    x, y, volumes, time, quantities = read_sww(sww_file)
+    print "x",x
+    print "y",y
+    
+    print "time", time
+    print "quantities", quantities
+
+    #Add the x and y together
+    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
+
+    #Will return the quantity values at the specified times and locations
+    interp =  Interpolation_interface(
+                 time,
+                 quantities,
+                 quantity_names = quantity_names,  
+                 vertex_coordinates = vertex_coordinates,
+                 triangles = volumes,
+                 interpolation_points = interpolation_points,
+                 verbose = verbose)
+
+
+def read_sww(file_name):
+    """
+    Read in an sww file.
+    
+    Input;
+    file_name - the sww file
+    
+    Output;
+    x - Vector of x values
+    y - Vector of y values
+    z - Vector of bed elevation
+    volumes - Array.  Each row has 3 values, representing
+    the vertices that define the volume
+    time - Vector of the times where there is stage information
+    stage - array with respect to time and vertices (x,y)
+    """
+    
+    #FIXME Have this reader as part of data_manager?
+    
+    from Scientific.IO.NetCDF import NetCDFFile     
+    import tempfile
+    import sys
+    import os
+    
+    #Check contents
+    #Get NetCDF
+    
+    # see if the file is there.  Throw a QUIET IO error if it isn't
+    # I don't think I could implement the above
+    
+    #throws prints to screen if file not present
+    junk = tempfile.mktemp(".txt")
+    fd = open(junk,'w')
+    stdout = sys.stdout
+    sys.stdout = fd
+    fid = NetCDFFile(file_name, 'r') 
+    sys.stdout = stdout
+    fd.close()
+    #clean up
+    os.remove(junk)     
+      
+    # Get the variables
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:] 
+    time = fid.variables['time'][:]
+
+    keys = fid.variables.keys()
+    keys.remove("x")
+    keys.remove("y")
+    keys.remove("volumes")
+    keys.remove("time")
+     #Turn NetCDF objects into Numeric arrays
+    quantities = {}
+    for name in keys:
+        quantities[name] = fid.variables[name][:]
+            
+    fid.close()
+    return x, y, volumes, time, quantities
+
 #-------------------------------------------------------------
 if __name__ == "__main__":

inundation/fit_interpolate/test_interpolate.py

@@ -2 +2 @@
 
 #TEST
+
+#import time, os
+
+
 import sys
+import os
 import unittest
 from math import sqrt
-
-
+import tempfile
+
+from Scientific.IO.NetCDF import NetCDFFile
+from Numeric import allclose, array, transpose, zeros, Float
+
+
+# ANUGA code imports
 from interpolate import *
-from Numeric import allclose, array, transpose
-
 from coordinate_transforms.geo_reference import Geo_reference
+from shallow_water import Domain, Transmissive_boundary #, mean
+from util import mean
+from data_manager import get_dataobject
 
 def distance(x, y):
@@ -23 +34 @@
 
     def setUp(self):
-        pass
+
+        import time
+        from mesh_factory import rectangular
+
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(2, 2)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order=2
+        domain.beta_h = 0
+
+
+        #Set some field values
+        domain.set_quantity('elevation', lambda x,y: -x)
+        domain.set_quantity('friction', 0.03)
+
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        ######################
+        #Initial condition - with jumps
+
+        bed = domain.quantities['elevation'].vertex_values
+        stage = zeros(bed.shape, Float)
+
+        h = 0.3
+        for i in range(stage.shape[0]):
+            if i % 2 == 0:
+                stage[i,:] = bed[i,:] + h
+            else:
+                stage[i,:] = bed[i,:]
+
+        domain.set_quantity('stage', stage)
+
+        domain.distribute_to_vertices_and_edges()
+
+
+        self.domain = domain
+
+        C = domain.get_vertex_coordinates()
+        self.X = C[:,0:6:2].copy()
+        self.Y = C[:,1:6:2].copy()
+
+        self.F = bed
+
+
 
     def tearDown(self):
@@ -530 +592 @@
         
 
+
+    def test_interpolation_interface_time_only(self):
+        """Test spatio-temporal interpolation
+        Test that spatio temporal function performs the correct
+        interpolations in both time and space
+        """
+
+
+        #Three timesteps
+        time = [1.0, 5.0, 6.0]
+        
+
+        #One quantity
+        Q = zeros( (3,6), Float )
+
+        #Linear in time and space
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        
+        for i, t in enumerate(time):
+            Q[i, :] = t*linear_function(points)
+
+            
+        #Check basic interpolation of one quantity using averaging
+        #(no interpolation points or spatial info)
+        from util import mean        
+        I = Interpolation_interface(time, [mean(Q[0,:]),
+                                          mean(Q[1,:]),
+                                          mean(Q[2,:])])
+
+
+
+        #Check temporal interpolation
+        for i in [0,1,2]:
+            assert allclose(I(time[i]), mean(Q[i,:]))
+
+        #Midway    
+        assert allclose(I( (time[0] + time[1])/2 ),
+                        (I(time[0]) + I(time[1]))/2 )
+
+        assert allclose(I( (time[1] + time[2])/2 ),
+                        (I(time[1]) + I(time[2]))/2 )
+
+        assert allclose(I( (time[0] + time[2])/2 ),
+                        (I(time[0]) + I(time[2]))/2 )                 
+
+        #1/3
+        assert allclose(I( (time[0] + time[2])/3 ),
+                        (I(time[0]) + I(time[2]))/3 )                         
+
+
+        #Out of bounds checks
+        try:
+            I(time[0]-1) 
+        except:
+            pass
+        else:
+            raise 'Should raise exception'
+
+        try:
+            I(time[-1]+1) 
+        except:
+            pass
+        else:
+            raise 'Should raise exception'        
+
+
+        
+
+    def test_interpolation_interface_spatial_only(self):
+        """Test spatio-temporal interpolation with constant time
+        """
+
+        #Three timesteps
+        time = [1.0, 5.0, 6.0]
+        
+        
+        #Setup mesh used to represent fitted function
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        #New datapoints where interpolated values are sought
+        interpolation_points = [[ 0.0, 0.0],
+                                [ 0.5, 0.5],
+                                [ 0.7, 0.7],
+                                [ 1.0, 0.5],
+                                [ 2.0, 0.4],
+                                [ 2.8, 1.2]]
+
+
+        #One quantity linear in space
+        Q = linear_function(points)
+
+
+        #Check interpolation of one quantity using interpolaton points
+        I = Interpolation_interface(time, Q,
+                                   vertex_coordinates = points,
+                                   triangles = triangles, 
+                                   interpolation_points = interpolation_points,
+                                   verbose = False)
+
+
+        answer = linear_function(interpolation_points)
+
+        t = time[0]
+        for j in range(50): #t in [1, 6]
+            for id in range(len(interpolation_points)):
+                assert allclose(I(t, id), answer[id])
+
+            t += 0.1    
+
+
+        try:    
+            I(1)
+        except:
+            pass
+        else:
+            raise 'Should raise exception'
+            
+
+
+    def test_interpolation_interface(self):
+        """Test spatio-temporal interpolation
+        Test that spatio temporal function performs the correct
+        interpolations in both time and space
+        """
+
+
+        #Three timesteps
+        time = [1.0, 5.0, 6.0]
+        
+
+        #Setup mesh used to represent fitted function
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        #New datapoints where interpolated values are sought
+        interpolation_points = [[ 0.0, 0.0],
+                                [ 0.5, 0.5],
+                                [ 0.7, 0.7],
+                                [ 1.0, 0.5],
+                                [ 2.0, 0.4],
+                                [ 2.8, 1.2]]
+
+
+        #One quantity
+        Q = zeros( (3,6), Float )
+
+        #Linear in time and space
+        for i, t in enumerate(time):
+            Q[i, :] = t*linear_function(points)
+
+
+        #Check interpolation of one quantity using interpolaton points)
+        I = Interpolation_interface(time, Q,
+                                   vertex_coordinates = points,
+                                   triangles = triangles, 
+                                   interpolation_points = interpolation_points,
+                                   verbose = False)
+
+
+        answer = linear_function(interpolation_points)
+        
+        t = time[0]
+        for j in range(50): #t in [1, 6]
+            for id in range(len(interpolation_points)):
+                assert allclose(I(t, id), t*answer[id])
+
+            t += 0.1    
+
+        try:    
+            I(1)
+        except:
+            pass
+        else:
+            raise 'Should raise exception'
+
+        
+    def BADtest_interpolate_sww(self):
+        """Not a unit test, rather a system test for interpolate_sww
+        This function is obsolete
+        """
+        
+        self.domain.filename = 'datatest' + str(time.time())
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.reduction = mean
+
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+        self.domain.time = 2.
+        sww.store_timestep('stage')
+
+        #print "self.domain.filename",self.domain.filename
+        interp = interpolate_sww(sww.filename, [0.0, 2.0],
+                                 [[0,1],[0.5,0.5]],
+                                 ['stage'])
+        #assert allclose(interp,[0.0,2.0])
+
+        #Cleanup
+        os.remove(sww.filename)
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":