"""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
import time


from pyvolution.least_squares import Interpolation
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 ztest_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=True,
                               verbose = False,
                               max_points_per_cell = 4)
        calc = interp.fit_points(point_atts )
        #print "interp.expanded_quad_searches", interp.expanded_quad_searches
        #print "calc", calc
        print 'That took %.2f seconds' %(time.time()-t0)

        
    def test_expand_search3(self):
        from random import seed, random
        from pmesh.mesh import Mesh

        num_of_points = 20000
        
        seed(2)
        m = Mesh()
        m.addUserVertex(0,0)
        m.addUserVertex(100,0)
        m.addUserVertex(0,100)
        m.addUserVertex(100,100)
        
        m.autoSegment(alpha = 100 )
        
        dict = {}
        dict['points'] = [[10,10],[90,20]]
        dict['segments'] = [[0,1]] 
        dict['segment_tags'] = ['wall1']   
        m.addVertsSegs(dict)
    
        dict = {}
        dict['points'] = [[10,90],[40,20]]
        dict['segments'] = [[0,1]] 
        dict['segment_tags'] = ['wall2']   
        m.addVertsSegs(dict)
        
        dict = {}
        dict['points'] = [[20,90],[60,60]]
        dict['segments'] = [[0,1]] 
        dict['segment_tags'] = ['wall3'] 
        m.addVertsSegs(dict)
        
        dict = {}
        dict['points'] = [[60,20],[90,90]]
        dict['segments'] = [[0,1]] 
        dict['segment_tags'] = ['wall4']   
        m.addVertsSegs(dict)
        
        m.addVertsSegs(dict)
        m.generateMesh(mode = "Q", maxArea = "1000")
        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=True,
                               verbose = True,
                               max_points_per_cell = 4)
        calc = interp.fit_points(point_atts )
        print "interp.expanded_quad_searches", interp.expanded_quad_searches
        #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)