source: inundation/fit_interpolate/benchmark_least_squares.py @ 2891

"""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 pyvolution.least_squares import Interpolation
from fit_interpolate.interpolate import Interpolate
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,
              use_least_squares=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

        mesh_dict = self._build_regular_mesh_dict(maxArea=maxArea,
                                                  save=save)
        points_dict = self._build_points_dict(num_of_points=num_of_points)
            
        #Initial time and memory
        t0 = time.time()
        #m0 = None on windows
        m0 = mem_usage()

        if use_least_squares is True:
            interp = Interpolation(mesh_dict['vertices'],
                                   mesh_dict['triangles'],
                                   points_dict['points'],
                                   expand_search=True,
                                   verbose = False,
                                   max_points_per_cell = max_points_per_cell) 
            if is_fit is True:
                calc = interp.fit_points(points_dict['point_attributes'])
                
            else:
                # run an interploate problem.
                print "Interpolate!"
                calc = interp.interpolate(mesh_dict['vertex_attributes'])
        else:
            # need to change to fit_interpolate code
            interp = Interpolate(mesh_dict['vertices'],
                                 mesh_dict['triangles'], 
                                 max_vertices_per_cell = max_points_per_cell) 
            if is_fit is True:
                #calc = interp.fit_points(point_atts)
                pass
            else:
                # run an interploate problem.
                print "Interpolate!"
                calc = interp.interpolate(mesh_dict['vertex_attributes']
                                          ,points_dict['points']
                                          ,start_blocking_len = 500000)
            
        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'])

    def _build_regular_mesh_dict(self,
                                 maxArea=1000,
                                 save=False):
      # 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.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.generateMesh(mode = "Q", maxArea = maxArea)       
        if save is True:
            m.export_mesh_file("aaaa.tsh")
        mesh_dict =  m.Mesh2IOTriangulationDict()

        #Add vert attribute info to the mesh
        mesh_dict['vertex_attributes'] = []
        # There has to be a better way of doing this..
        for vertex in mesh_dict['vertices']:
            mesh_dict['vertex_attributes'].append([10.0])

        return mesh_dict

    def _build_points_dict(self, num_of_points=20000):
        
        points_dict = {}
        points = []
        point_atts = []

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

        points_dict['points'] = points
        points_dict['point_attributes'] = point_atts
        return points_dict


#-------------------------------------------------------------
if __name__ == "__main__":
        b = BenchmarkLeastSquares()
        b._build_regular_mesh_dict()
        b._build_points_dict()