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 | |
---|
18 | import os |
---|
19 | import sys |
---|
20 | import time |
---|
21 | from random import seed, random |
---|
22 | import tempfile |
---|
23 | import profile , pstats |
---|
24 | import tempfile |
---|
25 | |
---|
26 | from anuga.fit_interpolate.interpolate import Interpolate |
---|
27 | from anuga.fit_interpolate.fit import Fit |
---|
28 | from anuga.pmesh.mesh import Mesh |
---|
29 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
30 | |
---|
31 | def 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 | |
---|
53 | class 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 | #------------------------------------------------------------- |
---|
254 | if __name__ == "__main__": |
---|
255 | b = BenchmarkLeastSquares() |
---|
256 | b._build_regular_mesh_dict() |
---|
257 | b._build_points_dict() |
---|