"""
This will produce a mesh which has vertex 10 that has no triangle.

Why?
   
   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
   Geoscience Australia, 2004.
"""


import unittest
from math import sqrt
import time


from anuga.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 test_expand_search3(self):
        from random import seed, random
        from pmesh.mesh import Mesh

        num_of_points = 100
        
        seed(2)
        m = Mesh()
        m.addUserVertex(0,0)
        m.addUserVertex(100,0)
        m.addUserVertex(0,100)
        m.addUserVertex(100,100)
        
        m.auto_segment(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")
        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.auto_segment(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)