Update isrMock for multi-component fringes.

Adding multi-component fringes to the mocks allows full tests of how
FringeTask handles these fits.  Return fringe solution and RMS values
for validation.

python/lsst/ip/isr/fringe.py

@@ -131,18 +131,17 @@ def run(self, exposure, fringes, seed=None):
             if self.config.pedestal:
                 self.removePedestal(fringe)
 
-        # Placeholder implementation for multiple fringe frames
-        # This needs to be revisited in DM-4441
         positions = self.generatePositions(fringes[0], rng)
         fluxes = numpy.ndarray([self.config.num, len(fringes)])
         for i, f in enumerate(fringes):
             fluxes[:, i] = self.measureExposure(f, positions, title="Fringe frame")
 
         expFringes = self.measureExposure(exposure, positions, title="Science")
-        solution = self.solve(expFringes, fluxes)
+        solution, rms = self.solve(expFringes, fluxes)
         self.subtract(exposure, fringes, solution)
         if display:
             afwDisplay.Display(frame=getFrame()).mtv(exposure, title="Fringe subtracted")
+        return solution, rms
 
     @timeMethod
     def runDataRef(self, exposure, dataRef, assembler=None):
@@ -160,7 +159,7 @@ def runDataRef(self, exposure, dataRef, assembler=None):
         if not self.checkFilter(exposure):
             return
         fringeStruct = self.readFringes(dataRef, assembler=assembler)
-        self.run(exposure, **fringeStruct.getDict())
+        return self.run(exposure, **fringeStruct.getDict())
 
     def checkFilter(self, exposure):
         """Check whether we should fringe-subtract the science exposure"""
@@ -248,7 +247,7 @@ def emptyResult(msg=""):
             out = [0]
             if len(fringes) > 1:
                 out = out*len(fringes)
-            return numpy.array(out)
+            return numpy.array(out), numpy.nan
 
         good = numpy.where(numpy.logical_and(numpy.isfinite(science), numpy.any(numpy.isfinite(fringes), 1)))
         science = science[good]
@@ -328,7 +327,7 @@ def emptyResult(msg=""):
         # Final solution without rejection
         solution = self._solve(science, fringes)
         self.log.info("Fringe solution: %s RMS: %f Good: %d/%d", solution, rms, len(science), origNum)
-        return solution
+        return solution, rms
 
     def _solve(self, science, fringes):
         """Solve for the scale factors

python/lsst/ip/isr/isrMock.py

@@ -19,6 +19,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+import copy
 import numpy as np
 import tempfile
 
@@ -136,10 +137,20 @@ class IsrMockConfig(pexConfig.Config):
         default=0.1,
         doc="Fractional flux drop due to flat from center to edge of detector along x-axis.",
     )
-    fringeScale = pexConfig.Field(
+    fringeScale = pexConfig.ListField(
         dtype=float,
-        default=200.0,
-        doc="Peak flux for the fringe ripple in DN.",
+        default=[200.0],
+        doc="Peak fluxes for the components of the fringe ripple in DN.",
+    )
+    fringeX0 = pexConfig.ListField(
+        dtype=float,
+        default=[-100],
+        doc="Center position for the fringe ripples.",
+    )
+    fringeY0 = pexConfig.ListField(
+        dtype=float,
+        default=[-0],
+        doc="Center position for the fringe ripples.",
     )
 
     # Inclusion parameters
@@ -156,7 +167,7 @@ class IsrMockConfig(pexConfig.Config):
     doAddCrosstalk = pexConfig.Field(
         dtype=bool,
         default=False,
-        doc="Apply simulated crosstalk to output image.  This cannot be corrected by ISR, "
+        doc="Apply simulated crosstalk to output image. This cannot be corrected by ISR, "
         "as detector.hasCrosstalk()==False.",
     )
     doAddOverscan = pexConfig.Field(
@@ -414,7 +425,9 @@ def makeImage(self):
                         self.amplifierAddSource(ampData, sourceFlux, sourceX, sourceY)
 
             if self.config.doAddFringe is True:
-                self.amplifierAddFringe(amp, ampData, self.config.fringeScale)
+                self.amplifierAddFringe(amp, ampData, np.array(self.config.fringeScale),
+                                        x0=np.array(self.config.fringeX0),
+                                        y0=np.array(self.config.fringeY0))
 
             if self.config.doAddFlat is True:
                 if ampData.getArray().sum() == 0.0:
@@ -662,18 +675,22 @@ def amplifierAddCT(self, ampDataSource, ampDataTarget, scale):
         ampDataTarget.array[:] = (ampDataTarget.array[:] +
                                   scale * ampDataSource.array[:])
 
-    # Functional form data values
-    def amplifierAddFringe(self, amp, ampData, scale):
+    # Functional form data values.
+    def amplifierAddFringe(self, amp, ampData, scale, x0=100, y0=0):
         """Add a fringe-like ripple pattern to an amplifier's image data.
 
         Parameters
         ----------
-        amp : `lsst.afw.ampInfo.AmpInfoRecord`
+        amp : `~lsst.afw.ampInfo.AmpInfoRecord`
             Amplifier to operate on. Needed for amp<->exp coordinate transforms.
         ampData : `lsst.afw.image.ImageF`
             Amplifier image to operate on.
-        scale : `float`
+        scale : `numpy.array` or `float`
             Peak intensity scaling for the ripple.
+        x0 : `numpy.array` or `float`, optional
+            Fringe center
+        y0 : `numpy.array` or `float`, optional
+            Fringe center
 
         Notes
         -----
@@ -686,8 +703,10 @@ def amplifierAddFringe(self, amp, ampData, scale):
         for x in range(0, ampData.getDimensions().getX()):
             for y in range(0, ampData.getDimensions().getY()):
                 (u, v) = self.localCoordToExpCoord(amp, x, y)
-                ampData.array[y][x] = (ampData.array[y][x] +
-                                       scale * np.sinc(((u - 100)/150)**2 + (v / 150)**2))
+                ampData.getArray()[y][x] = np.sum((ampData.getArray()[y][x] +
+                                                   scale *
+                                                   np.sinc(((u - x0) / 50)**2 +
+                                                           ((v - y0) / 50)**2)))
 
     def amplifierMultiplyFlat(self, amp, ampData, fracDrop, u0=100.0, v0=100.0):
         """Multiply an amplifier's image data by a flat-like pattern.
@@ -755,12 +774,17 @@ class CalibratedRawMock(RawMock):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
         self.config.isTrimmed = True
+        self.config.doGenerateImage = True
         self.config.doAddOverscan = False
-        self.config.doAddBias = True
+        self.config.doAddSky = True
+        self.config.doAddSource = True
+        self.config.doAddCrosstalk = False
+
+        self.config.doAddBias = False
         self.config.doAddDark = False
         self.config.doAddFlat = False
-        self.config.doAddFringe = False
-        self.config.doAddCrosstalk = True
+        self.config.doAddFringe = True
+
         self.config.biasLevel = 0.0
         self.config.readNoise = 10.0
 
@@ -981,3 +1005,88 @@ def put(self, exposure, filename):
             Base name of the output file.
         """
         exposure.writeFits(filename+".fits")
+
+
+class FringeDataRefMock(object):
+    """Simulated gen2 butler data ref.
+
+    Currently only supports get and put operations, which are most
+    likely to be called for data in ISR processing.
+
+    """
+    dataId = "isrMock Fake Data"
+    darkval = 2.  # e-/sec
+    oscan = 250.  # DN
+    gradient = .10
+    exptime = 15  # seconds
+    darkexptime = 40.  # seconds
+
+    def __init__(self, **kwargs):
+        if 'config' in kwargs.keys():
+            self.config = kwargs['config']
+        else:
+            self.config = IsrMockConfig()
+            self.config.isTrimmed = True
+            self.config.doAddFringe = True
+            self.config.readNoise = 10.0
+
+    def get(self, dataType, **kwargs):
+        """Return an appropriate data product.
+
+        Parameters
+        ----------
+        dataType : `str`
+            Type of data product to return.
+
+        Returns
+        -------
+        mock : IsrMock.run() result
+            The output product.
+        """
+        if "_filename" in dataType:
+            return tempfile.mktemp(), "mock"
+        elif 'transmission_' in dataType:
+            return TransmissionMock(config=self.config).run()
+        elif dataType == 'ccdExposureId':
+            return 20090913
+        elif dataType == 'camera':
+            return IsrMock(config=self.config).getCamera()
+        elif dataType == 'raw':
+            return CalibratedRawMock(config=self.config).run()
+        elif dataType == 'bias':
+            return BiasMock(config=self.config).run()
+        elif dataType == 'dark':
+            return DarkMock(config=self.config).run()
+        elif dataType == 'flat':
+            return FlatMock(config=self.config).run()
+        elif dataType == 'fringe':
+            fringes = []
+            configCopy = copy.deepcopy(self.config)
+            for scale, x, y in zip(self.config.fringeScale, self.config.fringeX0, self.config.fringeY0):
+                configCopy.fringeScale = [1.0]
+                configCopy.fringeX0 = [x]
+                configCopy.fringeY0 = [y]
+                fringes.append(FringeMock(config=configCopy).run())
+            return fringes
+        elif dataType == 'defects':
+            return DefectMock(config=self.config).run()
+        elif dataType == 'bfKernel':
+            return BfKernelMock(config=self.config).run()
+        elif dataType == 'linearizer':
+            return None
+        elif dataType == 'crosstalkSources':
+            return None
+        else:
+            return None
+
+    def put(self, exposure, filename):
+        """Write an exposure to a FITS file.
+
+        Parameters
+        ----------
+        exposure : `lsst.afw.image.Exposure`
+            Image data to write out.
+        filename : `str`
+            Base name of the output file.
+        """
+        exposure.writeFits(filename+".fits")

tests/test_fringes.py

@@ -265,6 +265,34 @@ def test_readFringes(self):
         result = task.readFringes(dataRef, assembler=None)
         self.assertIsInstance(result, pipeBase.Struct)
 
+    def test_multiFringes(self):
+        """Test that multi-fringe results are handled correctly by the task.
+        """
+        self.config.filters.append("_unknown_")
+        self.config.large = 16
+        task = FringeTask(name="multiFringeMock", config=self.config)
+
+        config = isrMock.IsrMockConfig()
+        config.fringeScale = [750.0, 240.0, 220.0]
+        config.fringeX0 = [100.0, 150.0, 200.0]
+        config.fringeY0 = [0.0, 200.0, 0.0]
+        dataRef = isrMock.FringeDataRefMock(config=config)
+
+        exp = dataRef.get("raw")
+        medianBefore = np.nanmedian(exp.getImage().getArray())
+        fringes = task.readFringes(dataRef, assembler=None)
+
+        solution, rms = task.run(exp, **fringes.getDict())
+        medianAfter = np.nanmedian(exp.getImage().getArray())
+        stdAfter = np.nanstd(exp.getImage().getArray())
+
+        self.assertLess(medianAfter, medianBefore)
+        self.assertFloatsAlmostEqual(medianAfter, 3002.233, atol=1e-4)
+        self.assertFloatsAlmostEqual(stdAfter, 3549.9375, atol=1e-4)
+
+        deviation = np.abs(solution - config.fringeScale)
+        self.assertTrue(np.all(deviation / rms < 1.0))
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass

tests/test_isrTask.py

@@ -147,7 +147,7 @@ def test_updateVariance(self):
         statAfter = computeImageMedianAndStd(self.inputExp.variance[self.amp.getBBox()])
         self.assertGreater(statAfter[0], statBefore[0])
         self.assertFloatsAlmostEqual(statBefore[0], 0.0, atol=1e-2)
-        self.assertFloatsAlmostEqual(statAfter[0], 8175.6885, atol=1e-2)
+        self.assertFloatsAlmostEqual(statAfter[0], 8170.0195, atol=1e-2)
 
     def test_darkCorrection(self):
         """Expect the median image value should decrease after this operation.
@@ -158,8 +158,8 @@ def test_darkCorrection(self):
         self.task.darkCorrection(self.inputExp, darkIm)
         statAfter = computeImageMedianAndStd(self.inputExp.image[self.amp.getBBox()])
         self.assertLess(statAfter[0], statBefore[0])
-        self.assertFloatsAlmostEqual(statBefore[0], 8075.6885, atol=1e-2)
-        self.assertFloatsAlmostEqual(statAfter[0], 8051.6206, atol=1e-2)
+        self.assertFloatsAlmostEqual(statBefore[0], 8070.0195, atol=1e-2)
+        self.assertFloatsAlmostEqual(statAfter[0], 8045.7773, atol=1e-2)
 
     def test_darkCorrection_noVisitInfo(self):
         """Expect the median image value should decrease after this operation.
@@ -171,8 +171,8 @@ def test_darkCorrection_noVisitInfo(self):
         self.task.darkCorrection(self.inputExp, darkIm)
         statAfter = computeImageMedianAndStd(self.inputExp.image[self.amp.getBBox()])
         self.assertLess(statAfter[0], statBefore[0])
-        self.assertFloatsAlmostEqual(statBefore[0], 8075.6885, atol=1e-2)
-        self.assertFloatsAlmostEqual(statAfter[0], 8051.6206, atol=1e-2)
+        self.assertFloatsAlmostEqual(statBefore[0], 8070.0195, atol=1e-2)
+        self.assertFloatsAlmostEqual(statAfter[0], 8045.7773, atol=1e-2)
 
     def test_flatCorrection(self):
         """Expect the image median should increase (divide by < 1).
@@ -183,8 +183,8 @@ def test_flatCorrection(self):
         self.task.flatCorrection(self.inputExp, flatIm)
         statAfter = computeImageMedianAndStd(self.inputExp.image[self.amp.getBBox()])
         self.assertGreater(statAfter[1], statBefore[1])
-        self.assertFloatsAlmostEqual(statAfter[1], 147417.36, atol=1e-2)
-        self.assertFloatsAlmostEqual(statBefore[1], 148.28436, atol=1e-2)
+        self.assertFloatsAlmostEqual(statAfter[1], 147407.02, atol=1e-2)
+        self.assertFloatsAlmostEqual(statBefore[1], 147.55304, atol=1e-2)
 
     def test_saturationDetection(self):
         """Expect the saturation level detection/masking to scale with
@@ -200,7 +200,7 @@ def test_saturationDetection(self):
 
         self.assertLessEqual(countBefore, countAfter)
         self.assertEqual(countBefore, 43)
-        self.assertEqual(countAfter, 141)
+        self.assertEqual(countAfter, 136)
 
     def test_measureBackground(self):
         """Expect the background measurement runs successfully and to save
@@ -222,7 +222,7 @@ def test_flatContext(self):
         mi = self.inputExp.getMaskedImage().clone()
         with self.task.flatContext(self.inputExp, flatExp, darkExp):
             contextStat = computeImageMedianAndStd(self.inputExp.getMaskedImage().getImage())
-            self.assertFloatsAlmostEqual(contextStat[0], 37269.914, atol=1e-2)
+            self.assertFloatsAlmostEqual(contextStat[0], 37165.594, atol=1e-2)
 
         self.assertMaskedImagesAlmostEqual(mi, self.inputExp.getMaskedImage())