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