Fix data handling in BFK code

python/lsst/cp/pipe/makeBrighterFatterKernel.py

@@ -203,8 +203,10 @@ def run(self, inputPtc, dummy, camera, inputDims):
             detectorFluxes = list()
 
         bfk = BrighterFatterKernel(camera=camera, detectorId=detector.getId(), level=self.config.level)
-        bfk.means = inputPtc.finalMeans  # ADU
-        bfk.variances = inputPtc.finalVars  # ADU^2
+        bfk.rawMeans = inputPtc.rawMeans  # ADU
+        bfk.rawVariances = inputPtc.rawVars  # ADU^2
+        bfk.expIdMask = inputPtc.expIdMask
+
         # Use the PTC covariances as the cross-correlations.  These
         # are scaled before the kernel is generated, which performs
         # the conversion.
@@ -219,30 +221,23 @@ def run(self, inputPtc, dummy, camera, inputDims):
         for amp in detector:
             ampName = amp.getName()
             gain = bfk.gain[ampName]
-
-            # Using the inputPtc.expIdMask works if the covariance
-            # array has the same length as the rawMeans/rawVars.  This
-            # isn't the case, as it's the same size as the
-            # finalMeans/finalVars.  However, these arrays (and the
-            # covariance) are padded with NAN values to match the
-            # longest amplifier vector.  We do not want to include
-            # these NAN values, so we construct a mask for all non-NAN
-            # values in finalMeans, and use that to filter finalVars
-            # and the covariances.
-            mask = np.isfinite(bfk.means[ampName])
-            fluxes = np.array(bfk.means[ampName])[mask]
-            variances = np.array(bfk.variances[ampName])[mask]
-            xCorrList = np.array([np.array(xcorr) for xcorr in bfk.rawXcorrs[ampName]])[mask]
-
+            mask = inputPtc.expIdMask[ampName]
             if gain <= 0:
                 # We've received very bad data.
-                self.log.warning("Impossible gain recieved from PTC for %s: %f.  Skipping amplifier.",
+                self.log.warning("Impossible gain recieved from PTC for %s: %f. Skipping bad amplifier.",
                                  ampName, gain)
                 bfk.meanXcorrs[ampName] = np.zeros(bfk.shape)
                 bfk.ampKernels[ampName] = np.zeros(bfk.shape)
+                bfk.rawXcorrs[ampName] = np.zeros((len(mask), inputPtc.covMatrixSide, inputPtc.covMatrixSide))
                 bfk.valid[ampName] = False
                 continue
 
+            # Use inputPtc.expIdMask to get the means, variances,
+            # and covariances that were not masked after PTC.
+            fluxes = np.array(bfk.rawMeans[ampName])[mask]
+            variances = np.array(bfk.rawVariances[ampName])[mask]
+            xCorrList = np.array([np.array(xcorr) for xcorr in bfk.rawXcorrs[ampName]])[mask]
+
             fluxes = np.array([flux*gain for flux in fluxes])  # Now in e^-
             variances = np.array([variance*gain*gain for variance in variances])  # Now in e^2-
 

tests/test_brighterFatterKernel.py

@@ -50,12 +50,13 @@ def setUp(self):
         self.defaultConfig = cpPipe.BrighterFatterKernelSolveConfig()
         self.ptc = ipIsr.PhotonTransferCurveDataset(ampNames=['amp 1'], ptcFitType='FULLCOVARIANCE',
                                                     covMatrixSide=3)
-        self.ptc.expIdMask['amp 1'] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
-        self.ptc.finalMeans['amp 1'] = [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000]
-        self.ptc.rawMeans['amp 1'] = self.ptc.finalMeans['amp 1']
-        self.ptc.finalVars['amp 1'] = 0.99 * np.array(self.ptc.finalMeans['amp 1'], dtype=float)
+        self.ptc.expIdMask['amp 1'] = np.array([False, True, True, True, True, True, True, True, True, True])
+        self.ptc.rawMeans['amp 1'] = np.array([1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000])
+        self.ptc.finalMeans['amp 1'] = self.ptc.rawMeans['amp 1'][self.ptc.expIdMask['amp 1']]
+        self.ptc.rawVars['amp 1'] = 0.99 * np.array(self.ptc.rawMeans['amp 1'], dtype=float)
+        self.ptc.finalVars['amp 1'] = self.ptc.rawVars['amp 1'][self.ptc.expIdMask['amp 1']]
         self.ptc.covariances['amp 1'] = []
-        for mean, variance in zip(self.ptc.finalMeans['amp 1'], self.ptc.finalVars['amp 1']):
+        for mean, variance in zip(self.ptc.rawMeans['amp 1'], self.ptc.rawVars['amp 1']):
             residual = mean - variance
             covariance = [[variance, 0.5 * residual, 0.1 * residual],
                           [0.2 * residual, 0.1 * residual, 0.05 * residual],
@@ -119,12 +120,13 @@ def test_quadratic(self):
 
         ptc = ipIsr.PhotonTransferCurveDataset(ampNames=['amp 1'], ptcFitType='FULLCOVARIANCE',
                                                covMatrixSide=3)
-        ptc.expIdMask['amp 1'] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
-        ptc.finalMeans['amp 1'] = [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000]
-        ptc.rawMeans['amp 1'] = ptc.finalMeans['amp 1']
-        ptc.finalVars['amp 1'] = 9e-5 * np.square(np.array(ptc.finalMeans['amp 1'], dtype=float))
+        ptc.expIdMask['amp 1'] = np.array([False, True, True, True, True, True, True, True, True, True])
+        ptc.rawMeans['amp 1'] = np.array([1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000])
+        ptc.finalMeans['amp 1'] = ptc.rawMeans['amp 1'][ptc.expIdMask['amp 1']]
+        ptc.rawVars['amp 1'] = 9e-5 * np.square(np.array(ptc.rawMeans['amp 1'], dtype=float))
+        ptc.finalVars['amp 1'] = ptc.rawVars['amp 1'][ptc.expIdMask['amp 1']]
         ptc.covariances['amp 1'] = []
-        for mean, variance in zip(ptc.finalMeans['amp 1'], ptc.finalVars['amp 1']):
+        for mean, variance in zip(ptc.rawMeans['amp 1'], ptc.rawVars['amp 1']):
             residual = variance
             covariance = [[variance, 0.5 * residual, 0.1 * residual],
                           [0.2 * residual, 0.1 * residual, 0.05 * residual],