DM-4957: Generate JSON output from validate_drp for inclusion in a test harness

README.md

@@ -4,7 +4,7 @@ and validate the astrometric and photometric repeatability of the results.
 
 Pre-requisites: install and declare the following
 
-1. `pipe_tasks` from the LSST DM stack (note that pipe_tasks is included with lsst_apps, which is the usual thing to install)
+1. `pipe_base` from the LSST DM stack (note that `pipe_base` is included with `lsst_apps`, which is the usual thing to install)
 2. `obs_decam` from https://github.com/lsst/obs_decam
 3. `obs_cfht` from https://github.com/lsst/obs_cfht
 4. `validation_data_cfht` from https://github.com/lsst/validation_data_cfht

python/lsst/validate/drp/base.py

@@ -1,5 +1,3 @@
-#!/usr/bin/env python
-
 # LSST Data Management System
 # Copyright 2008-2016 AURA/LSST.
 #
@@ -26,3 +24,13 @@
 class ValidateError(Exception):
     """Base classes for exceptions in validate_drp."""
     pass
+
+
+class ValidateErrorNoStars(ValidateError):
+    """To be returned by tests that find no stars satisfying a set of criteria.
+
+    Some example cases that might return such an error:
+    1. There are no stars between 19-21 arcmin apart.
+    2. There are no stars in the given magnitude range.
+    """
+    pass

python/lsst/validate/drp/calcSrd.py

@@ -1,5 +1,3 @@
-#!/usr/bin/env python
-
 # LSST Data Management System
 # Copyright 2008-2016 AURA/LSST.
 #
@@ -20,7 +18,7 @@
 # the GNU General Public License along with this program.  If not,
 # see <https://www.lsstcorp.org/LegalNotices/>.
 
-from __future__ import print_function, division
+from __future__ import print_function, division, absolute_import
 
 import math
 
@@ -29,30 +27,45 @@
 
 import lsst.pipe.base as pipeBase
 
+from .base import ValidateErrorNoStars
 from .util import averageRaFromCat, averageDecFromCat
-from .srdSpec import srdSpec
+from .srdSpec import srdSpec, getAstrometricSpec
 
 
-def calcPA1(groupView, magKey):
+def calcPA1(matches, magKey, numRandomShuffles=50):
     """Calculate the photometric repeatability of measurements across a set of observations.
 
     Parameters
     ----------
-    groupView : lsst.afw.table.GroupView
+    matches : lsst.afw.table.GroupView
          GroupView object of matched observations from MultiMatch.
     magKey : lookup key to a `schema`
          The lookup key of the field storing the magnitude of interest.
          E.g., `magKey = allMatches.schema.find("base_PsfFlux_mag").key`
-         where `allMatches` is a the result of lsst.afw.table.MultiMatch.finish()
+         where `allMatches` is the result of lsst.afw.table.MultiMatch.finish()
+    numRandomShuffles : int
+        Number of times to draw random pairs from the different observations.
 
     Returns
     -------
     pipeBase.Struct
-       The RMS, inter-quartile range,
-       differences between pairs of observations, mean mag of each object.
+        name -- str: 'PA1'.  Name of Key Performance Metric stored in this struct.
+        rms -- average RMS
+        iqr -- and average inter-quartile range (IQR)
+        rmsStd -- standard deviation of the RMS
+        iqrStd -- standard deviation of the IQR
+           These 4 quantities are derived from the `numRandomShuffles` trials.
+        magDiffs -- differences between pairs of observations.  Selected from one random trial.
+        magMean -- mean magnitude of each object.
 
     Notes
     -----
+    We calculate differences for 50 different random realizations
+    of the measurement pairs, to provide some estimate of the uncertainty
+    on our RMS estimates due to the random shuffling.
+    This estimate could be stated and calculated from a more formally
+    derived motivation but in practice 50 should be sufficient.
+
     The LSST Science Requirements Document (LPM-17), or SRD,
     characterizes the photometric repeatability by putting a requirement
     on the median RMS of measurements of non-variable bright stars.
@@ -76,6 +89,7 @@ def calcPA1(groupView, magKey):
 
     See Also
     --------
+    doCalcPA1 : Worker routine that calculates PA1 for one random sample.
     calcPA2 : Calculate photometric repeatability outlier fraction.
 
     Examples
@@ -99,16 +113,49 @@ def calcPA1(groupView, magKey):
     >>> allMatches = GroupView.build(matchCat)
     >>> allMatches
     >>> psfMagKey = allMatches.schema.find("base_PsfFlux_mag").key
-    >>> pa1_sample = calcPA1(allMatches, psfMagKey)
-    >>> print("The RMS was %.3f, the IQR was %.3f" % (pa1_sample.rms, pa1_sample.iqr))
+    >>> pa1 = calcPA1(allMatches, psfMagKey)
+    >>> print("LSST SRD Key Performance Metric %s" % pa1.name)
+    >>> print("The RMS was %.3f+-%.3f, the IQR was %.3f+-%.3f" % (pa1.rms, pa1.iqr, pa1.rmsStd, pa1.iqrStd))
     """
 
-    diffs = groupView.aggregate(getRandomDiffRmsInMas, field=magKey)
-    means = groupView.aggregate(np.mean, field=magKey)
-    rmsPA1, iqrPA1 = computeWidths(diffs)
+    pa1Samples = [doCalcPA1(matches, magKey) for n in range(numRandomShuffles)]
+    rmsPA1 = np.array([pa1.rms for pa1 in pa1Samples])
+    iqrPA1 = np.array([pa1.iqr for pa1 in pa1Samples])
+
+    return pipeBase.Struct(name='PA1',
+                           rms=np.mean(rmsPA1), iqr=np.mean(iqrPA1),
+                           rmsStd=np.std(rmsPA1), iqrStd=np.std(iqrPA1),
+                           rmsUnits='mmag', iqrUnits='mmag',
+                           magDiffs=pa1Samples[0].magDiffs, magMean=pa1Samples[0].magMean,
+                           magDiffsUnits='mmag', magMeanUnits='mag')
 
+
+def doCalcPA1(groupView, magKey):
+    """Calculate PA1 for one random realization.
+
+    Parameters
+    ----------
+    groupView : lsst.afw.table.GroupView
+         GroupView object of matched observations from MultiMatch.
+    magKey : lookup key to a `schema`
+         The lookup key of the field storing the magnitude of interest.
+         E.g., `magKey = allMatches.schema.find("base_PsfFlux_mag").key`
+         where `allMatches` is the result of lsst.afw.table.MultiMatch.finish()
+
+    Returns
+    -------
+    pipeBase.Struct
+       The RMS, inter-quartile range,
+       differences between pairs of observations, mean mag of each object.
+    """
+
+    magDiffs = groupView.aggregate(getRandomDiffRmsInMas, field=magKey)
+    magMean = groupView.aggregate(np.mean, field=magKey)
+    rmsPA1, iqrPA1 = computeWidths(magDiffs)
     return pipeBase.Struct(rms=rmsPA1, iqr=iqrPA1,
-                           diffs=diffs, means=means)
+                           rmsUnits='mmag', iqrUnits='mmag',
+                           magDiffs=magDiffs, magMean=magMean,
+                           magDiffsUnits='mmag', magMeanUnits='mag')
 
 
 def calcPA2(groupView, magKey):
@@ -129,7 +176,8 @@ def calcPA2(groupView, magKey):
     Returns
     -------
     pipeBase.Struct
-       Contains PA2, the millimags of varaition at the
+        name -- str: 'PA2'.  Name of Key Performance Metric stored in this struct.
+       Contains the results of calculating PA2, the millimags of variation at the
        `design`, `minimum`, and `stretch` fraction of outliers
        The specified fractions are also avialable as `PF1`.
 
@@ -159,11 +207,12 @@ def calcPA2(groupView, magKey):
     >>> allMatches
     >>> psfMagKey = allMatches.schema.find("base_PsfFlux_mag").key
     >>> pa2 = calcPA2(allMatches, psfMagKey)
-    >>> print("minimum: PF1=%2d%% of diffs more than PA2 = %4.2f mmag (target is PA2 < 15 mmag)" %
+    >>> print("LSST SRD Key Performance Metric %s" % pa2.name)
+    >>> print("minimum: PF1=%2d%% of magDiffs more than PA2 = %4.2f mmag (target is PA2 < 15 mmag)" %
     ...       (pa2.PF1['minimum'], pa2.minimum))
-    >>> print("design:  PF1=%2d%% of diffs more than PA2 = %4.2f mmag (target is PA2 < 15 mmag)" %
+    >>> print("design:  PF1=%2d%% of magDiffs more than PA2 = %4.2f mmag (target is PA2 < 15 mmag)" %
     ...       (pa2.PF1['design'], pa2.design))
-    >>> print("stretch: PF1=%2d%% of diffs more than PA2 = %4.2f mmag (target is PA2 < 10 mmag)" %
+    >>> print("stretch: PF1=%2d%% of magDiffs more than PA2 = %4.2f mmag (target is PA2 < 10 mmag)" %
     ...       (pa2.PF1['stretch'], pa2.stretch))
 
 
@@ -176,11 +225,13 @@ def calcPA2(groupView, magKey):
       following LPM-17 as of 2011-07-06, available at http://ls.st/LPM-17.
     """
 
-    diffs = groupView.aggregate(getRandomDiffRmsInMas, field=magKey)
+    magDiffs = groupView.aggregate(getRandomDiffRmsInMas, field=magKey)
     PF1 = {'minimum': 20, 'design': 10, 'stretch': 5}
     PF1_percentiles = 100 - np.asarray([PF1['minimum'], PF1['design'], PF1['stretch']])
-    minPA2, designPA2, stretchPA2 = np.percentile(np.abs(diffs), PF1_percentiles)
-    return pipeBase.Struct(design=designPA2, minimum=minPA2, stretch=stretchPA2, PF1=PF1)
+    minPA2, designPA2, stretchPA2 = np.percentile(np.abs(magDiffs), PF1_percentiles)
+    return pipeBase.Struct(name='PA2', pa2Units='mmag', pf1Units='%',
+                           design=designPA2, minimum=minPA2, stretch=stretchPA2,
+                           PF1=PF1)
 
 
 def getRandomDiffRmsInMas(array):
@@ -200,7 +251,7 @@ def getRandomDiffRmsInMas(array):
     -----
     The LSST SRD recommends computing repeatability from a histogram of
     magnitude differences for the same star measured on two visits
-    (using a median over the diffs to reject outliers).
+    (using a median over the magDiffs to reject outliers).
     Because we have N>=2 measurements for each star, we select a random
     pair of visits for each star.  We divide each difference by sqrt(2)
     to obtain RMS about the (unknown) mean magnitude,
@@ -319,8 +370,8 @@ def matchVisitComputeDistance(visit_obj1, ra_obj1, dec_obj1,
     For each visit shared between visit_obj1 and visit_obj2,
     calculate the spherical distance between the obj1 and obj2.
 
-    Inputs
-    ------
+    Parameters
+    ----------
     visit_obj1 : scalar, list, or numpy.array of int or str
         List of visits for object 1.
     ra_obj1 : scalar, list, or numpy.array of float
@@ -361,21 +412,22 @@ def calcAM1(*args, **kwargs):
     """Calculate the SRD definition of astrometric performance for AM1
 
     See `calcAMx` for more details."""
-    return calcAMx(*args, D=srdSpec.D1, width=2, **kwargs)
+    return calcAMx(*args, x=1, D=srdSpec.D1, width=2, **kwargs)
 
 def calcAM2(*args, **kwargs):
     """Calculate the SRD definition of astrometric performance for AM2
 
     See `calcAMx` for more details."""
-    return calcAMx(*args, D=srdSpec.D2, width=2, **kwargs)
+    return calcAMx(*args, x=2, D=srdSpec.D2, width=2, **kwargs)
 
 def calcAM3(*args, **kwargs):
     """Calculate the SRD definition of astrometric performance for AM3
 
     See `calcAMx` for more details."""
-    return calcAMx(*args, D=srdSpec.D3, width=2, **kwargs)
+    return calcAMx(*args, x=3, D=srdSpec.D3, width=2, **kwargs)
 
-def calcAMx(groupView, D=5, width=2, magRange=None):
+def calcAMx(groupView, D=5, width=2, magRange=None,
+            x=None, level="design"):
     """Calculate the SRD definition of astrometric performance
 
     Parameters
@@ -389,12 +441,20 @@ def calcAMx(groupView, D=5, width=2, magRange=None):
     magRange : 2-element list or tuple
         brighter, fainter limits of the magnitude range to include.
         E.g., `magRange=[17.5, 21.0]`
+    x : int
+        Which of AM1, AM2, AM3.  One of [1,2,3].
+    level : str
+        One of "minimum", "design", "stretch" indicating the level of the specification desired.
 
     Returns
     -------
-    list, list, list
+    pipeBase.Struct
         rmsDistMas, annulus, magRange
 
+    Raises
+    ------
+    ValueError if `x` isn't in `getAstrometricSpec` values of [1,2,3]
+
     Notes
     -----
     This table below is provided in this package in the `srdSpec.py` file.
@@ -436,6 +496,66 @@ def calcAMx(groupView, D=5, width=2, magRange=None):
     and to astrometric measurements performed in the r and i bands.
     """
 
+    AMx_spec, AFx_spec, ADx_spec = getAstrometricSpec(x=x, level=level)
+
+    annulus = D + (width/2)*np.array([-1, +1])
+
+    rmsDistances, annulus, magRange = \
+        calcRmsDistances(groupView, annulus, magRange=magRange)
+
+    if not list(rmsDistances):
+        raise ValidateErrorNoStars('No stars found that are %.1f--%.1f arcmin apart.' %
+                                   (annulus[0], annulus[1]))
+
+    rmsDistMas = radiansToMilliarcsec(rmsDistances)
+    AMx = np.median(rmsDistMas)
+    fractionOver = np.mean(np.asarray(rmsDistMas) > AMx_spec+ADx_spec)
+
+    return pipeBase.Struct(
+        name='AM%d' % x,
+        AMx=AMx,
+        amxUnits='mas',
+        rmsDistMas=rmsDistMas,
+        rmsUnits='mas',
+        fractionOver=fractionOver,
+        D=D,
+        DUnits='arcmin',
+        annulus=annulus,
+        annulusUnits='arcmin',
+        magRange=magRange,
+        magRangeUnits='mag',
+        x=x,
+        level=level,
+        AMx_spec=AMx_spec,
+        AFx_spec=AFx_spec,
+        ADx_spec=ADx_spec,
+        afxUnits='%',
+        adxUnits='mas',
+        )
+
+
+def calcRmsDistances(groupView, annulus, magRange=None):
+    """Calculate the RMS distance of a set of matched objects over visits.
+
+    Parameters
+    ----------
+    groupView : lsst.afw.table.GroupView
+        GroupView object of matched observations from MultiMatch.
+    annulus : 2-element list or tuple of float
+        Distance range in which to compare object.  [arcmin]
+        E.g., `annulus=[19, 21]` would consider all objects
+        separated from each other between 19 and 21 arcminutes.
+    magRange : 2-element list or tuple of float, optional
+        Magnitude range from which to select objects.
+        Default of `None` will result in all objects being considered.
+
+    Returns
+    -------
+    list
+        rmsDistMas
+
+    """
+
     # Default is specified here separately because defaults that are mutable
     # get overridden by previous calls of the function.
     if magRange is None:
@@ -469,7 +589,6 @@ def magInRange(cat):
     meanRa = groupViewInMagRange.aggregate(averageRaFromCat)
     meanDec = groupViewInMagRange.aggregate(averageDecFromCat)
 
-    annulus = D + (width/2)*np.array([-1, +1])
     annulusRadians = arcminToRadians(annulus)
 
     rmsDistances = list()
@@ -488,6 +607,4 @@ def magInRange(cat):
                 rmsDist = np.std(np.array(distances)[finiteEntries])
                 rmsDistances.append(rmsDist)
 
-    rmsDistMas = radiansToMilliarcsec(rmsDistances)
-
-    return rmsDistMas, annulus, magRange
+    return rmsDistances, annulus, magRange