Add param for PTC turn-off search

Fix bug in config file

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

@@ -27,6 +27,8 @@
 from lsst.cp.pipe.utils import (fitLeastSq, fitBootstrap, funcPolynomial, funcAstier)
 
 from scipy.optimize import least_squares
+from itertools import groupby
+from operator import itemgetter
 
 import lsst.pipe.base.connectionTypes as cT
 
@@ -109,10 +111,11 @@ class PhotonTransferCurveSolveConfig(pipeBase.PipelineTaskConfig,
         doc="Sigma cut for outlier rejection in PTC.",
         default=5.0,
     )
-    maxIterationsPtcOutliers = pexConfig.Field(
+    maxIterationsPtcOutliers = pexConfig.RangeField(
         dtype=int,
         doc="Maximum number of iterations for outlier rejection in PTC.",
         default=2,
+        min=0
     )
     minVarPivotSearch = pexConfig.Field(
         dtype=float,
@@ -122,6 +125,13 @@ class PhotonTransferCurveSolveConfig(pipeBase.PipelineTaskConfig,
             " should be sought.",
         default=10000,
     )
+    consecutivePointsVarDecreases = pexConfig.RangeField(
+        dtype=int,
+        doc="Required number of consecutive points/fluxes in the PTC where the variance "
+            "decreases in order to find a first estimate of the PTC turn-off. ",
+        default=2,
+        min=2
+    )
     doFitBootstrap = pexConfig.Field(
         dtype=bool,
         doc="Use bootstrap for the PTC fit parameters and errors?.",
@@ -395,7 +405,7 @@ def _boundsForAstier(initialPars, lowers=[], uppers=[]):
         return (lowers, uppers)
 
     @staticmethod
-    def _getInitialGoodPoints(means, variances, minVarPivotSearch):
+    def _getInitialGoodPoints(means, variances, minVarPivotSearch, consecutivePointsVarDecreases):
         """Return a boolean array to mask bad points.
 
         Parameters
@@ -407,6 +417,11 @@ def _getInitialGoodPoints(means, variances, minVarPivotSearch):
         minVarPivotSearch : `float`
             The variance (in ADU^2), above which, the point
             of decreasing variance should be sought.
+        consecutivePointsVarDecreases : `int`
+            Required number of consecutive points/fluxes
+            in the PTC where the variance
+            decreases in order to find a first
+            estimate of the PTC turn-off.
 
         Returns
         ------
@@ -416,20 +431,27 @@ def _getInitialGoodPoints(means, variances, minVarPivotSearch):
 
         Notes
         -----
-        Eliminate points beyond which the variance decreases
+        Eliminate points beyond which the variance decreases.
         """
         goodPoints = np.ones_like(means, dtype=bool)
         pivotList = np.where(np.array(np.diff(variances)) < 0)[0]
         if len(pivotList) > 0:
             # For small values, sometimes the variance decreases slightly
             # Only look when var > self.config.minVarPivotSearch
             pivotList = [p for p in pivotList if variances[p] > minVarPivotSearch]
+            # Require that the varince decreases during
+            # consecutivePointsVarDecreases
+            # consecutive points. This will give a first
+            # estimate of the PTC turn-off, which
+            # may be updated (reduced) further in the code.
             if len(pivotList) > 1:
-                # Require that the decrease in variance happen for two
-                # consecutive signal levels
-                pivotIndex = np.min(np.where(np.diff(pivotList) == 1)[0])
-                pivot = pivotList[pivotIndex]
-                goodPoints[pivot+1:] = False
+                for k, g in groupby(enumerate(pivotList), lambda x: x[0]-x[1]):
+                    group = (map(itemgetter(1), g))
+                    group = list(map(int, group))
+                    if len(group) >= consecutivePointsVarDecreases:
+                        pivotIndex = np.min(group)
+                        goodPoints[pivotIndex+1:] = False
+                        break
 
         return goodPoints
 
@@ -518,7 +540,8 @@ def errFunc(p, x, y):
             # Discard points when the variance starts to decrease after two
             # consecutive signal levels
             goodPoints = self._getInitialGoodPoints(meanVecOriginal, varVecOriginal,
-                                                    self.config.minVarPivotSearch)
+                                                    self.config.minVarPivotSearch,
+                                                    self.config.consecutivePointsVarDecreases)
             # Check if all points are bad from the 'cpExtractPtcTask'
             initialExpIdMask = np.ravel(np.array(dataset.expIdMask[ampName]))
 
@@ -545,8 +568,14 @@ def errFunc(p, x, y):
                 parsIniPtc = self._initialParsForPolynomial(self.config.polynomialFitDegree + 1)
                 bounds = self._boundsForPolynomial(parsIniPtc)
 
-            # Before bootstrap fit, do an iterative fit to get rid of outliers
+            # Before bootstrap fit, do an iterative fit to get rid of outliers.
+            # This further process of outlier rejection be skipped
+            # if self.config.maxIterationsPtcOutliers = 0.
+            # We already did some initial outlier rejection about in
+            # self._getInitialGoodPoints.
             count = 1
+            newMask = np.ones_like(meanVecOriginal)
+            pars = parsIniPtc
             while count <= maxIterationsPtcOutliers:
                 # Note that application of the mask actually shrinks the array
                 # to size rather than setting elements to zero (as we want) so

tests/test_ptc.py

@@ -377,12 +377,14 @@ def test_makeZeroSafe(self):
     def test_getInitialGoodPoints(self):
         xs = [1, 2, 3, 4, 5, 6]
         ys = [2*x for x in xs]
-        points = self.defaultTaskSolve._getInitialGoodPoints(xs, ys, minVarPivotSearch=0.)
+        points = self.defaultTaskSolve._getInitialGoodPoints(xs, ys, minVarPivotSearch=0.,
+                                                             consecutivePointsVarDecreases=2)
         assert np.all(points) == np.all(np.array([True for x in xs]))
 
         ys[4] = 7  # Variance decreases in two consecutive points after ys[3]=8
         ys[5] = 6
-        points = self.defaultTaskSolve._getInitialGoodPoints(xs, ys, minVarPivotSearch=0.)
+        points = self.defaultTaskSolve._getInitialGoodPoints(xs, ys, minVarPivotSearch=0.,
+                                                             consecutivePointsVarDecreases=2)
         assert np.all(points) == np.all(np.array([True, True, True, True, False]))
 
     def test_getExpIdsUsed(self):