Merge branch 'dev' of https://github.com/OpenWaterAnalytics/EPANET into dev

include/epanet2.bas

@@ -154,6 +154,13 @@ Public Const EN_INITFLOW = 10    ' Re-initialize flow flag
 Public Const EN_CONST_HP = 0      ' constant horsepower
 Public Const EN_POWER_FUNC = 1    ' power function
 Public Const EN_CUSTOM = 2        ' user-defined custom curve
+Public Const EN_NOCURVE = 3       ' no curve
+
+Public Const EN_V_CURVE = 0       ' volume curve
+Public Const EN_P_CURVE = 1       ' pump curve
+Public Const EN_E_CURVE = 2       ' efficiency curve
+Public Const EN_H_CURVE = 3       ' head loss curve
+Public Const EN_G_CURVE = 4       ' General\default curve
 
 'These are the external functions that comprise the DLL
 

include/epanet2.h

@@ -224,8 +224,6 @@ typedef enum {
   EN_TIMEOFDAY   = 3
 } EN_ControlType;
 
-
-
 typedef enum {
   EN_AVERAGE     = 1,   /* Time statistic types.    */
   EN_MINIMUM     = 2,   /* See TstatType in TYPES.H */
@@ -233,8 +231,6 @@ typedef enum {
   EN_RANGE       = 4
 } EN_StatisticType;
 
-
-
 typedef enum {
   EN_MIX1        = 0,   /* Tank mixing models */
   EN_MIX2        = 1,
@@ -242,8 +238,6 @@ typedef enum {
   EN_LIFO        = 3
 } EN_MixingModel;
 
-
-
 typedef enum {
   EN_NOSAVE        = 0,
   EN_SAVE          = 1,
@@ -251,16 +245,21 @@ typedef enum {
   EN_SAVE_AND_INIT = 11
 } EN_SaveOption;
 
-
-
 typedef enum {
   EN_CONST_HP    = 0,   /* constant horsepower       */
   EN_POWER_FUNC  = 1,   /* power function            */
-  EN_CUSTOM      = 2    /* user-defined custom curve */
+  EN_CUSTOM      = 2,   /* user-defined custom curve */
+  EN_NOCURVE     = 3    /* no curve                  */
+} EN_PumpType;
+
+typedef enum {
+  EN_V_CURVE     = 0,   /*    volume curve                      */
+  EN_P_CURVE     = 1,   /*    pump curve                        */
+  EN_E_CURVE     = 2,   /*    efficiency curve                  */
+  EN_H_CURVE     = 3,   /*    head loss curve                   */
+  EN_G_CURVE     = 4    /*    General\default curve             */
 } EN_CurveType;
 
-
-
 // --- Declare the EPANET toolkit functions
 #if defined(__cplusplus)
 extern "C" {

src/epanet.c

@@ -2933,7 +2933,7 @@ int DLLEXPORT EN_addcurve(EN_Project *p, char *id) {
 
   strcpy(tmpCur[n].ID, id);
   tmpCur[n].Npts = 1;
-  tmpCur[n].Type = -1;
+  tmpCur[n].Type = G_CURVE;
   tmpCur[n].X = (double *)calloc(tmpCur[n].Npts, sizeof(double));
   tmpCur[n].Y = (double *)calloc(tmpCur[n].Npts, sizeof(double));
   if (tmpCur[n].X == NULL)
@@ -2974,9 +2974,7 @@ int DLLEXPORT EN_addcurve(EN_Project *p, char *id) {
 int DLLEXPORT EN_setcurve(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y, int n) {
   
   EN_Network *net = &p->network;  
-  Scurve *Curve = net->Curve;
+  Scurve *Curve = net->Curve;  
-  
-  
   int j;
 
   /* Check for valid arguments */
@@ -3007,12 +3005,9 @@ int DLLEXPORT EN_setcurve(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API
 int DLLEXPORT EN_setcurvevalue(EN_Project *p, int index, int pnt, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y) {
   
   EN_Network *net = &p->network;
-  
   Scurve *Curve = net->Curve;
-  
   const int Ncurves = net->Ncurves;
   
-  
   if (!p->Openflag)
     return (102);
   if (index <= 0 || index > Ncurves)
@@ -3713,7 +3708,7 @@ int allocdata(EN_Project *p)
     }
     for (n = 0; n <= par->MaxCurves; n++) {
       net->Curve[n].Npts = 0;
-      net->Curve[n].Type = -1;
+      net->Curve[n].Type = G_CURVE;
       net->Curve[n].X = NULL;
       net->Curve[n].Y = NULL;
     }

src/types.h

@@ -151,8 +151,9 @@ typedef enum {
   V_CURVE,       /*    volume curve                      */
   P_CURVE,       /*    pump curve                        */
   E_CURVE,       /*    efficiency curve                  */
-  H_CURVE
+  H_CURVE,       /*    head loss curve                   */
-} CurveType;     /*    head loss curve                   */
+  G_CURVE        /*    General\default curve             */
+} CurveType;
 
 typedef enum {
   CONST_HP,      /*    constant horsepower              */

tools/nrtest-epanet/main.py

@@ -4,7 +4,9 @@ import time
 import cStringIO
 import itertools as it
 
-import epanet_reader as er
+# project import
+import nrtest_epanet.output_reader as er
+
 
 def result_compare(path_test, path_ref, comp_args):
 
@@ -15,29 +17,36 @@ def result_compare(path_test, path_ref, comp_args):
     total = 0
     output = cStringIO.StringIO()
     eps =  np.finfo(float).eps
+    min_cdd = 100.0
     
     start = time.time()
     
-    test_reader = er.reader(path_test)
+    test_reader = er.output_generator(path_test)
-    ref_reader = er.reader(path_ref)
+    ref_reader = er.output_generator(path_ref)
     
     for test, ref in it.izip(test_reader, ref_reader):
         total += 1
         if total%100000 == 0:
             print(total)
         
-        if test.size != ref.size:
+        if len(test[0]) != len(ref[0]):
             raise ValueError('Inconsistent lengths')
         
         # Skip results if they are zero or equal
-        if np.array_equal(test, ref):
+        #if np.array_equal(test, ref):
-            equal += 1
+        #    equal += 1
-            continue
+        #    continue
         else:
             try:
-                np.testing.assert_allclose(test, ref, 1.0e-06, 2*eps)
+                diff = np.fabs(np.subtract(test[0], ref[0]))
-                close += 1
+                idx = np.unravel_index(np.argmax(diff), diff.shape)
                 
+                if diff[idx] != 0.0:
+                    tmp = - np.log10(diff[idx])
+                    
+                    if tmp < min_cdd:
+                        min_cdd = tmp;
+                        
             except AssertionError as ae:
                 notclose += 1
                 output.write(str(ae))
@@ -49,8 +58,9 @@ def result_compare(path_test, path_ref, comp_args):
     print(output.getvalue())
     output.close()
     
-    print('equal: %d  close: %d  notclose: %d  total: %d  in %f (sec)\n' % 
+    print('mincdd: %d in %f (sec)' % (np.floor(min_cdd), (stop - start)))
-          (equal, close, notclose, total, (stop - start)))
+    #print('equal: %d  close: %d  notclose: %d  total: %d  in %f (sec)\n' % 
+    #      (equal, close, notclose, total, (stop - start)))
 
     if notclose > 0:
         print('%d differences found\n' % notclose)
@@ -120,6 +130,8 @@ def nrtest_execute(app_path, test_path, output_path):
     exit(not success)
 
 
+import report_diff as rd
+
 if __name__ == "__main__":
 #    app_path = "apps\\swmm-5.1.11.json"
 #    test_path = "tests\\examples\\example1.json"
@@ -130,7 +142,10 @@ if __name__ == "__main__":
 #    ref_path  = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2012"
 #    print(nrtest_compare(test_path, ref_path, (0.001, 0.0)))
 
-    
+    benchmark_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\michaeltryby\\epanet-lr\\nrtestsuite\\benchmarks\\"
-    path_test = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2011a\\Example_3\\example3.out"
+    path_test = benchmark_path + "epanet-220dev\\example2\\example2.out"
-    path_ref  = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2012\\Example_3\\example3.out"    
+    path_ref  = benchmark_path + "epanet-2012\\example2\\example2.out"
-    result_compare(path_test, path_ref, (0.001, 0.0))
+
+    #result_compare(path_test, path_ref, (0.001, 0.0))
+    rd.report_diff(path_test, path_ref, 2)
+    

tools/nrtest-epanet/nrtest_epanet/__init__.py

@@ -28,7 +28,7 @@ __copyright__ = "None"
 __credits__ = "Colleen Barr, Maurizio Cingi, Mark Gray, David Hall, Bryant McDonnell"
 __license__ = "CC0 1.0 Universal"
 
-__version__ = "0.4.0"
+__version__ = "0.5.0"
 __date__ = "September 6, 2017"
 
 __maintainer__ = "Michael Tryby"
@@ -38,7 +38,7 @@ __status  = "Development"
 
 def epanet_allclose_compare(path_test, path_ref, rtol, atol):
     ''' 
-    Compares results in two EPANET binary files. Using the comparison criteria 
+    Compares results in two EPANET binary files using the comparison criteria 
     described in the numpy assert_allclose documentation. 
             
         (test_value - ref_value) <= atol + rtol * abs(ref_value) 
@@ -67,22 +67,67 @@ def epanet_allclose_compare(path_test, path_ref, rtol, atol):
     for (test, ref) in it.izip(ordr.output_generator(path_test), 
                                ordr.output_generator(path_ref)):
         
-        if len(test) != len(ref):
+        if len(test[0]) != len(ref[0]):
             raise ValueError('Inconsistent lengths')
         
         # Skip over arrays that are equal
-        if np.array_equal(test, ref):
+        if np.array_equal(test[0], ref[0]):
             continue
         else:
-            np.testing.assert_allclose(test, ref, rtol, atol)
+            np.testing.assert_allclose(test[0], ref[0], rtol, atol)
 
     return True
 
-# def epanet_better_compare(path_test, path_ref, rtol, atol):
+
-#     '''
+def epanet_mincdd_compare(path_test, path_ref, rtol, atol):
-#     If you don't like assert_allclose you can add another function here. 
+    '''
-#     '''
+    Compares the results of two EPANET binary files using a correct decimal
-#     pass
+    digits (cdd) comparison criteria:
+
+           min cdd(test, ref) >= atol
+
+    Returns true if min cdd in the file is greater than or equal to atol,
+    otherwise an AssertionError is thrown.
+
+    Arguments:
+        path_test - path to result file being testedgit
+        path_ref  - path to reference result file
+        rtol - ignored
+        atol - minimum allowable cdd value (i.e. 3)
+
+    Returns:
+        True
+
+    Raises:
+        ValueError()
+        AssertionError()
+    '''
+    min_cdd = 100.0
+
+    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:
+            diff = np.fabs(np.subtract(test[0], ref[0]))
+            idx = np.unravel_index(np.argmax(diff), diff.shape)
+
+            if diff[idx] != 0.0:
+                tmp = - np.log10(diff[idx])
+
+                if tmp < min_cdd:
+                    min_cdd = tmp;
+
+    if np.floor(min_cdd) >= atol:
+        return True
+    else:
+        raise AssertionError('min_cdd=%d less than atol=%g' % (min_cdd, atol))
+
 
 def epanet_report_compare(path_test, path_ref, rtol, atol):
     '''

tools/nrtest-epanet/nrtest_epanet/output_reader.py

@@ -23,7 +23,8 @@ def output_generator(path_ref):
     yield element attributes. It is useful for comparing contents of binary 
     files for numerical regression testing. 
     
-    The generator yields a Python list containing element attributes. 
+    The generator yields a Python tuple containing an array of element 
+    attributes and a tuple containing the element type, period, and attribute. 
     
     Arguments: 
         path_ref - path to result file
@@ -38,7 +39,8 @@ def output_generator(path_ref):
             for element_type in oapi.ElementType:
                 for attribute in br.elementAttributes[element_type]:
         
-                    yield br.element_attribute(element_type, period_index, attribute)
+                    yield (br.element_attribute(element_type, period_index, attribute), 
+                           (element_type, period_index, attribute))
 
 
 class OutputReader():    

tools/nrtest-epanet/report_diff.py

@@ -0,0 +1,96 @@
+# -*- 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 _binary_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):
+    
+    idx_val = (idx[0], ref[1])
+    test_val = (test[0][idx[0]])
+    ref_val = (ref[0][idx[0]])
+    diff_val = (test_val - ref_val)
+    lre_val = (lre[idx[0]])
+    
+    print("Idx: %s\nSut: %f  Ref: %f  Diff: %f  LRE: %f\n" 
+          % (idx_val, test_val, ref_val, diff_val, lre_val))
+
+
+def report(args):
+    _binary_diff(args.test, args.ref, args.mincdd)
+
+
+if __name__ == '__main__':
+    from argparse import ArgumentParser
+    
+    parser = ArgumentParser(description='EPANET benchmark difference reporting')
+    parser.set_defaults(func=report)
+    parser.add_argument('-t', '--test', default=None, 
+                        help='Path to test benchmark')
+    parser.add_argument('-r', '--ref', default=None, 
+                        help='Path to reference benchmark')
+    parser.add_argument('-mc', '--mincdd', type=int, default=3, 
+                        help='Minimum correct decimal digits')
+    
+    args = parser.parse_args()
+    args.func(args)
+    

tools/nrtest-epanet/setup.py

@@ -17,6 +17,7 @@ except ImportError:
 entry_points = {
     'nrtest.compare': [
         'epanet allclose = nrtest_epanet:epanet_allclose_compare',
+        'epanet mincdd = nrtest_epanet:epanet_mincdd_compare',
         'epanet report = nrtest_epanet:epanet_report_compare',
         # Add entry point for new comparison functions here
     ]
@@ -24,7 +25,7 @@ entry_points = {
 
 setup(
     name='nrtest-epanet',
-    version='0.4.0',
+    version='0.5.0',
     description="EPANET extension for nrtest",
     
     author="Michael E. Tryby",

win_build/WinSDK/Makefile.bat

@@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp
 If %ERRORLEVEL% == 1 (
 	CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release
 	rem : create EPANET2.DLL
-	cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL 
+	cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL 
 	rem : create EPANET2.EXE
-	cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
+	cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
 	md "%Build_PATH%"\64bit
 	move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit
 	move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit
@@ -35,9 +35,9 @@ rem : 32 bit with DEF
 CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release
 echo "32 bit with epanet2.def mapping"
 rem : create EPANET2.DLL
-cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
+cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
 rem : create EPANET2.EXE
-cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
+cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
 md "%Build_PATH%"\32bit
 move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit
 move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit