source: inundation/ga/storm_surge/pyvolution/test_least_squares.py @ 342

#!/usr/bin/env python

import unittest
from math import sqrt

from least_squares import *
from Numeric import allclose, array

def distance(x, y):
    return sqrt( sum( (array(x)-array(y))**2 ))         

def linear_function(point):
    point = array(point)
    return point[:,0]+point[:,1]
    
        
class TestCase(unittest.TestCase):

    def setUp(self):
        pass
        
    def tearDown(self):
        pass


    def test_datapoint_at_centroid(self):
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        points = [a, b, c]
        vertices = [ [1,0,2] ]   #bac

        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point        
        
        interp = Interpolation(points, vertices, data)
        assert allclose(interp.A, [[1./3, 1./3, 1./3]]) 
        


    def test_datapoints_at_vertices(self):
        """Test that data points coinciding with vertices yield a diagonal matrix
        """
        
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        points = [a, b, c]
        vertices = [ [1,0,2] ]   #bac

        data = points #Use data at vertices        
        
        interp = Interpolation(points, vertices, data)
        assert allclose(interp.A, [[1., 0., 0.],
                                   [0., 1., 0.],
                                   [0., 0., 1.]]) 
        


    def test_datapoints_on_edge_midpoints(self):
        """Try datapoints midway on edges -
        each point should affect two matrix entries equally
        """
        
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        points = [a, b, c]
        vertices = [ [1,0,2] ]   #bac

        data = [ [0., 1.], [1., 0.], [1., 1.] ]
        
        interp = Interpolation(points, vertices, data)
        
        assert allclose(interp.A, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2


    def test_datapoints_on_edges(self):
        """Try datapoints on edges -
        each point should affect two matrix entries in proportion
        """
        
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        points = [a, b, c]
        vertices = [ [1,0,2] ]   #bac

        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
        
        interp = Interpolation(points, vertices, data)
        
        assert allclose(interp.A, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2

    def test_arbitrary_datapoints(self):
        """Try arbitrary datapoints
        """

        from Numeric import sum
        
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        points = [a, b, c]
        vertices = [ [1,0,2] ]   #bac

        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
        
        interp = Interpolation(points, vertices, data)

        assert allclose(sum(interp.A, axis=1), 1.0)
        

           
    def test_more_triangles(self):
        a = [-1.0, 0.0]
        b = [3.0, 4.0]
        c = [4.0,1.0]
        d = [-3.0, 2.0] #3
        e = [-1.0,-2.0]
        f = [1.0, -2.0] #5

        points = [a, b, c, d,e,f]
        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc

        #Data points
        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]

        interp = Interpolation(points, triangles, data) 
        
        answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],   #Affects point d     
                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f
        assert allclose(interp.A, answer)




    def test_smooth_attributes_to_mesh(self):
        a = [0.0, 0.0]
        b = [0.0, 5.0]
        c = [5.0, 0.0]
        points = [a, b, c]
        triangles = [ [1,0,2] ]   #bac
         
        d1 = [1.0, 1.0]
        d2 = [1.0, 3.0]
        d3 = [3.0,1.0]
        z1 = 2
        z2 = 4
        z3 = 4
        data_coords = [d1, d2, d3] 
        
        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
        z = [z1, z2, z3]
        f = interp.fit(z)
        answer = [0, 5., 5.]
        
        assert allclose(f, answer)
        
        
    def test_smooth_att_to_meshII(self):
        a = [0.0, 0.0]
        b = [0.0, 5.0]
        c = [5.0, 0.0]
        points = [a, b, c]
        triangles = [ [1,0,2] ]   #bac
         
        d1 = [1.0, 1.0]
        d2 = [1.0, 2.0]
        d3 = [3.0,1.0]
        data_coords = [d1, d2, d3]       
        z = linear_function(data_coords)
        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
        f = interp.fit(z)
        answer = linear_function(points)
        
        assert allclose(f, answer)
    
    def test_smooth_attributes_to_meshIII(self):
        a = [-1.0, 0.0]
        b = [3.0, 4.0]
        c = [4.0,1.0]
        d = [-3.0, 2.0] #3
        e = [-1.0,-2.0]
        f = [1.0, -2.0] #5

        vertices = [a, b, c, d,e,f]
        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc

        point_coords = [[-2.0, 2.0],
                        [-1.0, 1.0],
                        [0.0,2.0],
                        [1.0, 1.0],
                        [2.0, 1.0],
                        [0.0,0.0],
                        [1.0, 0.0],
                        [0.0, -1.0],
                        [-0.2,-0.5],
                        [-0.9, -1.5],
                        [0.5, -1.9],
                        [3.0,1.0]]
        
        z = linear_function(point_coords)
        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
        f = interp.fit(z)
        answer = linear_function(vertices)
        assert allclose(f, answer)
        

    def test_smooth_attributes_to_meshIV(self):
        """ Testing 2 attributes smoothed to the mesh
        """
        a = [0.0, 0.0]
        b = [0.0, 5.0]
        c = [5.0, 0.0]
        points = [a, b, c]
        triangles = [ [1,0,2] ]   #bac
         
        d1 = [1.0, 1.0]
        d2 = [1.0, 3.0]
        d3 = [3.0,1.0]
        z1 = [2,4]
        z2 = [4,8]
        z3 = [4,8]
        data_coords = [ d1, d2, d3] 
        
        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
        z = [z1, z2, z3]
        f =  interp.fit(z)
        answer = [[0,0], [5., 10.], [5., 10.]]
        assert allclose(f, answer)
        
    def test_interpolate_attributes_to_points(self):
        v0 = [0.0, 0.0]
        v1 = [0.0, 5.0]
        v2 = [5.0, 0.0]
        
        vertices = [v0, v1, v2]
        triangles = [ [1,0,2] ]   #bac
         
        d0 = [1.0, 1.0]
        d1 = [1.0, 2.0]
        d2 = [3.0,1.0]
        point_coords = [ d0, d1, d2]
        
        interp = Interpolation(vertices, triangles, point_coords)
        f = linear_function(vertices)
        z = interp.interpolate(f)
        answer = linear_function(point_coords)
        

        assert allclose(z, answer)

       
    def test_interpolate_attributes_to_pointsII(self):
        a = [-1.0, 0.0]
        b = [3.0, 4.0]
        c = [4.0,1.0]
        d = [-3.0, 2.0] #3
        e = [-1.0,-2.0]
        f = [1.0, -2.0] #5

        vertices = [a, b, c, d,e,f]
        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc

         
        point_coords = [[-2.0, 2.0],
                        [-1.0, 1.0],
                        [0.0,2.0],
                        [1.0, 1.0],
                        [2.0, 1.0],
                        [0.0,0.0],
                        [1.0, 0.0],
                        [0.0, -1.0],
                        [-0.2,-0.5],
                        [-0.9, -1.5],
                        [0.5, -1.9],
                        [3.0,1.0]]
        
        interp = Interpolation(vertices, triangles, point_coords)
        f = linear_function(vertices)
        z = interp.interpolate(f)
        answer = linear_function(point_coords)
        #print "z",z
        #print "answer",answer
        assert allclose(z, answer)

        
        
    def test_smooth_attributes_to_mesh_function(self):
        """ Testing 2 attributes smoothed to the mesh
        """
        """Test multiple attributes
        """
        
        a = [0.0, 0.0]
        b = [0.0, 5.0]
        c = [5.0, 0.0]
        points = [a, b, c]
        triangles = [ [1,0,2] ]   #bac
         
        d1 = [1.0, 1.0]
        d2 = [1.0, 3.0]
        d3 = [3.0,1.0]
        z1 = [2,4]
        z2 = [4,8]
        z3 = [4,8]
        data_coords = [ d1, d2, d3]
        z = [z1, z2, z3]
        
        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
        answer = [[0, 0], [5., 10.], [5., 10.]]
        
        assert allclose(f, answer)



    #Tests of smoothing matrix    
    def test_smoothing_matrix_one_triangle(self):
        from Numeric import dot
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        points = [a, b, c]
        
        vertices = [ [1,0,2] ]   #bac

        interp = Interpolation(points, vertices)

        assert allclose(interp.D, [[1, -0.5, -0.5],
                                   [-0.5, 0.5, 0],
                                   [-0.5, 0, 0.5]])

        #Define f(x,y) = x
        f = array([0,0,2]) #Value at global vertex 2

        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
        #           int 1 dx dy = area = 2
        assert dot(dot(f, interp.D), f) == 2
        
        #Define f(x,y) = y
        f = array([0,2,0])  #Value at global vertex 1

        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
        #           int 1 dx dy = area = 2
        assert dot(dot(f, interp.D), f) == 2

        #Define f(x,y) = x+y
        f = array([0,2,2])  #Values at global vertex 1 and 2

        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
        #           int 2 dx dy = 2*area = 4
        assert dot(dot(f, interp.D), f) == 4        
        


    def test_smoothing_matrix_more_triangles(self):
        from Numeric import dot
        
        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, daf, dae 
        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]

        interp = Interpolation(points, vertices)


        #assert allclose(interp.D, [[1, -0.5, -0.5],
        #                           [-0.5, 0.5, 0],
        #                           [-0.5, 0, 0.5]])

        #Define f(x,y) = x
        f = array([0,0,2,0,2,4]) #f evaluated at points a-f 

        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
        #           int 1 dx dy = total area = 8
        assert dot(dot(f, interp.D), f) == 8
        
        #Define f(x,y) = y
        f = array([0,2,0,4,2,0]) #f evaluated at points a-f 

        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
        #           int 1 dx dy = area = 8
        assert dot(dot(f, interp.D), f) == 8

        #Define f(x,y) = x+y
        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   

        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
        #           int 2 dx dy = 2*area = 16
        assert dot(dot(f, interp.D), f) == 16        


    def test_fit_and_interpolation(self):
        from mesh import Mesh
        
        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, daf, dae 
        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]

        #Get (enough) datapoints
        data_points = [[ 0.66666667, 0.66666667],
                       [ 1.33333333, 1.33333333],
                       [ 2.66666667, 0.66666667],
                       [ 0.66666667, 2.66666667],
                       [ 0.0, 1.0],
                       [ 0.0, 3.0],
                       [ 1.0, 0.0],
                       [ 1.0, 1.0],
                       [ 1.0, 2.0],
                       [ 1.0, 3.0],                                          
                       [ 2.0, 1.0],
                       [ 3.0, 0.0],
                       [ 3.0, 1.0]]                                          
        
        interp = Interpolation(points, triangles, data_points, alpha=0.0)

        z = linear_function(data_points)
        answer = linear_function(points)

        f = interp.fit(z)

        assert allclose(f, answer)

        #Map back
        z1 = interp.interpolate(f)
        assert allclose(z, z1)


    def test_smoothing_and_interpolation(self):
        
        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, daf, dae 
        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]

        #Get (too few!) datapoints
        data_points = [[ 0.66666667, 0.66666667],
                       [ 1.33333333, 1.33333333],
                       [ 2.66666667, 0.66666667],
                       [ 0.66666667, 2.66666667]]

        z = linear_function(data_points)
        answer = linear_function(points)        

        #Make interpolator with too few data points and no smoothing
        interp = Interpolation(points, triangles, data_points, alpha=0.0)
        #Must raise an exception        
        try:
            f = interp.fit(z)
        except:
            pass

        #Now try with smoothing parameter
        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
        
        f = interp.fit(z)

        #f will be different from answerr due to smoothing
        #assert allclose(f, answer)

        #Map back
        z1 = interp.interpolate(f)
        assert allclose(z, z1)
        
#-------------------------------------------------------------
if __name__ == "__main__":
    suite = unittest.makeSuite(TestCase,'test')
    runner = unittest.TextTestRunner()
    runner.run(suite)