source: inundation/fit_interpolate/spike_test_least_squares.py @ 2891

#!/usr/bin/env python

#TEST
import sys
import unittest
from math import sqrt


from spike_least_squares import *
from Numeric import allclose, array, transpose
from Numeric import zeros, take, compress, array, Float, Int, dot, transpose, concatenate, ArrayType
from utilities.sparse import Sparse, Sparse_CSR

from coordinate_transforms.geo_reference import Geo_reference

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 Test_Least_Squares(unittest.TestCase):

    def setUp(self):
        pass

    def tearDown(self):
        pass

    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]

        z = [z1, z2, z3]
        interp = Interpolation(points, triangles, z,data_coords, alpha=0)
        #print "interp.get_A()", interp.get_A()
        A = interp.A
        #print "A",A
        #print "z",z 
        Atz = A.trans_mult(z)
        #print "Atz",Atz
        
        f = interp.fit(z)
        answer = [0, 5., 5.]

        #print "f\n",f
        #print "answer\n",answer

        assert allclose(f, answer, atol=1e-7)

       
    def test_smooth_attributes_to_mesh_one_point(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]

        z = [z1]
        interp = Interpolation(points, triangles, z,data_coords, alpha=0)
        #print "interp.get_A()", interp.get_A()
        A = interp.A
        #print "A",A
        #print "z",z 
        Atz_actual = A.trans_mult(z)
        #Atz = interp.Atz.todense()

        
        #print "Atz",Atz
        #print "Atz_actual",Atz_actual 
        

        #assert allclose(Atz_actual, Atz, atol=1e-7)

    def ytest_chewin_the_fat(self):
        
        A = Sparse(2,2)
        A[0,0] = 1
        A[1,0] = 1
        A[0,1] = 0
        A[1,1] = 1
        z = [1,1]
        m = 2
        n = 2
        Atz = zeros((n), Float)
        r = A.trans_mult(z)
        for i in range(m):
            for j in range(n):
                Atz[i] += A[m-i-1,n-j-1]*z[i]
        print "A.trans_mult(z)",A.trans_mult(z)
        print "Atz",Atz 
        print "A*z", A*z
        

        #print "answer\n",answer

        assert allclose(f, answer, atol=1e-7)


    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)
        #print "z",z 

        interp = Interpolation(points, triangles, z,
                               data_coords, alpha=0.0)
        f = interp.fit(z)
        answer = linear_function(points)
        #print "f\n",f
        #print "answer\n",answer

        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, z,
                               point_coords, alpha=0.0)

        #print 'z',z
        f = interp.fit(z)
        answer = linear_function(vertices)
        #print "f\n",f
        #print "answer\n",answer
        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]

        z = [z1, z2, z3]
        interp = Interpolation(points, triangles, z,
                               data_coords, alpha=0.0)
        f =  interp.fit_points(z)
        answer = [[0,0], [5., 10.], [5., 10.]]
        assert allclose(f, answer)


    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]]

        z = linear_function(data_points)
        interp = Interpolation(points, triangles, z,
                               data_points, alpha=0.0)

        answer = linear_function(points)

        f = interp.fit(z)

        #print "f",f
        #print "answer",answer
        assert allclose(f, answer)

        
    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, z,
                               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, z,
                               data_points, alpha=1.0e-13)

        f = interp.fit(z)
        #f will be different from answer due to smoothing
        assert allclose(f, answer,atol=5)

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