580 def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
581 """Provide logging diagnostics on quality of spatial kernel fit
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
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)
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)
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)
610 self.metadata["spatialConditionNum"] = conditionNum
611 self.metadata["spatialKernelSum"] = kSum
613 # Look at how well the solution is constrained
614 nBasisKernels = spatialKernel.getNBasisKernels()
615 nKernelTerms = spatialKernel.getNSpatialParameters()
616 if nKernelTerms == 0: # order 0
620 nBgTerms = spatialBg.getNParameters()
622 if spatialBg.getParameters()[0] == 0.0:
628 for cell in kernelCellSet.getCellList():
629 for cand in cell.begin(False): # False = include bad candidates
631 if cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
633 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
636 self.log.info("Doing stats of kernel candidates used in the spatial fit.")
638 # Counting statistics
640 self.log.warning("Many more candidates rejected than accepted; %d total, %d rejected, %d used",
643 self.log.info("%d candidates total, %d rejected, %d used", nTot, nBad, nGood)
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")
655 self.log.info("Spatial kernel model well constrained; %d candidates, %d terms, %d bases",
656 nGood, nKernelTerms, nBasisKernels)
659 self.log.warning("Spatial background model underconstrained; %d candidates, %d terms",
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",
665 self.log.warning("Consider lowering the spatial order")
667 self.log.info("Spatial background model appears well constrained; %d candidates, %d terms",
670 def _displayDebug(self, kernelCellSet, spatialKernel, spatialBackground):
671 """Provide visualization of the inputs and ouputs to the Psf-matching code
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
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:
692 afwDisplay.setDefaultMaskTransparency(maskTransparency)
694 if displayCandidates:
695 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground,
696 frame=lsstDebug.frame,
697 showBadCandidates=showBadCandidates)
699 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground,
700 frame=lsstDebug.frame,
701 showBadCandidates=showBadCandidates,
704 diutils.showKernelCandidates(kernelCellSet, kernel=spatialKernel, background=spatialBackground,
705 frame=lsstDebug.frame,
706 showBadCandidates=showBadCandidates,
710 if displayKernelBasis:
711 diutils.showKernelBasis(spatialKernel, frame=lsstDebug.frame)
714 if displayKernelMosaic:
715 diutils.showKernelMosaic(kernelCellSet.getBBox(), spatialKernel, frame=lsstDebug.frame)
718 if plotKernelSpatialModel:
719 diutils.plotKernelSpatialModel(spatialKernel, kernelCellSet, showBadCandidates=showBadCandidates)
721 if plotKernelCoefficients:
722 diutils.plotKernelCoefficients(spatialKernel, kernelCellSet)
724 def _createPcaBasis(self, kernelCellSet, nStarPerCell, ps):
725 """Create Principal Component basis
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.
732 kernelCellSet : `lsst.afw.math.SpatialCellSet`
733 a SpatialCellSet containing KernelCandidates, from which components are derived
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
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
749 If the Eigenvalues sum to zero.
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()
759 eigenValues = imagePca.getEigenValues()
760 pcaBasisList = importStarVisitor.getEigenKernels()
762 eSum = np.sum(eigenValues)
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)
769 nToUse = min(nComponents, len(eigenValues))
771 for j in range(nToUse):
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])
778 # Put all the power in the first kernel, which will not vary spatially
779 spatialBasisList = diffimLib.renormalizeKernelList(trimBasisList)
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()
788 return nRejectedPca, spatialBasisList
797 def _solve(self, kernelCellSet, basisList):
798 """Solve for the PSF matching kernel
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)
811 psfMatchingKernel : `lsst.afw.math.LinearCombinationKernel`
812 Spatially varying Psf-matching kernel
813 backgroundModel : `lsst.afw.math.Function2D`
814 Spatially varying background-matching function
818 NoKernelCandidatesError :
819 If there are no useable kernel candidates.
823 display = lsstDebug.Info(__name__).display
825 maxSpatialIterations = self.kConfig.maxSpatialIterations
826 nStarPerCell = self.kConfig.nStarPerCell
827 usePcaForSpatialKernel = self.kConfig.usePcaForSpatialKernel
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)
834 singlekv = diffimLib.BuildSingleKernelVisitorF(basisList, ps)
836 # Visitor for the kernel sum rejection
837 ksv = diffimLib.KernelSumVisitorF(ps)
840 trace_loggers = [getTraceLogger(self.log.getChild("_solve"), i) for i in range(5)]
844 while (thisIteration < maxSpatialIterations):
846 # Make sure there are no uninitialized candidates as active occupants of Cell
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
857 # Reject outliers in kernel sum
859 ksv.setMode(diffimLib.KernelSumVisitorF.AGGREGATE)
860 kernelCellSet.visitCandidates(ksv, nStarPerCell, ignoreExceptions=True)
863 for cell in kernelCellSet.getCellList():
865 allCellsEmpty = False
868 raise NoKernelCandidatesError("All spatial cells are emtpy of candidates")
871 ksv.processKsumDistribution()
872 except RuntimeError as e:
873 raise NoKernelCandidatesError("Unable to calculate the kernel sum") from e
875 ksv.setMode(diffimLib.KernelSumVisitorF.REJECT)
876 kernelCellSet.visitCandidates(ksv, nStarPerCell, ignoreExceptions=True)
878 nRejectedKsum = ksv.getNRejected()
879 trace_loggers[1].debug(
880 "Iteration %d, rejected %d candidates due to kernel sum",
881 thisIteration, nRejectedKsum
884 # Do we jump back to the top without incrementing thisIteration?
885 if nRejectedKsum > 0:
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
894 if (usePcaForSpatialKernel):
895 trace_loggers[0].debug("Building Pca basis")
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
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.
911 # Don't count against thisIteration
912 if (nRejectedPca > 0):
916 spatialBasisList = basisList
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()
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
935 trace_loggers[1].debug("%d candidates used in fit", nGoodSpatial)
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")
941 if nRejectedSpatial == 0:
942 # Nothing rejected, finished with spatial fit
945 # Otherwise, iterate on...
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()
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")
964 spatialSolution = spatialkv.getKernelSolution()
967 trace_loggers[0].debug("Total time to compute the spatial kernel : %.2f s", (t1 - t0))
970 self._displayDebug(kernelCellSet, spatialKernel, spatialBackground)
972 self._diagnostic(kernelCellSet, spatialSolution, spatialKernel, spatialBackground)
974 return spatialSolution, spatialKernel, spatialBackground