source: anuga_work/development/friction_UA_flume_2006/calc_norm.py @ 4351

"""
Idea's
have it use cache!
"""


#----------------------------------------------------------------------------
# Import necessary modules
#----------------------------------------------------------------------------

# Standard modules
import os
from os import sep, path
import csv

# Related major packages
import project              
from anuga.abstract_2d_finite_volumes.util import file_function
from numerical_tools import  err   # norm, corr,  


def load_csv(file_name):
    """
    Load in the csv as a dic, title as key and column info as value, .
    Also, create a dic, title as key and column index as value,
    to keep track of the column order. 
    """
    #
    attribute_dic = {}
    title_index_dic = {}
    titles_stripped = [] # list of titles
    
    reader = csv.reader(file(file_name))

    # Read in and manipulate the title info
    titles = reader.next()
    for i,title in enumerate(titles):
        titles_stripped.append(title.strip())
        title_index_dic[title.strip()] = i
    title_count = len(titles_stripped)       
    #print "title_index_dic",title_index_dic

    
    #create a dic of colum values, indexed by column title
    for line in reader:
        if len(line) <> title_count:
            raise IOError #FIXME make this nicer
        for i, value in enumerate(line):
            attribute_dic.setdefault(titles_stripped[i],[]).append(value)
        
    return attribute_dic, title_index_dic


def write_norm_results( file_name, norms):
    fd = open(file_name, 'w')
    fd.write( 'sww_file, friction, 2nd normal' + '\n')
    for norm in norms:
        fd.write(norm[0] + ', ' + str(norm[1]) + ', ' + str(norm[2]) + '\n')
    fd.close()    
    
def calc_norm(flume_name,
              is_trial_run=False,
              outputdir_name=None,
              pro_instance=None,
              max_time=None):    
    if is_trial_run is True:
        outputdir_name += '_test'
    if pro_instance is None:
        pro_instance = project.Project(['data','flumes','dam_2007'],
                                       outputdir_name=outputdir_name)
    
    
    #  Read in the flume data
    flume_results = path.join(pro_instance.home, 'data', 'flumes',
                              'dam_UA_data_2006',
                              flume_name)
    flume_results, _ = load_csv(flume_results)
    #print "flume_results",flume_results #.keys()
    if max_time is not None:
        for i,time in enumerate(flume_results['time']):
            if float(time) > float(max_time):
                break
        for keys in flume_results.iterkeys():
            flume_results[keys] = flume_results[keys][0:i]
    
    #print "flume_results['time']", flume_results['time']
    #Get a list of the files
    files = os.listdir(pro_instance.outputdir)
    files = [x for x in files if x[-4:] == '.sww']
    files.sort()
    print "checking against files: ",files 

    quantities = ['stage', 'elevation']
    gauge_locations = [
        [5.1,0.225],  #0.5m from SWL
        [6.6,0.225],  #2m from SWL
        [6.95,0.255], #2.35m from SWL
        [7.6,0.255],  #3m from SWL
        [8.1,0.255],  #3.5m from SWL
        [9.1,0.255]  #4.5m from SWL
        ]

    norm_sums = []
    for sww_file in files:
        #print "sww_file",sww_file 
        file_instance = file_function(pro_instance.outputdir + sep + sww_file,
                                      quantities = quantities,
                                      interpolation_points = gauge_locations,
                                      verbose = True,
                                      use_cache = True)
        # file_instance(0.02, point_id=0) #This works
        norm_sum = 0
        for location in range(len(gauge_locations)):
            #print "location",location
            #print " gauge_locations[location]", gauge_locations[location]
            simulated_depths = []
            #print "flume_results['time']",flume_results['time'] 
            for time in flume_results['time']:
                quantities_slice = file_instance(float(time),
                                                 point_id=location)
                #print "quantities", quantities_slice 
                depth = quantities_slice[0] - quantities_slice[1]
                simulated_depths.append(depth)
            #print "simulated_depths",simulated_depths 
            # get the flume depth at this location
            flume_depths = flume_results[str(gauge_locations[location][0])]
            flume_depths = [float(i) for i in flume_depths]
            # calc the norm
            #print "map(None, simulated_depths, flume_depths)", \
            #      map(None, simulated_depths, flume_depths)
            norm = err(simulated_depths,
                       flume_depths, 2, relative = True)  # 2nd norm (rel. RMS)
            norm_sum += norm
        print "norm_sum",norm_sum
        norm_sums.append([sww_file, sww_file[2:-4], norm_sum])
        norm_file_name = pro_instance.outputdir + sep + \
                           flume_name[:-4]+'_norm' + flume_name[-4:]
        write_norm_results(norm_file_name, norm_sums)

#-------------------------------------------------------------
if __name__ == "__main__":
    calc_norm(
        'ensemble_average_h_smooth.csv',
        #'ensemble_average_h_rough.csv',
        is_trial_run = False,   # adds _test to outputdir_name
        outputdir_name='friction_UA_new_meshII',
        max_time=None)