# -*- coding: utf-8 -*-

#
#  report_diff.py 
# 
#  Date Created: July 11, 2018
#
#  Author:     Michael E. Tryby
#              US EPA - ORD/NRMRL
#

# system imports
import itertools as it

# third party imports
import numpy as np

# project imports
import nrtest_epanet.output_reader as ordr


def report_diff(path_test, path_ref, min_cdd):
    for (test, ref) in it.izip(ordr.output_generator(path_test), 
                               ordr.output_generator(path_ref)):
        
        if len(test[0]) != len(ref[0]):
            raise ValueError('Inconsistent lengths')
        
        # Skip over arrays that are equal
        if np.array_equal(test[0], ref[0]):
            continue
        else:
            lre = log_relative_error(test[0], ref[0])
            idx = np.unravel_index(np.argmin(lre), lre.shape)

            if lre[idx] < min_cdd: 
                print_diff(idx, lre, test, ref)

    return 


def log_relative_error(q, c):
    '''
    Computes log relative error, a measure of numerical accuracy. 
    
    Single precision machine epsilon is between 2^-24 and 2^-23.
    
    Reference:
    McCullough, B. D. "Assessing the Reliability of Statistical Software: Part I."
    The American Statistician, vol. 52, no. 4, 1998, pp. 358-366. 
    '''
    diff = np.subtract(q, c)
    tmp_c = np.copy(c)
    # If ref value is small compute absolute error
    tmp_c[np.fabs(tmp_c) <= 1.0e-6] = 1.0
    
    re = np.fabs(diff)/np.fabs(tmp_c)
    # If re is tiny set lre to number of digits
    re[re < 1.0e-7] = 1.0e-7
    # If re is very large set lre to zero
    re[re > 2.0] = 1.0

    return np.negative(np.log10(re))
    
    
def print_diff(idx, lre, test, ref):
    print("LRE: %f\nIdx: %s\nSut: %f\nRef: %f\n" 
          % ((lre[idx]),(idx[0], ref[1]),(test[0][idx[0]]),(ref[0][idx[0]])))