source: inundation/fit_interpolate/benchmark_least_squares.py @ 1982

"""Least squares smooting and interpolation.

   measure the speed of least squares.

   ________________________
   General comments
   
   The max_points_per_cell does effect the time spent solving a
   problem.  The best value to use is probably dependent on the number
   of triangles.  Maybe develop a simple imperical algorithm, based on
   test results.
   
   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
import os
import sys

def mem_usage():
    '''
    returns the rss.

  RSS  The total amount of physical memory used by  the  task,  in  kilo-
            bytes,  is  shown  here.  For ELF processes used library pages are
            counted here, for a.out processes not.
            
    Only works on nix systems.
    '''
    import string
    p=os.popen('ps uwp %s'%os.getpid()) 
    lines=p.readlines()
    #print "lines", lines
    status=p.close() 
    if status or len(lines)!=2 or sys.platform == 'win32': 
        return None 
    return int(string.split(lines[1])[4]) 




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

        # takes 7.73 seconds
        #num_of_points = 20000
        #maxArea = 1000
        #max_points_per_cell = 4

        # took 9.84 secs
        #num_of_points = 20000
        #maxArea = 100
        #max_points_per_cell = 4

        # 16.61 secs
        num_of_points = 20000
        maxArea = 10 #836 basis functions
        #max_points_per_cell = 4 #16.61 sec
        
        #max_points_per_cell = 30 #31.17 sec
        #max_points_per_cell = 2 #16.125 sec
        max_points_per_cell = 8 #17.66 sec

        
        # 16.61 secs
        num_of_points = 20000
        maxArea = 1 #7917 basis functions
        max_points_per_cell = 4 # 76.047 sec
        
        #max_points_per_cell = 30 #
        #max_points_per_cell = 2 #78.97 sec
        #max_points_per_cell = 16 # 82.33 sec
        #max_points_per_cell = 8 #77.984 sec
        
        # 8.44 secs
        #num_of_points = 20000
        #maxArea = 1
        
        
        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 = maxArea)
        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, #False,
                               max_points_per_cell = max_points_per_cell)
        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__":
    import os
    suite = unittest.makeSuite(Test_Least_Squares,'test')
    runner = unittest.TextTestRunner(verbosity=1)
    runner.run(suite)
    val = mem_usage()
    print 'RSS', val