source: branches/source_numpy_conversion/anuga/fit_interpolate/benchmark_least_squares.py @ 6768

Last change on this file since 6768 was 5905, checked in by rwilson, 16 years ago

NumPy? conversion.
File size: 12.6 KB
Line 
1"""Least squares smooting and interpolation.
2
3   measure the speed of least squares.
4
5   ________________________
6   General comments
7   
8   The max_points_per_cell does effect the time spent solving a
9   problem.  The best value to use is probably dependent on the number
10   of triangles.  Maybe develop a simple imperical algorithm, based on
11   test results.
12   
13   Duncan Gray
14   Geoscience Australia, 2004.
15"""
16
17
18import os
19import sys
20import time
21from random import seed, random
22import tempfile
23import profile , pstats
24##from math import sqrt
25##import numpy
26
27from anuga.fit_interpolate.search_functions import search_times, \
28     reset_search_times
29
30from anuga.fit_interpolate.interpolate import Interpolate
31from anuga.fit_interpolate.fit import Fit
32from anuga.pmesh.mesh import Mesh
33from anuga.geospatial_data.geospatial_data import Geospatial_data
34from anuga.shallow_water import Domain
35from anuga.fit_interpolate.fit import Fit, fit_to_mesh
36from anuga.fit_interpolate.interpolate import benchmark_interpolate
37from anuga.coordinate_transforms.geo_reference import Geo_reference
38from anuga.fit_interpolate.general_fit_interpolate import \
39     get_build_quadtree_time
40
41"""
42
43Code from the web;
44
45from ctypes import *
46from ctypes.wintypes import DWORD
47
48SIZE_T = c_ulong
49
50class _MEMORYSTATUS(Structure):
51_fields_ = [("dwLength", DWORD),
52("dwMemoryLength", DWORD),
53("dwTotalPhys", SIZE_T),
54("dwAvailPhys", SIZE_T),
55("dwTotalPageFile", SIZE_T),
56("dwAvailPageFile", SIZE_T),
57("dwTotalVirtual", SIZE_T),
58("dwAvailVirtualPhys", SIZE_T)]
59def show(self):
60for field_name, field_type in self._fields_:
61print field_name, getattr(self, field_name)
62
63memstatus = _MEMORYSTATUS()
64windll.kernel32.GlobalMemoryStatus(byref(memstatus ))
65memstatus.show()
66
67
68_______________________________
69
70from ctypes import *
71from ctypes.wintypes import *
72
73class MEMORYSTATUS(Structure):
74_fields_ = [
75('dwLength', DWORD),
76('dwMemoryLoad', DWORD),
77('dwTotalPhys', DWORD),
78('dwAvailPhys', DWORD),
79('dwTotalPageFile', DWORD),
80('dwAvailPageFile', DWORD),
81('dwTotalVirtual', DWORD),
82('dwAvailVirtual', DWORD),
83]
84
85def winmem():
86x = MEMORYSTATUS()
87windll.kernel32.GlobalMemoryStatus(byref(x))
88return x
89
90"""
91
92def mem_usage():
93    '''
94    returns the rss.
95
96  RSS  The total amount of physical memory used by  the  task,  in  kilo-
97            bytes,  is  shown  here.  For ELF processes used library pages are
98            counted here, for a.out processes not.
99           
100    Only works on nix systems.
101    '''
102    import string
103    p=os.popen('ps uwp %s'%os.getpid()) 
104    lines=p.readlines()
105    #print "lines", lines
106    status=p.close() 
107    if status or len(lines)!=2 or sys.platform == 'win32': 
108        return None 
109    return int(string.split(lines[1])[4]) 
110
111
112
113
114class BenchmarkLeastSquares:
115
116    """
117
118    Note(DSG-DSG): If you are interested in benchmarking fitting, before
119    and after blocking O:\1\dgray\before_blocking_subsandpit is before blocking
120
121    """
122
123    def __init__(self):
124        pass
125
126    def trial(self,
127              num_of_points=20000,
128              maxArea=1000,
129              max_points_per_cell=13,
130              is_fit=True,
131              use_file_type=None,
132              blocking_len=500000,
133              segments_in_mesh=True,
134              save=False,
135              verbose=False,
136              run_profile=False,
137              gridded=True,
138              geo_ref=True):
139        '''
140        num_of_points
141        '''
142        #print "num_of_points",num_of_points
143        #print "maxArea",maxArea
144        #print "max_points_per_cell", max_points_per_cell
145        if geo_ref is True:
146            geo = Geo_reference(xllcorner = 2.0, yllcorner = 2.0)
147        else:
148            geo = None
149        mesh_dict = self._build_regular_mesh_dict(maxArea=maxArea,
150                                                  is_segments=segments_in_mesh,
151                                                  save=save,
152                                                  geo=geo)
153        points_dict = self._build_points_dict(num_of_points=num_of_points,
154                                              gridded=gridded,
155                                              verbose=verbose)
156
157        #print "len(mesh_dict['triangles'])",len(mesh_dict['triangles'])
158        if is_fit is True:
159            op = "Fit_"
160        else:
161            op = "Interp_"
162        profile_file = op + "P" + str(num_of_points) + \
163                       "T" + str(len(mesh_dict['triangles'])) + \
164                       "PPC" + str(max_points_per_cell) + \
165                       ".txt"
166                   
167        # Apply the geo_ref to the points, so they are relative
168        # Pass in the geo_ref
169       
170        domain = Domain(mesh_dict['vertices'], mesh_dict['triangles'],
171                        use_cache=False, verbose=verbose,
172                                     geo_reference=geo)
173        #Initial time and memory
174        t0 = time.time()
175        #m0 = None on windows
176        m0 = mem_usage()
177       
178        # Apply the geo_ref to the points, so they are relative
179        # Pass in the geo_ref
180        geospatial = Geospatial_data(points_dict['points'],
181                                     points_dict['point_attributes'],
182                                     geo_reference=geo)
183        del points_dict
184        if is_fit is True:
185
186            # print "Fit in Fit"
187            if use_file_type == None:
188                points = geospatial
189                filename = None
190            else:
191                #FIXME (DSG) check that the type
192                fileName = tempfile.mktemp("." + use_file_type)
193                geospatial.export_points_file(fileName, absolute=True)
194                points = None
195                filename = fileName
196            if run_profile is True:
197                   
198                s = """domain.set_quantity('elevation',points,filename=filename,use_cache=False)"""
199                pobject = profile.Profile()
200                presult = pobject.runctx(s,
201                                         vars(sys.modules[__name__]),
202                                         vars())
203                prof_file = tempfile.mktemp(".prof")
204                presult.dump_stats(prof_file)
205                #
206                # Let process these results
207                S = pstats.Stats(prof_file)
208                saveout = sys.stdout
209                pfile = open(profile_file, "w")
210                sys.stdout = pfile
211                s = S.sort_stats('cumulative').print_stats(60)
212                sys.stdout = saveout
213                pfile.close()
214                os.remove(prof_file)
215            else:
216                domain.set_quantity('elevation',points,filename=filename,
217                                    use_cache=False, verbose=verbose)
218            if not use_file_type == None:
219                os.remove(fileName)
220                   
221        else:
222            # run an interploate problem.
223            #print "Interpolate!"
224           
225            if run_profile:
226                # pass in the geospatial points
227                # and the mesh origin
228                 
229                s="""benchmark_interpolate(mesh_dict['vertices'],mesh_dict['vertex_attributes'],mesh_dict['triangles'],geospatial,max_points_per_cell=max_points_per_cell,mesh_origin=geo)"""
230                pobject = profile.Profile()
231                presult = pobject.runctx(s,
232                                         vars(sys.modules[__name__]),
233                                         vars())
234                prof_file = tempfile.mktemp(".prof")
235                presult.dump_stats(prof_file)
236                #
237                # Let process these results
238                S = pstats.Stats(prof_file)
239                saveout = sys.stdout
240                pfile = open(profile_file, "w")
241                sys.stdout = pfile
242                s = S.sort_stats('cumulative').print_stats(60)
243                sys.stdout = saveout
244                pfile.close()
245                os.remove(prof_file)
246                   
247            else:
248                # pass in the geospatial points
249                 benchmark_interpolate(mesh_dict['vertices'],
250                                       mesh_dict['vertex_attributes'],
251                                       mesh_dict['triangles'],
252                                       geospatial,
253                                       mesh_origin=geo,
254                                       max_points_per_cell=max_points_per_cell,
255                                       verbose=verbose)
256        time_taken_sec = (time.time()-t0)
257        m1 = mem_usage()
258        if m0 is None or m1 is None:
259            memory_used = None
260        else:
261            memory_used = (m1 - m0)
262        #print 'That took %.2f seconds' %time_taken_sec
263
264        # return the times spent in first cell searching and
265        # backing up.
266       
267        search_one_cell_time, search_more_cells_time = search_times()
268        reset_search_times()
269        #print "bench - search_one_cell_time",search_one_cell_time
270        #print "bench - search_more_cells_time", search_more_cells_time
271        #print "bench - build_quadtree_time", get_build_quadtree_time()
272        return time_taken_sec, memory_used, len(mesh_dict['triangles']), \
273               search_one_cell_time, search_more_cells_time, \
274               get_build_quadtree_time()
275   
276
277    def _build_regular_mesh_dict(self,
278                                 maxArea=1000,
279                                 is_segments=True,
280                                 save=False,
281                                 geo=None):
282      # make a normalised mesh
283        # pretty regular size, with some segments thrown in.
284
285        # don't pass in the geo ref.
286        # it is applied in domain
287        m = Mesh() #geo_reference=geo)
288        m.addUserVertex(0,0)
289        m.addUserVertex(1.0,0)
290        m.addUserVertex(0,1.0)
291        m.addUserVertex(1.0,1.0)
292       
293        m.auto_segment(alpha = 100 )
294
295        if is_segments:
296            dict = {}
297            dict['points'] = [[.10,.10],[.90,.20]]
298            dict['segments'] = [[0,1]] 
299            dict['segment_tags'] = ['wall1']   
300            m.addVertsSegs(dict)
301   
302            dict = {}
303            dict['points'] = [[.10,.90],[.40,.20]]
304            dict['segments'] = [[0,1]] 
305            dict['segment_tags'] = ['wall2']   
306            m.addVertsSegs(dict)
307       
308            dict = {}
309            dict['points'] = [[.20,.90],[.60,.60]]
310            dict['segments'] = [[0,1]] 
311            dict['segment_tags'] = ['wall3'] 
312            m.addVertsSegs(dict)
313       
314            dict = {}
315            dict['points'] = [[.60,.20],[.90,.90]]
316            dict['segments'] = [[0,1]] 
317            dict['segment_tags'] = ['wall4']   
318            m.addVertsSegs(dict)
319
320        m.generateMesh(mode = "Q", maxArea = maxArea, minAngle=20.0)       
321        if save is True:
322            m.export_mesh_file("aaaa.tsh")
323        mesh_dict =  m.Mesh2IOTriangulationDict()
324
325        #Add vert attribute info to the mesh
326        mesh_dict['vertex_attributes'] = []
327        # There has to be a better way of doing this..
328        for vertex in mesh_dict['vertices']:
329            mesh_dict['vertex_attributes'].append([10.0])
330
331        return mesh_dict
332
333    def _build_points_dict(self, num_of_points=20000
334                           , gridded=True, verbose=False):
335       
336        points_dict = {}
337        points = []
338        point_atts = []
339
340        if gridded is True:
341            grid = int(sqrt(num_of_points))
342       
343        for point in range(num_of_points):
344            if gridded is True:
345
346                # point starts at 0.0
347                # the 2 and 0.25 is to make sure all points are in the
348                # range 0 - 1
349                points.append([float(point/grid)/float(grid*1.1)+0.0454,
350                               float(point%grid)/float(grid*1.1)+0.0454])
351            else:
352                points.append([random(), random()])
353            point_atts.append(10.0)
354
355        points_dict['points'] = points
356        points_dict['point_attributes'] = point_atts
357       
358        for point in points:
359            assert point[0] < 1.0
360            assert point[1] < 1.0
361            assert point[0] > 0.0
362            assert point[1] > 0.0
363           
364        if verbose is True:
365            pass
366            #print "points", points
367            #import sys; sys.exit()
368           
369                                   
370           
371        return points_dict
372
373
374#-------------------------------------------------------------
375if __name__ == "__main__":
376        b._build_points_dict()
Note: See TracBrowser for help on using the repository browser.