source: inundation-numpy-branch/fit_interpolate/benchmark_least_squares.py @ 3514

"""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.
   
   Duncan Gray
   Geoscience Australia, 2004.
"""


import os
import sys
import time
from random import seed, random

from anuga.pyvolution.least_squares import Interpolation
from pmesh.mesh import Mesh

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 BenchmarkLeastSquares:


    def __init__(self):
        pass

    def trial(self,
              num_of_points=20000,
              maxArea=1000,
              max_points_per_cell=4,
              is_fit=True,
              save=False):
        '''
        num_of_points 
        '''
        
        #print "num_of_points",num_of_points
        #print "maxArea",maxArea
        #print "max_points_per_cell", max_points_per_cell

        # make a mesh
        # pretty regular size, with some segments thrown in. 
        m = Mesh()
        m.addUserVertex(0,0)
        m.addUserVertex(1.0,0)
        m.addUserVertex(0,1.0)
        m.addUserVertex(1.0,1.0)
        
        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.generateMesh(mode = "Q", maxArea = maxArea)       
        if save is True:
            m.export_mesh_file("aaaa.tsh")
        mesh_dict =  m.Mesh2IOTriangulationDict()
        #print "mesh_dict",mesh_dict
        points = []
        point_atts = []
        vertex_atts = []

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

        # There has to be a better way of doing this..
        for vertex in mesh_dict['vertices']:
            vertex_atts.append(10.0)
            
        #Initial time and memory
        t0 = time.time()
        #m0 = None on windows
        m0 = mem_usage()

        interp = Interpolation(mesh_dict['vertices'],
                               mesh_dict['triangles'], points,
                               alpha=0.2, expand_search=True,
                               verbose = False,
                               max_points_per_cell = 4)
        if is_fit is True:
            calc = interp.fit_points(point_atts)

        else:
            # run an interploate problem.
            print "Interpolate!"
            calc = interp.interpolate(vertex_atts)
            
        time_taken_sec = (time.time()-t0)
        m1 = mem_usage()
        if m0 is None or m1 is None:
            memory_used = None
        else:
            memory_used = (m1 - m0)
        #print 'That took %.2f seconds' %time_taken_sec
        return time_taken_sec, memory_used, len(mesh_dict['triangles'])