Coverage for python / lsst / ip / diffim / psfMatch.py: 24%

1# This file is part of ip_diffim. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22__all__ = ["PsfMatchConfig", "PsfMatchConfigAL", "PsfMatchConfigDF", "PsfMatchTask"] 

23 

24import abc 

25import time 

26 

27import numpy as np 

28 

29import lsst.afw.image as afwImage 

30import lsst.pex.config as pexConfig 

31import lsst.afw.math as afwMath 

32import lsst.afw.display as afwDisplay 

33import lsst.pipe.base as pipeBase 

34from lsst.meas.algorithms import SubtractBackgroundConfig 

35from lsst.utils.logging import getTraceLogger 

36from lsst.utils.timer import timeMethod 

37from . import utils as diutils 

38from . import diffimLib 

39 

40 

41class NoKernelCandidatesError(pipeBase.AlgorithmError): 

42 """Raised if there are too few candidates to compute the PSF matching 

43 kernel. 

44 """ 

45 

46 @property 

47 def metadata(self) -> dict: 

48 return {} 

49 

50 

51class PsfMatchConfig(pexConfig.Config): 

52 """Base configuration for Psf-matching 

53 

54 The base configuration of the Psf-matching kernel, and of the warping, detection, 

55 and background modeling subTasks.""" 

56 

57 warpingConfig = pexConfig.ConfigField("Config for warping exposures to a common alignment", 

58 afwMath.WarperConfig) 

59 afwBackgroundConfig = pexConfig.ConfigField("Controlling the Afw background fitting", 

60 SubtractBackgroundConfig) 

61 

62 useAfwBackground = pexConfig.Field( 

63 dtype=bool, 

64 doc="Use afw background subtraction instead of ip_diffim", 

65 default=False, 

66 ) 

67 fitForBackground = pexConfig.Field( 

68 dtype=bool, 

69 doc="Include terms (including kernel cross terms) for background in ip_diffim", 

70 default=False, 

71 ) 

72 kernelBasisSet = pexConfig.ChoiceField( 

73 dtype=str, 

74 doc="Type of basis set for PSF matching kernel.", 

75 default="alard-lupton", 

76 allowed={ 

77 "alard-lupton": """Alard-Lupton sum-of-gaussians basis set, 

78 * The first term has no spatial variation 

79 * The kernel sum is conserved 

80 * You may want to turn off 'usePcaForSpatialKernel'""", 

81 "delta-function": """Delta-function kernel basis set, 

82 * You may enable the option useRegularization 

83 * You should seriously consider usePcaForSpatialKernel, which will also 

84 enable kernel sum conservation for the delta function kernels""" 

85 } 

86 ) 

87 kernelSize = pexConfig.Field( 

88 dtype=int, 

89 doc="""Number of rows/columns in the convolution kernel; should be odd-valued. 

90 Modified by kernelSizeFwhmScaling if scaleByFwhm = true""", 

91 default=21, 

92 ) 

93 scaleByFwhm = pexConfig.Field( 

94 dtype=bool, 

95 doc="Scale kernelSize, alardGaussians by input Fwhm", 

96 default=True, 

97 ) 

98 kernelSizeFwhmScaling = pexConfig.Field( 

99 dtype=float, 

100 doc="Multiplier of the largest AL Gaussian basis sigma to get the kernel bbox (pixel) size.", 

101 default=6.0, 

102 check=lambda x: x >= 1.0 

103 ) 

104 kernelSizeMin = pexConfig.Field( 

105 dtype=int, 

106 doc="Minimum kernel bbox (pixel) size.", 

107 default=21, 

108 ) 

109 kernelSizeMax = pexConfig.Field( 

110 dtype=int, 

111 doc="Maximum kernel bbox (pixel) size.", 

112 default=35, 

113 ) 

114 spatialModelType = pexConfig.ChoiceField( 

115 dtype=str, 

116 doc="Type of spatial functions for kernel and background", 

117 default="chebyshev1", 

118 allowed={ 

119 "chebyshev1": "Chebyshev polynomial of the first kind", 

120 "polynomial": "Standard x,y polynomial", 

121 } 

122 ) 

123 spatialKernelOrder = pexConfig.Field( 

124 dtype=int, 

125 doc="Spatial order of convolution kernel variation", 

126 default=2, 

127 check=lambda x: x >= 0 

128 ) 

129 spatialBgOrder = pexConfig.Field( 

130 dtype=int, 

131 doc="Spatial order of differential background variation", 

132 default=1, 

133 check=lambda x: x >= 0 

134 ) 

135 sizeCellX = pexConfig.Field( 

136 dtype=int, 

137 doc="Size (rows) in pixels of each SpatialCell for spatial modeling", 

138 default=128, 

139 check=lambda x: x >= 32 

140 ) 

141 sizeCellY = pexConfig.Field( 

142 dtype=int, 

143 doc="Size (columns) in pixels of each SpatialCell for spatial modeling", 

144 default=128, 

145 check=lambda x: x >= 32 

146 ) 

147 nStarPerCell = pexConfig.Field( 

148 dtype=int, 

149 doc="Maximum number of KernelCandidates in each SpatialCell to use in the spatial fitting. " 

150 "Set to -1 to use all candidates in each cell.", 

151 default=5, 

152 ) 

153 maxSpatialIterations = pexConfig.Field( 

154 dtype=int, 

155 doc="Maximum number of iterations for rejecting bad KernelCandidates in spatial fitting", 

156 default=3, 

157 check=lambda x: x >= 1 and x <= 5 

158 ) 

159 usePcaForSpatialKernel = pexConfig.Field( 

160 dtype=bool, 

161 doc="""Use Pca to reduce the dimensionality of the kernel basis sets. 

162 This is particularly useful for delta-function kernels. 

163 Functionally, after all Cells have their raw kernels determined, we run 

164 a Pca on these Kernels, re-fit the Cells using the eigenKernels and then 

165 fit those for spatial variation using the same technique as for Alard-Lupton kernels. 

166 If this option is used, the first term will have no spatial variation and the 

167 kernel sum will be conserved.""", 

168 default=False, 

169 ) 

170 subtractMeanForPca = pexConfig.Field( 

171 dtype=bool, 

172 doc="Subtract off the mean feature before doing the Pca", 

173 default=True, 

174 ) 

175 numPrincipalComponents = pexConfig.Field( 

176 dtype=int, 

177 doc="""Number of principal components to use for Pca basis, including the 

178 mean kernel if requested.""", 

179 default=5, 

180 check=lambda x: x >= 3 

181 ) 

182 singleKernelClipping = pexConfig.Field( 

183 dtype=bool, 

184 doc="Do sigma clipping on each raw kernel candidate", 

185 default=True, 

186 ) 

187 kernelSumClipping = pexConfig.Field( 

188 dtype=bool, 

189 doc="Do sigma clipping on the ensemble of kernel sums", 

190 default=True, 

191 ) 

192 spatialKernelClipping = pexConfig.Field( 

193 dtype=bool, 

194 doc="Do sigma clipping after building the spatial model", 

195 default=True, 

196 ) 

197 checkConditionNumber = pexConfig.Field( 

198 dtype=bool, 

199 doc="""Test for maximum condition number when inverting a kernel matrix. 

200 Anything above maxConditionNumber is not used and the candidate is set as BAD. 

201 Also used to truncate inverse matrix in estimateBiasedRisk. However, 

202 if you are doing any deconvolution you will want to turn this off, or use 

203 a large maxConditionNumber""", 

204 default=False, 

205 ) 

206 badMaskPlanes = pexConfig.ListField( 

207 dtype=str, 

208 doc="""Mask planes to ignore when calculating diffim statistics 

209 Options: NO_DATA EDGE SAT BAD CR INTRP""", 

210 default=("NO_DATA", "EDGE", "SAT") 

211 ) 

212 candidateResidualMeanMax = pexConfig.Field( 

213 dtype=float, 

214 doc="""Rejects KernelCandidates yielding bad difference image quality. 

215 Used by BuildSingleKernelVisitor, AssessSpatialKernelVisitor. 

216 Represents average over pixels of (image/sqrt(variance)).""", 

217 default=0.25, 

218 check=lambda x: x >= 0.0 

219 ) 

220 candidateResidualStdMax = pexConfig.Field( 

221 dtype=float, 

222 doc="""Rejects KernelCandidates yielding bad difference image quality. 

223 Used by BuildSingleKernelVisitor, AssessSpatialKernelVisitor. 

224 Represents stddev over pixels of (image/sqrt(variance)).""", 

225 default=1.50, 

226 check=lambda x: x >= 0.0 

227 ) 

228 useCoreStats = pexConfig.Field( 

229 dtype=bool, 

230 doc="""Use the core of the footprint for the quality statistics, instead of the entire footprint. 

231 WARNING: if there is deconvolution we probably will need to turn this off""", 

232 default=False, 

233 ) 

234 candidateCoreRadius = pexConfig.Field( 

235 dtype=int, 

236 doc="""Radius for calculation of stats in 'core' of KernelCandidate diffim. 

237 Total number of pixels used will be (2*radius)**2. 

238 This is used both for 'core' diffim quality as well as ranking of 

239 KernelCandidates by their total flux in this core""", 

240 default=3, 

241 check=lambda x: x >= 1 

242 ) 

243 maxKsumSigma = pexConfig.Field( 

244 dtype=float, 

245 doc="""Maximum allowed sigma for outliers from kernel sum distribution. 

246 Used to reject variable objects from the kernel model""", 

247 default=3.0, 

248 check=lambda x: x >= 0.0 

249 ) 

250 maxConditionNumber = pexConfig.Field( 

251 dtype=float, 

252 doc="Maximum condition number for a well conditioned matrix", 

253 default=5.0e7, 

254 check=lambda x: x >= 0.0 

255 ) 

256 conditionNumberType = pexConfig.ChoiceField( 

257 dtype=str, 

258 doc="Use singular values (SVD) or eigen values (EIGENVALUE) to determine condition number", 

259 default="EIGENVALUE", 

260 allowed={ 

261 "SVD": "Use singular values", 

262 "EIGENVALUE": "Use eigen values (faster)", 

263 } 

264 ) 

265 maxSpatialConditionNumber = pexConfig.Field( 

266 dtype=float, 

267 doc="Maximum condition number for a well conditioned spatial matrix", 

268 default=1.0e10, 

269 check=lambda x: x >= 0.0 

270 ) 

271 iterateSingleKernel = pexConfig.Field( 

272 dtype=bool, 

273 doc="""Remake KernelCandidate using better variance estimate after first pass? 

274 Primarily useful when convolving a single-depth image, otherwise not necessary.""", 

275 default=False, 

276 ) 

277 constantVarianceWeighting = pexConfig.Field( 

278 dtype=bool, 

279 doc="""Use constant variance weighting in single kernel fitting? 

280 In some cases this is better for bright star residuals.""", 

281 default=True, 

282 ) 

283 calculateKernelUncertainty = pexConfig.Field( 

284 dtype=bool, 

285 doc="""Calculate kernel and background uncertainties for each kernel candidate? 

286 This comes from the inverse of the covariance matrix. 

287 Warning: regularization can cause problems for this step.""", 

288 default=False, 

289 ) 

290 useBicForKernelBasis = pexConfig.Field( 

291 dtype=bool, 

292 doc="""Use Bayesian Information Criterion to select the number of bases going into the kernel""", 

293 default=False, 

294 ) 

295 

296 

297class PsfMatchConfigAL(PsfMatchConfig): 

298 """The parameters specific to the "Alard-Lupton" (sum-of-Gaussian) Psf-matching basis""" 

299 

300 def setDefaults(self): 

301 PsfMatchConfig.setDefaults(self) 

302 self.kernelBasisSet = "alard-lupton" 

303 self.maxConditionNumber = 5.0e7 

304 

305 alardNGauss = pexConfig.Field( 

306 dtype=int, 

307 doc="Number of base Gaussians in alard-lupton kernel basis function generation.", 

308 default=3, 

309 check=lambda x: x >= 1 

310 ) 

311 alardDegGauss = pexConfig.ListField( 

312 dtype=int, 

313 doc="Polynomial order of spatial modification of base Gaussians. " 

314 "List length must be `alardNGauss`.", 

315 default=(4, 2, 2), 

316 ) 

317 alardSigGauss = pexConfig.ListField( 

318 dtype=float, 

319 doc="Default sigma values in pixels of base Gaussians. " 

320 "List length must be `alardNGauss`." 

321 "Only used if the template and science image PSFs have equal size.", 

322 default=(0.7, 1.5, 3.0), 

323 ) 

324 alardGaussBeta = pexConfig.Field( 

325 dtype=float, 

326 doc="Used if `scaleByFwhm==True`, scaling multiplier of base " 

327 "Gaussian sigmas for automated sigma determination", 

328 default=2.0, 

329 check=lambda x: x >= 0.0, 

330 ) 

331 alardMinSig = pexConfig.Field( 

332 dtype=float, 

333 doc="Used if `scaleByFwhm==True`, minimum sigma (pixels) for base Gaussians", 

334 default=0.7, 

335 check=lambda x: x >= 0.25 

336 ) 

337 alardDegGaussDeconv = pexConfig.Field( 

338 dtype=int, 

339 doc="Used if `scaleByFwhm==True`, degree of spatial modification of ALL base Gaussians " 

340 "in AL basis during deconvolution", 

341 default=3, 

342 check=lambda x: x >= 1, 

343 deprecated="No longer used. Will be removed after v29" 

344 ) 

345 alardMinSigDeconv = pexConfig.Field( 

346 dtype=float, 

347 doc="Used if `scaleByFwhm==True`, minimum sigma (pixels) for base Gaussians during deconvolution; " 

348 "make smaller than `alardMinSig` as this is only indirectly used", 

349 default=0.4, 

350 check=lambda x: x >= 0.25 

351 ) 

352 alardNGaussDeconv = pexConfig.Field( 

353 dtype=int, 

354 doc="Used if `scaleByFwhm==True`, number of base Gaussians in AL basis during deconvolution", 

355 default=3, 

356 check=lambda x: x >= 1, 

357 deprecated="No longer used. Will be removed after v29" 

358 ) 

359 

360 

361class PsfMatchConfigDF(PsfMatchConfig): 

362 """The parameters specific to the delta-function (one basis per-pixel) Psf-matching basis""" 

363 

364 def setDefaults(self): 

365 PsfMatchConfig.setDefaults(self) 

366 self.kernelBasisSet = "delta-function" 

367 self.maxConditionNumber = 5.0e6 

368 self.usePcaForSpatialKernel = True 

369 self.subtractMeanForPca = True 

370 self.useBicForKernelBasis = False 

371 

372 useRegularization = pexConfig.Field( 

373 dtype=bool, 

374 doc="Use regularization to smooth the delta function kernels", 

375 default=True, 

376 ) 

377 regularizationType = pexConfig.ChoiceField( 

378 dtype=str, 

379 doc="Type of regularization.", 

380 default="centralDifference", 

381 allowed={ 

382 "centralDifference": "Penalize second derivative using 2-D stencil of central finite difference", 

383 "forwardDifference": "Penalize first, second, third derivatives using forward finite differeces" 

384 } 

385 ) 

386 centralRegularizationStencil = pexConfig.ChoiceField( 

387 dtype=int, 

388 doc="Type of stencil to approximate central derivative (for centralDifference only)", 

389 default=9, 

390 allowed={ 

391 5: "5-point stencil including only adjacent-in-x,y elements", 

392 9: "9-point stencil including diagonal elements" 

393 } 

394 ) 

395 forwardRegularizationOrders = pexConfig.ListField( 

396 dtype=int, 

397 doc="Array showing which order derivatives to penalize (for forwardDifference only)", 

398 default=(1, 2), 

399 itemCheck=lambda x: (x > 0) and (x < 4) 

400 ) 

401 regularizationBorderPenalty = pexConfig.Field( 

402 dtype=float, 

403 doc="Value of the penalty for kernel border pixels", 

404 default=3.0, 

405 check=lambda x: x >= 0.0 

406 ) 

407 lambdaType = pexConfig.ChoiceField( 

408 dtype=str, 

409 doc="How to choose the value of the regularization strength", 

410 default="absolute", 

411 allowed={ 

412 "absolute": "Use lambdaValue as the value of regularization strength", 

413 "relative": "Use lambdaValue as fraction of the default regularization strength (N.R. 18.5.8)", 

414 "minimizeBiasedRisk": "Minimize biased risk estimate", 

415 "minimizeUnbiasedRisk": "Minimize unbiased risk estimate", 

416 } 

417 ) 

418 lambdaValue = pexConfig.Field( 

419 dtype=float, 

420 doc="Value used for absolute determinations of regularization strength", 

421 default=0.2, 

422 ) 

423 lambdaScaling = pexConfig.Field( 

424 dtype=float, 

425 doc="Fraction of the default lambda strength (N.R. 18.5.8) to use. 1e-4 or 1e-5", 

426 default=1e-4, 

427 ) 

428 lambdaStepType = pexConfig.ChoiceField( 

429 dtype=str, 

430 doc="""If a scan through lambda is needed (minimizeBiasedRisk, minimizeUnbiasedRisk), 

431 use log or linear steps""", 

432 default="log", 

433 allowed={ 

434 "log": "Step in log intervals; e.g. lambdaMin, lambdaMax, lambdaStep = -1.0, 2.0, 0.1", 

435 "linear": "Step in linear intervals; e.g. lambdaMin, lambdaMax, lambdaStep = 0.1, 100, 0.1", 

436 } 

437 ) 

438 lambdaMin = pexConfig.Field( 

439 dtype=float, 

440 doc="""If scan through lambda needed (minimizeBiasedRisk, minimizeUnbiasedRisk), 

441 start at this value. If lambdaStepType = log:linear, suggest -1:0.1""", 

442 default=-1.0, 

443 ) 

444 lambdaMax = pexConfig.Field( 

445 dtype=float, 

446 doc="""If scan through lambda needed (minimizeBiasedRisk, minimizeUnbiasedRisk), 

447 stop at this value. If lambdaStepType = log:linear, suggest 2:100""", 

448 default=2.0, 

449 ) 

450 lambdaStep = pexConfig.Field( 

451 dtype=float, 

452 doc="""If scan through lambda needed (minimizeBiasedRisk, minimizeUnbiasedRisk), 

453 step in these increments. If lambdaStepType = log:linear, suggest 0.1:0.1""", 

454 default=0.1, 

455 ) 

456 

457 

458class PsfMatchTask(pipeBase.Task, abc.ABC): 

459 """Base class for Psf Matching; should not be called directly 

460 

461 Notes 

462 ----- 

463 PsfMatchTask is a base class that implements the core functionality for matching the 

464 Psfs of two images using a spatially varying Psf-matching `lsst.afw.math.LinearCombinationKernel`. 

465 The Task requires the user to provide an instance of an `lsst.afw.math.SpatialCellSet`, 

466 filled with `lsst.ip.diffim.KernelCandidate` instances, and a list of `lsst.afw.math.Kernels` 

467 of basis shapes that will be used for the decomposition. If requested, the Task 

468 also performs background matching and returns the differential background model as an 

469 `lsst.afw.math.Kernel.SpatialFunction`. 

470 

471 The initialization sets the Psf-matching kernel configuration using the 

472 value of self.config.kernel.active. If the kernel is requested with 

473 regularization to moderate the bias/variance tradeoff, currently only used 

474 when a delta function kernel basis is provided, it creates a 

475 regularization matrix stored as member variable self.hMat. 

476 

477 **Invoking the Task** 

478 

479 As a base class, this Task is not directly invoked. However, ``run()`` methods that are 

480 implemented on derived classes will make use of the core ``_solve()`` functionality, 

481 which defines a sequence of `lsst.afw.math.CandidateVisitor` classes that iterate 

482 through the KernelCandidates, first building up a per-candidate solution and then 

483 building up a spatial model from the ensemble of candidates. Sigma clipping is 

484 performed using the mean and standard deviation of all kernel sums (to reject 

485 variable objects), on the per-candidate substamp diffim residuals 

486 (to indicate a bad choice of kernel basis shapes for that particular object), 

487 and on the substamp diffim residuals using the spatial kernel fit (to indicate a bad 

488 choice of spatial kernel order, or poor constraints on the spatial model). The 

489 ``_diagnostic()`` method logs information on the quality of the spatial fit, and also 

490 modifies the Task metadata. 

491 

492 .. list-table:: Quantities set in Metadata 

493 :header-rows: 1 

494 

495 * - Parameter 

496 - Description 

497 * - ``spatialConditionNum`` 

498 - Condition number of the spatial kernel fit 

499 * - ``spatialKernelSum`` 

500 - Kernel sum (10^{-0.4 * ``Delta``; zeropoint}) of the spatial Psf-matching kernel 

501 * - ``ALBasisNGauss`` 

502 - If using sum-of-Gaussian basis, the number of gaussians used 

503 * - ``ALBasisDegGauss`` 

504 - If using sum-of-Gaussian basis, the deg of spatial variation of the Gaussians 

505 * - ``ALBasisSigGauss`` 

506 - If using sum-of-Gaussian basis, the widths (sigma) of the Gaussians 

507 * - ``ALKernelSize`` 

508 - If using sum-of-Gaussian basis, the kernel size 

509 * - ``NFalsePositivesTotal`` 

510 - Total number of diaSources 

511 * - ``NFalsePositivesRefAssociated`` 

512 - Number of diaSources that associate with the reference catalog 

513 * - ``NFalsePositivesRefAssociated`` 

514 - Number of diaSources that associate with the source catalog 

515 * - ``NFalsePositivesUnassociated`` 

516 - Number of diaSources that are orphans 

517 * - ``metric_MEAN`` 

518 - Mean value of substamp diffim quality metrics across all KernelCandidates, 

519 for both the per-candidate (LOCAL) and SPATIAL residuals 

520 * - ``metric_MEDIAN`` 

521 - Median value of substamp diffim quality metrics across all KernelCandidates, 

522 for both the per-candidate (LOCAL) and SPATIAL residuals 

523 * - ``metric_STDEV`` 

524 - Standard deviation of substamp diffim quality metrics across all KernelCandidates, 

525 for both the per-candidate (LOCAL) and SPATIAL residuals 

526 

527 **Debug variables** 

528 

529 The ``pipetask`` command line interface supports a 

530 flag --debug to import @b debug.py from your PYTHONPATH. The relevant contents of debug.py 

531 for this Task include: 

532 

533 .. code-block:: py 

534 

535 import sys 

536 import lsstDebug 

537 def DebugInfo(name): 

538 di = lsstDebug.getInfo(name) 

539 if name == "lsst.ip.diffim.psfMatch": 

540 # enable debug output 

541 di.display = True 

542 # display mask transparency 

543 di.maskTransparency = 80 

544 # show all the candidates and residuals 

545 di.displayCandidates = True 

546 # show kernel basis functions 

547 di.displayKernelBasis = False 

548 # show kernel realized across the image 

549 di.displayKernelMosaic = True 

550 # show coefficients of spatial model 

551 di.plotKernelSpatialModel = False 

552 # show fixed and spatial coefficients and coefficient histograms 

553 di.plotKernelCoefficients = True 

554 # show the bad candidates (red) along with good (green) 

555 di.showBadCandidates = True 

556 return di 

557 lsstDebug.Info = DebugInfo 

558 lsstDebug.frame = 1 

559 

560 Note that if you want additional logging info, you may add to your scripts: 

561 

562 .. code-block:: py 

563 

564 import lsst.utils.logging as logUtils 

565 logUtils.trace_set_at("lsst.ip.diffim", 4) 

566 """ 

567 ConfigClass = PsfMatchConfig 

568 _DefaultName = "psfMatch" 

569 

570 def __init__(self, *args, **kwargs): 

571 pipeBase.Task.__init__(self, *args, **kwargs) 

572 self.kConfig = self.config.kernel.active 

573 

574 if 'useRegularization' in self.kConfig: 

575 self.useRegularization = self.kConfig.useRegularization 

576 else: 

577 self.useRegularization = False 

578 

579 if self.useRegularization: 

580 self.hMat = diffimLib.makeRegularizationMatrix(pexConfig.makePropertySet(self.kConfig)) 

581 

582 def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg): 

583 """Provide logging diagnostics on quality of spatial kernel fit 

584 

585 Parameters 

586 ---------- 

587 kernelCellSet : `lsst.afw.math.SpatialCellSet` 

588 Cellset that contains the KernelCandidates used in the fitting 

589 spatialSolution : `lsst.ip.diffim.SpatialKernelSolution` 

590 KernelSolution of best-fit 

591 spatialKernel : `lsst.afw.math.LinearCombinationKernel` 

592 Best-fit spatial Kernel model 

593 spatialBg : `lsst.afw.math.Function2D` 

594 Best-fit spatial background model 

595 """ 

596 # What is the final kernel sum 

597 kImage = afwImage.ImageD(spatialKernel.getDimensions()) 

598 # Evaluate the kernel at the cell-set center, not at the world origin. 

599 xcen, ycen = kernelCellSet.getBBox().getCenter() 

600 kSum = spatialKernel.computeImage(kImage, doNormalize=False, x=xcen, y=ycen) 

601 self.log.info("Final spatial kernel sum %.3f", kSum) 

602 

603 # Look at how well conditioned the matrix is 

604 conditionNum = spatialSolution.getConditionNumber( 

605 getattr(diffimLib.KernelSolution, self.kConfig.conditionNumberType)) 

606 self.log.info("Spatial model condition number %.3e", conditionNum) 

607 

608 if conditionNum < 0.0: 

609 self.log.warning("Condition number is negative (%.3e)", conditionNum) 

610 if conditionNum > self.kConfig.maxSpatialConditionNumber: 

611 self.log.warning("Spatial solution exceeds max condition number (%.3e > %.3e)", 

612 conditionNum, self.kConfig.maxSpatialConditionNumber) 

613 

614 self.metadata["spatialConditionNum"] = conditionNum 

615 self.metadata["spatialKernelSum"] = kSum 

616 

617 # Look at how well the solution is constrained 

618 nBasisKernels = spatialKernel.getNBasisKernels() 

619 nKernelTerms = spatialKernel.getNSpatialParameters() 

620 if nKernelTerms == 0: # order 0 

621 nKernelTerms = 1 

622 

623 # Not fit for 

624 nBgTerms = spatialBg.getNParameters() 

625 if nBgTerms == 1: 

626 if spatialBg.getParameters()[0] == 0.0: 

627 nBgTerms = 0 

628 

629 nGood = 0 

630 nBad = 0 

631 nTot = 0 

632 for cell in kernelCellSet.getCellList(): 

633 for cand in cell.begin(False): # False = include bad candidates 

634 nTot += 1 

635 if cand.getStatus() == afwMath.SpatialCellCandidate.GOOD: 

636 nGood += 1 

637 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD: 

638 nBad += 1 

639 

640 self.log.info("Doing stats of kernel candidates used in the spatial fit.") 

641 

642 # Counting statistics 

643 if nBad > 2*nGood: 

644 self.log.warning("Many more candidates rejected than accepted; %d total, %d rejected, %d used", 

645 nTot, nBad, nGood) 

646 else: 

647 self.log.info("%d candidates total, %d rejected, %d used", nTot, nBad, nGood) 

648 

649 # Some judgements on the quality of the spatial models 

650 if nGood < nKernelTerms: 

651 self.log.warning("Spatial kernel model underconstrained; %d candidates, %d terms, %d bases", 

652 nGood, nKernelTerms, nBasisKernels) 

653 self.log.warning("Consider lowering the spatial order") 

654 elif nGood <= 2*nKernelTerms: 

655 self.log.warning("Spatial kernel model poorly constrained; %d candidates, %d terms, %d bases", 

656 nGood, nKernelTerms, nBasisKernels) 

657 self.log.warning("Consider lowering the spatial order") 

658 else: 

659 self.log.info("Spatial kernel model well constrained; %d candidates, %d terms, %d bases", 

660 nGood, nKernelTerms, nBasisKernels) 

661 

662 if nGood < nBgTerms: 

663 self.log.warning("Spatial background model underconstrained; %d candidates, %d terms", 

664 nGood, nBgTerms) 

665 self.log.warning("Consider lowering the spatial order") 

666 elif nGood <= 2*nBgTerms: 

667 self.log.warning("Spatial background model poorly constrained; %d candidates, %d terms", 

668 nGood, nBgTerms) 

669 self.log.warning("Consider lowering the spatial order") 

670 else: 

671 self.log.info("Spatial background model appears well constrained; %d candidates, %d terms", 

672 nGood, nBgTerms) 

673 

674 def _displayDebug(self, kernelCellSet, spatialKernel, spatialBackground): 

675 """Provide visualization of the inputs and ouputs to the Psf-matching code 

676 

677 Parameters 

678 ---------- 

679 kernelCellSet : `lsst.afw.math.SpatialCellSet` 

680 The SpatialCellSet used in determining the matching kernel and background 

681 spatialKernel : `lsst.afw.math.LinearCombinationKernel` 

682 Spatially varying Psf-matching kernel 

683 spatialBackground : `lsst.afw.math.Function2D` 

684 Spatially varying background-matching function 

685 """ 

686 import lsstDebug 

687 displayCandidates = lsstDebug.Info(__name__).displayCandidates 

688 displayKernelBasis = lsstDebug.Info(__name__).displayKernelBasis 

689 displayKernelMosaic = lsstDebug.Info(__name__).displayKernelMosaic 

690 plotKernelSpatialModel = lsstDebug.Info(__name__).plotKernelSpatialModel 

691 plotKernelCoefficients = lsstDebug.Info(__name__).plotKernelCoefficients 

692 showBadCandidates = lsstDebug.Info(__name__).showBadCandidates 

693 maskTransparency = lsstDebug.Info(__name__).maskTransparency 

694 if not maskTransparency: 

695 maskTransparency = 0 

696 afwDisplay.setDefaultMaskTransparency(maskTransparency) 

697 

698 if displayCandidates: 

699 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground, 

700 frame=lsstDebug.frame, 

701 showBadCandidates=showBadCandidates) 

702 lsstDebug.frame += 1 

703 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground, 

704 frame=lsstDebug.frame, 

705 showBadCandidates=showBadCandidates, 

706 kernels=True) 

707 lsstDebug.frame += 1 

708 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground, 

709 frame=lsstDebug.frame, 

710 showBadCandidates=showBadCandidates, 

711 resids=True) 

712 lsstDebug.frame += 1 

713 

714 if displayKernelBasis: 

715 diutils.showKernelBasis(spatialKernel, frame=lsstDebug.frame) 

716 lsstDebug.frame += 1 

717 

718 if displayKernelMosaic: 

719 diutils.showKernelMosaic(kernelCellSet.getBBox(), spatialKernel, frame=lsstDebug.frame) 

720 lsstDebug.frame += 1 

721 

722 if plotKernelSpatialModel: 

723 diutils.plotKernelSpatialModel(spatialKernel, kernelCellSet, showBadCandidates=showBadCandidates) 

724 

725 if plotKernelCoefficients: 

726 diutils.plotKernelCoefficients(spatialKernel, kernelCellSet) 

727 

728 def _createPcaBasis(self, kernelCellSet, nStarPerCell, ps): 

729 """Create Principal Component basis 

730 

731 If a principal component analysis is requested, typically when using a delta function basis, 

732 perform the PCA here and return a new basis list containing the new principal components. 

733 

734 Parameters 

735 ---------- 

736 kernelCellSet : `lsst.afw.math.SpatialCellSet` 

737 a SpatialCellSet containing KernelCandidates, from which components are derived 

738 nStarPerCell : `int` 

739 the number of stars per cell to visit when doing the PCA 

740 ps : `lsst.daf.base.PropertySet` 

741 input property set controlling the single kernel visitor 

742 

743 Returns 

744 ------- 

745 nRejectedPca : `int` 

746 number of KernelCandidates rejected during PCA loop 

747 spatialBasisList : `list` of `lsst.afw.math.kernel.FixedKernel` 

748 basis list containing the principal shapes as Kernels 

749 

750 Raises 

751 ------ 

752 RuntimeError 

753 If the Eigenvalues sum to zero. 

754 """ 

755 nComponents = self.kConfig.numPrincipalComponents 

756 imagePca = diffimLib.KernelPcaD() 

757 importStarVisitor = diffimLib.KernelPcaVisitorF(imagePca) 

758 kernelCellSet.visitCandidates(importStarVisitor, nStarPerCell) 

759 if self.kConfig.subtractMeanForPca: 

760 importStarVisitor.subtractMean() 

761 imagePca.analyze() 

762 

763 eigenValues = imagePca.getEigenValues() 

764 pcaBasisList = importStarVisitor.getEigenKernels() 

765 

766 eSum = np.sum(eigenValues) 

767 if eSum == 0.0: 

768 raise RuntimeError("Eigenvalues sum to zero") 

769 trace_logger = getTraceLogger(self.log.getChild("_solve"), 5) 

770 for j in range(len(eigenValues)): 

771 trace_logger.debug("Eigenvalue %d : %f (%f)", j, eigenValues[j], eigenValues[j]/eSum) 

772 

773 nToUse = min(nComponents, len(eigenValues)) 

774 trimBasisList = [] 

775 for j in range(nToUse): 

776 # Check for NaNs? 

777 kimage = afwImage.ImageD(pcaBasisList[j].getDimensions()) 

778 pcaBasisList[j].computeImage(kimage, doNormalize=False) 

779 if not (True in np.isnan(kimage.array)): 

780 trimBasisList.append(pcaBasisList[j]) 

781 

782 # Put all the power in the first kernel, which will not vary spatially 

783 spatialBasisList = diffimLib.renormalizeKernelList(trimBasisList) 

784 

785 # New Kernel visitor for this new basis list (no regularization explicitly) 

786 singlekvPca = diffimLib.BuildSingleKernelVisitorF(spatialBasisList, ps) 

787 singlekvPca.setSkipBuilt(False) 

788 kernelCellSet.visitCandidates(singlekvPca, nStarPerCell) 

789 singlekvPca.setSkipBuilt(True) 

790 nRejectedPca = singlekvPca.getNRejected() 

791 

792 return nRejectedPca, spatialBasisList 

793 

794 @abc.abstractmethod 

795 def _buildCellSet(self, *args): 

796 """Fill a SpatialCellSet with KernelCandidates for the Psf-matching process; 

797 override in derived classes""" 

798 return 

799 

800 @timeMethod 

801 def _solve(self, kernelCellSet, basisList): 

802 """Solve for the PSF matching kernel 

803 

804 Parameters 

805 ---------- 

806 kernelCellSet : `lsst.afw.math.SpatialCellSet` 

807 a SpatialCellSet to use in determining the matching kernel 

808 (typically as provided by _buildCellSet) 

809 basisList : `list` of `lsst.afw.math.kernel.FixedKernel` 

810 list of Kernels to be used in the decomposition of the spatially varying kernel 

811 (typically as provided by makeKernelBasisList) 

812 

813 Returns 

814 ------- 

815 psfMatchingKernel : `lsst.afw.math.LinearCombinationKernel` 

816 Spatially varying Psf-matching kernel 

817 backgroundModel : `lsst.afw.math.Function2D` 

818 Spatially varying background-matching function 

819 

820 Raises 

821 ------ 

822 NoKernelCandidatesError : 

823 If there are no useable kernel candidates. 

824 """ 

825 

826 import lsstDebug 

827 display = lsstDebug.Info(__name__).display 

828 

829 maxSpatialIterations = self.kConfig.maxSpatialIterations 

830 nStarPerCell = self.kConfig.nStarPerCell 

831 usePcaForSpatialKernel = self.kConfig.usePcaForSpatialKernel 

832 

833 # Visitor for the single kernel fit 

834 ps = pexConfig.makePropertySet(self.kConfig) 

835 if self.useRegularization: 

836 singlekv = diffimLib.BuildSingleKernelVisitorF(basisList, ps, self.hMat) 

837 else: 

838 singlekv = diffimLib.BuildSingleKernelVisitorF(basisList, ps) 

839 

840 # Visitor for the kernel sum rejection 

841 ksv = diffimLib.KernelSumVisitorF(ps) 

842 

843 # Main loop 

844 trace_loggers = [getTraceLogger(self.log.getChild("_solve"), i) for i in range(5)] 

845 t0 = time.time() 

846 totalIterations = 0 

847 thisIteration = 0 

848 while (thisIteration < maxSpatialIterations): 

849 

850 # Make sure there are no uninitialized candidates as active occupants of Cell 

851 nRejectedSkf = -1 

852 while (nRejectedSkf != 0): 

853 trace_loggers[1].debug("Building single kernels...") 

854 kernelCellSet.visitCandidates(singlekv, nStarPerCell, ignoreExceptions=True) 

855 nRejectedSkf = singlekv.getNRejected() 

856 trace_loggers[1].debug( 

857 "Iteration %d, rejected %d candidates due to initial kernel fit", 

858 thisIteration, nRejectedSkf 

859 ) 

860 

861 # Reject outliers in kernel sum 

862 ksv.resetKernelSum() 

863 ksv.setMode(diffimLib.KernelSumVisitorF.AGGREGATE) 

864 kernelCellSet.visitCandidates(ksv, nStarPerCell, ignoreExceptions=True) 

865 

866 allCellsEmpty = True 

867 for cell in kernelCellSet.getCellList(): 

868 if not cell.empty(): 

869 allCellsEmpty = False 

870 break 

871 if allCellsEmpty: 

872 raise NoKernelCandidatesError("All spatial cells are empty of candidates") 

873 

874 try: 

875 ksv.processKsumDistribution() 

876 except RuntimeError as e: 

877 raise NoKernelCandidatesError("Unable to calculate the kernel sum") from e 

878 

879 ksv.setMode(diffimLib.KernelSumVisitorF.REJECT) 

880 kernelCellSet.visitCandidates(ksv, nStarPerCell, ignoreExceptions=True) 

881 

882 nRejectedKsum = ksv.getNRejected() 

883 trace_loggers[1].debug( 

884 "Iteration %d, rejected %d candidates due to kernel sum", 

885 thisIteration, nRejectedKsum 

886 ) 

887 

888 # Do we jump back to the top without incrementing thisIteration? 

889 if nRejectedKsum > 0: 

890 totalIterations += 1 

891 continue 

892 

893 # At this stage we can either apply the spatial fit to 

894 # the kernels, or we run a PCA, use these as a *new* 

895 # basis set with lower dimensionality, and then apply 

896 # the spatial fit to these kernels 

897 

898 if (usePcaForSpatialKernel): 

899 trace_loggers[0].debug("Building Pca basis") 

900 

901 nRejectedPca, spatialBasisList = self._createPcaBasis(kernelCellSet, nStarPerCell, ps) 

902 trace_loggers[1].debug( 

903 "Iteration %d, rejected %d candidates due to Pca kernel fit", 

904 thisIteration, nRejectedPca 

905 ) 

906 

907 # We don't want to continue on (yet) with the 

908 # spatial modeling, because we have bad objects 

909 # contributing to the Pca basis. We basically 

910 # need to restart from the beginning of this loop, 

911 # since the cell-mates of those objects that were 

912 # rejected need their original Kernels built by 

913 # singleKernelFitter. 

914 

915 # Don't count against thisIteration 

916 if (nRejectedPca > 0): 

917 totalIterations += 1 

918 continue 

919 else: 

920 spatialBasisList = basisList 

921 

922 # We have gotten on to the spatial modeling part 

923 regionBBox = kernelCellSet.getBBox() 

924 spatialkv = diffimLib.BuildSpatialKernelVisitorF(spatialBasisList, regionBBox, ps) 

925 kernelCellSet.visitCandidates(spatialkv, nStarPerCell) 

926 spatialkv.solveLinearEquation() 

927 trace_loggers[2].debug("Spatial kernel built with %d candidates", spatialkv.getNCandidates()) 

928 spatialKernel, spatialBackground = spatialkv.getSolutionPair() 

929 

930 # Check the quality of the spatial fit (look at residuals) 

931 assesskv = diffimLib.AssessSpatialKernelVisitorF(spatialKernel, spatialBackground, ps)