source: inundation/fit_interpolate/benchmark_least_squares.py @ 1938

"""Least squares smooting and interpolation.

   measure the speed of least squares.
   
   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
   Geoscience Australia, 2004.
"""


import unittest
from math import sqrt


from least_squares import *
from Numeric import allclose, array, transpose

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_expand_search2(self):
        from random import seed, random
        from pmesh.mesh import Mesh
        
        seed(2)
        m = Mesh()
        verts_per_sector = 100 #100
        num_of_points = 300
        sec_side = 9.0

        vert_atts = []
        for vert in range(verts_per_sector):
            m.addUserVertex(random()*sec_side, random()*sec_side)
            vert_atts.append(random())
        for vert in range(verts_per_sector):
            m.addUserVertex(90.0+random()*sec_side, 90.0+random()*sec_side)
            vert_atts.append(random())
        for vert in range(verts_per_sector):
            m.addUserVertex(random()*sec_side, 90.0+random()*sec_side)
            vert_atts.append(random())
        for vert in range(verts_per_sector):
            m.addUserVertex(90.0+random()*sec_side, random()*sec_side)
            vert_atts.append(random())
        #for vert in range(verts_per_sector):
        #    m.addUserVertex(40.0+random()*sec_side, 40.0+random()*sec_side)
        #    vert_atts.append(random())
            
        m.autoSegment(alpha = 100 )
        m.generateMesh(mode = "Q")
        m.export_mesh_file("aaaa.tsh")
        mesh_dict =  m.Mesh2IOTriangulationDict()

        points = []
        point_atts = []
        for point in range(num_of_points):
            points.append([random()*100, random()*100])
            point_atts.append(10.0)
        m.autoSegment(alpha = 100 )

        #print 'points', points
        t0 = time.time()
        interp = Interpolation(mesh_dict['vertices'],
                               mesh_dict['triangles'], points,
                               alpha=0.2, expand_search=False, #True,
                               verbose = False,
                               max_points_per_cell = 4)
        calc = interp.fit_points(point_atts )
        #print "calc", calc
        print 'That took %.2f seconds' %(time.time()-t0)
        
#-------------------------------------------------------------
if __name__ == "__main__":
    suite = unittest.makeSuite(Test_Least_Squares,'test')
    runner = unittest.TextTestRunner(verbosity=1)
    runner.run(suite)