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

Last change on this file since 4648 was 4648, checked in by duncan, 17 years ago

comments and updating benchmark

File size: 9.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
24import tempfile
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
30
31def mem_usage():
32    '''
33    returns the rss.
34
35  RSS  The total amount of physical memory used by  the  task,  in  kilo-
36            bytes,  is  shown  here.  For ELF processes used library pages are
37            counted here, for a.out processes not.
38           
39    Only works on nix systems.
40    '''
41    import string
42    p=os.popen('ps uwp %s'%os.getpid()) 
43    lines=p.readlines()
44    #print "lines", lines
45    status=p.close() 
46    if status or len(lines)!=2 or sys.platform == 'win32': 
47        return None 
48    return int(string.split(lines[1])[4]) 
49
50
51
52
53class BenchmarkLeastSquares:
54
55    """
56
57    Note(DSG-DSG): If you are interested in benchmarking fitting, before
58    and after blocking O:\1\dgray\before_blocking_subsandpit is before blocking
59
60    """
61
62    def __init__(self):
63        pass
64
65    def trial(self,
66              num_of_points=20000,
67              maxArea=1000,
68              max_points_per_cell=4,
69              is_fit=True,
70              use_least_squares=False,
71              use_file_type=None,
72              blocking_len=500000,
73              segments_in_mesh=True,
74              save=False,
75              verbose=False,
76              run_profile=False):
77        '''
78        num_of_points
79        '''
80        #print "num_of_points",num_of_points
81        #print "maxArea",maxArea
82        #print "max_points_per_cell", max_points_per_cell
83
84        mesh_dict = self._build_regular_mesh_dict(maxArea=maxArea,
85                                                  is_segments=segments_in_mesh,
86                                                  save=save)
87        points_dict = self._build_points_dict(num_of_points=num_of_points)
88           
89        #Initial time and memory
90        t0 = time.time()
91        #m0 = None on windows
92        m0 = mem_usage()
93       
94        profile_file = "P" + str(num_of_points) + \
95                       "T" + str(len(mesh_dict['triangles'])) + \
96                       "PPC" + str(max_points_per_cell) + \
97                       ".txt"
98                       
99        if use_least_squares is True:
100            from anuga.where.least_squares import Interpolation
101            interp = Interpolation(mesh_dict['vertices'],
102                                   mesh_dict['triangles'],
103                                   points_dict['points'],
104                                   expand_search=True,
105                                   verbose = False,
106                                   max_points_per_cell = max_points_per_cell) 
107            if is_fit is True:
108                print "Fit in least squares"
109                calc = interp.fit_points(points_dict['point_attributes'])
110               
111            else:
112                # run an interploate problem.
113                print "Interpolate!"
114                calc = interp.interpolate(mesh_dict['vertex_attributes'])
115        else: 
116            if is_fit is True:
117                from anuga.fit_interpolate.fit import Fit
118                interp = Fit(mesh_dict['vertices'],
119                                 mesh_dict['triangles'], 
120                                 max_vertices_per_cell = max_points_per_cell)
121                print "Fit in Fit"
122                if use_file_type == None:
123                    calc = interp.fit(points_dict['points'],
124                                      points_dict['point_attributes'])
125                else:
126                    #check that the type
127                    fileName = tempfile.mktemp("." + use_file_type)
128                    G1 = Geospatial_data(points_dict['points'],
129                                         points_dict['point_attributes'])
130                    G1.export_points_file(fileName, absolute=True)
131
132                    if run_profile:
133                   
134                        s = """interp.fit(fileName, verbose=verbose)"""
135                        pobject = profile.Profile()
136                        presult = pobject.runctx(s,
137                                                 vars(sys.modules[__name__]),
138                                                 vars())
139                        prof_file = tempfile.mktemp(".prof")
140                        presult.dump_stats(prof_file)
141                        #
142                        # Let process these results
143                        S = pstats.Stats(prof_file)
144                        saveout = sys.stdout
145                        pfile = open(profile_file, "w")
146                        sys.stdout = pfile
147                        s = S.sort_stats('cumulative').print_stats(30)
148                        sys.stdout = saveout
149                        pfile.close()
150                        os.remove(prof_file)
151                    else:
152                        interp.fit(fileName, verbose=verbose)
153                    os.remove(fileName)
154                   
155            else:
156                # run an interploate problem.
157                print "Interpolate!"
158               
159                interp = Interpolate(mesh_dict['vertices'],
160                                     mesh_dict['triangles'], 
161                                 max_vertices_per_cell = max_points_per_cell)
162                s = """calc = interp.interpolate(mesh_dict['vertex_attributes']
163                                          ,points_dict['points']
164                                          ,start_blocking_len=blocking_len)"""
165               
166                fileName = tempfile.mktemp(".prof")
167                profile.run(s, fileName) #profile_file)
168               
169                S = pstats.Stats(fileName)
170                s = S.sort_stats('cumulative').print_stats(30)
171                print "***********"
172                print s
173                print "***********"
174                pfile = file.open(profile_file, "w")
175                pfile.write(s)
176                pfile.close()
177        time_taken_sec = (time.time()-t0)
178        m1 = mem_usage()
179        if m0 is None or m1 is None:
180            memory_used = None
181        else:
182            memory_used = (m1 - m0)
183        #print 'That took %.2f seconds' %time_taken_sec
184        return time_taken_sec, memory_used, len(mesh_dict['triangles'])
185
186    def _build_regular_mesh_dict(self,
187                                 maxArea=1000,
188                                 is_segments=True,
189                                 save=False):
190      # make a normalised mesh
191        # pretty regular size, with some segments thrown in.
192        m = Mesh()
193        m.addUserVertex(0,0)
194        m.addUserVertex(1.0,0)
195        m.addUserVertex(0,1.0)
196        m.addUserVertex(1.0,1.0)
197       
198        m.auto_segment(alpha = 100 )
199
200        if is_segments:
201            dict = {}
202            dict['points'] = [[.10,.10],[.90,.20]]
203            dict['segments'] = [[0,1]] 
204            dict['segment_tags'] = ['wall1']   
205            m.addVertsSegs(dict)
206   
207            dict = {}
208            dict['points'] = [[.10,.90],[.40,.20]]
209            dict['segments'] = [[0,1]] 
210            dict['segment_tags'] = ['wall2']   
211            m.addVertsSegs(dict)
212       
213            dict = {}
214            dict['points'] = [[.20,.90],[.60,.60]]
215            dict['segments'] = [[0,1]] 
216            dict['segment_tags'] = ['wall3'] 
217            m.addVertsSegs(dict)
218       
219            dict = {}
220            dict['points'] = [[.60,.20],[.90,.90]]
221            dict['segments'] = [[0,1]] 
222            dict['segment_tags'] = ['wall4']   
223            m.addVertsSegs(dict)
224
225        m.generateMesh(mode = "Q", maxArea = maxArea, minAngle=20.0)       
226        if save is True:
227            m.export_mesh_file("aaaa.tsh")
228        mesh_dict =  m.Mesh2IOTriangulationDict()
229
230        #Add vert attribute info to the mesh
231        mesh_dict['vertex_attributes'] = []
232        # There has to be a better way of doing this..
233        for vertex in mesh_dict['vertices']:
234            mesh_dict['vertex_attributes'].append([10.0])
235
236        return mesh_dict
237
238    def _build_points_dict(self, num_of_points=20000):
239       
240        points_dict = {}
241        points = []
242        point_atts = []
243
244        for point in range(num_of_points):
245            points.append([random(), random()])
246            point_atts.append(10.0)
247
248        points_dict['points'] = points
249        points_dict['point_attributes'] = point_atts
250        return points_dict
251
252
253#-------------------------------------------------------------
254if __name__ == "__main__":
255        b = BenchmarkLeastSquares()
256        b._build_regular_mesh_dict()
257        b._build_points_dict()
Note: See TracBrowser for help on using the repository browser.