source: anuga_validation/Hinwood_2008/calc_norm.py @ 5689

"""
This is getting really messy.

All err results are going into the same dir, and it can't really be changed.
"""


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

# Standard modules
import os
from os import sep, path
from csv import writer
from time import localtime, strftime

# Related major packages
from Numeric import arange, array, zeros, Float, where, greater, less, \
     compress, argmin, choose, searchsorted, sqrt, sum

import project          
from os import sep
from anuga.shallow_water.data_manager import csv2dict
from anuga.utilities.numerical_tools import ensure_numeric, err, norm

from anuga.utilities.interp import interp

def get_max_min_condition_array(min, max, vector):
    
    SMALL_MIN = -1e10  # Not that small, but small enough
    vector = ensure_numeric(vector)
    assert min > SMALL_MIN
    no_maxs = where(less(vector,max), vector, SMALL_MIN)
    #print "no_maxs", no_maxs
    band_condition = greater(no_maxs, min)
    return band_condition

    
def auto_rrms(outputdir_tag, scenarios, quantity='stage',
              y_location_tag=':0.0'):
    """
    Given a bunch of scenarios that have CSV guage files, calc the
    err, Number_of_samples and rmsd for all gauges in each scenario.
    Write this info to a file for each scenario.
    """
    for run_data in scenarios:           
        location_sims = []
        location_exps = []
        for gauge_x in run_data['gauge_x']:
            gauge_x = str(gauge_x)
            location_sims.append(gauge_x + y_location_tag)
            location_exps.append(gauge_x)
        
        id = run_data['scenario_id']
        outputdir_name = id + outputdir_tag
        file_sim = outputdir_name + '_' +  quantity + ".csv"
        file_exp = id + '_exp_' + quantity + '.csv'
        file_err = outputdir_name + "_" + quantity + "_err.csv"
        

        simulation, _ = csv2dict(file_sim)
        experiment, _ = csv2dict(file_exp)
        
        time_sim = [float(x) for x in simulation['time']]
        time_exp = [float(x) for x in experiment['Time']]
        time_sim = ensure_numeric(time_sim)
        time_exp = ensure_numeric(time_exp)
        #print "min(time_exp)", min(time_exp)
        #print "max(time_exp)", max(time_exp)
        
        condition = get_max_min_condition_array(run_data['wave_times'][0],
                                                run_data['wave_times'][1],
                                                time_exp)
        time_exp_cut = compress(condition, time_exp) #, axis=axis)
        #print "min(time_exp_cut)", min(time_exp_cut)
        #print "max(time_exp_cut)", max(time_exp_cut)
        
        #assert min(time_sim) < min(time_exp)
        
        print "Writing to ", file_err
        
        err_list = []
        points = []
        rmsd_list = []
        for location_sim, location_exp in map(None, location_sims,
                                              location_exps):
            quantity_sim = [float(x) for x in simulation[location_sim]]
            quantity_exp = [float(x) for x in experiment[location_exp]]

            quantity_exp_cut = compress(condition, quantity_exp)

            # Now let's do interpolation
            quantity_sim_interp = interp(quantity_sim, time_sim, time_exp_cut)

            assert len(quantity_sim_interp) == len(quantity_exp_cut)
            norm = err(quantity_sim_interp,
                       quantity_exp_cut,
                       2, relative = False)  # 2nd norm (rel. RMS)
            err_list.append(norm)
            points.append(len(quantity_sim_interp))
            rmsd_list.append(norm/sqrt(len(quantity_sim_interp)))              
            #print "norm", norm
            #for i in range(len(quantity_sim_interp)):
                
                #print "quantity_sim_interp", quantity_sim_interp[i]
                #print "quantity_exp_cut", quantity_exp_cut[i]
        assert len(location_exps) == len(err_list)

        # Writing the file out for one scenario
        a_writer = writer(file(file_err, "wb"))
        a_writer.writerow(["x location", "err", "Number_of_samples", "rmsd"])
        a_writer.writerows(map(None,
                               location_exps,
                               err_list,
                               points,
                               rmsd_list))



def load_sensors(quantity_file):
    """
    Load a csv file, where the first row is the column header and
    the first colum explains the rows.
    """
    #slope, _ = csv2dict(file_sim)
    
    # Read the depth file
    dfid = open(quantity_file)
    lines = dfid.readlines()
    dfid.close()

    title = lines.pop(0)
    n_time = len(lines)
    n_sensors = len(lines[0].split(','))-1  # -1 to remove time
    times = zeros(n_time, Float)  #Time
    depths = zeros(n_time, Float)  #
    sensors = zeros((n_time,n_sensors), Float)
    quantity_locations = title.split(',') #(',')
    quantity_locations.pop(0) # remove 'time'

    # Doing j.split(':')[0] drops the y location
    locations = [float(j.split(':')[0]) for j in quantity_locations]
    
    for i, line in enumerate(lines):
        fields = line.split(',') #(',')
        fields = [float(j) for j in fields]
        times[i] = fields[0]
        sensors[i] = fields[1:] # 1: to remove time

    #print "times",times
    #print "locations", locations
    #print "sensors", sensors
    return times, locations, sensors                  
    
        
        

# Return a bunch of lists
# The err files, for all scenarios
def err_files(scenarios, outputdir_tag, quantity='stage'):
    """
    The err files, for a list of scenarios
    """
    file_errs = []
    for scenario in scenarios:
        id = scenario['scenario_id']
        outputdir_name = id + outputdir_tag
        file_err =  outputdir_name + "_" + quantity + "_err.csv"
        file_errs.append(file_err)
    return file_errs
    

def compare_different_settings(outputdir_tag, scenarios, quantity='stage'):

    files = err_files(scenarios, outputdir_tag, quantity=quantity)
    err = 0.0
    number_of_samples = 0
    for run_data, file in map(None, scenarios, files):
        
        simulation, _ = csv2dict(file)
        err_list = [float(x) for x in simulation['err']]
        number_of_samples_list = [float(x) for x in \
                                  simulation['Number_of_samples']]
        
        if number_of_samples is not 0:
            err_list.append(err)
            number_of_samples_list.append(number_of_samples)
        err, number_of_samples = err_addition(err_list, number_of_samples_list)
    rmsd = err/sqrt(number_of_samples)
    print outputdir_tag + "   " + str(rmsd)
    
    
    
def err_addition(err_list, number_of_samples_list):
    """
    err1 is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values
    number_of_samples1 is the number of values associated with the err.
    
    If this function gets used alot, maybe pull this out and make it an object
    """
    err = norm(ensure_numeric(err_list))
    number_of_samples = sum(ensure_numeric(number_of_samples_list))

    return err, number_of_samples

                 
#-------------------------------------------------------------
if __name__ == "__main__":
    """ 
    """
    from scenarios import scenarios
    
    
    scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!

    
    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.5_mta_0.01"
    calc_norms = True
    #calc_norms = False
    if calc_norms:
        auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
    compare_different_settings(outputdir_tag, scenarios, "stage")