| 1 | """ |
|---|
| 2 | Idea's |
|---|
| 3 | have it use cache! |
|---|
| 4 | """ |
|---|
| 5 | |
|---|
| 6 | |
|---|
| 7 | #---------------------------------------------------------------------------- |
|---|
| 8 | # Import necessary modules |
|---|
| 9 | #---------------------------------------------------------------------------- |
|---|
| 10 | |
|---|
| 11 | # Standard modules |
|---|
| 12 | import os |
|---|
| 13 | from os import sep, path |
|---|
| 14 | import csv |
|---|
| 15 | |
|---|
| 16 | # Related major packages |
|---|
| 17 | import project |
|---|
| 18 | from anuga.abstract_2d_finite_volumes.util import file_function |
|---|
| 19 | from numerical_tools import err # norm, corr, |
|---|
| 20 | |
|---|
| 21 | |
|---|
| 22 | def load_csv(file_name): |
|---|
| 23 | """ |
|---|
| 24 | Load in the csv as a dic, title as key and column info as value, . |
|---|
| 25 | Also, create a dic, title as key and column index as value, |
|---|
| 26 | to keep track of the column order. |
|---|
| 27 | """ |
|---|
| 28 | # |
|---|
| 29 | attribute_dic = {} |
|---|
| 30 | title_index_dic = {} |
|---|
| 31 | titles_stripped = [] # list of titles |
|---|
| 32 | |
|---|
| 33 | reader = csv.reader(file(file_name)) |
|---|
| 34 | |
|---|
| 35 | # Read in and manipulate the title info |
|---|
| 36 | titles = reader.next() |
|---|
| 37 | for i,title in enumerate(titles): |
|---|
| 38 | titles_stripped.append(title.strip()) |
|---|
| 39 | title_index_dic[title.strip()] = i |
|---|
| 40 | title_count = len(titles_stripped) |
|---|
| 41 | #print "title_index_dic",title_index_dic |
|---|
| 42 | |
|---|
| 43 | |
|---|
| 44 | #create a dic of colum values, indexed by column title |
|---|
| 45 | for line in reader: |
|---|
| 46 | if len(line) <> title_count: |
|---|
| 47 | raise IOError #FIXME make this nicer |
|---|
| 48 | for i, value in enumerate(line): |
|---|
| 49 | attribute_dic.setdefault(titles_stripped[i],[]).append(value) |
|---|
| 50 | |
|---|
| 51 | return attribute_dic, title_index_dic |
|---|
| 52 | |
|---|
| 53 | |
|---|
| 54 | def write_norm_results( file_name, norms): |
|---|
| 55 | fd = open(file_name, 'w') |
|---|
| 56 | fd.write( 'sww_file, friction, 2nd normal' + '\n') |
|---|
| 57 | for norm in norms: |
|---|
| 58 | fd.write(norm[0] + ', ' + str(norm[1]) + ', ' + str(norm[2]) + '\n') |
|---|
| 59 | fd.close() |
|---|
| 60 | |
|---|
| 61 | def calc_norm(flume_name, |
|---|
| 62 | is_trial_run=False, |
|---|
| 63 | outputdir_name=None, |
|---|
| 64 | pro_instance=None): |
|---|
| 65 | if is_trial_run is True: |
|---|
| 66 | outputdir_name += '_test' |
|---|
| 67 | if pro_instance is None: |
|---|
| 68 | pro_instance = project.Project(['data','flumes','dam_2006'], |
|---|
| 69 | outputdir_name=outputdir_name) |
|---|
| 70 | |
|---|
| 71 | |
|---|
| 72 | # Read in the flume data |
|---|
| 73 | flume_results = path.join(pro_instance.home, 'data', 'flumes', |
|---|
| 74 | 'dam_UA_data_2006', |
|---|
| 75 | flume_name) |
|---|
| 76 | flume_results, _ = load_csv(flume_results) |
|---|
| 77 | #print "flume_results",flume_results #.keys() |
|---|
| 78 | |
|---|
| 79 | #print "flume_results['time']", flume_results['time'] |
|---|
| 80 | #Get a list of the files |
|---|
| 81 | files = os.listdir(pro_instance.outputdir) |
|---|
| 82 | files = [x for x in files if x[-4:] == '.sww'] |
|---|
| 83 | files.sort() |
|---|
| 84 | print "checking against files: ",files |
|---|
| 85 | |
|---|
| 86 | quantities = ['stage', 'elevation'] |
|---|
| 87 | gauge_locations = [ |
|---|
| 88 | [5.1,0.225], #0.5m from SWL |
|---|
| 89 | [6.6,0.225], #2m from SWL |
|---|
| 90 | [6.95,0.255], #2.35m from SWL |
|---|
| 91 | [7.6,0.255], #3m from SWL |
|---|
| 92 | [8.1,0.255], #3.5m from SWL |
|---|
| 93 | [9.1,0.255] #4.5m from SWL |
|---|
| 94 | ] |
|---|
| 95 | |
|---|
| 96 | norm_sums = [] |
|---|
| 97 | for sww_file in files: |
|---|
| 98 | #print "sww_file",sww_file |
|---|
| 99 | file_instance = file_function(pro_instance.outputdir + sep + sww_file, |
|---|
| 100 | quantities = quantities, |
|---|
| 101 | interpolation_points = gauge_locations, |
|---|
| 102 | verbose = True, |
|---|
| 103 | use_cache = True) |
|---|
| 104 | # file_instance(0.02, point_id=0) #This works |
|---|
| 105 | norm_sum = 0 |
|---|
| 106 | for location in range(len(gauge_locations)): |
|---|
| 107 | #print "location",location |
|---|
| 108 | #print " gauge_locations[location]", gauge_locations[location] |
|---|
| 109 | simulated_depths = [] |
|---|
| 110 | #print "flume_results['time']",flume_results['time'] |
|---|
| 111 | for time in flume_results['time']: |
|---|
| 112 | quantities_slice = file_instance(float(time), |
|---|
| 113 | point_id=location) |
|---|
| 114 | #print "quantities", quantities_slice |
|---|
| 115 | depth = quantities_slice[0] - quantities_slice[1] |
|---|
| 116 | simulated_depths.append(depth) |
|---|
| 117 | #print "simulated_depths",simulated_depths |
|---|
| 118 | # get the flume depth at this location |
|---|
| 119 | flume_depths = flume_results[str(gauge_locations[location][0])] |
|---|
| 120 | flume_depths = [float(i) for i in flume_depths] |
|---|
| 121 | # calc the norm |
|---|
| 122 | #print "map(None, simulated_depths, flume_depths)", \ |
|---|
| 123 | # map(None, simulated_depths, flume_depths) |
|---|
| 124 | norm = err(simulated_depths, |
|---|
| 125 | flume_depths, 2, relative = True) # 2nd norm (rel. RMS) |
|---|
| 126 | norm_sum += norm |
|---|
| 127 | print "norm_sum",norm_sum |
|---|
| 128 | norm_sums.append([sww_file, sww_file[2:-4], norm_sum]) |
|---|
| 129 | norm_file_name = pro_instance.outputdir + sep + \ |
|---|
| 130 | flume_name[:-4]+'_norm' + flume_name[-4:] |
|---|
| 131 | write_norm_results(norm_file_name, norm_sums) |
|---|
| 132 | |
|---|
| 133 | #------------------------------------------------------------- |
|---|
| 134 | if __name__ == "__main__": |
|---|
| 135 | calc_norm( |
|---|
| 136 | #'ensemble_average_h_smooth.csv', |
|---|
| 137 | 'ensemble_average_h_rough.csv', |
|---|
| 138 | is_trial_run = False, |
|---|
| 139 | outputdir_name='friction_UA') |
|---|
| 140 | |
|---|
