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 kSum = spatialKernel.computeImage(kImage, False) 

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

600 

601 # Look at how well conditioned the matrix is 

602 conditionNum = spatialSolution.getConditionNumber( 

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

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

605 

606 if conditionNum < 0.0: 

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

608 if conditionNum > self.kConfig.maxSpatialConditionNumber: 

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

610 conditionNum, self.kConfig.maxSpatialConditionNumber) 

611 

612 self.metadata["spatialConditionNum"] = conditionNum 

613 self.metadata["spatialKernelSum"] = kSum 

614 

615 # Look at how well the solution is constrained 

616 nBasisKernels = spatialKernel.getNBasisKernels() 

617 nKernelTerms = spatialKernel.getNSpatialParameters() 

618 if nKernelTerms == 0: # order 0 

619 nKernelTerms = 1 

620 

621 # Not fit for 

622 nBgTerms = spatialBg.getNParameters() 

623 if nBgTerms == 1: 

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

625 nBgTerms = 0 

626 

627 nGood = 0 

628 nBad = 0 

629 nTot = 0 

630 for cell in kernelCellSet.getCellList(): 

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

632 nTot += 1 

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

634 nGood += 1 

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

636 nBad += 1 

637 

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

639 

640 # Counting statistics 

641 if nBad > 2*nGood: 

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

643 nTot, nBad, nGood) 

644 else: 

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

646 

647 # Some judgements on the quality of the spatial models 

648 if nGood < nKernelTerms: 

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

650 nGood, nKernelTerms, nBasisKernels) 

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

652 elif nGood <= 2*nKernelTerms: 

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

654 nGood, nKernelTerms, nBasisKernels) 

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

656 else: 

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

658 nGood, nKernelTerms, nBasisKernels) 

659 

660 if nGood < nBgTerms: 

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

662 nGood, nBgTerms) 

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

664 elif nGood <= 2*nBgTerms: 

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

666 nGood, nBgTerms) 

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

668 else: 

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

670 nGood, nBgTerms) 

671 

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

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

674 

675 Parameters 

676 ---------- 

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

678 The SpatialCellSet used in determining the matching kernel and background 

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

680 Spatially varying Psf-matching kernel 

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

682 Spatially varying background-matching function 

683 """ 

684 import lsstDebug 

685 displayCandidates = lsstDebug.Info(__name__).displayCandidates 

686 displayKernelBasis = lsstDebug.Info(__name__).displayKernelBasis 

687 displayKernelMosaic = lsstDebug.Info(__name__).displayKernelMosaic 

688 plotKernelSpatialModel = lsstDebug.Info(__name__).plotKernelSpatialModel 

689 plotKernelCoefficients = lsstDebug.Info(__name__).plotKernelCoefficients 

690 showBadCandidates = lsstDebug.Info(__name__).showBadCandidates 

691 maskTransparency = lsstDebug.Info(__name__).maskTransparency 

692 if not maskTransparency: 

693 maskTransparency = 0 

694 afwDisplay.setDefaultMaskTransparency(maskTransparency) 

695 

696 if displayCandidates: 

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

698 frame=lsstDebug.frame, 

699 showBadCandidates=showBadCandidates) 

700 lsstDebug.frame += 1 

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

702 frame=lsstDebug.frame, 

703 showBadCandidates=showBadCandidates, 

704 kernels=True) 

705 lsstDebug.frame += 1 

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

707 frame=lsstDebug.frame, 

708 showBadCandidates=showBadCandidates, 

709 resids=True) 

710 lsstDebug.frame += 1 

711 

712 if displayKernelBasis: 

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

714 lsstDebug.frame += 1 

715 

716 if displayKernelMosaic: 

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

718 lsstDebug.frame += 1 

719 

720 if plotKernelSpatialModel: 

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

722 

723 if plotKernelCoefficients: 

724 diutils.plotKernelCoefficients(spatialKernel, kernelCellSet) 

725 

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

727 """Create Principal Component basis 

728 

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

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

731 

732 Parameters 

733 ---------- 

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

735 a SpatialCellSet containing KernelCandidates, from which components are derived 

736 nStarPerCell : `int` 

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

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

739 input property set controlling the single kernel visitor 

740 

741 Returns 

742 ------- 

743 nRejectedPca : `int` 

744 number of KernelCandidates rejected during PCA loop 

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

746 basis list containing the principal shapes as Kernels 

747 

748 Raises 

749 ------ 

750 RuntimeError 

751 If the Eigenvalues sum to zero. 

752 """ 

753 nComponents = self.kConfig.numPrincipalComponents 

754 imagePca = diffimLib.KernelPcaD() 

755 importStarVisitor = diffimLib.KernelPcaVisitorF(imagePca) 

756 kernelCellSet.visitCandidates(importStarVisitor, nStarPerCell) 

757 if self.kConfig.subtractMeanForPca: 

758 importStarVisitor.subtractMean() 

759 imagePca.analyze() 

760 

761 eigenValues = imagePca.getEigenValues() 

762 pcaBasisList = importStarVisitor.getEigenKernels() 

763 

764 eSum = np.sum(eigenValues) 

765 if eSum == 0.0: 

766 raise RuntimeError("Eigenvalues sum to zero") 

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

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

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

770 

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

772 trimBasisList = [] 

773 for j in range(nToUse): 

774 # Check for NaNs? 

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

776 pcaBasisList[j].computeImage(kimage, False) 

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

778 trimBasisList.append(pcaBasisList[j]) 

779 

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

781 spatialBasisList = diffimLib.renormalizeKernelList(trimBasisList) 

782 

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

784 singlekvPca = diffimLib.BuildSingleKernelVisitorF(spatialBasisList, ps) 

785 singlekvPca.setSkipBuilt(False) 

786 kernelCellSet.visitCandidates(singlekvPca, nStarPerCell) 

787 singlekvPca.setSkipBuilt(True) 

788 nRejectedPca = singlekvPca.getNRejected() 

789 

790 return nRejectedPca, spatialBasisList 

791 

792 @abc.abstractmethod 

793 def _buildCellSet(self, *args): 

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

795 override in derived classes""" 

796 return 

797 

798 @timeMethod 

799 def _solve(self, kernelCellSet, basisList): 

800 """Solve for the PSF matching kernel 

801 

802 Parameters 

803 ---------- 

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

805 a SpatialCellSet to use in determining the matching kernel 

806 (typically as provided by _buildCellSet) 

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

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

809 (typically as provided by makeKernelBasisList) 

810 

811 Returns 

812 ------- 

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

814 Spatially varying Psf-matching kernel 

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

816 Spatially varying background-matching function 

817 

818 Raises 

819 ------ 

820 NoKernelCandidatesError : 

821 If there are no useable kernel candidates. 

822 """ 

823 

824 import lsstDebug 

825 display = lsstDebug.Info(__name__).display 

826 

827 maxSpatialIterations = self.kConfig.maxSpatialIterations 

828 nStarPerCell = self.kConfig.nStarPerCell 

829 usePcaForSpatialKernel = self.kConfig.usePcaForSpatialKernel 

830 

831 # Visitor for the single kernel fit 

832 ps = pexConfig.makePropertySet(self.kConfig) 

833 if self.useRegularization: 

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

835 else: 

836 singlekv = diffimLib.BuildSingleKernelVisitorF(basisList, ps) 

837 

838 # Visitor for the kernel sum rejection 

839 ksv = diffimLib.KernelSumVisitorF(ps) 

840 

841 # Main loop 

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

843 t0 = time.time() 

844 totalIterations = 0 

845 thisIteration = 0 

846 while (thisIteration < maxSpatialIterations): 

847 

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

849 nRejectedSkf = -1 

850 while (nRejectedSkf != 0): 

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

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

853 nRejectedSkf = singlekv.getNRejected() 

854 trace_loggers[1].debug( 

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

856 thisIteration, nRejectedSkf 

857 ) 

858 

859 # Reject outliers in kernel sum 

860 ksv.resetKernelSum() 

861 ksv.setMode(diffimLib.KernelSumVisitorF.AGGREGATE) 

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

863 

864 allCellsEmpty = True 

865 for cell in kernelCellSet.getCellList(): 

866 if not cell.empty(): 

867 allCellsEmpty = False 

868 break 

869 if allCellsEmpty: 

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

871 

872 try: 

873 ksv.processKsumDistribution() 

874 except RuntimeError as e: 

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

876 

877 ksv.setMode(diffimLib.KernelSumVisitorF.REJECT) 

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

879 

880 nRejectedKsum = ksv.getNRejected() 

881 trace_loggers[1].debug( 

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

883 thisIteration, nRejectedKsum 

884 ) 

885 

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

887 if nRejectedKsum > 0: 

888 totalIterations += 1 

889 continue 

890 

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

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

893 # basis set with lower dimensionality, and then apply 

894 # the spatial fit to these kernels 

895 

896 if (usePcaForSpatialKernel): 

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

898 

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

900 trace_loggers[1].debug( 

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

902 thisIteration, nRejectedPca 

903 ) 

904 

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

906 # spatial modeling, because we have bad objects 

907 # contributing to the Pca basis. We basically 

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

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

910 # rejected need their original Kernels built by 

911 # singleKernelFitter. 

912 

913 # Don't count against thisIteration 

914 if (nRejectedPca > 0): 

915 totalIterations += 1 

916 continue 

917 else: 

918 spatialBasisList = basisList 

919 

920 # We have gotten on to the spatial modeling part 

921 regionBBox = kernelCellSet.getBBox() 

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

923 kernelCellSet.visitCandidates(spatialkv, nStarPerCell) 

924 spatialkv.solveLinearEquation() 

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

926 spatialKernel, spatialBackground = spatialkv.getSolutionPair() 

927 

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

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

930 kernelCellSet.visitCandidates(assesskv, nStarPerCell) 

931 nRejectedSpatial = assesskv.getNRejected()