DM-14806: bring demo to standards

.gitignore

@@ -1,5 +1,19 @@
+.cache
+pytest_session.txt
 output
 output_small
 detected-sources.txt
 detected-sources_small.txt
 astrometry_sdss.png
+bin/check_astrometry.py
+bin/compare.py
+bin/export-results.py
+
+.coverage
+.sconf_temp/
+.sconsign.dblite
+config.log
+python/lsst/lsst/dm/stack/demo/version.py
+tests/.tests/
+tests/__pycache__/
+ups/__pycache__/

.travis.yaml

@@ -0,0 +1,8 @@
+sudo: false
+language: python
+matrix:
+  include:
+    - python: '3.6'
+      install:
+      - pip install flake8
+      script: flake8

README

@@ -14,7 +14,7 @@ Usage
 To run the example, do::
 
 	$ source $LSST_HOME/loadLSST.sh
-	$ setup obs_sdss
+	$ setup -r .
 	$ ./bin/demo.sh  (or: "./bin/demo.sh --small"  if local memory is constrained)
 
 where $LSST_HOME is the directory where your LSST DM stack resides.
@@ -35,7 +35,7 @@ output; choose a result set appropriate for your system. The ./bin/compare
 script can then be used to check if your results are within the expected
 numerical tolerances::
 
-    $ python bin/compare detected-sources.txt
+    $ python bin/compare.py detected-sources.txt
     Ok.
 
 If you use the "--small" run option then the corresponding output file is
@@ -44,7 +44,7 @@ called "output_small".
 
 Check the astrometric relative RMS with::
 
-    $ python bin/check_astrometry output
+    $ python bin/check_astrometry.py output
 
 Included Data
 -------------

SConstruct

@@ -0,0 +1,3 @@
+# -*- python -*-
+from lsst.sconsUtils import scripts
+scripts.BasicSConstruct("lsst_dm_stack_demo", disableCc=True)

bin.src/SConscript

@@ -0,0 +1,3 @@
+# -*- python -*-
+from lsst.sconsUtils import scripts
+scripts.BasicSConscript.shebang()

bin/check_astrometry → bin.src/check_astrometry.py

@@ -26,7 +26,6 @@
 import os.path
 import sys
 
-import matplotlib.pylab as plt
 import numpy as np
 
 import lsst.daf.persistence as dafPersist
@@ -138,22 +137,26 @@ def loadAndMatchData(repo, visits, fields, ref, ref_field, camcol, filter):
             did = {'run': ref, 'filter': filter, 'field': ref_field, 'camcol': camcolRef}
             md = butler.get("calexp_md", did, immediate=True)
             calib = afwImage.Calib(md)
+            calib.setThrowOnNegativeFlux(False)
             # compute magnitude
             refMag = calib.getMagnitude(mRef.get('base_PsfFlux_flux'))
 
             mag.append(refMag)
             dist.append(ang)
 
     return pipeBase.Struct(
-        mag = mag,
-        dist = dist,
-        match = matchNum
+        mag=mag,
+        dist=dist,
+        match=matchNum
     )
 
 
 def plotAstrometry(mag, dist, match, good_mag_limit=19.5):
     """Plot angular distance between matched sources from different exposures."""
 
+    # Defer importing of matplotlib until we need it.
+    import matplotlib.pylab as plt
+
     plt.rcParams['axes.linewidth'] = 2
     plt.rcParams['mathtext.default'] = 'regular'
 
@@ -208,7 +211,9 @@ def checkAstrometry(mag, dist, match,
     @param medianRef  Median reference astrometric scatter in arcseconds.
     @param matchRef   Should match at least matchRef stars.
 
-    Return the astrometric scatter (RMS, arcsec) for all good stars.
+    Return a boolean indicating whether the astrometric scatter was less than
+    the supplied reference, and the astrometric scatter (RMS, arcsec) for all
+    good stars.
 
     Notes:
        The scatter and match defaults are appropriate to SDSS are stored here.
@@ -225,23 +230,36 @@ def checkAstrometry(mag, dist, match,
     if astromScatter > medianRef:
         print("Median astrometric scatter %.1f mas is larger than reference : %.1f mas " %
               (astromScatter, medianRef))
+        passed = False
+    else:
+        passed = True
     if match < matchRef:
-        print("Number of matched sources %d is too small (shoud be > %d)" % (match, matchRef))
+        print("Number of matched sources %d is too small (should be > %d)" % (match, matchRef))
+        passed = False
 
-    return astromScatter
+    return passed, astromScatter
 
 
-def main(repo, runs, fields, ref, ref_field, camcol, filter):
+def main(repo, runs, fields, ref, ref_field, camcol, filter, plot=False):
     """Main executable.
 
+    Returns True if the test passed, False otherwise.
     """
 
     struct = loadAndMatchData(repo, runs, fields, ref, ref_field, camcol, filter)
     mag = struct.mag
     dist = struct.dist
     match = struct.match
-    checkAstrometry(mag, dist, match)
-    plotAstrometry(mag, dist, match)
+
+    # Limit depends on filter
+    medianRef = 100
+    if filter == 'i':
+        medianRef = 105
+
+    passed, astromScatter = checkAstrometry(mag, dist, match, medianRef=medianRef)
+    if plot:
+        plotAstrometry(mag, dist, match)
+    return passed
 
 
 def defaultData(repo):
@@ -255,7 +273,7 @@ def defaultData(repo):
 
     # List of camcol to be considered (source calalogs will be concateneted)
     camcol = [4]
-    filter = 'r'
+    filter = 'i'
 
     return runs, fields, ref, ref_field, camcol, filter
 
@@ -273,4 +291,8 @@ def defaultData(repo):
         sys.exit(1)
 
     runs, fields, ref, ref_field, camcol, filter = defaultData(repo)
-    main(repo, runs, fields, ref, ref_field, camcol, filter)
+    passed = main(repo, runs, fields, ref, ref_field, camcol, filter)
+    if passed:
+        print("Ok.")
+    else:
+        sys.exit(1)

bin.src/compare.py

@@ -0,0 +1,185 @@
+#!/usr/bin/env python
+from __future__ import division
+from __future__ import print_function
+
+import argparse
+import os
+import sys
+import numpy as np
+
+# Data columns in the output
+# This needs to be kept in sync with the columns in export-results, until
+# this script and the expected results files are modified to run directly
+# on the binary table outputs.
+# Note that we load the flags as 'str', since that's what export-results
+# does, but that doesn't matter since we're just comparing for equality
+# for all non-float fields.
+DTYPE = np.dtype(
+    [("id", np.int64),
+     ("coord_ra", float),
+     ("coord_dec", float),
+     ("flags_negative", str, 5),
+     ("base_SdssCentroid_flag", str, 5),
+     ("base_PixelFlags_flag_edge", str, 5),
+     ("base_PixelFlags_flag_interpolated", str, 5),
+     ("base_PixelFlags_flag_interpolatedCenter", str, 5),
+     ("base_PixelFlags_flag_saturated", str, 5),
+     ("base_PixelFlags_flag_saturatedCenter", str, 5),
+     ("base_SdssCentroid_x", float),
+     ("base_SdssCentroid_y", float),
+     ("base_SdssCentroid_xSigma", float),
+     ("base_SdssCentroid_ySigma", float),
+     ("base_SdssShape_xx", float),
+     ("base_SdssShape_xy", float),
+     ("base_SdssShape_yy", float),
+     ("base_SdssShape_xxSigma", float),
+     ("base_SdssShape_xySigma", float),
+     ("base_SdssShape_yySigma", float),
+     ("base_SdssShape_flag", str, 5),
+     ("base_GaussianFlux_flux", float),
+     ("base_GaussianFlux_fluxSigma", float),
+     ("base_PsfFlux_flux", float),
+     ("base_PsfFlux_fluxSigma", float),
+     ("base_CircularApertureFlux_2_flux", float),
+     ("base_CircularApertureFlux_2_fluxSigma", float),
+     ("base_ClassificationExtendedness_value", float),
+     ])
+
+
+def get_array(filename):
+    with open(filename, 'r') as f:
+        array = np.loadtxt(f, dtype=DTYPE)
+    return array
+
+
+def difference(arr1, arr2):
+    """
+    Compute the relative and absolute differences of numpy arrays arr1 & arr2.
+
+    The relative difference R between numbers n1 and n2 is defined as per
+    numdiff (http://www.nongnu.org/numdiff/numdiff.html):
+    * R = 0 if n1 and n2 are equal,
+    * R = Inf if n2 differs from n1 and at least one of them is zero,
+    * R = A/ min(|n1|, |n2|) if n1 and n2 are both non zero and n2 differs from n1.
+    """
+    absDiff = np.abs(arr1 - arr2)
+
+    # If there is a difference between 0 and something else, the result is
+    # infinite.
+    absDiff = np.where((absDiff != 0) & ((arr1 == 0) | (arr2 == 0)), np.inf, absDiff)
+
+    # If both inputs are nan, the result is 0.
+    absDiff = np.where(np.isnan(arr1) & np.isnan(arr2), 0, absDiff)
+
+    # If one input is nan, the result is infinite.
+    absDiff = np.where(np.logical_xor(np.isnan(arr1), np.isnan(arr2)), np.inf, absDiff)
+
+    # Divide by the minimum of the inputs, unless 0 or nan.
+    # If the minimum is 0 or nan, then either both inputs are 0/nan (so there's no
+    # difference) or one is 0/nan (in which case the absolute difference is
+    # already inf).
+    divisor = np.where(np.minimum(arr1, arr2) == 0, 1, np.minimum(arr1, arr2))
+    divisor = np.where(np.isnan(divisor), 1, divisor)
+
+    return absDiff, absDiff/np.abs(divisor)
+
+
+def compareWithNumPy(filename, reference, tolerance):
+    """
+    Compare a generated data file to a reference using NumPy.
+
+    The comparison is successful if:
+    * The files contain the same data columns (both by name and by data type);
+    * Each numeric value is in the input is either:
+        a) Within ``tolerance`` of the corresponding value in the reference, or
+        b) The relative difference with the reference (defined as above) is
+           within ``tolerance``;
+    * Flags recorded in the input and the reference are identical.
+
+    @param filename  Path to input data file.
+    @param reference Path to reference file.
+    @param tolerance Tolerance.
+    """
+    with open(filename, 'r') as data, open(reference, 'r') as ref:
+        data_columns = data.readline().strip('#').split()
+        ref_columns = ref.readline().strip('#').split()
+    table1, table2 = get_array(filename), get_array(reference)
+    if (table1.dtype != table2.dtype) or (data_columns != ref_columns):
+        print("Files do not contain the same columns.")
+        return False
+    valid = True
+    for name in table1.dtype.names:
+        dtype, count = table1.dtype.fields[name]
+        if dtype == np.dtype(float):
+            absDiff, relDiff = difference(table1[name], table2[name])
+            for pos in np.where((relDiff > tolerance) & (absDiff > tolerance))[0]:
+                valid = False
+                print("Failed (absolute difference %g, relative difference %g over tolerance %g) "
+                      "in column %s." % (absDiff[pos], relDiff[pos], tolerance, name))
+        else:
+            if not np.all(table1[name] == table2[name]):
+                nTotal = len(table1[name])
+                nDiff = len(np.where(table1[name] != table2[name])[0])
+                print("Failed (%s of %s flags do not match) in column %s." % (str(nDiff), str(nTotal), name))
+                valid = False
+    return valid
+
+
+def determineFlavor():
+    """
+    Return a string representing the 'flavor' of the local system.
+
+    Based on the equivalent logic in EUPS, but without introducing an EUPS
+    dependency.
+    """
+    uname, machine = os.uname()[0:5:4]
+    if uname == "Linux":
+        if machine[-2:] == "64":
+            return "Linux64"
+        else:
+            return "Linux"
+    elif uname == "Darwin":
+        if machine in ("x86_64", "i686"):
+            return "DarwinX86"
+        else:
+            return "Darwin"
+    else:
+        raise RuntimeError("Unknown flavor: (%s, %s)" % (uname, machine))
+
+
+def extantFile(filename):
+    """
+    Raise ArgumentTypeError if ``filename`` does not exist.
+    """
+    if not os.path.isfile(filename):
+        raise argparse.ArgumentTypeError(filename + " is not a file.")
+    return filename
+
+
+def referenceFilename(checkFilename):
+    """
+    Attempt to guess the filename to compare our input against.
+    """
+    guess = os.path.join(os.path.split(os.path.dirname(__file__))[0], "expected",
+                         determineFlavor(), os.path.basename(checkFilename))
+    if os.path.isfile(guess):
+        return guess
+    else:
+        raise ValueError("Cannot find reference data (looked for %s)." % (guess,))
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument('filename', type=extantFile, help="Input data file.")
+    parser.add_argument('--tolerance', default=0, type=float, help="Tolerance for errors. "
+                        "The test will fail if both the relative and absolute errors exceed the tolerance.")
+    parser.add_argument('--reference', type=extantFile, help="Reference data for comparison.")
+    args = parser.parse_args()
+
+    if not args.reference:
+        args.reference = referenceFilename(args.filename)
+
+    if compareWithNumPy(args.filename, args.reference, args.tolerance):
+        print("Ok.")
+    else:
+        sys.exit(1)