Coverage for python / lsst / cp / pipe / cpCrosstalk.py: 13%
473 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-23 08:59 +0000
1# This file is part of cp_pipe
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <http://www.gnu.org/licenses/>.
21import itertools
22import numpy as np
24from collections import defaultdict
25from copy import copy
26from scipy.stats import median_abs_deviation
28import lsst.pipe.base as pipeBase
29import lsst.pipe.base.connectionTypes as cT
31from lsstDebug import getDebugFrame
32from lsst.afw.cameraGeom import ReadoutCorner
33from lsst.afw.detection import FootprintSet, Threshold
34from lsst.afw.display import getDisplay
35from lsst.pex.config import ConfigurableField, Field, ListField
36from lsst.ip.isr import CrosstalkCalib, IsrProvenance, growMasks
37from lsst.cp.pipe.utils import (ddict2dict, sigmaClipCorrection)
38from lsst.meas.algorithms import SubtractBackgroundTask
40__all__ = ["CrosstalkExtractConfig", "CrosstalkExtractTask",
41 "CrosstalkSolveTask", "CrosstalkSolveConfig",
42 "CrosstalkFilterTask", "CrosstalkFilterConfig"]
45class CrosstalkExtractConnections(pipeBase.PipelineTaskConnections,
46 dimensions=("instrument", "exposure", "detector")):
47 inputExp = cT.Input(
48 name="crosstalkInputs",
49 doc="Input post-ISR processed exposure to measure crosstalk from.",
50 storageClass="Exposure",
51 dimensions=("instrument", "exposure", "detector"),
52 multiple=False,
53 )
54 # TODO: Depends on DM-21904.
55 sourceExp = cT.Input(
56 name="crosstalkSource",
57 doc="Post-ISR exposure to measure for inter-chip crosstalk onto inputExp.",
58 storageClass="Exposure",
59 dimensions=("instrument", "exposure", "detector"),
60 multiple=True,
61 deferLoad=True,
62 # lookupFunction=None,
63 )
65 outputRatios = cT.Output(
66 name="crosstalkRatios",
67 doc="Extracted crosstalk pixel ratios.",
68 storageClass="StructuredDataDict",
69 dimensions=("instrument", "exposure", "detector"),
70 )
71 outputFluxes = cT.Output(
72 name="crosstalkFluxes",
73 doc="Source pixel fluxes used in ratios.",
74 storageClass="StructuredDataDict",
75 dimensions=("instrument", "exposure", "detector"),
76 )
78 def __init__(self, *, config=None):
79 super().__init__(config=config)
80 # Discard sourceExp until DM-21904 allows full interchip
81 # measurements.
82 self.inputs.discard("sourceExp")
85class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig,
86 pipelineConnections=CrosstalkExtractConnections):
87 """Configuration for the measurement of pixel ratios.
88 """
90 doMeasureInterchip = Field(
91 dtype=bool,
92 default=False,
93 doc="Measure inter-chip crosstalk as well?",
94 )
95 threshold = Field(
96 dtype=float,
97 default=30000,
98 doc="Minimum level of source pixels for which to measure crosstalk."
99 )
100 ignoreSaturatedPixels = Field(
101 dtype=bool,
102 default=False,
103 doc="Should saturated pixels be ignored?"
104 )
105 badMask = ListField(
106 dtype=str,
107 default=["BAD", "INTRP"],
108 doc="Mask planes to ignore when identifying source pixels."
109 )
110 isTrimmed = Field(
111 dtype=bool,
112 default=True,
113 doc="Is the input exposure trimmed?"
114 )
115 background = ConfigurableField(
116 target=SubtractBackgroundTask,
117 doc="Background estimation task.",
118 )
119 growMaskRadius = Field(
120 dtype=int,
121 default=20,
122 doc="Radius to grow CT_TEMP masks prior to background estimation."
123 )
125 def setDefaults(self):
126 super().setDefaults()
127 # is this really the best way to do this??
128 self.background.useApprox = False
129 self.background.ignoredPixelMask.append("CT_TEMP")
131 def validate(self):
132 super().validate()
134 # Ensure the handling of the SAT mask plane is consistent
135 # with the ignoreSaturatedPixels value.
136 if self.ignoreSaturatedPixels:
137 if 'SAT' not in self.badMask:
138 self.badMask.append('SAT')
139 else:
140 if 'SAT' in self.badMask:
141 self.badMask = [mask for mask in self.badMask if mask != 'SAT']
144class CrosstalkExtractTask(pipeBase.PipelineTask):
145 """Task to measure pixel ratios to find crosstalk.
146 """
148 ConfigClass = CrosstalkExtractConfig
149 _DefaultName = 'cpCrosstalkExtract'
151 def __init__(self, **kwargs):
152 super().__init__(**kwargs)
153 self.makeSubtask("background")
155 def run(self, inputExp, sourceExps=[]):
156 """Measure pixel ratios between amplifiers in inputExp.
158 Extract crosstalk ratios between different amplifiers.
160 For pixels above ``config.threshold``, we calculate the ratio
161 between each background-subtracted target amp and the source
162 amp. We return a list of ratios for each pixel for each
163 target/source combination, as nested dictionary containing the
164 ratio.
166 Parameters
167 ----------
168 inputExp : `lsst.afw.image.Exposure`
169 Input exposure to measure pixel ratios on.
170 sourceExp : `list` [`lsst.afw.image.Exposure`], optional
171 List of chips to use as sources to measure inter-chip
172 crosstalk.
174 Returns
175 -------
176 results : `lsst.pipe.base.Struct`
177 The results struct containing:
179 ``outputRatios``
180 A catalog of ratio lists. The dictionaries are
181 indexed such that:
182 outputRatios[targetChip][sourceChip][targetAmp][sourceAmp]
183 contains the ratio list for that combination (`dict`
184 [`dict` [`dict` [`dict` [`list`]]]]).
185 ``outputFluxes``
186 A catalog of flux lists. The dictionaries are
187 indexed such that:
188 outputFluxes[sourceChip][sourceAmp] contains the flux
189 list used in the outputRatios (`dict` [`dict`
190 [`list`]]).
191 """
192 outputRatios = defaultdict(lambda: defaultdict(dict))
193 outputFluxes = defaultdict(lambda: defaultdict(dict))
195 threshold = self.config.threshold
196 badPixels = list(self.config.badMask)
198 targetDetector = inputExp.getDetector()
199 targetChip = targetDetector.getName()
201 # Always look at the target chip first, then go to any other
202 # supplied exposures.
203 sourceExtractExps = [inputExp]
204 sourceExtractExps.extend(sourceExps)
206 self.log.info("Measuring full detector background for target: %s", targetChip)
207 targetIm = inputExp.getMaskedImage()
209 FootprintSet(targetIm, Threshold(threshold), "DETECTED")
210 detected = targetIm.getMask().getPlaneBitMask("DETECTED")
211 targetIm.getMask().addMaskPlane("CT_TEMP")
212 maskBit = targetIm.getMask().getPlaneBitMask("CT_TEMP")
213 bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"])
215 # Carry over over-threshold masked pixels to other amplifiers.
216 for amp in targetDetector:
217 ampIm = inputExp[amp.getBBox()]
218 ampName = amp.getName()
220 mask = ampIm.mask.array & detected
221 if np.sum(mask) == 0:
222 continue
224 newMask = np.where(mask & detected, maskBit, 0)
225 for ampToMask in targetDetector:
226 if ampName == ampToMask.getName():
227 # The amp we're considering already
228 continue
230 extractedAmp = inputExp[ampToMask.getBBox()]
231 # The mask needs to be flipped to match the target
232 flippedMask = self._flipMaskArray(newMask, amp, ampToMask)
234 extractedAmp.mask.array[:, :] |= flippedMask
235 # Optionally dilate these masks by some amount:
236 growMasks(inputExp.mask, radius=self.config.growMaskRadius,
237 maskNameList=["CT_TEMP"], maskValue="CT_TEMP")
239 # We've now masked the source pixels, and any potential CT
240 # pixels, so this should be just the
241 # background/reflections/etc.
242 backgroundModel = self.background.fitBackground(inputExp.maskedImage)
243 backgroundIm = backgroundModel.getImageF()
244 self.debugView('extract', inputExp)
246 # Begin search for bright pixels, and their associated crosstalk
247 # signals.
248 for sourceExp in sourceExtractExps:
249 # This loop exists to support future inter-chip searches.
250 sourceDetector = sourceExp.getDetector()
251 sourceChip = sourceDetector.getName()
252 sourceIm = sourceExp.getMaskedImage()
253 bad = sourceIm.getMask().getPlaneBitMask(badPixels)
254 self.log.info("Measuring crosstalk from source: %s", sourceChip)
256 if sourceExp != inputExp:
257 FootprintSet(sourceIm, Threshold(threshold), "DETECTED")
258 detected = sourceIm.getMask().getPlaneBitMask("DETECTED")
260 # The dictionary of amp-to-amp ratios for this pair of
261 # source->target detectors.
262 ratioDict = defaultdict(lambda: defaultdict(np.array))
263 extractedCount = 0
265 for sourceAmp in sourceDetector:
266 sourceAmpName = sourceAmp.getName()
267 sourceAmpBBox = sourceAmp.getBBox() if self.config.isTrimmed else sourceAmp.getRawDataBBox()
268 sourceAmpImage = sourceIm[sourceAmpBBox]
269 sourceMask = sourceAmpImage.mask.array
270 select = ((sourceMask & detected > 0)
271 & (sourceMask & bad == 0)
272 & np.isfinite(sourceAmpImage.image.array))
273 count = np.sum(select)
274 self.log.debug(" Source amplifier: %s", sourceAmpName)
276 outputFluxes[sourceChip][sourceAmpName] = sourceAmpImage.image.array[select].tolist()
278 for targetAmp in targetDetector:
279 # iterate over targetExposure
280 targetAmpName = targetAmp.getName()
281 if sourceAmpName == targetAmpName and sourceChip == targetChip:
282 ratioDict[targetAmpName][sourceAmpName] = []
283 continue
285 self.log.debug(" Target amplifier: %s", targetAmpName)
287 targetAmpImage = CrosstalkCalib.extractAmp(targetIm,
288 targetAmp, sourceAmp,
289 isTrimmed=self.config.isTrimmed)
290 targetBkgImage = CrosstalkCalib.extractAmp(backgroundIm,
291 targetAmp, sourceAmp,
292 isTrimmed=self.config.isTrimmed)
294 bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"])
296 ratios = ((targetAmpImage.image.array[select] - targetBkgImage.array[select])
297 / sourceAmpImage.image.array[select])
299 ratioDict[targetAmpName][sourceAmpName] = ratios.tolist()
300 self.log.info("Amp extracted %d pixels from %s -> %s",
301 count, sourceAmpName, targetAmpName)
302 extractedCount += count
304 self.debugPixels('pixels',
305 sourceAmpImage.image.array[select],
306 targetAmpImage.image.array[select] - bg,
307 sourceAmpName, targetAmpName)
309 self.log.info("Extracted %d pixels from %s -> %s (targetBG: %f)",
310 extractedCount, sourceChip, targetChip, bg)
311 outputRatios[targetChip][sourceChip] = ratioDict
313 return pipeBase.Struct(
314 outputRatios=ddict2dict(outputRatios),
315 outputFluxes=ddict2dict(outputFluxes)
316 )
318 def _flipMaskArray(self, maskArray, sourceAmp, targetAmp):
319 """Flip an array from a sourceAmp to match the readout order of
320 targetAmp.
322 Parameters
323 ----------
324 maskArray : `np.ndarray`
325 Mask data to flip.
326 sourceAmp : `lsst.afw.cameraGeom.Amplifier`
327 Amplifier corresponding to the maskArray.
328 targetAmp : `lsst.afw.cameraGeom.Amplifier`
329 Amplifier corresponding to the output mask.
331 Returns
332 -------
333 maskFlipped : `np.ndarray`
334 The flipped mask.
336 See Also
337 -----
338 lsst.ip.isr.CrosstalkCalib.extractAmp()
339 """
340 maskFlipped = copy(maskArray)
341 if sourceAmp.getReadoutCorner() == targetAmp.getReadoutCorner():
342 return maskFlipped
344 X_FLIP = {ReadoutCorner.LL: False,
345 ReadoutCorner.LR: True,
346 ReadoutCorner.UL: False,
347 ReadoutCorner.UR: True}
348 Y_FLIP = {ReadoutCorner.LL: False,
349 ReadoutCorner.LR: False,
350 ReadoutCorner.UL: True,
351 ReadoutCorner.UR: True}
353 sourceAmpCorner = sourceAmp.getReadoutCorner()
354 targetAmpCorner = targetAmp.getReadoutCorner()
356 # Flipping is necessary only if the desired configuration doesn't match
357 # what we currently have.
358 xFlip = X_FLIP[targetAmpCorner] ^ X_FLIP[sourceAmpCorner]
359 yFlip = Y_FLIP[targetAmpCorner] ^ Y_FLIP[sourceAmpCorner]
361 if xFlip:
362 maskFlipped = np.fliplr(maskFlipped)
363 if yFlip:
364 maskFlipped = np.flipud(maskFlipped)
366 return maskFlipped
368 def debugView(self, stepname, exposure):
370 """Utility function to examine the image being processed.
372 Parameters
373 ----------
374 stepname : `str`
375 State of processing to view.
376 exposure : `lsst.afw.image.Exposure`
377 Exposure to view.
378 """
379 frame = getDebugFrame(self._display, stepname)
380 if frame:
381 display = getDisplay(frame)
382 display.scale('asinh', 'zscale')
383 display.mtv(exposure)
385 prompt = "Press Enter to continue: "
386 while True:
387 ans = input(prompt).lower()
388 if ans in ("", "c",):
389 break
391 def debugPixels(self, stepname, pixelsIn, pixelsOut, sourceName, targetName):
392 """Utility function to examine the CT ratio pixel values.
394 Parameters
395 ----------
396 stepname : `str`
397 State of processing to view.
398 pixelsIn : `np.ndarray`, (N,)
399 Pixel values from the potential crosstalk source.
400 pixelsOut : `np.ndarray`, (N,)
401 Pixel values from the potential crosstalk target.
402 sourceName : `str`
403 Source amplifier name
404 targetName : `str`
405 Target amplifier name
406 """
407 frame = getDebugFrame(self._display, stepname)
408 if frame:
409 import matplotlib.pyplot as plt
410 figure = plt.figure(1)
411 figure.clear()
413 axes = figure.add_axes((0.1, 0.1, 0.8, 0.8))
414 axes.plot(pixelsIn, pixelsOut / pixelsIn, 'k+')
415 plt.xlabel("Source amplifier pixel value")
416 plt.ylabel("Measured pixel ratio")
417 plt.title(f"(Source {sourceName} -> Target {targetName}) median ratio: "
418 f"{(np.median(pixelsOut / pixelsIn))}")
419 figure.show()
421 prompt = "Press Enter to continue: "
422 while True:
423 ans = input(prompt).lower()
424 if ans in ("", "c",):
425 break
426 plt.close()
429class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections,
430 dimensions=("instrument", "detector")):
431 inputRatios = cT.Input(
432 name="crosstalkRatios",
433 doc="Ratios measured for an input exposure.",
434 storageClass="StructuredDataDict",
435 dimensions=("instrument", "exposure", "detector"),
436 multiple=True,
437 deferLoad=True,
438 )
439 inputFluxes = cT.Input(
440 name="crosstalkFluxes",
441 doc="Fluxes of CT source pixels, for nonlinear fits.",
442 storageClass="StructuredDataDict",
443 dimensions=("instrument", "exposure", "detector"),
444 multiple=True,
445 deferLoad=True,
446 )
447 camera = cT.PrerequisiteInput(
448 name="camera",
449 doc="Camera the input data comes from.",
450 storageClass="Camera",
451 dimensions=("instrument",),
452 isCalibration=True,
453 )
455 outputCrosstalk = cT.Output(
456 name="crosstalkProposal",
457 doc="Output proposed crosstalk calibration.",
458 storageClass="CrosstalkCalib",
459 dimensions=("instrument", "detector"),
460 multiple=False,
461 isCalibration=True,
462 )
464 def __init__(self, *, config=None):
465 super().__init__(config=config)
468class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig,
469 pipelineConnections=CrosstalkSolveConnections):
470 """Configuration for the solving of crosstalk from pixel ratios.
471 """
473 rejIter = Field(
474 dtype=int,
475 default=3,
476 doc="Number of rejection iterations for final coefficient calculation.",
477 )
478 rejSigma = Field(
479 dtype=float,
480 default=2.0,
481 doc="Rejection threshold (sigma) for final coefficient calculation.",
482 )
483 fluxOrder = Field(
484 dtype=int,
485 default=0,
486 doc="Order of source flux fit to crosstalk ratios. 0=simple linear; 1=first order non-linear.",
487 )
489 rejectNegativeSolutions = Field(
490 dtype=bool,
491 default=False,
492 doc="Should solutions with negative coefficients (which add flux to the target) be excluded?",
493 )
495 significanceLimit = Field(
496 dtype=float,
497 default=3.0,
498 doc="Sigma significance level to use in marking a coefficient valid.",
499 )
500 doSignificanceScaling = Field(
501 dtype=bool,
502 default=True,
503 doc="Scale error by 1/sqrt(N) in calculating significant coefficients?",
504 )
505 doFiltering = Field(
506 dtype=bool,
507 default=False,
508 doc="Filter generated crosstalk to remove marginal measurements?",
509 )
511 unitsAreElectrons = Field(
512 dtype=bool,
513 default=True,
514 doc="Crosstalk measurements have been done in electrons.",
515 )
518class CrosstalkSolveTask(pipeBase.PipelineTask):
519 """Task to solve crosstalk from pixel ratios.
520 """
522 ConfigClass = CrosstalkSolveConfig
523 _DefaultName = 'cpCrosstalkSolve'
525 def runQuantum(self, butlerQC, inputRefs, outputRefs):
526 """Ensure that the input and output dimensions are passed along.
528 Parameters
529 ----------
530 butlerQC : `lsst.daf.butler.QuantumContext`
531 Butler to operate on.
532 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
533 Input data refs to load.
534 ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection`
535 Output data refs to persist.
536 """
537 inputs = butlerQC.get(inputRefs)
539 # Use the dimensions to set calib/provenance information.
540 inputs['inputDims'] = [dict(exp.dataId.required) for exp in inputRefs.inputRatios]
541 inputs['outputDims'] = dict(outputRefs.outputCrosstalk.dataId.required)
543 outputs = self.run(**inputs)
544 butlerQC.put(outputs, outputRefs)
546 def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, outputDims=None):
547 """Combine ratios to produce crosstalk coefficients.
549 Parameters
550 ----------
551 inputRatios : `list` [`dict` [`dict` [`dict` [`dict` [`list`]]]]]
552 A list of nested dictionaries of ratios indexed by target
553 and source chip, then by target and source amplifier.
554 inputFluxes : `list` [`dict` [`dict` [`list`]]]
555 A list of nested dictionaries of source pixel fluxes, indexed
556 by source chip and amplifier.
557 camera : `lsst.afw.cameraGeom.Camera`
558 Input camera.
559 inputDims : `list` [`lsst.daf.butler.DataCoordinate`]
560 DataIds to use to construct provenance.
561 outputDims : `list` [`lsst.daf.butler.DataCoordinate`]
562 DataIds to use to populate the output calibration.
564 Returns
565 -------
566 results : `lsst.pipe.base.Struct`
567 The results struct containing:
569 ``outputCrosstalk``
570 Final crosstalk calibration
571 (`lsst.ip.isr.CrosstalkCalib`).
572 ``outputProvenance``
573 Provenance data for the new calibration
574 (`lsst.ip.isr.IsrProvenance`).
576 Raises
577 ------
578 RuntimeError
579 Raised if the input data contains multiple target detectors.
580 """
581 if outputDims:
582 calibChip = outputDims['detector']
583 instrument = outputDims['instrument']
584 else:
585 # calibChip needs to be set manually in Gen2.
586 calibChip = None
587 instrument = None
589 if camera and calibChip is not None:
590 calibDetector = camera[calibChip]
591 ordering = [amp.getName() for amp in calibDetector]
592 else:
593 calibDetector = None
594 ordering = None
596 self.log.info("Combining measurements from %d ratios and %d fluxes",
597 len(inputRatios), len(inputFluxes) if inputFluxes else 0)
599 if inputFluxes is None:
600 inputFluxes = [None for exp in inputRatios]
601 if inputDims is None:
602 inputDims = [{} for exp in inputRatios]
604 combinedRatios = defaultdict(lambda: defaultdict(list))
605 combinedFluxes = defaultdict(lambda: defaultdict(list))
607 for ratioRef, fluxRef, dimensions in zip(inputRatios, inputFluxes, inputDims):
608 ratioDict = ratioRef.get()
609 if fluxRef is not None:
610 fluxDict = fluxRef.get()
611 else:
612 fluxDict = None
613 for targetChip in ratioDict:
614 if calibChip and targetChip != calibChip and targetChip != calibDetector.getName():
615 raise RuntimeError(f"Target chip: {targetChip} does not match calibration dimension: "
616 f"{calibChip}, {calibDetector.getName()}!")
618 sourceChip = targetChip
619 if sourceChip in ratioDict[targetChip]:
620 ratios = ratioDict[targetChip][sourceChip]
622 for targetAmp in ratios:
623 for sourceAmp in ratios[targetAmp]:
624 if fluxDict:
625 if len(ratios[targetAmp][sourceAmp]) != len(fluxDict[sourceChip][sourceAmp]):
626 if targetAmp != sourceAmp:
627 # This usually triggers when
628 # targetAmp == sourceAmp. The
629 # sourceAmp has flux entries,
630 # but by definition cannot
631 # have ratios for itself. Only
632 # warn when this isn't the
633 # case.
634 self.log.warning(f"Length mismatch for ratios {len(ratios[targetAmp][sourceAmp])} " # noqa E501
635 f"and fluxes {len(fluxDict[sourceChip][sourceAmp])} for " # noqa E501
636 f"Source {sourceAmp} and target {targetAmp} "
637 f"Rejecting this {dimensions}")
638 continue
639 combinedFluxes[targetAmp][sourceAmp].extend(fluxDict[sourceChip][sourceAmp])
640 combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp])
641 # TODO: DM-21904
642 # Iterating over all other entries in
643 # ratioDict[targetChip] will yield inter-chip terms.
645 for targetAmp in combinedRatios:
646 for sourceAmp in combinedRatios[targetAmp]:
647 self.log.info("Read %d pixels for %s -> %s",
648 len(combinedRatios[targetAmp][sourceAmp]),
649 sourceAmp, targetAmp)
650 if len(combinedRatios[targetAmp][sourceAmp]) > 1:
651 self.debugRatios('reduce', combinedRatios, targetAmp, sourceAmp)
653 if self.config.fluxOrder < 2:
654 self.log.info("Fitting crosstalk coefficients with order {self.config.fluxOrder}")
656 calib = self.measureCrosstalkCoefficients(combinedRatios, ordering,
657 combinedFluxes,
658 self.config.rejIter, self.config.rejSigma)
659 else:
660 raise NotImplementedError("Higher order non-linear crosstalk terms are not yet supported.")
662 self.log.info("Number of valid coefficients: %d", np.sum(calib.coeffValid))
664 if self.config.doFiltering:
665 # This step will apply the calculated validity values to
666 # censor poorly measured coefficients.
667 self.log.info("Filtering measured crosstalk to remove invalid solutions.")
668 calib = self.filterCrosstalkCalib(calib)
670 # Populate the remainder of the calibration information.
671 calib.hasCrosstalk = True
672 calib.interChip = {}
673 calib.crosstalkRatiosUnits = 'electron' if self.config.unitsAreElectrons else 'adu'
674 calib.updateMetadata(
675 camera=camera,
676 detector=calibDetector,
677 setCalibId=True,
678 setCalibInfo=True,
679 setDate=True,
680 )
682 # Make an IsrProvenance().
683 provenance = IsrProvenance(calibType="CROSSTALK")
684 provenance._detectorName = calibChip
685 if inputDims:
686 provenance.fromDataIds(inputDims)
687 provenance._instrument = instrument
688 provenance.updateMetadata()
690 return pipeBase.Struct(
691 outputCrosstalk=calib,
692 outputProvenance=provenance,
693 )
695 def measureCrosstalkCoefficients(self, ratios, ordering, fluxes, rejIter, rejSigma):
696 """Measure crosstalk coefficients from the ratios.
698 Given a list of ratios for each target/source amp combination,
699 we measure a sigma clipped mean and error.
701 The coefficient errors returned are the standard deviation of
702 the final set of clipped input ratios.
704 Parameters
705 ----------
706 ratios : `dict` [`dict` [`numpy.ndarray`]]
707 Catalog of arrays of ratios. The ratio arrays are one-dimensional.
708 ordering : `list` [`str`] or None
709 List to use as a mapping between amplifier names (the
710 elements of the list) and their position in the output
711 calibration (the matching index of the list). If no
712 ordering is supplied, the order of the keys in the ratio
713 catalog is used.
714 fluxes : `dict` [`dict` [`numpy.ndarray`]]
715 Catalog of arrays of fluxes. The flux arrays are one-dimensional.
716 rejIter : `int`
717 Number of rejection iterations.
718 rejSigma : `float`
719 Rejection threshold (sigma).
721 Returns
722 -------
723 calib : `lsst.ip.isr.CrosstalkCalib`
724 The output crosstalk calibration.
726 Raises
727 ------
728 RuntimeError
729 Raised if filtering values and fluxes results in length mismatch.
730 """
731 calib = CrosstalkCalib(nAmp=len(ratios))
733 if ordering is None:
734 ordering = list(ratios.keys())
736 # Calibration stores coefficients as a numpy ndarray.
737 for ss, tt in itertools.product(range(calib.nAmp), range(calib.nAmp)):
738 if ss == tt:
739 values = []
740 myfluxes = []
741 else:
742 # ratios is ratios[Target][Source]
743 # use tt for Target, use ss for Source, to match ip_isr.
744 values = np.asarray(ratios[ordering[tt]][ordering[ss]])
745 # Discard unreasonable ratios: anything with
746 # abs(value) > 1.0 indicates that the target location
747 # is brighter than the source, so these target pixels
748 # have some additional contaminating flux (another
749 # spot, diffraction spike, etc).
750 good_values = np.abs(values) < 1.0
751 values = values[good_values]
753 myfluxes = np.asarray(fluxes[ordering[tt]][ordering[ss]])
754 if len(myfluxes) != 0:
755 myfluxes = myfluxes[good_values]
756 else:
757 myfluxes = np.ones_like(values)
759 if len(values) != len(myfluxes):
760 raise RuntimeError(
761 f"Flux and ratio length disagree after first filter: {len(values)} {len(myfluxes)}"
762 )
764 # Sigma clip using the inter-quartile distance and a
765 # normal distribution.
766 if ss != tt:
767 for rej in range(rejIter):
768 if len(values) == 0:
769 break
770 lo, med, hi = np.percentile(values, [25.0, 50.0, 75.0])
771 sigma = 0.741*(hi - lo)
772 good = np.abs(values - med) < rejSigma*sigma
773 if good.sum() == len(good) or good.sum() == 0:
774 break
775 values = values[good]
776 myfluxes = myfluxes[good]
777 if len(values) != len(myfluxes):
778 raise RuntimeError(
779 f"Flux and ratio length disagree after second filter: {len(values)} {len(myfluxes)}"
780 )
782 # Crosstalk calib is property[Source][Target].
783 calib.coeffNum[ss][tt] = len(values)
784 significanceThreshold = 0.0
785 if len(values) == 0:
786 self.log.warning("No values for matrix element %d,%d" % (ss, tt))
787 calib.coeffs[ss][tt] = np.nan
788 calib.coeffErr[ss][tt] = np.nan
789 calib.coeffValid[ss][tt] = False
790 polyfit = [0.0, 0.0]
791 else:
792 calib.coeffs[ss][tt] = np.mean(values)
793 polyfit = np.polyfit(myfluxes, values, 1)
795 if self.config.fluxOrder == 1:
796 # substitute polyfit solution.
797 calib.coeffs[ss][tt] = polyfit[1]
798 calib.coeffsSqr[ss][tt] = polyfit[0]
800 if self.config.rejectNegativeSolutions and calib.coeffs[ss][tt] < 0.0:
801 calib.coeffs[ss][tt] = 0.0
803 if calib.coeffNum[ss][tt] <= 1:
804 calib.coeffErr[ss][tt] = np.nan
805 calib.coeffSqr[ss][tt] = np.nan
806 calib.coeffValid[ss][tt] = False
807 else:
808 correctionFactor = sigmaClipCorrection(rejSigma)
809 calib.coeffErr[ss][tt] = np.std(values) * correctionFactor
811 # Use sample stdev.
812 significanceThreshold = self.config.significanceLimit * calib.coeffErr[ss][tt]
813 if self.config.doSignificanceScaling is True:
814 # Enabling this calculates the stdev of the mean.
815 significanceThreshold /= np.sqrt(calib.coeffNum[ss][tt])
816 calib.coeffValid[ss][tt] = np.abs(calib.coeffs[ss][tt]) > significanceThreshold
817 self.debugRatios('measure', ratios, ordering[ss], ordering[tt],
818 calib.coeffs[ss][tt], calib.coeffValid[ss][tt])
819 self.log.info("Measured %s -> %s Coeff: %e Err: %e N: %d Valid: %s Limit: %e Quadratic: %s",
820 ordering[ss], ordering[tt], calib.coeffs[ss][tt], calib.coeffErr[ss][tt],
821 calib.coeffNum[ss][tt], calib.coeffValid[ss][tt], significanceThreshold,
822 polyfit)
824 return calib
826 @staticmethod
827 def filterCrosstalkCalib(inCalib):
828 """Apply valid constraints to the measured values.
830 Any measured coefficient that is determined to be invalid is
831 set to zero, and has the error set to nan. The validation is
832 determined by checking that the measured coefficient is larger
833 than the calculated standard error of the mean.
835 Parameters
836 ----------
837 inCalib : `lsst.ip.isr.CrosstalkCalib`
838 Input calibration to filter.
840 Returns
841 -------
842 outCalib : `lsst.ip.isr.CrosstalkCalib`
843 Filtered calibration.
844 """
845 outCalib = CrosstalkCalib()
846 outCalib.nAmp = inCalib.nAmp
848 outCalib.coeffs = inCalib.coeffs
849 outCalib.coeffs[~inCalib.coeffValid] = 0.0
851 outCalib.coeffErr = inCalib.coeffErr
852 outCalib.coeffErr[~inCalib.coeffValid] = np.nan
854 outCalib.coeffNum = inCalib.coeffNum
855 outCalib.coeffValid = inCalib.coeffValid
857 outCalib.coeffsSqr = inCalib.coeffsSqr
858 outCalib.coeffsSqr[~inCalib.coeffValid] = 0.0
860 outCalib.coeffErrSqr = inCalib.coeffErrSqr
861 outCalib.coeffErrSqr[~inCalib.coeffValid] = np.nan
863 outCalib.ampGainRatios = inCalib.ampGainRatios
864 outCalib.crosstalkRatiosUnits = inCalib.crosstalkRatiosUnits
866 outCalib.fitGains = inCalib.fitGains
868 return outCalib
870 def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False):
871 """Utility function to examine the final CT ratio set.
873 Parameters
874 ----------
875 stepname : `str`
876 State of processing to view.
877 ratios : `dict` [`dict` [`numpy.ndarray`]]
878 Array of measured CT ratios, indexed by source/victim
879 amplifier. These arrays are one-dimensional.
880 i : `str`
881 Index of the target amplifier.
882 j : `str`
883 Index of the source amplifier.
884 coeff : `float`, optional
885 Coefficient calculated to plot along with the simple mean.
886 valid : `bool`, optional
887 Validity to be added to the plot title.
888 """
889 frame = getDebugFrame(self._display, stepname)
890 if frame:
891 if i == j or ratios is None or len(ratios) < 1:
892 pass
894 ratioList = ratios[i][j]
895 if ratioList is None or len(ratioList) < 1:
896 pass
898 mean = np.mean(ratioList)
899 std = np.std(ratioList)
900 import matplotlib.pyplot as plt
901 figure = plt.figure(1)
902 figure.clear()
903 plt.hist(x=ratioList, bins=len(ratioList),
904 cumulative=True, color='b', density=True, histtype='step')
905 plt.xlabel("Measured pixel ratio")
906 plt.ylabel(f"CDF: n={len(ratioList)}")
907 plt.xlim(np.percentile(ratioList, [1.0, 99]))
908 plt.axvline(x=mean, color="k")
909 plt.axvline(x=coeff, color='g')
910 plt.axvline(x=(std / np.sqrt(len(ratioList))), color='r')
911 plt.axvline(x=-(std / np.sqrt(len(ratioList))), color='r')
912 plt.title(f"(Source {j} -> Target {i}) mean: {mean:.2g} coeff: {coeff:.2g} valid: {valid}")
913 figure.show()
915 prompt = "Press Enter to continue: "
916 while True:
917 ans = input(prompt).lower()
918 if ans in ("", "c",):
919 break
920 elif ans in ("pdb", "p",):
921 import pdb
922 pdb.set_trace()
923 plt.close()
926class CrosstalkFilterConnections(pipeBase.PipelineTaskConnections,
927 dimensions=("instrument", )):
928 inputCrosstalk = cT.Input(
929 name="crosstalkProposal",
930 doc="Input crosstalk calibrations as measured by CrosstalkSolveTask.",
931 storageClass="IsrCalib",
932 dimensions=("instrument", "detector"),
933 isCalibration=True,
934 multiple=True,
935 )
937 camera = cT.PrerequisiteInput(
938 name="camera",
939 doc="Camera containing cameraGeom information.",
940 storageClass="Camera",
941 dimensions=("instrument", ),
942 isCalibration=True,
943 )
945 outputCrosstalk = cT.Output(
946 name="crosstalk",
947 doc="Filtered crosstalk solutions.",
948 storageClass="CrosstalkCalib",
949 dimensions=("instrument", "detector"),
950 isCalibration=True,
951 multiple=True,
952 )
955class CrosstalkFilterConfig(pipeBase.PipelineTaskConfig,
956 pipelineConnections=CrosstalkFilterConnections):
957 """Configuration for the filtering of measured crosstalk solutions."""
958 doFiltering = Field(
959 dtype=bool,
960 default=True,
961 doc="Do filtering? If false, then this task acts as a pass-through to rename dataset types.",
962 )
964 nSigmaCoeffClip = Field(
965 dtype=float,
966 default=3.0,
967 doc="Coefficient outlier clipping significance.",
968 )
969 nSigmaCoeffSqrClip = Field(
970 dtype=float,
971 default=6.0,
972 doc="Squared-term coefficient outlier clipping significance.",
973 )
974 unitsAreElectrons = Field(
975 dtype=bool,
976 default=True,
977 doc="Crosstalk measurements have been done in electrons.",
978 )
981class CrosstalkFilterTask(pipeBase.PipelineTask):
982 """Task to compare crosstalk solutions between detectors, to identify
983 and remove outliers.
984 """
986 ConfigClass = CrosstalkFilterConfig
987 _DefaultName = 'cpCrosstalkFilter'
989 def runQuantum(self, butlerQC, inputRefs, outputRefs):
990 """Ensure that the input and output dimensions are passed along.
992 Parameters
993 ----------
994 butlerQC : `lsst.daf.butler.QuantumContext`
995 Butler to operate on.
996 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
997 Input data refs to load.
998 ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection`
999 Output data refs to persist.
1000 """
1001 inputs = butlerQC.get(inputRefs)
1003 # Use the dimensions to set calib/provenance information.
1004 inputs['inputDims'] = [dict(inCT.dataId.required) for inCT in inputRefs.inputCrosstalk]
1005 inputs['outputDims'] = [dict(outCT.dataId.required) for outCT in outputRefs.outputCrosstalk]
1007 outputs = self.run(**inputs)
1008 butlerQC.put(outputs, outputRefs)
1010 def run(self, inputCrosstalk, camera, inputDims, outputDims):
1011 """Compare crosstalk solutions to produce filtered crosstalk
1012 calibrations.
1014 Parameters
1015 ----------
1016 inputCrosstalk : `list` [`lsst.ip.isr.CrosstalkCalib`]
1017 List of crosstalk solutions to filter.
1018 camera : `lsst.afw.cameraGeom.Camera`
1019 Input camera.
1020 inputDims : `list` [`lsst.daf.butler.DataCoordinate`]
1021 DataIds to use to construct provenance.
1022 outputDims : `list` [`lsst.daf.butler.DataCoordinate`]
1023 DataIds to use to populate the output calibration.
1025 Returns
1026 -------
1027 results : `lsst.pipe.base.Struct`
1028 The results struct containing:
1030 ``outputCrosstalk``
1031 Final crosstalk calibration
1032 (`lsst.ip.isr.CrosstalkCalib`).
1034 Raises
1035 ------
1036 RuntimeError
1037 Raised if something goes bad. CZW/Fix me.
1038 """
1039 # These will hold all of the input data.
1040 itl_c0 = []
1041 e2v_c0 = []
1042 itl_c1 = []
1043 e2v_c1 = []
1044 detector_map = {}
1045 itl_counter = 0
1046 e2v_counter = 0
1047 for inputCT, inputDim in zip(inputCrosstalk, inputDims):
1048 detId = inputDim['detector']
1049 detector = camera[detId]
1051 if detector.getPhysicalType() == 'ITL':
1052 itl_c0.append(inputCT.coeffs)
1053 itl_c1.append(inputCT.coeffsSqr)
1055 detector_map[detId] = itl_counter
1056 itl_counter += 1
1058 elif detector.getPhysicalType() == 'E2V':
1059 e2v_c0.append(inputCT.coeffs)
1060 e2v_c1.append(inputCT.coeffsSqr)
1062 detector_map[detId] = e2v_counter
1063 e2v_counter += 1
1064 else:
1065 # This is a wavefront sensor, and we don't want to
1066 # filter those.
1067 pass
1069 itl_c0 = np.asarray(itl_c0)
1070 itl_c1 = np.asarray(itl_c1)
1071 e2v_c0 = np.asarray(e2v_c0)
1072 e2v_c1 = np.asarray(e2v_c1)
1074 itl_outliers = self.find_outliers(itl_c0, itl_c1)
1075 e2v_outliers = self.find_outliers(e2v_c0, e2v_c1)
1077 if self.config.doFiltering:
1078 itl_final = self.replace_outliers(itl_c0, itl_c1,
1079 itl_outliers.isBad, itl_outliers.median0, itl_outliers.median1)
1080 e2v_final = self.replace_outliers(e2v_c0, e2v_c1,
1081 e2v_outliers.isBad, e2v_outliers.median0, e2v_outliers.median1)
1082 outputCrosstalkList = []
1083 for inputCT, inputDim, outputDim in zip(inputCrosstalk, inputDims, outputDims):
1084 if inputDim['detector'] != outputDim['detector']:
1085 raise RuntimeError("Inconsistent dimension records")
1087 detId = inputDim['detector']
1088 detector = camera[detId]
1090 outputCT = copy(inputCT)
1092 if detector.getPhysicalType() == 'ITL':
1093 itl_index = detector_map[detId]
1094 if (np.any(itl_final.new_matrix0[itl_index] != outputCT.coeffs)
1095 or np.any(itl_final.new_matrix1[itl_index] != outputCT.coeffsSqr)):
1096 outputCT.coeffs = itl_final.new_matrix0[itl_index]
1097 outputCT.coeffsSqr = itl_final.new_matrix1[itl_index]
1098 elif detector.getPhysicalType() == 'E2V':
1099 e2v_index = detector_map[detId]
1100 if (np.any(e2v_final.new_matrix0[e2v_index] != outputCT.coeffs)
1101 or np.any(e2v_final.new_matrix1[e2v_index] != outputCT.coeffsSqr)):
1102 outputCT.coeffs = e2v_final.new_matrix0[e2v_index]
1103 outputCT.coeffsSqr = e2v_final.new_matrix1[e2v_index]
1105 outputCT.crosstalkRatiosUnits = 'electron' if self.config.unitsAreElectrons else 'adu'
1106 outputCT.updateMetadata(
1107 camera=camera,
1108 detector=camera[detId],
1109 setCalibId=True,
1110 setCalibInfo=True,
1111 setDate=True,
1112 )
1113 outputCrosstalkList.append(outputCT)
1114 return pipeBase.Struct(
1115 outputCrosstalk=outputCrosstalkList,
1116 )
1118 def find_outliers(self, matrix0, matrix1):
1119 """Do checks to see if an element of the matrix is out-of-family.
1121 Parameters
1122 ----------
1123 matrix0 : `np.ndarray`, (Ndet, Namp, Namp)
1124 Matrix holding the 0th-order terms.
1125 matrix1 : `np.ndarray`, (Ndet, Namp, Namp)
1126 Matrix holding the 1st-order terms.
1128 Returns
1129 -------
1130 results : `lsst.pipe.base.Struct`
1131 Results struct containing
1133 ``median0``
1134 Median in-family value (`np.ndarray` (Namp, Namp)).
1135 ``stdev0``
1136 MAD effective sigma in-family value (`np.ndarray`
1137 (Namp, Namp)).
1138 ``median1``
1139 Median in-family value (`np.ndarray` (Namp, Namp)).
1140 ``stdev1``
1141 MAD effective sigma in-family value (`np.ndarray`
1142 (Namp, Namp)).
1143 ``isBad``
1144 Boolean indicator that an element has been replaced
1145 (`np.ndarray` (Ndet, Namp, Namp)).
1147 Raises
1148 ------
1149 ValueError :
1150 Raised if the inputs have a mismatch in size.
1151 """
1152 if matrix0.shape != matrix1.shape:
1153 raise ValueError("Shape disagreement!")
1154 if matrix0.shape[1] != matrix0.shape[2]:
1155 raise ValueError("Shape disagreement!")
1157 nAmp = matrix0.shape[1]
1159 median0 = np.nanmedian(matrix0, axis=0)
1160 median1 = np.nanmedian(matrix1, axis=0)
1162 stdev0 = median_abs_deviation(matrix0, axis=0,
1163 center=np.nanmedian, scale="normal", nan_policy="omit")
1164 stdev1 = median_abs_deviation(matrix1, axis=0,
1165 center=np.nanmedian, scale="normal", nan_policy="omit")
1167 isBad = np.full_like(matrix0, False, dtype=bool)
1169 for i in range(nAmp):
1170 for j in range(nAmp):
1171 m0 = median0[i][j]
1172 m1 = median1[i][j]
1173 s0 = stdev0[i][j]
1174 s1 = stdev1[i][j]
1176 min0 = m0 - self.config.nSigmaCoeffClip * s0
1177 max0 = m0 + self.config.nSigmaCoeffClip * s0
1179 min1 = m1 - self.config.nSigmaCoeffSqrClip * s1
1180 max1 = m1 + self.config.nSigmaCoeffSqrClip * s1
1182 bad, = np.where(
1183 (matrix0[:, i, j] < min0)
1184 | (matrix0[:, i, j] > max0)
1185 | (matrix1[:, i, j] < min1)
1186 | (matrix1[:, i, j] > max1)
1187 )
1189 if len(bad) > 0:
1190 for detIdx in bad:
1191 isBad[detIdx, i, j] = True
1192 return pipeBase.Struct(
1193 median0=median0,
1194 median1=median1,
1195 stdev0=stdev0,
1196 stdev1=stdev1,
1197 isBad=isBad
1198 )
1200 def replace_outliers(self, matrix0, matrix1, isBad, median0, median1):
1201 """Do checks to see if an element of the matrix is out-of-family.
1203 Parameters
1204 ----------
1205 matrix0 : `np.ndarray`, (Ndet, Namp, Namp)
1206 Matrix holding the 0th-order terms.
1207 matrix1 : `np.ndarray`, (Ndet, Namp, Namp)
1208 Matrix holding the 1st-order terms.
1209 isBad : `np.ndarray`, (Ndet, Namp, Namp)
1210 Matrix holding the boolean "is bad".
1211 median0 : `np.ndarray`, (Namp, Namp)
1212 Matrix of median 0th-order terms.
1213 median1 : `np.ndarray`, (Namp, Namp)
1214 Matrix of median 1st-order terms.
1216 Returns
1217 -------
1218 results : `lsst.pipe.base.Struct`
1219 Results struct containing
1221 ``new_matrix0``
1222 Replacement matrix0, with median substitutions.
1223 (`np.ndarray` (Ndet, Namp, Namp)).
1224 ``new_matrix1``
1225 Replacement matrix1, with median substitutions.
1226 (`np.ndarray` (Ndet, Namp, Namp)).
1228 Raises
1229 ------
1230 ValueError :
1231 Raised if the inputs have a mismatch in size.
1232 """
1233 if matrix0.shape != matrix1.shape:
1234 raise ValueError("Shape disagreement!")
1235 if matrix0.shape != isBad.shape:
1236 raise ValueError("Shape disagreement!")
1237 if median0.shape != median1.shape:
1238 raise ValueError("Shape disagreement!")
1240 out0 = np.full_like(matrix0, 0.0)
1241 out1 = np.full_like(matrix1, 0.0)
1243 for detIdx in range(matrix0.shape[0]):
1244 for srcIdx in range(matrix0.shape[1]):
1245 for tgtIdx in range(matrix0.shape[2]):
1246 if isBad[detIdx, srcIdx, tgtIdx]:
1247 self.log.info(f"Setting {detIdx} {srcIdx} {tgtIdx} from "
1248 f"{matrix0[detIdx, srcIdx, tgtIdx]} to "
1249 f"{median0[srcIdx, tgtIdx]} and "
1250 f"{matrix1[detIdx, srcIdx, tgtIdx]} to "
1251 f"{median1[srcIdx, tgtIdx]}")
1252 out0[detIdx, srcIdx, tgtIdx] = median0[srcIdx, tgtIdx]
1253 out1[detIdx, srcIdx, tgtIdx] = median1[srcIdx, tgtIdx]
1254 else:
1255 out0[detIdx, srcIdx, tgtIdx] = matrix0[detIdx, srcIdx, tgtIdx]
1256 out1[detIdx, srcIdx, tgtIdx] = matrix1[detIdx, srcIdx, tgtIdx]
1257 return pipeBase.Struct(
1258 new_matrix0=out0,
1259 new_matrix1=out1,
1260 )