1 | """Least squares smoothing 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 profile , pstats |
---|
23 | import tempfile |
---|
24 | import gc |
---|
25 | |
---|
26 | from Scientific.IO.NetCDF import NetCDFFile |
---|
27 | from Numeric import array, Float |
---|
28 | |
---|
29 | from anuga.fit_interpolate.interpolate import Interpolate |
---|
30 | from anuga.fit_interpolate.fit import Fit |
---|
31 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
32 | from anuga.shallow_water import Domain |
---|
33 | from anuga.fit_interpolate.fit import Fit, fit_to_mesh |
---|
34 | from anuga.fit_interpolate.interpolate import benchmark_interpolate |
---|
35 | from anuga.shallow_water.data_manager import Write_sww, export_grid |
---|
36 | |
---|
37 | def mem_usage(): |
---|
38 | ''' |
---|
39 | returns the rss. |
---|
40 | |
---|
41 | RSS The total amount of physical memory used by the task, in kilo- |
---|
42 | bytes, is shown here. For ELF processes used library pages are |
---|
43 | counted here, for a.out processes not. |
---|
44 | |
---|
45 | Only works on nix systems. |
---|
46 | ''' |
---|
47 | import string |
---|
48 | p=os.popen('ps uwp %s'%os.getpid()) |
---|
49 | lines=p.readlines() |
---|
50 | #print "lines", lines |
---|
51 | status=p.close() |
---|
52 | if status or len(lines)!=2 or sys.platform == 'win32': |
---|
53 | return None |
---|
54 | return int(string.split(lines[1])[4]) |
---|
55 | |
---|
56 | |
---|
57 | |
---|
58 | |
---|
59 | |
---|
60 | def trial(vert_rows, |
---|
61 | vert_columns, |
---|
62 | output_quantities=['elevation'], |
---|
63 | save=False, |
---|
64 | verbose=False, |
---|
65 | run_profile=False): |
---|
66 | import time |
---|
67 | |
---|
68 | sww_file_name = build_sww(vert_rows, vert_columns, save=save) |
---|
69 | |
---|
70 | #Initial time and memory |
---|
71 | t0 = time.time() |
---|
72 | #m0 = None on windows |
---|
73 | m0 = mem_usage() |
---|
74 | #print "m0", m0 |
---|
75 | |
---|
76 | fileout = export_grid(sww_file_name, quantities=output_quantities) |
---|
77 | |
---|
78 | time_taken_sec = (time.time()-t0) |
---|
79 | m1 = mem_usage() |
---|
80 | #print "m1", m1 |
---|
81 | if m0 is None or m1 is None: |
---|
82 | memory_used = None |
---|
83 | else: |
---|
84 | memory_used = (m1 - m0) |
---|
85 | return time_taken_sec, memory_used, m0, m1 |
---|
86 | |
---|
87 | def build_sww(vert_rows, vert_columns, save): |
---|
88 | """ |
---|
89 | Build an sww file we can read. |
---|
90 | verts_rows and vert_columns have to be greater that 1 |
---|
91 | """ |
---|
92 | from anuga.pmesh.mesh import Mesh |
---|
93 | m = Mesh() |
---|
94 | m.build_grid(vert_rows, vert_columns) |
---|
95 | if save is True: |
---|
96 | m.export_mesh_file("aaaa.tsh") |
---|
97 | mesh_dict = m.Mesh2IOTriangulationDict() |
---|
98 | |
---|
99 | sww_fileName = tempfile.mktemp(".sww" ) |
---|
100 | # sww_fileName = "aa.sww" |
---|
101 | elevation = array(range(len(mesh_dict["vertices"]))) |
---|
102 | stage = elevation |
---|
103 | ymomentum = elevation |
---|
104 | xmomentum = elevation |
---|
105 | |
---|
106 | # NetCDF file definition |
---|
107 | fid = NetCDFFile(sww_fileName, 'w') |
---|
108 | sww = Write_sww() |
---|
109 | sww.store_header(fid, 0, |
---|
110 | len(mesh_dict['triangles']), |
---|
111 | len(mesh_dict["vertices"]),sww_precision=Float) |
---|
112 | sww.store_triangulation(fid, |
---|
113 | mesh_dict["vertices"], mesh_dict['triangles'], |
---|
114 | elevation) |
---|
115 | |
---|
116 | for time_step in range(10): |
---|
117 | sww.store_quantities(fid, |
---|
118 | time=time_step, |
---|
119 | stage=stage, |
---|
120 | xmomentum=xmomentum, |
---|
121 | ymomentum=ymomentum) |
---|
122 | #Close |
---|
123 | fid.close() |
---|
124 | del mesh_dict |
---|
125 | del Mesh |
---|
126 | return sww_fileName |
---|
127 | |
---|
128 | #------------------------------------------------------------- |
---|
129 | if __name__ == "__main__": |
---|
130 | delimiter = ',' |
---|
131 | |
---|
132 | ofile = 'lbm_resultsII.csv' |
---|
133 | run_profile = False #True |
---|
134 | size_list = [[4,5],[50,40]] |
---|
135 | #size_list = [[5,4]] |
---|
136 | |
---|
137 | quantitiy_list = [['elevation'], |
---|
138 | ['elevation','xmomentum'], |
---|
139 | ['elevation','xmomentum','ymomentum']] |
---|
140 | |
---|
141 | fd = open(ofile,'a') |
---|
142 | # write the title line |
---|
143 | fd.write("size" + delimiter + |
---|
144 | "num_of_quantities" + delimiter + |
---|
145 | "is_profiling" + delimiter + |
---|
146 | "mem - diff" + delimiter + |
---|
147 | "mem - initial" + delimiter + |
---|
148 | "mem - final" + delimiter + |
---|
149 | "time" "\n") |
---|
150 | |
---|
151 | |
---|
152 | |
---|
153 | for size in size_list: |
---|
154 | for quantitiy in quantitiy_list: |
---|
155 | #sys.gettotalrefcount() |
---|
156 | #sys.getobjects() |
---|
157 | gc.collect() |
---|
158 | #print " gc.get_objects() ", gc.get_objects() |
---|
159 | time, mem, m0, m1 = trial(size[0] |
---|
160 | ,size[1] |
---|
161 | ,output_quantities=quantitiy |
---|
162 | ,run_profile=run_profile |
---|
163 | ) |
---|
164 | print "time",time |
---|
165 | print "mem", mem |
---|
166 | gc.collect() |
---|
167 | print "mem", mem |
---|
168 | #print "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", |
---|
169 | #print "gc.get_referrers()", gc.get_referrers() |
---|
170 | #sys.gettotalrefcount() |
---|
171 | fd.write(str(size[0]*size[1]) + delimiter + |
---|
172 | str(len(quantitiy)) + delimiter + |
---|
173 | str(run_profile) + delimiter + |
---|
174 | str(mem) + delimiter + |
---|
175 | str(m0) + delimiter + |
---|
176 | str(m1) + delimiter + |
---|
177 | str(time) + delimiter + "\n") |
---|
178 | print " gc.get_objects() ", gc.get_objects() |
---|
179 | fd.close() |
---|