Coverage for python/lsst/cp/pipe/cpBfk.py: 11%
296 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-19 02:31 -0800
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 <https://www.gnu.org/licenses/>.
21#
22"""Calculation of brighter-fatter effect correlations and kernels."""
24__all__ = ['BrighterFatterKernelSolveTask',
25 'BrighterFatterKernelSolveConfig']
27import numpy as np
29import lsst.afw.math as afwMath
30import lsst.pex.config as pexConfig
31import lsst.pipe.base as pipeBase
32import lsst.pipe.base.connectionTypes as cT
34from lsst.ip.isr import (BrighterFatterKernel)
35from .utils import (funcPolynomial, irlsFit, extractCalibDate)
36from .cpLinearitySolve import ptcLookup
39class BrighterFatterKernelSolveConnections(pipeBase.PipelineTaskConnections,
40 dimensions=("instrument", "detector")):
41 dummy = cT.Input(
42 name="raw",
43 doc="Dummy exposure.",
44 storageClass='Exposure',
45 dimensions=("instrument", "exposure", "detector"),
46 multiple=True,
47 deferLoad=True,
48 )
49 camera = cT.PrerequisiteInput(
50 name="camera",
51 doc="Camera associated with this data.",
52 storageClass="Camera",
53 dimensions=("instrument", ),
54 isCalibration=True,
55 )
56 inputPtc = cT.PrerequisiteInput(
57 name="ptc",
58 doc="Photon transfer curve dataset.",
59 storageClass="PhotonTransferCurveDataset",
60 dimensions=("instrument", "detector"),
61 isCalibration=True,
62 lookupFunction=ptcLookup,
63 )
64 inputBfkPtc = cT.Input(
65 name="bfkPtc",
66 doc="Input BFK PTC dataset.",
67 storageClass="PhotonTransferCurveDataset",
68 dimensions=("instrument", "detector"),
69 isCalibration=True,
70 )
72 outputBFK = cT.Output(
73 name="brighterFatterKernel",
74 doc="Output measured brighter-fatter kernel.",
75 storageClass="BrighterFatterKernel",
76 dimensions=("instrument", "detector"),
77 isCalibration=True,
78 )
80 def __init__(self, *, config=None):
81 if config.useBfkPtc:
82 del self.inputPtc
83 del self.dummy
84 else:
85 del self.inputBfkPtc
88class BrighterFatterKernelSolveConfig(pipeBase.PipelineTaskConfig,
89 pipelineConnections=BrighterFatterKernelSolveConnections):
90 level = pexConfig.ChoiceField(
91 doc="The level at which to calculate the brighter-fatter kernels",
92 dtype=str,
93 default="AMP",
94 allowed={
95 "AMP": "Every amplifier treated separately",
96 "DETECTOR": "One kernel per detector",
97 }
98 )
99 ignoreAmpsForAveraging = pexConfig.ListField(
100 dtype=str,
101 doc="List of amp names to ignore when averaging the amplifier kernels into the detector"
102 " kernel. Only relevant for level = DETECTOR",
103 default=[]
104 )
105 xcorrCheckRejectLevel = pexConfig.Field(
106 dtype=float,
107 doc="Rejection level for the sum of the input cross-correlations. Arrays which "
108 "sum to greater than this are discarded before the clipped mean is calculated.",
109 default=2.0
110 )
111 nSigmaClip = pexConfig.Field(
112 dtype=float,
113 doc="Number of sigma to clip when calculating means for the cross-correlation",
114 default=5
115 )
116 forceZeroSum = pexConfig.Field(
117 dtype=bool,
118 doc="Force the correlation matrix to have zero sum by adjusting the (0,0) value?"
119 "Defaults to true bsed on recommendation of Broughton et al. 2024.",
120 default=True,
121 )
122 useAmatrix = pexConfig.Field(
123 dtype=bool,
124 doc="Use the PTC 'a' matrix (Astier et al. 2019 equation 20) "
125 "instead of the average of measured covariances?",
126 default=False,
127 )
129 useCovModelSample = pexConfig.Field(
130 dtype=bool,
131 doc="Use the covariance matrix sampled from the full covariance model "
132 "(Astier et al. 2019 equation 20) instead of the average measured covariances?",
133 default=False,
134 )
136 covModelFluxSample = pexConfig.DictField(
137 keytype=str,
138 itemtype=float,
139 doc="Flux level in electrons at which to sample the full covariance"
140 "model if useCovModelSample=True. The same level is applied to all"
141 "amps if this parameter [`dict`] is passed as {'ALL_AMPS': value}",
142 default={'ALL_AMPS': 25000.0},
143 )
144 maxIterSuccessiveOverRelaxation = pexConfig.Field(
145 dtype=int,
146 doc="The maximum number of iterations allowed for the successive over-relaxation method",
147 default=10000
148 )
149 eLevelSuccessiveOverRelaxation = pexConfig.Field(
150 dtype=float,
151 doc="The target residual error for the successive over-relaxation method",
152 default=5.0e-14
153 )
154 correlationQuadraticFit = pexConfig.Field(
155 dtype=bool,
156 doc="Use a quadratic fit to find the correlations instead of simple averaging?",
157 default=False,
158 )
159 correlationModelRadius = pexConfig.Field(
160 dtype=int,
161 doc="Build a model of the correlation coefficients for radii larger than this value in pixels?",
162 default=100,
163 )
164 correlationModelSlope = pexConfig.Field(
165 dtype=float,
166 doc="Slope of the correlation model for radii larger than correlationModelRadius",
167 default=-1.35,
168 )
170 nSigmaTolForValidEdge = pexConfig.Field(
171 dtype=float,
172 doc="A valid kernel will have all pixels in a 3px picture frame (edges) "
173 "be within +/- nSigmaTolForValidEdge of zero.",
174 default=5.0,
175 )
176 fractionK00ForKernelMax = pexConfig.Field(
177 dtype=float,
178 doc="All kernel values must be less than value*abs(K_00) to be valid.",
179 default=0.025,
180 )
182 useBfkPtc = pexConfig.Field(
183 dtype=bool,
184 doc="Use a BFK ptc in a single pipeline?",
185 default=False,
186 )
188 doCheckValidity = pexConfig.Field(
189 dtype=bool,
190 doc="Check the AMP kernels for basic validity criteria? "
191 "Will set bfk.valid for each amp.",
192 default=True,
193 )
196class BrighterFatterKernelSolveTask(pipeBase.PipelineTask):
197 """Measure appropriate Brighter-Fatter Kernel from the PTC dataset.
198 """
200 ConfigClass = BrighterFatterKernelSolveConfig
201 _DefaultName = 'cpBfkMeasure'
203 def runQuantum(self, butlerQC, inputRefs, outputRefs):
204 """Ensure that the input and output dimensions are passed along.
206 Parameters
207 ----------
208 butlerQC : `lsst.daf.butler.QuantumContext`
209 Butler to operate on.
210 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
211 Input data refs to load.
212 ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection`
213 Output data refs to persist.
214 """
215 inputs = butlerQC.get(inputRefs)
217 # Use the dimensions to set calib/provenance information.
219 if self.config.useBfkPtc:
220 inputs["inputDims"] = dict(inputRefs.inputBfkPtc.dataId.required)
221 inputs["inputPtc"] = inputs["inputBfkPtc"]
222 inputs["dummy"] = []
223 del inputs["inputBfkPtc"]
224 else:
225 inputs["inputDims"] = dict(inputRefs.inputPtc.dataId.required)
227 # Add calibration provenance info to header.
228 kwargs = dict()
229 reference = getattr(inputRefs, "inputPtc", None)
231 if reference is not None and hasattr(reference, "run"):
232 runKey = "PTC_RUN"
233 runValue = reference.run
234 idKey = "PTC_UUID"
235 idValue = str(reference.id)
236 dateKey = "PTC_DATE"
237 calib = inputs.get("inputPtc", None)
238 dateValue = extractCalibDate(calib)
240 kwargs[runKey] = runValue
241 kwargs[idKey] = idValue
242 kwargs[dateKey] = dateValue
244 self.log.info("Using " + str(reference.run))
246 outputs = self.run(**inputs)
247 outputs.outputBFK.updateMetadata(setDate=False, **kwargs)
249 butlerQC.put(outputs, outputRefs)
251 def run(self, inputPtc, dummy, camera, inputDims):
252 """Combine covariance information from PTC into brighter-fatter
253 kernels.
255 Parameters
256 ----------
257 inputPtc : `lsst.ip.isr.PhotonTransferCurveDataset`
258 PTC data containing per-amplifier covariance measurements.
259 dummy : `lsst.afw.image.Exposure`
260 The exposure used to select the appropriate PTC dataset.
261 In almost all circumstances, one of the input exposures
262 used to generate the PTC dataset is the best option.
263 camera : `lsst.afw.cameraGeom.Camera`
264 Camera to use for camera geometry information.
265 inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
266 DataIds to use to populate the output calibration.
268 Returns
269 -------
270 results : `lsst.pipe.base.Struct`
271 The resulst struct containing:
273 ``outputBfk``
274 Resulting Brighter-Fatter Kernel
275 (`lsst.ip.isr.BrighterFatterKernel`).
276 """
277 detector = camera[inputDims['detector']]
278 detName = detector.getName()
280 if self.config.level == 'DETECTOR':
281 detectorCorrList = list()
282 detectorFluxes = list()
284 if not inputPtc.ptcFitType == "FULLCOVARIANCE" and self.config.useCovModelSample:
285 raise ValueError("ptcFitType must be FULLCOVARIANCE if useCovModelSample=True.")
287 # Get flux sample dictionary
288 fluxSampleDict = {ampName: 0.0 for ampName in inputPtc.ampNames}
289 for ampName in inputPtc.ampNames:
290 if 'ALL_AMPS' in self.config.covModelFluxSample:
291 fluxSampleDict[ampName] = self.config.covModelFluxSample['ALL_AMPS']
292 elif ampName in self.config.covModelFluxSample:
293 fluxSampleDict[ampName] = self.config.covModelFluxSample[ampName]
295 bfk = BrighterFatterKernel(camera=camera, detectorId=detector.getId(), level=self.config.level)
296 bfk.rawMeans = inputPtc.rawMeans # ADU
297 bfk.rawVariances = inputPtc.rawVars # ADU^2
298 bfk.expIdMask = inputPtc.expIdMask
300 # Use the PTC covariances as the cross-correlations. These
301 # are scaled before the kernel is generated, which performs
302 # the conversion. The input covariances are in (x, y) index
303 # ordering, as is the aMatrix.
304 bfk.rawXcorrs = inputPtc.covariances # ADU^2
305 bfk.badAmps = inputPtc.badAmps
306 bfk.shape = (inputPtc.covMatrixSide*2 + 1, inputPtc.covMatrixSide*2 + 1)
307 bfk.gain = inputPtc.gain
308 bfk.noise = inputPtc.noise
309 bfk.meanXcorrs = dict()
310 bfk.valid = dict()
311 bfk.updateMetadataFromExposures([inputPtc])
313 for amp in detector:
314 ampName = amp.getName()
315 gain = bfk.gain[ampName]
316 noiseMatrix = inputPtc.noiseMatrix[ampName]
317 mask = inputPtc.expIdMask[ampName]
318 if gain <= 0:
319 # We've received very bad data.
320 self.log.warning("Impossible gain recieved from PTC for %s: %f. Skipping bad amplifier.",
321 ampName, gain)
322 bfk.meanXcorrs[ampName] = np.zeros(bfk.shape)
323 bfk.ampKernels[ampName] = np.zeros(bfk.shape)
324 bfk.rawXcorrs[ampName] = np.zeros((len(mask), inputPtc.covMatrixSide, inputPtc.covMatrixSide))
325 bfk.valid[ampName] = False
326 continue
328 # Use inputPtc.expIdMask to get the means, variances, and
329 # covariances that were not masked after PTC. The
330 # covariances may now have the mask already applied.
331 fluxes = np.array(bfk.rawMeans[ampName])[mask]
332 variances = np.array(bfk.rawVariances[ampName])[mask]
333 covModelList = np.array(inputPtc.covariancesModel[ampName])
335 xCorrList = np.array([np.array(xcorr) for xcorr in bfk.rawXcorrs[ampName]])
336 if np.sum(mask) < len(xCorrList):
337 # Only apply the mask if needed.
338 xCorrList = xCorrList[mask]
340 fluxes = np.array([flux*gain for flux in fluxes]) # Now in e^-
341 variances = np.array([variance*gain*gain for variance in variances]) # Now in e^2-
343 # This should duplicate Coulton et al. 2017 Equation 22-29
344 # (arxiv:1711.06273)
345 scaledCorrList = list()
346 corrList = list()
347 truncatedFluxes = list()
348 for xcorrNum, (xcorr, flux, var) in enumerate(zip(xCorrList, fluxes, variances), 1):
349 q = np.array(xcorr) * gain * gain # xcorr now in e^-
350 q *= 2.0 # Remove factor of 1/2 applied in PTC.
351 self.log.info("Amp: %s %d/%d Flux: %f Var: %f Q(0,0): %g Q(1,0): %g Q(0,1): %g",
352 ampName, xcorrNum, len(xCorrList), flux, var, q[0][0], q[1][0], q[0][1])
354 # Normalize by the flux, which removes the (0,0)
355 # component attributable to Poisson noise. This
356 # contains the two "t I delta(x - x')" terms in
357 # Coulton et al. 2017 equation 29
358 q[0][0] -= 2.0*(flux)
360 if q[0][0] > 0.0:
361 self.log.warning("Amp: %s %d skipped due to value of (variance-mean)=%f",
362 ampName, xcorrNum, q[0][0])
363 # If we drop an element of ``scaledCorrList``
364 # (which is what this does), we need to ensure we
365 # drop the flux entry as well.
366 continue
368 # This removes the "t (I_a^2 + I_b^2)" factor in
369 # Coulton et al. 2017 equation 29.
370 # The quadratic fit option needs the correlations unscaled
371 q /= -2.0
372 unscaled = self._tileArray(q)
373 q /= flux**2
374 scaled = self._tileArray(q)
375 xcorrCheck = np.abs(np.sum(scaled))/np.sum(np.abs(scaled))
376 if (xcorrCheck > self.config.xcorrCheckRejectLevel) or not (np.isfinite(xcorrCheck)):
377 self.log.warning("Amp: %s %d skipped due to value of triangle-inequality sum %f",
378 ampName, xcorrNum, xcorrCheck)
379 continue
381 scaledCorrList.append(scaled)
382 corrList.append(unscaled)
383 truncatedFluxes.append(flux)
384 self.log.info("Amp: %s %d/%d Final: %g XcorrCheck: %f",
385 ampName, xcorrNum, len(xCorrList), q[0][0], xcorrCheck)
387 fluxes = np.array(truncatedFluxes)
389 if len(scaledCorrList) == 0:
390 self.log.warning("Amp: %s All inputs rejected for amp!", ampName)
391 bfk.meanXcorrs[ampName] = np.zeros(bfk.shape)
392 bfk.ampKernels[ampName] = np.zeros(bfk.shape)
393 bfk.valid[ampName] = False
394 continue
396 if self.config.useAmatrix:
397 # Use the aMatrix, ignoring the meanXcorr generated above.
398 preKernel = np.pad(self._tileArray(-1.0 * np.array(inputPtc.aMatrix[ampName])), ((1, 1)))
399 elif self.config.correlationQuadraticFit:
400 # Use a quadratic fit to the correlations as a
401 # function of flux.
402 preKernel = self.quadraticCorrelations(corrList, fluxes, f"Amp: {ampName}")
403 elif self.config.useCovModelSample:
404 # Sample the full covariance model at a given flux.
405 # Use the non-truncated fluxes for this
406 mu = bfk.rawMeans[ampName]
407 covTilde = self.sampleCovModel(mu, noiseMatrix, gain,
408 covModelList, fluxSampleDict[ampName],
409 f"Amp: {ampName}")
410 preKernel = np.pad(self._tileArray(-1.0 * covTilde), ((1, 1)))
411 else:
412 # Use a simple average of the measured correlations.
413 preKernel = self.averageCorrelations(scaledCorrList, f"Amp: {ampName}")
415 center = int((bfk.shape[0] - 1) / 2)
417 if self.config.forceZeroSum:
418 totalSum = np.sum(preKernel)
420 if self.config.correlationModelRadius < (preKernel.shape[0] - 1) / 2:
421 # Assume a correlation model of
422 # Corr(r) = -preFactor * r^(2 * slope)
423 preFactor = np.sqrt(preKernel[center, center + 1] * preKernel[center + 1, center])
424 slopeFactor = 2.0 * np.abs(self.config.correlationModelSlope)
425 totalSum += 2.0*np.pi*(preFactor / (slopeFactor*(center + 0.5))**slopeFactor)
427 preKernel[center, center] -= totalSum
428 self.log.info("%s Zero-Sum Scale: %g", ampName, totalSum)
430 finalSum = np.sum(preKernel)
431 bfk.meanXcorrs[ampName] = preKernel
433 postKernel = self.successiveOverRelax(preKernel)
434 bfk.ampKernels[ampName] = postKernel
435 if self.config.level == 'DETECTOR' and ampName not in self.config.ignoreAmpsForAveraging:
436 detectorCorrList.extend(scaledCorrList)
437 detectorFluxes.extend(fluxes)
439 # Check for validity
440 if self.config.doCheckValidity:
441 # We will check 3 boundary conditions to test the
442 # validitity of the kernel:
443 # (1) The edge of the kernel should be consistent with zero.
444 # (2) The kernel should be negative
445 # (3) The center of the kernal should be the absolute minimum.
447 # Since there are typically too few pixels along the edge, we
448 # will use the stddev of the 3px picture frame of the kernel to
449 # set the tolerances of the validity conditions.
450 kernelSize = bfk.ampKernels[ampName].shape[0]
451 kernelCenterValue = bfk.ampKernels[ampName][kernelSize//2, kernelSize//2]
452 xv, yv = np.meshgrid(range(kernelSize), range(kernelSize))
454 # Get a mask for a 3px picture frame
455 valid = True
456 if kernelSize >= 7:
457 # Last 1px
458 edge = (xv == 0)
459 edge += (xv == kernelSize - 1)
460 edge += (yv == 0)
461 edge += (yv == kernelSize - 1)
463 # Last 3 px
464 edges = (xv < 3)
465 edges += (xv > kernelSize - 3 - 1)
466 edges += (yv < 3)
467 edges += (yv > kernelSize - 3 - 1)
469 kernelEdge = bfk.ampKernels[ampName][edge].ravel()
470 kernelEdges = bfk.ampKernels[ampName][edges].ravel()
471 kernelEdgeStd = np.std(kernelEdges)
473 threshold = self.config.nSigmaTolForValidEdge*kernelEdgeStd
475 # Edge should converge to zero
476 valid = np.all(np.isclose(
477 kernelEdge, 0,
478 rtol=threshold,
479 ))
480 if not valid:
481 self.log.warning("%s: the kernel edges did not converge to 0 within "
482 "+/- %s sigma.", ampName, threshold)
484 # The center should be the minimum
485 valid *= np.all(bfk.ampKernels[ampName] >= kernelCenterValue)
486 if not np.all(bfk.ampKernels[ampName] >= kernelCenterValue):
487 self.log.warning("%s: the kernel center is not the absolute minimum.", ampName)
489 # Kernel should be negative
490 negativityThreshold = self.config.fractionK00ForKernelMax*np.abs(kernelCenterValue)
491 valid *= np.all(bfk.ampKernels[ampName] <= negativityThreshold)
492 if not np.all(bfk.ampKernels[ampName] <= negativityThreshold):
493 self.log.warning("%s: the kernel is not negative", ampName)
494 else:
495 raise RuntimeError("%s: The kernel is too small for validity check.", ampName)
497 bfk.valid[ampName] = valid
498 else:
499 # The kernel at this point will be valid
500 bfk.valid[ampName] = True
502 self.log.info("Amp: %s Sum: %g Center Info Pre: %g Post: %g",
503 ampName, finalSum, preKernel[center, center], postKernel[center, center])
505 # Assemble a detector kernel?
506 if self.config.level == 'DETECTOR':
507 if self.config.correlationQuadraticFit:
508 preKernel = self.quadraticCorrelations(detectorCorrList, detectorFluxes, f"Amp: {ampName}")
509 else:
510 preKernel = self.averageCorrelations(detectorCorrList, f"Det: {detName}")
511 finalSum = np.sum(preKernel)
512 center = int((bfk.shape[0] - 1) / 2)
514 postKernel = self.successiveOverRelax(preKernel)
515 bfk.detKernels[detName] = postKernel
516 self.log.info("Det: %s Sum: %g Center Info Pre: %g Post: %g",
517 detName, finalSum, preKernel[center, center], postKernel[center, center])
519 return pipeBase.Struct(
520 outputBFK=bfk,
521 )
523 def averageCorrelations(self, xCorrList, name):
524 """Average input correlations.
526 Parameters
527 ----------
528 xCorrList : `list` [`numpy.array`]
529 List of cross-correlations. These are expected to be
530 square arrays.
531 name : `str`
532 Name for log messages.
534 Returns
535 -------
536 meanXcorr : `numpy.array`, (N, N)
537 The averaged cross-correlation.
538 """
539 meanXcorr = np.zeros_like(xCorrList[0])
540 xCorrList = np.array(xCorrList)
542 sctrl = afwMath.StatisticsControl()
543 sctrl.setNumSigmaClip(self.config.nSigmaClip)
544 for i in range(np.shape(meanXcorr)[0]):
545 for j in range(np.shape(meanXcorr)[1]):
546 meanXcorr[i, j] = afwMath.makeStatistics(xCorrList[:, i, j],
547 afwMath.MEANCLIP, sctrl).getValue()
549 # To match previous definitions, pad by one element.
550 meanXcorr = np.pad(meanXcorr, ((1, 1)))
552 return meanXcorr
554 def quadraticCorrelations(self, xCorrList, fluxList, name):
555 """Measure a quadratic correlation model.
557 Parameters
558 ----------
559 xCorrList : `list` [`numpy.array`]
560 List of cross-correlations. These are expected to be
561 square arrays.
562 fluxList : `numpy.array`, (Nflux,)
563 Associated list of fluxes.
564 name : `str`
565 Name for log messages.
567 Returns
568 -------
569 meanXcorr : `numpy.array`, (N, N)
570 The averaged cross-correlation.
571 """
572 meanXcorr = np.zeros_like(xCorrList[0])
573 fluxList = np.square(fluxList)
574 xCorrList = np.array(xCorrList)
576 for i in range(np.shape(meanXcorr)[0]):
577 for j in range(np.shape(meanXcorr)[1]):
578 # Fit corrlation_i(x, y) = a0 + a1 * (flux_i)^2 We do
579 # not want to transpose, so use (i, j) without
580 # inversion.
581 linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 1e-4], fluxList,
582 xCorrList[:, i, j], funcPolynomial,
583 scaleResidual=False)
584 meanXcorr[i, j] = linearFit[1] # Discard the intercept.
585 self.log.info("Quad fit meanXcorr[%d,%d] = %g", i, j, linearFit[1])
587 # To match previous definitions, pad by one element.
588 meanXcorr = np.pad(meanXcorr, ((1, 1)))
590 return meanXcorr
592 def sampleCovModel(self, fluxes, noiseMatrix, gain, covModelList, flux, name):
593 """Sample the correlation model and measure
594 widetile{C}_{ij} from Broughton et al. 2023 (eq. 4)
596 Parameters
597 ----------
598 fluxes : `list` [`float`]
599 List of fluxes (in ADU)
600 noiseMatrix : `numpy.array`, (N, N)
601 Noise matrix
602 gain : `float`
603 Amplifier gain
604 covModelList : `numpy.array`, (N, N)
605 List of covariance model matrices. These are
606 expected to be square arrays.
607 flux : `float`
608 Flux in electrons at which to sample the
609 covariance model.
610 name : `str`
611 Name for log messages.
613 Returns
614 -------
615 covTilde : `numpy.array`, (N, N)
616 The calculated C-tilde from Broughton et al. 2023 (eq. 4).
617 """
619 # Get the index of the flux sample
620 # (this must be done in electron units)
621 ix = np.argmin((fluxes*gain - flux)**2)
622 assert len(fluxes) == len(covModelList)
624 # Find the nearest measured flux level
625 # and the full covariance model at that point
626 nearestFlux = fluxes[ix]
627 covModelSample = covModelList[ix]
629 # Calculate flux sample
630 # covTilde returned in ADU units
631 covTilde = (covModelSample - noiseMatrix/gain**2)/(nearestFlux**2)
632 covTilde[0][0] -= (nearestFlux/gain)/(nearestFlux**2)
634 return covTilde
636 @staticmethod
637 def _tileArray(in_array):
638 """Given an input quarter-image, tile/mirror it and return full image.
640 Given a square input of side-length n, of the form
642 input = array([[1, 2, 3],
643 [4, 5, 6],
644 [7, 8, 9]])
646 return an array of size 2n-1 as
648 output = array([[ 9, 8, 7, 8, 9],
649 [ 6, 5, 4, 5, 6],
650 [ 3, 2, 1, 2, 3],
651 [ 6, 5, 4, 5, 6],
652 [ 9, 8, 7, 8, 9]])
654 Parameters
655 ----------
656 input : `np.array`, (N, N)
657 The square input quarter-array
659 Returns
660 -------
661 output : `np.array`, (2*N + 1, 2*N + 1)
662 The full, tiled array
663 """
664 assert in_array.shape[0] == in_array.shape[1]
665 length = in_array.shape[0] - 1
666 output = np.zeros((2*length + 1, 2*length + 1))
668 for i in range(length + 1):
669 for j in range(length + 1):
670 output[i + length, j + length] = in_array[i, j]
671 output[-i + length, j + length] = in_array[i, j]
672 output[i + length, -j + length] = in_array[i, j]
673 output[-i + length, -j + length] = in_array[i, j]
674 return output
676 def successiveOverRelax(self, source, maxIter=None, eLevel=None):
677 """An implementation of the successive over relaxation (SOR) method.
679 A numerical method for solving a system of linear equations
680 with faster convergence than the Gauss-Seidel method.
682 Parameters
683 ----------
684 source : `numpy.ndarray`, (N, N)
685 The input array.
686 maxIter : `int`, optional
687 Maximum number of iterations to attempt before aborting.
688 eLevel : `float`, optional
689 The target error level at which we deem convergence to have
690 occurred.
692 Returns
693 -------
694 output : `numpy.ndarray`, (N, N)
695 The solution.
696 """
697 if not maxIter:
698 maxIter = self.config.maxIterSuccessiveOverRelaxation
699 if not eLevel:
700 eLevel = self.config.eLevelSuccessiveOverRelaxation
702 assert source.shape[0] == source.shape[1], "Input array must be square"
703 # initialize, and set boundary conditions
704 func = np.zeros([source.shape[0] + 2, source.shape[1] + 2])
705 resid = np.zeros([source.shape[0] + 2, source.shape[1] + 2])
706 rhoSpe = np.cos(np.pi/source.shape[0]) # Here a square grid is assumed
708 # Calculate the initial error
709 for i in range(1, func.shape[0] - 1):
710 for j in range(1, func.shape[1] - 1):
711 resid[i, j] = (func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
712 + func[i + 1, j] - 4*func[i, j] - source[i - 1, j - 1])
713 inError = np.sum(np.abs(resid))
715 # Iterate until convergence
716 # We perform two sweeps per cycle,
717 # updating 'odd' and 'even' points separately
718 nIter = 0
719 omega = 1.0
720 dx = 1.0
721 while nIter < maxIter*2:
722 outError = 0
723 if nIter%2 == 0:
724 for i in range(1, func.shape[0] - 1, 2):
725 for j in range(1, func.shape[1] - 1, 2):
726 resid[i, j] = float(func[i, j-1] + func[i, j + 1] + func[i - 1, j]
727 + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
728 func[i, j] += omega*resid[i, j]*.25
729 for i in range(2, func.shape[0] - 1, 2):
730 for j in range(2, func.shape[1] - 1, 2):
731 resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
732 + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
733 func[i, j] += omega*resid[i, j]*.25
734 else:
735 for i in range(1, func.shape[0] - 1, 2):
736 for j in range(2, func.shape[1] - 1, 2):
737 resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
738 + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
739 func[i, j] += omega*resid[i, j]*.25
740 for i in range(2, func.shape[0] - 1, 2):
741 for j in range(1, func.shape[1] - 1, 2):
742 resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
743 + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
744 func[i, j] += omega*resid[i, j]*.25
745 outError = np.sum(np.abs(resid))
746 if outError < inError*eLevel:
747 break
748 if nIter == 0:
749 omega = 1.0/(1 - rhoSpe*rhoSpe/2.0)
750 else:
751 omega = 1.0/(1 - rhoSpe*rhoSpe*omega/4.0)
752 nIter += 1
754 if nIter >= maxIter*2:
755 self.log.warning("Failure: SuccessiveOverRelaxation did not converge in %s iterations."
756 "\noutError: %s, inError: %s,", nIter//2, outError, inError*eLevel)
757 else:
758 self.log.info("Success: SuccessiveOverRelaxation converged in %s iterations."
759 "\noutError: %s, inError: %s", nIter//2, outError, inError*eLevel)
760 return func[1: -1, 1: -1]