LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
psfMatch.py
Go to the documentation of this file.
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 def metadata(self) -> dict:
46 return {}
47
48
49class PsfMatchConfig(pexConfig.Config):
50 """Base configuration for Psf-matching
51
52 The base configuration of the Psf-matching kernel, and of the warping, detection,
53 and background modeling subTasks."""
54
55 warpingConfig = pexConfig.ConfigField("Config for warping exposures to a common alignment",
57 afwBackgroundConfig = pexConfig.ConfigField("Controlling the Afw background fitting",
58 SubtractBackgroundConfig)
59
60 useAfwBackground = pexConfig.Field(
61 dtype=bool,
62 doc="Use afw background subtraction instead of ip_diffim",
63 default=False,
64 )
65 fitForBackground = pexConfig.Field(
66 dtype=bool,
67 doc="Include terms (including kernel cross terms) for background in ip_diffim",
68 default=False,
69 )
70 kernelBasisSet = pexConfig.ChoiceField(
71 dtype=str,
72 doc="Type of basis set for PSF matching kernel.",
73 default="alard-lupton",
74 allowed={
75 "alard-lupton": """Alard-Lupton sum-of-gaussians basis set,
76 * The first term has no spatial variation
77 * The kernel sum is conserved
78 * You may want to turn off 'usePcaForSpatialKernel'""",
79 "delta-function": """Delta-function kernel basis set,
80 * You may enable the option useRegularization
81 * You should seriously consider usePcaForSpatialKernel, which will also
82 enable kernel sum conservation for the delta function kernels"""
83 }
84 )
85 kernelSize = pexConfig.Field(
86 dtype=int,
87 doc="""Number of rows/columns in the convolution kernel; should be odd-valued.
88 Modified by kernelSizeFwhmScaling if scaleByFwhm = true""",
89 default=21,
90 )
91 scaleByFwhm = pexConfig.Field(
92 dtype=bool,
93 doc="Scale kernelSize, alardGaussians by input Fwhm",
94 default=True,
95 )
96 kernelSizeFwhmScaling = pexConfig.Field(
97 dtype=float,
98 doc="Multiplier of the largest AL Gaussian basis sigma to get the kernel bbox (pixel) size.",
99 default=6.0,
100 check=lambda x: x >= 1.0
101 )
102 kernelSizeMin = pexConfig.Field(
103 dtype=int,
104 doc="Minimum kernel bbox (pixel) size.",
105 default=21,
106 )
107 kernelSizeMax = pexConfig.Field(
108 dtype=int,
109 doc="Maximum kernel bbox (pixel) size.",
110 default=35,
111 )
112 spatialModelType = pexConfig.ChoiceField(
113 dtype=str,
114 doc="Type of spatial functions for kernel and background",
115 default="chebyshev1",
116 allowed={
117 "chebyshev1": "Chebyshev polynomial of the first kind",
118 "polynomial": "Standard x,y polynomial",
119 }
120 )
121 spatialKernelOrder = pexConfig.Field(
122 dtype=int,
123 doc="Spatial order of convolution kernel variation",
124 default=2,
125 check=lambda x: x >= 0
126 )
127 spatialBgOrder = pexConfig.Field(
128 dtype=int,
129 doc="Spatial order of differential background variation",
130 default=1,
131 check=lambda x: x >= 0
132 )
133 sizeCellX = pexConfig.Field(
134 dtype=int,
135 doc="Size (rows) in pixels of each SpatialCell for spatial modeling",
136 default=128,
137 check=lambda x: x >= 32
138 )
139 sizeCellY = pexConfig.Field(
140 dtype=int,
141 doc="Size (columns) in pixels of each SpatialCell for spatial modeling",
142 default=128,
143 check=lambda x: x >= 32
144 )
145 nStarPerCell = pexConfig.Field(
146 dtype=int,
147 doc="Maximum number of KernelCandidates in each SpatialCell to use in the spatial fitting. "
148 "Set to -1 to use all candidates in each cell.",
149 default=5,
150 )
151 maxSpatialIterations = pexConfig.Field(
152 dtype=int,
153 doc="Maximum number of iterations for rejecting bad KernelCandidates in spatial fitting",
154 default=3,
155 check=lambda x: x >= 1 and x <= 5
156 )
157 usePcaForSpatialKernel = pexConfig.Field(
158 dtype=bool,
159 doc="""Use Pca to reduce the dimensionality of the kernel basis sets.
160 This is particularly useful for delta-function kernels.
161 Functionally, after all Cells have their raw kernels determined, we run
162 a Pca on these Kernels, re-fit the Cells using the eigenKernels and then
163 fit those for spatial variation using the same technique as for Alard-Lupton kernels.
164 If this option is used, the first term will have no spatial variation and the
165 kernel sum will be conserved.""",
166 default=False,
167 )
168 subtractMeanForPca = pexConfig.Field(
169 dtype=bool,
170 doc="Subtract off the mean feature before doing the Pca",
171 default=True,
172 )
173 numPrincipalComponents = pexConfig.Field(
174 dtype=int,
175 doc="""Number of principal components to use for Pca basis, including the
176 mean kernel if requested.""",
177 default=5,
178 check=lambda x: x >= 3
179 )
180 singleKernelClipping = pexConfig.Field(
181 dtype=bool,
182 doc="Do sigma clipping on each raw kernel candidate",
183 default=True,
184 )
185 kernelSumClipping = pexConfig.Field(
186 dtype=bool,
187 doc="Do sigma clipping on the ensemble of kernel sums",
188 default=True,
189 )
190 spatialKernelClipping = pexConfig.Field(
191 dtype=bool,
192 doc="Do sigma clipping after building the spatial model",
193 default=True,
194 )
195 checkConditionNumber = pexConfig.Field(
196 dtype=bool,
197 doc="""Test for maximum condition number when inverting a kernel matrix.
198 Anything above maxConditionNumber is not used and the candidate is set as BAD.
199 Also used to truncate inverse matrix in estimateBiasedRisk. However,
200 if you are doing any deconvolution you will want to turn this off, or use
201 a large maxConditionNumber""",
202 default=False,
203 )
204 badMaskPlanes = pexConfig.ListField(
205 dtype=str,
206 doc="""Mask planes to ignore when calculating diffim statistics
207 Options: NO_DATA EDGE SAT BAD CR INTRP""",
208 default=("NO_DATA", "EDGE", "SAT")
209 )
210 candidateResidualMeanMax = pexConfig.Field(
211 dtype=float,
212 doc="""Rejects KernelCandidates yielding bad difference image quality.
213 Used by BuildSingleKernelVisitor, AssessSpatialKernelVisitor.
214 Represents average over pixels of (image/sqrt(variance)).""",
215 default=0.25,
216 check=lambda x: x >= 0.0
217 )
218 candidateResidualStdMax = pexConfig.Field(
219 dtype=float,
220 doc="""Rejects KernelCandidates yielding bad difference image quality.
221 Used by BuildSingleKernelVisitor, AssessSpatialKernelVisitor.
222 Represents stddev over pixels of (image/sqrt(variance)).""",
223 default=1.50,
224 check=lambda x: x >= 0.0
225 )
226 useCoreStats = pexConfig.Field(
227 dtype=bool,
228 doc="""Use the core of the footprint for the quality statistics, instead of the entire footprint.
229 WARNING: if there is deconvolution we probably will need to turn this off""",
230 default=False,
231 )
232 candidateCoreRadius = pexConfig.Field(
233 dtype=int,
234 doc="""Radius for calculation of stats in 'core' of KernelCandidate diffim.
235 Total number of pixels used will be (2*radius)**2.
236 This is used both for 'core' diffim quality as well as ranking of
237 KernelCandidates by their total flux in this core""",
238 default=3,
239 check=lambda x: x >= 1
240 )
241 maxKsumSigma = pexConfig.Field(
242 dtype=float,
243 doc="""Maximum allowed sigma for outliers from kernel sum distribution.
244 Used to reject variable objects from the kernel model""",
245 default=3.0,
246 check=lambda x: x >= 0.0
247 )
248 maxConditionNumber = pexConfig.Field(
249 dtype=float,
250 doc="Maximum condition number for a well conditioned matrix",
251 default=5.0e7,
252 check=lambda x: x >= 0.0
253 )
254 conditionNumberType = pexConfig.ChoiceField(
255 dtype=str,
256 doc="Use singular values (SVD) or eigen values (EIGENVALUE) to determine condition number",
257 default="EIGENVALUE",
258 allowed={
259 "SVD": "Use singular values",
260 "EIGENVALUE": "Use eigen values (faster)",
261 }
262 )
263 maxSpatialConditionNumber = pexConfig.Field(
264 dtype=float,
265 doc="Maximum condition number for a well conditioned spatial matrix",
266 default=1.0e10,
267 check=lambda x: x >= 0.0
268 )
269 iterateSingleKernel = pexConfig.Field(
270 dtype=bool,
271 doc="""Remake KernelCandidate using better variance estimate after first pass?
272 Primarily useful when convolving a single-depth image, otherwise not necessary.""",
273 default=False,
274 )
275 constantVarianceWeighting = pexConfig.Field(
276 dtype=bool,
277 doc="""Use constant variance weighting in single kernel fitting?
278 In some cases this is better for bright star residuals.""",
279 default=True,
280 )
281 calculateKernelUncertainty = pexConfig.Field(
282 dtype=bool,
283 doc="""Calculate kernel and background uncertainties for each kernel candidate?
284 This comes from the inverse of the covariance matrix.
285 Warning: regularization can cause problems for this step.""",
286 default=False,
287 )
288 useBicForKernelBasis = pexConfig.Field(
289 dtype=bool,
290 doc="""Use Bayesian Information Criterion to select the number of bases going into the kernel""",
291 default=False,
292 )
293
294
295class PsfMatchConfigAL(PsfMatchConfig):
296 """The parameters specific to the "Alard-Lupton" (sum-of-Gaussian) Psf-matching basis"""
297
298 def setDefaults(self):
299 PsfMatchConfig.setDefaults(self)
300 self.kernelBasisSet = "alard-lupton"
301 self.maxConditionNumber = 5.0e7
302
303 alardNGauss = pexConfig.Field(
304 dtype=int,
305 doc="Number of base Gaussians in alard-lupton kernel basis function generation.",
306 default=3,
307 check=lambda x: x >= 1
308 )
309 alardDegGauss = pexConfig.ListField(
310 dtype=int,
311 doc="Polynomial order of spatial modification of base Gaussians. "
312 "List length must be `alardNGauss`.",
313 default=(4, 2, 2),
314 )
315 alardSigGauss = pexConfig.ListField(
316 dtype=float,
317 doc="Default sigma values in pixels of base Gaussians. "
318 "List length must be `alardNGauss`."
319 "Only used if the template and science image PSFs have equal size.",
320 default=(0.7, 1.5, 3.0),
321 )
322 alardGaussBeta = pexConfig.Field(
323 dtype=float,
324 doc="Used if `scaleByFwhm==True`, scaling multiplier of base "
325 "Gaussian sigmas for automated sigma determination",
326 default=2.0,
327 check=lambda x: x >= 0.0,
328 )
329 alardMinSig = pexConfig.Field(
330 dtype=float,
331 doc="Used if `scaleByFwhm==True`, minimum sigma (pixels) for base Gaussians",
332 default=0.7,
333 check=lambda x: x >= 0.25
334 )
335 alardDegGaussDeconv = pexConfig.Field(
336 dtype=int,
337 doc="Used if `scaleByFwhm==True`, degree of spatial modification of ALL base Gaussians "
338 "in AL basis during deconvolution",
339 default=3,
340 check=lambda x: x >= 1,
341 deprecated="No longer used. Will be removed after v29"
342 )
343 alardMinSigDeconv = pexConfig.Field(
344 dtype=float,
345 doc="Used if `scaleByFwhm==True`, minimum sigma (pixels) for base Gaussians during deconvolution; "
346 "make smaller than `alardMinSig` as this is only indirectly used",
347 default=0.4,
348 check=lambda x: x >= 0.25
349 )
350 alardNGaussDeconv = pexConfig.Field(
351 dtype=int,
352 doc="Used if `scaleByFwhm==True`, number of base Gaussians in AL basis during deconvolution",
353 default=3,
354 check=lambda x: x >= 1,
355 deprecated="No longer used. Will be removed after v29"
356 )
357
358
359class PsfMatchConfigDF(PsfMatchConfig):
360 """The parameters specific to the delta-function (one basis per-pixel) Psf-matching basis"""
361
362 def setDefaults(self):
363 PsfMatchConfig.setDefaults(self)
364 self.kernelBasisSet = "delta-function"
365 self.maxConditionNumber = 5.0e6
366 self.usePcaForSpatialKernel = True
367 self.subtractMeanForPca = True
368 self.useBicForKernelBasis = False
369
370 useRegularization = pexConfig.Field(
371 dtype=bool,
372 doc="Use regularization to smooth the delta function kernels",
373 default=True,
374 )
375 regularizationType = pexConfig.ChoiceField(
376 dtype=str,
377 doc="Type of regularization.",
378 default="centralDifference",
379 allowed={
380 "centralDifference": "Penalize second derivative using 2-D stencil of central finite difference",
381 "forwardDifference": "Penalize first, second, third derivatives using forward finite differeces"
382 }
383 )
384 centralRegularizationStencil = pexConfig.ChoiceField(
385 dtype=int,
386 doc="Type of stencil to approximate central derivative (for centralDifference only)",
387 default=9,
388 allowed={
389 5: "5-point stencil including only adjacent-in-x,y elements",
390 9: "9-point stencil including diagonal elements"
391 }
392 )
393 forwardRegularizationOrders = pexConfig.ListField(
394 dtype=int,
395 doc="Array showing which order derivatives to penalize (for forwardDifference only)",
396 default=(1, 2),
397 itemCheck=lambda x: (x > 0) and (x < 4)
398 )
399 regularizationBorderPenalty = pexConfig.Field(
400 dtype=float,
401 doc="Value of the penalty for kernel border pixels",
402 default=3.0,
403 check=lambda x: x >= 0.0
404 )
405 lambdaType = pexConfig.ChoiceField(
406 dtype=str,
407 doc="How to choose the value of the regularization strength",
408 default="absolute",
409 allowed={
410 "absolute": "Use lambdaValue as the value of regularization strength",
411 "relative": "Use lambdaValue as fraction of the default regularization strength (N.R. 18.5.8)",
412 "minimizeBiasedRisk": "Minimize biased risk estimate",
413 "minimizeUnbiasedRisk": "Minimize unbiased risk estimate",
414 }
415 )
416 lambdaValue = pexConfig.Field(
417 dtype=float,
418 doc="Value used for absolute determinations of regularization strength",
419 default=0.2,
420 )
421 lambdaScaling = pexConfig.Field(
422 dtype=float,
423 doc="Fraction of the default lambda strength (N.R. 18.5.8) to use. 1e-4 or 1e-5",
424 default=1e-4,
425 )
426 lambdaStepType = pexConfig.ChoiceField(
427 dtype=str,
428 doc="""If a scan through lambda is needed (minimizeBiasedRisk, minimizeUnbiasedRisk),
429 use log or linear steps""",
430 default="log",
431 allowed={
432 "log": "Step in log intervals; e.g. lambdaMin, lambdaMax, lambdaStep = -1.0, 2.0, 0.1",
433 "linear": "Step in linear intervals; e.g. lambdaMin, lambdaMax, lambdaStep = 0.1, 100, 0.1",
434 }
435 )
436 lambdaMin = pexConfig.Field(
437 dtype=float,
438 doc="""If scan through lambda needed (minimizeBiasedRisk, minimizeUnbiasedRisk),
439 start at this value. If lambdaStepType = log:linear, suggest -1:0.1""",
440 default=-1.0,
441 )
442 lambdaMax = pexConfig.Field(
443 dtype=float,
444 doc="""If scan through lambda needed (minimizeBiasedRisk, minimizeUnbiasedRisk),
445 stop at this value. If lambdaStepType = log:linear, suggest 2:100""",
446 default=2.0,
447 )
448 lambdaStep = pexConfig.Field(
449 dtype=float,
450 doc="""If scan through lambda needed (minimizeBiasedRisk, minimizeUnbiasedRisk),
451 step in these increments. If lambdaStepType = log:linear, suggest 0.1:0.1""",
452 default=0.1,
453 )
454
455
456class PsfMatchTask(pipeBase.Task, abc.ABC):
457 """Base class for Psf Matching; should not be called directly
458
459 Notes
460 -----
461 PsfMatchTask is a base class that implements the core functionality for matching the
462 Psfs of two images using a spatially varying Psf-matching `lsst.afw.math.LinearCombinationKernel`.
463 The Task requires the user to provide an instance of an `lsst.afw.math.SpatialCellSet`,
464 filled with `lsst.ip.diffim.KernelCandidate` instances, and a list of `lsst.afw.math.Kernels`
465 of basis shapes that will be used for the decomposition. If requested, the Task
466 also performs background matching and returns the differential background model as an
467 `lsst.afw.math.Kernel.SpatialFunction`.
468
469 The initialization sets the Psf-matching kernel configuration using the
470 value of self.config.kernel.active. If the kernel is requested with
471 regularization to moderate the bias/variance tradeoff, currently only used
472 when a delta function kernel basis is provided, it creates a
473 regularization matrix stored as member variable self.hMat.
474
475 **Invoking the Task**
476
477 As a base class, this Task is not directly invoked. However, ``run()`` methods that are
478 implemented on derived classes will make use of the core ``_solve()`` functionality,
479 which defines a sequence of `lsst.afw.math.CandidateVisitor` classes that iterate
480 through the KernelCandidates, first building up a per-candidate solution and then
481 building up a spatial model from the ensemble of candidates. Sigma clipping is
482 performed using the mean and standard deviation of all kernel sums (to reject
483 variable objects), on the per-candidate substamp diffim residuals
484 (to indicate a bad choice of kernel basis shapes for that particular object),
485 and on the substamp diffim residuals using the spatial kernel fit (to indicate a bad
486 choice of spatial kernel order, or poor constraints on the spatial model). The
487 ``_diagnostic()`` method logs information on the quality of the spatial fit, and also
488 modifies the Task metadata.
489
490 .. list-table:: Quantities set in Metadata
491 :header-rows: 1
492
493 * - Parameter
494 - Description
495 * - ``spatialConditionNum``
496 - Condition number of the spatial kernel fit
497 * - ``spatialKernelSum``
498 - Kernel sum (10^{-0.4 * ``Delta``; zeropoint}) of the spatial Psf-matching kernel
499 * - ``ALBasisNGauss``
500 - If using sum-of-Gaussian basis, the number of gaussians used
501 * - ``ALBasisDegGauss``
502 - If using sum-of-Gaussian basis, the deg of spatial variation of the Gaussians
503 * - ``ALBasisSigGauss``
504 - If using sum-of-Gaussian basis, the widths (sigma) of the Gaussians
505 * - ``ALKernelSize``
506 - If using sum-of-Gaussian basis, the kernel size
507 * - ``NFalsePositivesTotal``
508 - Total number of diaSources
509 * - ``NFalsePositivesRefAssociated``
510 - Number of diaSources that associate with the reference catalog
511 * - ``NFalsePositivesRefAssociated``
512 - Number of diaSources that associate with the source catalog
513 * - ``NFalsePositivesUnassociated``
514 - Number of diaSources that are orphans
515 * - ``metric_MEAN``
516 - Mean value of substamp diffim quality metrics across all KernelCandidates,
517 for both the per-candidate (LOCAL) and SPATIAL residuals
518 * - ``metric_MEDIAN``
519 - Median value of substamp diffim quality metrics across all KernelCandidates,
520 for both the per-candidate (LOCAL) and SPATIAL residuals
521 * - ``metric_STDEV``
522 - Standard deviation of substamp diffim quality metrics across all KernelCandidates,
523 for both the per-candidate (LOCAL) and SPATIAL residuals
524
525 **Debug variables**
526
527 The ``pipetask`` command line interface supports a
528 flag --debug to import @b debug.py from your PYTHONPATH. The relevant contents of debug.py
529 for this Task include:
530
531 .. code-block:: py
532
533 import sys
534 import lsstDebug
535 def DebugInfo(name):
536 di = lsstDebug.getInfo(name)
537 if name == "lsst.ip.diffim.psfMatch":
538 # enable debug output
539 di.display = True
540 # display mask transparency
541 di.maskTransparency = 80
542 # show all the candidates and residuals
543 di.displayCandidates = True
544 # show kernel basis functions
545 di.displayKernelBasis = False
546 # show kernel realized across the image
547 di.displayKernelMosaic = True
548 # show coefficients of spatial model
549 di.plotKernelSpatialModel = False
550 # show fixed and spatial coefficients and coefficient histograms
551 di.plotKernelCoefficients = True
552 # show the bad candidates (red) along with good (green)
553 di.showBadCandidates = True
554 return di
555 lsstDebug.Info = DebugInfo
556 lsstDebug.frame = 1
557
558 Note that if you want additional logging info, you may add to your scripts:
559
560 .. code-block:: py
561
562 import lsst.utils.logging as logUtils
563 logUtils.trace_set_at("lsst.ip.diffim", 4)
564 """
565 ConfigClass = PsfMatchConfig
566 _DefaultName = "psfMatch"
567
568 def __init__(self, *args, **kwargs):
569 pipeBase.Task.__init__(self, *args, **kwargs)
570 self.kConfig = self.config.kernel.active
571
572 if 'useRegularization' in self.kConfig:
573 self.useRegularization = self.kConfig.useRegularization
574 else:
575 self.useRegularization = False
576
577 if self.useRegularization:
578 self.hMat = diffimLib.makeRegularizationMatrix(pexConfig.makePropertySet(self.kConfig))
579
580 def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
581 """Provide logging diagnostics on quality of spatial kernel fit
582
583 Parameters
584 ----------
585 kernelCellSet : `lsst.afw.math.SpatialCellSet`
586 Cellset that contains the KernelCandidates used in the fitting
587 spatialSolution : `lsst.ip.diffim.SpatialKernelSolution`
588 KernelSolution of best-fit
589 spatialKernel : `lsst.afw.math.LinearCombinationKernel`
590 Best-fit spatial Kernel model
591 spatialBg : `lsst.afw.math.Function2D`
592 Best-fit spatial background model
593 """
594 # What is the final kernel sum
595 kImage = afwImage.ImageD(spatialKernel.getDimensions())
596 kSum = spatialKernel.computeImage(kImage, False)
597 self.log.info("Final spatial kernel sum %.3f", kSum)
598
599 # Look at how well conditioned the matrix is
600 conditionNum = spatialSolution.getConditionNumber(
601 getattr(diffimLib.KernelSolution, self.kConfig.conditionNumberType))
602 self.log.info("Spatial model condition number %.3e", conditionNum)
603
604 if conditionNum < 0.0:
605 self.log.warning("Condition number is negative (%.3e)", conditionNum)
606 if conditionNum > self.kConfig.maxSpatialConditionNumber:
607 self.log.warning("Spatial solution exceeds max condition number (%.3e > %.3e)",
608 conditionNum, self.kConfig.maxSpatialConditionNumber)
609
610 self.metadata["spatialConditionNum"] = conditionNum
611 self.metadata["spatialKernelSum"] = kSum
612
613 # Look at how well the solution is constrained
614 nBasisKernels = spatialKernel.getNBasisKernels()
615 nKernelTerms = spatialKernel.getNSpatialParameters()
616 if nKernelTerms == 0: # order 0
617 nKernelTerms = 1
618
619 # Not fit for
620 nBgTerms = spatialBg.getNParameters()
621 if nBgTerms == 1:
622 if spatialBg.getParameters()[0] == 0.0:
623 nBgTerms = 0
624
625 nGood = 0
626 nBad = 0
627 nTot = 0
628 for cell in kernelCellSet.getCellList():
629 for cand in cell.begin(False): # False = include bad candidates
630 nTot += 1
631 if cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
632 nGood += 1
633 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
634 nBad += 1
635
636 self.log.info("Doing stats of kernel candidates used in the spatial fit.")
637
638 # Counting statistics
639 if nBad > 2*nGood:
640 self.log.warning("Many more candidates rejected than accepted; %d total, %d rejected, %d used",
641 nTot, nBad, nGood)
642 else:
643 self.log.info("%d candidates total, %d rejected, %d used", nTot, nBad, nGood)
644
645 # Some judgements on the quality of the spatial models
646 if nGood < nKernelTerms:
647 self.log.warning("Spatial kernel model underconstrained; %d candidates, %d terms, %d bases",
648 nGood, nKernelTerms, nBasisKernels)
649 self.log.warning("Consider lowering the spatial order")
650 elif nGood <= 2*nKernelTerms:
651 self.log.warning("Spatial kernel model poorly constrained; %d candidates, %d terms, %d bases",
652 nGood, nKernelTerms, nBasisKernels)
653 self.log.warning("Consider lowering the spatial order")
654 else:
655 self.log.info("Spatial kernel model well constrained; %d candidates, %d terms, %d bases",
656 nGood, nKernelTerms, nBasisKernels)
657
658 if nGood < nBgTerms:
659 self.log.warning("Spatial background model underconstrained; %d candidates, %d terms",
660 nGood, nBgTerms)
661 self.log.warning("Consider lowering the spatial order")
662 elif nGood <= 2*nBgTerms:
663 self.log.warning("Spatial background model poorly constrained; %d candidates, %d terms",
664 nGood, nBgTerms)
665 self.log.warning("Consider lowering the spatial order")
666 else:
667 self.log.info("Spatial background model appears well constrained; %d candidates, %d terms",
668 nGood, nBgTerms)
669
670 def _displayDebug(self, kernelCellSet, spatialKernel, spatialBackground):
671 """Provide visualization of the inputs and ouputs to the Psf-matching code
672
673 Parameters
674 ----------
675 kernelCellSet : `lsst.afw.math.SpatialCellSet`
676 The SpatialCellSet used in determining the matching kernel and background
677 spatialKernel : `lsst.afw.math.LinearCombinationKernel`
678 Spatially varying Psf-matching kernel
679 spatialBackground : `lsst.afw.math.Function2D`
680 Spatially varying background-matching function
681 """
682 import lsstDebug
683 displayCandidates = lsstDebug.Info(__name__).displayCandidates
684 displayKernelBasis = lsstDebug.Info(__name__).displayKernelBasis
685 displayKernelMosaic = lsstDebug.Info(__name__).displayKernelMosaic
686 plotKernelSpatialModel = lsstDebug.Info(__name__).plotKernelSpatialModel
687 plotKernelCoefficients = lsstDebug.Info(__name__).plotKernelCoefficients
688 showBadCandidates = lsstDebug.Info(__name__).showBadCandidates
689 maskTransparency = lsstDebug.Info(__name__).maskTransparency
690 if not maskTransparency:
691 maskTransparency = 0
692 afwDisplay.setDefaultMaskTransparency(maskTransparency)
693
694 if displayCandidates:
695 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground,
696 frame=lsstDebug.frame,
697 showBadCandidates=showBadCandidates)
698 lsstDebug.frame += 1
699 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground,
700 frame=lsstDebug.frame,
701 showBadCandidates=showBadCandidates,
702 kernels=True)
703 lsstDebug.frame += 1
704 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground,
705 frame=lsstDebug.frame,
706 showBadCandidates=showBadCandidates,
707 resids=True)
708 lsstDebug.frame += 1
709
710 if displayKernelBasis:
711 diutils.showKernelBasis(spatialKernel, frame=lsstDebug.frame)
712 lsstDebug.frame += 1
713
714 if displayKernelMosaic:
715 diutils.showKernelMosaic(kernelCellSet.getBBox(), spatialKernel, frame=lsstDebug.frame)
716 lsstDebug.frame += 1
717
718 if plotKernelSpatialModel:
719 diutils.plotKernelSpatialModel(spatialKernel, kernelCellSet, showBadCandidates=showBadCandidates)
720
721 if plotKernelCoefficients:
722 diutils.plotKernelCoefficients(spatialKernel, kernelCellSet)
723
724 def _createPcaBasis(self, kernelCellSet, nStarPerCell, ps):
725 """Create Principal Component basis
726
727 If a principal component analysis is requested, typically when using a delta function basis,
728 perform the PCA here and return a new basis list containing the new principal components.
729
730 Parameters
731 ----------
732 kernelCellSet : `lsst.afw.math.SpatialCellSet`
733 a SpatialCellSet containing KernelCandidates, from which components are derived
734 nStarPerCell : `int`
735 the number of stars per cell to visit when doing the PCA
736 ps : `lsst.daf.base.PropertySet`
737 input property set controlling the single kernel visitor
738
739 Returns
740 -------
741 nRejectedPca : `int`
742 number of KernelCandidates rejected during PCA loop
743 spatialBasisList : `list` of `lsst.afw.math.kernel.FixedKernel`
744 basis list containing the principal shapes as Kernels
745
746 Raises
747 ------
748 RuntimeError
749 If the Eigenvalues sum to zero.
750 """
751 nComponents = self.kConfig.numPrincipalComponents
752 imagePca = diffimLib.KernelPcaD()
753 importStarVisitor = diffimLib.KernelPcaVisitorF(imagePca)
754 kernelCellSet.visitCandidates(importStarVisitor, nStarPerCell)
755 if self.kConfig.subtractMeanForPca:
756 importStarVisitor.subtractMean()
757 imagePca.analyze()
758
759 eigenValues = imagePca.getEigenValues()
760 pcaBasisList = importStarVisitor.getEigenKernels()
761
762 eSum = np.sum(eigenValues)
763 if eSum == 0.0:
764 raise RuntimeError("Eigenvalues sum to zero")
765 trace_logger = getTraceLogger(self.log.getChild("_solve"), 5)
766 for j in range(len(eigenValues)):
767 trace_logger.debug("Eigenvalue %d : %f (%f)", j, eigenValues[j], eigenValues[j]/eSum)
768
769 nToUse = min(nComponents, len(eigenValues))
770 trimBasisList = []
771 for j in range(nToUse):
772 # Check for NaNs?
773 kimage = afwImage.ImageD(pcaBasisList[j].getDimensions())
774 pcaBasisList[j].computeImage(kimage, False)
775 if not (True in np.isnan(kimage.array)):
776 trimBasisList.append(pcaBasisList[j])
777
778 # Put all the power in the first kernel, which will not vary spatially
779 spatialBasisList = diffimLib.renormalizeKernelList(trimBasisList)
780
781 # New Kernel visitor for this new basis list (no regularization explicitly)
782 singlekvPca = diffimLib.BuildSingleKernelVisitorF(spatialBasisList, ps)
783 singlekvPca.setSkipBuilt(False)
784 kernelCellSet.visitCandidates(singlekvPca, nStarPerCell)
785 singlekvPca.setSkipBuilt(True)
786 nRejectedPca = singlekvPca.getNRejected()
787
788 return nRejectedPca, spatialBasisList
789
790 @abc.abstractmethod
791 def _buildCellSet(self, *args):
792 """Fill a SpatialCellSet with KernelCandidates for the Psf-matching process;
793 override in derived classes"""
794 return
795
796 @timeMethod
797 def _solve(self, kernelCellSet, basisList):
798 """Solve for the PSF matching kernel
799
800 Parameters
801 ----------
802 kernelCellSet : `lsst.afw.math.SpatialCellSet`
803 a SpatialCellSet to use in determining the matching kernel
804 (typically as provided by _buildCellSet)
805 basisList : `list` of `lsst.afw.math.kernel.FixedKernel`
806 list of Kernels to be used in the decomposition of the spatially varying kernel
807 (typically as provided by makeKernelBasisList)
808
809 Returns
810 -------
811 psfMatchingKernel : `lsst.afw.math.LinearCombinationKernel`
812 Spatially varying Psf-matching kernel
813 backgroundModel : `lsst.afw.math.Function2D`
814 Spatially varying background-matching function
815
816 Raises
817 ------
818 NoKernelCandidatesError :
819 If there are no useable kernel candidates.
820 """
821
822 import lsstDebug
823 display = lsstDebug.Info(__name__).display
824
825 maxSpatialIterations = self.kConfig.maxSpatialIterations
826 nStarPerCell = self.kConfig.nStarPerCell
827 usePcaForSpatialKernel = self.kConfig.usePcaForSpatialKernel
828
829 # Visitor for the single kernel fit
830 ps = pexConfig.makePropertySet(self.kConfig)
831 if self.useRegularization:
832 singlekv = diffimLib.BuildSingleKernelVisitorF(basisList, ps, self.hMat)
833 else:
834 singlekv = diffimLib.BuildSingleKernelVisitorF(basisList, ps)
835
836 # Visitor for the kernel sum rejection
837 ksv = diffimLib.KernelSumVisitorF(ps)
838
839 # Main loop
840 trace_loggers = [getTraceLogger(self.log.getChild("_solve"), i) for i in range(5)]
841 t0 = time.time()
842 totalIterations = 0
843 thisIteration = 0
844 while (thisIteration < maxSpatialIterations):
845
846 # Make sure there are no uninitialized candidates as active occupants of Cell
847 nRejectedSkf = -1
848 while (nRejectedSkf != 0):
849 trace_loggers[1].debug("Building single kernels...")
850 kernelCellSet.visitCandidates(singlekv, nStarPerCell, ignoreExceptions=True)
851 nRejectedSkf = singlekv.getNRejected()
852 trace_loggers[1].debug(
853 "Iteration %d, rejected %d candidates due to initial kernel fit",
854 thisIteration, nRejectedSkf
855 )
856
857 # Reject outliers in kernel sum
858 ksv.resetKernelSum()
859 ksv.setMode(diffimLib.KernelSumVisitorF.AGGREGATE)
860 kernelCellSet.visitCandidates(ksv, nStarPerCell, ignoreExceptions=True)
861
862 allCellsEmpty = True
863 for cell in kernelCellSet.getCellList():
864 if not cell.empty():
865 allCellsEmpty = False
866 break
867 if allCellsEmpty:
868 raise NoKernelCandidatesError("All spatial cells are emtpy of candidates")
869
870 try:
871 ksv.processKsumDistribution()
872 except RuntimeError as e:
873 raise NoKernelCandidatesError("Unable to calculate the kernel sum") from e
874
875 ksv.setMode(diffimLib.KernelSumVisitorF.REJECT)
876 kernelCellSet.visitCandidates(ksv, nStarPerCell, ignoreExceptions=True)
877
878 nRejectedKsum = ksv.getNRejected()
879 trace_loggers[1].debug(
880 "Iteration %d, rejected %d candidates due to kernel sum",
881 thisIteration, nRejectedKsum
882 )
883
884 # Do we jump back to the top without incrementing thisIteration?
885 if nRejectedKsum > 0:
886 totalIterations += 1
887 continue
888
889 # At this stage we can either apply the spatial fit to
890 # the kernels, or we run a PCA, use these as a *new*
891 # basis set with lower dimensionality, and then apply
892 # the spatial fit to these kernels
893
894 if (usePcaForSpatialKernel):
895 trace_loggers[0].debug("Building Pca basis")
896
897 nRejectedPca, spatialBasisList = self._createPcaBasis(kernelCellSet, nStarPerCell, ps)
898 trace_loggers[1].debug(
899 "Iteration %d, rejected %d candidates due to Pca kernel fit",
900 thisIteration, nRejectedPca
901 )
902
903 # We don't want to continue on (yet) with the
904 # spatial modeling, because we have bad objects
905 # contributing to the Pca basis. We basically
906 # need to restart from the beginning of this loop,
907 # since the cell-mates of those objects that were
908 # rejected need their original Kernels built by
909 # singleKernelFitter.
910
911 # Don't count against thisIteration
912 if (nRejectedPca > 0):
913 totalIterations += 1
914 continue
915 else:
916 spatialBasisList = basisList
917
918 # We have gotten on to the spatial modeling part
919 regionBBox = kernelCellSet.getBBox()
920 spatialkv = diffimLib.BuildSpatialKernelVisitorF(spatialBasisList, regionBBox, ps)
921 kernelCellSet.visitCandidates(spatialkv, nStarPerCell)
922 spatialkv.solveLinearEquation()
923 trace_loggers[2].debug("Spatial kernel built with %d candidates", spatialkv.getNCandidates())
924 spatialKernel, spatialBackground = spatialkv.getSolutionPair()
925
926 # Check the quality of the spatial fit (look at residuals)
927 assesskv = diffimLib.AssessSpatialKernelVisitorF(spatialKernel, spatialBackground, ps)
928 kernelCellSet.visitCandidates(assesskv, nStarPerCell)
929 nRejectedSpatial = assesskv.getNRejected()
930 nGoodSpatial = assesskv.getNGood()
931 trace_loggers[1].debug(
932 "Iteration %d, rejected %d candidates due to spatial kernel fit",
933 thisIteration, nRejectedSpatial
934 )
935 trace_loggers[1].debug("%d candidates used in fit", nGoodSpatial)
936
937 # If only nGoodSpatial == 0, might be other candidates in the cells
938 if nGoodSpatial == 0 and nRejectedSpatial == 0:
939 raise NoKernelCandidatesError("No kernel candidates for spatial fit")
940
941 if nRejectedSpatial == 0:
942 # Nothing rejected, finished with spatial fit
943 break
944
945 # Otherwise, iterate on...
946 thisIteration += 1
947
948 # Final fit if above did not converge
949 if (nRejectedSpatial > 0) and (thisIteration == maxSpatialIterations):
950 trace_loggers[1].debug("Final spatial fit")
951 if (usePcaForSpatialKernel):
952 nRejectedPca, spatialBasisList = self._createPcaBasis(kernelCellSet, nStarPerCell, ps)
953 regionBBox = kernelCellSet.getBBox()
954 spatialkv = diffimLib.BuildSpatialKernelVisitorF(spatialBasisList, regionBBox, ps)
955 kernelCellSet.visitCandidates(spatialkv, nStarPerCell)
956 spatialkv.solveLinearEquation()
957 trace_loggers[2].debug("Spatial kernel built with %d candidates", spatialkv.getNCandidates())
958 spatialKernel, spatialBackground = spatialkv.getSolutionPair()
959
960 # Run after the final fit to use the updated kernel visitor `spatialkv`
961 if spatialkv.getNCandidates() < 1:
962 raise NoKernelCandidatesError("No kernel candidates remain after max iterations")
963
964 spatialSolution = spatialkv.getKernelSolution()
965
966 t1 = time.time()
967 trace_loggers[0].debug("Total time to compute the spatial kernel : %.2f s", (t1 - t0))
968
969 if display:
970 self._displayDebug(kernelCellSet, spatialKernel, spatialBackground)
971
972 self._diagnostic(kernelCellSet, spatialSolution, spatialKernel, spatialBackground)
973
974 return spatialSolution, spatialKernel, spatialBackground
975
976
977PsfMatch = PsfMatchTask