Coverage for python/lsst/cp/pipe/utils.py: 7%
763 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-05 01:07 -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#
23__all__ = ['ddict2dict', 'CovFastFourierTransform', 'getReadNoise', 'ampOffsetGainRatioFixup']
25from astropy.table import Table
26import galsim
27import logging
28import numpy as np
29import itertools
30import numpy.polynomial.polynomial as poly
31import warnings
33from scipy.interpolate import Akima1DInterpolator
34from scipy.optimize import leastsq, minimize
35from scipy.stats import median_abs_deviation, norm
37from lsst.ip.isr import isrMock
38import lsst.afw.cameraGeom
39import lsst.afw.image
40import lsst.afw.math
43def sigmaClipCorrection(nSigClip):
44 """Correct measured sigma to account for clipping.
46 If we clip our input data and then measure sigma, then the
47 measured sigma is smaller than the true value because real
48 points beyond the clip threshold have been removed. This is a
49 small (1.5% at nSigClip=3) effect when nSigClip >~ 3, but the
50 default parameters for measure crosstalk use nSigClip=2.0.
51 This causes the measured sigma to be about 15% smaller than
52 real. This formula corrects the issue, for the symmetric case
53 (upper clip threshold equal to lower clip threshold).
55 Parameters
56 ----------
57 nSigClip : `float`
58 Number of sigma the measurement was clipped by.
60 Returns
61 -------
62 scaleFactor : `float`
63 Scale factor to increase the measured sigma by.
64 """
65 varFactor = 1.0 - (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))
66 return 1.0 / np.sqrt(varFactor)
69def calculateWeightedReducedChi2(measured, model, weightsMeasured, nData, nParsModel):
70 """Calculate weighted reduced chi2.
72 Parameters
73 ----------
74 measured : `list`
75 List with measured data.
76 model : `list`
77 List with modeled data.
78 weightsMeasured : `list`
79 List with weights for the measured data.
80 nData : `int`
81 Number of data points.
82 nParsModel : `int`
83 Number of parameters in the model.
85 Returns
86 -------
87 redWeightedChi2 : `float`
88 Reduced weighted chi2.
89 """
90 wRes = (measured - model)*weightsMeasured
91 return ((wRes*wRes).sum())/(nData-nParsModel)
94def makeMockFlats(expTime, gain=1.0, readNoiseElectrons=5, fluxElectrons=1000,
95 randomSeedFlat1=1984, randomSeedFlat2=666, powerLawBfParams=[],
96 expId1=0, expId2=1, ampNames=[]):
97 """Create a pair or mock flats with isrMock.
99 Parameters
100 ----------
101 expTime : `float`
102 Exposure time of the flats.
103 gain : `float`, optional
104 Gain, in e/ADU.
105 readNoiseElectrons : `float`, optional
106 Read noise rms, in electrons.
107 fluxElectrons : `float`, optional
108 Flux of flats, in electrons per second.
109 randomSeedFlat1 : `int`, optional
110 Random seed for the normal distrubutions for the mean signal
111 and noise (flat1).
112 randomSeedFlat2 : `int`, optional
113 Random seed for the normal distrubutions for the mean signal
114 and noise (flat2).
115 powerLawBfParams : `list`, optional
116 Parameters for `galsim.cdmodel.PowerLawCD` to simulate the
117 brightter-fatter effect.
118 expId1 : `int`, optional
119 Exposure ID for first flat.
120 expId2 : `int`, optional
121 Exposure ID for second flat.
122 ampNames : `list` [`str`], optional
123 Names of amplifiers for filling in header information.
125 Returns
126 -------
127 flatExp1 : `lsst.afw.image.exposure.ExposureF`
128 First exposure of flat field pair.
129 flatExp2 : `lsst.afw.image.exposure.ExposureF`
130 Second exposure of flat field pair.
132 Notes
133 -----
134 The parameters of `galsim.cdmodel.PowerLawCD` are `n, r0, t0, rx,
135 tx, r, t, alpha`. For more information about their meaning, see
136 the Galsim documentation
137 https://galsim-developers.github.io/GalSim/_build/html/_modules/galsim/cdmodel.html # noqa: W505
138 and Gruen+15 (1501.02802).
140 Example: galsim.cdmodel.PowerLawCD(8, 1.1e-7, 1.1e-7, 1.0e-8,
141 1.0e-8, 1.0e-9, 1.0e-9, 2.0)
142 """
143 flatFlux = fluxElectrons # e/s
144 flatMean = flatFlux*expTime # e
145 readNoise = readNoiseElectrons # e
147 mockImageConfig = isrMock.IsrMock.ConfigClass()
149 mockImageConfig.flatDrop = 0.99999
150 mockImageConfig.isTrimmed = True
152 flatExp1 = isrMock.FlatMock(config=mockImageConfig).run()
153 flatExp2 = flatExp1.clone()
154 (shapeY, shapeX) = flatExp1.getDimensions()
155 flatWidth = np.sqrt(flatMean)
157 rng1 = np.random.RandomState(randomSeedFlat1)
158 flatData1 = rng1.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng1.normal(0.0, readNoise,
159 (shapeX, shapeY))
160 rng2 = np.random.RandomState(randomSeedFlat2)
161 flatData2 = rng2.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng2.normal(0.0, readNoise,
162 (shapeX, shapeY))
163 # Simulate BF with power law model in galsim
164 if len(powerLawBfParams):
165 if not len(powerLawBfParams) == 8:
166 raise RuntimeError("Wrong number of parameters for `galsim.cdmodel.PowerLawCD`. "
167 f"Expected 8; passed {len(powerLawBfParams)}.")
168 cd = galsim.cdmodel.PowerLawCD(*powerLawBfParams)
169 tempFlatData1 = galsim.Image(flatData1)
170 temp2FlatData1 = cd.applyForward(tempFlatData1)
172 tempFlatData2 = galsim.Image(flatData2)
173 temp2FlatData2 = cd.applyForward(tempFlatData2)
175 flatExp1.image.array[:] = temp2FlatData1.array/gain # ADU
176 flatExp2.image.array[:] = temp2FlatData2.array/gain # ADU
177 else:
178 flatExp1.image.array[:] = flatData1/gain # ADU
179 flatExp2.image.array[:] = flatData2/gain # ADU
181 visitInfoExp1 = lsst.afw.image.VisitInfo(exposureTime=expTime)
182 visitInfoExp2 = lsst.afw.image.VisitInfo(exposureTime=expTime)
184 flatExp1.info.id = expId1
185 flatExp1.getInfo().setVisitInfo(visitInfoExp1)
186 flatExp2.info.id = expId2
187 flatExp2.getInfo().setVisitInfo(visitInfoExp2)
189 # Set ISR metadata.
190 flatExp1.metadata["LSST ISR UNITS"] = "adu"
191 flatExp2.metadata["LSST ISR UNITS"] = "adu"
192 for ampName in ampNames:
193 key = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"
194 value = readNoiseElectrons / gain
196 flatExp1.metadata[key] = value
197 flatExp2.metadata[key] = value
199 key = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
200 flatExp1.metadata[key] = 25000.0 # adu
201 flatExp2.metadata[key] = 25000.0 # adu
203 key = f"LSST ISR AMPOFFSET PEDESTAL {ampName}"
204 value = 0.0
206 flatExp1.metadata[key] = value
207 flatExp2.metadata[key] = value
209 return flatExp1, flatExp2
212def irlsFit(initialParams, dataX, dataY, function, weightsY=None, weightType='Cauchy', scaleResidual=True):
213 """Iteratively reweighted least squares fit.
215 This uses the `lsst.cp.pipe.utils.fitLeastSq`, but applies weights
216 based on the Cauchy distribution by default. Other weight options
217 are implemented. See e.g. Holland and Welsch, 1977,
218 doi:10.1080/03610927708827533
220 Parameters
221 ----------
222 initialParams : `list` [`float`]
223 Starting parameters.
224 dataX : `numpy.array`, (N,)
225 Abscissa data.
226 dataY : `numpy.array`, (N,)
227 Ordinate data.
228 function : callable
229 Function to fit.
230 weightsY : `numpy.array`, (N,)
231 Weights to apply to the data.
232 weightType : `str`, optional
233 Type of weighting to use. One of Cauchy, Anderson, bisquare,
234 box, Welsch, Huber, logistic, or Fair.
235 scaleResidual : `bool`, optional
236 If true, the residual is scaled by the sqrt of the Y values.
238 Returns
239 -------
240 polyFit : `list` [`float`]
241 Final best fit parameters.
242 polyFitErr : `list` [`float`]
243 Final errors on fit parameters.
244 chiSq : `float`
245 Reduced chi squared.
246 weightsY : `list` [`float`]
247 Final weights used for each point.
249 Raises
250 ------
251 RuntimeError :
252 Raised if an unknown weightType string is passed.
253 """
254 if not weightsY:
255 weightsY = np.ones_like(dataX)
257 polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
258 for iteration in range(10):
259 resid = np.abs(dataY - function(polyFit, dataX))
260 if scaleResidual:
261 resid = resid / np.sqrt(dataY)
262 if weightType == 'Cauchy':
263 # Use Cauchy weighting. This is a soft weight.
264 # At [2, 3, 5, 10] sigma, weights are [.59, .39, .19, .05].
265 Z = resid / 2.385
266 weightsY = 1.0 / (1.0 + np.square(Z))
267 elif weightType == 'Anderson':
268 # Anderson+1972 weighting. This is a hard weight.
269 # At [2, 3, 5, 10] sigma, weights are [.67, .35, 0.0, 0.0].
270 Z = resid / (1.339 * np.pi)
271 weightsY = np.where(Z < 1.0, np.sinc(Z), 0.0)
272 elif weightType == 'bisquare':
273 # Beaton and Tukey (1974) biweight. This is a hard weight.
274 # At [2, 3, 5, 10] sigma, weights are [.81, .59, 0.0, 0.0].
275 Z = resid / 4.685
276 weightsY = np.where(Z < 1.0, 1.0 - np.square(Z), 0.0)
277 elif weightType == 'box':
278 # Hinich and Talwar (1975). This is a hard weight.
279 # At [2, 3, 5, 10] sigma, weights are [1.0, 0.0, 0.0, 0.0].
280 weightsY = np.where(resid < 2.795, 1.0, 0.0)
281 elif weightType == 'Welsch':
282 # Dennis and Welsch (1976). This is a hard weight.
283 # At [2, 3, 5, 10] sigma, weights are [.64, .36, .06, 1e-5].
284 Z = resid / 2.985
285 weightsY = np.exp(-1.0 * np.square(Z))
286 elif weightType == 'Huber':
287 # Huber (1964) weighting. This is a soft weight.
288 # At [2, 3, 5, 10] sigma, weights are [.67, .45, .27, .13].
289 Z = resid / 1.345
290 weightsY = np.where(Z < 1.0, 1.0, 1 / Z)
291 elif weightType == 'logistic':
292 # Logistic weighting. This is a soft weight.
293 # At [2, 3, 5, 10] sigma, weights are [.56, .40, .24, .12].
294 Z = resid / 1.205
295 weightsY = np.tanh(Z) / Z
296 elif weightType == 'Fair':
297 # Fair (1974) weighting. This is a soft weight.
298 # At [2, 3, 5, 10] sigma, weights are [.41, .32, .22, .12].
299 Z = resid / 1.4
300 weightsY = (1.0 / (1.0 + (Z)))
301 else:
302 raise RuntimeError(f"Unknown weighting type: {weightType}")
303 polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
305 return polyFit, polyFitErr, chiSq, weightsY
308def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None):
309 """Do a fit and estimate the parameter errors using using
310 scipy.optimize.leastq.
312 optimize.leastsq returns the fractional covariance matrix. To
313 estimate the standard deviation of the fit parameters, multiply
314 the entries of this matrix by the unweighted reduced chi squared
315 and take the square root of the diagonal elements.
317 Parameters
318 ----------
319 initialParams : `list` [`float`]
320 initial values for fit parameters.
321 dataX : `numpy.array`, (N,)
322 Data in the abscissa axis.
323 dataY : `numpy.array`, (N,)
324 Data in the ordinate axis.
325 function : callable object (function)
326 Function to fit the data with.
327 weightsY : `numpy.array`, (N,)
328 Weights of the data in the ordinate axis.
330 Return
331 ------
332 pFitSingleLeastSquares : `list` [`float`]
333 List with fitted parameters.
334 pErrSingleLeastSquares : `list` [`float`]
335 List with errors for fitted parameters.
337 reducedChiSqSingleLeastSquares : `float`
338 Reduced chi squared, unweighted if weightsY is not provided.
339 """
340 if weightsY is None:
341 weightsY = np.ones(len(dataX))
343 def errFunc(p, x, y, weightsY=None):
344 if weightsY is None:
345 weightsY = np.ones(len(x))
346 return (function(p, x) - y)*weightsY
348 pFit, pCov, infoDict, errMessage, success = leastsq(errFunc, initialParams,
349 args=(dataX, dataY, weightsY), full_output=1,
350 epsfcn=0.0001)
352 if (len(dataY) > len(initialParams)) and pCov is not None:
353 reducedChiSq = calculateWeightedReducedChi2(dataY, function(pFit, dataX), weightsY, len(dataY),
354 len(initialParams))
355 pCov *= reducedChiSq
356 else:
357 pCov = np.zeros((len(initialParams), len(initialParams)))
358 pCov[:, :] = np.nan
359 reducedChiSq = np.nan
361 errorVec = []
362 for i in range(len(pFit)):
363 errorVec.append(np.fabs(pCov[i][i])**0.5)
365 pFitSingleLeastSquares = pFit
366 pErrSingleLeastSquares = np.array(errorVec)
368 return pFitSingleLeastSquares, pErrSingleLeastSquares, reducedChiSq
371def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.):
372 """Do a fit using least squares and bootstrap to estimate parameter errors.
374 The bootstrap error bars are calculated by fitting 100 random data sets.
376 Parameters
377 ----------
378 initialParams : `list` [`float`]
379 initial values for fit parameters.
380 dataX : `numpy.array`, (N,)
381 Data in the abscissa axis.
382 dataY : `numpy.array`, (N,)
383 Data in the ordinate axis.
384 function : callable object (function)
385 Function to fit the data with.
386 weightsY : `numpy.array`, (N,), optional.
387 Weights of the data in the ordinate axis.
388 confidenceSigma : `float`, optional.
389 Number of sigmas that determine confidence interval for the
390 bootstrap errors.
392 Return
393 ------
394 pFitBootstrap : `list` [`float`]
395 List with fitted parameters.
396 pErrBootstrap : `list` [`float`]
397 List with errors for fitted parameters.
398 reducedChiSqBootstrap : `float`
399 Reduced chi squared, unweighted if weightsY is not provided.
400 """
401 if weightsY is None:
402 weightsY = np.ones(len(dataX))
404 def errFunc(p, x, y, weightsY):
405 if weightsY is None:
406 weightsY = np.ones(len(x))
407 return (function(p, x) - y)*weightsY
409 # Fit first time
410 pFit, _ = leastsq(errFunc, initialParams, args=(dataX, dataY, weightsY), full_output=0)
412 # Get the stdev of the residuals
413 residuals = errFunc(pFit, dataX, dataY, weightsY)
414 # 100 random data sets are generated and fitted
415 pars = []
416 for i in range(100):
417 randomDelta = np.random.normal(0., np.fabs(residuals), len(dataY))
418 randomDataY = dataY + randomDelta
419 randomFit, _ = leastsq(errFunc, initialParams,
420 args=(dataX, randomDataY, weightsY), full_output=0)
421 pars.append(randomFit)
422 pars = np.array(pars)
423 meanPfit = np.mean(pars, 0)
425 # confidence interval for parameter estimates
426 errPfit = confidenceSigma*np.std(pars, 0)
427 pFitBootstrap = meanPfit
428 pErrBootstrap = errPfit
430 reducedChiSq = calculateWeightedReducedChi2(dataY, function(pFitBootstrap, dataX), weightsY, len(dataY),
431 len(initialParams))
432 return pFitBootstrap, pErrBootstrap, reducedChiSq
435def funcPolynomial(pars, x):
436 """Polynomial function definition
437 Parameters
438 ----------
439 params : `list`
440 Polynomial coefficients. Its length determines the polynomial order.
442 x : `numpy.array`, (N,)
443 Abscisa array.
445 Returns
446 -------
447 y : `numpy.array`, (N,)
448 Ordinate array after evaluating polynomial of order
449 len(pars)-1 at `x`.
450 """
451 return poly.polyval(x, [*pars])
454def funcAstier(pars, x):
455 """Single brighter-fatter parameter model for PTC; Equation 16 of
456 Astier+19.
458 Model:
460 C_{00}(mu) = frac{1}{2 a_{00} g**2} * [exp(2 a_{00} mu g ) - 1]
461 + n_{00} / g**2
463 Parameters
464 ----------
465 params : `list`
466 Parameters of the model: a00 (brightter-fatter), gain (e/ADU),
467 and noise (e^2).
468 x : `numpy.array`, (N,)
469 Signal mu (ADU).
471 Returns
472 -------
473 y : `numpy.array`, (N,)
474 C_00 (variance) in ADU^2.
475 """
476 a00, gain, noiseSquared = pars
477 return 0.5/(a00*gain*gain)*(np.exp(2*a00*x*gain)-1) + noiseSquared/(gain*gain) # C_00
480def ptcRolloffModel(params, mu):
481 """Piece-wise exponential saturation roll-off model of the PTC.
483 Parameters
484 ----------
485 params : `list`
486 Parameters of the model: muTurnoff (adu), tau (rolloff sharpness).
487 mu : `numpy.array`, (N,)
488 Signal mu (ADU).
490 Returns
491 -------
492 y : `numpy.array`, (N,)
493 Difference in variance in ADU^2.
494 """
495 muTurnoff, tau = params
496 return np.where(mu < muTurnoff, 0, np.exp(-(mu - muTurnoff) / tau) - 1)
499def funcAstierWithRolloff(pars, x):
500 """Single brighter-fatter parameter model for PTC; Equation 16 of
501 Astier+19 with an piece-wise exponential model for the PTC roll-off
502 of the PTC caused by saturation.
504 The nominal turnoff is calculated beforehand, and we extend the PTC
505 fit to include signal values up to 5% above the nominally computed
506 turnoff.
508 Model:
510 C_{00}(mu) = funcAstier - np.where(
511 x < muTurnoff,
512 0,
513 np.exp(-(x - muTurnoff) / tau) - 1,
514 )
516 Parameters
517 ----------
518 params : `list`
519 Parameters of the model: a00 (brightter-fatter), gain (e/ADU),
520 and noise (e^2), muTurnoff (adu), tau (rolloff sharpness).
521 x : `numpy.array`, (N,)
522 Signal mu (ADU).
524 Returns
525 -------
526 y : `numpy.array`, (N,)
527 C_00 (variance) in ADU^2.
528 """
529 # Initial computation of Astier+19 (Eqn 19).
530 originalModelPars = pars[:-2]
531 ptcRolloffModelPars = pars[-2:]
532 model = funcAstier(originalModelPars, x)
534 return model - ptcRolloffModel(ptcRolloffModelPars, x) # C_00
537def arrangeFlatsByExpTime(exposureList, exposureIdList, log=None):
538 """Arrange exposures by exposure time.
540 Parameters
541 ----------
542 exposureList : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
543 Input list of exposure references.
544 exposureIdList : `list` [`int`]
545 List of exposure ids as obtained by dataId[`exposure`].
546 log : `lsst.utils.logging.LsstLogAdapter`, optional
547 Log object.
549 Returns
550 ------
551 flatsAtExpTime : `dict` [`float`,
552 `list`[(`lsst.pipe.base.connections.DeferredDatasetRef`,
553 `int`)]]
554 Dictionary that groups references to flat-field exposures
555 (and their IDs) that have the same exposure time (seconds).
556 """
557 flatsAtExpTime = {}
558 assert len(exposureList) == len(exposureIdList), "Different lengths for exp. list and exp. ID lists"
559 for expRef, expId in zip(exposureList, exposureIdList):
560 expTime = expRef.get(component='visitInfo').exposureTime
561 if not np.isfinite(expTime) and log is not None:
562 log.warning("Exposure %d has non-finite exposure time.", expId)
563 listAtExpTime = flatsAtExpTime.setdefault(expTime, [])
564 listAtExpTime.append((expRef, expId))
566 return flatsAtExpTime
569def arrangeFlatsByExpFlux(exposureList, exposureIdList, fluxKeyword, log=None):
570 """Arrange exposures by exposure flux.
572 Parameters
573 ----------
574 exposureList : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
575 Input list of exposure references.
576 exposureIdList : `list` [`int`]
577 List of exposure ids as obtained by dataId[`exposure`].
578 fluxKeyword : `str`
579 Header keyword that contains the flux per exposure.
580 log : `lsst.utils.logging.LsstLogAdapter`, optional
581 Log object.
583 Returns
584 -------
585 flatsAtFlux : `dict` [`float`,
586 `list`[(`lsst.pipe.base.connections.DeferredDatasetRef`,
587 `int`)]]
588 Dictionary that groups references to flat-field exposures
589 (and their IDs) that have the same flux.
590 """
591 flatsAtExpFlux = {}
592 assert len(exposureList) == len(exposureIdList), "Different lengths for exp. list and exp. ID lists"
593 for expRef, expId in zip(exposureList, exposureIdList):
594 # Get flux from header, assuming it is in the metadata.
595 try:
596 expFlux = expRef.get().getMetadata()[fluxKeyword]
597 except KeyError:
598 # If it's missing from the header, continue; it will
599 # be caught and rejected when pairing exposures.
600 expFlux = None
601 if expFlux is None:
602 if log is not None:
603 log.warning("Exposure %d does not have valid header keyword %s.", expId, fluxKeyword)
604 expFlux = np.nan
605 listAtExpFlux = flatsAtExpFlux.setdefault(expFlux, [])
606 listAtExpFlux.append((expRef, expId))
608 return flatsAtExpFlux
611def arrangeFlatsByExpId(exposureList, exposureIdList):
612 """Arrange exposures by exposure ID.
614 There is no guarantee that this will properly group exposures,
615 but allows a sequence of flats that have different
616 illumination (despite having the same exposure time) to be
617 processed.
619 Parameters
620 ----------
621 exposureList : `list`[`lsst.pipe.base.connections.DeferredDatasetRef`]
622 Input list of exposure references.
623 exposureIdList : `list`[`int`]
624 List of exposure ids as obtained by dataId[`exposure`].
626 Returns
627 ------
628 flatsAtExpId : `dict` [`float`,
629 `list`[(`lsst.pipe.base.connections.DeferredDatasetRef`,
630 `int`)]]
631 Dictionary that groups references to flat-field exposures (and their
632 IDs) sequentially by their exposure id.
634 Notes
635 -----
637 This algorithm sorts the input exposure references by their exposure
638 id, and then assigns each pair of exposure references (exp_j, exp_{j+1})
639 to pair k, such that 2*k = j, where j is the python index of one of the
640 exposure references (starting from zero). By checking for the IndexError
641 while appending, we can ensure that there will only ever be fully
642 populated pairs.
643 """
644 flatsAtExpId = {}
645 assert len(exposureList) == len(exposureIdList), "Different lengths for exp. list and exp. ID lists"
646 # Sort exposures by expIds, which are in the second list `exposureIdList`.
647 sortedExposures = sorted(zip(exposureList, exposureIdList), key=lambda pair: pair[1])
649 for jPair, expTuple in enumerate(sortedExposures):
650 if (jPair + 1) % 2:
651 kPair = jPair // 2
652 listAtExpId = flatsAtExpId.setdefault(kPair, [])
653 try:
654 listAtExpId.append(expTuple)
655 listAtExpId.append(sortedExposures[jPair + 1])
656 except IndexError:
657 pass
659 return flatsAtExpId
662def extractCalibDate(calib):
663 """Extract common calibration metadata values that will be written to
664 output header.
666 Parameters
667 ----------
668 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
669 Calibration to pull date information from.
671 Returns
672 -------
673 dateString : `str`
674 Calibration creation date string to add to header.
675 """
676 if hasattr(calib, "getMetadata"):
677 if 'CALIB_CREATION_DATE' in calib.getMetadata():
678 return " ".join((calib.getMetadata().get("CALIB_CREATION_DATE", "Unknown"),
679 calib.getMetadata().get("CALIB_CREATION_TIME", "Unknown")))
680 else:
681 return " ".join((calib.getMetadata().get("CALIB_CREATE_DATE", "Unknown"),
682 calib.getMetadata().get("CALIB_CREATE_TIME", "Unknown")))
683 else:
684 return "Unknown Unknown"
687class CovFastFourierTransform:
688 """A class to compute (via FFT) the nearby pixels correlation function.
690 Implements appendix of Astier+19.
692 Parameters
693 ----------
694 diff : `numpy.array`
695 Image where to calculate the covariances (e.g., the difference
696 image of two flats).
697 w : `numpy.array`
698 Weight image (mask): it should consist of 1's (good pixel) and
699 0's (bad pixels).
700 fftShape : `tuple`
701 2d-tuple with the shape of the FFT
702 maxRangeCov : `int`
703 Maximum range for the covariances.
704 """
706 def __init__(self, diff, w, fftShape, maxRangeCov):
707 # check that the zero padding implied by "fft_shape"
708 # is large enough for the required correlation range
709 assert fftShape[0] > diff.shape[0]+maxRangeCov+1
710 assert fftShape[1] > diff.shape[1]+maxRangeCov+1
711 # for some reason related to numpy.fft.rfftn,
712 # the second dimension should be even, so
713 if fftShape[1]%2 == 1:
714 fftShape = (fftShape[0], fftShape[1]+1)
715 tIm = np.fft.rfft2(diff*w, fftShape)
716 tMask = np.fft.rfft2(w, fftShape)
717 # sum of "squares"
718 self.pCov = np.fft.irfft2(tIm*tIm.conjugate())
719 # sum of values
720 self.pMean = np.fft.irfft2(tIm*tMask.conjugate())
721 # number of w!=0 pixels.
722 self.pCount = np.fft.irfft2(tMask*tMask.conjugate())
724 def cov(self, dx, dy):
725 """Covariance for dx,dy averaged with dx,-dy if both non zero.
727 Implements appendix of Astier+19.
729 Parameters
730 ----------
731 dx : `int`
732 Lag in x
733 dy : `int`
734 Lag in y
736 Returns
737 -------
738 0.5*(cov1+cov2) : `float`
739 Covariance at (dx, dy) lag
740 npix1+npix2 : `int`
741 Number of pixels used in covariance calculation.
743 Raises
744 ------
745 ValueError if number of pixels for a given lag is 0.
746 """
747 # compensate rounding errors
748 nPix1 = int(round(self.pCount[dy, dx]))
749 if nPix1 == 0:
750 raise ValueError(f"Could not compute covariance term {dy}, {dx}, as there are no good pixels.")
751 cov1 = self.pCov[dy, dx]/nPix1-self.pMean[dy, dx]*self.pMean[-dy, -dx]/(nPix1*nPix1)
752 if (dx == 0 or dy == 0):
753 return cov1, nPix1
754 nPix2 = int(round(self.pCount[-dy, dx]))
755 if nPix2 == 0:
756 raise ValueError("Could not compute covariance term {dy}, {dx} as there are no good pixels.")
757 cov2 = self.pCov[-dy, dx]/nPix2-self.pMean[-dy, dx]*self.pMean[dy, -dx]/(nPix2*nPix2)
758 return 0.5*(cov1+cov2), nPix1+nPix2
760 def reportCovFastFourierTransform(self, maxRange):
761 """Produce a list of tuples with covariances.
763 Implements appendix of Astier+19.
765 Parameters
766 ----------
767 maxRange : `int`
768 Maximum range of covariances.
770 Returns
771 -------
772 tupleVec : `list`
773 List with covariance tuples.
774 """
775 tupleVec = []
776 # (dy,dx) = (0,0) has to be first
777 for dy in range(maxRange+1):
778 for dx in range(maxRange+1):
779 cov, npix = self.cov(dx, dy)
780 if (dx == 0 and dy == 0):
781 var = cov
782 tupleVec.append((dx, dy, var, cov, npix))
783 return tupleVec
786def getFitDataFromCovariances(i, j, mu, fullCov, fullCovModel, fullCovSqrtWeights, gain=1.0,
787 divideByMu=False, returnMasked=False):
788 """Get measured signal and covariance, cov model, weigths, and mask at
789 covariance lag (i, j).
791 Parameters
792 ----------
793 i : `int`
794 Lag for covariance matrix.
795 j : `int`
796 Lag for covariance matrix.
797 mu : `list`
798 Mean signal values.
799 fullCov : `list` of `numpy.array`
800 Measured covariance matrices at each mean signal level in mu.
801 fullCovSqrtWeights : `list` of `numpy.array`
802 List of square root of measured covariances at each mean
803 signal level in mu.
804 fullCovModel : `list` of `numpy.array`
805 List of modeled covariances at each mean signal level in mu.
806 gain : `float`, optional
807 Gain, in e-/ADU. If other than 1.0 (default), the returned
808 quantities will be in electrons or powers of electrons.
809 divideByMu : `bool`, optional
810 Divide returned covariance, model, and weights by the mean
811 signal mu?
812 returnMasked : `bool`, optional
813 Use mask (based on weights) in returned arrays (mu,
814 covariance, and model)?
816 Returns
817 -------
818 mu : `numpy.array`
819 list of signal values at (i, j).
820 covariance : `numpy.array`
821 Covariance at (i, j) at each mean signal mu value (fullCov[:, i, j]).
822 covarianceModel : `numpy.array`
823 Covariance model at (i, j).
824 weights : `numpy.array`
825 Weights at (i, j).
826 maskFromWeights : `numpy.array`, optional
827 Boolean mask of the covariance at (i,j), where the weights
828 differ from 0.
829 """
830 mu = np.array(mu)
831 fullCov = np.array(fullCov)
832 fullCovModel = np.array(fullCovModel)
833 fullCovSqrtWeights = np.array(fullCovSqrtWeights)
834 covariance = fullCov[:, i, j]*(gain**2)
835 covarianceModel = fullCovModel[:, i, j]*(gain**2)
836 weights = fullCovSqrtWeights[:, i, j]/(gain**2)
838 maskFromWeights = weights != 0
839 if returnMasked:
840 weights = weights[maskFromWeights]
841 covarianceModel = covarianceModel[maskFromWeights]
842 mu = mu[maskFromWeights]
843 covariance = covariance[maskFromWeights]
845 if divideByMu:
846 covariance /= mu
847 covarianceModel /= mu
848 weights *= mu
849 return mu, covariance, covarianceModel, weights, maskFromWeights
852def symmetrize(inputArray):
853 """ Copy array over 4 quadrants prior to convolution.
855 Parameters
856 ----------
857 inputarray : `numpy.array`
858 Input array to symmetrize.
860 Returns
861 -------
862 aSym : `numpy.array`
863 Symmetrized array.
864 """
865 targetShape = list(inputArray.shape)
866 r1, r2 = inputArray.shape[-1], inputArray.shape[-2]
867 targetShape[-1] = 2*r1-1
868 targetShape[-2] = 2*r2-1
869 aSym = np.ndarray(tuple(targetShape))
870 aSym[..., r2-1:, r1-1:] = inputArray
871 aSym[..., r2-1:, r1-1::-1] = inputArray
872 aSym[..., r2-1::-1, r1-1::-1] = inputArray
873 aSym[..., r2-1::-1, r1-1:] = inputArray
875 return aSym
878def ddict2dict(d):
879 """Convert nested default dictionaries to regular dictionaries.
881 This is needed to prevent yaml persistence issues.
883 Parameters
884 ----------
885 d : `defaultdict`
886 A possibly nested set of `defaultdict`.
888 Returns
889 -------
890 dict : `dict`
891 A possibly nested set of `dict`.
892 """
893 for k, v in d.items():
894 if isinstance(v, dict):
895 d[k] = ddict2dict(v)
896 return dict(d)
899class Pol2D:
900 """2D Polynomial Regression.
902 Parameters
903 ----------
904 x : numpy.ndarray
905 Input array for the x-coordinate.
906 y : numpy.ndarray
907 Input array for the y-coordinate.
908 z : numpy.ndarray
909 Input array for the dependent variable.
910 order : int
911 Order of the polynomial.
912 w : numpy.ndarray, optional
913 Weight array for weighted regression. Default is None.
915 Notes
916 -----
917 Ported from by https://gitlab.in2p3.fr/astier/bfptc P. Astier.
919 Example:
920 >>> x = np.array([1, 2, 3])
921 >>> y = np.array([4, 5, 6])
922 >>> z = np.array([7, 8, 9])
923 >>> order = 2
924 >>> poly_reg = Pol2D(x, y, z, order)
925 >>> result = poly_reg.eval(2.5, 5.5)
926 """
927 def __init__(self, x, y, z, order, w=None):
928 """
929 orderx : `int`
930 Effective order in the x-direction.
931 ordery : `int`
932 Effective order in the y-direction.
933 coeff : `numpy.ndarray`
934 Coefficients of the polynomial regression.
935 """
936 self.orderx = min(order, x.shape[0] - 1)
937 self.ordery = min(order, x.shape[1] - 1)
938 G = self.monomials(x.ravel(), y.ravel())
939 if w is None:
940 self.coeff, _, rank, _ = np.linalg.lstsq(G, z.ravel(), rcond=None)
941 else:
942 self.coeff, _, rank, _ = np.linalg.lstsq((w.ravel() * G.T).T, z.ravel() * w.ravel(), rcond=None)
944 def monomials(self, x, y):
945 """
946 Generate the monomials matrix for the given x and y.
948 Parameters
949 ----------
950 x : numpy.ndarray
951 Input array for the x-coordinate.
952 y : numpy.ndarray
953 Input array for the y-coordinate.
955 Returns
956 -------
957 G : numpy.ndarray
958 Monomials matrix.
959 """
960 ncols = (self.orderx + 1) * (self.ordery + 1)
961 G = np.zeros(x.shape + (ncols,))
962 ij = itertools.product(range(self.orderx + 1), range(self.ordery + 1))
963 for k, (i, j) in enumerate(ij):
964 G[..., k] = x**i * y**j
965 return G
967 def eval(self, x, y):
968 """
969 Evaluate the polynomial at the given x and y coordinates.
971 Parameters
972 ----------
973 x : `float`
974 x-coordinate for evaluation.
975 y : `float`
976 y-coordinate for evaluation.
978 Returns
979 -------
980 result : `float`
981 Result of the polynomial evaluation.
982 """
983 G = self.monomials(x, y)
984 return np.dot(G, self.coeff)
987class AstierSplineLinearityFitter:
988 """Class to fit the Astier spline linearity model.
990 This is a spline fit with photodiode data based on a model
991 from Pierre Astier, referenced in June 2023 from
992 https://me.lsst.eu/astier/bot/7224D/model_nonlin.py
994 This model fits a spline with (optional) nuisance parameters
995 to allow for different linearity coefficients with different
996 photodiode settings. The minimization is a least-squares
997 fit with the residual of
998 Sum[(S(mu_i) + mu_i)/(k_j * D_i - O) - 1]**2, where S(mu_i)
999 is an Akima Spline function of mu_i, the observed flat-pair
1000 mean; D_j is the photo-diode measurement corresponding to
1001 that flat-pair; and k_j is a constant of proportionality
1002 which is over index j as it is allowed to
1003 be different based on different photodiode settings (e.g.
1004 CCOBCURR); and O is a constant offset to allow for light
1005 leaks (and is only fit if fit_offset=True). In
1006 addition, if config.doSplineFitTemperature is True then
1007 the fit will adjust mu such that
1008 mu = mu_input*(1 + alpha*(temperature_scaled))
1009 and temperature_scaled = T - T_ref. Finally, if
1010 config.doSplineFitTemporal is True then the fit will
1011 further adjust mu such that
1012 mu = mu_input*(1 + beta*(mjd_scaled))
1013 and mjd_scaled = mjd - mjd_ref.
1015 The fit has additional constraints to ensure that the spline
1016 goes through the (0, 0) point, as well as a normalization
1017 condition so that the average of the spline over the full
1018 range is 0. The normalization ensures that the spline only
1019 fits deviations from linearity, rather than the linear
1020 function itself which is degenerate with the gain.
1022 Parameters
1023 ----------
1024 nodes : `np.ndarray` (N,)
1025 Array of spline node locations.
1026 grouping_values : `np.ndarray` (M,)
1027 Array of values to group values for different proportionality
1028 constants (e.g. CCOBCURR).
1029 pd : `np.ndarray` (M,)
1030 Array of photodiode measurements.
1031 mu : `np.ndarray` (M,)
1032 Array of flat mean values.
1033 mask : `np.ndarray` (M,), optional
1034 Input mask (True is good point, False is bad point).
1035 log : `logging.logger`, optional
1036 Logger object to use for logging.
1037 fit_offset : `bool`, optional
1038 Fit a constant offset to allow for light leaks?
1039 fit_weights : `bool`, optional
1040 Fit for the weight parameters?
1041 weight_pars_start : `list` [ `float` ]
1042 Iterable of 2 weight parameters for weighed fit. These will
1043 be used as input if the weight parameters are not fit.
1044 fit_temperature : `bool`, optional
1045 Fit for temperature scaling?
1046 temperature_scaled : `np.ndarray` (M,), optional
1047 Input scaled temperature values (T - T_ref).
1048 max_signal_nearly_linear : `float`, optional
1049 Maximum signal that we are confident the input is nearly
1050 linear. This is used both for regularization, and for
1051 fitting the raw slope. Usually set to the ptc turnoff,
1052 above which we allow the spline to significantly deviate
1053 and do not demand the deviation to average to zero.
1054 fit_temporal : `bool`, optional
1055 Fit for temporal scaling?
1056 mjd_scaled : `np.ndarray` (M,), optional
1057 Input scaled mjd values (mjd - mjd_ref).
1058 max_correction : `float`, optional
1059 Maximum fractional correction.
1060 """
1061 def __init__(
1062 self,
1063 nodes,
1064 grouping_values,
1065 pd,
1066 mu,
1067 mask=None,
1068 log=None,
1069 fit_offset=True,
1070 fit_weights=False,
1071 weight_pars_start=[1.0, 0.0],
1072 fit_temperature=False,
1073 temperature_scaled=None,
1074 max_signal_nearly_linear=None,
1075 fit_temporal=False,
1076 mjd_scaled=None,
1077 max_correction=0.25,
1078 ):
1079 self._pd = pd
1080 self._mu = mu
1081 self._grouping_values = grouping_values
1082 self.log = log if log else logging.getLogger(__name__)
1083 self._fit_offset = fit_offset
1084 self._fit_weights = fit_weights
1085 self._weight_pars_start = weight_pars_start
1086 self._fit_temperature = fit_temperature
1087 self._fit_temporal = fit_temporal
1088 self._max_correction = max_correction
1090 self._nodes = nodes
1091 if nodes[0] != 0.0:
1092 raise ValueError("First node must be 0.0")
1093 if not np.all(np.diff(nodes) > 0):
1094 raise ValueError("Nodes must be sorted with no repeats.")
1096 # Find the group indices.
1097 u_group_values = np.unique(self._grouping_values)
1098 self.ngroup = len(u_group_values)
1100 self.group_indices = []
1101 for i in range(self.ngroup):
1102 self.group_indices.append(np.where(self._grouping_values == u_group_values[i])[0])
1104 # Weight values. Outliers will be set to 0.
1105 if mask is None:
1106 _mask = np.ones(len(mu), dtype=np.bool_)
1107 else:
1108 _mask = mask
1109 self._w = self.compute_weights(self._weight_pars_start, self._mu, _mask)
1111 if temperature_scaled is None:
1112 temperature_scaled = np.zeros(len(self._mu))
1113 else:
1114 if len(np.atleast_1d(temperature_scaled)) != len(self._mu):
1115 raise ValueError("temperature_scaled must be the same length as input mu.")
1116 self._temperature_scaled = temperature_scaled
1118 if mjd_scaled is None:
1119 mjd_scaled = np.zeros(len(self._mu))
1120 else:
1121 if len(np.atleast_1d(mjd_scaled)) != len(self._mu):
1122 raise ValueError("mjd_scaled must be the same length as input mu.")
1123 self._mjd_scaled = mjd_scaled
1125 # Values to regularize spline fit.
1126 if max_signal_nearly_linear is None:
1127 max_signal_nearly_linear = self._mu[self.mask].max()
1128 self._max_signal_nearly_linear = max_signal_nearly_linear
1129 self._x_regularize = np.linspace(0.0, self._max_signal_nearly_linear, 100)
1131 # Set up the indices for the fit parameters.
1132 self.par_indices = {
1133 "values": np.arange(len(self._nodes)),
1134 "groups": len(self._nodes) + np.arange(self.ngroup),
1135 "offset": np.zeros(0, dtype=np.int64),
1136 "weight_pars": np.zeros(0, dtype=np.int64),
1137 "temperature_coeff": np.zeros(0, dtype=np.int64),
1138 "temporal_coeff": np.zeros(0, dtype=np.int64),
1139 }
1140 if self._fit_offset:
1141 self.par_indices["offset"] = np.arange(1) + (
1142 len(self.par_indices["values"])
1143 + len(self.par_indices["groups"])
1144 )
1145 if self._fit_weights:
1146 self.par_indices["weight_pars"] = np.arange(2) + (
1147 len(self.par_indices["values"])
1148 + len(self.par_indices["groups"])
1149 + len(self.par_indices["offset"])
1150 )
1151 if self._fit_temperature:
1152 self.par_indices["temperature_coeff"] = np.arange(1) + (
1153 len(self.par_indices["values"])
1154 + len(self.par_indices["groups"])
1155 + len(self.par_indices["offset"])
1156 + len(self.par_indices["weight_pars"])
1157 )
1158 if self._fit_temporal:
1159 self.par_indices["temporal_coeff"] = np.arange(1) + (
1160 len(self.par_indices["values"])
1161 + len(self.par_indices["groups"])
1162 + len(self.par_indices["offset"])
1163 + len(self.par_indices["weight_pars"])
1164 + len(self.par_indices["temperature_coeff"])
1165 )
1167 @staticmethod
1168 def compute_weights(weight_pars, mu, mask):
1169 w = 1./np.sqrt(weight_pars[0]**2. + weight_pars[1]**2./mu)
1170 w[~mask] = 0.0
1172 return w
1174 def estimate_p0(self):
1175 """Estimate initial fit parameters.
1177 Returns
1178 -------
1179 p0 : `np.ndarray`
1180 Parameter array, with spline values (one for each node) followed
1181 by proportionality constants (one for each group); one extra
1182 for the offset O (if fit_offset was set to True); two extra
1183 for the weights (if fit_weights was set to True); one
1184 extra for the temperature coefficient (if fit_temperature was
1185 set to True); and one extra for the temporal coefficient (if
1186 fit_temporal was set to True).
1187 """
1188 npt = (len(self.par_indices["values"])
1189 + len(self.par_indices["groups"])
1190 + len(self.par_indices["offset"])
1191 + len(self.par_indices["weight_pars"])
1192 + len(self.par_indices["temperature_coeff"])
1193 + len(self.par_indices["temporal_coeff"]))
1194 p0 = np.zeros(npt)
1196 # Do a simple linear fit for each group.
1197 for i, indices in enumerate(self.group_indices):
1198 mask = self.mask[indices]
1199 mu = self._mu[indices][mask]
1200 pd = self._pd[indices][mask]
1201 to_fit = (mu < self._max_signal_nearly_linear)
1202 linfit = np.polyfit(pd[to_fit], mu[to_fit], 1)
1203 p0[self.par_indices["groups"][i]] = linfit[0]
1205 # Look at the residuals...
1206 ratio_model = self.compute_ratio_model(
1207 self._nodes,
1208 self.group_indices,
1209 self.par_indices,
1210 p0,
1211 self._pd,
1212 self._mu,
1213 self._temperature_scaled,
1214 self._mjd_scaled,
1215 )
1216 # ...and adjust the linear parameters accordingly.
1217 p0[self.par_indices["groups"]] *= np.median(ratio_model[self.mask])
1219 # Re-compute the residuals.
1220 ratio_model2 = self.compute_ratio_model(
1221 self._nodes,
1222 self.group_indices,
1223 self.par_indices,
1224 p0,
1225 self._pd,
1226 self._mu,
1227 self._temperature_scaled,
1228 self._mjd_scaled,
1229 )
1231 # And compute a first guess of the spline nodes.
1232 bins = np.clip(np.searchsorted(self._nodes, self._mu[self.mask]), 0, len(self._nodes) - 1)
1233 tot_arr = np.zeros(len(self._nodes))
1234 n_arr = np.zeros(len(self._nodes), dtype=int)
1235 np.add.at(tot_arr, bins, ratio_model2[self.mask])
1236 np.add.at(n_arr, bins, 1)
1238 ratio = np.ones(len(self._nodes))
1239 ratio[n_arr > 0] = tot_arr[n_arr > 0]/n_arr[n_arr > 0]
1240 ratio[0] = 1.0
1241 p0[self.par_indices["values"]] = (ratio - 1) * self._nodes
1243 if self._fit_offset:
1244 p0[self.par_indices["offset"]] = 0.0
1246 if self._fit_weights:
1247 p0[self.par_indices["weight_pars"]] = self._weight_pars_start
1249 return p0
1251 @staticmethod
1252 def compute_ratio_model(
1253 nodes,
1254 group_indices,
1255 par_indices,
1256 pars,
1257 pd,
1258 mu,
1259 temperature_scaled,
1260 mjd_scaled,
1261 return_spline=False,
1262 ):
1263 """Compute the ratio model values.
1265 Parameters
1266 ----------
1267 nodes : `np.ndarray` (M,)
1268 Array of node positions.
1269 group_indices : `list` [`np.ndarray`]
1270 List of group indices, one array for each group.
1271 par_indices : `dict`
1272 Dictionary showing which indices in the pars belong to
1273 each set of fit values.
1274 pars : `np.ndarray`
1275 Parameter array, with spline values (one for each node) followed
1276 by proportionality constants (one for each group); one extra
1277 for the offset O (if fit_offset was set to True); two extra
1278 for the weights (if fit_weights was set to True); and one
1279 extra for the temperature coefficient (if fit_temperature was
1280 set to True).
1281 pd : `np.ndarray` (N,)
1282 Array of photodiode measurements.
1283 mu : `np.ndarray` (N,)
1284 Array of flat means.
1285 temperature_scaled : `np.ndarray` (N,)
1286 Array of scaled temperature values.
1287 mjd_scaled : `np.ndarray` (N,)
1288 Array of scaled mjd values.
1289 return_spline : `bool`, optional
1290 Return the spline interpolation as well as the model ratios?
1292 Returns
1293 -------
1294 ratio_models : `np.ndarray` (N,)
1295 Model ratio, (mu_i - S(mu_i) - O)/(k_j * D_i)
1296 spl : `scipy.interpolate.Akima1DInterpolator`
1297 Spline interpolator (returned if return_spline=True).
1298 """
1299 spl = Akima1DInterpolator(nodes, pars[par_indices["values"]], method="akima")
1301 # Check if we want to do just the left or both with temp scale.
1302 if len(par_indices["temperature_coeff"]) == 1:
1303 mu_corr = mu*(1. + pars[par_indices["temperature_coeff"]]*temperature_scaled)
1304 else:
1305 mu_corr = mu
1307 if len(par_indices["temporal_coeff"]) == 1:
1308 mu_corr = mu_corr*(1. + pars[par_indices["temporal_coeff"]]*mjd_scaled)
1310 numerator = mu_corr - spl(np.clip(mu_corr, nodes[0], nodes[-1]))
1311 if len(par_indices["offset"]) == 1:
1312 numerator -= pars[par_indices["offset"][0]]
1313 denominator = pd.copy()
1314 ngroup = len(group_indices)
1315 kj = pars[par_indices["groups"]]
1316 for j in range(ngroup):
1317 denominator[group_indices[j]] *= kj[j]
1319 if return_spline:
1320 return numerator / denominator, spl
1321 else:
1322 return numerator / denominator
1324 def fit(self, p0, min_iter=3, max_iter=20, max_rejection_per_iteration=5, n_sigma_clip=5.0):
1325 """
1326 Perform iterative fit for linear + spline model with offsets.
1328 Parameters
1329 ----------
1330 p0 : `np.ndarray`
1331 Initial parameter array, with spline values (one for each node)
1332 followed by proportionality constants (one for each group); one
1333 extra for the offset O (if fit_offset was set to True); two extra
1334 for the weights (if fit_weights was set to True); one
1335 extra for the temperature coefficient (if fit_temperature was
1336 set to True); and one extra for the temporal coefficient (if
1337 fit_temporal was set to True).
1338 min_iter : `int`, optional
1339 Minimum number of fit iterations.
1340 max_iter : `int`, optional
1341 Maximum number of fit iterations.
1342 max_rejection_per_iteration : `int`, optional
1343 Maximum number of points to reject per iteration.
1344 n_sigma_clip : `float`, optional
1345 Number of sigma to do clipping in each iteration.
1346 """
1347 init_params = p0
1348 for k in range(max_iter):
1349 params, cov_params, _, msg, ierr = leastsq(
1350 self,
1351 init_params,
1352 full_output=True,
1353 ftol=1e-5,
1354 maxfev=12000,
1355 )
1356 init_params = params.copy()
1358 # We need to cut off the constraints at the end (there are more
1359 # residuals than data points.)
1360 res = self(params)[: len(self._w)]
1361 std_res = median_abs_deviation(res[self.good_points], scale="normal")
1362 sample = len(self.good_points)
1364 # We don't want to reject too many outliers at once.
1365 if sample > max_rejection_per_iteration:
1366 sres = np.sort(np.abs(res))
1367 cut = max(sres[-max_rejection_per_iteration], std_res*n_sigma_clip)
1368 else:
1369 cut = std_res*n_sigma_clip
1371 outliers = np.abs(res) > cut
1372 self._w[outliers] = 0
1373 if outliers.sum() != 0:
1374 self.log.info(
1375 "After iteration %d there are %d outliers (of %d).",
1376 k,
1377 outliers.sum(),
1378 sample,
1379 )
1380 elif k >= min_iter:
1381 self.log.info("After iteration %d there are no more outliers.", k)
1382 break
1384 return params
1386 @property
1387 def mask(self):
1388 return (self._w > 0)
1390 @property
1391 def good_points(self):
1392 return self.mask.nonzero()[0]
1394 def compute_chisq_dof(self, pars):
1395 """Compute the chi-squared per degree of freedom for a set of pars.
1397 Parameters
1398 ----------
1399 pars : `np.ndarray`
1400 Parameter array.
1402 Returns
1403 -------
1404 chisq_dof : `float`
1405 Chi-squared per degree of freedom.
1406 """
1407 resids = self(pars)[0: len(self.mask)]
1408 chisq = np.sum(resids[self.mask]**2.)
1409 dof = self.mask.sum() - self.ngroup
1410 if self._fit_temporal:
1411 dof -= 1
1412 if self._fit_temperature:
1413 dof -= 1
1414 if self._fit_offset:
1415 dof -= 1
1416 if self._fit_weights:
1417 dof -= 2
1419 return chisq/dof
1421 def __call__(self, pars):
1423 ratio_model, spl = self.compute_ratio_model(
1424 self._nodes,
1425 self.group_indices,
1426 self.par_indices,
1427 pars,
1428 self._pd,
1429 self._mu,
1430 self._temperature_scaled,
1431 self._mjd_scaled,
1432 return_spline=True,
1433 )
1435 _mask = self.mask
1436 # Update the weights if we are fitting them.
1437 if self._fit_weights:
1438 self._w = self.compute_weights(pars[self.par_indices["weight_pars"]], self._mu, _mask)
1439 resid = self._w*(ratio_model - 1.0)
1441 # Ensure masked points have 0 residual.
1442 resid[~_mask] = 0.0
1444 constraint = [1e3 * np.mean(spl(np.clip(self._x_regularize, self._nodes[0], self._nodes[-1])))]
1445 # 0 should transform to 0
1446 constraint.append(spl(0)*1e10)
1447 # Use a Jeffreys prior on the weight if we are fitting it.
1448 if self._fit_weights:
1449 # This factor ensures that log(fact * w) is negative.
1450 fact = 1e-3 / self._w.max()
1451 # We only add non-zero weights to the constraint array.
1452 log_w = np.sqrt(-2.*np.log(fact*self._w[self._w > 0]))
1453 constraint = np.hstack([constraint, log_w])
1455 # Don't let it get to >5% correction.
1456 values = pars[self.par_indices["values"]]
1457 if np.abs(values[-1])/self._nodes[-1] > self._max_correction:
1458 extra_constraint = 1e10
1459 else:
1460 extra_constraint = 0
1462 return np.hstack([resid, constraint, extra_constraint])
1465def getReadNoise(exposure, ampName, taskMetadata=None, log=None):
1466 """Gets readout noise for an amp from ISR metadata.
1468 If possible, this attempts to get the now-standard headers
1469 added to the exposure itself. If not found there, the ISR
1470 TaskMetadata is searched. If neither of these has the value,
1471 warn and set the read noise to NaN.
1473 Parameters
1474 ----------
1475 exposure : `lsst.afw.image.Exposure`
1476 Exposure to check for read noise first.
1477 ampName : `str`
1478 Amplifier name.
1479 taskMetadata : `lsst.pipe.base.TaskMetadata`, optional
1480 List of exposures metadata from ISR for this exposure.
1481 log : `logging.logger`, optional
1482 Log for messages.
1484 Returns
1485 -------
1486 readNoise : `float`
1487 The read noise for this set of exposure/amplifier.
1488 """
1489 exposureMetadata = exposure.getMetadata()
1491 # Try from the exposure first.
1492 expectedKey = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"
1493 if expectedKey in exposureMetadata:
1494 return exposureMetadata[expectedKey]
1496 # If not, try getting it from the task metadata.
1497 if taskMetadata:
1498 expectedKey = f"RESIDUAL STDEV {ampName}"
1499 if "isr" in taskMetadata:
1500 if expectedKey in taskMetadata["isr"]:
1501 return taskMetadata["isr"][expectedKey]
1503 log = log if log else logging.getLogger(__name__)
1504 log.warning("Median readout noise from ISR metadata for amp %s "
1505 "could not be found.", ampName)
1506 return np.nan
1509def ampOffsetGainRatioFixup(ptc, minAdu, maxAdu, log=None):
1510 """Apply gain ratio fixup using amp offsets.
1512 Parameters
1513 ----------
1514 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
1515 Input PTC. Will be modified in place.
1516 minAdu : `float`
1517 Minimum number of ADU in mean of amplifier to use for
1518 computing gain ratio fixup.
1519 maxAdu : `float`
1520 Maximum number of ADU in mean of amplifier to use for
1521 computing gain ratio fixup.
1522 log : `lsst.utils.logging.LsstLogAdapter`, optional
1523 Log object.
1524 """
1525 # We need to find the reference amplifier.
1526 # Find an amp near the middle to use as a pivot.
1527 if log is None:
1528 log = logging.getLogger(__name__)
1530 # We first check for which amps have gain measurements
1531 # (fully bad amps are filled with NaN for gain.)
1532 gainArray = np.zeros(len(ptc.ampNames))
1533 for i, ampName in enumerate(ptc.ampNames):
1534 gainArray[i] = ptc.gain[ampName]
1535 good, = np.where(np.isfinite(gainArray))
1537 if len(good) > 1:
1538 # This only works with more than 1 good amp.
1540 # We sort the gains and take the one that is closest
1541 # to the median to use as the reference amplifier for
1542 # gain ratios.
1543 st = np.argsort(gainArray[good])
1544 midAmp = good[st[int(0.5*len(good))]]
1545 midAmpName = ptc.ampNames[midAmp]
1547 log.info("Using amplifier %s as the reference for doLinearityGainRatioFixup.", midAmpName)
1549 # First pass, we need to compute the corrections.
1550 corrections = {}
1551 for ampName in ptc.ampNames:
1552 if not np.isfinite(ptc.gain[ampName]) or ampName == midAmpName:
1553 continue
1555 ratioPtc = ptc.gain[ampName] / ptc.gain[midAmpName]
1557 deltas = ptc.ampOffsets[ampName] - ptc.ampOffsets[midAmpName]
1558 use = (
1559 (ptc.expIdMask[ampName])
1560 & (np.isfinite(deltas))
1561 & (ptc.finalMeans[ampName] >= minAdu)
1562 & (ptc.finalMeans[ampName] <= maxAdu)
1563 & (np.isfinite(ptc.finalMeans[midAmpName]))
1564 & (ptc.expIdMask[midAmpName])
1565 )
1566 if use.sum() < 3:
1567 log.warning("Not enough good amp offset measurements to fix up amp %s "
1568 "gains from amp ratios.", ampName)
1569 continue
1571 ratios = 1. / (deltas / ptc.finalMeans[midAmpName] + 1.0)
1572 ratio = np.median(ratios[use])
1573 corrections[ampName] = ratio / ratioPtc
1575 # For the final correction, we need to make sure that the
1576 # reference amplifier is included. By definition, it has a
1577 # correction factor of 1.0 before any final fix.
1578 corrections[midAmpName] = 1.0
1580 # Adjust the median correction to be 1.0 so we do not
1581 # change the gain of the detector on average.
1582 # This is needed in case the reference amplifier is
1583 # skewed in terms of offsets even though it has the median
1584 # gain.
1585 medCorrection = np.median([corrections[key] for key in corrections])
1587 for ampName in ptc.ampNames:
1588 if ampName not in corrections:
1589 continue
1591 correction = corrections[ampName] / medCorrection
1592 newGain = ptc.gain[ampName] * correction
1593 log.info(
1594 "Adjusting gain from amplifier %s by factor of %.5f (from %.5f to %.5f)",
1595 ampName,
1596 correction,
1597 ptc.gain[ampName],
1598 newGain,
1599 )
1600 # Copying the value should not be necessary, but we record
1601 # it just in case.
1602 ptc.gainUnadjusted[ampName] = ptc.gain[ampName]
1603 ptc.gain[ampName] = newGain
1604 else:
1605 log.warning("Cannot apply ampOffsetGainRatioFixup with fewer than 2 good amplifiers.")
1608class FlatGradientFitter:
1609 """Class to fit various large-scale flat-field gradients.
1611 This fitter will take arrays of x/y/value and (by default) fit a radial
1612 gradient, using a spline function at the specified nodes. The fitter
1613 will also fit out a nuisance parameter for the ratio of the ITL/E2V
1614 throughput (though in general this could work for any such focal plane).
1615 The focal plane origin is set by fp_centroid_x and fp_centroid_y which
1616 is the center of the radial gradient, though it may be modified if
1617 fit_centroid is True. In addition, the fitter may fit a linear gradient
1618 in x/y if fit_gradient is True. The "pivot" of the gradient is at
1619 fp_centroid_x/fp_centroid_y.
1621 Parameters
1622 ----------
1623 nodes : `np.ndarray`
1624 Array of spline nodes for radial spline.
1625 x : `np.ndarray`
1626 Array of x values for points to fit (mm).
1627 y : `np.ndarray`
1628 Array of y values for points to fit (mm).
1629 value : `np.ndarray`
1630 Array of values describing the flat field to fit the gradient.
1631 itl_indices : `np.ndarray`
1632 Array of indices corresponding to the ITL detectors.
1633 constrain_zero : `bool`, optional
1634 Constrain the outermost radial spline value to 0?
1635 fit_centroid : `bool`, optional
1636 Fit an additional centroid offset?
1637 fit_gradient : `bool`, optional
1638 Fit an additional plane gradient?
1639 fp_centroid_x : `float`, optional
1640 Focal plane centroid x (mm).
1641 fp_centroid_y : `float`, optional
1642 Focal plane centroid y (mm).
1643 """
1644 def __init__(
1645 self,
1646 nodes,
1647 x,
1648 y,
1649 value,
1650 itl_indices,
1651 constrain_zero=True,
1652 fit_centroid=False,
1653 fit_gradient=False,
1654 fp_centroid_x=0.0,
1655 fp_centroid_y=0.0
1656 ):
1657 self._nodes = nodes
1658 self._x = x
1659 self._y = y
1660 self._value = value
1661 self._itl_indices = itl_indices
1662 self._fp_centroid_x = fp_centroid_x
1663 self._fp_centroid_y = fp_centroid_y
1665 self._constrain_zero = constrain_zero
1667 self._fit_centroid = fit_centroid
1668 self._fit_gradient = fit_gradient
1670 self.indices = {"spline": np.arange(len(nodes))}
1671 npar = len(nodes)
1673 self._fit_itl_ratio = False
1674 if len(itl_indices) > 0:
1675 self._fit_itl_ratio = True
1676 self.indices["itl_ratio"] = npar
1677 npar += 1
1679 radius = np.sqrt((self._x - self._fp_centroid_x)**2. + (self._y - self._fp_centroid_y)**2.)
1681 if fit_centroid:
1682 self.indices["centroid_delta"] = np.arange(2) + npar
1683 npar += 2
1684 else:
1685 self._radius = radius
1687 if fit_gradient:
1688 self.indices["gradient"] = np.arange(2) + npar
1689 npar += 2
1691 self._npar = npar
1693 self._bounds = [(None, None)]*npar
1695 if self._constrain_zero:
1696 self._bounds[self.indices["spline"][-1]] = (0.0, 0.0)
1698 def compute_p0(self, itl_ratio=None):
1699 """Compute initial guess for fit parameters.
1701 Returns
1702 -------
1703 pars : `np.ndarray`
1704 Array of first guess fit parameters.
1705 """
1706 pars = np.zeros(self._npar)
1707 value = self._value.copy()
1709 if itl_ratio is not None and self._fit_itl_ratio:
1710 value[self._itl_indices] /= itl_ratio
1712 # Initial spline values
1713 radius = np.sqrt((self._x - self._fp_centroid_x)**2. + (self._y - self._fp_centroid_y)**2.)
1714 for i, index in enumerate(self.indices["spline"]):
1715 if i == 0:
1716 low = self._nodes[i]
1717 else:
1718 low = (self._nodes[i - 1] + self._nodes[i])/2.
1719 if i == (len(self._nodes) - 1):
1720 high = self._nodes[i]
1721 else:
1722 high = (self._nodes[i] + self._nodes[i + 1])/2.
1723 u = ((radius > low) & (radius < high))
1724 if u.sum() == 0:
1725 pars[index] = 0.0
1726 else:
1727 pars[index] = np.median(value[u])
1729 if self._constrain_zero:
1730 pars[self.indices["spline"][-1]] = 0.0
1732 spl = Akima1DInterpolator(self._nodes, pars[self.indices["spline"]], method="akima")
1733 model = spl(radius)
1734 resid_ratio = value / model
1736 if self._fit_itl_ratio:
1737 if itl_ratio is not None:
1738 pars[self.indices["itl_ratio"]] = itl_ratio
1739 else:
1740 e2v_indices = np.delete(np.arange(len(self._value)), self._itl_indices)
1742 itl_inner = radius[self._itl_indices] < 0.8*np.max(radius)
1743 e2v_inner = radius[e2v_indices] < 0.8*np.max(radius)
1745 itl_median = np.nanmedian(resid_ratio[self._itl_indices][itl_inner])
1746 e2v_median = np.nanmedian(resid_ratio[e2v_indices][e2v_inner])
1748 pars[self.indices["itl_ratio"]] = itl_median / e2v_median
1750 if self._fit_centroid:
1751 pars[self.indices["centroid_delta"]] = [0.0, 0.0]
1753 if self._fit_gradient:
1754 resid_ratio = self._value / model
1756 ok = (np.isfinite(resid_ratio) & (model > 0.5))
1758 fit = np.polyfit(self._x[ok], resid_ratio[ok], 1)
1759 pars[self.indices["gradient"][0]] = fit[0]
1761 fit = np.polyfit(self._y[ok], resid_ratio[ok], 1)
1762 pars[self.indices["gradient"][1]] = fit[0]
1764 return pars
1766 def fit(self, p0, freeze_itl_ratio=False, fit_eps=1e-8, fit_gtol=1e-10):
1767 """Do a non-linear minimization to fit the parameters.
1769 Parameters
1770 ----------
1771 p0 : `np.ndarray`
1772 Array of initial parameter estimates.
1773 freeze_itl_ratio : `bool`, optional
1774 Freeze the ITL ratio in the fit?
1775 fit_eps : `float`, optional
1776 Value of ``eps`` to send to the scipy minimizer.
1777 fit_gtol : `float`, optional
1778 Value of ``gtol`` to send to the scipy minimizer.
1780 Returns
1781 -------
1782 pars : `np.ndarray`
1783 Array of parameters. Use ``fitter.indices`` for the
1784 dictionary to map parameters to subsets.
1785 """
1786 bounds = self._bounds
1787 if freeze_itl_ratio and self._fit_itl_ratio:
1788 ind = self.indices["itl_ratio"]
1789 par = p0[ind]
1790 bounds[ind] = (par, par)
1792 res = minimize(
1793 self,
1794 p0,
1795 method="L-BFGS-B",
1796 jac=False,
1797 bounds=bounds,
1798 options={
1799 "maxfun": 10000,
1800 "maxiter": 10000,
1801 "maxcor": 20,
1802 "eps": fit_eps,
1803 "gtol": fit_gtol,
1804 },
1805 callback=None,
1806 )
1807 pars = res.x
1809 return pars
1811 def compute_model(self, pars):
1812 """Compute the model given a set of parameters.
1814 Parameters
1815 ----------
1816 pars : `np.ndarray`
1817 Parameter array to compute model.
1819 Returns
1820 -------
1821 model_array : `np.ndarray`
1822 Array of model parameters at the input x/y.
1823 """
1824 spl = Akima1DInterpolator(self._nodes, pars[self.indices["spline"]], method="akima")
1825 if self._fit_centroid:
1826 centroid_delta = pars[self.indices["centroid_delta"]]
1827 centroid_x = self._fp_centroid_x + centroid_delta[0]
1828 centroid_y = self._fp_centroid_y + centroid_delta[1]
1829 radius = np.sqrt((self._x - centroid_x)**2. + (self._y - centroid_y)**2.)
1830 else:
1831 radius = self._radius
1833 model = spl(np.clip(radius, self._nodes[0], self._nodes[-1]))
1834 if self._fit_itl_ratio:
1835 model[self._itl_indices] *= pars[self.indices["itl_ratio"]]
1837 if self._fit_gradient:
1838 a, b = pars[self.indices["gradient"]]
1839 gradient = 1 + a*(self._x - self._fp_centroid_x) + b*(self._y - self._fp_centroid_y)
1840 model /= gradient
1842 return model
1844 def __call__(self, pars):
1845 """Compute the cost function for a set of parameters.
1847 Parameters
1848 ----------
1849 pars : `np.ndarray`
1850 Parameter array to compute model.
1852 Returns
1853 -------
1854 cost : `float`
1855 Cost value computed from the absolute deviation.
1856 """
1858 model = self.compute_model(pars)
1860 absdev = np.abs(self._value - model)
1861 t = np.sum(absdev.astype(np.float64))
1863 return t
1866class FlatGainRatioFitter:
1867 """Class to fit amplifier gain ratios, removing a simple gradient.
1869 This fitter will take arrays of x/y/amp_num/value and fit amplifier
1870 gain ratios, using one amplifier as the fixed point. The fitter
1871 uses a low-order chebyshev polynomial to fit out the gradient, with
1872 amp ratios on top of this.
1874 Parameters
1875 ----------
1876 bbox : `lsst.geom.Box2I`
1877 Bounding box for Chebyshev polynomial gradient.
1878 order : `int`
1879 Order of Chebyshev polynomial.
1880 x : `np.ndarray`
1881 Array of x values for points to fit (detector pixels).
1882 y : `np.ndarray`
1883 Array of y values for points to fit (detector pixels).
1884 amp_index : `np.ndarray`
1885 Array of amp numbers associated with each x/y pair.
1886 value : `np.ndarray`
1887 Flat value at each x/y pair.
1888 amps : `np.ndarray`
1889 Array of unique amplifier numbers that will be parameterized.
1890 Any of these amps that does not have any data with the same
1891 amp_index will be set to 1.0.
1892 fixed_amp_index : `int`
1893 Amplifier number to keep fixed.
1894 """
1895 def __init__(self, bbox, order, x, y, amp_index, value, amps, fixed_amp_index):
1896 self._bbox = bbox
1897 self._order = order
1898 self._x = x.astype(np.float64)
1899 self._y = y.astype(np.float64)
1900 self._amp_index = amp_index
1901 self._value = value.astype(np.float64)
1902 self._fixed_amp_index = fixed_amp_index
1904 self.indices = {"chebyshev": np.arange((order + 1) * (order + 1))}
1905 npar = len(self.indices["chebyshev"])
1907 self._amps = amps
1908 self._n_amp = len(self._amps)
1910 self.indices["amp_pars"] = np.arange(self._n_amp) + npar
1911 npar += self._n_amp
1913 self._npar = npar
1915 self._amp_indices = {}
1916 for i in range(self._n_amp):
1917 amp_index = self._amps[i]
1918 self._amp_indices[amp_index] = (self._amp_index == amp_index)
1920 def fit(self, n_iter=10):
1921 """Fit the amp ratio parameters.
1923 This uses an iterative fit, where it fits a Chebyshev gradient,
1924 computes amp ratios, and re-fits the gradient.
1926 Returns
1927 -------
1928 pars : `np.ndarray`
1929 Chebyshev parameters and amp offset parameters.
1930 """
1931 value = self._value.copy()
1933 control = lsst.afw.math.ChebyshevBoundedFieldControl()
1934 control.orderX = self._order
1935 control.orderY = self._order
1936 control.triangular = False
1938 pars = np.zeros(self._npar)
1940 pars[self.indices["amp_pars"]] = 1.0
1942 for i in range(n_iter):
1943 field = lsst.afw.math.ChebyshevBoundedField.fit(
1944 self._bbox,
1945 self._x,
1946 self._y,
1947 value,
1948 control,
1949 )
1951 pars[self.indices["chebyshev"][:]] = field.getCoefficients().ravel()
1953 ratio = self._value / field.evaluate(self._x, self._y)
1955 fixed_med = np.median(ratio[self._amp_indices[self._fixed_amp_index]])
1956 ratio /= fixed_med
1958 pars[self.indices["amp_pars"][self._fixed_amp_index]] = 1.0
1960 value = self._value.copy()
1962 for j in range(self._n_amp):
1963 amp_index = self._amps[j]
1965 if np.sum(self._amp_indices[amp_index]) == 0:
1966 continue
1968 pars[self.indices["amp_pars"][j]] = np.median(ratio[self._amp_indices[amp_index]])
1969 value[self._amp_indices[amp_index]] *= pars[self.indices["amp_pars"][j]]
1971 return pars
1973 def compute_model(self, pars):
1974 """Compute the gradient/amp ratio model for a given set of parameters.
1976 Parameters
1977 ----------
1978 pars : `np.ndarray`
1979 Chebyshev parameters and amp offset parameters.
1981 Returns
1982 -------
1983 model : `np.ndarray`
1984 The model at each x/y pair.
1985 """
1986 field = lsst.afw.math.ChebyshevBoundedField(
1987 self._bbox,
1988 pars[self.indices["chebyshev"]].reshape(self._order + 1, self._order + 1),
1989 )
1990 model = field.evaluate(self._x, self._y)
1992 for i in range(self._n_amp):
1993 amp_index = self._amps[i]
1995 model[self._amp_indices[amp_index]] *= pars[self.indices["amp_pars"][i]]
1997 return model
1999 def __call__(self, pars):
2000 """Compute the cost function for a set of parameters.
2002 Parameters
2003 ----------
2004 pars : `np.ndarray`
2005 Chebyshev parameters and amp offset parameters.
2007 Returns
2008 -------
2009 cost : `float`
2010 Cost value computed from the absolute deviation.
2011 """
2012 model = self.compute_model(pars)
2014 absdev = np.abs(self._value - model)
2015 t = np.sum(absdev.astype(np.float64))
2017 return t
2020def bin_focal_plane(
2021 exposure_handle_dict,
2022 detector_boundary,
2023 bin_factor,
2024 defect_handle_dict={},
2025 include_itl_flag=True,
2026):
2027 """Bin all the detectors into the full focal plane.
2029 This function reads in images; takes a simple average if there
2030 are more than one input per detector; excludes detector edges;
2031 and bins according to the bin factor. The output is a struct
2032 with focal plane coordinates, detector numbers, and a flag
2033 if the detector is an ITL detector.
2035 Parameters
2036 ----------
2037 exposure_handle_dict : `dict`
2038 Dict keyed by detector (`int`), each element is a list
2039 of `lsst.daf.butler.DeferredDatasetHandle` that will be averaged.
2040 detector_boundary : `int`
2041 Boundary of the detector to ignore (pixels).
2042 bin_factor : `int`
2043 Binning factor. Detectors will be cropped (after applying the
2044 ``detector_boundary``) such that there are no partially
2045 covered binned pixels.
2046 defect_handle_dict : `dict`, optional
2047 Dict keyed by detector (`int`), each element is a
2048 `lsst.daf.butler.DeferredDatasetHandle` for defects to be applied.
2049 include_itl_flag : `bool`, optional
2050 Include a flag for which detectors are ITL?
2052 Returns
2053 -------
2054 binned : `astropy.table.Table`
2055 Table with focal plane positions at the center of each bin
2056 (``xf``, ``yf``); average image values (``value``); and detector
2057 number (``detector``).
2058 """
2059 xf_arrays = []
2060 yf_arrays = []
2061 value_arrays = []
2062 detector_arrays = []
2063 itl_arrays = []
2065 for det in exposure_handle_dict.keys():
2066 flat = exposure_handle_dict[det].get()
2067 defect_handle = defect_handle_dict.get(det, None)
2068 if defect_handle is not None:
2069 defects = defect_handle.get()
2070 else:
2071 defects = None
2073 detector = flat.getDetector()
2075 # Mask out defects if we have them.
2076 if defects is not None:
2077 for defect in defects:
2078 flat.image[defect.getBBox()].array[:, :] = np.nan
2080 # Mask NO_DATA pixels if we have them.
2081 no_data = ((flat.mask.array[:, :] & flat.mask.getPlaneBitMask("NO_DATA")) > 0)
2082 flat.image.array[no_data] = np.nan
2084 # Bin the image, avoiding the boundary and the masked pixels.
2085 # We also make sure we are using an integral number of
2086 # steps to avoid partially covered binned pixels.
2088 arr = flat.image.array
2090 n_step_y = (arr.shape[0] - (2 * detector_boundary)) // bin_factor
2091 y_min = detector_boundary
2092 y_max = bin_factor * n_step_y + y_min
2093 n_step_x = (arr.shape[1] - (2 * detector_boundary)) // bin_factor
2094 x_min = detector_boundary
2095 x_max = bin_factor * n_step_x + x_min
2097 arr = arr[y_min: y_max, x_min: x_max]
2098 binned = arr.reshape((n_step_y, bin_factor, n_step_x, bin_factor))
2099 with warnings.catch_warnings():
2100 warnings.filterwarnings("ignore", r"Mean of empty")
2101 binned = np.nanmean(binned, axis=1)
2102 binned = np.nanmean(binned, axis=2)
2104 xx = np.arange(binned.shape[1]) * bin_factor + bin_factor / 2. + x_min
2105 yy = np.arange(binned.shape[0]) * bin_factor + bin_factor / 2. + y_min
2106 x, y = np.meshgrid(xx, yy)
2107 x = x.ravel()
2108 y = y.ravel()
2109 value = binned.ravel()
2111 # Transform to focal plane coordinates.
2112 transform = detector.getTransform(lsst.afw.cameraGeom.PIXELS, lsst.afw.cameraGeom.FOCAL_PLANE)
2113 xy = np.vstack((x, y))
2114 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)
2115 xf = xf.ravel()
2116 yf = yf.ravel()
2118 if include_itl_flag:
2119 is_itl = np.zeros(len(value), dtype=np.bool_)
2120 # We use this check so that ITL matches ITL science detectors,
2121 # ITL_WF wavefront detectors, and pseudoITL test detectors.
2122 is_itl[:] = ("ITL" in detector.getPhysicalType())
2124 xf_arrays.append(xf)
2125 yf_arrays.append(yf)
2126 value_arrays.append(value)
2127 detector_arrays.append(np.full_like(xf, det, dtype=np.int32))
2128 if include_itl_flag:
2129 itl_arrays.append(is_itl)
2131 xf = np.concatenate(xf_arrays)
2132 yf = np.concatenate(yf_arrays)
2133 value = np.concatenate(value_arrays)
2134 detector = np.concatenate(detector_arrays)
2136 binned = Table(
2137 np.zeros(
2138 len(xf),
2139 dtype=[
2140 ("xf", "f8"),
2141 ("yf", "f8"),
2142 ("value", "f8"),
2143 ("detector", "i4"),
2144 ],
2145 )
2146 )
2147 binned["xf"] = xf
2148 binned["yf"] = yf
2149 binned["value"] = value
2150 binned["detector"] = detector
2152 if include_itl_flag:
2153 binned["itl"] = np.concatenate(itl_arrays).astype(np.bool_)
2155 return binned
2158def bin_flat(ptc, exposure, bin_factor=8, amp_boundary=20, apply_gains=True, gain_ratios=None):
2159 """Bin a flat image, being careful with amplifier edges.
2161 This will optionally apply gains, and apply any gain
2162 ratios.
2164 Parameters
2165 ----------
2166 ptc : `lsst.ip.isr.PhotonTransferCurveDatasets`
2167 PTC dataset with relevant gains.
2168 exposure : `lsst.afw.image.Exposure`
2169 Exposure to bin.
2170 bin_factor : `int`, optional
2171 Binning factor.
2172 amp_boundary : `int`, optional
2173 Boundary around each amp to ignore in binning.
2174 apply_gains : `bool`, optional
2175 Apply gains before binning?
2176 gain_ratios : `np.ndarray`, optional
2177 Array of gain ratios to apply.
2179 Returns
2180 -------
2181 binned : `astropy.table.Table`
2182 Table with detector coordinates at the center of each bin
2183 (``xd``, ``yd``); average image values (``value``); and amplifier
2184 index (``amp_index``).
2185 """
2186 detector = exposure.getDetector()
2188 for i, amp_name in enumerate(ptc.ampNames):
2189 bbox = detector[amp_name].getBBox()
2190 if amp_name in ptc.badAmps:
2191 exposure[bbox].image.array[:, :] = np.nan
2192 continue
2194 if apply_gains:
2195 exposure[bbox].image.array[:, :] *= ptc.gainUnadjusted[amp_name]
2196 if gain_ratios is not None:
2197 exposure[bbox].image.array[:, :] /= gain_ratios[i]
2199 # Next we bin the detector, avoiding amp edges.
2200 xd_arrays = []
2201 yd_arrays = []
2202 value_arrays = []
2203 amp_arrays = []
2205 for i, amp in enumerate(detector):
2206 arr = exposure[amp.getBBox()].image.array
2208 n_step_y = (arr.shape[0] - (2 * amp_boundary)) // bin_factor
2209 y_min = amp_boundary
2210 y_max = bin_factor * n_step_y + y_min
2211 n_step_x = (arr.shape[1] - (2 * amp_boundary)) // bin_factor
2212 x_min = amp_boundary
2213 x_max = bin_factor * n_step_x + x_min
2215 arr = arr[y_min: y_max, x_min: x_max]
2216 binned = arr.reshape((n_step_y, bin_factor, n_step_x, bin_factor))
2217 with warnings.catch_warnings():
2218 warnings.filterwarnings("ignore", r"Mean of empty")
2219 binned = np.nanmean(binned, axis=1)
2220 binned = np.nanmean(binned, axis=2)
2222 xx = np.arange(binned.shape[1]) * bin_factor + bin_factor / 2. + x_min
2223 yy = np.arange(binned.shape[0]) * bin_factor + bin_factor / 2. + y_min
2224 x, y = np.meshgrid(xx, yy)
2225 x = x.ravel() + amp.getBBox().getBeginX()
2226 y = y.ravel() + amp.getBBox().getBeginY()
2227 value = binned.ravel()
2229 xd_arrays.append(x)
2230 yd_arrays.append(y)
2231 value_arrays.append(value)
2232 amp_arrays.append(np.full(len(x), i))
2234 xd = np.concatenate(xd_arrays)
2235 yd = np.concatenate(yd_arrays)
2236 value = np.concatenate(value_arrays)
2237 amp_index = np.concatenate(amp_arrays)
2239 binned = np.zeros(
2240 len(xd),
2241 dtype=[
2242 ("xd", "f8"),
2243 ("yd", "f8"),
2244 ("value", "f8"),
2245 ("amp_index", "i4"),
2246 ],
2247 )
2248 binned["xd"] = xd
2249 binned["yd"] = yd
2250 binned["value"] = value
2251 binned["amp_index"] = amp_index
2253 return Table(data=binned)