Merge pull request #275 from lsst/tickets/DM-46750

DM-46750: Add linearityTurnoff and linearityMaxSignal fitting to cpLinearizerSolve.

python/lsst/cp/pipe/cpLinearitySolve.py

@@ -179,12 +179,12 @@ class LinearitySolveConfig(pipeBase.PipelineTaskConfig,
     )
     maxLinearAdu = pexConfig.Field(
         dtype=float,
-        doc="Maximum DN value to use to estimate linear term.",
+        doc="Maximum adu value to use to estimate linear term; not used with spline fits.",
         default=20000.0,
     )
     minLinearAdu = pexConfig.Field(
         dtype=float,
-        doc="Minimum DN value to use to estimate linear term.",
+        doc="Minimum adu value to use to estimate linear term.",
         default=30.0,
     )
     nSigmaClipLinear = pexConfig.Field(
@@ -196,6 +196,19 @@ class LinearitySolveConfig(pipeBase.PipelineTaskConfig,
         dtype=bool,
         doc="Ignore the expIdMask set by the PTC solver?",
         default=False,
+        deprecated="This field is no longer used. Will be removed after v28.",
+    )
+    maxFracLinearityDeviation = pexConfig.Field(
+        dtype=float,
+        doc="Maximum fraction deviation from raw linearity to compute "
+            "linearityTurnoff and linearityMaxSignal.",
+        # TODO: DM-46811 investigate if this can be raised to 0.05.
+        default=0.01,
+    )
+    minSignalFitLinearityTurnoff = pexConfig.Field(
+        dtype=float,
+        doc="Minimum signal to compute raw linearity slope for linearityTurnoff.",
+        default=1000.0,
     )
     usePhotodiode = pexConfig.Field(
         dtype=bool,
@@ -444,12 +457,8 @@ def run(self, inputPtc, dummy, camera, inputDims,
                     "If you really know what you are doing, you may reduce "
                     "config.splineGroupingMinPoints.")
 
-            if (len(inputPtc.expIdMask[ampName]) == 0) or self.config.ignorePtcMask:
-                self.log.warning("Mask not found for %s in detector %s in fit. Using all points.",
-                                 ampName, detector.getName())
-                mask = np.ones(len(inputPtc.expIdMask[ampName]), dtype=bool)
-            else:
-                mask = inputPtc.expIdMask[ampName].copy()
+            # We start with all finite values.
+            mask = np.isfinite(inputPtc.rawMeans[ampName])
 
             if self.config.linearityType == "Spline" and temperatureValues is not None:
                 mask &= np.isfinite(temperatureValues)
@@ -477,9 +486,25 @@ def run(self, inputPtc, dummy, camera, inputDims,
             else:
                 inputAbscissa = inputPtc.rawExpTimes[ampName].copy()
 
+            # Compute linearityTurnoff and linearitySignalMax.
+            turnoff, maxSignal = self._computeTurnoffAndMax(
+                inputAbscissa[mask],
+                inputPtc.rawMeans[ampName][mask],
+                inputPtc.expIdMask[ampName][mask],
+                ampName=ampName,
+            )
+            linearizer.linearityTurnoff[ampName] = turnoff
+            linearizer.linearityMaxSignal[ampName] = maxSignal
+
             inputOrdinate = inputPtc.rawMeans[ampName].copy()
 
-            mask &= (inputOrdinate < self.config.maxLinearAdu)
+            if self.config.linearityType != 'Spline':
+                mask &= (inputOrdinate < self.config.maxLinearAdu)
+            else:
+                # For spline fits, cut above the turnoff.
+                self.log.info("Using linearityTurnoff of %.4f adu for amplifier %s", turnoff, ampName)
+                mask &= (inputOrdinate <= turnoff)
+
             mask &= (inputOrdinate > self.config.minLinearAdu)
 
             if mask.sum() < 2:
@@ -579,7 +604,7 @@ def run(self, inputPtc, dummy, camera, inputDims,
                 nodes = np.linspace(0.0, np.max(inputOrdinate[mask]), self.config.splineKnots)
 
                 if temperatureValues is not None:
-                    temperatureValuesScaled = temperatureValues - np.median(temperatureValues[~mask])
+                    temperatureValuesScaled = temperatureValues - np.median(temperatureValues[mask])
                 else:
                     temperatureValuesScaled = None
 
@@ -595,6 +620,7 @@ def run(self, inputPtc, dummy, camera, inputDims,
                     weight_pars_start=self.config.splineFitWeightParsStart,
                     fit_temperature=self.config.doSplineFitTemperature,
                     temperature_scaled=temperatureValuesScaled,
+                    max_signal_nearly_linear=inputPtc.ptcTurnoff[ampName],
                 )
                 p0 = fitter.estimate_p0()
                 pars = fitter.fit(
@@ -756,6 +782,8 @@ def fillBadAmp(self, linearizer, fitOrder, inputPtc, amp):
         linearizer.fitResiduals[ampName] = np.zeros(len(inputPtc.expIdMask[ampName]))
         linearizer.fitResidualsSigmaMad[ampName] = np.nan
         linearizer.linearFit[ampName] = np.zeros(2)
+        linearizer.linearityTurnoff[ampName] = np.nan
+        linearizer.linearityMaxSignal[ampName] = np.nan
         return linearizer
 
     def fixupBadAmps(self, linearizer):
@@ -778,6 +806,65 @@ def fixupBadAmps(self, linearizer):
             elif len(linearizer.fitParams[ampName]) != fitParamsMaxLen:
                 raise RuntimeError("Linearity has mismatched fitParams; check code/data.")
 
+    def _computeTurnoffAndMax(self, abscissa, ordinate, initialMask, ampName="UNKNOWN"):
+        """Compute the turnoff and max signal.
+
+        Parameters
+        ----------
+        abscissa : `np.ndarray`
+            Input x values, either photoCharges or exposure times.
+            These should be cleaned of any non-finite values.
+        ordinate : `np.ndarray`
+            Input y values, the raw mean values for the amp.
+            These should be cleaned of any non-finite values.
+        initialMask : `np.ndarray`
+            Mask to use for initial fit (usually from PTC).
+        ampName : `str`, optional
+            Amplifier name (used for logging).
+
+        Returns
+        -------
+        turnoff : `float`
+            Fit turnoff value.
+        maxSignal : `float`
+            Fit maximum signal value.
+        """
+        # Follow eo_pipe:
+        # https://github.com/lsst-camera-dh/eo_pipe/blob/6afa546569f622b8d604921e248200481c445730/python/lsst/eo/pipe/linearityPlotsTask.py#L50
+        # Replacing flux with abscissa, Ne with ordinate.
+
+        # Fit a line with the y-intercept fixed to zero, using the
+        # signal counts Ne as the variance in the chi-square, i.e.,
+        # chi2 = sum( (ordinate - aa*abscissa)**2/ordinate )
+        # Minimizing chi2 wrt aa, gives
+        # aa = sum(abscissa) / sum(abscissa**2/ordinate)
+
+        fitMask = initialMask.copy()
+        fitMask[ordinate < self.config.minSignalFitLinearityTurnoff] = False
+
+        aa = sum(abscissa[fitMask])/sum(abscissa[fitMask]**2/ordinate[fitMask])
+
+        # Use the residuals to compute the turnoff.
+        residuals = (ordinate - aa*abscissa)/ordinate
+        turnoffIndex = np.argmax(ordinate[np.abs(residuals) < self.config.maxFracLinearityDeviation])
+        turnoff = ordinate[turnoffIndex]
+
+        # Fit the maximum signal.
+        if turnoffIndex == (len(residuals) - 1):
+            # This is the last point; we can't do a fit.
+            self.log.warning(
+                "No linearity turnoff detected for amplifier %s; try to increase the signal range.",
+                ampName,
+            )
+            maxSignal = ordinate[turnoffIndex]
+        else:
+            maxSignalInitial = np.max(ordinate)
+
+            highFluxPoints = (ordinate > (1.0 - self.config.maxFracLinearityDeviation)*maxSignalInitial)
+            maxSignal = np.median(ordinate[highFluxPoints])
+
+        return turnoff, maxSignal
+
     def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName):
         """Debug method for linearity fitting.
 

python/lsst/cp/pipe/utils.py

@@ -973,6 +973,12 @@ class AstierSplineLinearityFitter:
         Fit for temperature scaling?
     temperature_scaled : `np.ndarray` (M,), optional
         Input scaled temperature values (T - T_ref).
+    max_signal_nearly_linear : `float`, optional
+        Maximum signal that we are confident the input is nearly
+        linear. This is used both for regularization, and for
+        fitting the raw slope. Usually set to the ptc turnoff,
+        above which we allow the spline to significantly deviate
+        and do not demand the deviation to average to zero.
     """
     def __init__(
         self,
@@ -987,6 +993,7 @@ def __init__(
         weight_pars_start=[1.0, 0.0],
         fit_temperature=False,
         temperature_scaled=None,
+        max_signal_nearly_linear=None,
     ):
         self._pd = pd
         self._mu = mu
@@ -1026,7 +1033,10 @@ def __init__(
         self._temperature_scaled = temperature_scaled
 
         # Values to regularize spline fit.
-        self._x_regularize = np.linspace(0.0, self._mu[self.mask].max(), 100)
+        if max_signal_nearly_linear is None:
+            max_signal_nearly_linear = self._mu[self.mask].max()
+        self._max_signal_nearly_linear = max_signal_nearly_linear
+        self._x_regularize = np.linspace(0.0, self._max_signal_nearly_linear, 100)
 
         # Set up the indices for the fit parameters.
         self.par_indices = {
@@ -1083,7 +1093,8 @@ def estimate_p0(self):
         p0 = np.zeros(npt)
 
         # Do a simple linear fit and set all the constants to this.
-        linfit = np.polyfit(self._pd[self.mask], self._mu[self.mask], 1)
+        to_fit = (self._mu[self.mask] < self._max_signal_nearly_linear)
+        linfit = np.polyfit(self._pd[self.mask][to_fit], self._mu[self.mask][to_fit], 1)
         p0[self.par_indices["groups"]] = linfit[0]
 
         # Look at the residuals...