source: trunk/anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py @ 8050

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

Cleaned up unit tests to use API.

File size: 11.9 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.shallow_water_domain 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 - build_quadtree_time", get_build_quadtree_time()
264        return time_taken_sec, memory_used, len(mesh_dict['triangles']), \
265               get_build_quadtree_time()
266   
267
268    def _build_regular_mesh_dict(self,
269                                 maxArea=1000,
270                                 is_segments=True,
271                                 save=False,
272                                 geo=None):
273      # make a normalised mesh
274        # pretty regular size, with some segments thrown in.
275
276        # don't pass in the geo ref.
277        # it is applied in domain
278        m = Mesh() #geo_reference=geo)
279        m.addUserVertex(0,0)
280        m.addUserVertex(1.0,0)
281        m.addUserVertex(0,1.0)
282        m.addUserVertex(1.0,1.0)
283       
284        m.auto_segment(alpha = 100 )
285
286        if is_segments:
287            dict = {}
288            dict['points'] = [[.10,.10],[.90,.20]]
289            dict['segments'] = [[0,1]] 
290            dict['segment_tags'] = ['wall1']   
291            m.addVertsSegs(dict)
292   
293            dict = {}
294            dict['points'] = [[.10,.90],[.40,.20]]
295            dict['segments'] = [[0,1]] 
296            dict['segment_tags'] = ['wall2']   
297            m.addVertsSegs(dict)
298       
299            dict = {}
300            dict['points'] = [[.20,.90],[.60,.60]]
301            dict['segments'] = [[0,1]] 
302            dict['segment_tags'] = ['wall3'] 
303            m.addVertsSegs(dict)
304       
305            dict = {}
306            dict['points'] = [[.60,.20],[.90,.90]]
307            dict['segments'] = [[0,1]] 
308            dict['segment_tags'] = ['wall4']   
309            m.addVertsSegs(dict)
310
311        m.generateMesh(mode = "Q", maxArea = maxArea, minAngle=20.0)       
312        if save is True:
313            m.export_mesh_file("aaaa.tsh")
314        mesh_dict =  m.Mesh2IOTriangulationDict()
315
316        #Add vert attribute info to the mesh
317        mesh_dict['vertex_attributes'] = []
318        # There has to be a better way of doing this..
319        for vertex in mesh_dict['vertices']:
320            mesh_dict['vertex_attributes'].append([10.0])
321
322        return mesh_dict
323
324    def _build_points_dict(self, num_of_points=20000
325                           , gridded=True, verbose=False):
326       
327        points_dict = {}
328        points = []
329        point_atts = []
330
331        if gridded is True:
332            grid = int(sqrt(num_of_points))
333       
334        for point in range(num_of_points):
335            if gridded is True:
336
337                # point starts at 0.0
338                # the 2 and 0.25 is to make sure all points are in the
339                # range 0 - 1
340                points.append([float(point/grid)/float(grid*1.1)+0.0454,
341                               float(point%grid)/float(grid*1.1)+0.0454])
342            else:
343                points.append([random(), random()])
344            point_atts.append(10.0)
345
346        points_dict['points'] = points
347        points_dict['point_attributes'] = point_atts
348       
349        for point in points:
350            assert point[0] < 1.0
351            assert point[1] < 1.0
352            assert point[0] > 0.0
353            assert point[1] > 0.0
354           
355        if verbose is True:
356            pass
357            #print "points", points
358            #import sys; sys.exit()
359           
360                                   
361           
362        return points_dict
363
364
365#-------------------------------------------------------------
366if __name__ == "__main__":
367        b._build_points_dict()
Note: See TracBrowser for help on using the repository browser.