Coverage for python / lsst / cp / pipe / cpCrosstalk.py: 15%
293 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:07 +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
26import lsst.pipe.base as pipeBase
27import lsst.pipe.base.connectionTypes as cT
29from lsstDebug import getDebugFrame
30from lsst.afw.detection import FootprintSet, Threshold
31from lsst.afw.display import getDisplay
32from lsst.pex.config import ConfigurableField, Field, ListField
33from lsst.ip.isr import CrosstalkCalib, IsrProvenance
34from lsst.cp.pipe.utils import (ddict2dict, sigmaClipCorrection)
35from lsst.meas.algorithms import SubtractBackgroundTask
37__all__ = ["CrosstalkExtractConfig", "CrosstalkExtractTask",
38 "CrosstalkSolveTask", "CrosstalkSolveConfig"]
41class CrosstalkExtractConnections(pipeBase.PipelineTaskConnections,
42 dimensions=("instrument", "exposure", "detector")):
43 inputExp = cT.Input(
44 name="crosstalkInputs",
45 doc="Input post-ISR processed exposure to measure crosstalk from.",
46 storageClass="Exposure",
47 dimensions=("instrument", "exposure", "detector"),
48 multiple=False,
49 )
50 # TODO: Depends on DM-21904.
51 sourceExp = cT.Input(
52 name="crosstalkSource",
53 doc="Post-ISR exposure to measure for inter-chip crosstalk onto inputExp.",
54 storageClass="Exposure",
55 dimensions=("instrument", "exposure", "detector"),
56 multiple=True,
57 deferLoad=True,
58 # lookupFunction=None,
59 )
61 outputRatios = cT.Output(
62 name="crosstalkRatios",
63 doc="Extracted crosstalk pixel ratios.",
64 storageClass="StructuredDataDict",
65 dimensions=("instrument", "exposure", "detector"),
66 )
67 outputFluxes = cT.Output(
68 name="crosstalkFluxes",
69 doc="Source pixel fluxes used in ratios.",
70 storageClass="StructuredDataDict",
71 dimensions=("instrument", "exposure", "detector"),
72 )
74 def __init__(self, *, config=None):
75 super().__init__(config=config)
76 # Discard sourceExp until DM-21904 allows full interchip
77 # measurements.
78 self.inputs.discard("sourceExp")
81class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig,
82 pipelineConnections=CrosstalkExtractConnections):
83 """Configuration for the measurement of pixel ratios.
84 """
86 doMeasureInterchip = Field(
87 dtype=bool,
88 default=False,
89 doc="Measure inter-chip crosstalk as well?",
90 )
91 threshold = Field(
92 dtype=float,
93 default=30000,
94 doc="Minimum level of source pixels for which to measure crosstalk."
95 )
96 ignoreSaturatedPixels = Field(
97 dtype=bool,
98 default=False,
99 doc="Should saturated pixels be ignored?"
100 )
101 badMask = ListField(
102 dtype=str,
103 default=["BAD", "INTRP"],
104 doc="Mask planes to ignore when identifying source pixels."
105 )
106 isTrimmed = Field(
107 dtype=bool,
108 default=True,
109 doc="Is the input exposure trimmed?"
110 )
111 background = ConfigurableField(
112 target=SubtractBackgroundTask,
113 doc="Background estimation task.",
114 )
116 def validate(self):
117 super().validate()
119 # Ensure the handling of the SAT mask plane is consistent
120 # with the ignoreSaturatedPixels value.
121 if self.ignoreSaturatedPixels:
122 if 'SAT' not in self.badMask:
123 self.badMask.append('SAT')
124 else:
125 if 'SAT' in self.badMask:
126 self.badMask = [mask for mask in self.badMask if mask != 'SAT']
129class CrosstalkExtractTask(pipeBase.PipelineTask):
130 """Task to measure pixel ratios to find crosstalk.
131 """
133 ConfigClass = CrosstalkExtractConfig
134 _DefaultName = 'cpCrosstalkExtract'
136 def __init__(self, **kwargs):
137 super().__init__(**kwargs)
138 self.makeSubtask("background")
140 def run(self, inputExp, sourceExps=[]):
141 """Measure pixel ratios between amplifiers in inputExp.
143 Extract crosstalk ratios between different amplifiers.
145 For pixels above ``config.threshold``, we calculate the ratio
146 between each background-subtracted target amp and the source
147 amp. We return a list of ratios for each pixel for each
148 target/source combination, as nested dictionary containing the
149 ratio.
151 Parameters
152 ----------
153 inputExp : `lsst.afw.image.Exposure`
154 Input exposure to measure pixel ratios on.
155 sourceExp : `list` [`lsst.afw.image.Exposure`], optional
156 List of chips to use as sources to measure inter-chip
157 crosstalk.
159 Returns
160 -------
161 results : `lsst.pipe.base.Struct`
162 The results struct containing:
164 ``outputRatios``
165 A catalog of ratio lists. The dictionaries are
166 indexed such that:
167 outputRatios[targetChip][sourceChip][targetAmp][sourceAmp]
168 contains the ratio list for that combination (`dict`
169 [`dict` [`dict` [`dict` [`list`]]]]).
170 ``outputFluxes``
171 A catalog of flux lists. The dictionaries are
172 indexed such that:
173 outputFluxes[sourceChip][sourceAmp] contains the flux
174 list used in the outputRatios (`dict` [`dict`
175 [`list`]]).
176 """
177 outputRatios = defaultdict(lambda: defaultdict(dict))
178 outputFluxes = defaultdict(lambda: defaultdict(dict))
180 threshold = self.config.threshold
181 badPixels = list(self.config.badMask)
183 targetDetector = inputExp.getDetector()
184 targetChip = targetDetector.getName()
186 # Always look at the target chip first, then go to any other
187 # supplied exposures.
188 sourceExtractExps = [inputExp]
189 sourceExtractExps.extend(sourceExps)
191 self.log.info("Measuring full detector background for target: %s", targetChip)
192 targetIm = inputExp.getMaskedImage()
193 FootprintSet(targetIm, Threshold(threshold), "DETECTED")
194 detected = targetIm.getMask().getPlaneBitMask("DETECTED")
195 bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"])
196 backgroundModel = self.background.fitBackground(inputExp.maskedImage)
197 backgroundIm = backgroundModel.getImageF()
199 self.debugView('extract', inputExp)
201 for sourceExp in sourceExtractExps:
202 sourceDetector = sourceExp.getDetector()
203 sourceChip = sourceDetector.getName()
204 sourceIm = sourceExp.getMaskedImage()
205 bad = sourceIm.getMask().getPlaneBitMask(badPixels)
206 self.log.info("Measuring crosstalk from source: %s", sourceChip)
208 if sourceExp != inputExp:
209 FootprintSet(sourceIm, Threshold(threshold), "DETECTED")
210 detected = sourceIm.getMask().getPlaneBitMask("DETECTED")
212 # The dictionary of amp-to-amp ratios for this pair of
213 # source->target detectors.
214 ratioDict = defaultdict(lambda: defaultdict(list))
215 extractedCount = 0
217 for sourceAmp in sourceDetector:
218 sourceAmpName = sourceAmp.getName()
219 sourceAmpBBox = sourceAmp.getBBox() if self.config.isTrimmed else sourceAmp.getRawDataBBox()
220 sourceAmpImage = sourceIm[sourceAmpBBox]
221 sourceMask = sourceAmpImage.mask.array
222 select = ((sourceMask & detected > 0)
223 & (sourceMask & bad == 0)
224 & np.isfinite(sourceAmpImage.image.array))
225 count = np.sum(select)
226 self.log.debug(" Source amplifier: %s", sourceAmpName)
228 outputFluxes[sourceChip][sourceAmpName] = sourceAmpImage.image.array[select].tolist()
230 for targetAmp in targetDetector:
231 # iterate over targetExposure
232 targetAmpName = targetAmp.getName()
233 if sourceAmpName == targetAmpName and sourceChip == targetChip:
234 ratioDict[targetAmpName][sourceAmpName] = []
235 continue
237 self.log.debug(" Target amplifier: %s", targetAmpName)
239 targetAmpImage = CrosstalkCalib.extractAmp(targetIm,
240 targetAmp, sourceAmp,
241 isTrimmed=self.config.isTrimmed)
242 targetBkgImage = CrosstalkCalib.extractAmp(backgroundIm,
243 targetAmp, sourceAmp,
244 isTrimmed=self.config.isTrimmed)
246 bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"])
248 ratios = ((targetAmpImage.image.array[select] - targetBkgImage.array[select])
249 / sourceAmpImage.image.array[select])
251 ratioDict[targetAmpName][sourceAmpName] = ratios.tolist()
252 self.log.info("Amp extracted %d pixels from %s -> %s",
253 count, sourceAmpName, targetAmpName)
254 extractedCount += count
256 self.debugPixels('pixels',
257 sourceAmpImage.image.array[select],
258 targetAmpImage.image.array[select] - bg,
259 sourceAmpName, targetAmpName)
261 self.log.info("Extracted %d pixels from %s -> %s (targetBG: %f)",
262 extractedCount, sourceChip, targetChip, bg)
263 outputRatios[targetChip][sourceChip] = ratioDict
265 return pipeBase.Struct(
266 outputRatios=ddict2dict(outputRatios),
267 outputFluxes=ddict2dict(outputFluxes)
268 )
270 def debugView(self, stepname, exposure):
271 """Utility function to examine the image being processed.
273 Parameters
274 ----------
275 stepname : `str`
276 State of processing to view.
277 exposure : `lsst.afw.image.Exposure`
278 Exposure to view.
279 """
280 frame = getDebugFrame(self._display, stepname)
281 if frame:
282 display = getDisplay(frame)
283 display.scale('asinh', 'zscale')
284 display.mtv(exposure)
286 prompt = "Press Enter to continue: "
287 while True:
288 ans = input(prompt).lower()
289 if ans in ("", "c",):
290 break
292 def debugPixels(self, stepname, pixelsIn, pixelsOut, sourceName, targetName):
293 """Utility function to examine the CT ratio pixel values.
295 Parameters
296 ----------
297 stepname : `str`
298 State of processing to view.
299 pixelsIn : `np.ndarray`, (N,)
300 Pixel values from the potential crosstalk source.
301 pixelsOut : `np.ndarray`, (N,)
302 Pixel values from the potential crosstalk target.
303 sourceName : `str`
304 Source amplifier name
305 targetName : `str`
306 Target amplifier name
307 """
308 frame = getDebugFrame(self._display, stepname)
309 if frame:
310 import matplotlib.pyplot as plt
311 figure = plt.figure(1)
312 figure.clear()
314 axes = figure.add_axes((0.1, 0.1, 0.8, 0.8))
315 axes.plot(pixelsIn, pixelsOut / pixelsIn, 'k+')
316 plt.xlabel("Source amplifier pixel value")
317 plt.ylabel("Measured pixel ratio")
318 plt.title(f"(Source {sourceName} -> Target {targetName}) median ratio: "
319 f"{(np.median(pixelsOut / pixelsIn))}")
320 figure.show()
322 prompt = "Press Enter to continue: "
323 while True:
324 ans = input(prompt).lower()
325 if ans in ("", "c",):
326 break
327 plt.close()
330class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections,
331 dimensions=("instrument", "detector")):
332 inputRatios = cT.Input(
333 name="crosstalkRatios",
334 doc="Ratios measured for an input exposure.",
335 storageClass="StructuredDataDict",
336 dimensions=("instrument", "exposure", "detector"),
337 multiple=True,
338 )
339 inputFluxes = cT.Input(
340 name="crosstalkFluxes",
341 doc="Fluxes of CT source pixels, for nonlinear fits.",
342 storageClass="StructuredDataDict",
343 dimensions=("instrument", "exposure", "detector"),
344 multiple=True,
345 )
346 camera = cT.PrerequisiteInput(
347 name="camera",
348 doc="Camera the input data comes from.",
349 storageClass="Camera",
350 dimensions=("instrument",),
351 isCalibration=True,
352 )
354 outputCrosstalk = cT.Output(
355 name="crosstalk",
356 doc="Output proposed crosstalk calibration.",
357 storageClass="CrosstalkCalib",
358 dimensions=("instrument", "detector"),
359 multiple=False,
360 isCalibration=True,
361 )
363 def __init__(self, *, config=None):
364 super().__init__(config=config)
366 if config.fluxOrder == 0:
367 self.inputs.discard("inputFluxes")
370class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig,
371 pipelineConnections=CrosstalkSolveConnections):
372 """Configuration for the solving of crosstalk from pixel ratios.
373 """
375 rejIter = Field(
376 dtype=int,
377 default=3,
378 doc="Number of rejection iterations for final coefficient calculation.",
379 )
380 rejSigma = Field(
381 dtype=float,
382 default=2.0,
383 doc="Rejection threshold (sigma) for final coefficient calculation.",
384 )
385 fluxOrder = Field(
386 dtype=int,
387 default=0,
388 doc="Polynomial order in source flux to fit crosstalk.",
389 )
391 rejectNegativeSolutions = Field(
392 dtype=bool,
393 default=True,
394 doc="Should solutions with negative coefficients (which add flux to the target) be excluded?",
395 )
397 significanceLimit = Field(
398 dtype=float,
399 default=3.0,
400 doc="Sigma significance level to use in marking a coefficient valid.",
401 )
402 doSignificanceScaling = Field(
403 dtype=bool,
404 default=True,
405 doc="Scale error by 1/sqrt(N) in calculating significant coefficients?",
406 )
407 doFiltering = Field(
408 dtype=bool,
409 default=False,
410 doc="Filter generated crosstalk to remove marginal measurements?",
411 )
414class CrosstalkSolveTask(pipeBase.PipelineTask):
415 """Task to solve crosstalk from pixel ratios.
416 """
418 ConfigClass = CrosstalkSolveConfig
419 _DefaultName = 'cpCrosstalkSolve'
421 def runQuantum(self, butlerQC, inputRefs, outputRefs):
422 """Ensure that the input and output dimensions are passed along.
424 Parameters
425 ----------
426 butlerQC : `lsst.daf.butler.QuantumContext`
427 Butler to operate on.
428 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
429 Input data refs to load.
430 ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection`
431 Output data refs to persist.
432 """
433 inputs = butlerQC.get(inputRefs)
435 # Use the dimensions to set calib/provenance information.
436 inputs['inputDims'] = [dict(exp.dataId.required) for exp in inputRefs.inputRatios]
437 inputs['outputDims'] = dict(outputRefs.outputCrosstalk.dataId.required)
439 outputs = self.run(**inputs)
440 butlerQC.put(outputs, outputRefs)
442 def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, outputDims=None):
443 """Combine ratios to produce crosstalk coefficients.
445 Parameters
446 ----------
447 inputRatios : `list` [`dict` [`dict` [`dict` [`dict` [`list`]]]]]
448 A list of nested dictionaries of ratios indexed by target
449 and source chip, then by target and source amplifier.
450 inputFluxes : `list` [`dict` [`dict` [`list`]]]
451 A list of nested dictionaries of source pixel fluxes, indexed
452 by source chip and amplifier.
453 camera : `lsst.afw.cameraGeom.Camera`
454 Input camera.
455 inputDims : `list` [`lsst.daf.butler.DataCoordinate`]
456 DataIds to use to construct provenance.
457 outputDims : `list` [`lsst.daf.butler.DataCoordinate`]
458 DataIds to use to populate the output calibration.
460 Returns
461 -------
462 results : `lsst.pipe.base.Struct`
463 The results struct containing:
465 ``outputCrosstalk``
466 Final crosstalk calibration
467 (`lsst.ip.isr.CrosstalkCalib`).
468 ``outputProvenance``
469 Provenance data for the new calibration
470 (`lsst.ip.isr.IsrProvenance`).
472 Raises
473 ------
474 RuntimeError
475 Raised if the input data contains multiple target detectors.
476 """
477 if outputDims:
478 calibChip = outputDims['detector']
479 instrument = outputDims['instrument']
480 else:
481 # calibChip needs to be set manually in Gen2.
482 calibChip = None
483 instrument = None
485 if camera and calibChip is not None:
486 calibDetector = camera[calibChip]
487 ordering = [amp.getName() for amp in calibDetector]
488 else:
489 calibDetector = None
490 ordering = None
492 self.log.info("Combining measurements from %d ratios and %d fluxes",
493 len(inputRatios), len(inputFluxes) if inputFluxes else 0)
495 if inputFluxes is None:
496 inputFluxes = [None for exp in inputRatios]
498 combinedRatios = defaultdict(lambda: defaultdict(list))
499 combinedFluxes = defaultdict(lambda: defaultdict(list))
501 for ratioDict, fluxDict in zip(inputRatios, inputFluxes):
502 for targetChip in ratioDict:
503 if calibChip and targetChip != calibChip and targetChip != calibDetector.getName():
504 raise RuntimeError(f"Target chip: {targetChip} does not match calibration dimension: "
505 f"{calibChip}, {calibDetector.getName()}!")
507 sourceChip = targetChip
508 if sourceChip in ratioDict[targetChip]:
509 ratios = ratioDict[targetChip][sourceChip]
511 for targetAmp in ratios:
512 for sourceAmp in ratios[targetAmp]:
513 combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp])
514 if fluxDict:
515 combinedFluxes[targetAmp][sourceAmp].extend(fluxDict[sourceChip][sourceAmp])
516 # TODO: DM-21904
517 # Iterating over all other entries in
518 # ratioDict[targetChip] will yield inter-chip terms.
520 for targetAmp in combinedRatios:
521 for sourceAmp in combinedRatios[targetAmp]:
522 self.log.info("Read %d pixels for %s -> %s",
523 len(combinedRatios[targetAmp][sourceAmp]),
524 sourceAmp, targetAmp)
525 if len(combinedRatios[targetAmp][sourceAmp]) > 1:
526 self.debugRatios('reduce', combinedRatios, targetAmp, sourceAmp)
528 if self.config.fluxOrder == 0:
529 self.log.info("Fitting crosstalk coefficients.")
531 calib = self.measureCrosstalkCoefficients(combinedRatios, ordering,
532 self.config.rejIter, self.config.rejSigma)
533 else:
534 raise NotImplementedError("Non-linear crosstalk terms are not yet supported.")
536 self.log.info("Number of valid coefficients: %d", np.sum(calib.coeffValid))
538 if self.config.doFiltering:
539 # This step will apply the calculated validity values to
540 # censor poorly measured coefficients.
541 self.log.info("Filtering measured crosstalk to remove invalid solutions.")
542 calib = self.filterCrosstalkCalib(calib)
544 # Populate the remainder of the calibration information.
545 calib.hasCrosstalk = True
546 calib.interChip = {}
548 calib.updateMetadata(camera=camera, detector=calibDetector)
549 calib.updateMetadata(setCalibId=True, setDate=True)
551 # Make an IsrProvenance().
552 provenance = IsrProvenance(calibType="CROSSTALK")
553 provenance._detectorName = calibChip
554 if inputDims:
555 provenance.fromDataIds(inputDims)
556 provenance._instrument = instrument
557 provenance.updateMetadata()
559 return pipeBase.Struct(
560 outputCrosstalk=calib,
561 outputProvenance=provenance,
562 )
564 def measureCrosstalkCoefficients(self, ratios, ordering, rejIter, rejSigma):
565 """Measure crosstalk coefficients from the ratios.
567 Given a list of ratios for each target/source amp combination,
568 we measure a sigma clipped mean and error.
570 The coefficient errors returned are the standard deviation of
571 the final set of clipped input ratios.
573 Parameters
574 ----------
575 ratios : `dict` [`dict` [`numpy.ndarray`]]
576 Catalog of arrays of ratios. The ratio arrays are one-dimensional
577 ordering : `list` [`str`] or None
578 List to use as a mapping between amplifier names (the
579 elements of the list) and their position in the output
580 calibration (the matching index of the list). If no
581 ordering is supplied, the order of the keys in the ratio
582 catalog is used.
583 rejIter : `int`
584 Number of rejection iterations.
585 rejSigma : `float`
586 Rejection threshold (sigma).
588 Returns
589 -------
590 calib : `lsst.ip.isr.CrosstalkCalib`
591 The output crosstalk calibration.
592 """
593 calib = CrosstalkCalib(nAmp=len(ratios))
595 if ordering is None:
596 ordering = list(ratios.keys())
598 # Calibration stores coefficients as a numpy ndarray.
599 for ss, tt in itertools.product(range(calib.nAmp), range(calib.nAmp)):
600 if ss == tt:
601 values = [0.0]
602 else:
603 # ratios is ratios[Target][Source]
604 # use tt for Target, use ss for Source, to match ip_isr.
605 values = np.array(ratios[ordering[tt]][ordering[ss]])
606 values = values[np.abs(values) < 1.0] # Discard unreasonable values
608 # Sigma clip using the inter-quartile distance and a
609 # normal distribution.
610 if ss != tt:
611 for rej in range(rejIter):
612 if len(values) == 0:
613 break
614 lo, med, hi = np.percentile(values, [25.0, 50.0, 75.0])
615 sigma = 0.741*(hi - lo)
616 good = np.abs(values - med) < rejSigma*sigma
617 if good.sum() == len(good) or good.sum() == 0:
618 break
619 values = values[good]
621 # Crosstalk calib is property[Source][Target].
622 calib.coeffNum[ss][tt] = len(values)
623 significanceThreshold = 0.0
624 if len(values) == 0:
625 self.log.warning("No values for matrix element %d,%d" % (ss, tt))
626 calib.coeffs[ss][tt] = np.nan
627 calib.coeffErr[ss][tt] = np.nan
628 calib.coeffValid[ss][tt] = False
629 else:
630 calib.coeffs[ss][tt] = np.mean(values)
631 if self.config.rejectNegativeSolutions and calib.coeffs[ss][tt] < 0.0:
632 calib.coeffs[ss][tt] = 0.0
634 if calib.coeffNum[ss][tt] == 1:
635 calib.coeffErr[ss][tt] = np.nan
636 calib.coeffValid[ss][tt] = False
637 else:
638 correctionFactor = sigmaClipCorrection(rejSigma)
639 calib.coeffErr[ss][tt] = np.std(values) * correctionFactor
641 # Use sample stdev.
642 significanceThreshold = self.config.significanceLimit * calib.coeffErr[ss][tt]
643 if self.config.doSignificanceScaling is True:
644 # Enabling this calculates the stdev of the mean.
645 significanceThreshold /= np.sqrt(calib.coeffNum[ss][tt])
646 calib.coeffValid[ss][tt] = np.abs(calib.coeffs[ss][tt]) > significanceThreshold
647 self.debugRatios('measure', ratios, ordering[ss], ordering[tt],
648 calib.coeffs[ss][tt], calib.coeffValid[ss][tt])
649 self.log.info("Measured %s -> %s Coeff: %e Err: %e N: %d Valid: %s Limit: %e",
650 ordering[ss], ordering[tt], calib.coeffs[ss][tt], calib.coeffErr[ss][tt],
651 calib.coeffNum[ss][tt], calib.coeffValid[ss][tt], significanceThreshold)
653 return calib
655 @staticmethod
656 def filterCrosstalkCalib(inCalib):
657 """Apply valid constraints to the measured values.
659 Any measured coefficient that is determined to be invalid is
660 set to zero, and has the error set to nan. The validation is
661 determined by checking that the measured coefficient is larger
662 than the calculated standard error of the mean.
664 Parameters
665 ----------
666 inCalib : `lsst.ip.isr.CrosstalkCalib`
667 Input calibration to filter.
669 Returns
670 -------
671 outCalib : `lsst.ip.isr.CrosstalkCalib`
672 Filtered calibration.
673 """
674 outCalib = CrosstalkCalib()
675 outCalib.nAmp = inCalib.nAmp
677 outCalib.coeffs = inCalib.coeffs
678 outCalib.coeffs[~inCalib.coeffValid] = 0.0
680 outCalib.coeffErr = inCalib.coeffErr
681 outCalib.coeffErr[~inCalib.coeffValid] = np.nan
683 outCalib.coeffNum = inCalib.coeffNum
684 outCalib.coeffValid = inCalib.coeffValid
686 outCalib.coeffsSqr = inCalib.coeffsSqr
687 outCalib.coeffsSqr[~inCalib.coeffValid] = 0.0
689 outCalib.coeffErrSqr = inCalib.coeffErrSqr
690 outCalib.coeffErrSqr[~inCalib.coeffValid] = np.nan
692 outCalib.ampGainRatios = inCalib.ampGainRatios
693 outCalib.crosstalkRatiosUnits = inCalib.crosstalkRatiosUnits
695 outCalib.fitGains = inCalib.fitGains
697 return outCalib
699 def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False):
700 """Utility function to examine the final CT ratio set.
702 Parameters
703 ----------
704 stepname : `str`
705 State of processing to view.
706 ratios : `dict` [`dict` [`numpy.ndarray`]]
707 Array of measured CT ratios, indexed by source/victim
708 amplifier. These arrays are one-dimensional.
709 i : `str`
710 Index of the target amplifier.
711 j : `str`
712 Index of the source amplifier.
713 coeff : `float`, optional
714 Coefficient calculated to plot along with the simple mean.
715 valid : `bool`, optional
716 Validity to be added to the plot title.
717 """
718 frame = getDebugFrame(self._display, stepname)
719 if frame:
720 if i == j or ratios is None or len(ratios) < 1:
721 pass
723 ratioList = ratios[i][j]
724 if ratioList is None or len(ratioList) < 1:
725 pass
727 mean = np.mean(ratioList)
728 std = np.std(ratioList)
729 import matplotlib.pyplot as plt
730 figure = plt.figure(1)
731 figure.clear()
732 plt.hist(x=ratioList, bins=len(ratioList),
733 cumulative=True, color='b', density=True, histtype='step')
734 plt.xlabel("Measured pixel ratio")
735 plt.ylabel(f"CDF: n={len(ratioList)}")
736 plt.xlim(np.percentile(ratioList, [1.0, 99]))
737 plt.axvline(x=mean, color="k")
738 plt.axvline(x=coeff, color='g')
739 plt.axvline(x=(std / np.sqrt(len(ratioList))), color='r')
740 plt.axvline(x=-(std / np.sqrt(len(ratioList))), color='r')
741 plt.title(f"(Source {j} -> Target {i}) mean: {mean:.2g} coeff: {coeff:.2g} valid: {valid}")
742 figure.show()
744 prompt = "Press Enter to continue: "
745 while True:
746 ans = input(prompt).lower()
747 if ans in ("", "c",):
748 break
749 elif ans in ("pdb", "p",):
750 import pdb
751 pdb.set_trace()
752 plt.close()