DM-27185: ptc.py fails with ptcFitType=FULLCOVARIANCE

python/lsst/cp/pipe/astierCovPtcFit.py

@@ -27,6 +27,7 @@
 from scipy.optimize import leastsq
 from .astierCovFitParameters import FitParameters
 
+import lsst.log as lsstLog
 
 __all__ = ["CovFit"]
 
@@ -237,8 +238,10 @@ class CovFit:
 
     def __init__(self, inputTuple, maxRangeFromTuple=8):
         self.cov, self.vcov, self.mu = makeCovArray(inputTuple, maxRangeFromTuple)
-        self.sqrtW = 1./np.sqrt(self.vcov)
+        # make it nan safe, replacing nan's with 0 in weights
+        self.sqrtW = np.nan_to_num(1./np.sqrt(self.vcov))
         self.r = self.cov.shape[1]
+        self.logger = lsstLog.Log.getDefaultLogger()
 
     def subtractDistantOffset(self, maxLag=8, startLag=5, polDegree=1):
         """Subtract a background/offset to the measured covariances.
@@ -507,7 +510,7 @@ def setAandB(self, a, b):
         return
 
     def chi2(self):
-        """Calculate weighte chi2 of full-model fit."""
+        """Calculate weighted chi2 of full-model fit."""
         return(self.weightedRes()**2).sum()
 
     def wres(self, params=None):
@@ -570,6 +573,7 @@ def fitFullModel(self, pInit=None, nSigma=5.0, maxFitIter=3):
         nOutliers = 1
         counter = 0
         while nOutliers != 0:
+            # If fit fails, None values are returned and caught in getOutputPtcDataCovAstier of ptc.py
             params, paramsCov, _, mesg, ierr = leastsq(self.weightedRes, pInit, full_output=True)
             wres = self.weightedRes(params)
             # Do not count the outliers as significant
@@ -579,10 +583,9 @@ def fitFullModel(self, pInit=None, nSigma=5.0, maxFitIter=3):
             nOutliers = mask.sum()
             counter += 1
             if counter == maxFitIter:
+                self.logger.warn(f"Max number of iterations,{maxFitIter}, reached in fitFullModel")
                 break
 
-        if ierr not in [1, 2, 3, 4]:
-            raise RuntimeError("Minimization failed: " + mesg)
         self.covParams = paramsCov
         return params
 

python/lsst/cp/pipe/astierCovPtcUtils.py

@@ -436,5 +436,4 @@ def getFitDataFromCovariances(i, j, mu, fullCov, fullCovModel, fullCovSqrtWeight
         covariance /= mu
         covarianceModel /= mu
         weights *= mu
-
     return mu, covariance, covarianceModel, weights, mask

python/lsst/cp/pipe/linearity.py

@@ -30,7 +30,7 @@
 from lsstDebug import getDebugFrame
 from lsst.ip.isr import (Linearizer, IsrProvenance)
 
-from .utils import (fitLeastSq, funcPolynomial)
+from .utils import (funcPolynomial, irlsFit)
 
 
 __all__ = ["LinearitySolveTask", "LinearitySolveConfig", "MeasureLinearityTask"]
@@ -205,8 +205,8 @@ def run(self, inputPtc, camera, inputDims):
             linearAbscissa = inputAbscissa[fluxMask]
             linearOrdinate = inputOrdinate[fluxMask]
 
-            linearFit, linearFitErr, chiSq, weights = self.irlsFit([0.0, 100.0], linearAbscissa,
-                                                                   linearOrdinate, funcPolynomial)
+            linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 100.0], linearAbscissa,
+                                                              linearOrdinate, funcPolynomial)
             # Convert this proxy-to-flux fit into an expected linear flux
             linearOrdinate = linearFit[0] + linearFit[1] * inputAbscissa
 
@@ -221,8 +221,8 @@ def run(self, inputPtc, camera, inputDims):
             if self.config.linearityType in ['Polynomial', 'Squared', 'LookupTable']:
                 polyFit = np.zeros(fitOrder + 1)
                 polyFit[1] = 1.0
-                polyFit, polyFitErr, chiSq, weights = self.irlsFit(polyFit, linearOrdinate,
-                                                                   fitOrdinate, funcPolynomial)
+                polyFit, polyFitErr, chiSq, weights = irlsFit(polyFit, linearOrdinate,
+                                                              fitOrdinate, funcPolynomial)
 
                 # Truncate the polynomial fit
                 k1 = polyFit[1]
@@ -329,50 +329,6 @@ def run(self, inputPtc, camera, inputDims):
             outputProvenance=provenance,
         )
 
-    def irlsFit(self, initialParams, dataX, dataY, function, weightsY=None):
-        """Iteratively reweighted least squares fit.
-
-        This uses the `lsst.cp.pipe.utils.fitLeastSq`, but applies
-        weights based on the Cauchy distribution to the fitter.  See
-        e.g. Holland and Welsch, 1977, doi:10.1080/03610927708827533
-
-        Parameters
-        ----------
-        initialParams : `list` [`float`]
-            Starting parameters.
-        dataX : `numpy.array` [`float`]
-            Abscissa data.
-        dataY : `numpy.array` [`float`]
-            Ordinate data.
-        function : callable
-            Function to fit.
-        weightsY : `numpy.array` [`float`]
-            Weights to apply to the data.
-
-        Returns
-        -------
-        polyFit : `list` [`float`]
-            Final best fit parameters.
-        polyFitErr : `list` [`float`]
-            Final errors on fit parameters.
-        chiSq : `float`
-            Reduced chi squared.
-        weightsY : `list` [`float`]
-            Final weights used for each point.
-
-        """
-        if not weightsY:
-            weightsY = np.ones_like(dataX)
-
-        polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
-        for iteration in range(10):
-            # Use Cauchy weights
-            resid = np.abs(dataY - function(polyFit, dataX)) / np.sqrt(dataY)
-            weightsY = 1.0 / (1.0 + np.sqrt(resid / 2.385))
-            polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
-
-        return polyFit, polyFitErr, chiSq, weightsY
-
     def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName):
         """Debug method for linearity fitting.
 

python/lsst/cp/pipe/plotPtc.py

@@ -293,7 +293,7 @@ def plotCovariances(mu, covs, covsModel, covsWeights, covsNoB, covsModelNoB, cov
                                                               axCov01.flatten(), axCov10.flatten())):
 
             muAmp, cov, model, weight = mu[amp], covs[amp], covsModel[amp], covsWeights[amp]
-            if not np.isnan(np.array(cov)).all():  # If all the entries ara np.nan, this is a bad amp.
+            if not np.isnan(np.array(cov)).all():  # If all the entries are np.nan, this is a bad amp.
                 aCoeffs, bCoeffs = np.array(aDict[amp]), np.array(bDict[amp])
                 gain, noise = gainDict[amp], noiseDict[amp]
                 (meanVecOriginal, varVecOriginal, varVecModelOriginal,

python/lsst/cp/pipe/ptc.py

@@ -319,26 +319,27 @@ def runDataRef(self, dataRefList):
             expId2 = exp2.getInfo().getVisitInfo().getExposureId()
             tupleRows = []
             nAmpsNan = 0
+            tags = ['mu', 'i', 'j', 'var', 'cov', 'npix', 'ext', 'expTime', 'ampName']
             for ampNumber, amp in enumerate(detector):
                 ampName = amp.getName()
                 # covAstier: (i, j, var (cov[0,0]), cov, npix)
                 doRealSpace = self.config.covAstierRealSpace
                 muDiff, varDiff, covAstier = self.measureMeanVarCov(exp1, exp2, region=amp.getBBox(),
                                                                     covAstierRealSpace=doRealSpace)
-                datasetPtc.rawExpTimes[ampName].append(expTime)
-                datasetPtc.rawMeans[ampName].append(muDiff)
-                datasetPtc.rawVars[ampName].append(varDiff)
 
                 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},"
                            f" {expId2} of detector {detNum}.")
                     self.log.warn(msg)
                     nAmpsNan += 1
                     continue
-                tags = ['mu', 'i', 'j', 'var', 'cov', 'npix', 'ext', 'expTime', 'ampName']
                 if (muDiff <= minMeanSignalDict[ampName]) or (muDiff >= maxMeanSignalDict[ampName]):
                     continue
 
+                datasetPtc.rawExpTimes[ampName].append(expTime)
+                datasetPtc.rawMeans[ampName].append(muDiff)
+                datasetPtc.rawVars[ampName].append(varDiff)
+
                 tupleRows += [(muDiff, ) + covRow + (ampNumber, expTime, ampName) for covRow in covAstier]
             if nAmpsNan == len(ampNames):
                 msg = f"NaN mean in all amps of exposure pair {expId1}, {expId2} of detector {detNum}."
@@ -347,6 +348,12 @@ def runDataRef(self, dataRefList):
             allTags += tags
             tupleRecords += tupleRows
         covariancesWithTags = np.core.records.fromrecords(tupleRecords, names=allTags)
+        # Sort raw vectors by rawMeans index
+        for ampName in datasetPtc.ampNames:
+            index = np.argsort(datasetPtc.rawMeans[ampName])
+            datasetPtc.rawExpTimes[ampName] = np.array(datasetPtc.rawExpTimes[ampName])[index]
+            datasetPtc.rawMeans[ampName] = np.array(datasetPtc.rawMeans[ampName])[index]
+            datasetPtc.rawVars[ampName] = np.array(datasetPtc.rawVars[ampName])[index]
 
         if self.config.ptcFitType in ["FULLCOVARIANCE", ]:
             # Calculate covariances and fit them, including the PTC, to Astier+19 full model (Eq. 20)
@@ -442,7 +449,8 @@ def makePairs(self, dataRefList):
                 self.log.warn(f"Only one exposure found at expTime {key}. Dropping exposure "
                               f"{flatPairs[key][0].getInfo().getVisitInfo().getExposureId()}.")
                 flatPairs.pop(key)
-        return flatPairs
+        sortedFlatPairs = {k: flatPairs[k] for k in sorted(flatPairs)}
+        return sortedFlatPairs
 
     def fitCovariancesAstier(self, dataset, covariancesWithTagsArray):
         """Fit measured flat covariances to full model in Astier+19.
@@ -510,7 +518,8 @@ def getOutputPtcDataCovAstier(self, dataset, covFits, covFitsNoB):
             dataset.ptcFitPars[amp] = np.nan
             dataset.ptcFitParsError[amp] = np.nan
             dataset.ptcFitChiSq[amp] = np.nan
-            if amp in covFits:
+            if (amp in covFits and (covFits[amp].covParams is not None) and
+                    (covFitsNoB[amp].covParams is not None)):
                 fit = covFits[amp]
                 fitNoB = covFitsNoB[amp]
                 # Save full covariances, covariances models, and their weights
@@ -527,6 +536,8 @@ def getOutputPtcDataCovAstier(self, dataset, covFits, covFitsNoB):
                 (meanVecFinal, varVecFinal, varVecModel,
                     wc, varMask) = fit.getFitData(0, 0, divideByMu=False, returnMasked=True)
                 gain = fit.getGain()
+                # adjust mask to original size of rawExpTimes
+                # so far, only the min/max signal cut is in dataset.expIdMask
                 dataset.expIdMask[amp] = varMask
                 dataset.gain[amp] = gain
                 dataset.gainErr[amp] = fit.getGainErr()
@@ -800,16 +811,20 @@ def _initialParsForPolynomial(order):
         return pars
 
     @staticmethod
-    def _boundsForPolynomial(initialPars):
-        lowers = [np.NINF for p in initialPars]
-        uppers = [np.inf for p in initialPars]
+    def _boundsForPolynomial(initialPars, lowers=[], uppers=[]):
+        if not len(lowers):
+            lowers = [np.NINF for p in initialPars]
+        if not len(uppers):
+            uppers = [np.inf for p in initialPars]
         lowers[1] = 0  # no negative gains
         return (lowers, uppers)
 
     @staticmethod
-    def _boundsForAstier(initialPars):
-        lowers = [np.NINF for p in initialPars]
-        uppers = [np.inf for p in initialPars]
+    def _boundsForAstier(initialPars, lowers=[], uppers=[]):
+        if not len(lowers):
+            lowers = [np.NINF for p in initialPars]
+        if not len(uppers):
+            uppers = [np.inf for p in initialPars]
         return (lowers, uppers)
 
     @staticmethod
@@ -971,7 +986,9 @@ def errFunc(p, x, y):
             if ptcFitType == 'EXPAPPROXIMATION':
                 ptcFunc = funcAstier
                 parsIniPtc = [-1e-9, 1.0, 10.]  # a00, gain, noise
-                bounds = self._boundsForAstier(parsIniPtc)
+                # lowers and uppers obtained from studies by C. Lage (UC Davis, 11/2020).
+                bounds = self._boundsForAstier(parsIniPtc, lowers=[-1e-4, 0.5, -100],
+                                               uppers=[1e-4, 2.5, 100])
             if ptcFitType == 'POLYNOMIAL':
                 ptcFunc = funcPolynomial
                 parsIniPtc = self._initialParsForPolynomial(self.config.polynomialFitDegree + 1)

python/lsst/cp/pipe/utils.py

@@ -282,6 +282,51 @@ def makeDataRefList(self, namespace):
             self.refList += refList
 
 
+def irlsFit(initialParams, dataX, dataY, function, weightsY=None):
+    """Iteratively reweighted least squares fit.
+
+    This uses the `lsst.cp.pipe.utils.fitLeastSq`, but applies
+    weights based on the Cauchy distribution to the fitter.  See
+    e.g. Holland and Welsch, 1977, doi:10.1080/03610927708827533
+
+    Parameters
+    ----------
+    initialParams : `list` [`float`]
+        Starting parameters.
+    dataX : `numpy.array` [`float`]
+        Abscissa data.
+    dataY : `numpy.array` [`float`]
+        Ordinate data.
+    function : callable
+        Function to fit.
+    weightsY : `numpy.array` [`float`]
+        Weights to apply to the data.
+
+    Returns
+    -------
+    polyFit : `list` [`float`]
+        Final best fit parameters.
+    polyFitErr : `list` [`float`]
+        Final errors on fit parameters.
+    chiSq : `float`
+        Reduced chi squared.
+    weightsY : `list` [`float`]
+        Final weights used for each point.
+
+    """
+    if not weightsY:
+        weightsY = np.ones_like(dataX)
+
+    polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
+    for iteration in range(10):
+        # Use Cauchy weights
+        resid = np.abs(dataY - function(polyFit, dataX)) / np.sqrt(dataY)
+        weightsY = 1.0 / (1.0 + np.sqrt(resid / 2.385))
+        polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
+
+    return polyFit, polyFitErr, chiSq, weightsY
+
+
 def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None):
     """Do a fit and estimate the parameter errors using using scipy.optimize.leastq.
 
@@ -428,13 +473,13 @@ def funcPolynomial(pars, x):
         Polynomial coefficients. Its length determines the polynomial order.
 
     x : `numpy.array`
-        Signal mu (ADU).
+        Abscisa array.
 
     Returns
     -------
-    C_00 (variance) in ADU^2.
+    Ordinate array after evaluating polynomial of order len(pars)-1 at `x`.
     """
-    return poly.polyval(x, [*pars])  # C_00
+    return poly.polyval(x, [*pars])
 
 
 def funcAstier(pars, x):