source: anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py @ 7317

Last change on this file since 7317 was 7317, checked in by rwilson, 15 years ago

Replaced 'print' statements with log.critical() calls.
File size: 12.2 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
24from math import sqrt
25
26from anuga.fit_interpolate.search_functions import search_times, \
27     reset_search_times
28
29from anuga.fit_interpolate.interpolate import Interpolate
30from anuga.fit_interpolate.fit import Fit
31from anuga.pmesh.mesh import Mesh
32from anuga.geospatial_data.geospatial_data import Geospatial_data
33from anuga.shallow_water import Domain
34from anuga.fit_interpolate.fit import Fit, fit_to_mesh
35from anuga.fit_interpolate.interpolate import benchmark_interpolate
36from anuga.coordinate_transforms.geo_reference import Geo_reference
37from anuga.fit_interpolate.general_fit_interpolate import \
38     get_build_quadtree_time
39
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
103    # GC stuff should be uncommented when debugging memory 'shift' problems,
104    # when (if?) the amount of memory used takes a sudden leap upwards.
105    #import gc
106    import string
107
108    #gc.collect()
109    #print('Ran a garbage collection')
110    p=os.popen('ps uwp %s'%os.getpid()) 
111    lines=p.readlines()
112    status=p.close() 
113    if status or len(lines)!=2 or sys.platform == 'win32': 
114        return None 
115    return int(string.split(lines[1])[4]) 
116
117
118class BenchmarkLeastSquares:
119    """
120
121    Note(DSG-DSG): If you are interested in benchmarking fitting, before
122    and after blocking O:\1\dgray\before_blocking_subsandpit is before blocking
123
124    """
125
126    def __init__(self):
127        pass
128
129    def trial(self,
130              num_of_points=20000,
131              maxArea=1000,
132              max_points_per_cell=13,
133              is_fit=True,
134              use_file_type=None,
135              blocking_len=500000,
136              segments_in_mesh=True,
137              save=False,
138              verbose=False,
139              run_profile=False,
140              gridded=True,
141              geo_ref=True):
142        '''
143        num_of_points
144        '''
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        if is_fit is True:
158            op = "Fit_"
159        else:
160            op = "Interp_"
161        profile_file = op + "P" + str(num_of_points) + \
162                       "T" + str(len(mesh_dict['triangles'])) + \
163                       "PPC" + str(max_points_per_cell) + \
164                       ".txt"
165                   
166        # Apply the geo_ref to the points, so they are relative
167        # Pass in the geo_ref
168       
169        domain = Domain(mesh_dict['vertices'], mesh_dict['triangles'],
170                        use_cache=False, verbose=verbose,
171                                     geo_reference=geo)
172        #Initial time and memory
173        t0 = time.time()
174        #m0 = None on windows
175        m0 = mem_usage()
176       
177        # Apply the geo_ref to the points, so they are relative
178        # Pass in the geo_ref
179        geospatial = Geospatial_data(points_dict['points'],
180                                     points_dict['point_attributes'],
181                                     geo_reference=geo)
182        del points_dict
183        if is_fit is True:
184
185            if use_file_type == None:
186                points = geospatial
187                filename = None
188            else:
189                #FIXME (DSG) check that the type
190                fileName = tempfile.mktemp("." + use_file_type)
191                geospatial.export_points_file(fileName, absolute=True)
192                points = None
193                filename = fileName
194            if run_profile is True:
195                   
196                s = """domain.set_quantity('elevation',points,filename=filename,use_cache=False)"""
197                pobject = profile.Profile()
198                presult = pobject.runctx(s,
199                                         vars(sys.modules[__name__]),
200                                         vars())
201                prof_file = tempfile.mktemp(".prof")
202                presult.dump_stats(prof_file)
203                #
204                # Let process these results
205                S = pstats.Stats(prof_file)
206                saveout = sys.stdout
207                pfile = open(profile_file, "w")
208                sys.stdout = pfile
209                s = S.sort_stats('cumulative').print_stats(60)
210                sys.stdout = saveout
211                pfile.close()
212                os.remove(prof_file)
213            else:
214                domain.set_quantity('elevation',points,filename=filename,
215                                    use_cache=False, verbose=verbose)
216            if not use_file_type == None:
217                os.remove(fileName)
218                   
219        else:
220            # run an interploate problem.
221           
222            if run_profile:
223                # pass in the geospatial points
224                # and the mesh origin
225                 
226                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)"""
227                pobject = profile.Profile()
228                presult = pobject.runctx(s,
229                                         vars(sys.modules[__name__]),
230                                         vars())
231                prof_file = tempfile.mktemp(".prof")
232                presult.dump_stats(prof_file)
233                #
234                # Let process these results
235                S = pstats.Stats(prof_file)
236                saveout = sys.stdout
237                pfile = open(profile_file, "w")
238                sys.stdout = pfile
239                s = S.sort_stats('cumulative').print_stats(60)
240                sys.stdout = saveout
241                pfile.close()
242                os.remove(prof_file)
243                   
244            else:
245                # pass in the geospatial points
246                 benchmark_interpolate(mesh_dict['vertices'],
247                                       mesh_dict['vertex_attributes'],
248                                       mesh_dict['triangles'],
249                                       geospatial,
250                                       mesh_origin=geo,
251                                       max_points_per_cell=max_points_per_cell,
252                                       verbose=verbose)
253        time_taken_sec = (time.time()-t0)
254        m1 = mem_usage()
255        if m0 is None or m1 is None:
256            memory_used = None
257        else:
258            memory_used = (m1 - m0)
259        #print 'That took %.2f seconds' %time_taken_sec
260
261        # return the times spent in first cell searching and
262        # backing up.
263       
264        search_one_cell_time, search_more_cells_time = search_times()
265        reset_search_times()
266        #print "bench - search_one_cell_time",search_one_cell_time
267        #print "bench - search_more_cells_time", search_more_cells_time
268        #print "bench - build_quadtree_time", get_build_quadtree_time()
269        return time_taken_sec, memory_used, len(mesh_dict['triangles']), \
270               search_one_cell_time, search_more_cells_time, \
271               get_build_quadtree_time()
272   
273
274    def _build_regular_mesh_dict(self,
275                                 maxArea=1000,
276                                 is_segments=True,
277                                 save=False,
278                                 geo=None):
279      # make a normalised mesh
280        # pretty regular size, with some segments thrown in.
281
282        # don't pass in the geo ref.
283        # it is applied in domain
284        m = Mesh() #geo_reference=geo)
285        m.addUserVertex(0,0)
286        m.addUserVertex(1.0,0)
287        m.addUserVertex(0,1.0)
288        m.addUserVertex(1.0,1.0)
289       
290        m.auto_segment(alpha = 100 )
291
292        if is_segments:
293            dict = {}
294            dict['points'] = [[.10,.10],[.90,.20]]
295            dict['segments'] = [[0,1]] 
296            dict['segment_tags'] = ['wall1']   
297            m.addVertsSegs(dict)
298   
299            dict = {}
300            dict['points'] = [[.10,.90],[.40,.20]]
301            dict['segments'] = [[0,1]] 
302            dict['segment_tags'] = ['wall2']   
303            m.addVertsSegs(dict)
304       
305            dict = {}
306            dict['points'] = [[.20,.90],[.60,.60]]
307            dict['segments'] = [[0,1]] 
308            dict['segment_tags'] = ['wall3'] 
309            m.addVertsSegs(dict)
310       
311            dict = {}
312            dict['points'] = [[.60,.20],[.90,.90]]
313            dict['segments'] = [[0,1]] 
314            dict['segment_tags'] = ['wall4']   
315            m.addVertsSegs(dict)
316
317        m.generateMesh(mode = "Q", maxArea = maxArea, minAngle=20.0)       
318        if save is True:
319            m.export_mesh_file("aaaa.tsh")
320        mesh_dict =  m.Mesh2IOTriangulationDict()
321
322        #Add vert attribute info to the mesh
323        mesh_dict['vertex_attributes'] = []
324        # There has to be a better way of doing this..
325        for vertex in mesh_dict['vertices']:
326            mesh_dict['vertex_attributes'].append([10.0])
327
328        return mesh_dict
329
330    def _build_points_dict(self, num_of_points=20000
331                           , gridded=True, verbose=False):
332       
333        points_dict = {}
334        points = []
335        point_atts = []
336
337        if gridded is True:
338            grid = int(sqrt(num_of_points))
339       
340        for point in range(num_of_points):
341            if gridded is True:
342
343                # point starts at 0.0
344                # the 2 and 0.25 is to make sure all points are in the
345                # range 0 - 1
346                points.append([float(point/grid)/float(grid*1.1)+0.0454,
347                               float(point%grid)/float(grid*1.1)+0.0454])
348            else:
349                points.append([random(), random()])
350            point_atts.append(10.0)
351
352        points_dict['points'] = points
353        points_dict['point_attributes'] = point_atts
354       
355        for point in points:
356            assert point[0] < 1.0
357            assert point[1] < 1.0
358            assert point[0] > 0.0
359            assert point[1] > 0.0
360           
361        if verbose is True:
362            pass
363            #print "points", points
364            #import sys; sys.exit()
365           
366                                   
367           
368        return points_dict
369
370
371#-------------------------------------------------------------
372if __name__ == "__main__":
373        b._build_points_dict()
Note: See TracBrowser for help on using the repository browser.