source: branches/numpy/anuga/fit_interpolate/benchmark_least_squares.py @ 7239

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

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