Coverage for python/lsst/cp/pipe/utils.py: 7%

1# This file is part of cp_pipe.

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.

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__ = [

24 "ddict2dict",

25 "CovFastFourierTransform",

26 "getReadNoise",

27 "ampOffsetGainRatioFixup",

28 "ElectrostaticFit",

29 "BoundaryShifts",

30 "ElectrostaticCcdGeom",

31]

33from astropy.table import Table

34import copy

35import galsim

36import logging

37import numpy as np

38import itertools

39import numpy.polynomial.polynomial as poly

40import warnings

42from scipy.interpolate import Akima1DInterpolator

43from scipy.optimize import leastsq, minimize

44from scipy.stats import median_abs_deviation, norm

45from lmfit import Minimizer

47from lsst.ip.isr import isrMock

48import lsst.afw.cameraGeom

49import lsst.afw.image

50import lsst.afw.math

52from numpy.polynomial.legendre import leggauss

55def sigmaClipCorrection(nSigClip):

56 """Correct measured sigma to account for clipping.

58 If we clip our input data and then measure sigma, then the

59 measured sigma is smaller than the true value because real

60 points beyond the clip threshold have been removed. This is a

61 small (1.5% at nSigClip=3) effect when nSigClip >~ 3, but the

62 default parameters for measure crosstalk use nSigClip=2.0.

63 This causes the measured sigma to be about 15% smaller than

64 real. This formula corrects the issue, for the symmetric case

65 (upper clip threshold equal to lower clip threshold).

67 Parameters

68 ----------

69 nSigClip : `float`

70 Number of sigma the measurement was clipped by.

72 Returns

73 -------

74 scaleFactor : `float`

75 Scale factor to increase the measured sigma by.

76 """

77 varFactor = 1.0 - (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))

78 return 1.0 / np.sqrt(varFactor)

81def calculateWeightedReducedChi2(measured, model, weightsMeasured, nData, nParsModel):

82 """Calculate weighted reduced chi2.

84 Parameters

85 ----------

86 measured : `list`

87 List with measured data.

88 model : `list`

89 List with modeled data.

90 weightsMeasured : `list`

91 List with weights for the measured data.

92 nData : `int`

93 Number of data points.

94 nParsModel : `int`

95 Number of parameters in the model.

97 Returns

98 -------

99 redWeightedChi2 : `float`

100 Reduced weighted chi2.

101 """

102 wRes = (measured - model)*weightsMeasured

103 return ((wRes*wRes).sum())/(nData-nParsModel)

104

105

106def makeMockFlats(expTime, gain=1.0, readNoiseElectrons=5, fluxElectrons=1000,

107 randomSeedFlat1=1984, randomSeedFlat2=666, powerLawBfParams=[],

108 expId1=0, expId2=1, ampNames=[]):

109 """Create a pair or mock flats with isrMock.

110

111 Parameters

112 ----------

113 expTime : `float`

114 Exposure time of the flats.

115 gain : `float`, optional

116 Gain, in e/ADU.

117 readNoiseElectrons : `float`, optional

118 Read noise rms, in electrons.

119 fluxElectrons : `float`, optional

120 Flux of flats, in electrons per second.

121 randomSeedFlat1 : `int`, optional

122 Random seed for the normal distrubutions for the mean signal

123 and noise (flat1).

124 randomSeedFlat2 : `int`, optional

125 Random seed for the normal distrubutions for the mean signal

126 and noise (flat2).

127 powerLawBfParams : `list`, optional

128 Parameters for `galsim.cdmodel.PowerLawCD` to simulate the

129 brightter-fatter effect.

130 expId1 : `int`, optional

131 Exposure ID for first flat.

132 expId2 : `int`, optional

133 Exposure ID for second flat.

134 ampNames : `list` [`str`], optional

135 Names of amplifiers for filling in header information.

136

137 Returns

138 -------

139 flatExp1 : `lsst.afw.image.exposure.ExposureF`

140 First exposure of flat field pair.

141 flatExp2 : `lsst.afw.image.exposure.ExposureF`

142 Second exposure of flat field pair.

143

144 Notes

145 -----

146 The parameters of `galsim.cdmodel.PowerLawCD` are `n, r0, t0, rx,

147 tx, r, t, alpha`. For more information about their meaning, see

148 the Galsim documentation

149 https://galsim-developers.github.io/GalSim/_build/html/_modules/galsim/cdmodel.html # noqa: W505

150 and Gruen+15 (1501.02802).

151

152 Example: galsim.cdmodel.PowerLawCD(8, 1.1e-7, 1.1e-7, 1.0e-8,

153 1.0e-8, 1.0e-9, 1.0e-9, 2.0)

154 """

155 flatFlux = fluxElectrons # e/s

156 flatMean = flatFlux*expTime # e

157 readNoise = readNoiseElectrons # e

158

159 mockImageConfig = isrMock.IsrMock.ConfigClass()

160

161 mockImageConfig.flatDrop = 0.99999

162 mockImageConfig.isTrimmed = True

163

164 flatExp1 = isrMock.FlatMock(config=mockImageConfig).run()

165 flatExp2 = flatExp1.clone()

166 (shapeY, shapeX) = flatExp1.getDimensions()

167 flatWidth = np.sqrt(flatMean)

168

169 rng1 = np.random.RandomState(randomSeedFlat1)

170 flatData1 = rng1.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng1.normal(0.0, readNoise,

171 (shapeX, shapeY))

172 rng2 = np.random.RandomState(randomSeedFlat2)

173 flatData2 = rng2.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng2.normal(0.0, readNoise,

174 (shapeX, shapeY))

175 # Simulate BF with power law model in galsim

176 if len(powerLawBfParams):

177 if not len(powerLawBfParams) == 8:

178 raise RuntimeError("Wrong number of parameters for `galsim.cdmodel.PowerLawCD`. "

179 f"Expected 8; passed {len(powerLawBfParams)}.")

180 cd = galsim.cdmodel.PowerLawCD(*powerLawBfParams)

181 tempFlatData1 = galsim.Image(flatData1)

182 temp2FlatData1 = cd.applyForward(tempFlatData1)

183

184 tempFlatData2 = galsim.Image(flatData2)

185 temp2FlatData2 = cd.applyForward(tempFlatData2)

186

187 flatExp1.image.array[:] = temp2FlatData1.array/gain # ADU

188 flatExp2.image.array[:] = temp2FlatData2.array/gain # ADU

189 else:

190 flatExp1.image.array[:] = flatData1/gain # ADU

191 flatExp2.image.array[:] = flatData2/gain # ADU

192

193 visitInfoExp1 = lsst.afw.image.VisitInfo(exposureTime=expTime)

194 visitInfoExp2 = lsst.afw.image.VisitInfo(exposureTime=expTime)

195

196 flatExp1.info.id = expId1

197 flatExp1.getInfo().setVisitInfo(visitInfoExp1)

198 flatExp2.info.id = expId2

199 flatExp2.getInfo().setVisitInfo(visitInfoExp2)

200

201 # Set ISR metadata.

202 flatExp1.metadata["LSST ISR UNITS"] = "adu"

203 flatExp2.metadata["LSST ISR UNITS"] = "adu"

204 for ampName in ampNames:

205 key = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"

206 value = readNoiseElectrons / gain

207

208 flatExp1.metadata[key] = value

209 flatExp2.metadata[key] = value

210

211 key = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"

212 flatExp1.metadata[key] = 25000.0 # adu

213 flatExp2.metadata[key] = 25000.0 # adu

214

215 key = f"LSST ISR AMPOFFSET PEDESTAL {ampName}"

216 value = 0.0

217

218 flatExp1.metadata[key] = value

219 flatExp2.metadata[key] = value

220

221 return flatExp1, flatExp2

222

223

224def irlsFit(initialParams, dataX, dataY, function, weightsY=None, weightType='Cauchy', scaleResidual=True):

225 """Iteratively reweighted least squares fit.

226

227 This uses the `lsst.cp.pipe.utils.fitLeastSq`, but applies weights

228 based on the Cauchy distribution by default. Other weight options

229 are implemented. See e.g. Holland and Welsch, 1977,

230 doi:10.1080/03610927708827533

231

232 Parameters

233 ----------

234 initialParams : `list` [`float`]

235 Starting parameters.

236 dataX : `numpy.array`, (N,)

237 Abscissa data.

238 dataY : `numpy.array`, (N,)

239 Ordinate data.

240 function : callable

241 Function to fit.

242 weightsY : `numpy.array`, (N,)

243 Weights to apply to the data.

244 weightType : `str`, optional

245 Type of weighting to use. One of Cauchy, Anderson, bisquare,

246 box, Welsch, Huber, logistic, or Fair.

247 scaleResidual : `bool`, optional

248 If true, the residual is scaled by the sqrt of the Y values.

249

250 Returns

251 -------

252 polyFit : `list` [`float`]

253 Final best fit parameters.

254 polyFitErr : `list` [`float`]

255 Final errors on fit parameters.

256 chiSq : `float`

257 Reduced chi squared.

258 weightsY : `list` [`float`]

259 Final weights used for each point.

260

261 Raises

262 ------

263 RuntimeError :

264 Raised if an unknown weightType string is passed.

265 """

266 if not weightsY:

267 weightsY = np.ones_like(dataX)

268

269 polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)

270 for iteration in range(10):

271 resid = np.abs(dataY - function(polyFit, dataX))

272 if scaleResidual:

273 resid = resid / np.sqrt(dataY)

274 if weightType == 'Cauchy':

275 # Use Cauchy weighting. This is a soft weight.

276 # At [2, 3, 5, 10] sigma, weights are [.59, .39, .19, .05].

277 Z = resid / 2.385

278 weightsY = 1.0 / (1.0 + np.square(Z))

279 elif weightType == 'Anderson':

280 # Anderson+1972 weighting. This is a hard weight.

281 # At [2, 3, 5, 10] sigma, weights are [.67, .35, 0.0, 0.0].

282 Z = resid / (1.339 * np.pi)

283 weightsY = np.where(Z < 1.0, np.sinc(Z), 0.0)

284 elif weightType == 'bisquare':

285 # Beaton and Tukey (1974) biweight. This is a hard weight.

286 # At [2, 3, 5, 10] sigma, weights are [.81, .59, 0.0, 0.0].

287 Z = resid / 4.685

288 weightsY = np.where(Z < 1.0, 1.0 - np.square(Z), 0.0)

289 elif weightType == 'box':

290 # Hinich and Talwar (1975). This is a hard weight.

291 # At [2, 3, 5, 10] sigma, weights are [1.0, 0.0, 0.0, 0.0].

292 weightsY = np.where(resid < 2.795, 1.0, 0.0)

293 elif weightType == 'Welsch':

294 # Dennis and Welsch (1976). This is a hard weight.

295 # At [2, 3, 5, 10] sigma, weights are [.64, .36, .06, 1e-5].

296 Z = resid / 2.985

297 weightsY = np.exp(-1.0 * np.square(Z))

298 elif weightType == 'Huber':

299 # Huber (1964) weighting. This is a soft weight.

300 # At [2, 3, 5, 10] sigma, weights are [.67, .45, .27, .13].

301 Z = resid / 1.345

302 weightsY = np.where(Z < 1.0, 1.0, 1 / Z)

303 elif weightType == 'logistic':

304 # Logistic weighting. This is a soft weight.

305 # At [2, 3, 5, 10] sigma, weights are [.56, .40, .24, .12].

306 Z = resid / 1.205

307 weightsY = np.tanh(Z) / Z

308 elif weightType == 'Fair':

309 # Fair (1974) weighting. This is a soft weight.

310 # At [2, 3, 5, 10] sigma, weights are [.41, .32, .22, .12].

311 Z = resid / 1.4

312 weightsY = (1.0 / (1.0 + (Z)))

313 else:

314 raise RuntimeError(f"Unknown weighting type: {weightType}")

315 polyFit, polyFitErr, chiSq = fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)

316

317 return polyFit, polyFitErr, chiSq, weightsY

318

319

320def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None):

321 """Do a fit and estimate the parameter errors using using

322 scipy.optimize.leastq.

323

324 optimize.leastsq returns the fractional covariance matrix. To

325 estimate the standard deviation of the fit parameters, multiply

326 the entries of this matrix by the unweighted reduced chi squared

327 and take the square root of the diagonal elements.

328

329 Parameters

330 ----------

331 initialParams : `list` [`float`]

332 initial values for fit parameters.

333 dataX : `numpy.array`, (N,)

334 Data in the abscissa axis.

335 dataY : `numpy.array`, (N,)

336 Data in the ordinate axis.

337 function : callable object (function)

338 Function to fit the data with.

339 weightsY : `numpy.array`, (N,)

340 Weights of the data in the ordinate axis.

341

342 Return

343 ------

344 pFitSingleLeastSquares : `list` [`float`]

345 List with fitted parameters.

346 pErrSingleLeastSquares : `list` [`float`]

347 List with errors for fitted parameters.

348

349 reducedChiSqSingleLeastSquares : `float`

350 Reduced chi squared, unweighted if weightsY is not provided.

351 """

352 if weightsY is None:

353 weightsY = np.ones(len(dataX))

354

355 def errFunc(p, x, y, weightsY=None):

356 if weightsY is None:

357 weightsY = np.ones(len(x))

358 return (function(p, x) - y)*weightsY

359

360 pFit, pCov, infoDict, errMessage, success = leastsq(errFunc, initialParams,

361 args=(dataX, dataY, weightsY), full_output=1,

362 epsfcn=0.0001)

363

364 if (len(dataY) > len(initialParams)) and pCov is not None:

365 reducedChiSq = calculateWeightedReducedChi2(dataY, function(pFit, dataX), weightsY, len(dataY),

366 len(initialParams))

367 pCov *= reducedChiSq

368 else:

369 pCov = np.zeros((len(initialParams), len(initialParams)))

370 pCov[:, :] = np.nan

371 reducedChiSq = np.nan

372

373 errorVec = []

374 for i in range(len(pFit)):

375 errorVec.append(np.fabs(pCov[i][i])**0.5)

376

377 pFitSingleLeastSquares = pFit

378 pErrSingleLeastSquares = np.array(errorVec)

379

380 return pFitSingleLeastSquares, pErrSingleLeastSquares, reducedChiSq

381

382

383def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.):

384 """Do a fit using least squares and bootstrap to estimate parameter errors.

385

386 The bootstrap error bars are calculated by fitting 100 random data sets.

387

388 Parameters

389 ----------

390 initialParams : `list` [`float`]

391 initial values for fit parameters.

392 dataX : `numpy.array`, (N,)

393 Data in the abscissa axis.

394 dataY : `numpy.array`, (N,)

395 Data in the ordinate axis.

396 function : callable object (function)

397 Function to fit the data with.

398 weightsY : `numpy.array`, (N,), optional.

399 Weights of the data in the ordinate axis.

400 confidenceSigma : `float`, optional.

401 Number of sigmas that determine confidence interval for the

402 bootstrap errors.

403

404 Return

405 ------

406 pFitBootstrap : `list` [`float`]

407 List with fitted parameters.

408 pErrBootstrap : `list` [`float`]

409 List with errors for fitted parameters.

410 reducedChiSqBootstrap : `float`

411 Reduced chi squared, unweighted if weightsY is not provided.

412 """

413 if weightsY is None:

414 weightsY = np.ones(len(dataX))

415

416 def errFunc(p, x, y, weightsY):

417 if weightsY is None:

418 weightsY = np.ones(len(x))

419 return (function(p, x) - y)*weightsY

420

421 # Fit first time

422 pFit, _ = leastsq(errFunc, initialParams, args=(dataX, dataY, weightsY), full_output=0)

423

424 # Get the stdev of the residuals

425 residuals = errFunc(pFit, dataX, dataY, weightsY)

426 # 100 random data sets are generated and fitted

427 pars = []

428 for i in range(100):

429 randomDelta = np.random.normal(0., np.fabs(residuals), len(dataY))

430 randomDataY = dataY + randomDelta

431 randomFit, _ = leastsq(errFunc, initialParams,

432 args=(dataX, randomDataY, weightsY), full_output=0)

433 pars.append(randomFit)

434 pars = np.array(pars)

435 meanPfit = np.mean(pars, 0)

436

437 # confidence interval for parameter estimates

438 errPfit = confidenceSigma*np.std(pars, 0)

439 pFitBootstrap = meanPfit

440 pErrBootstrap = errPfit

441

442 reducedChiSq = calculateWeightedReducedChi2(dataY, function(pFitBootstrap, dataX), weightsY, len(dataY),

443 len(initialParams))

444 return pFitBootstrap, pErrBootstrap, reducedChiSq

445

446

447def funcPolynomial(pars, x):

448 """Polynomial function definition

449 Parameters

450 ----------

451 params : `list`

452 Polynomial coefficients. Its length determines the polynomial order.

453

454 x : `numpy.array`, (N,)

455 Abscisa array.

456

457 Returns

458 -------

459 y : `numpy.array`, (N,)

460 Ordinate array after evaluating polynomial of order

461 len(pars)-1 at `x`.

462 """

463 return poly.polyval(x, [*pars])

464

465

466def funcAstier(pars, x):

467 """Single brighter-fatter parameter model for PTC; Equation 16 of

468 Astier+19.

469

470 Model:

471

472 C_{00}(mu) = frac{1}{2 a_{00} g**2} * [exp(2 a_{00} mu g ) - 1]

473 + n_{00} / g**2

474

475 Parameters

476 ----------

477 params : `list`

478 Parameters of the model: a00 (brightter-fatter), gain (e/ADU),

479 and noise (e^2).

480 x : `numpy.array`, (N,)

481 Signal mu (ADU).

482

483 Returns

484 -------

485 y : `numpy.array`, (N,)

486 C_00 (variance) in ADU^2.

487 """

488 a00, gain, noiseSquared = pars

489 return 0.5/(a00*gain*gain)*(np.exp(2*a00*x*gain)-1) + noiseSquared/(gain*gain) # C_00

490

491

492def ptcRolloffModel(params, mu):

493 """Piece-wise exponential saturation roll-off model of the PTC.

494

495 Parameters

496 ----------

497 params : `list`

498 Parameters of the model: muTurnoff (adu), tau (rolloff sharpness).

499 mu : `numpy.array`, (N,)

500 Signal mu (ADU).

501

502 Returns

503 -------

504 y : `numpy.array`, (N,)

505 Difference in variance in ADU^2.

506 """

507 muTurnoff, tau = params

508 return np.where(mu < muTurnoff, 0, np.exp(-(mu - muTurnoff) / tau) - 1)

509

510

511def funcAstierWithRolloff(pars, x):

512 """Single brighter-fatter parameter model for PTC; Equation 16 of

513 Astier+19 with an piece-wise exponential model for the PTC roll-off

514 of the PTC caused by saturation.

515

516 The nominal turnoff is calculated beforehand, and we extend the PTC

517 fit to include signal values up to 5% above the nominally computed

518 turnoff.

519

520 Model:

521

522 C_{00}(mu) = funcAstier - np.where(

523 x < muTurnoff,

524 0,

525 np.exp(-(x - muTurnoff) / tau) - 1,

526 )

527

528 Parameters

529 ----------

530 params : `list`

531 Parameters of the model: a00 (brightter-fatter), gain (e/ADU),

532 and noise (e^2), muTurnoff (adu), tau (rolloff sharpness).

533 x : `numpy.array`, (N,)

534 Signal mu (ADU).

535

536 Returns

537 -------

538 y : `numpy.array`, (N,)

539 C_00 (variance) in ADU^2.

540 """

541 # Initial computation of Astier+19 (Eqn 19).

542 originalModelPars = pars[:-2]

543 ptcRolloffModelPars = pars[-2:]

544 model = funcAstier(originalModelPars, x)

545

546 return model - ptcRolloffModel(ptcRolloffModelPars, x) # C_00

547

548

549def arrangeFlatsByExpTime(exposureList, exposureIdList, log=None):

550 """Arrange exposures by exposure time.

551

552 Parameters

553 ----------

554 exposureList : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]

555 Input list of exposure references.

556 exposureIdList : `list` [`int`]

557 List of exposure ids as obtained by dataId[`exposure`].

558 log : `lsst.utils.logging.LsstLogAdapter`, optional

559 Log object.

560

561 Returns

562 ------

563 flatsAtExpTime : `dict` [`float`,

564 `list`[(`lsst.pipe.base.connections.DeferredDatasetRef`,

565 `int`)]]

566 Dictionary that groups references to flat-field exposures

567 (and their IDs) that have the same exposure time (seconds).

568 """

569 flatsAtExpTime = {}

570 assert len(exposureList) == len(exposureIdList), "Different lengths for exp. list and exp. ID lists"

571 for expRef, expId in zip(exposureList, exposureIdList):

572 expTime = expRef.get(component='visitInfo').exposureTime

573 if not np.isfinite(expTime) and log is not None:

574 log.warning("Exposure %d has non-finite exposure time.", expId)

575 listAtExpTime = flatsAtExpTime.setdefault(expTime, [])

576 listAtExpTime.append((expRef, expId))

577

578 return flatsAtExpTime

579

580

581def arrangeFlatsByExpFlux(exposureList, exposureIdList, fluxKeyword, log=None):

582 """Arrange exposures by exposure flux.

583

584 Parameters

585 ----------

586 exposureList : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]

587 Input list of exposure references.

588 exposureIdList : `list` [`int`]

589 List of exposure ids as obtained by dataId[`exposure`].

590 fluxKeyword : `str`

591 Header keyword that contains the flux per exposure.

592 log : `lsst.utils.logging.LsstLogAdapter`, optional

593 Log object.

594

595 Returns

596 -------

597 flatsAtFlux : `dict` [`float`,

598 `list`[(`lsst.pipe.base.connections.DeferredDatasetRef`,

599 `int`)]]

600 Dictionary that groups references to flat-field exposures

601 (and their IDs) that have the same flux.

602 """

603 flatsAtExpFlux = {}

604 assert len(exposureList) == len(exposureIdList), "Different lengths for exp. list and exp. ID lists"

605 for expRef, expId in zip(exposureList, exposureIdList):

606 # Get flux from header, assuming it is in the metadata.

607 try:

608 expFlux = expRef.get().getMetadata()[fluxKeyword]

609 except KeyError:

610 # If it's missing from the header, continue; it will

611 # be caught and rejected when pairing exposures.

612 expFlux = None

613 if expFlux is None:

614 if log is not None:

615 log.warning("Exposure %d does not have valid header keyword %s.", expId, fluxKeyword)

616 expFlux = np.nan

617 listAtExpFlux = flatsAtExpFlux.setdefault(expFlux, [])

618 listAtExpFlux.append((expRef, expId))

619

620 return flatsAtExpFlux

621

622

623def arrangeFlatsByExpId(exposureList, exposureIdList):

624 """Arrange exposures by exposure ID.

625

626 There is no guarantee that this will properly group exposures,

627 but allows a sequence of flats that have different

628 illumination (despite having the same exposure time) to be

629 processed.

630

631 Parameters

632 ----------

633 exposureList : `list`[`lsst.pipe.base.connections.DeferredDatasetRef`]

634 Input list of exposure references.

635 exposureIdList : `list`[`int`]

636 List of exposure ids as obtained by dataId[`exposure`].

637

638 Returns

639 ------

640 flatsAtExpId : `dict` [`float`,

641 `list`[(`lsst.pipe.base.connections.DeferredDatasetRef`,

642 `int`)]]

643 Dictionary that groups references to flat-field exposures (and their

644 IDs) sequentially by their exposure id.

645

646 Notes

647 -----

648

649 This algorithm sorts the input exposure references by their exposure

650 id, and then assigns each pair of exposure references (exp_j, exp_{j+1})

651 to pair k, such that 2*k = j, where j is the python index of one of the

652 exposure references (starting from zero). By checking for the IndexError

653 while appending, we can ensure that there will only ever be fully

654 populated pairs.

655 """

656 flatsAtExpId = {}

657 assert len(exposureList) == len(exposureIdList), "Different lengths for exp. list and exp. ID lists"

658 # Sort exposures by expIds, which are in the second list `exposureIdList`.

659 sortedExposures = sorted(zip(exposureList, exposureIdList), key=lambda pair: pair[1])

660

661 for jPair, expTuple in enumerate(sortedExposures):

662 if (jPair + 1) % 2:

663 kPair = jPair // 2

664 listAtExpId = flatsAtExpId.setdefault(kPair, [])

665 try:

666 listAtExpId.append(expTuple)

667 listAtExpId.append(sortedExposures[jPair + 1])

668 except IndexError:

669 pass

670

671 return flatsAtExpId

672

673

674def extractCalibDate(calib):

675 """Extract common calibration metadata values that will be written to

676 output header.

677

678 Parameters

679 ----------

680 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`

681 Calibration to pull date information from.

682

683 Returns

684 -------

685 dateString : `str`

686 Calibration creation date string to add to header.

687 """

688 if hasattr(calib, "getMetadata"):

689 if 'CALIB_CREATION_DATE' in calib.getMetadata():

690 return " ".join((calib.getMetadata().get("CALIB_CREATION_DATE", "Unknown"),

691 calib.getMetadata().get("CALIB_CREATION_TIME", "Unknown")))

692 else:

693 return " ".join((calib.getMetadata().get("CALIB_CREATE_DATE", "Unknown"),

694 calib.getMetadata().get("CALIB_CREATE_TIME", "Unknown")))

695 else:

696 return "Unknown Unknown"

697

698

699class CovFastFourierTransform:

700 """A class to compute (via FFT) the nearby pixels correlation function.

701

702 Implements appendix of Astier+19.

703

704 Parameters

705 ----------

706 diff : `numpy.array`

707 Image where to calculate the covariances (e.g., the difference

708 image of two flats).

709 w : `numpy.array`

710 Weight image (mask): it should consist of 1's (good pixel) and

711 0's (bad pixels).

712 fftShape : `tuple`

713 2d-tuple with the shape of the FFT

714 maxRangeCov : `int`

715 Maximum range for the covariances.

716 """

717

718 def __init__(self, diff, w, fftShape, maxRangeCov):

719 # check that the zero padding implied by "fft_shape"

720 # is large enough for the required correlation range

721 assert fftShape[0] > diff.shape[0]+maxRangeCov+1

722 assert fftShape[1] > diff.shape[1]+maxRangeCov+1

723 # for some reason related to numpy.fft.rfftn,

724 # the second dimension should be even, so

725 if fftShape[1]%2 == 1:

726 fftShape = (fftShape[0], fftShape[1]+1)

727 tIm = np.fft.rfft2(diff*w, fftShape)

728 tMask = np.fft.rfft2(w, fftShape)

729 # sum of "squares"

730 self.pCov = np.fft.irfft2(tIm*tIm.conjugate())

731 # sum of values

732 self.pMean = np.fft.irfft2(tIm*tMask.conjugate())

733 # number of w!=0 pixels.

734 self.pCount = np.fft.irfft2(tMask*tMask.conjugate())

735

736 def cov(self, dx, dy):

737 """Covariance for dx,dy averaged with dx,-dy if both non zero.

738

739 Implements appendix of Astier+19.

740

741 Parameters

742 ----------

743 dx : `int`

744 Lag in x

745 dy : `int`

746 Lag in y

747

748 Returns

749 -------

750 0.5*(cov1+cov2): `float`

751 Covariance at (dx, dy) lag

752 npix1+npix2 : `int`

753 Number of pixels used in covariance calculation.

754

755 Raises

756 ------

757 ValueError if number of pixels for a given lag is 0.

758 """

759 # compensate rounding errors

760 nPix1 = int(round(self.pCount[dy, dx]))

761 if nPix1 == 0:

762 raise ValueError(f"Could not compute covariance term {dy}, {dx}, as there are no good pixels.")

763 cov1 = self.pCov[dy, dx]/nPix1-self.pMean[dy, dx]*self.pMean[-dy, -dx]/(nPix1*nPix1)

764 if (dx == 0 or dy == 0):

765 return cov1, nPix1

766 nPix2 = int(round(self.pCount[-dy, dx]))

767 if nPix2 == 0:

768 raise ValueError("Could not compute covariance term {dy}, {dx} as there are no good pixels.")

769 cov2 = self.pCov[-dy, dx]/nPix2-self.pMean[-dy, dx]*self.pMean[dy, -dx]/(nPix2*nPix2)

770 return 0.5*(cov1+cov2), nPix1+nPix2

771

772 def reportCovFastFourierTransform(self, maxRange):

773 """Produce a list of tuples with covariances.

774

775 Implements appendix of Astier+19.

776

777 Parameters

778 ----------

779 maxRange : `int`

780 Maximum range of covariances.

781

782 Returns

783 -------

784 tupleVec : `list`

785 List with covariance tuples.

786 """

787 tupleVec = []

788 # (dy,dx) = (0,0) has to be first

789 for dy in range(maxRange+1):

790 for dx in range(maxRange+1):

791 cov, npix = self.cov(dx, dy)

792 if (dx == 0 and dy == 0):

793 var = cov

794 tupleVec.append((dx, dy, var, cov, npix))

795 return tupleVec

796

797

798def getFitDataFromCovariances(i, j, mu, fullCov, fullCovModel, fullCovSqrtWeights, gain=1.0,

799 divideByMu=False, returnMasked=False):

800 """Get measured signal and covariance, cov model, weigths, and mask at

801 covariance lag (i, j).

802

803 Parameters

804 ----------

805 i : `int`

806 Lag for covariance matrix.

807 j : `int`

808 Lag for covariance matrix.

809 mu : `list`

810 Mean signal values.

811 fullCov : `list` of `numpy.array`

812 Measured covariance matrices at each mean signal level in mu.

813 fullCovSqrtWeights : `list` of `numpy.array`

814 List of square root of measured covariances at each mean

815 signal level in mu.

816 fullCovModel : `list` of `numpy.array`

817 List of modeled covariances at each mean signal level in mu.

818 gain : `float`, optional

819 Gain, in e-/ADU. If other than 1.0 (default), the returned

820 quantities will be in electrons or powers of electrons.

821 divideByMu : `bool`, optional

822 Divide returned covariance, model, and weights by the mean

823 signal mu?

824 returnMasked : `bool`, optional

825 Use mask (based on weights) in returned arrays (mu,

826 covariance, and model)?

827

828 Returns

829 -------

830 mu : `numpy.array`

831 list of signal values at (i, j).

832 covariance : `numpy.array`

833 Covariance at (i, j) at each mean signal mu value (fullCov[:, i, j]).

834 covarianceModel : `numpy.array`

835 Covariance model at (i, j).

836 weights : `numpy.array`

837 Weights at (i, j).

838 maskFromWeights : `numpy.array`, optional

839 Boolean mask of the covariance at (i,j), where the weights

840 differ from 0.

841 """

842 mu = np.array(mu)

843 fullCov = np.array(fullCov)

844 fullCovModel = np.array(fullCovModel)

845 fullCovSqrtWeights = np.array(fullCovSqrtWeights)

846 covariance = fullCov[:, i, j]*(gain**2)

847 covarianceModel = fullCovModel[:, i, j]*(gain**2)

848 weights = fullCovSqrtWeights[:, i, j]/(gain**2)

849

850 maskFromWeights = weights != 0

851 if returnMasked:

852 weights = weights[maskFromWeights]

853 covarianceModel = covarianceModel[maskFromWeights]

854 mu = mu[maskFromWeights]

855 covariance = covariance[maskFromWeights]

856

857 if divideByMu:

858 covariance /= mu

859 covarianceModel /= mu

860 weights *= mu

861 return mu, covariance, covarianceModel, weights, maskFromWeights

862

863

864def symmetrize(inputArray):

865 """ Copy array over 4 quadrants prior to convolution.

866

867 Parameters

868 ----------

869 inputarray : `numpy.array`

870 Input array to symmetrize.

871

872 Returns

873 -------

874 aSym : `numpy.array`

875 Symmetrized array.

876 """

877 targetShape = list(inputArray.shape)

878 r1, r2 = inputArray.shape[-1], inputArray.shape[-2]

879 targetShape[-1] = 2*r1-1

880 targetShape[-2] = 2*r2-1

881 aSym = np.ndarray(tuple(targetShape))

882 aSym[..., r2-1:, r1-1:] = inputArray

883 aSym[..., r2-1:, r1-1::-1] = inputArray

884 aSym[..., r2-1::-1, r1-1::-1] = inputArray

885 aSym[..., r2-1::-1, r1-1:] = inputArray

886

887 return aSym

888

889

890def ddict2dict(d):

891 """Convert nested default dictionaries to regular dictionaries.

892

893 This is needed to prevent yaml persistence issues.

894

895 Parameters

896 ----------

897 d : `defaultdict`

898 A possibly nested set of `defaultdict`.

899

900 Returns

901 -------

902 dict : `dict`

903 A possibly nested set of `dict`.

904 """

905 for k, v in d.items():

906 if isinstance(v, dict):

907 d[k] = ddict2dict(v)

908 return dict(d)

909

910

911class Pol2D:

912 """2D Polynomial Regression.

913

914 Parameters

915 ----------

916 x : numpy.ndarray

917 Input array for the x-coordinate.

918 y : numpy.ndarray

919 Input array for the y-coordinate.

920 z : numpy.ndarray

921 Input array for the dependent variable.

922 order : int

923 Order of the polynomial.

924 w : numpy.ndarray, optional

925 Weight array for weighted regression. Default is None.

926

927 Notes

928 -----

929 Ported from by https://gitlab.in2p3.fr/astier/bfptc P. Astier.

930

931 Example:

932 >>> x = np.array([1, 2, 3])

933 >>> y = np.array([4, 5, 6])

934 >>> z = np.array([7, 8, 9])

935 >>> order = 2

936 >>> poly_reg = Pol2D(x, y, z, order)

937 >>> result = poly_reg.eval(2.5, 5.5)

938 """

939 def __init__(self, x, y, z, order, w=None):

940 """

941 orderx : `int`

942 Effective order in the x-direction.

943 ordery : `int`

944 Effective order in the y-direction.

945 coeff : `numpy.ndarray`

946 Coefficients of the polynomial regression.

947 """

948 self.orderx = min(order, x.shape[0] - 1)

949 self.ordery = min(order, x.shape[1] - 1)

950 G = self.monomials(x.ravel(), y.ravel())

951 if w is None:

952 self.coeff, _, rank, _ = np.linalg.lstsq(G, z.ravel(), rcond=None)

953 else:

954 self.coeff, _, rank, _ = np.linalg.lstsq((w.ravel() * G.T).T, z.ravel() * w.ravel(), rcond=None)

955

956 def monomials(self, x, y):

957 """

958 Generate the monomials matrix for the given x and y.

959

960 Parameters

961 ----------

962 x : numpy.ndarray

963 Input array for the x-coordinate.

964 y : numpy.ndarray

965 Input array for the y-coordinate.

966

967 Returns

968 -------

969 G : numpy.ndarray

970 Monomials matrix.

971 """

972 ncols = (self.orderx + 1) * (self.ordery + 1)

973 G = np.zeros(x.shape + (ncols,))

974 ij = itertools.product(range(self.orderx + 1), range(self.ordery + 1))

975 for k, (i, j) in enumerate(ij):

976 G[..., k] = x**i * y**j

977 return G

978

979 def eval(self, x, y):

980 """

981 Evaluate the polynomial at the given x and y coordinates.

982

983 Parameters

984 ----------

985 x : `float`

986 x-coordinate for evaluation.

987 y : `float`

988 y-coordinate for evaluation.

989

990 Returns

991 -------

992 result : `float`

993 Result of the polynomial evaluation.

994 """

995 G = self.monomials(x, y)

996 return np.dot(G, self.coeff)

997

998

999class AstierSplineLinearityFitter:

1000 """Class to fit the Astier spline linearity model.

1001

1002 This is a spline fit with photodiode data based on a model

1003 from Pierre Astier, referenced in June 2023 from

1004 https://me.lsst.eu/astier/bot/7224D/model_nonlin.py

1005

1006 This model fits a spline with (optional) nuisance parameters

1007 to allow for different linearity coefficients with different

1008 photodiode settings. The minimization is a least-squares

1009 fit with the residual of

1010 Sum[(S(mu_i) + mu_i)/(k_j * D_i - O) - 1]**2, where S(mu_i)

1011 is an Akima Spline function of mu_i, the observed flat-pair

1012 mean; D_j is the photo-diode measurement corresponding to

1013 that flat-pair; and k_j is a constant of proportionality

1014 which is over index j as it is allowed to

1015 be different based on different photodiode settings (e.g.

1016 CCOBCURR); and O is a constant offset to allow for light

1017 leaks (and is only fit if fit_offset=True). In

1018 addition, if config.doSplineFitTemperature is True then

1019 the fit will adjust mu such that

1020 mu = mu_input*(1 + alpha*(temperature_scaled))

1021 and temperature_scaled = T - T_ref. Finally, if

1022 config.doSplineFitTemporal is True then the fit will

1023 further adjust mu such that

1024 mu = mu_input*(1 + beta*(mjd_scaled))

1025 and mjd_scaled = mjd - mjd_ref.

1026

1027 The fit has additional constraints to ensure that the spline

1028 goes through the (0, 0) point, as well as a normalization

1029 condition so that the average of the spline over the full

1030 range is 0. The normalization ensures that the spline only

1031 fits deviations from linearity, rather than the linear

1032 function itself which is degenerate with the gain.

1033

1034 Parameters

1035 ----------

1036 nodes : `np.ndarray` (N,)

1037 Array of spline node locations.

1038 grouping_values : `np.ndarray` (M,)

1039 Array of values to group values for different proportionality

1040 constants (e.g. CCOBCURR).

1041 pd : `np.ndarray` (M,)

1042 Array of photodiode measurements.

1043 mu : `np.ndarray` (M,)

1044 Array of flat mean values.

1045 mask : `np.ndarray` (M,), optional

1046 Input mask (True is good point, False is bad point).

1047 log : `logging.logger`, optional

1048 Logger object to use for logging.

1049 fit_offset : `bool`, optional

1050 Fit a constant offset to allow for light leaks?

1051 fit_weights : `bool`, optional

1052 Fit for the weight parameters?

1053 weight_pars_start : `list` [ `float` ]

1054 Iterable of 2 weight parameters for weighed fit. These will

1055 be used as input if the weight parameters are not fit.

1056 fit_temperature : `bool`, optional

1057 Fit for temperature scaling?

1058 temperature_scaled : `np.ndarray` (M,), optional

1059 Input scaled temperature values (T - T_ref).

1060 max_signal_nearly_linear : `float`, optional

1061 Maximum signal that we are confident the input is nearly

1062 linear. This is used both for regularization, and for

1063 fitting the raw slope. Usually set to the ptc turnoff,

1064 above which we allow the spline to significantly deviate

1065 and do not demand the deviation to average to zero.

1066 fit_temporal : `bool`, optional

1067 Fit for temporal scaling?

1068 mjd_scaled : `np.ndarray` (M,), optional

1069 Input scaled mjd values (mjd - mjd_ref).

1070 max_frac_correction : `float`, optional

1071 Maximum fractional correction.

1072 max_correction : `float`, optional

1073 Maximum correction (unscaled).

1074 """

1075 def __init__(

1076 self,

1077 nodes,

1078 grouping_values,

1079 pd,

1080 mu,

1081 mask=None,

1082 log=None,

1083 fit_offset=True,

1084 fit_weights=False,

1085 weight_pars_start=[1.0, 0.0],

1086 fit_temperature=False,

1087 temperature_scaled=None,

1088 max_signal_nearly_linear=None,

1089 fit_temporal=False,

1090 mjd_scaled=None,

1091 max_frac_correction=0.25,

1092 max_correction=np.inf,

1093 ):

1094 self._pd = pd

1095 self._mu = mu

1096 self._grouping_values = grouping_values

1097 self.log = log if log else logging.getLogger(__name__)

1098 self._fit_offset = fit_offset

1099 self._fit_weights = fit_weights

1100 self._weight_pars_start = weight_pars_start

1101 self._fit_temperature = fit_temperature

1102 self._fit_temporal = fit_temporal

1103 self._max_frac_correction = max_frac_correction

1104 self._max_correction = max_correction

1105

1106 self._nodes = nodes

1107 if nodes[0] != 0.0:

1108 raise ValueError("First node must be 0.0")

1109 if not np.all(np.diff(nodes) > 0):

1110 raise ValueError("Nodes must be sorted with no repeats.")

1111

1112 # Find the group indices.

1113 u_group_values = np.unique(self._grouping_values)

1114 self.ngroup = len(u_group_values)

1115

1116 self.group_indices = []

1117 for i in range(self.ngroup):

1118 self.group_indices.append(np.where(self._grouping_values == u_group_values[i])[0])

1119

1120 # Weight values. Outliers will be set to 0.

1121 if mask is None:

1122 _mask = np.ones(len(mu), dtype=np.bool_)

1123 else:

1124 _mask = mask

1125 self._w = self.compute_weights(self._weight_pars_start, self._mu, _mask)

1126

1127 if temperature_scaled is None:

1128 temperature_scaled = np.zeros(len(self._mu))

1129 else:

1130 if len(np.atleast_1d(temperature_scaled)) != len(self._mu):

1131 raise ValueError("temperature_scaled must be the same length as input mu.")

1132 self._temperature_scaled = temperature_scaled

1133

1134 if mjd_scaled is None:

1135 mjd_scaled = np.zeros(len(self._mu))

1136 else:

1137 if len(np.atleast_1d(mjd_scaled)) != len(self._mu):

1138 raise ValueError("mjd_scaled must be the same length as input mu.")

1139 self._mjd_scaled = mjd_scaled

1140

1141 # Values to regularize spline fit.

1142 if max_signal_nearly_linear is None:

1143 max_signal_nearly_linear = self._mu[self.mask].max()

1144 self._max_signal_nearly_linear = max_signal_nearly_linear

1145 self._x_regularize = np.linspace(0.0, self._max_signal_nearly_linear, 100)

1146

1147 # Set up the indices for the fit parameters.

1148 self.par_indices = {

1149 "values": np.arange(len(self._nodes)),

1150 "groups": len(self._nodes) + np.arange(self.ngroup),

1151 "offset": np.zeros(0, dtype=np.int64),

1152 "weight_pars": np.zeros(0, dtype=np.int64),

1153 "temperature_coeff": np.zeros(0, dtype=np.int64),

1154 "temporal_coeff": np.zeros(0, dtype=np.int64),

1155 }

1156 if self._fit_offset:

1157 self.par_indices["offset"] = np.arange(1) + (

1158 len(self.par_indices["values"])

1159 + len(self.par_indices["groups"])

1160 )

1161 if self._fit_weights:

1162 self.par_indices["weight_pars"] = np.arange(2) + (

1163 len(self.par_indices["values"])

1164 + len(self.par_indices["groups"])

1165 + len(self.par_indices["offset"])

1166 )

1167 if self._fit_temperature:

1168 self.par_indices["temperature_coeff"] = np.arange(1) + (

1169 len(self.par_indices["values"])

1170 + len(self.par_indices["groups"])

1171 + len(self.par_indices["offset"])

1172 + len(self.par_indices["weight_pars"])

1173 )

1174 if self._fit_temporal:

1175 self.par_indices["temporal_coeff"] = np.arange(1) + (

1176 len(self.par_indices["values"])

1177 + len(self.par_indices["groups"])

1178 + len(self.par_indices["offset"])

1179 + len(self.par_indices["weight_pars"])

1180 + len(self.par_indices["temperature_coeff"])

1181 )

1182

1183 self._npt = (

1184 len(self.par_indices["values"])

1185 + len(self.par_indices["groups"])

1186 + len(self.par_indices["offset"])

1187 + len(self.par_indices["weight_pars"])

1188 + len(self.par_indices["temperature_coeff"])

1189 + len(self.par_indices["temporal_coeff"])

1190 )

1191

1192 @staticmethod

1193 def compute_weights(weight_pars, mu, mask):

1194 w = 1./np.sqrt(weight_pars[0]**2. + weight_pars[1]**2./mu)

1195 w[~mask] = 0.0

1196

1197 return w

1198

1199 def estimate_p0(self, use_all_for_normalization=False):

1200 """Estimate initial fit parameters.

1201

1202 Parameters

1203 ----------

1204 use_all_for_normalization : `bool`, optional

1205 Use all the points (not just below turnoff) for normalization;

1206 this is for compatibility with the old linearizer fits.

1207

1208 Returns

1209 -------

1210 p0 : `np.ndarray`

1211 Parameter array, with spline values (one for each node) followed

1212 by proportionality constants (one for each group); one extra

1213 for the offset O (if fit_offset was set to True); two extra

1214 for the weights (if fit_weights was set to True); one

1215 extra for the temperature coefficient (if fit_temperature was

1216 set to True); and one extra for the temporal coefficient (if

1217 fit_temporal was set to True).

1218 """

1219 p0 = np.zeros(self._npt)

1220

1221 # Do a simple linear fit for each group.

1222 for i, indices in enumerate(self.group_indices):

1223 mask = self.mask[indices]

1224 if mask.sum() == 0:

1225 # There are no points in this group; we can ignore it.

1226 p0[self.par_indices["groups"][i]] = 1.0

1227 continue

1228 mu = self._mu[indices][mask]

1229 pd = self._pd[indices][mask]

1230 to_fit = (mu < self._max_signal_nearly_linear)

1231 # If we have no points now (low PTC turnoff) then

1232 # we just use the first 25% of the points in the group.

1233 if to_fit.sum() == 0:

1234 to_fit = (mu < np.nanpercentile(mu, 25.0))

1235 linfit = np.polyfit(pd[to_fit], mu[to_fit], 1)

1236 p0[self.par_indices["groups"][i]] = linfit[0]

1237

1238 # This re-fit will be skipped for compatibility.

1239 # if use_all_for_normalization is set.

1240 if not use_all_for_normalization:

1241 params, cov_params, _, _, _ = leastsq(

1242 self._group_minfunc,

1243 p0[self.par_indices["groups"]],

1244 full_output=True,

1245 ftol=1e-5,

1246 maxfev=12000,

1247 )

1248

1249 p0[self.par_indices["groups"]] = params

1250

1251 # Look at the residuals...

1252 ratio_model = self.compute_ratio_model(

1253 self._nodes,

1254 self.group_indices,

1255 self.par_indices,

1256 p0,

1257 self._pd,

1258 self._mu,

1259 self._temperature_scaled,

1260 self._mjd_scaled,

1261 )

1262 # ...and adjust the linear parameters accordingly.

1263 if use_all_for_normalization:

1264 ok = self.mask

1265 else:

1266 ok = (self.mask & (self._mu < self._max_signal_nearly_linear))

1267 p0[self.par_indices["groups"]] *= np.median(ratio_model[ok])

1268

1269 # Re-compute the residuals.

1270 ratio_model2 = self.compute_ratio_model(

1271 self._nodes,

1272 self.group_indices,

1273 self.par_indices,

1274 p0,

1275 self._pd,

1276 self._mu,

1277 self._temperature_scaled,

1278 self._mjd_scaled,

1279 )

1280

1281 # And compute a first guess of the spline nodes.

1282 bins = np.clip(np.searchsorted(self._nodes, self._mu[self.mask]), 0, len(self._nodes) - 1)

1283 tot_arr = np.zeros(len(self._nodes))

1284 n_arr = np.zeros(len(self._nodes), dtype=int)

1285 np.add.at(tot_arr, bins, ratio_model2[self.mask])

1286 np.add.at(n_arr, bins, 1)

1287

1288 ratio = np.ones(len(self._nodes))

1289 ratio[n_arr > 0] = tot_arr[n_arr > 0]/n_arr[n_arr > 0]

1290 ratio[0] = 1.0

1291 p0[self.par_indices["values"]] = (ratio - 1) * self._nodes

1292

1293 if self._fit_offset:

1294 p0[self.par_indices["offset"]] = 0.0

1295

1296 if self._fit_weights:

1297 p0[self.par_indices["weight_pars"]] = self._weight_pars_start

1298

1299 return p0

1300

1301 @staticmethod

1302 def compute_ratio_model(

1303 nodes,

1304 group_indices,

1305 par_indices,

1306 pars,

1307 pd,

1308 mu,

1309 temperature_scaled,

1310 mjd_scaled,

1311 return_spline=False,

1312 ):

1313 """Compute the ratio model values.

1314

1315 Parameters

1316 ----------

1317 nodes : `np.ndarray` (M,)

1318 Array of node positions.

1319 group_indices : `list` [`np.ndarray`]

1320 List of group indices, one array for each group.

1321 par_indices : `dict`

1322 Dictionary showing which indices in the pars belong to

1323 each set of fit values.

1324 pars : `np.ndarray`

1325 Parameter array, with spline values (one for each node) followed

1326 by proportionality constants (one for each group); one extra

1327 for the offset O (if fit_offset was set to True); two extra

1328 for the weights (if fit_weights was set to True); and one

1329 extra for the temperature coefficient (if fit_temperature was

1330 set to True).

1331 pd : `np.ndarray` (N,)

1332 Array of photodiode measurements.

1333 mu : `np.ndarray` (N,)

1334 Array of flat means.

1335 temperature_scaled : `np.ndarray` (N,)

1336 Array of scaled temperature values.

1337 mjd_scaled : `np.ndarray` (N,)

1338 Array of scaled mjd values.

1339 return_spline : `bool`, optional

1340 Return the spline interpolation as well as the model ratios?

1341

1342 Returns

1343 -------

1344 ratio_models : `np.ndarray` (N,)

1345 Model ratio, (mu_i - S(mu_i) - O)/(k_j * D_i)

1346 spl : `scipy.interpolate.Akima1DInterpolator`

1347 Spline interpolator (returned if return_spline=True).

1348 """

1349 spl = Akima1DInterpolator(nodes, pars[par_indices["values"]], method="akima")

1350

1351 # Check if we want to do just the left or both with temp scale.

1352 if len(par_indices["temperature_coeff"]) == 1:

1353 mu_corr = mu*(1. + pars[par_indices["temperature_coeff"]]*temperature_scaled)

1354 else:

1355 mu_corr = mu

1356

1357 if len(par_indices["temporal_coeff"]) == 1:

1358 mu_corr = mu_corr*(1. + pars[par_indices["temporal_coeff"]]*mjd_scaled)

1359

1360 numerator = mu_corr - spl(np.clip(mu_corr, nodes[0], nodes[-1]))

1361 if len(par_indices["offset"]) == 1:

1362 numerator -= pars[par_indices["offset"][0]]

1363 denominator = pd.copy()

1364 ngroup = len(group_indices)

1365 kj = pars[par_indices["groups"]]

1366 for j in range(ngroup):

1367 denominator[group_indices[j]] *= kj[j]

1368

1369 if return_spline:

1370 return numerator / denominator, spl

1371 else:

1372 return numerator / denominator

1373

1374 def fit(self, p0, min_iter=3, max_iter=20, max_rejection_per_iteration=5, n_sigma_clip=5.0):

1375 """

1376 Perform iterative fit for linear + spline model with offsets.

1377

1378 Parameters

1379 ----------

1380 p0 : `np.ndarray`

1381 Initial parameter array, with spline values (one for each node)

1382 followed by proportionality constants (one for each group); one

1383 extra for the offset O (if fit_offset was set to True); two extra

1384 for the weights (if fit_weights was set to True); one

1385 extra for the temperature coefficient (if fit_temperature was

1386 set to True); and one extra for the temporal coefficient (if

1387 fit_temporal was set to True).

1388 min_iter : `int`, optional

1389 Minimum number of fit iterations.

1390 max_iter : `int`, optional

1391 Maximum number of fit iterations.

1392 max_rejection_per_iteration : `int`, optional

1393 Maximum number of points to reject per iteration.

1394 n_sigma_clip : `float`, optional

1395 Number of sigma to do clipping in each iteration.

1396 """

1397 init_params = p0

1398 for k in range(max_iter):

1399 params, cov_params, _, msg, ierr = leastsq(

1400 self,

1401 init_params,

1402 full_output=True,

1403 ftol=1e-5,

1404 maxfev=12000,

1405 )

1406 init_params = params.copy()

1407

1408 # We need to cut off the constraints at the end (there are more

1409 # residuals than data points.)

1410 res = self(params)[: len(self._w)]

1411 std_res = median_abs_deviation(res[self.good_points], scale="normal")

1412 sample = len(self.good_points)

1413

1414 # We don't want to reject too many outliers at once.

1415 if sample > max_rejection_per_iteration:

1416 sres = np.sort(np.abs(res))

1417 cut = max(sres[-max_rejection_per_iteration], std_res*n_sigma_clip)

1418 else:

1419 cut = std_res*n_sigma_clip

1420

1421 outliers = np.abs(res) > cut

1422 self._w[outliers] = 0

1423 if outliers.sum() != 0:

1424 self.log.info(

1425 "After iteration %d there are %d outliers (of %d).",

1426 k,

1427 outliers.sum(),

1428 sample,

1429 )

1430 elif k >= min_iter:

1431 self.log.info("After iteration %d there are no more outliers.", k)

1432 break

1433

1434 return params

1435

1436 @property

1437 def mask(self):

1438 return (self._w > 0)

1439

1440 @property

1441 def good_points(self):

1442 return self.mask.nonzero()[0]

1443

1444 def compute_chisq_dof(self, pars):

1445 """Compute the chi-squared per degree of freedom for a set of pars.

1446

1447 Parameters

1448 ----------

1449 pars : `np.ndarray`

1450 Parameter array.

1451

1452 Returns

1453 -------

1454 chisq_dof : `float`

1455 Chi-squared per degree of freedom.

1456 """

1457 resids = self(pars)[0: len(self.mask)]

1458 chisq = np.sum(resids[self.mask]**2.)

1459 dof = self.mask.sum() - self.ngroup

1460 if self._fit_temporal:

1461 dof -= 1

1462 if self._fit_temperature:

1463 dof -= 1

1464 if self._fit_offset:

1465 dof -= 1

1466 if self._fit_weights:

1467 dof -= 2

1468

1469 return chisq/dof

1470

1471 def __call__(self, pars):

1472

1473 ratio_model, spl = self.compute_ratio_model(

1474 self._nodes,

1475 self.group_indices,

1476 self.par_indices,

1477 pars,

1478 self._pd,

1479 self._mu,

1480 self._temperature_scaled,

1481 self._mjd_scaled,

1482 return_spline=True,

1483 )

1484

1485 _mask = self.mask

1486 # Update the weights if we are fitting them.

1487 if self._fit_weights:

1488 self._w = self.compute_weights(pars[self.par_indices["weight_pars"]], self._mu, _mask)

1489 resid = self._w*(ratio_model - 1.0)

1490

1491 # Ensure masked points have 0 residual.

1492 resid[~_mask] = 0.0

1493

1494 constraint = [1e3 * np.mean(spl(np.clip(self._x_regularize, self._nodes[0], self._nodes[-1])))]

1495 # 0 should transform to 0

1496 constraint.append(spl(0)*1e10)

1497 # Use a Jeffreys prior on the weight if we are fitting it.

1498 if self._fit_weights:

1499 # This factor ensures that log(fact * w) is negative.

1500 fact = 1e-3 / self._w.max()

1501 # We only add non-zero weights to the constraint array.

1502 log_w = np.sqrt(-2.*np.log(fact*self._w[self._w > 0]))

1503 constraint = np.hstack([constraint, log_w])

1504

1505 # Don't let it get to >5% correction.

1506 values = pars[self.par_indices["values"]]

1507 if np.abs(values[-1])/self._nodes[-1] > self._max_frac_correction \

1508 or np.abs(values[-1]) > self._max_correction:

1509 extra_constraint = 1e10

1510 else:

1511 extra_constraint = 0

1512

1513 return np.hstack([resid, constraint, extra_constraint])

1514

1515 def _group_minfunc(self, group_pars):

1516 """Minimization function for initial group parameters.

1517

1518 Parameters

1519 ----------

1520 group_pars : `np.ndarray`

1521 Array of slope parameters.

1522 """

1523 full_pars = np.zeros(self._npt)

1524 full_pars[self.par_indices["groups"]] = group_pars

1525

1526 ratio_model = self.compute_ratio_model(

1527 self._nodes,

1528 self.group_indices,

1529 self.par_indices,

1530 full_pars,

1531 self._pd,

1532 self._mu,

1533 self._temperature_scaled,

1534 self._mjd_scaled,

1535 return_spline=False,

1536 )

1537 wt = self._w.copy()

1538 wt[self._mu > self._max_signal_nearly_linear] = 0.0

1539 return wt*(ratio_model - 1.0)

1540

1541

1542def getReadNoise(exposure, ampName, taskMetadata=None, log=None):

1543 """Gets readout noise for an amp from ISR metadata.

1544

1545 If possible, this attempts to get the now-standard headers

1546 added to the exposure itself. If not found there, the ISR

1547 TaskMetadata is searched. If neither of these has the value,

1548 warn and set the read noise to NaN.

1549

1550 Parameters

1551 ----------

1552 exposure : `lsst.afw.image.Exposure`

1553 Exposure to check for read noise first.

1554 ampName : `str`

1555 Amplifier name.

1556 taskMetadata : `lsst.pipe.base.TaskMetadata`, optional

1557 List of exposures metadata from ISR for this exposure.

1558 log : `logging.logger`, optional

1559 Log for messages.

1560

1561 Returns

1562 -------

1563 readNoise : `float`

1564 The read noise for this set of exposure/amplifier.

1565 """

1566 exposureMetadata = exposure.getMetadata()

1567

1568 # Try from the exposure first.

1569 expectedKey = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"

1570 if expectedKey in exposureMetadata:

1571 return exposureMetadata[expectedKey]

1572

1573 # If not, try getting it from the task metadata.

1574 if taskMetadata:

1575 expectedKey = f"RESIDUAL STDEV {ampName}"

1576 if "isr" in taskMetadata:

1577 if expectedKey in taskMetadata["isr"]:

1578 return taskMetadata["isr"][expectedKey]

1579

1580 log = log if log else logging.getLogger(__name__)

1581 log.warning("Median readout noise from ISR metadata for amp %s "

1582 "could not be found.", ampName)

1583 return np.nan

1584

1585

1586def ampOffsetGainRatioFixup(ptc, minAdu, maxAdu, log=None):

1587 """Apply gain ratio fixup using amp offsets.

1588

1589 Parameters

1590 ----------

1591 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`

1592 Input PTC. Will be modified in place.

1593 minAdu : `float`

1594 Minimum number of ADU in mean of amplifier to use for

1595 computing gain ratio fixup.

1596 maxAdu : `float`

1597 Maximum number of ADU in mean of amplifier to use for

1598 computing gain ratio fixup.

1599 log : `lsst.utils.logging.LsstLogAdapter`, optional

1600 Log object.

1601 """

1602 # We need to find the reference amplifier.

1603 # Find an amp near the middle to use as a pivot.

1604 if log is None:

1605 log = logging.getLogger(__name__)

1606

1607 # We first check for which amps have gain measurements

1608 # (fully bad amps are filled with NaN for gain.)

1609 gainArray = np.zeros(len(ptc.ampNames))

1610 for i, ampName in enumerate(ptc.ampNames):

1611 gainArray[i] = ptc.gain[ampName]

1612 good, = np.where(np.isfinite(gainArray))

1613

1614 if len(good) > 1:

1615 # This only works with more than 1 good amp.

1616

1617 # We sort the gains and take the one that is closest

1618 # to the median to use as the reference amplifier for

1619 # gain ratios.

1620 st = np.argsort(gainArray[good])

1621 midAmp = good[st[int(0.5*len(good))]]

1622 midAmpName = ptc.ampNames[midAmp]

1623

1624 log.info("Using amplifier %s as the reference for doLinearityGainRatioFixup.", midAmpName)

1625

1626 # First pass, we need to compute the corrections.

1627 corrections = {}

1628 for ampName in ptc.ampNames:

1629 if not np.isfinite(ptc.gain[ampName]) or ampName == midAmpName:

1630 continue

1631

1632 ratioPtc = ptc.gain[ampName] / ptc.gain[midAmpName]

1633

1634 deltas = ptc.ampOffsets[ampName] - ptc.ampOffsets[midAmpName]

1635 use = (

1636 (ptc.expIdMask[ampName])

1637 & (np.isfinite(deltas))

1638 & (ptc.finalMeans[ampName] >= minAdu)

1639 & (ptc.finalMeans[ampName] <= maxAdu)

1640 & (np.isfinite(ptc.finalMeans[midAmpName]))

1641 & (ptc.expIdMask[midAmpName])

1642 )

1643 if use.sum() < 3:

1644 log.warning("Not enough good amp offset measurements to fix up amp %s "

1645 "gains from amp ratios.", ampName)

1646 continue

1647

1648 ratios = 1. / (deltas / ptc.finalMeans[midAmpName] + 1.0)

1649 ratio = np.median(ratios[use])

1650 corrections[ampName] = ratio / ratioPtc

1651

1652 # For the final correction, we need to make sure that the

1653 # reference amplifier is included. By definition, it has a

1654 # correction factor of 1.0 before any final fix.

1655 corrections[midAmpName] = 1.0

1656

1657 # Adjust the median correction to be 1.0 so we do not

1658 # change the gain of the detector on average.

1659 # This is needed in case the reference amplifier is

1660 # skewed in terms of offsets even though it has the median

1661 # gain.

1662 medCorrection = np.median([corrections[key] for key in corrections])

1663

1664 for ampName in ptc.ampNames:

1665 if ampName not in corrections:

1666 continue

1667

1668 correction = corrections[ampName] / medCorrection

1669 newGain = ptc.gain[ampName] * correction

1670 log.info(

1671 "Adjusting gain from amplifier %s by factor of %.5f (from %.5f to %.5f)",

1672 ampName,

1673 correction,

1674 ptc.gain[ampName],

1675 newGain,

1676 )

1677 # Copying the value should not be necessary, but we record

1678 # it just in case.

1679 ptc.gainUnadjusted[ampName] = ptc.gain[ampName]

1680 ptc.gain[ampName] = newGain

1681 else:

1682 log.warning("Cannot apply ampOffsetGainRatioFixup with fewer than 2 good amplifiers.")

1683

1684

1685class FlatGradientFitter:

1686 """Class to fit various large-scale flat-field gradients.

1687

1688 This fitter will take arrays of x/y/value and (by default) fit a radial

1689 gradient, using a spline function at the specified nodes. The fitter

1690 will also fit out a nuisance parameter for the ratio of the ITL/E2V

1691 throughput (though in general this could work for any such focal plane).

1692 The focal plane origin is set by fp_centroid_x and fp_centroid_y which

1693 is the center of the radial gradient, though it may be modified if

1694 fit_centroid is True. In addition, the fitter may fit a linear gradient

1695 in x/y if fit_gradient is True. The "pivot" of the gradient is at

1696 fp_centroid_x/fp_centroid_y.

1697

1698 Parameters

1699 ----------

1700 nodes : `np.ndarray`

1701 Array of spline nodes for radial spline.

1702 x : `np.ndarray`

1703 Array of x values for points to fit (mm).

1704 y : `np.ndarray`

1705 Array of y values for points to fit (mm).

1706 value : `np.ndarray`

1707 Array of values describing the flat field to fit the gradient.

1708 itl_indices : `np.ndarray`

1709 Array of indices corresponding to the ITL detectors.

1710 constrain_zero : `bool`, optional

1711 Constrain the outermost radial spline value to 0?

1712 fit_centroid : `bool`, optional

1713 Fit an additional centroid offset?

1714 fit_gradient : `bool`, optional

1715 Fit an additional plane gradient?

1716 fp_centroid_x : `float`, optional

1717 Focal plane centroid x (mm).

1718 fp_centroid_y : `float`, optional

1719 Focal plane centroid y (mm).

1720 """

1721 def __init__(

1722 self,

1723 nodes,

1724 x,

1725 y,

1726 value,

1727 itl_indices,

1728 constrain_zero=True,

1729 fit_centroid=False,

1730 fit_gradient=False,

1731 fp_centroid_x=0.0,

1732 fp_centroid_y=0.0

1733 ):

1734 self._nodes = nodes

1735 self._x = x

1736 self._y = y

1737 self._value = value

1738 self._itl_indices = itl_indices

1739 self._fp_centroid_x = fp_centroid_x

1740 self._fp_centroid_y = fp_centroid_y

1741

1742 self._constrain_zero = constrain_zero

1743

1744 self._fit_centroid = fit_centroid

1745 self._fit_gradient = fit_gradient

1746

1747 self.indices = {"spline": np.arange(len(nodes))}

1748 npar = len(nodes)

1749

1750 self._fit_itl_ratio = False

1751 if len(itl_indices) > 0:

1752 self._fit_itl_ratio = True

1753 self.indices["itl_ratio"] = npar

1754 npar += 1

1755

1756 radius = np.sqrt((self._x - self._fp_centroid_x)**2. + (self._y - self._fp_centroid_y)**2.)

1757

1758 if fit_centroid:

1759 self.indices["centroid_delta"] = np.arange(2) + npar

1760 npar += 2

1761 else:

1762 self._radius = radius

1763

1764 if fit_gradient:

1765 self.indices["gradient"] = np.arange(2) + npar

1766 npar += 2

1767

1768 self._npar = npar

1769

1770 self._bounds = [(None, None)]*npar

1771

1772 if self._constrain_zero:

1773 self._bounds[self.indices["spline"][-1]] = (0.0, 0.0)

1774

1775 def compute_p0(self, itl_ratio=None):

1776 """Compute initial guess for fit parameters.

1777

1778 Returns

1779 -------

1780 pars : `np.ndarray`

1781 Array of first guess fit parameters.

1782 """

1783 pars = np.zeros(self._npar)

1784 value = self._value.copy()

1785

1786 if itl_ratio is not None and self._fit_itl_ratio:

1787 value[self._itl_indices] /= itl_ratio

1788

1789 # Initial spline values

1790 radius = np.sqrt((self._x - self._fp_centroid_x)**2. + (self._y - self._fp_centroid_y)**2.)

1791 for i, index in enumerate(self.indices["spline"]):

1792 if i == 0:

1793 low = self._nodes[i]

1794 else:

1795 low = (self._nodes[i - 1] + self._nodes[i])/2.

1796 if i == (len(self._nodes) - 1):

1797 high = self._nodes[i]

1798 else:

1799 high = (self._nodes[i] + self._nodes[i + 1])/2.

1800 u = ((radius > low) & (radius < high))

1801 if u.sum() == 0:

1802 pars[index] = 0.0

1803 else:

1804 pars[index] = np.median(value[u])

1805

1806 if self._constrain_zero:

1807 pars[self.indices["spline"][-1]] = 0.0

1808

1809 spl = Akima1DInterpolator(self._nodes, pars[self.indices["spline"]], method="akima")

1810 model = spl(radius)

1811 resid_ratio = value / model

1812

1813 if self._fit_itl_ratio:

1814 if itl_ratio is not None:

1815 pars[self.indices["itl_ratio"]] = itl_ratio

1816 else:

1817 e2v_indices = np.delete(np.arange(len(self._value)), self._itl_indices)

1818

1819 itl_inner = radius[self._itl_indices] < 0.8*np.max(radius)

1820 e2v_inner = radius[e2v_indices] < 0.8*np.max(radius)

1821

1822 itl_median = np.nanmedian(resid_ratio[self._itl_indices][itl_inner])

1823 e2v_median = np.nanmedian(resid_ratio[e2v_indices][e2v_inner])

1824

1825 pars[self.indices["itl_ratio"]] = itl_median / e2v_median

1826

1827 if self._fit_centroid:

1828 pars[self.indices["centroid_delta"]] = [0.0, 0.0]

1829

1830 if self._fit_gradient:

1831 resid_ratio = self._value / model

1832

1833 ok = (np.isfinite(resid_ratio) & (model > 0.5))

1834

1835 fit = np.polyfit(self._x[ok], resid_ratio[ok], 1)

1836 pars[self.indices["gradient"][0]] = fit[0]

1837

1838 fit = np.polyfit(self._y[ok], resid_ratio[ok], 1)

1839 pars[self.indices["gradient"][1]] = fit[0]

1840

1841 return pars

1842

1843 def fit(self, p0, freeze_itl_ratio=False, fit_eps=1e-8, fit_gtol=1e-10):

1844 """Do a non-linear minimization to fit the parameters.

1845

1846 Parameters

1847 ----------

1848 p0 : `np.ndarray`

1849 Array of initial parameter estimates.

1850 freeze_itl_ratio : `bool`, optional

1851 Freeze the ITL ratio in the fit?

1852 fit_eps : `float`, optional

1853 Value of ``eps`` to send to the scipy minimizer.

1854 fit_gtol : `float`, optional

1855 Value of ``gtol`` to send to the scipy minimizer.

1856

1857 Returns

1858 -------

1859 pars : `np.ndarray`

1860 Array of parameters. Use ``fitter.indices`` for the

1861 dictionary to map parameters to subsets.

1862 """

1863 bounds = self._bounds

1864 if freeze_itl_ratio and self._fit_itl_ratio:

1865 ind = self.indices["itl_ratio"]

1866 par = p0[ind]

1867 bounds[ind] = (par, par)

1868

1869 res = minimize(

1870 self,

1871 p0,

1872 method="L-BFGS-B",

1873 jac=False,

1874 bounds=bounds,

1875 options={

1876 "maxfun": 10000,

1877 "maxiter": 10000,

1878 "maxcor": 20,

1879 "eps": fit_eps,

1880 "gtol": fit_gtol,

1881 },

1882 callback=None,

1883 )

1884 pars = res.x

1885

1886 return pars

1887

1888 def compute_model(self, pars):

1889 """Compute the model given a set of parameters.

1890

1891 Parameters

1892 ----------

1893 pars : `np.ndarray`

1894 Parameter array to compute model.

1895

1896 Returns

1897 -------

1898 model_array : `np.ndarray`

1899 Array of model parameters at the input x/y.

1900 """

1901 spl = Akima1DInterpolator(self._nodes, pars[self.indices["spline"]], method="akima")

1902 if self._fit_centroid:

1903 centroid_delta = pars[self.indices["centroid_delta"]]

1904 centroid_x = self._fp_centroid_x + centroid_delta[0]

1905 centroid_y = self._fp_centroid_y + centroid_delta[1]

1906 radius = np.sqrt((self._x - centroid_x)**2. + (self._y - centroid_y)**2.)

1907 else:

1908 radius = self._radius

1909

1910 model = spl(np.clip(radius, self._nodes[0], self._nodes[-1]))

1911 if self._fit_itl_ratio:

1912 model[self._itl_indices] *= pars[self.indices["itl_ratio"]]

1913

1914 if self._fit_gradient:

1915 a, b = pars[self.indices["gradient"]]

1916 gradient = 1 + a*(self._x - self._fp_centroid_x) + b*(self._y - self._fp_centroid_y)

1917 model /= gradient

1918

1919 return model

1920

1921 def __call__(self, pars):

1922 """Compute the cost function for a set of parameters.

1923

1924 Parameters

1925 ----------

1926 pars : `np.ndarray`

1927 Parameter array to compute model.

1928

1929 Returns

1930 -------

1931 cost : `float`

1932 Cost value computed from the absolute deviation.

1933 """

1934

1935 model = self.compute_model(pars)

1936

1937 absdev = np.abs(self._value - model)

1938 t = np.sum(absdev.astype(np.float64))

1939

1940 return t

1941

1942

1943class FlatGainRatioFitter:

1944 """Class to fit amplifier gain ratios, removing a simple gradient.

1945

1946 This fitter will take arrays of x/y/amp_num/value and fit amplifier

1947 gain ratios, using one amplifier as the fixed point. The fitter

1948 uses a low-order chebyshev polynomial to fit out the gradient, with

1949 amp ratios on top of this.

1950

1951 Parameters

1952 ----------

1953 bbox : `lsst.geom.Box2I`

1954 Bounding box for Chebyshev polynomial gradient.

1955 order : `int`

1956 Order of Chebyshev polynomial.

1957 x : `np.ndarray`

1958 Array of x values for points to fit (detector pixels).

1959 y : `np.ndarray`

1960 Array of y values for points to fit (detector pixels).

1961 amp_index : `np.ndarray`

1962 Array of amp numbers associated with each x/y pair.

1963 value : `np.ndarray`

1964 Flat value at each x/y pair.

1965 amps : `np.ndarray`

1966 Array of unique amplifier numbers that will be parameterized.

1967 Any of these amps that does not have any data with the same

1968 amp_index will be set to 1.0.

1969 fixed_amp_index : `int`

1970 Amplifier number to keep fixed.

1971 """

1972 def __init__(self, bbox, order, x, y, amp_index, value, amps, fixed_amp_index):

1973 self._bbox = bbox

1974 self._order = order

1975 self._x = x.astype(np.float64)

1976 self._y = y.astype(np.float64)

1977 self._amp_index = amp_index

1978 self._value = value.astype(np.float64)

1979 self._fixed_amp_index = fixed_amp_index

1980

1981 self.indices = {"chebyshev": np.arange((order + 1) * (order + 1))}

1982 npar = len(self.indices["chebyshev"])

1983

1984 self._amps = amps

1985 self._n_amp = len(self._amps)

1986

1987 self.indices["amp_pars"] = np.arange(self._n_amp) + npar

1988 npar += self._n_amp

1989

1990 self._npar = npar

1991

1992 self._amp_indices = {}

1993 for i in range(self._n_amp):

1994 amp_index = self._amps[i]

1995 self._amp_indices[amp_index] = (self._amp_index == amp_index)

1996

1997 def fit(self, n_iter=10, nsig_clip=5.0):

1998 """Fit the amp ratio parameters.

1999

2000 This uses an iterative fit, where it fits a Chebyshev gradient,

2001 computes amp ratios, and re-fits the gradient.

2002

2003 Parameters

2004 ----------

2005 n_iter : `int`, optional

2006 Number of iterations for fit.

2007 nsig_clip : `float`, optional

2008 Number of sigma in gain correction distribution to clip when

2009 fitting out the Chebyshev gradient.

2010

2011 Returns

2012 -------

2013 pars : `np.ndarray`

2014 Chebyshev parameters and amp offset parameters.

2015 """

2016 value = self._value.copy()

2017

2018 control = lsst.afw.math.ChebyshevBoundedFieldControl()

2019 control.orderX = self._order

2020 control.orderY = self._order

2021 control.triangular = False

2022

2023 pars = np.zeros(self._npar)

2024

2025 pars[self.indices["amp_pars"]] = 1.0

2026

2027 field_use = np.ones(len(value), dtype=np.bool_)

2028

2029 for i in range(n_iter):

2030 field = lsst.afw.math.ChebyshevBoundedField.fit(

2031 self._bbox,

2032 self._x[field_use],

2033 self._y[field_use],

2034 value[field_use],

2035 control,

2036 )

2037

2038 pars[self.indices["chebyshev"][:]] = field.getCoefficients().ravel()

2039

2040 ratio = self._value / field.evaluate(self._x, self._y)

2041

2042 fixed_med = np.median(ratio[self._amp_indices[self._fixed_amp_index]])

2043 ratio /= fixed_med

2044

2045 pars[self.indices["amp_pars"][self._fixed_amp_index]] = 1.0

2046

2047 value = self._value.copy()

2048

2049 for j in range(self._n_amp):

2050 amp_index = self._amps[j]

2051

2052 if np.sum(self._amp_indices[amp_index]) == 0:

2053 continue

2054

2055 pars[self.indices["amp_pars"][j]] = np.median(ratio[self._amp_indices[amp_index]])

2056

2057 amp_pars = pars[self.indices["amp_pars"]]

2058 med = np.median(amp_pars)

2059 sig = median_abs_deviation(amp_pars, scale="normal")

2060 amp_gradient_bad, = np.where(np.abs(amp_pars - med) > nsig_clip * sig)

2061 field_use[:] = True

2062 for agb in amp_gradient_bad:

2063 field_use[self._amp_indices[agb]] = False

2064

2065 return pars

2066

2067 def compute_model(self, pars):

2068 """Compute the gradient/amp ratio model for a given set of parameters.

2069

2070 Parameters

2071 ----------

2072 pars : `np.ndarray`

2073 Chebyshev parameters and amp offset parameters.

2074

2075 Returns

2076 -------

2077 model : `np.ndarray`

2078 The model at each x/y pair.

2079 """

2080 field = lsst.afw.math.ChebyshevBoundedField(

2081 self._bbox,

2082 pars[self.indices["chebyshev"]].reshape(self._order + 1, self._order + 1),

2083 )

2084 model = field.evaluate(self._x, self._y)

2085

2086 for i in range(self._n_amp):

2087 amp_index = self._amps[i]

2088

2089 model[self._amp_indices[amp_index]] *= pars[self.indices["amp_pars"][i]]

2090

2091 return model

2092

2093 def __call__(self, pars):

2094 """Compute the cost function for a set of parameters.

2095

2096 Parameters

2097 ----------

2098 pars : `np.ndarray`

2099 Chebyshev parameters and amp offset parameters.

2100

2101 Returns

2102 -------

2103 cost : `float`

2104 Cost value computed from the absolute deviation.

2105 """

2106 model = self.compute_model(pars)

2107

2108 absdev = np.abs(self._value - model)

2109 t = np.sum(absdev.astype(np.float64))

2110

2111 return t

2112

2113

2114def bin_focal_plane(

2115 exposure_handle_dict,

2116 detector_boundary,

2117 bin_factor,

2118 defect_handle_dict={},

2119 include_itl_flag=True,

2120):

2121 """Bin all the detectors into the full focal plane.

2122

2123 This function reads in images; takes a simple average if there

2124 are more than one input per detector; excludes detector edges;

2125 and bins according to the bin factor. The output is a struct

2126 with focal plane coordinates, detector numbers, and a flag

2127 if the detector is an ITL detector.

2128

2129 Parameters

2130 ----------

2131 exposure_handle_dict : `dict`

2132 Dict keyed by detector (`int`), each element is a list

2133 of `lsst.daf.butler.DeferredDatasetHandle` that will be averaged.

2134 detector_boundary : `int`

2135 Boundary of the detector to ignore (pixels).

2136 bin_factor : `int`

2137 Binning factor. Detectors will be cropped (after applying the

2138 ``detector_boundary``) such that there are no partially

2139 covered binned pixels.

2140 defect_handle_dict : `dict`, optional

2141 Dict keyed by detector (`int`), each element is a

2142 `lsst.daf.butler.DeferredDatasetHandle` for defects to be applied.

2143 include_itl_flag : `bool`, optional

2144 Include a flag for which detectors are ITL?

2145

2146 Returns

2147 -------

2148 binned : `astropy.table.Table`

2149 Table with focal plane positions at the center of each bin

2150 (``xf``, ``yf``); average image values (``value``); and detector

2151 number (``detector``).

2152 """

2153 xf_arrays = []

2154 yf_arrays = []

2155 value_arrays = []

2156 detector_arrays = []

2157 itl_arrays = []

2158

2159 for det in exposure_handle_dict.keys():

2160 flat = exposure_handle_dict[det].get()

2161 defect_handle = defect_handle_dict.get(det, None)

2162 if defect_handle is not None:

2163 defects = defect_handle.get()

2164 else:

2165 defects = None

2166

2167 detector = flat.getDetector()

2168

2169 # Mask out defects if we have them.

2170 if defects is not None:

2171 for defect in defects:

2172 flat.image[defect.getBBox()].array[:, :] = np.nan

2173

2174 # Mask NO_DATA pixels if we have them.

2175 no_data = ((flat.mask.array[:, :] & flat.mask.getPlaneBitMask("NO_DATA")) > 0)

2176 flat.image.array[no_data] = np.nan

2177

2178 # Bin the image, avoiding the boundary and the masked pixels.

2179 # We also make sure we are using an integral number of

2180 # steps to avoid partially covered binned pixels.

2181

2182 arr = flat.image.array

2183

2184 n_step_y = (arr.shape[0] - (2 * detector_boundary)) // bin_factor

2185 y_min = detector_boundary

2186 y_max = bin_factor * n_step_y + y_min

2187 n_step_x = (arr.shape[1] - (2 * detector_boundary)) // bin_factor

2188 x_min = detector_boundary

2189 x_max = bin_factor * n_step_x + x_min

2190

2191 arr = arr[y_min: y_max, x_min: x_max]

2192 binned = arr.reshape((n_step_y, bin_factor, n_step_x, bin_factor))

2193 with warnings.catch_warnings():

2194 warnings.filterwarnings("ignore", r"Mean of empty")

2195 binned = np.nanmean(binned, axis=1)

2196 binned = np.nanmean(binned, axis=2)

2197

2198 xx = np.arange(binned.shape[1]) * bin_factor + bin_factor / 2. + x_min

2199 yy = np.arange(binned.shape[0]) * bin_factor + bin_factor / 2. + y_min

2200 x, y = np.meshgrid(xx, yy)

2201 x = x.ravel()

2202 y = y.ravel()

2203 value = binned.ravel()

2204

2205 # Transform to focal plane coordinates.

2206 transform = detector.getTransform(lsst.afw.cameraGeom.PIXELS, lsst.afw.cameraGeom.FOCAL_PLANE)

2207 xy = np.vstack((x, y))

2208 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)

2209 xf = xf.ravel()

2210 yf = yf.ravel()

2211

2212 if include_itl_flag:

2213 is_itl = np.zeros(len(value), dtype=np.bool_)

2214 # We use this check so that ITL matches ITL science detectors,

2215 # ITL_WF wavefront detectors, and pseudoITL test detectors.

2216 is_itl[:] = ("ITL" in detector.getPhysicalType())

2217

2218 xf_arrays.append(xf)

2219 yf_arrays.append(yf)

2220 value_arrays.append(value)

2221 detector_arrays.append(np.full_like(xf, det, dtype=np.int32))

2222 if include_itl_flag:

2223 itl_arrays.append(is_itl)

2224

2225 xf = np.concatenate(xf_arrays)

2226 yf = np.concatenate(yf_arrays)

2227 value = np.concatenate(value_arrays)

2228 detector = np.concatenate(detector_arrays)

2229

2230 binned = Table(

2231 np.zeros(

2232 len(xf),

2233 dtype=[

2234 ("xf", "f8"),

2235 ("yf", "f8"),

2236 ("value", "f8"),

2237 ("detector", "i4"),

2238 ],

2239 )

2240 )

2241 binned["xf"] = xf

2242 binned["yf"] = yf

2243 binned["value"] = value

2244 binned["detector"] = detector

2245

2246 if include_itl_flag:

2247 binned["itl"] = np.concatenate(itl_arrays).astype(np.bool_)

2248

2249 return binned

2250

2251

2252def bin_flat(ptc, exposure, bin_factor=8, amp_boundary=20, apply_gains=True, gain_ratios=None):

2253 """Bin a flat image, being careful with amplifier edges.

2254

2255 This will optionally apply gains, and apply any gain

2256 ratios.

2257

2258 Parameters

2259 ----------

2260 ptc : `lsst.ip.isr.PhotonTransferCurveDatasets`

2261 PTC dataset with relevant gains.

2262 exposure : `lsst.afw.image.Exposure`

2263 Exposure to bin.

2264 bin_factor : `int`, optional

2265 Binning factor.

2266 amp_boundary : `int`, optional

2267 Boundary around each amp to ignore in binning.

2268 apply_gains : `bool`, optional

2269 Apply gains before binning?

2270 gain_ratios : `np.ndarray`, optional

2271 Array of gain ratios to apply.

2272

2273 Returns

2274 -------

2275 binned : `astropy.table.Table`

2276 Table with detector coordinates at the center of each bin

2277 (``xd``, ``yd``); average image values (``value``); and amplifier

2278 index (``amp_index``).

2279 """

2280 detector = exposure.getDetector()

2281

2282 for i, amp_name in enumerate(ptc.ampNames):

2283 bbox = detector[amp_name].getBBox()

2284 if amp_name in ptc.badAmps:

2285 exposure[bbox].image.array[:, :] = np.nan

2286 continue

2287

2288 if apply_gains:

2289 exposure[bbox].image.array[:, :] *= ptc.gainUnadjusted[amp_name]

2290 if gain_ratios is not None:

2291 exposure[bbox].image.array[:, :] /= gain_ratios[i]

2292

2293 # Next we bin the detector, avoiding amp edges.

2294 xd_arrays = []

2295 yd_arrays = []

2296 value_arrays = []

2297 amp_arrays = []

2298

2299 for i, amp in enumerate(detector):

2300 arr = exposure[amp.getBBox()].image.array

2301

2302 n_step_y = (arr.shape[0] - (2 * amp_boundary)) // bin_factor

2303 y_min = amp_boundary

2304 y_max = bin_factor * n_step_y + y_min

2305 n_step_x = (arr.shape[1] - (2 * amp_boundary)) // bin_factor

2306 x_min = amp_boundary

2307 x_max = bin_factor * n_step_x + x_min

2308

2309 arr = arr[y_min: y_max, x_min: x_max]

2310 binned = arr.reshape((n_step_y, bin_factor, n_step_x, bin_factor))

2311 with warnings.catch_warnings():

2312 warnings.filterwarnings("ignore", r"Mean of empty")

2313 binned = np.nanmean(binned, axis=1)

2314 binned = np.nanmean(binned, axis=2)

2315

2316 xx = np.arange(binned.shape[1]) * bin_factor + bin_factor / 2. + x_min

2317 yy = np.arange(binned.shape[0]) * bin_factor + bin_factor / 2. + y_min

2318 x, y = np.meshgrid(xx, yy)

2319 x = x.ravel() + amp.getBBox().getBeginX()

2320 y = y.ravel() + amp.getBBox().getBeginY()

2321 value = binned.ravel()

2322

2323 xd_arrays.append(x)

2324 yd_arrays.append(y)

2325 value_arrays.append(value)

2326 amp_arrays.append(np.full(len(x), i))

2327

2328 xd = np.concatenate(xd_arrays)

2329 yd = np.concatenate(yd_arrays)

2330 value = np.concatenate(value_arrays)

2331 amp_index = np.concatenate(amp_arrays)

2332

2333 binned = np.zeros(

2334 len(xd),

2335 dtype=[

2336 ("xd", "f8"),

2337 ("yd", "f8"),

2338 ("value", "f8"),

2339 ("amp_index", "i4"),

2340 ],

2341 )

2342 binned["xd"] = xd

2343 binned["yd"] = yd

2344 binned["value"] = value

2345 binned["amp_index"] = amp_index

2346

2347 return Table(data=binned)

2348

2349

2350class ElectrostaticFit():

2351 """

2352 Class to handle the electrostatic fit of area coefficients.

2353

2354 This class manages the fitting of electrostatic model parameters to

2355 measured area change coefficients. The actual electrostatic calculations

2356 are performed by the ElectrostaticCcdGeom class. The fit is performed

2357 using the lmfit Minimizer, and the model includes normalization

2358 parameters (alpha, beta) as well as physical CCD parameters.

2359

2360 Attributes

2361 ----------

2362 fitRange : int

2363 Range of pixels to fit.

2364 inputRange : int

2365 Range of input data.

2366 doFitNormalizationOffset : bool

2367 Whether to fit an offset parameter.

2368 nImageChargePairs : int

2369 Number of image charge pairs used in the electrostatic calculation.

2370 fitMethod : str

2371 Minimization method used by lmfit.

2372 aMatrix : np.ndarray

2373 Measured area change matrix.

2374 aMatrixSigma : np.ndarray

2375 Uncertainty matrix for area changes.

2376 sqrtWeights : np.ndarray

2377 Weights for fitting, inverse of aMatrixSigma.

2378 params : lmfit.Parameters

2379 Fit parameters.

2380 """

2381

2382 def __init__(

2383 self,

2384 initialParams,

2385 fitMethod,

2386 aMatrix,

2387 aMatrixSigma,

2388 fitRange,

2389 doFitNormalizationOffset,

2390 nImageChargePairs,

2391 ):

2392 """

2393 Initialize the ElectrostaticFit class.

2394

2395 Parameters

2396 ----------

2397 initialParams : lmfit.Parameters

2398 Initial fit parameters.

2399 fitMethod : str

2400 Minimization method for lmfit.

2401 aMatrix : np.ndarray

2402 Measured area change matrix.

2403 aMatrixSigma : np.ndarray

2404 Uncertainty matrix for area changes.

2405 fitRange : int

2406 Range of pixels to fit.

2407 doFitNormalizationOffset : bool

2408 Whether to fit an offset parameter.

2409 nImageChargePairs : int

2410 Number of image charge pairs for electrostatic calculation.

2411 """

2412 self.fitRange = fitRange

2413 self.inputRange = aMatrix.shape[0]

2414 self.doFitNormalizationOffset = doFitNormalizationOffset

2415 self.nImageChargePairs = nImageChargePairs

2416 self.fitMethod = fitMethod

2417

2418 if self.fitRange != 0 and self.inputRange < self.fitRange:

2419 print("INFO : truncating input data at %d" % self.inputRange)

2420 self.fitRange = self.inputRange

2421

2422 self.aMatrix = aMatrix[0:fitRange, 0:fitRange]

2423 self.aMatrixSigma = aMatrixSigma[0:fitRange, 0:fitRange]

2424 self.sqrtWeights = 1.0 / aMatrixSigma

2425

2426 self.params = initialParams

2427

2428 def fit(self):

2429 minner = Minimizer(

2430 self.computeWeightedResidual,

2431 self.params,

2432 )

2433

2434 result = minner.minimize(

2435 method=self.fitMethod,

2436 max_nfev=20000,

2437 epsfcn=1.0e-8,

2438 ftol=1.0e-8,

2439 xtol=1.0e-8,

2440 gtol=0.0,

2441 )

2442

2443 return result

2444

2445 def getParamsDict(self):

2446 """

2447 Return a copy of the free parameters vector as a dictionary.

2448 """

2449 return self.params.valuesdict()

2450

2451 def computeWeightedResidual(self, params=None):

2452 self.params = params

2453 m = self.model(params)

2454 w = self.sqrtWeights[:self.fitRange, :self.fitRange]

2455 y = self.aMatrix[:self.fitRange, :self.fitRange]

2456 # Multiply two 2-d arrays term by term:

2457 weightedResiduals = (w * (m - y))

2458 # The result has the same size as both of them.

2459 return weightedResiduals

2460

2461 def model(self, params):

2462 m = self.rawModel(params)

2463 alpha = params['alpha']

2464 beta = params['beta']

2465 return alpha * m + beta

2466

2467 def rawModel(self, params):

2468 # Get all parameters as a dictionary:

2469 parameterDict = self.getParamsDict()

2470 del parameterDict['alpha']

2471 del parameterDict['beta']

2472 del parameterDict['zsh_minus_zq']

2473 del parameterDict['zsv_minus_zq']

2474

2475 # Push them into the electrostatic calculator

2476 c = ElectrostaticCcdGeom(parameterDict)

2477

2478 # Compute the observables:

2479 # If you change the routine called here, also

2480 # change it in BoundaryShifts.__init__()

2481 m = c.evalAreaChangeSidesFast(

2482 self.fitRange,

2483 nImageChargePairs=self.nImageChargePairs,

2484 )

2485 m = m[:self.fitRange, :self.fitRange]

2486

2487 return m

2488

2489 def normalizeModel(self, m):

2490 """

2491 Compute optimal normalization and offset for the model.

2492

2493 The overall normalization is a linear parameter. This

2494 method computes the optimal value given the other

2495 parameters.

2496 """

2497 fr = m.shape[0]

2498 sqrtw = self.sqrtWeights[:fr, :fr]

2499 w = sqrtw ** 2

2500 y = self.aMatrix[:fr, :fr]

2501 if self.doFitNormalizationOffset:

2502 sxx = (w * m * m).sum()

2503 sx = (w * m).sum()

2504 s1 = w.sum()

2505 sxy = (w * m * y).sum()

2506 sy = (w * y).sum()

2507 d = sxx * s1 - sx * sx

2508 a = (s1 * sxy - sx * sy) / d

2509 b = (-sx * sxy + sxx * sy) / d

2510 else:

2511 # Just scale

2512 a = (w * y * m).sum() / (w * m * m).sum()

2513 b = 0

2514 return a, b

2515

2516 def computePixelDistortions(self, conversionWeights=None):

2517 """

2518 Compute pixel distortions using a probability distribution of

2519 conversion depths.

2520

2521 If provided, conversionWeights is expected to be a tuple of

2522 (depth, probability). The routine computes the model corresponding

2523 to this probability distribution. If conversionWeights is not

2524 provided, then [(0, 1.)] is used as the distribution.

2525

2526 Parameters

2527 ----------

2528 conversionWeights : tuple of (np.ndarray, np.ndarray), optional

2529 Tuple containing depths and their associated probabilities.

2530

2531 Returns

2532 -------

2533 BoundaryShifts

2534 The computed boundary shifts for the pixel.

2535 """

2536 zf = self.params["thickness"].value

2537 zsh = self.params["zsh"].value

2538 zsv = self.params["zsv"].value

2539 if conversionWeights is None:

2540 conversionWeights = (np.array([0.0]), np.array([1.0]))

2541

2542 r = None

2543 (d, p) = conversionWeights

2544

2545 # Zero depths lower than the end of drift

2546 too_low = ((zf - d) < zsh) | ((zf - d) < zsv)

2547 p[too_low] = 0

2548 p /= p.sum() # Normalize to 1.

2549 for (depth, prob) in zip(d, p):

2550 if prob == 0:

2551 continue

2552 if r is None:

2553 r = prob * BoundaryShifts(

2554 electrostaticFit=self,

2555 zf=zf - depth

2556 )

2557 else:

2558 r = r + prob * BoundaryShifts(

2559 electrostaticFit=self,

2560 zf=zf - depth

2561 )

2562

2563 return r

2564

2565

2566class BoundaryShifts:

2567 """

2568 Class to compute and store boundary shift values for CCD pixels based on

2569 electrostatic field models.

2570

2571 This class calculates the shifts at the North, South, East, and West

2572 boundaries of each pixel, as well as the area change at the pixel sides,

2573 using parameters from an electrostatic fit and a specified integration

2574 depth `zf`. The shifts are computed using fast integration methods from

2575 the ElectrostaticCcdGeom model.

2576

2577 Attributes

2578 ----------

2579 aN : np.ndarray

2580 Array of North boundary shifts for each pixel.

2581 aS : np.ndarray

2582 Array of South boundary shifts for each pixel.

2583 aE : np.ndarray

2584 Array of East boundary shifts for each pixel.

2585 aW : np.ndarray

2586 Array of West boundary shifts for each pixel.

2587 ath : np.ndarray

2588 Array of area changes at the pixel sides, including beta.

2589 athMinusBeta : np.ndarray

2590 Array of area changes at the pixel sides, excluding beta.

2591 """

2592

2593 def __init__(self, electrostaticFit, zf):

2594 if zf <= 0:

2595 raise RuntimeError(

2596 f"Cannot integrate below the bottom of the pixel, zf={zf}."

2597 )

2598

2599 paramDict = electrostaticFit.getParamsDict()

2600 del paramDict['alpha']

2601 del paramDict['beta']

2602 del paramDict['zsh_minus_zq']

2603 del paramDict['zsv_minus_zq']

2604

2605 c = ElectrostaticCcdGeom(paramDict)

2606 fr = electrostaticFit.fitRange

2607 ii, jj = np.meshgrid(list(range(fr)), list(range(fr)))

2608 ii = ii.flatten()

2609 jj = jj.flatten()

2610 self.aN = np.ndarray((fr, fr))

2611 self.aS = np.zeros_like(self.aN)

2612 self.aE = np.zeros_like(self.aN)

2613 self.aW = np.zeros_like(self.aN)

2614 self.ath = np.zeros_like(self.aN)

2615 alpha = electrostaticFit.params['alpha'].value

2616 beta = electrostaticFit.params['beta'].value

2617

2618 # alpha*rawModel+beta is the description of the measurements

2619 # We should not apply beta to the outcome, because beta is meant to

2620 # accomodate some long-range contamination (non-electrostatic)

2621 # of the covariance measurements.

2622 iMax, jMax = fr, fr

2623 self.aN = -alpha * c.integrateEyFieldFast(iMax, jMax, 1, zf=zf)

2624 self.aS = -alpha * c.integrateEyFieldFast(iMax, jMax, -1, zf=zf)

2625 self.aW = -alpha * c.integrateExFieldFast(iMax, jMax, -1, zf=zf)

2626 self.aE = -alpha * c.integrateExFieldFast(iMax, jMax, 1, zf=zf)

2627 self.ath = alpha * c.evalAreaChangeSidesFast(fr, zf=zf) + beta

2628

2629 self.athMinusBeta = self.ath - beta

2630

2631 def __rmul__(self, factor):

2632 """

2633 """

2634 res = copy.deepcopy(self)

2635 res.aN *= factor

2636 res.aS *= factor

2637 res.aE *= factor

2638 res.aW *= factor

2639 res.ath *= factor

2640 res.athMinusBeta *= factor

2641 return res

2642

2643 def __add__(self, other):

2644 """

2645 """

2646 res = copy.deepcopy(self)

2647 res.aN += other.aN

2648 res.aS += other.aS

2649 res.aE += other.aE

2650 res.aW += other.aW

2651 res.ath += other.ath

2652 res.athMinusBeta += other.athMinusBeta

2653 return res

2654

2655

2656def calcEFieldCoulomb(xv, xqv):

2657 """

2658 xv = where, xqv = charge location.

2659 both should be numpy arrays.

2660 if qv is multi-d, the routine assumes that the

2661 physical coordinates (x,y,z) are patrolled by

2662 the last index

2663 """

2664 d = xv - xqv

2665 r3 = np.power((d**2).sum(axis=-1), 1.5)

2666 # anything more clever ?

2667 # of course d/r3 does not work

2668 return (d.T / r3.T).T

2669

2670

2671class ElectrostaticCcdGeom():

2672 def __init__(self, paramDict):

2673 """

2674 Dictionary of parameters containing the following keys:

2675 parameters : (all in microns)

2676 zq : altitude of the burried channel (microns)

2677 zsh : vertex altitude for horizontal boundaries

2678 zsv : vertex altitude for vertical boundaries

2679 a, b : size of the rectangular charge source

2680 thickness : thickness

2681 pixelsize : pixel size

2682 """

2683 zq = paramDict['zq']

2684 zsh = paramDict['zsh']

2685 zsv = paramDict['zsv']

2686 a = paramDict['a']

2687 b = paramDict['b']

2688 thickness = paramDict['thickness']

2689 pixelsize = paramDict['pixelsize']

2690

2691 self.zq = np.fabs(zq)

2692 self.zsh = np.fabs(zsh)

2693 self.zsv = np.fabs(zsv)

2694 self.b = np.fabs(b)

2695 self.a = np.fabs(a)

2696 self.thickness = np.fabs(thickness)

2697 self.pixelsize = pixelsize

2698

2699 self.nStepsZ = 100

2700 self.nStepsXY = 20

2701

2702 # yields a ~ 1% precision of the field at z~10

2703 # if compared to the uniform sheet model (eField

2704 self.nChargeElements = 3

2705

2706 # Set up weights

2707 x, w = leggauss(self.nStepsXY)

2708 self.integrationWeights = w*0.5

2709

2710 # Abcissa refer to [-1, 1], we want [0, self.pixelsize]

2711 self.xyOffsets = (x+1) * 0.5 * self.pixelsize

2712

2713 def calcEFieldCoulombChargeSheet(self, xv, xqv):

2714 """

2715 xv = where (the last index should address x,y,z.

2716 xqv = charge location

2717 Both Should be numpy arrays.

2718 if xv is multi-d, the routine assumes

2719 that the physical coordinates (x,y,z) are patrolled by the last index.

2720 Returns the electric field from a unitely charged horizontal rectangle

2721 centered at x_q of size 2a * 2b.

2722 The returned electric field assumes 4*pi*epsilon=1.

2723 """

2724 # Use Durand page 244 tome 1

2725 # four corners :

2726 x1 = xqv + np.array([self.a, self.b, 0])

2727 x2 = xqv + np.array([-self.a, self.b, 0])

2728 x3 = xqv + np.array([-self.a, -self.b, 0])

2729 x4 = xqv + np.array([self.a, -self.b, 0])

2730

2731 # Distances to the four corners

2732 d1 = np.sqrt(((xv-x1)**2).sum(axis=-1))

2733 d2 = np.sqrt(((xv-x2)**2).sum(axis=-1))

2734 d3 = np.sqrt(((xv-x3)**2).sum(axis=-1))

2735 d4 = np.sqrt(((xv-x4)**2).sum(axis=-1))

2736

2737 # Reserve the returned array

2738 r = np.ndarray(xv.shape)

2739 xvOffset = xv[..., 0] - xqv[0]

2740 yvOffset = xv[..., 1] - xqv[1]

2741

2742 # Ex

2743 # note : if a or b goes to 0, the log is 0 and the denominator (last

2744 # line) is zero as well. So some expansion would be required

2745 ao = yvOffset + self.b

2746 bo = yvOffset - self.b

2747 co = xvOffset + self.a

2748 do = xvOffset - self.a

2749 # Ex (eq 105)

2750 r[..., 0] = np.log((d4+ao) * (d2+bo) / (d3+ao) / (d1+bo))

2751 # Ey (eq 106)

2752 r[..., 1] = np.log((d2+co) * (d4+do) / (d3+co) / (d1+do))

2753 # point source approximation

2754 # ret[...,2] = (4*self.a*self.b)*calcEFieldCoulomb(X,X_q)[...,2]

2755 # full expression for ez : p 244

2756 # ez (eq 111 is only valid if the x and y are "inside")

2757 # there is a discussion of the general case around Fig VI-18.

2758 z = xv[..., 2] - xqv[2]

2759 # it ressembles equation 111 but I flipped two signs

2760 r[..., 2] = (np.arctan(do * bo / z / d1)

2761 - np.arctan(bo * co / z / d2)

2762 + np.arctan(co * ao / z / d3)

2763 - np.arctan(ao * do / z / d4))

2764

2765 # seems OK both "inside" and "outside"

2766 return r / (4 * self.a * self.b)

2767

2768 def integrateEField(self, xv, component, z0, zf, nImageChargePairs=11):

2769 """

2770 Integrate transverse E Field along Z at point X (2 coordinates, last

2771 index). The coordinate of the field is given by component (0,

2772 or 1). at point X from the point charge The computation uses

2773 the dipole series trick. The number of dipoles is an optional

2774 argument. Odd numbers are better for what we are doing here.

2775

2776 """

2777 # The integral of the field (x/r^3 dz from z1 to z2) reads

2778 # x/rho**2*(z2/r2-z1/r1) with rho2 = x**2+y**2

2779 # x/rho2 does not change when going through image sources

2780 # so we use them as arguments, dz1 and dz2 z{begin,end}--Xq[2]

2781 # just for test: if z0==zf, then return the field value

2782 if z0 != zf:

2783 def integral(rho2, x, dz1, dz2):

2784 r1 = np.sqrt(rho2 + dz1**2)

2785 r2 = np.sqrt(rho2 + dz2**2)

2786 return x * (dz2/r2 - dz1/r1) / rho2

2787 else: # see the comment above: return the value, not the integral.

2788 def integral(rho2, x, dz1, dz2):

2789 """

2790 x/r**3

2791 """

2792 r = np.sqrt(rho2 + dz1**2)

2793 vals = x / r**3

2794 return vals

2795

2796 # Check validity

2797 if component not in [0, 1]:

2798 raise RuntimeError("Can only integrate electric field along component 0 or 1, "

2799 f"not {component}.")

2800

2801 # Reserve the result array

2802 result = np.zeros(xv.shape[:-1])

2803 zqp = self.zq

2804 zqm = -zqp

2805

2806 # 1st dipole

2807 # Generate a set of point charges to emulate

2808 # an extended distribution (size 2a*2b)

2809 xStep = (2 * self.a) / self.nChargeElements

2810 yStep = (2 * self.b) / self.nChargeElements

2811 xqs = -self.a + (np.linspace(0, self.nChargeElements-1, self.nChargeElements) + 0.5) * xStep

2812 yqs = -self.b + (np.linspace(0, self.nChargeElements-1, self.nChargeElements) + 0.5) * yStep

2813 for xq in xqs:

2814 for yq in yqs:

2815 dx = xv[..., 0] - xq

2816 dy = xv[..., 1] - yq

2817 dxv = [dx, dy]

2818 rho2 = dx**2 + dy**2

2819 result += integral(rho2, dxv[component], z0 - zqp, zf - zqp)

2820 # Construct image charge, change sign of z and q

2821 result -= integral(rho2, dxv[component], z0 - zqm, zf - zqm)

2822 result /= (self.nChargeElements * self.nChargeElements)

2823

2824 # Remaining dipoles : no more extended charge

2825 # The (x,y) charge coordinates are 0, and common to all images:

2826 rho2 = xv[..., 0]**2 + xv[..., 1]**2

2827 xi = xv[..., component]

2828 for i in range(1, nImageChargePairs):

2829 odd = i%2

2830 if odd:

2831 ztmp = (2*self.thickness) - zqm

2832 zqm = (2*self.thickness) - zqp

2833 zqp = ztmp

2834 else:

2835 ztmp = -zqm

2836 zqm = -zqp

2837 zqp = ztmp

2838 result += integral(rho2, xi, z0 - zqp, zf - zqp)

2839 result -= integral(rho2, xi, z0 - zqm, zf - zqm)

2840

2841 # 8.85418781e-12 (F/m) * 1e-6 (microns/m) / 1.602e-19 (C/electron)

2842 # = 55.26958

2843 # eps_r_Si = 12, so eps = 55.27*12 = 663.23 el/V/um

2844 # This routine hence returns the field sourced by -1 electron

2845 result *= 1 / (4*np.pi*660)

2846

2847 return result

2848

2849 def eField(self, xv, nImageChargePairs=11):

2850 """

2851 Field at point xv from the point charge

2852 if xv is multi-dimensional, x,y,z should be represented

2853 by the last index ([0:3]).

2854 The computation uses the dipole series trick. The number

2855 of dipoles is an optional argument. Odd numbers are better

2856 for what we are doing here.

2857 """

2858 # Put the center of the source pixel at (x, y) = (0, 0)

2859 # This assumption is relied on in Eval_Eth and Eval_Etv

2860 xqv = np.array([0, 0, self.zq])

2861

2862 # Image charge w.r.t. the parallel clock lines

2863 xqvImage = np.array([xqv[0], xqv[1], -xqv[2]])

2864

2865 # Split the calculation in 2 parts: approximation

2866 # when far from the source, image method when near.

2867 rho = np.sqrt(xv[..., 0]**2 + xv[..., 1]**2)

2868

2869 # Field near the charge cloud

2870 # First dipole

2871 # The points lower than half the thickness of the pixel

2872 # are considered "close"

2873 pointsNearChargeCloudIdxs = (rho/self.thickness < 2)

2874 xvNear = xv[pointsNearChargeCloudIdxs]

2875 eFieldNear = self.calcEFieldCoulombChargeSheet(xvNear, xqv)

2876 eFieldNear -= self.calcEFieldCoulombChargeSheet(xvNear, xqvImage)

2877

2878 # Next dipoles

2879 for i in range(1, nImageChargePairs):

2880 if (i % 2): # Is odd?

2881 xqv[2] = 2 * self.thickness - xqv[2]

2882 xqvImage[2] = 2 * self.thickness - xqvImage[2]

2883 eFieldNear += calcEFieldCoulomb(xvNear, xqvImage)

2884 eFieldNear -= calcEFieldCoulomb(xvNear, xqv)

2885 else:

2886 xqv[2] = -xqv[2]

2887 xqvImage[2] = -xqvImage[2]

2888 eFieldNear += calcEFieldCoulomb(xvNear, xqv)

2889 eFieldNear -= calcEFieldCoulomb(xvNear, xqvImage)

2890

2891 # Field far from the charge cloud

2892 xvFar = xv[~pointsNearChargeCloudIdxs]

2893 rhoFar = rho[~pointsNearChargeCloudIdxs]

2894

2895 # Jon Pumplin, Am. Jour. Phys. 37,7 (1969), eq 7

2896 # When changing coordinate system (shift along z), cos -> sin.

2897 # And since this only applies far from the source, the latter

2898 # can be regarded as point-like.

2899 # I checked the continuity over the separation point.

2900 factor = -np.sin(np.pi * self.zq / self.thickness)

2901 factor *= np.sin(np.pi * xvFar[..., 2] / self.thickness)

2902 factor *= np.exp(-np.pi * rhoFar / self.thickness)

2903 factor *= np.sqrt(8 / rhoFar / self.thickness)

2904 factor *= (-np.pi / self.thickness - 0.5/rhoFar)

2905 eFieldFar = np.zeros_like(xvFar)

2906 eFieldFar[..., 0] = xvFar[..., 0] / (rhoFar*factor)

2907 eFieldFar[..., 1] = xvFar[..., 1] / (rhoFar*factor)

2908

2909 # aggregate the results

2910 eField = np.zeros_like(xv)

2911 eField[pointsNearChargeCloudIdxs] = eFieldNear

2912 eField[~pointsNearChargeCloudIdxs] = eFieldFar

2913

2914 # epsilon0 = 55 el/V/micron

2915 # 8.85418781e-12 (F/m) * 1e-6 (microns/m) / 1.602e-19 (C/electron)

2916 # = 55.26958

2917 # eps_r_Si = 12, so eps = 55.27*12 = 663.23 el/V/um

2918 # This routine hence returns the field sourced by -1 electron

2919 eField *= 1 / (4*np.pi*660)

2920

2921 return eField

2922

2923 def integrateExFieldFast(self, iMax, jMax, leftOrRight, zf=None, nImageChargePairs=11):

2924 """

2925 Computes the integrals of Ex along z from

2926 self.zsh to zf. The returned array is 2d of

2927 shape (iMax, jMax).

2928 """

2929 zf = self.thickness if zf is None else zf

2930 assert zf > self.zsv

2931 ii, jj = np.indices((iMax, jMax))

2932 yy = (jj - 0.5)[:, :, np.newaxis] * self.pixelsize + self.xyOffsets[np.newaxis, np.newaxis, :]

2933 xx = (ii + 0.5 * leftOrRight) * self.pixelsize

2934 xx = np.broadcast_to(xx[..., None], yy.shape)

2935 xv = np.stack([xx, yy], axis=-1)

2936

2937 # By definition of zsh, we integrate from zsv to zf,

2938 # and divide by the pixel size to be consistent with Eval_ET{v,h}

2939 w = np.broadcast_to(self.integrationWeights, yy.shape)

2940 integral = np.sum(

2941 w * self.integrateEField(xv, 0, self.zsv, zf,

2942 nImageChargePairs=nImageChargePairs),

2943 axis=2,

2944 )

2945 return leftOrRight * integral/self.pixelsize

2946

2947 def integrateEyFieldFast(self, iMax, jMax, topOrBottom, zf=None, nImageChargePairs=11):

2948 """

2949 Computes the integral of Ey along z from

2950 self.zsh to zf. The returned array is

2951 2d [0:iMax, 0:jMax].

2952 """

2953 zf = self.thickness if zf is None else zf

2954 # assert zf > self.zsh

2955 ii, jj = np.indices((iMax, jMax))

2956 xx = (ii - 0.5)[:, :, np.newaxis] * self.pixelsize + self.xyOffsets[np.newaxis, np.newaxis, :]

2957 yy = (jj + 0.5*topOrBottom) * self.pixelsize

2958 yy = np.broadcast_to(yy[..., None], xx.shape)

2959 xv = np.stack([xx, yy], axis=-1)

2960

2961 # By definition of zsh, we integrate from zsh to zf,

2962 # and divide by the pixel size to be consistent with Eval_ET{v,h}

2963 # integrate

2964 w = np.broadcast_to(self.integrationWeights, yy.shape) # Add leading dimensions

2965 integral = np.sum(

2966 w * self.integrateEField(xv, 1, self.zsh, zf,

2967 nImageChargePairs=nImageChargePairs),

2968 axis=2,

2969 )

2970 return topOrBottom * integral / self.pixelsize

2971

2972 def eFieldTransverseToHorizontalPixelEdgeGrid(self, i, j, topOrBottom, zf=None):

2973 """

2974 Returns the field transverse to the horizontal pixel boundary.

2975 return a 2d array of shifts at evenly spaced points in x and z.

2976 normalized in units of pixel size, for a unit charge.

2977 The returned array has 3 indices:

2978 [along the pixel side, along the drift,E-field coordinate].

2979 The resutl is multiplied by the z- and x- steps, so that the

2980 sum is the averge over x, divided by the pixel size.

2981 """

2982 assert np.abs(topOrBottom) == 1

2983 zf = self.thickness if zf is None else zf

2984

2985 # by definition of zsh, we integrate from zsh to zf:

2986 zStep = (zf - self.zsh) / self.nStepsZ

2987 xyStep = self.pixelsize / self.nStepsXY

2988

2989 z = self.zsh + (np.linspace(0, self.nStepsZ - 1, num=self.nStepsZ) + 0.5) * zStep

2990 x = (i - 0.5)*self.pixelsize + (np.linspace(0, self.nStepsXY - 1, self.nStepsXY) + 0.5) * xyStep

2991

2992 xx, zz = np.meshgrid(x, z)

2993 yy = np.ones(xx.shape) * (j + 0.5*topOrBottom) * self.pixelsize

2994 xv = np.array([xx, yy, zz]).T

2995

2996 return self.eField(xv) * zStep * xyStep

2997

2998 def averageHorizontalBoundaryShift(self, i, j, topOrBottom, zf=None):

2999 """

3000 Integrate the field transverse to the horizontal pixel boundary

3001 """

3002 sum = self.eFieldTransverseToHorizontalPixelEdgeGrid(i, j, topOrBottom, zf)[..., 1].sum()

3003 # we want the integral over z and the average over x,

3004 # divided by the pixel size, with a sign that defines

3005 # if in moves inside or outside. Here is it:

3006 return sum * (topOrBottom / (self.pixelsize**2))

3007

3008 def corner_shift_h2(self, i, j, topOrBottom, zf=None): # Do we even need this??

3009 E = self.eFieldTransverseToHorizontalPixelEdgeGrid(i, j, topOrBottom, zf)[..., 1]

3010 # Integrate over z (multiplication by the

3011 # step done in the calling routine)

3012 intz = E.sum(axis=1)

3013 x = range(intz.shape[0])

3014 p = np.polyfit(x, intz, 1)

3015 return p[0] * -0.5 + p[1], p[0] * (x[-1] + 0.5) + p[1]

3016

3017 def eFieldTransverseToVerticalPixelEdgeGrid(self, i, j, leftOrRight, zf=None):

3018 """

3019 Returns the field transverse to the vertical pixel boundary.

3020 return a 3d array of shifts at evenly spaced points in y and z.

3021 Ex,y,z is indexed by the last index.

3022 """

3023 # Determine the depth of the photon conversion

3024 zf = self.thickness if zf is None else zf

3025

3026 # The source charge is at x,y=0

3027 zStep = (zf-self.zsv) / (self.nStepsZ)

3028 xyStep = self.pixelsize / (self.nStepsXY)

3029

3030 z = self.zsv + (np.linspace(0, self.nStepsZ - 1, self.nStepsZ) + 0.5) * zStep

3031 y = (j - 0.5)*self.pixelsize + (np.linspace(0, self.nStepsXY - 1, self.nStepsXY) + 0.5) * xyStep

3032

3033 yy, zz = np.meshgrid(y, z)

3034 xx = np.ones(yy.shape) * (i + 0.5*leftOrRight) * self.pixelsize

3035 xv = np.array([xx, yy, zz]).T

3036

3037 return self.eField(xv) * zStep * xyStep

3038

3039 def averageVerticalBoundaryShift(self, i, j, leftOrRight, zf=None):

3040 """

3041 Average shift of the vertical boundary of pixel (i j).

3042 """

3043 sum = self.eFieldTransverseToVerticalPixelEdgeGrid(i, j, leftOrRight, zf)[..., 0].sum()

3044 # we want the integral over z and the average over x,

3045 # divided by the pixel size, with a sign that defines

3046 # if in moves inside or outside. Here is it:

3047 return sum * (leftOrRight / (self.pixelsize**2))

3048

3049 def dxdy(self, iMax, zf=None):

3050 """

3051 corner shifts calculations. The returned array are larger by 1

3052 than iMax, because there are more corners than pixels

3053 """

3054 horiztonalShifts = np.ndarray((iMax, iMax))

3055 verticalShifts = np.ndarray((iMax, iMax))

3056 for i in range(iMax):

3057 for j in range(iMax):

3058 horiztonalShifts[i, j] = self.averageHorizontalBoundaryShift(i + 0.5, j, 1, zf)

3059 verticalShifts[i, j] = self.averageVerticalBoundaryShift(i, j + 0.5, 1, zf)

3060

3061 # Parametrize the corner shifts of iMax pixels: iMax+1

3062 # corners in each direction

3063 dx = np.zeros((iMax+1, iMax+1))

3064 dy = dx + 0.

3065 dx[1:, 1:] = verticalShifts

3066 dx[0, 1:] = -verticalShifts[0, :] # leftmost column

3067 dx[:, 0] = dx[:, 1] # bottom row

3068 dy[1:, 1:] = horiztonalShifts

3069 dy[1:, 0] = -horiztonalShifts[:, 0] # bottom row

3070 dy[0, :] = dy[1, :] # leftmost column

3071

3072 return dx, dy

3073

3074 def evalAreaChangeCorners(self, iMax, zf=None):

3075 """

3076 pixel area alterations computed through corner shifts

3077 """

3078 dx, dy = self.dxdy(iMax, zf)

3079 areaChange = dx[1:, 1:] - dx[:-1, 1:] + dx[1:, :-1] - dx[:-1, :-1]

3080 areaChange += dy[1:, 1:] - dy[1:, :-1] + dy[:-1, 1:] - dy[:-1, :-1]

3081 return -0.5*areaChange

3082

3083 def evalAreaChangeSidesFast(self, iMax, zf=None, nImageChargePairs=11):

3084 """

3085 Same as EvalAreaChangeSides, but uses direct integration.

3086 it evaluates the divergence of the discrete boundary

3087 displacement field.

3088 This routine groups the calls to the field computing routines.

3089 """

3090 horiztonalShifts = np.ndarray((iMax, iMax+1))

3091 verticalShifts = np.zeros_like(horiztonalShifts.T)

3092 horiztonalShifts[:, 1:] = self.integrateEyFieldFast(iMax, iMax, 1, zf,

3093 nImageChargePairs=nImageChargePairs)

3094 verticalShifts[1:, :] = self.integrateExFieldFast(iMax, iMax, 1, zf,

3095 nImageChargePairs=nImageChargePairs)

3096

3097 # special case for [0,j] and [i,0] (they have two opposite values)

3098 horiztonalShifts[:, 0] = -horiztonalShifts[:, 1]

3099 verticalShifts[0, :] = -verticalShifts[1, :]

3100

3101 # the divergence

3102 areaChange = verticalShifts[1:, :]

3103 areaChange -= verticalShifts[:-1, :]

3104 areaChange += horiztonalShifts[:, 1:]

3105 areaChange -= horiztonalShifts[:, :-1]

3106 areaChange *= -1 # Pixel area is decreased

3107

3108 return areaChange

Coverage for python / lsst / cp / pipe / utils.py: 7%

1126 statements