From b7a42ba087540e6aea5268b7a36edfecc295b0a9 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 1 Aug 2025 14:17:19 -0700 Subject: [PATCH 1/8] Add crosstalk pipeline definition. --- bps/templates/bps_crosstalk.yaml | 19 ++++++++++++ pipelines/LSSTCam/cpCrosstalk.yaml | 4 +++ pipelines/_ingredients/cpCrosstalkLSST.yaml | 32 +++++++++++++++++++++ 3 files changed, 55 insertions(+) create mode 100644 bps/templates/bps_crosstalk.yaml create mode 100644 pipelines/LSSTCam/cpCrosstalk.yaml create mode 100644 pipelines/_ingredients/cpCrosstalkLSST.yaml diff --git a/bps/templates/bps_crosstalk.yaml b/bps/templates/bps_crosstalk.yaml new file mode 100644 index 000000000..059a614ba --- /dev/null +++ b/bps/templates/bps_crosstalk.yaml @@ -0,0 +1,19 @@ +pipelineYaml: "${CP_PIPE_DIR}/pipelines/${INSTRUMENT}/cpCrosstalk.yaml" + +# Necessary to cluster by exposure if there are many input exposures +clusterAlgorithm: lsst.ctrl.bps.quantum_clustering_funcs.dimension_clustering +cluster: + cpBiasIsrByExposure: + pipetasks: cpCrosstalkIsr + dimensions: exposure + +project: "${TICKET}" +campaign: "${TICKET}" +submitPath: "${SCRATCH}/submit/${outputRun}" + +payload: + payloadName: "${INSTRUMENT}_${TICKET}_crosstalk" + output: "${USER_CALIB_PREFIX}${INSTRUMENT}/calib/${TICKET}/${TAG}/crosstalkGen.${RERUN}" + butlerConfig: "${REPO}" + inCollection: "${USER_CALIB_PREFIX}${INSTRUMENT}/calib/${TICKET},${RAW_COLLECTION},${CALIB_COLLECTIONS}" + dataQuery: "${SELECTION_CROSSTALK}" diff --git a/pipelines/LSSTCam/cpCrosstalk.yaml b/pipelines/LSSTCam/cpCrosstalk.yaml new file mode 100644 index 000000000..b9e6c19d1 --- /dev/null +++ b/pipelines/LSSTCam/cpCrosstalk.yaml @@ -0,0 +1,4 @@ +description: cp_pipe LSSTCam crosstalk calibration construction +instrument: lsst.obs.lsst.LsstCam +imports: + - location: $CP_PIPE_DIR/pipelines/_ingredients/cpCrosstalkLSST.yaml diff --git a/pipelines/_ingredients/cpCrosstalkLSST.yaml b/pipelines/_ingredients/cpCrosstalkLSST.yaml new file mode 100644 index 000000000..f55866bc1 --- /dev/null +++ b/pipelines/_ingredients/cpCrosstalkLSST.yaml @@ -0,0 +1,32 @@ +description: cp_pipe CROSSTALK calibration construction. +tasks: + cpCrosstalkIsr: + class: lsst.ip.isr.IsrTaskLSST + config: + connections.ccdExposure: "raw" + connections.outputExposure: "cpCrosstalkIsrExp" + python: | + from lsst.cp.pipe import configureIsrTaskLSSTForCalibrations + configureIsrTaskLSSTForCalibrations(config) + + config.doCrosstalk = False + config.doApplyGains = True + cpCrosstalkExtract: + class: lsst.cp.pipe.CrosstalkExtractTask + config: + connections.inputExp: "cpCrosstalkIsrExp" + connections.outputRatios: "cpCrosstalkRatio" + cpCrosstalkSolve: + class: lsst.cp.pipe.CrosstalkSolveTask + config: + connections.inputRatios: "cpCrosstalkRatio" + connections.outputCrosstalk: "crosstalk" +subsets: + cpCrosstalk: + subset: + - cpCrosstalkExtract + - cpCrosstalkSolve + description: Just the extract and solve. +contracts: + - cpCrosstalkExtract.connections.inputExp == cpCrosstalkIsr.connections.outputExposure + - cpCrosstalkSolve.connections.inputRatios == cpCrosstalkExtract.connections.outputRatios From a777d4bb881c0ad4912133d2c6bae2cb444b7f21 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Thu, 14 Aug 2025 15:25:02 -0700 Subject: [PATCH 2/8] Improvements for crosstalk: a) check that flux and ratio length match b) pass fluxes to measurement code, to allow for preliminary non-linear fitting c) log clarifications d) Add stage subsets. --- pipelines/LSSTCam/cpCrosstalk.yaml | 11 ++ pipelines/_ingredients/cpCrosstalk.yaml | 2 + pipelines/_ingredients/cpCrosstalkLSST.yaml | 6 + python/lsst/cp/pipe/cpCrosstalk.py | 170 +++++++++++++++++--- tests/test_measureCrosstalk.py | 2 + 5 files changed, 173 insertions(+), 18 deletions(-) diff --git a/pipelines/LSSTCam/cpCrosstalk.yaml b/pipelines/LSSTCam/cpCrosstalk.yaml index b9e6c19d1..cf4498fc8 100644 --- a/pipelines/LSSTCam/cpCrosstalk.yaml +++ b/pipelines/LSSTCam/cpCrosstalk.yaml @@ -2,3 +2,14 @@ description: cp_pipe LSSTCam crosstalk calibration construction instrument: lsst.obs.lsst.LsstCam imports: - location: $CP_PIPE_DIR/pipelines/_ingredients/cpCrosstalkLSST.yaml +tasks: + cpCrosstalkIsr: + class: lsst.ip.isr.IsrTaskLSST + config: + doDefect: true + cpCrosstalkExtract: + class: lsst.cp.pipe.CrosstalkExtractTask + cpCrosstalkSolve: + class: lsst.cp.pipe.CrosstalkSolveTask + config: + fluxOrder: 1 diff --git a/pipelines/_ingredients/cpCrosstalk.yaml b/pipelines/_ingredients/cpCrosstalk.yaml index f1674bad8..db4b9f30c 100644 --- a/pipelines/_ingredients/cpCrosstalk.yaml +++ b/pipelines/_ingredients/cpCrosstalk.yaml @@ -26,6 +26,8 @@ tasks: config: connections.inputExp: "cpCrosstalkIsrExp" connections.outputRatios: "cpCrosstalkRatio" + background.useApprox: false + background.ignoredPixelMask: ["BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA", "CT_TEMP"] cpCrosstalkSolve: class: lsst.cp.pipe.CrosstalkSolveTask config: diff --git a/pipelines/_ingredients/cpCrosstalkLSST.yaml b/pipelines/_ingredients/cpCrosstalkLSST.yaml index f55866bc1..1a2c3e651 100644 --- a/pipelines/_ingredients/cpCrosstalkLSST.yaml +++ b/pipelines/_ingredients/cpCrosstalkLSST.yaml @@ -21,12 +21,18 @@ tasks: config: connections.inputRatios: "cpCrosstalkRatio" connections.outputCrosstalk: "crosstalk" + rejectNegativeSolutions: false subsets: cpCrosstalk: subset: - cpCrosstalkExtract - cpCrosstalkSolve description: Just the extract and solve. + cpCrosstalkProc: + subset: + - cpCrosstalkIsr + - cpCrosstalkExtract + description: Handle the per-exposure analysis contracts: - cpCrosstalkExtract.connections.inputExp == cpCrosstalkIsr.connections.outputExposure - cpCrosstalkSolve.connections.inputRatios == cpCrosstalkExtract.connections.outputRatios diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index ac80ffd5b..67a423a1a 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -22,11 +22,13 @@ import numpy as np from collections import defaultdict +from copy import copy import lsst.pipe.base as pipeBase import lsst.pipe.base.connectionTypes as cT from lsstDebug import getDebugFrame +from lsst.afw.cameraGeom import ReadoutCorner from lsst.afw.detection import FootprintSet, Threshold from lsst.afw.display import getDisplay from lsst.pex.config import ConfigurableField, Field, ListField @@ -113,6 +115,12 @@ class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig, doc="Background estimation task.", ) + def setDefaults(self): + super().setDefaults() + # is this really the best way to do this?? + self.background.useApprox = False + self.background.ignoredPixelMask.append("CT_TEMP") + def validate(self): super().validate() @@ -190,15 +198,45 @@ def run(self, inputExp, sourceExps=[]): self.log.info("Measuring full detector background for target: %s", targetChip) targetIm = inputExp.getMaskedImage() + FootprintSet(targetIm, Threshold(threshold), "DETECTED") detected = targetIm.getMask().getPlaneBitMask("DETECTED") + targetIm.getMask().addMaskPlane("CT_TEMP") + maskBit = targetIm.getMask().getPlaneBitMask("CT_TEMP") bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"]) + + # Carry over over-threshold masked pixels to other amplifiers. + for amp in targetDetector: + ampIm = inputExp[amp.getBBox()] + ampName = amp.getName() + + mask = ampIm.mask.array & detected + if np.sum(mask) == 0: + continue + + newMask = np.where(np.bitwise_and(mask, detected), maskBit, 0) + for ampToMask in targetDetector: + if ampName == ampToMask.getName(): + # The amp we're considering already + continue + + extractedAmp = inputExp[ampToMask.getBBox()] + # The mask needs to be flipped to match the target + flippedMask = self._flipMask(newMask, amp, ampToMask) + + extractedAmp.mask.array[:, :] |= flippedMask + + # We've now masked the source pixels, and any potential CT + # pixels, so this should be just the + # background/reflections/etc. backgroundModel = self.background.fitBackground(inputExp.maskedImage) backgroundIm = backgroundModel.getImageF() - self.debugView('extract', inputExp) + # Begin search for bright pixels, and their associated crosstalk + # signals. for sourceExp in sourceExtractExps: + # This loop exists to support future inter-chip searches. sourceDetector = sourceExp.getDetector() sourceChip = sourceDetector.getName() sourceIm = sourceExp.getMaskedImage() @@ -211,7 +249,7 @@ def run(self, inputExp, sourceExps=[]): # The dictionary of amp-to-amp ratios for this pair of # source->target detectors. - ratioDict = defaultdict(lambda: defaultdict(list)) + ratioDict = defaultdict(lambda: defaultdict(np.array)) extractedCount = 0 for sourceAmp in sourceDetector: @@ -267,7 +305,58 @@ def run(self, inputExp, sourceExps=[]): outputFluxes=ddict2dict(outputFluxes) ) + def _flipMask(self, maskArray, sourceAmp, targetAmp): + """Flip an array from a sourceAmp to match the readout order of + targetAmp. + + Parameters + ---------- + maskArray : `np.ndarray` + Mask data to flip. + sourceAmp : `lsst.afw.cameraGeom.Amplifier` + Amplifier corresponding to the maskArray. + targetAmp : `lsst.afw.cameraGeom.Amplifier` + Amplifier corresponding to the output mask. + + Returns + ------- + maskFlipped : `np.ndarray` + The flipped mask. + + See Also + ----- + lsst.ip.isr.CrosstalkCalib.extractAmp() + """ + maskFlipped = copy(maskArray) + if sourceAmp.getReadoutCorner() == targetAmp.getReadoutCorner(): + return maskFlipped + + X_FLIP = {ReadoutCorner.LL: False, + ReadoutCorner.LR: True, + ReadoutCorner.UL: False, + ReadoutCorner.UR: True} + Y_FLIP = {ReadoutCorner.LL: False, + ReadoutCorner.LR: False, + ReadoutCorner.UL: True, + ReadoutCorner.UR: True} + + sourceAmpCorner = sourceAmp.getReadoutCorner() + targetAmpCorner = targetAmp.getReadoutCorner() + + # Flipping is necessary only if the desired configuration doesn't match + # what we currently have. + xFlip = X_FLIP[targetAmpCorner] ^ X_FLIP[sourceAmpCorner] + yFlip = Y_FLIP[targetAmpCorner] ^ Y_FLIP[sourceAmpCorner] + + if xFlip: + maskFlipped = np.fliplr(maskFlipped) + if yFlip: + maskFlipped = np.flipud(maskFlipped) + + return maskFlipped + def debugView(self, stepname, exposure): + """Utility function to examine the image being processed. Parameters @@ -335,6 +424,7 @@ class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections, storageClass="StructuredDataDict", dimensions=("instrument", "exposure", "detector"), multiple=True, + deferLoad=True, ) inputFluxes = cT.Input( name="crosstalkFluxes", @@ -342,6 +432,7 @@ class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections, storageClass="StructuredDataDict", dimensions=("instrument", "exposure", "detector"), multiple=True, + deferLoad=True, ) camera = cT.PrerequisiteInput( name="camera", @@ -363,8 +454,8 @@ class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections, def __init__(self, *, config=None): super().__init__(config=config) - if config.fluxOrder == 0: - self.inputs.discard("inputFluxes") + # if config.fluxOrder == 0 and False: + # self.inputs.discard("inputFluxes") class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig, @@ -385,12 +476,12 @@ class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig, fluxOrder = Field( dtype=int, default=0, - doc="Polynomial order in source flux to fit crosstalk.", + doc="Order of source flux fit to crosstalk. 0=simple linear; 1=first order non-linear.", ) rejectNegativeSolutions = Field( dtype=bool, - default=True, + default=False, doc="Should solutions with negative coefficients (which add flux to the target) be excluded?", ) @@ -494,11 +585,18 @@ def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, output if inputFluxes is None: inputFluxes = [None for exp in inputRatios] + if inputDims is None: + inputDims = [{} for exp in inputRatios] combinedRatios = defaultdict(lambda: defaultdict(list)) combinedFluxes = defaultdict(lambda: defaultdict(list)) - for ratioDict, fluxDict in zip(inputRatios, inputFluxes): + for ratioRef, fluxRef, dimensions in zip(inputRatios, inputFluxes, inputDims): + ratioDict = ratioRef.get() + if fluxRef is not None: + fluxDict = fluxRef.get() + else: + fluxDict = None for targetChip in ratioDict: if calibChip and targetChip != calibChip and targetChip != calibDetector.getName(): raise RuntimeError(f"Target chip: {targetChip} does not match calibration dimension: " @@ -510,9 +608,23 @@ def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, output for targetAmp in ratios: for sourceAmp in ratios[targetAmp]: - combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp]) if fluxDict: + if len(ratios[targetAmp][sourceAmp]) != len(fluxDict[sourceChip][sourceAmp]): + if targetAmp != sourceAmp: + # This usually triggers when + # targetAmp == sourceAmp. The + # sourceAmp has flux entries, + # but by definition cannot + # have ratios for itself. Only + # warn when this isn't the + # case. + self.log.warning(f"Length mismatch for ratios {len(ratios[targetAmp][sourceAmp])} " # noqa E501 + f"and fluxes {len(fluxDict[sourceChip][sourceAmp])} for " # noqa E501 + f"Source {sourceAmp} and target {targetAmp} " + f"Rejecting this {dimensions}") + continue combinedFluxes[targetAmp][sourceAmp].extend(fluxDict[sourceChip][sourceAmp]) + combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp]) # TODO: DM-21904 # Iterating over all other entries in # ratioDict[targetChip] will yield inter-chip terms. @@ -525,13 +637,14 @@ def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, output if len(combinedRatios[targetAmp][sourceAmp]) > 1: self.debugRatios('reduce', combinedRatios, targetAmp, sourceAmp) - if self.config.fluxOrder == 0: - self.log.info("Fitting crosstalk coefficients.") + if self.config.fluxOrder < 2: + self.log.info("Fitting crosstalk coefficients with order {self.config.fluxOrder}") calib = self.measureCrosstalkCoefficients(combinedRatios, ordering, + combinedFluxes, self.config.rejIter, self.config.rejSigma) else: - raise NotImplementedError("Non-linear crosstalk terms are not yet supported.") + raise NotImplementedError("Higher order non-linear crosstalk terms are not yet supported.") self.log.info("Number of valid coefficients: %d", np.sum(calib.coeffValid)) @@ -561,7 +674,7 @@ def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, output outputProvenance=provenance, ) - def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma): + def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSigma): """Measure crosstalk coefficients from the ratios. Given a list of ratios for each target/source amp combination, @@ -573,13 +686,15 @@ def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma): Parameters ---------- ratios : `dict` [`dict` [`numpy.ndarray`]] - Catalog of arrays of ratios. The ratio arrays are one-dimensional + Catalog of arrays of ratios. The ratio arrays are one-dimensional. ordering : `list` [`str`] or None List to use as a mapping between amplifier names (the elements of the list) and their position in the output calibration (the matching index of the list). If no ordering is supplied, the order of the keys in the ratio catalog is used. + fluxes : `dict` [`dict` [`numpy.ndarray`]] + Catalog of arrays of fluxes. The flux arrays are one-dimensional. rejIter : `int` Number of rejection iterations. rejSigma : `float` @@ -598,12 +713,18 @@ def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma): # Calibration stores coefficients as a numpy ndarray. for ss, tt in itertools.product(range(calib.nAmp), range(calib.nAmp)): if ss == tt: - values = [0.0] + values = [] + myfluxes = [] else: # ratios is ratios[Target][Source] # use tt for Target, use ss for Source, to match ip_isr. values = np.array(ratios[ordering[tt]][ordering[ss]]) - values = values[np.abs(values) < 1.0] # Discard unreasonable values + good_values = np.abs(values) < 1.0 # Discard unreasonable values + values = values[good_values] + myfluxes = np.array(fluxes[ordering[tt]][ordering[ss]]) + myfluxes = myfluxes[good_values] + if len(values) != len(myfluxes): + self.log.warning(f"Flux and ratio length disagree after first filter: {len(values)} {len(myfluxes)}") # noqa E501 # Sigma clip using the inter-quartile distance and a # normal distribution. @@ -617,6 +738,9 @@ def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma): if good.sum() == len(good) or good.sum() == 0: break values = values[good] + myfluxes = myfluxes[good] + if len(values) != len(myfluxes): + self.log.warning(f"Flux and ratio length disagree after second filter: {len(values)} {len(myfluxes)}") # noqa E501 # Crosstalk calib is property[Source][Target]. calib.coeffNum[ss][tt] = len(values) @@ -626,13 +750,22 @@ def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma): calib.coeffs[ss][tt] = np.nan calib.coeffErr[ss][tt] = np.nan calib.coeffValid[ss][tt] = False + polyfit = [0.0, 0.0] else: calib.coeffs[ss][tt] = np.mean(values) + polyfit = np.polyfit(myfluxes, values, 1) + + if self.config.fluxOrder == 1: + # substitute polyfit solution. + calib.coeffs[ss][tt] = polyfit[1] + calib.coeffsSqr[ss][tt] = polyfit[0] + if self.config.rejectNegativeSolutions and calib.coeffs[ss][tt] < 0.0: calib.coeffs[ss][tt] = 0.0 - if calib.coeffNum[ss][tt] == 1: + if calib.coeffNum[ss][tt] <= 1: calib.coeffErr[ss][tt] = np.nan + calib.coeffSqr[ss][tt] = np.nan calib.coeffValid[ss][tt] = False else: correctionFactor = sigmaClipCorrection(rejSigma) @@ -646,9 +779,10 @@ def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma): calib.coeffValid[ss][tt] = np.abs(calib.coeffs[ss][tt]) > significanceThreshold self.debugRatios('measure', ratios, ordering[ss], ordering[tt], calib.coeffs[ss][tt], calib.coeffValid[ss][tt]) - self.log.info("Measured %s -> %s Coeff: %e Err: %e N: %d Valid: %s Limit: %e", + self.log.info("Measured %s -> %s Coeff: %e Err: %e N: %d Valid: %s Limit: %e Quadratic: %s", ordering[ss], ordering[tt], calib.coeffs[ss][tt], calib.coeffErr[ss][tt], - calib.coeffNum[ss][tt], calib.coeffValid[ss][tt], significanceThreshold) + calib.coeffNum[ss][tt], calib.coeffValid[ss][tt], significanceThreshold, + polyfit) return calib diff --git a/tests/test_measureCrosstalk.py b/tests/test_measureCrosstalk.py index 03c2875eb..88dbb3cca 100644 --- a/tests/test_measureCrosstalk.py +++ b/tests/test_measureCrosstalk.py @@ -86,6 +86,8 @@ def setup_measureCrosstalk(self, isTrimmed=False, nSources=8): result = ctex.run(exposure) fullResult.append(result.outputRatios) + # This test needs fixing: CZW: + return [True, True, True] # Generate the final measured CT ratios, uncertainties, pixel counts. ctsConfig = CrosstalkSolveConfig() cts = CrosstalkSolveTask(config=ctsConfig) From cb4a7e5027b7fab382ab60ad348277364c97540b Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Wed, 12 Nov 2025 12:22:13 -0800 Subject: [PATCH 3/8] Add crosstalk filtering step. --- pipelines/LSSTCam/cpCrosstalk.yaml | 2 + pipelines/_ingredients/cpCrosstalkLSST.yaml | 8 +- python/lsst/cp/pipe/cpCrosstalk.py | 345 +++++++++++++++++++- 3 files changed, 351 insertions(+), 4 deletions(-) diff --git a/pipelines/LSSTCam/cpCrosstalk.yaml b/pipelines/LSSTCam/cpCrosstalk.yaml index cf4498fc8..9768be939 100644 --- a/pipelines/LSSTCam/cpCrosstalk.yaml +++ b/pipelines/LSSTCam/cpCrosstalk.yaml @@ -9,6 +9,8 @@ tasks: doDefect: true cpCrosstalkExtract: class: lsst.cp.pipe.CrosstalkExtractTask + config: + growMaskRadius: 20 cpCrosstalkSolve: class: lsst.cp.pipe.CrosstalkSolveTask config: diff --git a/pipelines/_ingredients/cpCrosstalkLSST.yaml b/pipelines/_ingredients/cpCrosstalkLSST.yaml index 1a2c3e651..e1f977fd6 100644 --- a/pipelines/_ingredients/cpCrosstalkLSST.yaml +++ b/pipelines/_ingredients/cpCrosstalkLSST.yaml @@ -20,8 +20,13 @@ tasks: class: lsst.cp.pipe.CrosstalkSolveTask config: connections.inputRatios: "cpCrosstalkRatio" - connections.outputCrosstalk: "crosstalk" + connections.outputCrosstalk: "crosstalkProposal" rejectNegativeSolutions: false + cpCrosstalkFilter: + class: lsst.cp.pipe.CrosstalkFilterTask + config: + connections.inputCrosstalk: "crosstalkProposal" + connections.outputCrosstalk: "crosstalk" subsets: cpCrosstalk: subset: @@ -36,3 +41,4 @@ subsets: contracts: - cpCrosstalkExtract.connections.inputExp == cpCrosstalkIsr.connections.outputExposure - cpCrosstalkSolve.connections.inputRatios == cpCrosstalkExtract.connections.outputRatios + - cpCrosstalkSolve.connections.outputCrosstalk == cpCrosstalkFilter.connections.inputCrosstalk diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index 67a423a1a..590e1c3a2 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -23,6 +23,7 @@ from collections import defaultdict from copy import copy +from scipy.stats import median_abs_deviation import lsst.pipe.base as pipeBase import lsst.pipe.base.connectionTypes as cT @@ -32,12 +33,13 @@ from lsst.afw.detection import FootprintSet, Threshold from lsst.afw.display import getDisplay from lsst.pex.config import ConfigurableField, Field, ListField -from lsst.ip.isr import CrosstalkCalib, IsrProvenance +from lsst.ip.isr import CrosstalkCalib, IsrProvenance, growMasks from lsst.cp.pipe.utils import (ddict2dict, sigmaClipCorrection) from lsst.meas.algorithms import SubtractBackgroundTask __all__ = ["CrosstalkExtractConfig", "CrosstalkExtractTask", - "CrosstalkSolveTask", "CrosstalkSolveConfig"] + "CrosstalkSolveTask", "CrosstalkSolveConfig", + "CrosstalkFilterTask", "CrosstalkFilterConfig"] class CrosstalkExtractConnections(pipeBase.PipelineTaskConnections, @@ -114,6 +116,11 @@ class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig, target=SubtractBackgroundTask, doc="Background estimation task.", ) + growMaskRadius = Field( + dtype=int, + default=0, + doc="Radius to grow CT_TEMP masks prior to background estimation." + ) def setDefaults(self): super().setDefaults() @@ -225,6 +232,9 @@ def run(self, inputExp, sourceExps=[]): flippedMask = self._flipMask(newMask, amp, ampToMask) extractedAmp.mask.array[:, :] |= flippedMask + # Optionally dilate these masks by some amount: + growMasks(inputExp.mask, radius=self.config.growMaskRadius, + maskNameList=["CT_TEMP"], maskValue="CT_TEMP") # We've now masked the source pixels, and any potential CT # pixels, so this should be just the @@ -443,7 +453,7 @@ class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections, ) outputCrosstalk = cT.Output( - name="crosstalk", + name="crosstalkProposal", doc="Output proposed crosstalk calibration.", storageClass="CrosstalkCalib", dimensions=("instrument", "detector"), @@ -884,3 +894,332 @@ def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False): import pdb pdb.set_trace() plt.close() + + +class CrosstalkFilterConnections(pipeBase.PipelineTaskConnections, + dimensions=("instrument", )): + inputCrosstalk = cT.Input( + name="crosstalkProposal", + doc="Input crosstalk calibrations as measured by CrosstalkSolveTask.", + storageClass="IsrCalib", + dimensions=("instrument", "detector"), + isCalibration=True, + multiple=True, + ) + + camera = cT.PrerequisiteInput( + name="camera", + doc="Camera containing cameraGeom information.", + storageClass="Camera", + dimensions=("instrument", ), + isCalibration=True, + ) + + outputCrosstalk = cT.Output( + name="crosstalk", + doc="Filtered crosstalk solutions.", + storageClass="CrosstalkCalib", + dimensions=("instrument", "detector"), + isCalibration=True, + multiple=True, + ) + + +class CrosstalkFilterConfig(pipeBase.PipelineTaskConfig, + pipelineConnections=CrosstalkFilterConnections): + """Configuration for the filtering of measured crosstalk solutions.""" + doFiltering = Field( + dtype=bool, + default=True, + doc="Do filtering? If false, then this task acts as a pass-through to rename dataset types.", + ) + + nSigmaCoeffClip = Field( + dtype=float, + default=3.0, + doc="Coefficient outlier clipping significance.", + ) + nSigmaCoeffSqrClip = Field( + dtype=float, + default=6.0, + doc="Squared-term coefficient outlier clipping significance.", + ) + + +class CrosstalkFilterTask(pipeBase.PipelineTask): + """Task to compare crosstalk solutions between detectors, to identify + and remove outliers. + """ + + ConfigClass = CrosstalkFilterConfig + _DefaultName = 'cpCrosstalkFilter' + + def runQuantum(self, butlerQC, inputRefs, outputRefs): + """Ensure that the input and output dimensions are passed along. + + Parameters + ---------- + butlerQC : `lsst.daf.butler.QuantumContext` + Butler to operate on. + inputRefs : `lsst.pipe.base.InputQuantizedConnection` + Input data refs to load. + ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection` + Output data refs to persist. + """ + inputs = butlerQC.get(inputRefs) + + # Use the dimensions to set calib/provenance information. + inputs['inputDims'] = [dict(inCT.dataId.required) for inCT in inputRefs.inputCrosstalk] + inputs['outputDims'] = [dict(outCT.dataId.required) for outCT in outputRefs.outputCrosstalk] + + outputs = self.run(**inputs) + butlerQC.put(outputs, outputRefs) + + def run(self, inputCrosstalk, camera, inputDims, outputDims): + """Compare crosstalk solutions to produce filtered crosstalk + calibrations. + + Parameters + ---------- + inputCrosstalk : `list` [`lsst.ip.isr.CrosstalkCalib`] + List of crosstalk solutions to filter. + camera : `lsst.afw.cameraGeom.Camera` + Input camera. + inputDims : `list` [`lsst.daf.butler.DataCoordinate`] + DataIds to use to construct provenance. + outputDims : `list` [`lsst.daf.butler.DataCoordinate`] + DataIds to use to populate the output calibration. + + Returns + ------- + results : `lsst.pipe.base.Struct` + The results struct containing: + + ``outputCrosstalk`` + Final crosstalk calibration + (`lsst.ip.isr.CrosstalkCalib`). + + Raises + ------ + RuntimeError + Raised if something goes bad. CZW/Fix me. + """ + # These will hold all of the input data. + itl_c0 = [] + e2v_c0 = [] + itl_c1 = [] + e2v_c1 = [] + detector_map = {} + itl_counter = 0 + e2v_counter = 0 + for inputCT, inputDim in zip(inputCrosstalk, inputDims): + detId = inputDim['detector'] + detector = camera[detId] + + if detector.getPhysicalType() == 'ITL': + itl_c0.append(inputCT.coeffs) + itl_c1.append(inputCT.coeffsSqr) + + detector_map[detId] = itl_counter + itl_counter += 1 + + elif detector.getPhysicalType() == 'E2V': + e2v_c0.append(inputCT.coeffs) + e2v_c1.append(inputCT.coeffsSqr) + + detector_map[detId] = e2v_counter + e2v_counter += 1 + else: + # This is a wavefront sensor, and we don't want to + # filter those. + pass + + itl_c0 = np.array(itl_c0) + itl_c1 = np.array(itl_c1) + e2v_c0 = np.array(e2v_c0) + e2v_c1 = np.array(e2v_c1) + + itl_outliers = self.find_outliers(itl_c0, itl_c1) + e2v_outliers = self.find_outliers(e2v_c0, e2v_c1) + + if self.config.doFiltering: + itl_final = self.replace_outliers(itl_c0, itl_c1, + itl_outliers.isBad, itl_outliers.median0, itl_outliers.median1) + e2v_final = self.replace_outliers(e2v_c0, e2v_c1, + e2v_outliers.isBad, e2v_outliers.median0, e2v_outliers.median1) + outputCrosstalkList = [] + for inputCT, inputDim, outputDim in zip(inputCrosstalk, inputDims, outputDims): + if inputDim['detector'] != outputDim['detector']: + raise RuntimeError("Inconsistent dimension records") + + detId = inputDim['detector'] + detector = camera[detId] + + outputCT = copy(inputCT) + + if detector.getPhysicalType() == 'ITL': + itl_index = detector_map[detId] + if (np.any(itl_final.new_matrix0[itl_index] != outputCT.coeffs) + or np.any(itl_final.new_matrix1[itl_index] != outputCT.coeffsSqr)): + outputCT.coeffs = itl_final.new_matrix0[itl_index] + outputCT.coeffsSqr = itl_final.new_matrix1[itl_index] + elif detector.getPhysicalType() == 'E2V': + e2v_index = detector_map[detId] + if (np.any(e2v_final.new_matrix0[e2v_index] != outputCT.coeffs) + or np.any(e2v_final.new_matrix1[e2v_index] != outputCT.coeffsSqr)): + outputCT.coeffs = e2v_final.new_matrix0[e2v_index] + outputCT.coeffsSqr = e2v_final.new_matrix1[e2v_index] + + outputCT.updateMetadata( + camera=camera, + detector=camera[detId], + setCalibId=True, + setCalibInfo=True, + setDate=True, + ) + outputCrosstalkList.append(outputCT) + return pipeBase.Struct( + outputCrosstalk=outputCrosstalkList, + ) + + def find_outliers(self, matrix0, matrix1): + """Do checks to see if an element of the matrix is out-of-family. + + Parameters + ---------- + matrix0 : `np.array`, (Ndet, Namp, Namp) + Matrix holding the 0th-order terms. + matrix1 : `np.array`, (Ndet, Namp, Namp) + Matrix holding the 1st-order terms. + + Returns + ------- + results : `lsst.pipe.base.Struct` + Results struct containing + + ``median0`` + Median in-family value (`np.array` (Namp, Namp)). + ``stdev0`` + MAD effective sigma in-family value (`np.array` (Namp, Namp)). + ``median1`` + Median in-family value (`np.array` (Namp, Namp)). + ``stdev1`` + MAD effective sigma in-family value (`np.array` (Namp, Namp)). + ``isBad`` + Boolean indicator that an element has been replaced + (`np.array` (Ndet, Namp, Namp)). + + Raises + ------ + ValueError : + Raised if the inputs have a mismatch in size. + """ + if matrix0.shape != matrix1.shape: + raise ValueError("Shape disagreement!") + if matrix0.shape[1] != matrix0.shape[2]: + raise ValueError("Shape disagreement!") + + nAmp = matrix0.shape[1] + + median0 = np.nanmedian(matrix0, axis=0) + median1 = np.nanmedian(matrix1, axis=0) + + stdev0 = median_abs_deviation(matrix0, axis=0, + center=np.nanmedian, scale="normal", nan_policy="omit") + stdev1 = median_abs_deviation(matrix1, axis=0, + center=np.nanmedian, scale="normal", nan_policy="omit") + + isBad = np.full_like(matrix0, False, dtype=bool) + + for i in range(nAmp): + for j in range(nAmp): + m0 = median0[i][j] + m1 = median1[i][j] + s0 = stdev0[i][j] + s1 = stdev1[i][j] + + min0 = m0 - self.config.nSigmaCoeffClip * s0 + max0 = m0 + self.config.nSigmaCoeffClip * s0 + + min1 = m1 - self.config.nSigmaCoeffSqrClip * s1 + max1 = m1 + self.config.nSigmaCoeffSqrClip * s1 + + bad, = np.where( + (matrix0[:, i, j] < min0) + | (matrix0[:, i, j] > max0) + | (matrix1[:, i, j] < min1) + | (matrix1[:, i, j] > max1) + ) + + if len(bad) > 0: + for detIdx in bad: + isBad[detIdx, i, j] = True + return pipeBase.Struct( + median0=median0, + median1=median1, + stdev0=stdev0, + stdev1=stdev1, + isBad=isBad + ) + + def replace_outliers(self, matrix0, matrix1, isBad, median0, median1): + """Do checks to see if an element of the matrix is out-of-family. + + Parameters + ---------- + matrix0 : `np.array`, (Ndet, Namp, Namp) + Matrix holding the 0th-order terms. + matrix1 : `np.array`, (Ndet, Namp, Namp) + Matrix holding the 1st-order terms. + isBad : `np.array`, (Ndet, Namp, Namp) + Matrix holding the boolean "is bad". + median0 : `np.array`, (Namp, Namp) + Matrix of median 0th-order terms. + median1 : `np.array`, (Namp, Namp) + Matrix of median 1st-order terms. + + Returns + ------- + results : `lsst.pipe.base.Struct` + Results struct containing + + ``new_matrix0`` + Replacement matrix0, with median substitutions. + (`np.array` (Ndet, Namp, Namp)). + ``new_matrix1`` + Replacement matrix1, with median substitutions. + (`np.array` (Ndet, Namp, Namp)). + + Raises + ------ + ValueError : + Raised if the inputs have a mismatch in size. + """ + if matrix0.shape != matrix1.shape: + raise ValueError("Shape disagreement!") + if matrix0.shape != isBad.shape: + raise ValueError("Shape disagreement!") + if median0.shape != median1.shape: + raise ValueError("Shape disagreement!") + + out0 = np.full_like(matrix0, 0.0) + out1 = np.full_like(matrix1, 0.0) + + for detIdx in range(matrix0.shape[0]): + for srcIdx in range(matrix0.shape[1]): + for tgtIdx in range(matrix0.shape[2]): + if isBad[detIdx, srcIdx, tgtIdx]: + self.log.info(f"Setting {detIdx} {srcIdx} {tgtIdx} from " + f"{matrix0[detIdx, srcIdx, tgtIdx]} to " + f"{median0[srcIdx, tgtIdx]} and " + f"{matrix1[detIdx, srcIdx, tgtIdx]} to " + f"{median1[srcIdx, tgtIdx]}") + out0[detIdx, srcIdx, tgtIdx] = median0[srcIdx, tgtIdx] + out1[detIdx, srcIdx, tgtIdx] = median1[srcIdx, tgtIdx] + else: + out0[detIdx, srcIdx, tgtIdx] = matrix0[detIdx, srcIdx, tgtIdx] + out1[detIdx, srcIdx, tgtIdx] = matrix1[detIdx, srcIdx, tgtIdx] + return pipeBase.Struct( + new_matrix0=out0, + new_matrix1=out1, + ) From c3eddb5aeec4b69230574651585ab0b2538a3baf Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 21 Nov 2025 14:32:39 -0800 Subject: [PATCH 4/8] Add code to set the crosstalk units correctly. --- pipelines/_ingredients/cpCrosstalkLSST.yaml | 4 ++++ python/lsst/cp/pipe/cpCrosstalk.py | 23 ++++++++++++++++++--- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/pipelines/_ingredients/cpCrosstalkLSST.yaml b/pipelines/_ingredients/cpCrosstalkLSST.yaml index e1f977fd6..ced7f7810 100644 --- a/pipelines/_ingredients/cpCrosstalkLSST.yaml +++ b/pipelines/_ingredients/cpCrosstalkLSST.yaml @@ -22,11 +22,13 @@ tasks: connections.inputRatios: "cpCrosstalkRatio" connections.outputCrosstalk: "crosstalkProposal" rejectNegativeSolutions: false + unitsAreElectrons: true cpCrosstalkFilter: class: lsst.cp.pipe.CrosstalkFilterTask config: connections.inputCrosstalk: "crosstalkProposal" connections.outputCrosstalk: "crosstalk" + unitsAreElectrons: true subsets: cpCrosstalk: subset: @@ -42,3 +44,5 @@ contracts: - cpCrosstalkExtract.connections.inputExp == cpCrosstalkIsr.connections.outputExposure - cpCrosstalkSolve.connections.inputRatios == cpCrosstalkExtract.connections.outputRatios - cpCrosstalkSolve.connections.outputCrosstalk == cpCrosstalkFilter.connections.inputCrosstalk + - cpCrosstalkIsr.config.doApplyGains == cpCrosstalkSolve.config.unitsAreElectrons + - cpCrosstalkSolve.config.unitsAreElectrons == cpCrosstalkFilter.config.unitsAreElectrons diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index 590e1c3a2..a0aae6132 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -511,6 +511,12 @@ class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig, doc="Filter generated crosstalk to remove marginal measurements?", ) + unitsAreElectrons = Field( + dtype=bool, + default=True, + doc="Crosstalk measurements have been done in electrons.", + ) + class CrosstalkSolveTask(pipeBase.PipelineTask): """Task to solve crosstalk from pixel ratios. @@ -667,9 +673,14 @@ def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, output # Populate the remainder of the calibration information. calib.hasCrosstalk = True calib.interChip = {} - - calib.updateMetadata(camera=camera, detector=calibDetector) - calib.updateMetadata(setCalibId=True, setDate=True) + calib.crosstalkRatiosUnits = 'electron' if self.config.unitsAreElectrons else 'adu' + calib.updateMetadata( + camera=camera, + detector=calibDetector, + setCalibId=True, + setCalibInfo=True, + setDate=True, + ) # Make an IsrProvenance(). provenance = IsrProvenance(calibType="CROSSTALK") @@ -944,6 +955,11 @@ class CrosstalkFilterConfig(pipeBase.PipelineTaskConfig, default=6.0, doc="Squared-term coefficient outlier clipping significance.", ) + unitsAreElectrons = Field( + dtype=bool, + default=True, + doc="Crosstalk measurements have been done in electrons.", + ) class CrosstalkFilterTask(pipeBase.PipelineTask): @@ -1070,6 +1086,7 @@ def run(self, inputCrosstalk, camera, inputDims, outputDims): outputCT.coeffs = e2v_final.new_matrix0[e2v_index] outputCT.coeffsSqr = e2v_final.new_matrix1[e2v_index] + outputCT.crosstalkRatiosUnits = 'electron' if self.config.unitsAreElectrons else 'adu' outputCT.updateMetadata( camera=camera, detector=camera[detId], From c98d0ddc682e344736cf1e5a898519e6bf7359a0 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Tue, 13 Jan 2026 16:44:43 -0800 Subject: [PATCH 5/8] Fix unit test error and underlying length bug. --- python/lsst/cp/pipe/cpCrosstalk.py | 7 ++++++- tests/test_measureCrosstalk.py | 27 +++++++++++++++++++++++---- 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index a0aae6132..b93dd502f 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -742,8 +742,13 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig values = np.array(ratios[ordering[tt]][ordering[ss]]) good_values = np.abs(values) < 1.0 # Discard unreasonable values values = values[good_values] + myfluxes = np.array(fluxes[ordering[tt]][ordering[ss]]) - myfluxes = myfluxes[good_values] + if len(myfluxes) != 0: + myfluxes = myfluxes[good_values] + else: + myfluxes = np.ones_like(values) + if len(values) != len(myfluxes): self.log.warning(f"Flux and ratio length disagree after first filter: {len(values)} {len(myfluxes)}") # noqa E501 diff --git a/tests/test_measureCrosstalk.py b/tests/test_measureCrosstalk.py index 88dbb3cca..021be693c 100644 --- a/tests/test_measureCrosstalk.py +++ b/tests/test_measureCrosstalk.py @@ -29,6 +29,19 @@ import lsst.ip.isr.isrMock as isrMock +class PretendRef: + "A class to act as a mock butler reference" + + def __init__(self, exposure): + self.exp = exposure + + def get(self, component=None): + if component: + raise RuntimeError("Components not supported in CT PretendRef.") + else: + return self.exp + + class MeasureCrosstalkTaskCases(lsst.utils.tests.TestCase): def setup_measureCrosstalk(self, isTrimmed=False, nSources=8): @@ -68,7 +81,7 @@ def setup_measureCrosstalk(self, isTrimmed=False, nSources=8): ctexConfig.threshold = 4000 ctexConfig.isTrimmed = isTrimmed ctex = CrosstalkExtractTask(config=ctexConfig) - fullResult = [] + measuredRatios = [] mockTask.config.isTrimmed = isTrimmed # Generate simulated set of exposures. @@ -84,14 +97,20 @@ def setup_measureCrosstalk(self, isTrimmed=False, nSources=8): exposure = mockTask.run() result = ctex.run(exposure) - fullResult.append(result.outputRatios) + measuredRatios.append(PretendRef(result.outputRatios)) + camera = mockTask.getCamera() # This test needs fixing: CZW: - return [True, True, True] + # return [True, True, True] # Generate the final measured CT ratios, uncertainties, pixel counts. + outputDims = { + 'detector': list(measuredRatios[0].get().keys())[0], + 'instrument': camera.getName(), + } + ctsConfig = CrosstalkSolveConfig() cts = CrosstalkSolveTask(config=ctsConfig) - finalResult = cts.run(fullResult) + finalResult = cts.run(inputRatios=measuredRatios, camera=camera, outputDims=outputDims) calib = finalResult.outputCrosstalk # Needed because measureCrosstalk cannot find coefficients equal to 0.0 From a6185dd934fe33a241046e9cc9484abbc06b7654 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 23 Jan 2026 12:22:57 -0800 Subject: [PATCH 6/8] Cleanup and docstring fixes. --- python/lsst/cp/pipe/cpCrosstalk.py | 55 ++++++++++++++++-------------- 1 file changed, 29 insertions(+), 26 deletions(-) diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index b93dd502f..a797a57a9 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -221,7 +221,7 @@ def run(self, inputExp, sourceExps=[]): if np.sum(mask) == 0: continue - newMask = np.where(np.bitwise_and(mask, detected), maskBit, 0) + newMask = np.where(mask & detected, maskBit, 0) for ampToMask in targetDetector: if ampName == ampToMask.getName(): # The amp we're considering already @@ -464,9 +464,6 @@ class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections, def __init__(self, *, config=None): super().__init__(config=config) - # if config.fluxOrder == 0 and False: - # self.inputs.discard("inputFluxes") - class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CrosstalkSolveConnections): @@ -739,18 +736,20 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig else: # ratios is ratios[Target][Source] # use tt for Target, use ss for Source, to match ip_isr. - values = np.array(ratios[ordering[tt]][ordering[ss]]) + values = np.asarray(ratios[ordering[tt]][ordering[ss]]) good_values = np.abs(values) < 1.0 # Discard unreasonable values values = values[good_values] - myfluxes = np.array(fluxes[ordering[tt]][ordering[ss]]) + myfluxes = np.asarray(fluxes[ordering[tt]][ordering[ss]]) if len(myfluxes) != 0: myfluxes = myfluxes[good_values] else: myfluxes = np.ones_like(values) if len(values) != len(myfluxes): - self.log.warning(f"Flux and ratio length disagree after first filter: {len(values)} {len(myfluxes)}") # noqa E501 + self.log.warning( + f"Flux and ratio length disagree after first filter: {len(values)} {len(myfluxes)}" + ) # Sigma clip using the inter-quartile distance and a # normal distribution. @@ -766,7 +765,9 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig values = values[good] myfluxes = myfluxes[good] if len(values) != len(myfluxes): - self.log.warning(f"Flux and ratio length disagree after second filter: {len(values)} {len(myfluxes)}") # noqa E501 + self.log.warning( + f"Flux and ratio length disagree after second filter: {len(values)} {len(myfluxes)}" + ) # Crosstalk calib is property[Source][Target]. calib.coeffNum[ss][tt] = len(values) @@ -1055,10 +1056,10 @@ def run(self, inputCrosstalk, camera, inputDims, outputDims): # filter those. pass - itl_c0 = np.array(itl_c0) - itl_c1 = np.array(itl_c1) - e2v_c0 = np.array(e2v_c0) - e2v_c1 = np.array(e2v_c1) + itl_c0 = np.asarray(itl_c0) + itl_c1 = np.asarray(itl_c1) + e2v_c0 = np.asarray(e2v_c0) + e2v_c1 = np.asarray(e2v_c1) itl_outliers = self.find_outliers(itl_c0, itl_c1) e2v_outliers = self.find_outliers(e2v_c0, e2v_c1) @@ -1109,9 +1110,9 @@ def find_outliers(self, matrix0, matrix1): Parameters ---------- - matrix0 : `np.array`, (Ndet, Namp, Namp) + matrix0 : `np.ndarray`, (Ndet, Namp, Namp) Matrix holding the 0th-order terms. - matrix1 : `np.array`, (Ndet, Namp, Namp) + matrix1 : `np.ndarray`, (Ndet, Namp, Namp) Matrix holding the 1st-order terms. Returns @@ -1120,16 +1121,18 @@ def find_outliers(self, matrix0, matrix1): Results struct containing ``median0`` - Median in-family value (`np.array` (Namp, Namp)). + Median in-family value (`np.ndarray` (Namp, Namp)). ``stdev0`` - MAD effective sigma in-family value (`np.array` (Namp, Namp)). + MAD effective sigma in-family value (`np.ndarray` + (Namp, Namp)). ``median1`` - Median in-family value (`np.array` (Namp, Namp)). + Median in-family value (`np.ndarray` (Namp, Namp)). ``stdev1`` - MAD effective sigma in-family value (`np.array` (Namp, Namp)). + MAD effective sigma in-family value (`np.ndarray` + (Namp, Namp)). ``isBad`` Boolean indicator that an element has been replaced - (`np.array` (Ndet, Namp, Namp)). + (`np.ndarray` (Ndet, Namp, Namp)). Raises ------ @@ -1189,15 +1192,15 @@ def replace_outliers(self, matrix0, matrix1, isBad, median0, median1): Parameters ---------- - matrix0 : `np.array`, (Ndet, Namp, Namp) + matrix0 : `np.ndarray`, (Ndet, Namp, Namp) Matrix holding the 0th-order terms. - matrix1 : `np.array`, (Ndet, Namp, Namp) + matrix1 : `np.ndarray`, (Ndet, Namp, Namp) Matrix holding the 1st-order terms. - isBad : `np.array`, (Ndet, Namp, Namp) + isBad : `np.ndarray`, (Ndet, Namp, Namp) Matrix holding the boolean "is bad". - median0 : `np.array`, (Namp, Namp) + median0 : `np.ndarray`, (Namp, Namp) Matrix of median 0th-order terms. - median1 : `np.array`, (Namp, Namp) + median1 : `np.ndarray`, (Namp, Namp) Matrix of median 1st-order terms. Returns @@ -1207,10 +1210,10 @@ def replace_outliers(self, matrix0, matrix1, isBad, median0, median1): ``new_matrix0`` Replacement matrix0, with median substitutions. - (`np.array` (Ndet, Namp, Namp)). + (`np.ndarray` (Ndet, Namp, Namp)). ``new_matrix1`` Replacement matrix1, with median substitutions. - (`np.array` (Ndet, Namp, Namp)). + (`np.ndarray` (Ndet, Namp, Namp)). Raises ------ From 0e9c156c7d8bbb053676d1e632685fe485e775c8 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Thu, 12 Mar 2026 15:06:21 -0700 Subject: [PATCH 7/8] Clarify names and comments. --- pipelines/LSSTCam/cpCrosstalk.yaml | 2 -- python/lsst/cp/pipe/cpCrosstalk.py | 17 +++++++++++------ 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/pipelines/LSSTCam/cpCrosstalk.yaml b/pipelines/LSSTCam/cpCrosstalk.yaml index 9768be939..cf4498fc8 100644 --- a/pipelines/LSSTCam/cpCrosstalk.yaml +++ b/pipelines/LSSTCam/cpCrosstalk.yaml @@ -9,8 +9,6 @@ tasks: doDefect: true cpCrosstalkExtract: class: lsst.cp.pipe.CrosstalkExtractTask - config: - growMaskRadius: 20 cpCrosstalkSolve: class: lsst.cp.pipe.CrosstalkSolveTask config: diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index a797a57a9..8a61f0c2f 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -118,7 +118,7 @@ class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig, ) growMaskRadius = Field( dtype=int, - default=0, + default=20, doc="Radius to grow CT_TEMP masks prior to background estimation." ) @@ -229,7 +229,7 @@ def run(self, inputExp, sourceExps=[]): extractedAmp = inputExp[ampToMask.getBBox()] # The mask needs to be flipped to match the target - flippedMask = self._flipMask(newMask, amp, ampToMask) + flippedMask = self._flipMaskArray(newMask, amp, ampToMask) extractedAmp.mask.array[:, :] |= flippedMask # Optionally dilate these masks by some amount: @@ -315,7 +315,7 @@ def run(self, inputExp, sourceExps=[]): outputFluxes=ddict2dict(outputFluxes) ) - def _flipMask(self, maskArray, sourceAmp, targetAmp): + def _flipMaskArray(self, maskArray, sourceAmp, targetAmp): """Flip an array from a sourceAmp to match the readout order of targetAmp. @@ -483,7 +483,7 @@ class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig, fluxOrder = Field( dtype=int, default=0, - doc="Order of source flux fit to crosstalk. 0=simple linear; 1=first order non-linear.", + doc="Order of source flux fit to crosstalk ratios. 0=simple linear; 1=first order non-linear.", ) rejectNegativeSolutions = Field( @@ -722,6 +722,11 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig ------- calib : `lsst.ip.isr.CrosstalkCalib` The output crosstalk calibration. + + Raises + ------ + RuntimeError + Raised if filtering values and fluxes results in length mismatch. """ calib = CrosstalkCalib(nAmp=len(ratios)) @@ -747,7 +752,7 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig myfluxes = np.ones_like(values) if len(values) != len(myfluxes): - self.log.warning( + raise RuntimeError( f"Flux and ratio length disagree after first filter: {len(values)} {len(myfluxes)}" ) @@ -765,7 +770,7 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig values = values[good] myfluxes = myfluxes[good] if len(values) != len(myfluxes): - self.log.warning( + raise RuntimeError( f"Flux and ratio length disagree after second filter: {len(values)} {len(myfluxes)}" ) From 2a06b403b82950613f78fc1437028412856f0bb4 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Wed, 18 Mar 2026 14:45:17 -0700 Subject: [PATCH 8/8] Clarify comment. --- python/lsst/cp/pipe/cpCrosstalk.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/python/lsst/cp/pipe/cpCrosstalk.py b/python/lsst/cp/pipe/cpCrosstalk.py index 8a61f0c2f..6931eba9f 100644 --- a/python/lsst/cp/pipe/cpCrosstalk.py +++ b/python/lsst/cp/pipe/cpCrosstalk.py @@ -742,7 +742,12 @@ def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSig # ratios is ratios[Target][Source] # use tt for Target, use ss for Source, to match ip_isr. values = np.asarray(ratios[ordering[tt]][ordering[ss]]) - good_values = np.abs(values) < 1.0 # Discard unreasonable values + # Discard unreasonable ratios: anything with + # abs(value) > 1.0 indicates that the target location + # is brighter than the source, so these target pixels + # have some additional contaminating flux (another + # spot, diffraction spike, etc). + good_values = np.abs(values) < 1.0 values = values[good_values] myfluxes = np.asarray(fluxes[ordering[tt]][ordering[ss]])