DM-26453: Add sqrt(var) as weight to EXPAPPROXIMATION and POLYNOMIAL fit residual in ptc.py

python/lsst/cp/pipe/plotPtc.py

@@ -787,6 +787,7 @@ def _plotStandardPtc(self, dataset, ptcFitType, pdfPages):
             meanVecOutliers = meanVecOriginal[np.invert(mask)]
             varVecOutliers = varVecOriginal[np.invert(mask)]
             pars, parsErr = dataset.ptcFitPars[amp], dataset.ptcFitParsError[amp]
+            ptcRedChi2 = dataset.ptcFitReducedChiSquared[amp]
             if ptcFitType == 'EXPAPPROXIMATION':
                 if len(meanVecFinal):
                     ptcA00, ptcA00error = pars[0], parsErr[0]
@@ -795,15 +796,17 @@ def _plotStandardPtc(self, dataset, ptcFitType, pdfPages):
                     ptcNoiseError = 0.5*(parsErr[2]/np.fabs(pars[2]))*np.sqrt(np.fabs(pars[2]))
                     stringLegend = (f"a00: {ptcA00:.2e}+/-{ptcA00error:.2e} 1/e"
                                     f"\n Gain: {ptcGain:.4}+/-{ptcGainError:.2e} e/DN"
-                                    f"\n Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n")
+                                    f"\n Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n"
+                                    r"$\chi^2_{\rm{red}}$: " + f"{ptcRedChi2:.4}")
 
             if ptcFitType == 'POLYNOMIAL':
                 if len(meanVecFinal):
                     ptcGain, ptcGainError = 1./pars[1], np.fabs(1./pars[1])*(parsErr[1]/pars[1])
                     ptcNoise = np.sqrt((pars[0]))*ptcGain
                     ptcNoiseError = (0.5*(parsErr[0]/np.fabs(pars[0]))*(np.sqrt(np.fabs(pars[0]))))*ptcGain
                     stringLegend = (f"Gain: {ptcGain:.4}+/-{ptcGainError:.2e} e/DN \n"
-                                    f"Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n")
+                                    f"Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n"
+                                    r"$\chi^2_{\rm{red}}$: " + f"{ptcRedChi2:.4}")
 
                 a.set_xlabel(r'Mean signal ($\mu$, DN)', fontsize=labelFontSize)
                 a.set_ylabel(r'Variance (DN$^2$)', fontsize=labelFontSize)
@@ -906,8 +909,10 @@ def _plotLinearizer(self, dataset, linearizer, pdfPages):
                 k0, k0Error = pars[0], parsErr[0]
                 k1, k1Error = pars[1], parsErr[1]
                 k2, k2Error = pars[2], parsErr[2]
+                linRedChi2 = linearizer.linearityFitReducedChiSquared[amp]
                 stringLegend = (f"k0: {k0:.4}+/-{k0Error:.2e} DN\n k1: {k1:.4}+/-{k1Error:.2e} DN/t"
-                                f"\n k2: {k2:.2e}+/-{k2Error:.2e} DN/t^2 \n")
+                                f"\n k2: {k2:.2e}+/-{k2Error:.2e} DN/t^2 \n"
+                                r"$\chi^2_{\rm{red}}$: " + f"{linRedChi2:.4}")
                 a.scatter(timeVecFinal, meanVecFinal)
                 a.plot(timeVecFinal, funcPolynomial(pars, timeVecFinal), color='red')
                 a.text(0.03, 0.75, stringLegend, transform=a.transAxes, fontsize=legendFontSize)

python/lsst/cp/pipe/ptc.py

@@ -625,8 +625,9 @@ def measureMeanVarCov(self, exposure1, exposure2, region=None, covAstierRealSpac
             im1Area = exposure1.maskedImage
             im2Area = exposure2.maskedImage
 
-        im1Area = afwMath.binImage(im1Area, self.config.binSize)
-        im2Area = afwMath.binImage(im2Area, self.config.binSize)
+        if self.config.binSize > 1:
+            im1Area = afwMath.binImage(im1Area, self.config.binSize)
+            im2Area = afwMath.binImage(im2Area, self.config.binSize)
 
         im1MaskVal = exposure1.getMask().getPlaneBitMask(self.config.maskNameList)
         im1StatsCtrl = afwMath.StatisticsControl(self.config.nSigmaClipPtc,
@@ -940,15 +941,19 @@ def calculateLinearityResidualAndLinearizers(self, exposureTimeVector, meanSigna
         # Lookup table linearizer
         parsIniNonLinearity = self._initialParsForPolynomial(self.config.polynomialFitDegreeNonLinearity + 1)
         if self.config.doFitBootstrap:
-            parsFit, parsFitErr, reducedChiSquaredNonLinearityFit = fitBootstrap(parsIniNonLinearity,
-                                                                                 exposureTimeVector,
-                                                                                 meanSignalVector,
-                                                                                 funcPolynomial)
+            (parsFit, parsFitErr,
+             reducedChiSquaredNonLinearityFit) = fitBootstrap(parsIniNonLinearity,
+                                                              exposureTimeVector,
+                                                              meanSignalVector,
+                                                              funcPolynomial,
+                                                              weightsY=1./np.sqrt(meanSignalVector))
         else:
-            parsFit, parsFitErr, reducedChiSquaredNonLinearityFit = fitLeastSq(parsIniNonLinearity,
-                                                                               exposureTimeVector,
-                                                                               meanSignalVector,
-                                                                               funcPolynomial)
+            (parsFit, parsFitErr,
+             reducedChiSquaredNonLinearityFit) = fitLeastSq(parsIniNonLinearity,
+                                                            exposureTimeVector,
+                                                            meanSignalVector,
+                                                            funcPolynomial,
+                                                            weightsY=1./np.sqrt(meanSignalVector))
 
         # LinearizeLookupTable:
         # Use linear part to get time at wich signal is maxAduForLookupTableLinearizer DN
@@ -1243,10 +1248,12 @@ def errFunc(p, x, y):
             # Fit the PTC
             if self.config.doFitBootstrap:
                 parsFit, parsFitErr, reducedChiSqPtc = fitBootstrap(parsIniPtc, meanVecFinal,
-                                                                    varVecFinal, ptcFunc)
+                                                                    varVecFinal, ptcFunc,
+                                                                    weightsY=1./np.sqrt(varVecFinal))
             else:
                 parsFit, parsFitErr, reducedChiSqPtc = fitLeastSq(parsIniPtc, meanVecFinal,
-                                                                  varVecFinal, ptcFunc)
+                                                                  varVecFinal, ptcFunc,
+                                                                  weightsY=1./np.sqrt(varVecFinal))
             dataset.ptcFitPars[ampName] = parsFit
             dataset.ptcFitParsError[ampName] = parsFitErr
             dataset.ptcFitReducedChiSquared[ampName] = reducedChiSqPtc

python/lsst/cp/pipe/utils.py

@@ -282,7 +282,7 @@ def makeDataRefList(self, namespace):
             self.refList += refList
 
 
-def fitLeastSq(initialParams, dataX, dataY, function):
+def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None):
     """Do a fit and estimate the parameter errors using using scipy.optimize.leastq.
 
     optimize.leastsq returns the fractional covariance matrix. To estimate the
@@ -304,6 +304,9 @@ def fitLeastSq(initialParams, dataX, dataY, function):
     function : callable object (function)
         Function to fit the data with.
 
+    weightsY : `numpy.array` of `float`
+        Weights of the data in the ordinate axis.
+
     Return
     ------
     pFitSingleLeastSquares : `list` of `float`
@@ -313,17 +316,23 @@ def fitLeastSq(initialParams, dataX, dataY, function):
         List with errors for fitted parameters.
 
     reducedChiSqSingleLeastSquares : `float`
-        Unweighted reduced chi squared
+        Reduced chi squared, unweighted if weightsY is not provided.
     """
+    if weightsY is None:
+        weightsY = np.ones(len(dataX))
 
-    def errFunc(p, x, y):
-        return function(p, x) - y
+    def errFunc(p, x, y, weightsY=None):
+        if weightsY is None:
+            weightsY = np.ones(len(x))
+        return (function(p, x) - y)*weightsY
 
     pFit, pCov, infoDict, errMessage, success = leastsq(errFunc, initialParams,
-                                                        args=(dataX, dataY), full_output=1, epsfcn=0.0001)
+                                                        args=(dataX, dataY, weightsY), full_output=1,
+                                                        epsfcn=0.0001)
 
     if (len(dataY) > len(initialParams)) and pCov is not None:
-        reducedChiSq = (errFunc(pFit, dataX, dataY)**2).sum()/(len(dataY)-len(initialParams))
+        reducedChiSq = calculateWeightedReducedChi2(dataY, function(pFit, dataX), weightsY, len(dataY),
+                                                    len(initialParams))
         pCov *= reducedChiSq
     else:
         pCov = np.zeros((len(initialParams), len(initialParams)))
@@ -340,7 +349,7 @@ def errFunc(p, x, y):
     return pFitSingleLeastSquares, pErrSingleLeastSquares, reducedChiSq
 
 
-def fitBootstrap(initialParams, dataX, dataY, function, confidenceSigma=1.):
+def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.):
     """Do a fit using least squares and bootstrap to estimate parameter errors.
 
     The bootstrap error bars are calculated by fitting 100 random data sets.
@@ -360,7 +369,10 @@ def fitBootstrap(initialParams, dataX, dataY, function, confidenceSigma=1.):
     function : callable object (function)
         Function to fit the data with.
 
-    confidenceSigma : `float`
+    weightsY : `numpy.array` of `float`, optional.
+        Weights of the data in the ordinate axis.
+
+    confidenceSigma : `float`, optional.
         Number of sigmas that determine confidence interval for the bootstrap errors.
 
     Return
@@ -372,37 +384,39 @@ def fitBootstrap(initialParams, dataX, dataY, function, confidenceSigma=1.):
         List with errors for fitted parameters.
 
     reducedChiSqBootstrap : `float`
-        Reduced chi squared.
+        Reduced chi squared, unweighted if weightsY is not provided.
     """
+    if weightsY is None:
+        weightsY = np.ones(len(dataX))
 
-    def errFunc(p, x, y):
-        return function(p, x) - y
+    def errFunc(p, x, y, weightsY):
+        if weightsY is None:
+            weightsY = np.ones(len(x))
+        return (function(p, x) - y)*weightsY
 
     # Fit first time
-    pFit, _ = leastsq(errFunc, initialParams, args=(dataX, dataY), full_output=0)
+    pFit, _ = leastsq(errFunc, initialParams, args=(dataX, dataY, weightsY), full_output=0)
 
     # Get the stdev of the residuals
-    residuals = errFunc(pFit, dataX, dataY)
-    sigmaErrTotal = np.std(residuals)
-
+    residuals = errFunc(pFit, dataX, dataY, weightsY)
     # 100 random data sets are generated and fitted
     pars = []
     for i in range(100):
-        randomDelta = np.random.normal(0., sigmaErrTotal, len(dataY))
+        randomDelta = np.random.normal(0., np.fabs(residuals), len(dataY))
         randomDataY = dataY + randomDelta
         randomFit, _ = leastsq(errFunc, initialParams,
-                               args=(dataX, randomDataY), full_output=0)
+                               args=(dataX, randomDataY, weightsY), full_output=0)
         pars.append(randomFit)
     pars = np.array(pars)
     meanPfit = np.mean(pars, 0)
 
     # confidence interval for parameter estimates
-    nSigma = confidenceSigma
-    errPfit = nSigma*np.std(pars, 0)
+    errPfit = confidenceSigma*np.std(pars, 0)
     pFitBootstrap = meanPfit
     pErrBootstrap = errPfit
 
-    reducedChiSq = (errFunc(pFitBootstrap, dataX, dataY)**2).sum()/(len(dataY)-len(initialParams))
+    reducedChiSq = calculateWeightedReducedChi2(dataY, function(pFitBootstrap, dataX), weightsY, len(dataY),
+                                                len(initialParams))
     return pFitBootstrap, pErrBootstrap, reducedChiSq
 
 

tests/test_ptc.py

@@ -136,9 +136,16 @@ def test_covAstier(self):
                 v2 *= (0.75**2)  # convert to electrons
                 self.assertAlmostEqual(v1/v2, 1.0, places=1)
 
-    def ptcFitAndCheckPtc(self, order=None, fitType='', doTableArray=False):
+    def ptcFitAndCheckPtc(self, order=None, fitType='', doTableArray=False, doFitBootstrap=False):
         localDataset = copy.copy(self.dataset)
         config = copy.copy(self.defaultConfig)
+        placesTests = 6
+        if doFitBootstrap:
+            config.doFitBootstrap = True
+            # Bootstrap method in cp_pipe/utils.py does multiple fits in the precense of noise.
+            # Allow for more margin of error.
+            placesTests = 3
+
         if fitType == 'POLYNOMIAL':
             if order not in [2, 3]:
                 RuntimeError("Enter a valid polynomial order for this test: 2 or 3")
@@ -187,7 +194,8 @@ def ptcFitAndCheckPtc(self, order=None, fitType='', doTableArray=False):
                 linearizerTableRow = signalIdeal - signalUncorrected
                 self.assertEqual(len(linearizerTableRow), len(lookupTableArray[i, :]))
                 for j in np.arange(len(linearizerTableRow)):
-                    self.assertAlmostEqual(linearizerTableRow[j], lookupTableArray[i, :][j], places=6)
+                    self.assertAlmostEqual(linearizerTableRow[j], lookupTableArray[i, :][j],
+                                           places=placesTests)
 
         # check entries in localDataset, which was modified by the function
         for ampName in self.ampNames:
@@ -215,31 +223,43 @@ def ptcFitAndCheckPtc(self, order=None, fitType='', doTableArray=False):
             inputFracNonLinearityResiduals = 100*(linearPart - finalMuVec)/linearPart
 
             # Nonlinearity fit parameters
-            self.assertAlmostEqual(0.0, linDataset[ampName].meanSignalVsTimePolyFitPars[0])
-            self.assertAlmostEqual(self.flux, linDataset[ampName].meanSignalVsTimePolyFitPars[1])
-            self.assertAlmostEqual(self.k2NonLinearity, linDataset[ampName].meanSignalVsTimePolyFitPars[2])
+            self.assertLess(np.fabs(linDataset[ampName].meanSignalVsTimePolyFitPars[0]), 0.01)
+            self.assertAlmostEqual(self.flux, linDataset[ampName].meanSignalVsTimePolyFitPars[1],
+                                   places=placesTests)
+            self.assertAlmostEqual(self.k2NonLinearity, linDataset[ampName].meanSignalVsTimePolyFitPars[2],
+                                   places=placesTests)
 
             # Non-linearity coefficient for linearizer
             self.assertAlmostEqual(-self.k2NonLinearity/(self.flux**2),
-                                   linDataset[ampName].quadraticPolynomialLinearizerCoefficient)
+                                   linDataset[ampName].quadraticPolynomialLinearizerCoefficient,
+                                   places=placesTests)
 
             linearPartModel = linDataset[ampName].meanSignalVsTimePolyFitPars[1]*finalTimeVec
             outputFracNonLinearityResiduals = 100*(linearPartModel - finalMuVec)/linearPartModel
             # Fractional nonlinearity residuals
             self.assertEqual(len(outputFracNonLinearityResiduals), len(inputFracNonLinearityResiduals))
             for calc, truth in zip(outputFracNonLinearityResiduals, inputFracNonLinearityResiduals):
-                self.assertAlmostEqual(calc, truth)
+                self.assertAlmostEqual(calc, truth, places=placesTests)
 
             # check calls to calculateLinearityResidualAndLinearizers
             datasetLinResAndLinearizers = task.calculateLinearityResidualAndLinearizers(
                 localDataset.rawExpTimes[ampName], localDataset.rawMeans[ampName])
 
             self.assertAlmostEqual(-self.k2NonLinearity/(self.flux**2),
-                                   datasetLinResAndLinearizers.quadraticPolynomialLinearizerCoefficient)
-            self.assertAlmostEqual(0.0, datasetLinResAndLinearizers.meanSignalVsTimePolyFitPars[0])
-            self.assertAlmostEqual(self.flux, datasetLinResAndLinearizers.meanSignalVsTimePolyFitPars[1])
+                                   datasetLinResAndLinearizers.quadraticPolynomialLinearizerCoefficient,
+                                   places=placesTests)
+            self.assertAlmostEqual(0.0, datasetLinResAndLinearizers.meanSignalVsTimePolyFitPars[0],
+                                   places=placesTests)
+            self.assertAlmostEqual(self.flux, datasetLinResAndLinearizers.meanSignalVsTimePolyFitPars[1],
+                                   places=placesTests)
             self.assertAlmostEqual(self.k2NonLinearity,
-                                   datasetLinResAndLinearizers.meanSignalVsTimePolyFitPars[2])
+                                   datasetLinResAndLinearizers.meanSignalVsTimePolyFitPars[2],
+                                   places=placesTests)
+
+    def test_ptcFitBootstrap(self):
+        """Test the bootstrap fit option for the PTC"""
+        for (fitType, order) in [('POLYNOMIAL', 2), ('POLYNOMIAL', 3), ('EXPAPPROXIMATION', None)]:
+            self.ptcFitAndCheckPtc(fitType=fitType, order=order, doTableArray=False, doFitBootstrap=True)
 
     def test_ptcFit(self):
         for createArray in [True, False]: