Changeset 7649


Ignore:
Timestamp:
Mar 2, 2010, 5:19:13 PM (14 years ago)
Author:
ole
Message:

Added variable tolerance and more measures to Patong auto validation

Location:
anuga_validation/automated_validation_tests/patong_beach_validation
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_validation/automated_validation_tests/patong_beach_validation/compare_model_timeseries.py

    r7645 r7649  
    3636def usage():
    3737    print 'Usage:'
    38     print 'python compare_model_timeseries.py swwfile1 swwfile2'
     38    print 'python compare_model_timeseries.py swwfile1 swwfile2 epsilon'
    3939
    4040
     
    7474                       timeseries2,
    7575                       name='',
    76                        legend = ('Time series 1', 'Time series 2')):
     76                       legend = ('Time series 1', 'Time series 2'),
     77                       eps=1.0e-6):
    7778    """Compare and plot two timeseries
    7879    """
     
    8485    assert timeseries1.shape[0] == N   
    8586    assert timeseries2.shape[0] == N       
    86    
    87     eps = 1.0e-3
    88    
     87
     88    print 'epsilon', eps
    8989    # Covariance measure   
    90     #res = cov(timeseries1, timeseries2)
    91     #msg = 'Covariance between timeseries was too large: %e' % res
    92     #assert res < 0.01, msg
     90    res = cov(timeseries1-timeseries2)
     91    print '2 norm diff', res       
     92    msg = '2-norm of timeseries difference was too large: %e' % res
     93    assert res < eps, msg
    9394     
    94     # Difference measures
    95     print timeseries1
    96     print timeseries1-timeseries2
    97    
    98    
    9995    # Maximum norm
    10096    nominator = max(abs(timeseries1-timeseries2))   
    10197    denominator = max(abs(timeseries1))
    102     print nominator, denominator, N
    10398    if denominator > 0:
    10499        # Relative measure
     
    107102        # Absolute measure
    108103        res = nominator/N
    109                
     104
     105    print 'Max abs diff', res   
    110106    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
    111     #assert res < eps, msg
     107    assert res < eps, msg
    112108   
    113109   
    114110    nominator = sum(abs(timeseries1-timeseries2))   
    115111    denominator = sum(abs(timeseries1))
    116     print nominator, denominator, N
    117112    if denominator > 0:
    118113        # Relative measure
     
    121116        # Absolute measure
    122117        res = nominator/N
    123                
     118
     119    print 'Sum abs diff', res   
    124120    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
    125     #assert res < eps, msg
     121    assert res < eps, msg
    126122
    127123    # Extrema
    128     #res = max(model)
    129     #report_difference('Maximum', res, expected_maximum[name], rtol, atol)
     124    max1 = max(timeseries1)
     125    max2 = max(timeseries2)
     126    print 'Max diff', abs(max1-max2)
     127
     128    min1 = min(timeseries1)
     129    min2 = min(timeseries2)
     130    print 'Min diff', abs(min1-min2)   
     131   
     132    # Locations of extrema
     133    i1 = num.argmax(timeseries1)
     134    i2 = num.argmax(timeseries2)
     135   
     136    res = abs(reference_time[i1]-reference_time[i2])
     137    print 'Max loc diff', res
     138
     139
     140    i1 = num.argmin(timeseries1)
     141    i2 = num.argmin(timeseries2)
     142   
     143    res = abs(reference_time[i1]-reference_time[i2])
     144    print 'Min loc diff', res   
    130145   
    131146   
     
    155170   
    156171
    157 def compare_all_timeseries(swwfile1, swwfile2):
     172def compare_all_timeseries(swwfile1, swwfile2, eps=1.0e-6):
    158173   
    159174    timevector1, timeseries1 = get_timeseries(swwfile1, gauges)
     
    177192                           timeseries1[name],
    178193                           timeseries2[name],
    179                            name=name)
     194                           name=name,
     195                           eps=eps)
    180196       
    181197       
     
    189205if __name__ == '__main__':
    190206
    191     if len(sys.argv) != 3:
     207    if len(sys.argv) != 4:
    192208        usage()
    193209        import sys; sys.exit(-1)
    194210
    195     #res = compare_all_timeseries(sys.argv[1], sys.argv[2])
     211    eps = float(sys.argv[3]) # Get tolerance to be used in comparisons
    196212    try:
    197         res = compare_all_timeseries(sys.argv[1], sys.argv[2])
     213        res = compare_all_timeseries(sys.argv[1], sys.argv[2], eps=eps)
    198214    except:
    199215        sys.exit(-1)
  • anuga_validation/automated_validation_tests/patong_beach_validation/validate.py

    r7645 r7649  
    5353                         'patong.sww.FINAL.tgz'
    5454                        )
    55 Optional_Data_Objects = ('patong.sww.BASIC.tgz',)
     55Optional_Data_Objects = ('patong.sww.TRIAL.tgz',)
     56
     57# Associated tolerances to be used in comparisons (would depend on discretisation errors)
     58epsilon = {'patong.sww.TRIAL.tgz': 1.0e-3,
     59           'patong.sww.BASIC.tgz': 1.0e-4,
     60           'patong.sww.FINAL.tgz': 1.0e-5}
     61
    5662
    5763# path to the local data directory
     
    314320    return res == 0
    315321
    316 def check_that_output_is_as_expected(expected_sww, valid_sww):
     322def check_that_output_is_as_expected(expected_sww, valid_sww, epsilon):
    317323    '''Check that validation output is as required.'''
    318324
     
    348354    new_output_sww = os.path.join(output_directory, expected_sww)
    349355    #cmd = 'python cmpsww.py %s %s > cmpsww.stdout' % (local_sww, new_output_sww)
    350     cmd = 'python compare_model_timeseries.py %s %s > compare_model_timeseries.stdout' % (local_sww, new_output_sww)
     356    cmd = 'python compare_model_timeseries.py %s %s %e > compare_model_timeseries.stdout' %\
     357          (local_sww, new_output_sww, epsilon)
    351358    print '-------------------------------------'
    352359    print cmd
     
    449456        (valid_sww, _) = odo.rsplit('.', 1)
    450457        (expected_sww, _) = valid_sww.rsplit('.', 1)
    451         check_that_output_is_as_expected(expected_sww, valid_sww)
     458        check_that_output_is_as_expected(expected_sww, valid_sww, epsilon[odo])
    452459
    453460    stop_time = time.time()
Note: See TracChangeset for help on using the changeset viewer.