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

Last change on this file since 7718 was 7718, checked in by hudson, 14 years ago

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