Coverage for python/lsst/cp/pipe/utils.py: 7%
773 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-24 02:10 -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_frac_correction : `float`, optional
1059 Maximum fractional correction.
1060 max_correction : `float`, optional
1061 Maximum correction (unscaled).
1062 """
1063 def __init__(
1064 self,
1065 nodes,
1066 grouping_values,
1067 pd,
1068 mu,
1069 mask=None,
1070 log=None,
1071 fit_offset=True,
1072 fit_weights=False,
1073 weight_pars_start=[1.0, 0.0],
1074 fit_temperature=False,
1075 temperature_scaled=None,
1076 max_signal_nearly_linear=None,
1077 fit_temporal=False,
1078 mjd_scaled=None,
1079 max_frac_correction=0.25,
1080 max_correction=np.inf,
1081 ):
1082 self._pd = pd
1083 self._mu = mu
1084 self._grouping_values = grouping_values
1085 self.log = log if log else logging.getLogger(__name__)
1086 self._fit_offset = fit_offset
1087 self._fit_weights = fit_weights
1088 self._weight_pars_start = weight_pars_start
1089 self._fit_temperature = fit_temperature
1090 self._fit_temporal = fit_temporal
1091 self._max_frac_correction = max_frac_correction
1092 self._max_correction = max_correction
1094 self._nodes = nodes
1095 if nodes[0] != 0.0:
1096 raise ValueError("First node must be 0.0")
1097 if not np.all(np.diff(nodes) > 0):
1098 raise ValueError("Nodes must be sorted with no repeats.")
1100 # Find the group indices.
1101 u_group_values = np.unique(self._grouping_values)
1102 self.ngroup = len(u_group_values)
1104 self.group_indices = []
1105 for i in range(self.ngroup):
1106 self.group_indices.append(np.where(self._grouping_values == u_group_values[i])[0])
1108 # Weight values. Outliers will be set to 0.
1109 if mask is None:
1110 _mask = np.ones(len(mu), dtype=np.bool_)
1111 else:
1112 _mask = mask
1113 self._w = self.compute_weights(self._weight_pars_start, self._mu, _mask)
1115 if temperature_scaled is None:
1116 temperature_scaled = np.zeros(len(self._mu))
1117 else:
1118 if len(np.atleast_1d(temperature_scaled)) != len(self._mu):
1119 raise ValueError("temperature_scaled must be the same length as input mu.")
1120 self._temperature_scaled = temperature_scaled
1122 if mjd_scaled is None:
1123 mjd_scaled = np.zeros(len(self._mu))
1124 else:
1125 if len(np.atleast_1d(mjd_scaled)) != len(self._mu):
1126 raise ValueError("mjd_scaled must be the same length as input mu.")
1127 self._mjd_scaled = mjd_scaled
1129 # Values to regularize spline fit.
1130 if max_signal_nearly_linear is None:
1131 max_signal_nearly_linear = self._mu[self.mask].max()
1132 self._max_signal_nearly_linear = max_signal_nearly_linear
1133 self._x_regularize = np.linspace(0.0, self._max_signal_nearly_linear, 100)
1135 # Set up the indices for the fit parameters.
1136 self.par_indices = {
1137 "values": np.arange(len(self._nodes)),
1138 "groups": len(self._nodes) + np.arange(self.ngroup),
1139 "offset": np.zeros(0, dtype=np.int64),
1140 "weight_pars": np.zeros(0, dtype=np.int64),
1141 "temperature_coeff": np.zeros(0, dtype=np.int64),
1142 "temporal_coeff": np.zeros(0, dtype=np.int64),
1143 }
1144 if self._fit_offset:
1145 self.par_indices["offset"] = np.arange(1) + (
1146 len(self.par_indices["values"])
1147 + len(self.par_indices["groups"])
1148 )
1149 if self._fit_weights:
1150 self.par_indices["weight_pars"] = np.arange(2) + (
1151 len(self.par_indices["values"])
1152 + len(self.par_indices["groups"])
1153 + len(self.par_indices["offset"])
1154 )
1155 if self._fit_temperature:
1156 self.par_indices["temperature_coeff"] = np.arange(1) + (
1157 len(self.par_indices["values"])
1158 + len(self.par_indices["groups"])
1159 + len(self.par_indices["offset"])
1160 + len(self.par_indices["weight_pars"])
1161 )
1162 if self._fit_temporal:
1163 self.par_indices["temporal_coeff"] = np.arange(1) + (
1164 len(self.par_indices["values"])
1165 + len(self.par_indices["groups"])
1166 + len(self.par_indices["offset"])
1167 + len(self.par_indices["weight_pars"])
1168 + len(self.par_indices["temperature_coeff"])
1169 )
1171 @staticmethod
1172 def compute_weights(weight_pars, mu, mask):
1173 w = 1./np.sqrt(weight_pars[0]**2. + weight_pars[1]**2./mu)
1174 w[~mask] = 0.0
1176 return w
1178 def estimate_p0(self):
1179 """Estimate initial fit parameters.
1181 Returns
1182 -------
1183 p0 : `np.ndarray`
1184 Parameter array, with spline values (one for each node) followed
1185 by proportionality constants (one for each group); one extra
1186 for the offset O (if fit_offset was set to True); two extra
1187 for the weights (if fit_weights was set to True); one
1188 extra for the temperature coefficient (if fit_temperature was
1189 set to True); and one extra for the temporal coefficient (if
1190 fit_temporal was set to True).
1191 """
1192 npt = (len(self.par_indices["values"])
1193 + len(self.par_indices["groups"])
1194 + len(self.par_indices["offset"])
1195 + len(self.par_indices["weight_pars"])
1196 + len(self.par_indices["temperature_coeff"])
1197 + len(self.par_indices["temporal_coeff"]))
1198 p0 = np.zeros(npt)
1200 # Do a simple linear fit for each group.
1201 for i, indices in enumerate(self.group_indices):
1202 mask = self.mask[indices]
1203 mu = self._mu[indices][mask]
1204 pd = self._pd[indices][mask]
1205 to_fit = (mu < self._max_signal_nearly_linear)
1206 # If we have no points (low PTC turnoff) then
1207 # we just use the first half for our initial value.
1208 if to_fit.sum() == 0:
1209 to_fit = (mu < np.median(mu))
1210 linfit = np.polyfit(pd[to_fit], mu[to_fit], 1)
1211 p0[self.par_indices["groups"][i]] = linfit[0]
1213 # Look at the residuals...
1214 ratio_model = self.compute_ratio_model(
1215 self._nodes,
1216 self.group_indices,
1217 self.par_indices,
1218 p0,
1219 self._pd,
1220 self._mu,
1221 self._temperature_scaled,
1222 self._mjd_scaled,
1223 )
1224 # ...and adjust the linear parameters accordingly.
1225 p0[self.par_indices["groups"]] *= np.median(ratio_model[self.mask])
1227 # Re-compute the residuals.
1228 ratio_model2 = self.compute_ratio_model(
1229 self._nodes,
1230 self.group_indices,
1231 self.par_indices,
1232 p0,
1233 self._pd,
1234 self._mu,
1235 self._temperature_scaled,
1236 self._mjd_scaled,
1237 )
1239 # And compute a first guess of the spline nodes.
1240 bins = np.clip(np.searchsorted(self._nodes, self._mu[self.mask]), 0, len(self._nodes) - 1)
1241 tot_arr = np.zeros(len(self._nodes))
1242 n_arr = np.zeros(len(self._nodes), dtype=int)
1243 np.add.at(tot_arr, bins, ratio_model2[self.mask])
1244 np.add.at(n_arr, bins, 1)
1246 ratio = np.ones(len(self._nodes))
1247 ratio[n_arr > 0] = tot_arr[n_arr > 0]/n_arr[n_arr > 0]
1248 ratio[0] = 1.0
1249 p0[self.par_indices["values"]] = (ratio - 1) * self._nodes
1251 if self._fit_offset:
1252 p0[self.par_indices["offset"]] = 0.0
1254 if self._fit_weights:
1255 p0[self.par_indices["weight_pars"]] = self._weight_pars_start
1257 return p0
1259 @staticmethod
1260 def compute_ratio_model(
1261 nodes,
1262 group_indices,
1263 par_indices,
1264 pars,
1265 pd,
1266 mu,
1267 temperature_scaled,
1268 mjd_scaled,
1269 return_spline=False,
1270 ):
1271 """Compute the ratio model values.
1273 Parameters
1274 ----------
1275 nodes : `np.ndarray` (M,)
1276 Array of node positions.
1277 group_indices : `list` [`np.ndarray`]
1278 List of group indices, one array for each group.
1279 par_indices : `dict`
1280 Dictionary showing which indices in the pars belong to
1281 each set of fit values.
1282 pars : `np.ndarray`
1283 Parameter array, with spline values (one for each node) followed
1284 by proportionality constants (one for each group); one extra
1285 for the offset O (if fit_offset was set to True); two extra
1286 for the weights (if fit_weights was set to True); and one
1287 extra for the temperature coefficient (if fit_temperature was
1288 set to True).
1289 pd : `np.ndarray` (N,)
1290 Array of photodiode measurements.
1291 mu : `np.ndarray` (N,)
1292 Array of flat means.
1293 temperature_scaled : `np.ndarray` (N,)
1294 Array of scaled temperature values.
1295 mjd_scaled : `np.ndarray` (N,)
1296 Array of scaled mjd values.
1297 return_spline : `bool`, optional
1298 Return the spline interpolation as well as the model ratios?
1300 Returns
1301 -------
1302 ratio_models : `np.ndarray` (N,)
1303 Model ratio, (mu_i - S(mu_i) - O)/(k_j * D_i)
1304 spl : `scipy.interpolate.Akima1DInterpolator`
1305 Spline interpolator (returned if return_spline=True).
1306 """
1307 spl = Akima1DInterpolator(nodes, pars[par_indices["values"]], method="akima")
1309 # Check if we want to do just the left or both with temp scale.
1310 if len(par_indices["temperature_coeff"]) == 1:
1311 mu_corr = mu*(1. + pars[par_indices["temperature_coeff"]]*temperature_scaled)
1312 else:
1313 mu_corr = mu
1315 if len(par_indices["temporal_coeff"]) == 1:
1316 mu_corr = mu_corr*(1. + pars[par_indices["temporal_coeff"]]*mjd_scaled)
1318 numerator = mu_corr - spl(np.clip(mu_corr, nodes[0], nodes[-1]))
1319 if len(par_indices["offset"]) == 1:
1320 numerator -= pars[par_indices["offset"][0]]
1321 denominator = pd.copy()
1322 ngroup = len(group_indices)
1323 kj = pars[par_indices["groups"]]
1324 for j in range(ngroup):
1325 denominator[group_indices[j]] *= kj[j]
1327 if return_spline:
1328 return numerator / denominator, spl
1329 else:
1330 return numerator / denominator
1332 def fit(self, p0, min_iter=3, max_iter=20, max_rejection_per_iteration=5, n_sigma_clip=5.0):
1333 """
1334 Perform iterative fit for linear + spline model with offsets.
1336 Parameters
1337 ----------
1338 p0 : `np.ndarray`
1339 Initial parameter array, with spline values (one for each node)
1340 followed by proportionality constants (one for each group); one
1341 extra for the offset O (if fit_offset was set to True); two extra
1342 for the weights (if fit_weights was set to True); one
1343 extra for the temperature coefficient (if fit_temperature was
1344 set to True); and one extra for the temporal coefficient (if
1345 fit_temporal was set to True).
1346 min_iter : `int`, optional
1347 Minimum number of fit iterations.
1348 max_iter : `int`, optional
1349 Maximum number of fit iterations.
1350 max_rejection_per_iteration : `int`, optional
1351 Maximum number of points to reject per iteration.
1352 n_sigma_clip : `float`, optional
1353 Number of sigma to do clipping in each iteration.
1354 """
1355 init_params = p0
1356 for k in range(max_iter):
1357 params, cov_params, _, msg, ierr = leastsq(
1358 self,
1359 init_params,
1360 full_output=True,
1361 ftol=1e-5,
1362 maxfev=12000,
1363 )
1364 init_params = params.copy()
1366 # We need to cut off the constraints at the end (there are more
1367 # residuals than data points.)
1368 res = self(params)[: len(self._w)]
1369 std_res = median_abs_deviation(res[self.good_points], scale="normal")
1370 sample = len(self.good_points)
1372 # We don't want to reject too many outliers at once.
1373 if sample > max_rejection_per_iteration:
1374 sres = np.sort(np.abs(res))
1375 cut = max(sres[-max_rejection_per_iteration], std_res*n_sigma_clip)
1376 else:
1377 cut = std_res*n_sigma_clip
1379 outliers = np.abs(res) > cut
1380 self._w[outliers] = 0
1381 if outliers.sum() != 0:
1382 self.log.info(
1383 "After iteration %d there are %d outliers (of %d).",
1384 k,
1385 outliers.sum(),
1386 sample,
1387 )
1388 elif k >= min_iter:
1389 self.log.info("After iteration %d there are no more outliers.", k)
1390 break
1392 return params
1394 @property
1395 def mask(self):
1396 return (self._w > 0)
1398 @property
1399 def good_points(self):
1400 return self.mask.nonzero()[0]
1402 def compute_chisq_dof(self, pars):
1403 """Compute the chi-squared per degree of freedom for a set of pars.
1405 Parameters
1406 ----------
1407 pars : `np.ndarray`
1408 Parameter array.
1410 Returns
1411 -------
1412 chisq_dof : `float`
1413 Chi-squared per degree of freedom.
1414 """
1415 resids = self(pars)[0: len(self.mask)]
1416 chisq = np.sum(resids[self.mask]**2.)
1417 dof = self.mask.sum() - self.ngroup
1418 if self._fit_temporal:
1419 dof -= 1
1420 if self._fit_temperature:
1421 dof -= 1
1422 if self._fit_offset:
1423 dof -= 1
1424 if self._fit_weights:
1425 dof -= 2
1427 return chisq/dof
1429 def __call__(self, pars):
1431 ratio_model, spl = self.compute_ratio_model(
1432 self._nodes,
1433 self.group_indices,
1434 self.par_indices,
1435 pars,
1436 self._pd,
1437 self._mu,
1438 self._temperature_scaled,
1439 self._mjd_scaled,
1440 return_spline=True,
1441 )
1443 _mask = self.mask
1444 # Update the weights if we are fitting them.
1445 if self._fit_weights:
1446 self._w = self.compute_weights(pars[self.par_indices["weight_pars"]], self._mu, _mask)
1447 resid = self._w*(ratio_model - 1.0)
1449 # Ensure masked points have 0 residual.
1450 resid[~_mask] = 0.0
1452 constraint = [1e3 * np.mean(spl(np.clip(self._x_regularize, self._nodes[0], self._nodes[-1])))]
1453 # 0 should transform to 0
1454 constraint.append(spl(0)*1e10)
1455 # Use a Jeffreys prior on the weight if we are fitting it.
1456 if self._fit_weights:
1457 # This factor ensures that log(fact * w) is negative.
1458 fact = 1e-3 / self._w.max()
1459 # We only add non-zero weights to the constraint array.
1460 log_w = np.sqrt(-2.*np.log(fact*self._w[self._w > 0]))
1461 constraint = np.hstack([constraint, log_w])
1463 # Don't let it get to >5% correction.
1464 values = pars[self.par_indices["values"]]
1465 if np.abs(values[-1])/self._nodes[-1] > self._max_frac_correction \
1466 or np.abs(values[-1]) > self._max_correction:
1467 extra_constraint = 1e10
1468 else:
1469 extra_constraint = 0
1471 return np.hstack([resid, constraint, extra_constraint])
1474def getReadNoise(exposure, ampName, taskMetadata=None, log=None):
1475 """Gets readout noise for an amp from ISR metadata.
1477 If possible, this attempts to get the now-standard headers
1478 added to the exposure itself. If not found there, the ISR
1479 TaskMetadata is searched. If neither of these has the value,
1480 warn and set the read noise to NaN.
1482 Parameters
1483 ----------
1484 exposure : `lsst.afw.image.Exposure`
1485 Exposure to check for read noise first.
1486 ampName : `str`
1487 Amplifier name.
1488 taskMetadata : `lsst.pipe.base.TaskMetadata`, optional
1489 List of exposures metadata from ISR for this exposure.
1490 log : `logging.logger`, optional
1491 Log for messages.
1493 Returns
1494 -------
1495 readNoise : `float`
1496 The read noise for this set of exposure/amplifier.
1497 """
1498 exposureMetadata = exposure.getMetadata()
1500 # Try from the exposure first.
1501 expectedKey = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"
1502 if expectedKey in exposureMetadata:
1503 return exposureMetadata[expectedKey]
1505 # If not, try getting it from the task metadata.
1506 if taskMetadata:
1507 expectedKey = f"RESIDUAL STDEV {ampName}"
1508 if "isr" in taskMetadata:
1509 if expectedKey in taskMetadata["isr"]:
1510 return taskMetadata["isr"][expectedKey]
1512 log = log if log else logging.getLogger(__name__)
1513 log.warning("Median readout noise from ISR metadata for amp %s "
1514 "could not be found.", ampName)
1515 return np.nan
1518def ampOffsetGainRatioFixup(ptc, minAdu, maxAdu, log=None):
1519 """Apply gain ratio fixup using amp offsets.
1521 Parameters
1522 ----------
1523 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
1524 Input PTC. Will be modified in place.
1525 minAdu : `float`
1526 Minimum number of ADU in mean of amplifier to use for
1527 computing gain ratio fixup.
1528 maxAdu : `float`
1529 Maximum number of ADU in mean of amplifier to use for
1530 computing gain ratio fixup.
1531 log : `lsst.utils.logging.LsstLogAdapter`, optional
1532 Log object.
1533 """
1534 # We need to find the reference amplifier.
1535 # Find an amp near the middle to use as a pivot.
1536 if log is None:
1537 log = logging.getLogger(__name__)
1539 # We first check for which amps have gain measurements
1540 # (fully bad amps are filled with NaN for gain.)
1541 gainArray = np.zeros(len(ptc.ampNames))
1542 for i, ampName in enumerate(ptc.ampNames):
1543 gainArray[i] = ptc.gain[ampName]
1544 good, = np.where(np.isfinite(gainArray))
1546 if len(good) > 1:
1547 # This only works with more than 1 good amp.
1549 # We sort the gains and take the one that is closest
1550 # to the median to use as the reference amplifier for
1551 # gain ratios.
1552 st = np.argsort(gainArray[good])
1553 midAmp = good[st[int(0.5*len(good))]]
1554 midAmpName = ptc.ampNames[midAmp]
1556 log.info("Using amplifier %s as the reference for doLinearityGainRatioFixup.", midAmpName)
1558 # First pass, we need to compute the corrections.
1559 corrections = {}
1560 for ampName in ptc.ampNames:
1561 if not np.isfinite(ptc.gain[ampName]) or ampName == midAmpName:
1562 continue
1564 ratioPtc = ptc.gain[ampName] / ptc.gain[midAmpName]
1566 deltas = ptc.ampOffsets[ampName] - ptc.ampOffsets[midAmpName]
1567 use = (
1568 (ptc.expIdMask[ampName])
1569 & (np.isfinite(deltas))
1570 & (ptc.finalMeans[ampName] >= minAdu)
1571 & (ptc.finalMeans[ampName] <= maxAdu)
1572 & (np.isfinite(ptc.finalMeans[midAmpName]))
1573 & (ptc.expIdMask[midAmpName])
1574 )
1575 if use.sum() < 3:
1576 log.warning("Not enough good amp offset measurements to fix up amp %s "
1577 "gains from amp ratios.", ampName)
1578 continue
1580 ratios = 1. / (deltas / ptc.finalMeans[midAmpName] + 1.0)
1581 ratio = np.median(ratios[use])
1582 corrections[ampName] = ratio / ratioPtc
1584 # For the final correction, we need to make sure that the
1585 # reference amplifier is included. By definition, it has a
1586 # correction factor of 1.0 before any final fix.
1587 corrections[midAmpName] = 1.0
1589 # Adjust the median correction to be 1.0 so we do not
1590 # change the gain of the detector on average.
1591 # This is needed in case the reference amplifier is
1592 # skewed in terms of offsets even though it has the median
1593 # gain.
1594 medCorrection = np.median([corrections[key] for key in corrections])
1596 for ampName in ptc.ampNames:
1597 if ampName not in corrections:
1598 continue
1600 correction = corrections[ampName] / medCorrection
1601 newGain = ptc.gain[ampName] * correction
1602 log.info(
1603 "Adjusting gain from amplifier %s by factor of %.5f (from %.5f to %.5f)",
1604 ampName,
1605 correction,
1606 ptc.gain[ampName],
1607 newGain,
1608 )
1609 # Copying the value should not be necessary, but we record
1610 # it just in case.
1611 ptc.gainUnadjusted[ampName] = ptc.gain[ampName]
1612 ptc.gain[ampName] = newGain
1613 else:
1614 log.warning("Cannot apply ampOffsetGainRatioFixup with fewer than 2 good amplifiers.")
1617class FlatGradientFitter:
1618 """Class to fit various large-scale flat-field gradients.
1620 This fitter will take arrays of x/y/value and (by default) fit a radial
1621 gradient, using a spline function at the specified nodes. The fitter
1622 will also fit out a nuisance parameter for the ratio of the ITL/E2V
1623 throughput (though in general this could work for any such focal plane).
1624 The focal plane origin is set by fp_centroid_x and fp_centroid_y which
1625 is the center of the radial gradient, though it may be modified if
1626 fit_centroid is True. In addition, the fitter may fit a linear gradient
1627 in x/y if fit_gradient is True. The "pivot" of the gradient is at
1628 fp_centroid_x/fp_centroid_y.
1630 Parameters
1631 ----------
1632 nodes : `np.ndarray`
1633 Array of spline nodes for radial spline.
1634 x : `np.ndarray`
1635 Array of x values for points to fit (mm).
1636 y : `np.ndarray`
1637 Array of y values for points to fit (mm).
1638 value : `np.ndarray`
1639 Array of values describing the flat field to fit the gradient.
1640 itl_indices : `np.ndarray`
1641 Array of indices corresponding to the ITL detectors.
1642 constrain_zero : `bool`, optional
1643 Constrain the outermost radial spline value to 0?
1644 fit_centroid : `bool`, optional
1645 Fit an additional centroid offset?
1646 fit_gradient : `bool`, optional
1647 Fit an additional plane gradient?
1648 fp_centroid_x : `float`, optional
1649 Focal plane centroid x (mm).
1650 fp_centroid_y : `float`, optional
1651 Focal plane centroid y (mm).
1652 """
1653 def __init__(
1654 self,
1655 nodes,
1656 x,
1657 y,
1658 value,
1659 itl_indices,
1660 constrain_zero=True,
1661 fit_centroid=False,
1662 fit_gradient=False,
1663 fp_centroid_x=0.0,
1664 fp_centroid_y=0.0
1665 ):
1666 self._nodes = nodes
1667 self._x = x
1668 self._y = y
1669 self._value = value
1670 self._itl_indices = itl_indices
1671 self._fp_centroid_x = fp_centroid_x
1672 self._fp_centroid_y = fp_centroid_y
1674 self._constrain_zero = constrain_zero
1676 self._fit_centroid = fit_centroid
1677 self._fit_gradient = fit_gradient
1679 self.indices = {"spline": np.arange(len(nodes))}
1680 npar = len(nodes)
1682 self._fit_itl_ratio = False
1683 if len(itl_indices) > 0:
1684 self._fit_itl_ratio = True
1685 self.indices["itl_ratio"] = npar
1686 npar += 1
1688 radius = np.sqrt((self._x - self._fp_centroid_x)**2. + (self._y - self._fp_centroid_y)**2.)
1690 if fit_centroid:
1691 self.indices["centroid_delta"] = np.arange(2) + npar
1692 npar += 2
1693 else:
1694 self._radius = radius
1696 if fit_gradient:
1697 self.indices["gradient"] = np.arange(2) + npar
1698 npar += 2
1700 self._npar = npar
1702 self._bounds = [(None, None)]*npar
1704 if self._constrain_zero:
1705 self._bounds[self.indices["spline"][-1]] = (0.0, 0.0)
1707 def compute_p0(self, itl_ratio=None):
1708 """Compute initial guess for fit parameters.
1710 Returns
1711 -------
1712 pars : `np.ndarray`
1713 Array of first guess fit parameters.
1714 """
1715 pars = np.zeros(self._npar)
1716 value = self._value.copy()
1718 if itl_ratio is not None and self._fit_itl_ratio:
1719 value[self._itl_indices] /= itl_ratio
1721 # Initial spline values
1722 radius = np.sqrt((self._x - self._fp_centroid_x)**2. + (self._y - self._fp_centroid_y)**2.)
1723 for i, index in enumerate(self.indices["spline"]):
1724 if i == 0:
1725 low = self._nodes[i]
1726 else:
1727 low = (self._nodes[i - 1] + self._nodes[i])/2.
1728 if i == (len(self._nodes) - 1):
1729 high = self._nodes[i]
1730 else:
1731 high = (self._nodes[i] + self._nodes[i + 1])/2.
1732 u = ((radius > low) & (radius < high))
1733 if u.sum() == 0:
1734 pars[index] = 0.0
1735 else:
1736 pars[index] = np.median(value[u])
1738 if self._constrain_zero:
1739 pars[self.indices["spline"][-1]] = 0.0
1741 spl = Akima1DInterpolator(self._nodes, pars[self.indices["spline"]], method="akima")
1742 model = spl(radius)
1743 resid_ratio = value / model
1745 if self._fit_itl_ratio:
1746 if itl_ratio is not None:
1747 pars[self.indices["itl_ratio"]] = itl_ratio
1748 else:
1749 e2v_indices = np.delete(np.arange(len(self._value)), self._itl_indices)
1751 itl_inner = radius[self._itl_indices] < 0.8*np.max(radius)
1752 e2v_inner = radius[e2v_indices] < 0.8*np.max(radius)
1754 itl_median = np.nanmedian(resid_ratio[self._itl_indices][itl_inner])
1755 e2v_median = np.nanmedian(resid_ratio[e2v_indices][e2v_inner])
1757 pars[self.indices["itl_ratio"]] = itl_median / e2v_median
1759 if self._fit_centroid:
1760 pars[self.indices["centroid_delta"]] = [0.0, 0.0]
1762 if self._fit_gradient:
1763 resid_ratio = self._value / model
1765 ok = (np.isfinite(resid_ratio) & (model > 0.5))
1767 fit = np.polyfit(self._x[ok], resid_ratio[ok], 1)
1768 pars[self.indices["gradient"][0]] = fit[0]
1770 fit = np.polyfit(self._y[ok], resid_ratio[ok], 1)
1771 pars[self.indices["gradient"][1]] = fit[0]
1773 return pars
1775 def fit(self, p0, freeze_itl_ratio=False, fit_eps=1e-8, fit_gtol=1e-10):
1776 """Do a non-linear minimization to fit the parameters.
1778 Parameters
1779 ----------
1780 p0 : `np.ndarray`
1781 Array of initial parameter estimates.
1782 freeze_itl_ratio : `bool`, optional
1783 Freeze the ITL ratio in the fit?
1784 fit_eps : `float`, optional
1785 Value of ``eps`` to send to the scipy minimizer.
1786 fit_gtol : `float`, optional
1787 Value of ``gtol`` to send to the scipy minimizer.
1789 Returns
1790 -------
1791 pars : `np.ndarray`
1792 Array of parameters. Use ``fitter.indices`` for the
1793 dictionary to map parameters to subsets.
1794 """
1795 bounds = self._bounds
1796 if freeze_itl_ratio and self._fit_itl_ratio:
1797 ind = self.indices["itl_ratio"]
1798 par = p0[ind]
1799 bounds[ind] = (par, par)
1801 res = minimize(
1802 self,
1803 p0,
1804 method="L-BFGS-B",
1805 jac=False,
1806 bounds=bounds,
1807 options={
1808 "maxfun": 10000,
1809 "maxiter": 10000,
1810 "maxcor": 20,
1811 "eps": fit_eps,
1812 "gtol": fit_gtol,
1813 },
1814 callback=None,
1815 )
1816 pars = res.x
1818 return pars
1820 def compute_model(self, pars):
1821 """Compute the model given a set of parameters.
1823 Parameters
1824 ----------
1825 pars : `np.ndarray`
1826 Parameter array to compute model.
1828 Returns
1829 -------
1830 model_array : `np.ndarray`
1831 Array of model parameters at the input x/y.
1832 """
1833 spl = Akima1DInterpolator(self._nodes, pars[self.indices["spline"]], method="akima")
1834 if self._fit_centroid:
1835 centroid_delta = pars[self.indices["centroid_delta"]]
1836 centroid_x = self._fp_centroid_x + centroid_delta[0]
1837 centroid_y = self._fp_centroid_y + centroid_delta[1]
1838 radius = np.sqrt((self._x - centroid_x)**2. + (self._y - centroid_y)**2.)
1839 else:
1840 radius = self._radius
1842 model = spl(np.clip(radius, self._nodes[0], self._nodes[-1]))
1843 if self._fit_itl_ratio:
1844 model[self._itl_indices] *= pars[self.indices["itl_ratio"]]
1846 if self._fit_gradient:
1847 a, b = pars[self.indices["gradient"]]
1848 gradient = 1 + a*(self._x - self._fp_centroid_x) + b*(self._y - self._fp_centroid_y)
1849 model /= gradient
1851 return model
1853 def __call__(self, pars):
1854 """Compute the cost function for a set of parameters.
1856 Parameters
1857 ----------
1858 pars : `np.ndarray`
1859 Parameter array to compute model.
1861 Returns
1862 -------
1863 cost : `float`
1864 Cost value computed from the absolute deviation.
1865 """
1867 model = self.compute_model(pars)
1869 absdev = np.abs(self._value - model)
1870 t = np.sum(absdev.astype(np.float64))
1872 return t
1875class FlatGainRatioFitter:
1876 """Class to fit amplifier gain ratios, removing a simple gradient.
1878 This fitter will take arrays of x/y/amp_num/value and fit amplifier
1879 gain ratios, using one amplifier as the fixed point. The fitter
1880 uses a low-order chebyshev polynomial to fit out the gradient, with
1881 amp ratios on top of this.
1883 Parameters
1884 ----------
1885 bbox : `lsst.geom.Box2I`
1886 Bounding box for Chebyshev polynomial gradient.
1887 order : `int`
1888 Order of Chebyshev polynomial.
1889 x : `np.ndarray`
1890 Array of x values for points to fit (detector pixels).
1891 y : `np.ndarray`
1892 Array of y values for points to fit (detector pixels).
1893 amp_index : `np.ndarray`
1894 Array of amp numbers associated with each x/y pair.
1895 value : `np.ndarray`
1896 Flat value at each x/y pair.
1897 amps : `np.ndarray`
1898 Array of unique amplifier numbers that will be parameterized.
1899 Any of these amps that does not have any data with the same
1900 amp_index will be set to 1.0.
1901 fixed_amp_index : `int`
1902 Amplifier number to keep fixed.
1903 """
1904 def __init__(self, bbox, order, x, y, amp_index, value, amps, fixed_amp_index):
1905 self._bbox = bbox
1906 self._order = order
1907 self._x = x.astype(np.float64)
1908 self._y = y.astype(np.float64)
1909 self._amp_index = amp_index
1910 self._value = value.astype(np.float64)
1911 self._fixed_amp_index = fixed_amp_index
1913 self.indices = {"chebyshev": np.arange((order + 1) * (order + 1))}
1914 npar = len(self.indices["chebyshev"])
1916 self._amps = amps
1917 self._n_amp = len(self._amps)
1919 self.indices["amp_pars"] = np.arange(self._n_amp) + npar
1920 npar += self._n_amp
1922 self._npar = npar
1924 self._amp_indices = {}
1925 for i in range(self._n_amp):
1926 amp_index = self._amps[i]
1927 self._amp_indices[amp_index] = (self._amp_index == amp_index)
1929 def fit(self, n_iter=10, nsig_clip=5.0):
1930 """Fit the amp ratio parameters.
1932 This uses an iterative fit, where it fits a Chebyshev gradient,
1933 computes amp ratios, and re-fits the gradient.
1935 Parameters
1936 ----------
1937 n_iter : `int`, optional
1938 Number of iterations for fit.
1939 nsig_clip : `float`, optional
1940 Number of sigma in gain correction distribution to clip when
1941 fitting out the Chebyshev gradient.
1943 Returns
1944 -------
1945 pars : `np.ndarray`
1946 Chebyshev parameters and amp offset parameters.
1947 """
1948 value = self._value.copy()
1950 control = lsst.afw.math.ChebyshevBoundedFieldControl()
1951 control.orderX = self._order
1952 control.orderY = self._order
1953 control.triangular = False
1955 pars = np.zeros(self._npar)
1957 pars[self.indices["amp_pars"]] = 1.0
1959 field_use = np.ones(len(value), dtype=np.bool_)
1961 for i in range(n_iter):
1962 field = lsst.afw.math.ChebyshevBoundedField.fit(
1963 self._bbox,
1964 self._x[field_use],
1965 self._y[field_use],
1966 value[field_use],
1967 control,
1968 )
1970 pars[self.indices["chebyshev"][:]] = field.getCoefficients().ravel()
1972 ratio = self._value / field.evaluate(self._x, self._y)
1974 fixed_med = np.median(ratio[self._amp_indices[self._fixed_amp_index]])
1975 ratio /= fixed_med
1977 pars[self.indices["amp_pars"][self._fixed_amp_index]] = 1.0
1979 value = self._value.copy()
1981 for j in range(self._n_amp):
1982 amp_index = self._amps[j]
1984 if np.sum(self._amp_indices[amp_index]) == 0:
1985 continue
1987 pars[self.indices["amp_pars"][j]] = np.median(ratio[self._amp_indices[amp_index]])
1989 amp_pars = pars[self.indices["amp_pars"]]
1990 med = np.median(amp_pars)
1991 sig = median_abs_deviation(amp_pars, scale="normal")
1992 amp_gradient_bad, = np.where(np.abs(amp_pars - med) > nsig_clip * sig)
1993 field_use[:] = True
1994 for agb in amp_gradient_bad:
1995 field_use[self._amp_indices[agb]] = False
1997 return pars
1999 def compute_model(self, pars):
2000 """Compute the gradient/amp ratio model for a given set of parameters.
2002 Parameters
2003 ----------
2004 pars : `np.ndarray`
2005 Chebyshev parameters and amp offset parameters.
2007 Returns
2008 -------
2009 model : `np.ndarray`
2010 The model at each x/y pair.
2011 """
2012 field = lsst.afw.math.ChebyshevBoundedField(
2013 self._bbox,
2014 pars[self.indices["chebyshev"]].reshape(self._order + 1, self._order + 1),
2015 )
2016 model = field.evaluate(self._x, self._y)
2018 for i in range(self._n_amp):
2019 amp_index = self._amps[i]
2021 model[self._amp_indices[amp_index]] *= pars[self.indices["amp_pars"][i]]
2023 return model
2025 def __call__(self, pars):
2026 """Compute the cost function for a set of parameters.
2028 Parameters
2029 ----------
2030 pars : `np.ndarray`
2031 Chebyshev parameters and amp offset parameters.
2033 Returns
2034 -------
2035 cost : `float`
2036 Cost value computed from the absolute deviation.
2037 """
2038 model = self.compute_model(pars)
2040 absdev = np.abs(self._value - model)
2041 t = np.sum(absdev.astype(np.float64))
2043 return t
2046def bin_focal_plane(
2047 exposure_handle_dict,
2048 detector_boundary,
2049 bin_factor,
2050 defect_handle_dict={},
2051 include_itl_flag=True,
2052):
2053 """Bin all the detectors into the full focal plane.
2055 This function reads in images; takes a simple average if there
2056 are more than one input per detector; excludes detector edges;
2057 and bins according to the bin factor. The output is a struct
2058 with focal plane coordinates, detector numbers, and a flag
2059 if the detector is an ITL detector.
2061 Parameters
2062 ----------
2063 exposure_handle_dict : `dict`
2064 Dict keyed by detector (`int`), each element is a list
2065 of `lsst.daf.butler.DeferredDatasetHandle` that will be averaged.
2066 detector_boundary : `int`
2067 Boundary of the detector to ignore (pixels).
2068 bin_factor : `int`
2069 Binning factor. Detectors will be cropped (after applying the
2070 ``detector_boundary``) such that there are no partially
2071 covered binned pixels.
2072 defect_handle_dict : `dict`, optional
2073 Dict keyed by detector (`int`), each element is a
2074 `lsst.daf.butler.DeferredDatasetHandle` for defects to be applied.
2075 include_itl_flag : `bool`, optional
2076 Include a flag for which detectors are ITL?
2078 Returns
2079 -------
2080 binned : `astropy.table.Table`
2081 Table with focal plane positions at the center of each bin
2082 (``xf``, ``yf``); average image values (``value``); and detector
2083 number (``detector``).
2084 """
2085 xf_arrays = []
2086 yf_arrays = []
2087 value_arrays = []
2088 detector_arrays = []
2089 itl_arrays = []
2091 for det in exposure_handle_dict.keys():
2092 flat = exposure_handle_dict[det].get()
2093 defect_handle = defect_handle_dict.get(det, None)
2094 if defect_handle is not None:
2095 defects = defect_handle.get()
2096 else:
2097 defects = None
2099 detector = flat.getDetector()
2101 # Mask out defects if we have them.
2102 if defects is not None:
2103 for defect in defects:
2104 flat.image[defect.getBBox()].array[:, :] = np.nan
2106 # Mask NO_DATA pixels if we have them.
2107 no_data = ((flat.mask.array[:, :] & flat.mask.getPlaneBitMask("NO_DATA")) > 0)
2108 flat.image.array[no_data] = np.nan
2110 # Bin the image, avoiding the boundary and the masked pixels.
2111 # We also make sure we are using an integral number of
2112 # steps to avoid partially covered binned pixels.
2114 arr = flat.image.array
2116 n_step_y = (arr.shape[0] - (2 * detector_boundary)) // bin_factor
2117 y_min = detector_boundary
2118 y_max = bin_factor * n_step_y + y_min
2119 n_step_x = (arr.shape[1] - (2 * detector_boundary)) // bin_factor
2120 x_min = detector_boundary
2121 x_max = bin_factor * n_step_x + x_min
2123 arr = arr[y_min: y_max, x_min: x_max]
2124 binned = arr.reshape((n_step_y, bin_factor, n_step_x, bin_factor))
2125 with warnings.catch_warnings():
2126 warnings.filterwarnings("ignore", r"Mean of empty")
2127 binned = np.nanmean(binned, axis=1)
2128 binned = np.nanmean(binned, axis=2)
2130 xx = np.arange(binned.shape[1]) * bin_factor + bin_factor / 2. + x_min
2131 yy = np.arange(binned.shape[0]) * bin_factor + bin_factor / 2. + y_min
2132 x, y = np.meshgrid(xx, yy)
2133 x = x.ravel()
2134 y = y.ravel()
2135 value = binned.ravel()
2137 # Transform to focal plane coordinates.
2138 transform = detector.getTransform(lsst.afw.cameraGeom.PIXELS, lsst.afw.cameraGeom.FOCAL_PLANE)
2139 xy = np.vstack((x, y))
2140 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)
2141 xf = xf.ravel()
2142 yf = yf.ravel()
2144 if include_itl_flag:
2145 is_itl = np.zeros(len(value), dtype=np.bool_)
2146 # We use this check so that ITL matches ITL science detectors,
2147 # ITL_WF wavefront detectors, and pseudoITL test detectors.
2148 is_itl[:] = ("ITL" in detector.getPhysicalType())
2150 xf_arrays.append(xf)
2151 yf_arrays.append(yf)
2152 value_arrays.append(value)
2153 detector_arrays.append(np.full_like(xf, det, dtype=np.int32))
2154 if include_itl_flag:
2155 itl_arrays.append(is_itl)
2157 xf = np.concatenate(xf_arrays)
2158 yf = np.concatenate(yf_arrays)
2159 value = np.concatenate(value_arrays)
2160 detector = np.concatenate(detector_arrays)
2162 binned = Table(
2163 np.zeros(
2164 len(xf),
2165 dtype=[
2166 ("xf", "f8"),
2167 ("yf", "f8"),
2168 ("value", "f8"),
2169 ("detector", "i4"),
2170 ],
2171 )
2172 )
2173 binned["xf"] = xf
2174 binned["yf"] = yf
2175 binned["value"] = value
2176 binned["detector"] = detector
2178 if include_itl_flag:
2179 binned["itl"] = np.concatenate(itl_arrays).astype(np.bool_)
2181 return binned
2184def bin_flat(ptc, exposure, bin_factor=8, amp_boundary=20, apply_gains=True, gain_ratios=None):
2185 """Bin a flat image, being careful with amplifier edges.
2187 This will optionally apply gains, and apply any gain
2188 ratios.
2190 Parameters
2191 ----------
2192 ptc : `lsst.ip.isr.PhotonTransferCurveDatasets`
2193 PTC dataset with relevant gains.
2194 exposure : `lsst.afw.image.Exposure`
2195 Exposure to bin.
2196 bin_factor : `int`, optional
2197 Binning factor.
2198 amp_boundary : `int`, optional
2199 Boundary around each amp to ignore in binning.
2200 apply_gains : `bool`, optional
2201 Apply gains before binning?
2202 gain_ratios : `np.ndarray`, optional
2203 Array of gain ratios to apply.
2205 Returns
2206 -------
2207 binned : `astropy.table.Table`
2208 Table with detector coordinates at the center of each bin
2209 (``xd``, ``yd``); average image values (``value``); and amplifier
2210 index (``amp_index``).
2211 """
2212 detector = exposure.getDetector()
2214 for i, amp_name in enumerate(ptc.ampNames):
2215 bbox = detector[amp_name].getBBox()
2216 if amp_name in ptc.badAmps:
2217 exposure[bbox].image.array[:, :] = np.nan
2218 continue
2220 if apply_gains:
2221 exposure[bbox].image.array[:, :] *= ptc.gainUnadjusted[amp_name]
2222 if gain_ratios is not None:
2223 exposure[bbox].image.array[:, :] /= gain_ratios[i]
2225 # Next we bin the detector, avoiding amp edges.
2226 xd_arrays = []
2227 yd_arrays = []
2228 value_arrays = []
2229 amp_arrays = []
2231 for i, amp in enumerate(detector):
2232 arr = exposure[amp.getBBox()].image.array
2234 n_step_y = (arr.shape[0] - (2 * amp_boundary)) // bin_factor
2235 y_min = amp_boundary
2236 y_max = bin_factor * n_step_y + y_min
2237 n_step_x = (arr.shape[1] - (2 * amp_boundary)) // bin_factor
2238 x_min = amp_boundary
2239 x_max = bin_factor * n_step_x + x_min
2241 arr = arr[y_min: y_max, x_min: x_max]
2242 binned = arr.reshape((n_step_y, bin_factor, n_step_x, bin_factor))
2243 with warnings.catch_warnings():
2244 warnings.filterwarnings("ignore", r"Mean of empty")
2245 binned = np.nanmean(binned, axis=1)
2246 binned = np.nanmean(binned, axis=2)
2248 xx = np.arange(binned.shape[1]) * bin_factor + bin_factor / 2. + x_min
2249 yy = np.arange(binned.shape[0]) * bin_factor + bin_factor / 2. + y_min
2250 x, y = np.meshgrid(xx, yy)
2251 x = x.ravel() + amp.getBBox().getBeginX()
2252 y = y.ravel() + amp.getBBox().getBeginY()
2253 value = binned.ravel()
2255 xd_arrays.append(x)
2256 yd_arrays.append(y)
2257 value_arrays.append(value)
2258 amp_arrays.append(np.full(len(x), i))
2260 xd = np.concatenate(xd_arrays)
2261 yd = np.concatenate(yd_arrays)
2262 value = np.concatenate(value_arrays)
2263 amp_index = np.concatenate(amp_arrays)
2265 binned = np.zeros(
2266 len(xd),
2267 dtype=[
2268 ("xd", "f8"),
2269 ("yd", "f8"),
2270 ("value", "f8"),
2271 ("amp_index", "i4"),
2272 ],
2273 )
2274 binned["xd"] = xd
2275 binned["yd"] = yd
2276 binned["value"] = value
2277 binned["amp_index"] = amp_index
2279 return Table(data=binned)