Skip to content

Commit

Permalink
Simplify amp offset correction logging and config
Browse files Browse the repository at this point in the history
  • Loading branch information
@leeskelvin
leeskelvin committed Oct 30, 2024
1 parent bb4f38f commit 5410abc
Show file tree
Hide file tree
Showing 4 changed files with 118 additions and 49 deletions.
145 changes: 108 additions & 37 deletions python/lsst/ip/isr/ampOffset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.",
)


Expand Down Expand Up @@ -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."
Expand Down Expand Up @@ -243,31 +251,32 @@ def run(self, exposure):

metadata = exposure.getMetadata() # Exposure metadata.
self.metadata["AMPOFFSET_PEDESTALS"] = {} # Task metadata.
for amp, pedestal in zip(amps, pedestals):
ampName = amp.getName()
ampNames = [amp.getName() for amp in amps]
for ampName, amp, pedestal in zip(ampNames, amps, pedestals):
# Add the amp pedestal to the exposure metadata.
metadata.set(
f"LSST ISR AMPOFFSET PEDESTAL {ampName}",
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])}")
ampPedestalReport = ", ".join(
[f"{ampName}: {ampPedestal:.4f}" for (ampName, ampPedestal) in zip(ampNames, pedestals)]
)
self.log.info(f"Amp pedestal values ({status}): {ampPedestalReport}")

return Struct(pedestals=pedestals)

Expand Down Expand Up @@ -344,8 +353,9 @@ def getAmpAssociations(self, amps):
if not np.all(ampAssociations == ampAssociations.T):
raise RuntimeError("The `ampAssociations` is not symmetric about the diagonal.")

self.log.debug("amp associations:\n%s", ampAssociations)
self.log.debug("amp sides:\n%s", ampSides)
with np.printoptions(linewidth=200):
self.log.debug("amp associations:\n%s", ampAssociations)
self.log.debug("amp sides:\n%s", ampSides)

return ampAssociations, ampSides

Expand Down Expand Up @@ -408,8 +418,14 @@ def getAmpOffsets(self, im, amps, associations, sides):
"""
ampsOffsets = np.zeros(len(amps))
ampsEdges = self.getAmpEdges(im, amps, sides)
ampsNames = [amp.getName() for amp in amps]
interfaceOffsetLookup = {}

interfaceIds = []
interfaceOffsetOriginals = []
ampEdgeGoodFracs = []
minFracFails = []
maxOffsetFails = []
for ampId, ampAssociations in enumerate(associations):
ampNeighbors = np.ravel(np.where(ampAssociations < 0))
for ampNeighbor in ampNeighbors:
Expand All @@ -420,11 +436,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(ampsNames[ampId], ampsNames[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 offset quartile summary (min, Q1, Q2, Q3, max): %.4f, %.4f, %.4f, %.4f, %.4f",
*quartile_summary,
)
log_fn = self.log.info 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(
"Absolute amp offsets exceed the configured maximum (%s) and have been set to zero for the "
"following amp interfaces: %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):
Expand Down Expand Up @@ -474,16 +538,16 @@ def getAmpEdges(self, im, amps, ampSides):
ampEdges[ampId][ampSide] = np.nanmedian(strip, axis=ampSide % 2) # 1D medianified strip
return ampEdges

def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB):
def getInterfaceOffset(self, ampNameA, ampNameB, edgeA, edgeB):
"""Calculate the amp offset for a given interface between two
amplifiers.
Parameters
----------
ampIdA : int
ID of the first amplifier.
ampIdB : int
ID of the second amplifier.
ampNameA : int
Name of the first amplifier.
ampNameB : int
Name of the second amplifier.
edgeA : numpy.ndarray
Amp edge for the first amplifier.
edgeB : numpy.ndarray
Expand All @@ -494,8 +558,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"{ampNameA}:{ampNameB}"
sctrl = StatisticsControl()
# NOTE: Taking the difference with the order below fixes the sign flip
# in the B matrix.
Expand All @@ -512,30 +587,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}"
f"viable edge difference frac = {ampEdgeGoodFrac:.2f}, "
f"amp offset = {interfaceOffsetOriginal:.3f}"
)
return (
interfaceId,
interfaceOffset,
interfaceOffsetOriginal,
ampEdgeGoodFrac,
minFracFail,
maxOffsetFail,
)
return interfaceOffset
10 changes: 5 additions & 5 deletions python/lsst/ip/isr/isrTask.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -987,8 +987,8 @@ def validate(self):
self.maskListToInterpolate.remove(self.saturatedMaskName)
if self.doNanInterpolation and "UNMASKEDNAN" not in self.maskListToInterpolate:
self.maskListToInterpolate.append("UNMASKEDNAN")
if self.ampOffset.doApplyAmpOffset and not self.doAmpOffset:
raise ValueError("ampOffset.doApplyAmpOffset requires doAmpOffset to be True.")
if self.ampOffset.measureOnly and not self.doAmpOffset:
raise ValueError("ampOffset.measureOnly requires doAmpOffset to be True.")


class IsrTask(pipeBase.PipelineTask):
Expand Down Expand Up @@ -1730,10 +1730,10 @@ def run(self, ccdExposure, *, camera=None, bias=None, linearizer=None,

# Calculate amp offset corrections within the CCD.
if self.config.doAmpOffset:
if self.config.ampOffset.doApplyAmpOffset:
self.log.info("Calculating and applying amp offset corrections.")
if self.config.ampOffset.measureOnly:
self.log.info("Measuring amp offset corrections only, without applying them.")
else:
self.log.info("Calculating amp offset corrections without applying them.")
self.log.info("Measuring and applying amp offset corrections.")
self.ampOffset.run(ccdExposure)

if self.config.doMeasureBackground:
Expand Down
6 changes: 3 additions & 3 deletions python/lsst/ip/isr/isrTaskLSST.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1825,10 +1825,10 @@ def run(self, ccdExposure, *, dnlLUT=None, bias=None, deferredChargeCalib=None,

# Calculate amp offset corrections within the CCD.
if self.config.doAmpOffset:
if self.config.ampOffset.doApplyAmpOffset:
self.log.info("Calculating and applying amp offset corrections.")
if self.config.ampOffset.measureOnly:
self.log.info("Measuring amp offset corrections only, without applying them.")
else:
self.log.info("Calculating amp offset corrections without applying them.")
self.log.info("Measuring and applying amp offset corrections.")
self.ampOffset.run(ccdExposure)

# Calculate standard image quality statistics
Expand Down
6 changes: 2 additions & 4 deletions tests/test_ampOffset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,6 @@ def runAmpOffsetWithBackground(self, valueType, rampBackground=False):
config.doDetection = True
config.ampEdgeWidth = 12
config.applyWeights = applyWeights
config.doApplyAmpOffset = True # Updates the exposure in place.
if valueType == "random":
# For this specific case, the fraction of unmasked pixels for
# amp interface 01 is unusually small.
Expand Down Expand Up @@ -376,13 +375,13 @@ def testAmpOffsetEffectOnExposure(self):

# Configure to not apply amp offset to the exposure and run the task.
# Verify that the exposure remains unchanged.
config.doApplyAmpOffset = False
config.measureOnly = True
AmpOffsetTask(config=config).run(exp)
self.assertFloatsEqual(exp0.image.array, exp.image.array)

# Configure to apply amp offset to the exposure and run the task.
# Verify that the exposure is updated.
config.doApplyAmpOffset = True
config.measureOnly = False
AmpOffsetTask(config=config).run(exp)
self.assertFloatsNotEqual(exp0.image.array, exp.image.array)

Expand All @@ -394,7 +393,6 @@ def testAmpOffset(self, valueType):
config.doBackground = False
config.doDetection = False
config.ampEdgeWidth = 12 # Given 100x51 amps in our mock detector.
config.doApplyAmpOffset = True # Updates the exposure in place.
if valueType == "artificial":
# For this extreme case, we expect the interface offsets to be
# unusually large.
Expand Down

0 comments on commit 5410abc

Please sign in to comment.