compare_MS takes data selection as extra argument

The data selection in python/compare_MS_numerically.py was hardcoded for
ES085A - but the script should work for any pair of measurementsets.

So the data selection was added as command line argument.
In test_ES085A.py the data selection for ES085 was added to the script
invocation.

ES085A/test_ES085A.py

@@ -199,12 +199,18 @@ class TTC(unittest.TestCase):
     # Compare the generated lisfile MS with the known-good one
     # We use the jplotter helper for this
     #
+    
+    # Select amp/phas vs freq from three different scans and sum them up.
+    # We expect the sum to be numerically identical between the known-good dataset 
+    # and the newly created one
+    __es085a_comparison__ = "indexr; scan 18 19 33; pt anpfreq; avt sum; solint none;"
     def testAA(self):
         rv = None
-        cmd = "{python} {compare} {gold_ms} es085a_cont.ms".format(
-                python  = options.binaries.python,
-                compare = os.path.join(TTC.__rootdir__, "../python", "compare_MS_numerically.py"),
-                gold_ms = os.path.join(TTC.__rootdir__, "es085a_cont.ms") )
+        cmd = "{python} {compare} '{selection}' {gold_ms} es085a_cont.ms".format(
+                python    = options.binaries.python,
+                compare   = os.path.join(TTC.__rootdir__, "../python", "compare_MS_numerically.py"),
+                selection = TTC.__es085a_comparison__,
+                gold_ms   = os.path.join(TTC.__rootdir__, "es085a_cont.ms") )
         # we need to prepare the environment
         jplEnv = copy.deepcopy(os.environ)
         # Make sure the pythonpath to jplotter is added
@@ -257,10 +263,14 @@ class TTC(unittest.TestCase):
 
     def testBA(self):
         rv = None
-        cmd = "{python} {compare} {gold_ms} es085a.ms".format(
-                python  = options.binaries.python,
-                compare = os.path.join(TTC.__rootdir__, "../python", "compare_MS_numerically.py"),
-                #gold_ms = os.path.join(TTC.__rootdir__, "es085a-hack.ms") )
+        cmd = "{python} {compare} '{selection}' {gold_ms} es085a.ms".format(
+                python    = options.binaries.python,
+                compare   = os.path.join(TTC.__rootdir__, "../python", "compare_MS_numerically.py"),
+                selection = TTC.__es085a_comparison__,
+                # this "-hack" MeasurementSet has one visibility modified (row 0, datapoint 0) by ~1x10-6
+                # and later had one row duplicated - i.e. one visibility extra on one baseline for
+                # one source for one spectral window
+                #gold_ms   = os.path.join(TTC.__rootdir__, "es085a-hack.ms") )
                 gold_ms = os.path.join(TTC.__rootdir__, "es085a.ms") )
         # we need to prepare the environment
         jplEnv = copy.deepcopy(os.environ)
@@ -326,6 +336,13 @@ class TTC(unittest.TestCase):
                 --ms es085a_cont.ms --ms es085a.ms --idi ES085A_CONT.IDI --idi ES085A.IDI".format(
                     python   = options.binaries.python,
                     verify   = os.path.join(TTC.__rootdir__, "../jive-toolchain-verify/compare-ms-idi.py"),
+                    # the "-hack" MeasurementSet has one visibility modified (row 0, datapoint 0) by ~1x10-6
+                    # and later had one row duplicated - i.e. one 2second visibility extra on one baseline for
+                    # one source for one spectral window.
+                    # This test should fail with a report showing 2s of data missing from the produced
+                    # data sets compared to the "gold" MS on Baseline McIb source 3c345 ...
+                    # (and reporting different number of integrations 156 in stead of 155)
+                    #gold_ms  = os.path.join(TTC.__rootdir__, "es085a-hack.ms"),
                     gold_ms  = os.path.join(TTC.__rootdir__, "es085a_cont.ms"),
                     gold_idi = os.path.join(TTC.__rootdir__, "ES085A_CONT.IDI"),
               )

python/compare_MS_numerically.py

@@ -9,15 +9,15 @@
 from __future__ import print_function
 import          __main__, os, sys, command, jplotter, shutil, tempfile, traceback, pyrap
 
-if len(sys.argv)!=3:
-    raise RuntimeError("Usage: {0} path/to/reference.ms path/to/test.ms".format( sys.argv[0] ))
+if len(sys.argv)!=4:
+    raise RuntimeError("Usage: {0} '<jplotter data selection>' path/to/reference.ms path/to/test.ms".format( sys.argv[0] ))
 
 # Since jplotter is very forgiving when the user attempts to open a non-existing MS
 # it doesn't generate an error in that case. So we verify here that both measurement sets at least _seem_ to exist
-with pyrap.tables.table(sys.argv[1]) as test:
-    pass
 with pyrap.tables.table(sys.argv[2]) as test:
     pass
+with pyrap.tables.table(sys.argv[3]) as test:
+    pass
 
 # add the path to this module such that Python can load the compare_ms_data postprocessing module
 sys.path.insert(0, os.path.dirname(__main__.__file__))
@@ -27,9 +27,11 @@ mk_plots = """
 refile {tempfile}
 postprocess ./compare_ms_data.compare_ms_data
 ms {reference}
-indexr; scan 18 19 33; pt anpfreq; avt sum; solint none; y0 local; y1 local; multi t; new src t; sort src bl; pl;
+{select_data};
+pl;
 ms {test}
-indexr; scan 18 19 33; pt anpfreq; avt sum; solint none; y0 local; y1 local; multi t; new src t; sort src bl; pl;
+{select_data};
+pl;
 """.format
 
 # We need to capture errors without exiting Python
@@ -63,7 +65,7 @@ try:
     sys.exit         = myexitfn
     sys.excepthook   = myexcepthook
     cmd              = mk_plots(tempfile=os.path.join(tmpdir,'plotfile.ps'),
-                                reference=sys.argv[1], test=sys.argv[2])
+                                select_data=sys.argv[1], reference=sys.argv[2], test=sys.argv[3])
     print("====> jplotter script")
     print(cmd)
     print("<==== jplotter script")

python/compare_ms_data.py

@@ -1,5 +1,5 @@
 from   __future__ import print_function
-import sys, numpy, plots, math, operator, copy
+import sys, numpy, plots, math, operator, copy, collections
 
 COPY = copy.deepcopy
 # define a post-processing operation that always remembers the last data set and 
@@ -40,25 +40,20 @@ def compare_ms_data(plotar, ms2mappings):
         new_keys = set(new_dataset.keys())
         # if the sets are not equal, the data will also not compare equal!
         if old_keys == new_keys:
-            # inspect x and y separately, add up all the diffs
-            dx, dy = 0, 0
+            # inspect x's and the y's separately, add up all the diffs per type
+            d = collections.defaultdict(int)
             for k in old_keys:
                 ods = last_dataset[ k ]
                 nds = new_dataset[ k ]
-                dx  = numpy.add(numpy.abs( ods.xval - nds.xval ), dx)
-                dy  = numpy.add(numpy.abs( ods.yval - nds.yval ), dy)
-            #print("dx=", dx)
-            #print("dy=", dy)
-            if numpy.any( dx>abs(tolerance) ):
-                print(">>> compare_data: total diffs in x exceed tolerance")
-                print("    tolerance=", tolerance)
-                print("    accum. dx=", dx)
-                error += "Datasets mismatch in X according to tolerance. "
-            if numpy.any( dy>abs(tolerance) ):
-                print(">>> compare_data: total diffs in y exceed tolerance")
-                print("    tolerance=", tolerance)
-                print("    accum. dy=", dy)
-                error = "Datasets mismatch in Y according to tolerance. "
+                d['x']    = numpy.add(numpy.abs(ods.xval - nds.xval), d['x'])
+                d[k.TYPE] = numpy.add(numpy.abs(ods.yval - nds.yval), d[k.TYPE])
+            # ok sum-of-all-diffs computes, now check if any of'em are > tolerance
+            for k,v in d.items():
+                if numpy.any( v>abs(tolerance) ):
+                    print(">>> compare_data: total diffs in", k, "exceed tolerance")
+                    print("    tolerance=", tolerance)
+                    print("    accum. d=", v)
+                    error += "Datasets mismatch in {0} according to tolerance. ".format(k)
         else:
             print(">>> compare_data: The datasets to compare have different data content?!")
             common = old_keys & new_keys