DM-5660: Add motivated model fits to validate_drp photometric and astrometric scatter/repeatability analysis and plots

README.md

@@ -1,20 +1,21 @@
 Validate an LSST DM processCcd.py output repository
 against a set of LSST Science Requirements Document Key Performance Metrics.
+and expected performance following the LSST Overview paper.
 
 ```
 setup validate_drp
 validateDrp.py CFHT/output
 ```
 
-will produces output, plots, JSON files analyzing the processed data in `CFHT/output`.  Metrics will be separately calculated for each filter represented in the repository.  Replace `CFHT/output` with your favorite processed data repository and you will get reasonable output.
+will produces output, plots, JSON files analyzing the processed data in `CFHT/output`.  Metrics will be separately calculated for each filter represented in the repository.  
+
+Replace `CFHT/output` with your favorite processed data repository and you will get reasonable output.
 
 One can run `validateDrp.py` in any of the following modes:
 1. use no configuration file (as above)
 2. pass a configuration file with just validation parameters (brightSnr, number of expected matches, ...) but no dataId specifications
 3. pass a configuration file that specifies validation parameters and the dataIds to process.  See examples below for use with a `--configFile`
 
-Caveat:  Will likely not successfully run on more than 500 catalogs per band due to memory limits and inefficiencies in the current matching approach.
-
 ------
 This package also includes examples that run processCcd task on some 
 CFHT data and DECam data
@@ -145,6 +146,14 @@ Once these basic steps are completed, then you can run any of the following:
 
 Note that the example validation test selects several of the CCDs and will fail if you just pass it a repository with 1 visit or just 1 CCD.
 
+Notes
+-----
+* Will likely not successfully run on more than 500 catalogs per band due to memory limits and inefficiencies in the current matching approach.
+* The astrometric and photometric error models are formally valid for individual images.  However, they are being applied here to the results from the set of images, which is implicitly looking at some sort of mean performance.
+E.g., the expected astrometric uncertainty is intimately related to the seeing of the image.  For collections of images where most have a similar seeing, these estimates are useful and reasonable.  However, if the data set analyzed consisted of a set of images distributed across a wide range of seeing values, then the fits here have less direct meaning.
+
+
+
 Files of Interest:
 ------------------
 * `bin.src/validateDrp.py`  : Analyze output data produced by processCcd.py

python/lsst/validate/drp/check.py

@@ -21,8 +21,10 @@
 from __future__ import print_function, division
 
 import numpy as np
+from scipy.optimize import curve_fit
 
 import lsst.afw.geom as afwGeom
+import lsst.pipe.base as pipeBase
 
 from .util import averageRaDecFromCat
 
@@ -65,6 +67,44 @@ def magNormDiff(cat):
     return normDiff
 
 
+def fitExp(x, y, y_err, deg=2):
+    """Fit an exponential quadratic to x, y, y_err.
+    """
+    fit_params, fit_param_covariance = \
+        curve_fit(expModel, x, y, p0=[1, 0.02, 5], sigma=y_err)
+
+    return fit_params
+
+
+def fitAstromErrModel(snr, dist):
+    fit_params, fit_param_covariance = \
+        curve_fit(astromErrModel, snr, dist, p0=[1, 0.01])
+
+    return fit_params
+
+
+def fitPhotErrModel(mag, mmag_err):
+    """Fit model of photometric error from LSST Overview paper
+
+    Parameters
+    ----------
+    mag : list or numpy.array
+        Magnitude
+    mmag_err : list or numpy.array
+        Magnitude uncertainty or variation in *mmag*.
+
+    Returns
+    -------
+    float, float, float
+        sigmaSys, gamma, m5 fit parameters.
+    """
+    mag_err = mmag_err / 1000
+    fit_params, fit_param_covariance = \
+        curve_fit(photErrModel, mag, mag_err, p0=[0.01, 0.039, 24.35])
+
+    return fit_params
+
+
 def positionRms(cat):
     """Calculate the RMS for RA, Dec for a set of observations an object.
 
@@ -114,8 +154,11 @@ def checkAstrometry(snr, dist, match,
 
     Returns
     -------
-    float
-        The astrometric scatter (RMS, milliarcsec) for all good stars.
+    pipeBase.Struct
+        name -- str: "checkAstrometry"
+        astromScatter -- float: The astrometric scatter (RMS, milliarcsec) for all good stars.
+        astromErrModel -- str: Description of astrometric error model
+        astromErrFitParams -- list or numpy.array of float:  Fit parameters.
 
     Notes
     -----
@@ -132,16 +175,21 @@ def checkAstrometry(snr, dist, match,
     print("Astrometric scatter (median) - snr > %.1f : %.1f %s" %
           (brightSnr, astromScatter, "mas"))
 
+    fit_params = fitAstromErrModel(snr[bright], dist[bright])
+
     if astromScatter > medianRef:
         print("Median astrometric scatter %.1f %s is larger than reference : %.1f %s " %
               (astromScatter, "mas", medianRef, "mas"))
     if match < matchRef:
         print("Number of matched sources %d is too small (shoud be > %d)" % (match, matchRef))
 
-    return astromScatter
+    return pipeBase.Struct(name="checkAstrometry",
+                           astromRmsScatter=astromScatter, 
+                           astromErrModel=astromErrModel.__doc__,
+                           astromErrFitParams=fit_params)
 
 
-def checkPhotometry(snr, mag, mmagrms, dist, match,
+def checkPhotometry(snr, mag, mmagErr, mmagrms, dist, match,
                     brightSnr=100,
                     medianRef=100, matchRef=500):
     """Print out the astrometric scatter for all stars, and for good stars.
@@ -151,9 +199,11 @@ def checkPhotometry(snr, mag, mmagrms, dist, match,
     snr : list or numpy.array
         Median SNR of PSF flux
     mag : list or numpy.array
-        Average magnitudes of each star
+        Average magnitudes of each star.  [mag]
+    mmagErr : list or numpy.array
+        Uncertainties in magnitudes of each star.  [mmag]
     mmagrms ; list or numpy.array
-        Magnitude RMS of the multiple observation of each star.
+        Magnitude RMS of the multiple observation of each star. [mmag]
     dist : list or numpy.array
         Distances between successive measurements of one star
     match : int
@@ -167,8 +217,10 @@ def checkPhotometry(snr, mag, mmagrms, dist, match,
 
     Returns
     -------
-    float
-        The photometry scatter (RMS, millimag) for all good stars.
+    pipeBase.Struct
+        name -- str: "checkPhotometry"
+        photScatter -- float: The photometric scatter (RMS, mmag) for all good star stars.
+        photFitParams -- list or numpy.array of float:  Fit parameters.
 
     Notes
     -----
@@ -181,15 +233,80 @@ def checkPhotometry(snr, mag, mmagrms, dist, match,
           (np.median(mmagrms), "mmag"))
 
     bright = np.where(np.asarray(snr) > brightSnr)
-    photoScatter = np.median(np.asarray(mmagrms)[bright])
+    photScatter = np.median(np.asarray(mmagrms)[bright])
     print("Photometric scatter (median) - SNR > %.1f : %.1f %s" %
-          (brightSnr, photoScatter, "mmag"))
+          (brightSnr, photScatter, "mmag"))
+
+    fit_params = fitPhotErrModel(mag[bright], mmagErr[bright])
 
-    if photoScatter > medianRef:
+    if photScatter > medianRef:
         print("Median photometric scatter %.3f %s is larger than reference : %.3f %s "
-              % (photoScatter, "mmag", medianRef, "mag"))
+              % (photScatter, "mmag", medianRef, "mag"))
     if match < matchRef:
         print("Number of matched sources %d is too small (shoud be > %d)"
               % (match, matchRef))
 
-    return photoScatter
+    return pipeBase.Struct(name="checkPhotometry",
+                           photRmsScatter=photScatter, 
+                           photErrModel=photErrModel.__doc__,
+                           photErrFitParams=fit_params)
+
+
+def astromErrModel(snr, theta=1, sigma_sys=0.01, C=1):
+    """Calculate expected astrometric uncertainty based on SNR.
+
+    mas = C*theta/SNR
+
+    Parameters
+    ----------
+    snr : list or numpy.array
+        S/N of photometric measurements
+    theta : float or numpy.array, optional
+        Seeing [mas]
+    sigma_sys : float
+        Systematic error floor [mas]
+    C : float
+        Scaling factor
+    verbose : bool, optional
+        Produce extra output to STDOUT
+
+    Returns
+    -------
+    np.array
+        Expected astrometric uncertainty. [mas]
+    """
+    return C*theta/snr + sigma_sys
+
+
+def photErrModel(mag, sigmaSys, gamma, m5):
+    """Fit model of photometric error from LSST Overview paper
+    http://arxiv.org/abs/0805.2366v4
+
+    Photometric errors described by
+    Eq. 4
+    sigma_1^2 = sigma_sys^2 + sigma_rand^2
+
+    Eq. 5
+    sigma_(rand) = (0.04 - gamma) * x + gamma * x^2 [mag^2]
+    where x = 10**(0.4*(m-m_5))
+
+    Parameters
+    ----------
+    mag : list or numpy.array
+        Magnitude
+    sigmaSq : float
+        Limiting systematics floor [mag]
+    gamma : float
+        proxy for sky brightness and readout noise
+    m5 : float
+        5-sigma depth [mag]
+
+    Returns
+    -------
+    numpy.array
+        Result of noise estimation function
+    """
+    x = 10**(0.4*(mag - m5))
+    sigmaRandSq = (0.04 - gamma) * x + gamma * x**2
+    sigmaSq = sigmaSys**2 + sigmaRandSq
+    return np.sqrt(sigmaSq)