diff --git a/python/lsst/ip/isr/ampOffset.py b/python/lsst/ip/isr/ampOffset.py index b91f76a5..27c0c14c 100644 --- a/python/lsst/ip/isr/ampOffset.py +++ b/python/lsst/ip/isr/ampOffset.py @@ -120,10 +120,18 @@ def setDefaults(self): dtype=bool, default=True, ) + measureOnly = Field( + doc="Measure amp offsets without applying them to the input exposure?", + dtype=bool, + default=False, + ) + # TODO: Remove this deprecated config after v29.0 on DM-XXXXX. doApplyAmpOffset = Field( doc="Apply amp offset corrections to the input exposure?", dtype=bool, - default=False, + default=True, + deprecated="Please use the `measureOnly` argument to measure without application. This config will " + "be deprecated after v29.0.", ) @@ -214,7 +222,7 @@ def run(self, exposure): # Safety check: do any pixels remain for amp offset estimation? if (exp.mask.array & bitMask).all(): - log_fn = self.log.warning if self.config.doApplyAmpOffset else self.log.debug + log_fn = self.log.debug if self.config.measureOnly else self.log.warning log_fn( "All pixels masked: cannot calculate any amp offset corrections. All pedestals are being set " "to zero." @@ -251,23 +259,21 @@ def run(self, exposure): float(pedestal), f"Pedestal level calculated for amp {ampName}", ) - if self.config.doApplyAmpOffset: + if not self.config.measureOnly: ampIm = exposure.image[amp.getBBox()].array ampIm -= pedestal # Add the amp pedestal to the "Task" metadata as well. # Needed for Sasquatch/Chronograf! self.metadata["AMPOFFSET_PEDESTALS"][ampName] = float(pedestal) - if self.config.doApplyAmpOffset: + if not self.config.measureOnly: status = "subtracted from exposure" - metadata.set( - "LSST ISR AMPOFFSET PEDESTAL SUBTRACTED", True, "Amp pedestals have been subtracted" - ) + metadata.set("LSST ISR AMPOFFSET PEDESTAL SUBTRACTED", True, "Amp pedestals have been subtracted") else: status = "not subtracted from exposure" metadata.set( "LSST ISR AMPOFFSET PEDESTAL SUBTRACTED", False, "Amp pedestals have not been subtracted" ) - self.log.info(f"amp pedestal values ({status}): {', '.join([f'{x:.4f}' for x in pedestals])}") + self.log.info(f"Amp pedestal values ({status}): {', '.join([f'{x:.4f}' for x in pedestals])}") return Struct(pedestals=pedestals) @@ -410,6 +416,11 @@ def getAmpOffsets(self, im, amps, associations, sides): ampsEdges = self.getAmpEdges(im, amps, sides) interfaceOffsetLookup = {} + interfaceIds = [] + interfaceOffsetOriginals = [] + ampEdgeGoodFracs = [] + minFracFails = [] + maxOffsetFails = [] for ampId, ampAssociations in enumerate(associations): ampNeighbors = np.ravel(np.where(ampAssociations < 0)) for ampNeighbor in ampNeighbors: @@ -420,11 +431,59 @@ def getAmpOffsets(self, im, amps, associations, sides): edgeA = ampsEdges[ampId][ampSide] edgeB = ampsEdges[ampNeighbor][(ampSide + 2) % 4] if ampId < ampNeighbor: - interfaceOffset = self.getInterfaceOffset(ampId, ampNeighbor, edgeA, edgeB) - interfaceOffsetLookup[f"{ampId}{ampNeighbor}"] = interfaceOffset + ( + interfaceId, + interfaceOffset, + interfaceOffsetOriginal, + ampEdgeGoodFrac, + minFracFail, + maxOffsetFail, + ) = self.getInterfaceOffset(ampId, ampNeighbor, edgeA, edgeB) + interfaceIds.append(interfaceId) + interfaceOffsetOriginals.append(interfaceOffsetOriginal) + ampEdgeGoodFracs.append(ampEdgeGoodFrac) + minFracFails.append(minFracFail) + maxOffsetFails.append(maxOffsetFail) + interfaceOffsetLookup[f"{ampId:02d}:{ampNeighbor:02d}"] = interfaceOffset else: - interfaceOffset = -interfaceOffsetLookup[f"{ampNeighbor}{ampId}"] + interfaceOffset = -interfaceOffsetLookup[f"{ampNeighbor:02d}:{ampId:02d}"] ampsOffsets[ampId] += interfaceWeight * interfaceOffset + quartile_summary = np.percentile(interfaceOffsetOriginals, [0, 25, 50, 75, 100]) + self.log.info( + "Amp interface offset quartile summary (min, Q1, median, Q3, max): %.4f, %.4f, %.4f, %.4f, %.4f", + *quartile_summary, + ) + log_fn = self.log.debug if self.config.measureOnly else self.log.warning + if any(minFracFails): + log_fn( + "The fraction of unmasked edge pixels for the following amp interfaces is below the " + "configured threshold (%s): %s", + self.config.ampEdgeMinFrac, + ", ".join( + [ + f"{interfaceId} ({ampEdgeGoodFrac:0.2f})" + for interfaceId, ampEdgeGoodFrac, minFracFail in zip( + interfaceIds, ampEdgeGoodFracs, minFracFails + ) + if minFracFail + ] + ), + ) + if any(maxOffsetFails): + log_fn( + "The absolute amp edge offset value for the following amp interfaces exceeds the configured " + "limit (%s): %s", + self.config.ampEdgeMaxOffset, + ", ".join( + [ + f"{interfaceId} ({np.abs(interfaceOffset):0.2f})" + for interfaceId, interfaceOffset, maxOffsetFail in zip( + interfaceIds, interfaceOffsetOriginals, maxOffsetFails + ) + if maxOffsetFail + ] + ), + ) return ampsOffsets def getAmpEdges(self, im, amps, ampSides): @@ -494,8 +553,19 @@ def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB): interfaceOffset : float The calculated amp offset value for the given interface between amps A and B. + interfaceOffsetOriginal : float + The original calculated amp offset value for the given interface + between amps A and B. + ampEdgeGoodFrac : float + Fraction of viable pixel rows along the amp edge. + minFracFail : bool + True if the fraction of unmasked pixel rows is below the + ampEdgeMinFrac threshold. + maxOffsetFail : bool + True if the absolute offset value exceeds the ampEdgeMaxOffset + threshold. """ - interfaceId = f"{ampIdA}{ampIdB}" + interfaceId = f"{ampIdA:02d}:{ampIdB:02d}" sctrl = StatisticsControl() # NOTE: Taking the difference with the order below fixes the sign flip # in the B matrix. @@ -512,30 +582,26 @@ def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB): edgeDiffAvg[np.isnan(edgeDiff)] = np.nan # Take clipped mean of rolling average data as amp offset value. interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue() + interfaceOffsetOriginal = interfaceOffset ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg)) + # Perform a couple of do-no-harm safety checks: # a) The fraction of unmasked pixel rows is > ampEdgeMinFrac, # b) The absolute offset ADU value is < ampEdgeMaxOffset. minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset if minFracFail or maxOffsetFail: - log_fn = self.log.warning if self.config.doApplyAmpOffset else self.log.debug - if minFracFail: - log_fn( - f"The fraction of unmasked pixels for amp interface {interfaceId} is below the " - f"threshold ({ampEdgeGoodFrac:.2f} < {self.config.ampEdgeMinFrac}). Resetting the " - f"interface offset from {interfaceOffset} to 0." - ) - if maxOffsetFail: - log_fn( - "The absolute offset value exceeds the limit " - f"({np.abs(interfaceOffset):.2f} > {self.config.ampEdgeMaxOffset} ADU). Resetting " - f"the interface offset from {interfaceOffset} to 0." - ) interfaceOffset = 0 self.log.debug( f"amp interface {interfaceId} : " f"viable edge difference frac = {ampEdgeGoodFrac}, " f"interface offset = {interfaceOffset:.3f}" ) - return interfaceOffset + return ( + interfaceId, + interfaceOffset, + interfaceOffsetOriginal, + ampEdgeGoodFrac, + minFracFail, + maxOffsetFail, + )