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

Last change on this file since 4839 was 4839, checked in by duncan, 16 years ago

speeding up slow functions in general mesh. Fix for ticket 201. Also checking in GA validation test for profiling fitting

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