source: anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py @ 4339

"""Verify that simulation produced by ANUGA compares to published
validation timeseries ch5, ch7 and ch9 as well as the boundary timeseries.

RMS norm is printed and plots are produced.
"""

from Numeric import allclose, argmin, argmax
from Scientific.IO.NetCDF import NetCDFFile

from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.utilities.numerical_tools import\
     ensure_numeric, cov, get_machine_precision

import project

try:
    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
except:
    print 'Could not import pylab'
    plotting = False
else:
    plotting = True

# No plotting for automated unittest (unless one wishes to run
# this manually and look at the outputs)
plotting = False


#-------------------------
# Basic data
#-------------------------

finaltime = 22.5
timestep = 0.05

gauge_locations = [[0.000, 1.696]] # Boundary gauge
gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']

validation_data = {}
for key in gauge_names:
    validation_data[key] = []


# Expected values
expected_covariance = {'Boundary': 5.269569575007607815e-05, 
                       'ch5': 1.166277999581819919e-04,
                       'ch7': 1.127136457890861503e-04,
                       'ch9': 1.250659477418482129e-04}

expected_difference = {'Boundary': 8.350712673810733924e-04,
                       'ch5': 3.405426180525532483e-03,
                       'ch7': 2.852870417368218517e-03,
                       'ch9': 3.248778982037564891e-03}

expected_maximum = {'Boundary': 1.611749508386188523e-02,  
                    'ch5': 3.551308418158714147e-02,
                    'ch7': 3.858418457126511908e-02,
                    'ch9': 4.317962986578308127e-02}

expected_minimum = {'Boundary': -1.164547474575844919e-02, 
                    'ch5': -8.664439185502026408e-03,
                    'ch7': -2.726335488279797541e-03,
                    'ch9': -5.977581218447349659e-03}

expected_argmax = {'Boundary': 1.255000000000000071e+01, 
                   'ch5': 1.839999999999999858e+01,
                   'ch7': 1.700000000000000000e+01,
                   'ch9': 1.685000000000000142e+01}

expected_argmin = {'Boundary': 2.064999999999999858e+01, 
                   'ch5': 1.459999999999999964e+01,
                   'ch7': 1.230000000000000071e+01,
                   'ch9': 1.315000000000000036e+01}

#-------------------------
# Read validation dataa
#-------------------------

print 'Reading', project.boundary_filename
fid = NetCDFFile(project.boundary_filename, 'r')
input_time = fid.variables['time'][:]
validation_data['Boundary'] = fid.variables['stage'][:]

reference_time = []
fid = open(project.validation_filename)
lines = fid.readlines()
fid.close()

for i, line in enumerate(lines[1:]):
    if i == len(input_time): break
    
    fields = line.split()

    reference_time.append(float(fields[0]))    # Record reference time
    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
        value = float(fields[1:][j])           # Omit time 
        validation_data[key].append(value/100) # Convert cm2m


# Checks
assert reference_time[0] == 0.0
assert reference_time[-1] == finaltime
assert allclose(reference_time, input_time)

for key in gauge_names:
    validation_data[key] = ensure_numeric(validation_data[key])

#--------------------------------------------------
# Read and interpolate model output
#--------------------------------------------------

import sys
if len(sys.argv) > 1:
    sww_filename = sys.argv[1]
else:    
    sww_filename = project.output_filename
    
f = file_function(sww_filename,
                  quantities='stage',
                  interpolation_points=gauge_locations,
                  use_cache=True,
                  verbose=True)


def report_difference(name, computed_value, reference_value, rtol, atol):

    if abs(reference_value) > 0:
        msg = '%s (expected, computed):\n  (%.18e, %.18e):\n  Relative error=%.18e'\
              %(name, reference_value, computed_value,
                abs(reference_value-computed_value)/reference_value)
        print msg
        

    msg = '  Absolute error=%.18e'\
          %(abs(reference_value-computed_value))        
    print msg

    
    #print 'Allclose:', allclose(reference_value, computed_value,
    #                            rtol=rtol, atol=atol)
    assert allclose(reference_value, computed_value, rtol=rtol, atol=atol), msg
    


#--------------------------------------------------
# Compare model output to validation data
#--------------------------------------------------


#eps = get_machine_precision()

# Windows tolerances
rtol = 1.0e-2
atol = 1.0e-3
print 'Precisions used: rtol=%e, atol=%e' %(rtol, atol)

for k, name in enumerate(gauge_names):

    sqsum = 0
    denom = 0
    model = []
    print 
    print 'Validating ' + name
    observed_timeseries = validation_data[name]
    for i, t in enumerate(reference_time):
        model.append(f(t, point_id=k)[0])

    # Covariance measure    
    res = cov(observed_timeseries, model)
    report_difference('Covariance', res, expected_covariance[name], rtol, atol)
     
    # Difference measures    
    res = sum(abs(observed_timeseries-model))/len(model)
    report_difference('Accumulated difference', res,
                      expected_difference[name], rtol, atol)    

    # Extrema
    res = max(model)
    report_difference('Maximum', res, expected_maximum[name], rtol, atol)
    
    res = min(model)
    report_difference('Minimum', res, expected_minimum[name], rtol, atol)    

    # Locations of extrema
    #i0 = argmax(observed_timeseries)
    i1 = argmax(model)
    res = reference_time[i1]
    report_difference('Location of maximum', res, expected_argmax[name], rtol, atol)    
    

    if name <> 'ch9':
        # Minimum of ch9 is very flat and hard to pinpoint
        i1 = argmin(model)
        res = reference_time[i1]
        report_difference('Location of minimum', res, expected_argmin[name],
                          rtol, atol)        


    if plotting is True:
        ion()
        hold(False)
    
        plot(reference_time, validation_data[name], 'r-',
             reference_time, model, 'k-')
        title('Gauge %s' %name)
        xlabel('time(s)')
        ylabel('stage (m)')    
        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
        savefig(name, dpi = 300)        

        raw_input('Next')