DM-24762: Add option to PTC task to correct for sigma clipping bias

python/lsst/cp/pipe/measureCrosstalk.py

@@ -20,7 +20,6 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 import itertools
 import numpy as np
-from scipy.stats import norm
 
 from collections import defaultdict
 
@@ -33,7 +32,7 @@
 from lsst.pex.config import Config, Field, ListField, ConfigurableField
 from lsst.ip.isr import CrosstalkCalib, IsrProvenance
 from lsst.pipe.tasks.getRepositoryData import DataRefListRunner
-from lsst.cp.pipe.utils import ddict2dict
+from lsst.cp.pipe.utils import (ddict2dict, sigmaClipCorrection)
 
 from ._lookupStaticCalibration import lookupStaticCalibration
 
@@ -604,7 +603,7 @@ def measureCrosstalkCoefficients(self, ratios, rejIter, rejSigma):
                 if calib.coeffNum[ii][jj] == 1:
                     calib.coeffErr[ii][jj] = np.nan
                 else:
-                    correctionFactor = self.sigmaClipCorrection(rejSigma)
+                    correctionFactor = sigmaClipCorrection(rejSigma)
                     calib.coeffErr[ii][jj] = np.std(values) * correctionFactor
                 calib.coeffValid[ii][jj] = (np.abs(calib.coeffs[ii][jj])
                                             > calib.coeffErr[ii][jj] / np.sqrt(calib.coeffNum[ii][jj]))
@@ -615,33 +614,6 @@ def measureCrosstalkCoefficients(self, ratios, rejIter, rejSigma):
 
         return calib
 
-    @staticmethod
-    def sigmaClipCorrection(nSigClip):
-        """Correct measured sigma to account for clipping.
-
-        If we clip our input data and then measure sigma, then the
-        measured sigma is smaller than the true value because real
-        points beyond the clip threshold have been removed.  This is a
-        small (1.5% at nSigClip=3) effect when nSigClip >~ 3, but the
-        default parameters for measure crosstalk use nSigClip=2.0.
-        This causes the measured sigma to be about 15% smaller than
-        real.  This formula corrects the issue, for the symmetric case
-        (upper clip threshold equal to lower clip threshold).
-
-        Parameters
-        ----------
-        nSigClip : `float`
-            Number of sigma the measurement was clipped by.
-
-        Returns
-        -------
-        scaleFactor : `float`
-            Scale factor to increase the measured sigma by.
-
-        """
-        varFactor = 1.0 + (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))
-        return 1.0 / np.sqrt(varFactor)
-
     @staticmethod
     def filterCrosstalkCalib(inCalib):
         """Apply valid constraints to the measured values.

python/lsst/cp/pipe/ptc/astierCovPtcFit.py

@@ -115,12 +115,6 @@ def makeCovArray(inputTuple, maxRangeFromTuple=8):
     cov[ind, i, j] = c
     var[ind, i, j] = v**2/n
     var[:, 0, 0] *= 2  # var(v) = 2*v**2/N
-    # compensate for loss of variance and covariance due to outlier elimination(sigma clipping)
-    # when computing variances(cut to 4 sigma): 1 per mill for variances and twice as
-    # much for covariances:
-    fact = 1.0  # 1.00107
-    cov *= fact*fact
-    cov[:, 0, 0] /= fact
 
     return cov, var, muVals
 

python/lsst/cp/pipe/ptc/cpExtractPtcTask.py

@@ -24,7 +24,8 @@
 import lsst.afw.math as afwMath
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
-from lsst.cp.pipe.utils import arrangeFlatsByExpTime, arrangeFlatsByExpId
+from lsst.cp.pipe.utils import (arrangeFlatsByExpTime, arrangeFlatsByExpId,
+                                sigmaClipCorrection)
 
 import lsst.pipe.base.connectionTypes as cT
 
@@ -279,6 +280,11 @@ def run(self, inputExp, inputDims):
                 # in {maxLag, maxLag}^2]
                 muDiff, varDiff, covAstier = self.measureMeanVarCov(exp1, exp2, region=region,
                                                                     covAstierRealSpace=doRealSpace)
+                # Correction factor for sigma clipping. Function returns 1/sqrt(varFactor),
+                # so it needs to be squared. varDiff is calculated via afwMath.VARIANCECLIP.
+                varFactor = sigmaClipCorrection(self.config.nSigmaClipPtc)**2
+                varDiff *= varFactor
+
                 expIdMask = True
                 if np.isnan(muDiff) or np.isnan(varDiff) or (covAstier is None):
                     msg = (f"NaN mean or var, or None cov in amp {ampName} in exposure pair {expId1},"
@@ -301,6 +307,12 @@ def run(self, inputExp, inputDims):
                                                      self.config.maximumRangeCovariancesAstier)
                     covSqrtWeights = np.nan_to_num(1./np.sqrt(vcov))
 
+                # Correct covArray for sigma clipping:
+                # 1) Apply varFactor twice for the whole covariance matrix
+                covArray *= varFactor**2
+                # 2) But, only once for the variance element of the matrix, covArray[0,0]
+                covArray[0, 0] /= varFactor
+
                 partialPtcDataset.setAmpValues(ampName, rawExpTime=[expTime], rawMean=[muDiff],
                                                rawVar=[varDiff], inputExpIdPair=[(expId1, expId2)],
                                                expIdMask=[expIdMask], covArray=covArray,

python/lsst/cp/pipe/utils.py

@@ -28,6 +28,7 @@
 import numpy as np
 from scipy.optimize import leastsq
 import numpy.polynomial.polynomial as poly
+from scipy.stats import norm
 
 import lsst.pipe.base as pipeBase
 import lsst.ip.isr as ipIsr
@@ -38,6 +39,33 @@
 import galsim
 
 
+def sigmaClipCorrection(nSigClip):
+    """Correct measured sigma to account for clipping.
+
+    If we clip our input data and then measure sigma, then the
+    measured sigma is smaller than the true value because real
+    points beyond the clip threshold have been removed.  This is a
+    small (1.5% at nSigClip=3) effect when nSigClip >~ 3, but the
+    default parameters for measure crosstalk use nSigClip=2.0.
+    This causes the measured sigma to be about 15% smaller than
+    real.  This formula corrects the issue, for the symmetric case
+    (upper clip threshold equal to lower clip threshold).
+
+    Parameters
+    ----------
+    nSigClip : `float`
+        Number of sigma the measurement was clipped by.
+
+    Returns
+    -------
+    scaleFactor : `float`
+        Scale factor to increase the measured sigma by.
+
+    """
+    varFactor = 1.0 - (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))
+    return 1.0 / np.sqrt(varFactor)
+
+
 def calculateWeightedReducedChi2(measured, model, weightsMeasured, nData, nParsModel):
     """Calculate weighted reduced chi2.