DM-19443: Extract and visualize the local and the spatial AL kernel solution coefficients

doc/lsst.ip.diffim/AL_implementation.rst

@@ -0,0 +1,148 @@
+################################
+Alard-Lupton implementation code
+################################
+
+This page focuses on key implementation details of the Alard-Lupton
+(AL) image differencing in the LSST pipeline.
+
+Introduction
+------------
+
+In image differencing, our goal is to detect variable brightness
+features between two images by removing their static brightness parts
+by subtraction. Before performing pixel by pixel subtraction, the
+images should be transformed to match their PSFs. A suitable
+transformation is convolution of one of the images by a small image, a
+convolution kernel. The optimal convolution kernel is called the
+*matching kernel*. Given the images :math:`R` and :math:`I`, we want
+to minimise the following expression in the least squares sense by
+finding :math:`K`:
+
+.. math::
+   R\otimes K - I
+
+In the AL approach the matching kernel is searched in the form of a
+linear combination of given basis functions:
+
+.. math::
+   K = \sum_i a_i K_i
+
+As the PSFs are not uniform all over the images, the coefficients
+:math:`a_i` are considered to be functions of image coordinates. In
+the AL algorithm, first local kernel solutions (spatially not varying
+values of :math:`a_i`-s) are determined around selected sources. Then
+polynomials are fitted to the local sets of :math:`a_i` solutions
+to obtain the spatially dependent coefficients for the image pair.
+
+Basis functions
+---------------
+
+Solving for the matching kernel starts with the generation of the
+basis functions. The number of basis functions and the properties of
+the Gaussian bases (e.g. widths in pixels) are determined in
+`lsst.ip.diffim.makeKernelBasisList.generateAlardLuptonBasisList`. The
+decision logic depends on configuration and psf fwhm values of the
+images, some details is given in the function API documentation.
+
+.. TODO More details about the basis function generation
+
+Implementation of the AL solution
+---------------------------------
+
+In the LSST pipeline, the following main steps can be distinguished:
+
+* A subset of detected or external catalog sources are selected on a
+  spatial grid in the image. *Kernel candidates* (``KernelCandidate``
+  in :numref:`figclass`) are created to perform the determination of
+  the local matching kernel around each source. The entry point for
+  the kernel optimization is at `lsst.ip.diffim.PsfMatchTask._solve()`
+  in the Python codebase.
+
+* For each kernel candidate, local (spatially non-varying)
+  coefficients are determined of the basis kernels that minimizes the
+  image difference for the image stamp around this source.
+
+* If principal component analysis (PCA) is selected, new (fewer) basis
+  functions are calculated from the local solutions, then the kernel
+  candidates are solved again with the new set of basis functions.
+
+* For the whole image, spatially varying coefficients are fitted onto
+  the local coefficient solutions. This is an iterative process
+  rejecting outlying local solutions. The spatial variation is either
+  generic 2nd order polynomials, or Chebyshev polynomials of the first
+  kind (``config.spatialModelType``).
+
+:numref:`figclass` shows the main functional relation among the
+implementation classes. These classes are implemented in the C++
+codebase. Each arrow represents a functional relationship, pointing
+from the subject to its target (object) class. The targets are
+created, contained or accessed (used) by the subject classes. Numbers
+show the possible multiplicities of the target instances in each
+action. This diagram show the functionally most important components
+only. It does not detail the abstraction hierarchy of the classes, all
+classes shown have more general base classes. There are also
+additional enums and classes declared to represent internal statuses
+and to select and configure numerical algorithms. Classes for the
+statistical analysis of the local and spatial solutions and for
+supporting PCA (``KernelSumVisitor``, ``AssessSpatialKernelVisitor``,
+``KernelPca``, ``KernelPcaVisitor``) are also not detailed
+here. Please refer to the `C++ API
+<http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/index.html>`_
+documentation for these details.
+
+.. _figclass:
+
+.. figure:: figures/AL_kernel_solution_classes.svg
+    :align: center
+    :alt: Kernel solution classes
+
+    Main Alard-Lupton kernel solution classes in the C++ codebase and
+    their relationships.
+
+    Arrows point from the class that performs the arrow action on the
+    *pointed* class. Numbers mark target multiplicity, "0.." marks
+    zero or more, "1.." marks one or more target instances. All
+    containment or storing relations mean holding of *shared_ptrs* of
+    the object instances, there is no exclusive ownership.  Figure
+    source :download:`here
+    <figures/AL_kernel_solution_classes.drawio>`.
+
+In the C++ codebase, the solution steps are implemented using the
+*visitor concept*. The visitor classes visit each ``KernelCandidate``
+instance, perform their operation on the kernel candidate itself by
+changing the state of the *visited instance* while also updating their
+own *visitor* state with data about the whole visiting process like
+the number of good or bad numerical solutions.
+
+``BuildSingleKernelVisitor`` visits all kernel candidates in each
+``SpatialCell`` and calls their ``build()`` method to solve for
+the coefficients of the basis kernels. ``KernelCandidate`` owns the
+knowledge of how to initialize ``KernelSolution``, the solution knows
+how to solve itself and how to turn that into an output kernel. The
+solution is stored in the ``KernelCandidate`` instance
+itself. Following the solution of each kernel candidate with the
+initial basis kernels, a principal component analysis (PCA) can be
+performed to reduce the number of basis functions for the spatial
+solution ( ``config.usePcaForSpatialKernel`` ), typically in the case
+of the delta-function basis kernels. This needs the calculation of the
+PCA basis of the initial local solutions, recalculation of the local
+solutions using the new PCA basis kernels and then solving for the
+spatial coefficients of the PCA basis. To support PCA,
+``KernelCandidate`` can store one *original* and one *PCA* kernel
+solution. See the C++ API docs of ``KernelCandidate.build()``,
+``KernelPca`` and ``KernelPcaVisitor`` for more details.
+
+Following the determination of the kernel solutions for each kernel
+candidate, the spatial solution is determined by
+``BuildSpatialKernelVisitor``. ``BuildSpatialKernelVisitor`` visits
+the kernel candidates and collects their local solution and
+initializes ``SpatialKernelSolution``. Then ``SpatialKernelSolution``
+solves itself.
+
+Both the *local* and the whole image *spatial* solutions are stored as
+``LinearCombinationKernel`` instances. Note, the python API of
+``LinearCombinationKernel`` currently does not support the evaluation
+of parameters for an arbitrary x,y position. Workarounds are to get
+and evaluate the parameter functions themselves directly or to compute
+a kernel image that updates the last parameter values retrievable from
+the instance.

doc/lsst.ip.diffim/code_notes.rst

@@ -3,7 +3,7 @@ Package usage and technical notes
 #################################
 
 This page is a collection of usage and code related notes about the
-image differencing implmentation. We do not summarise the Alard-Lupton
+image differencing implementation. We do not summarise the Alard-Lupton
 (AL) [AL_1998]_ and Zackay, Ofek, Gal-Yam (ZOGY) [ZOGY2016]_ papers
 themselves here. 
 
@@ -30,13 +30,15 @@ steps included in ``ImageDifferenceTask`` around the actual
 subtraction operation of the images. Including the subtraction
 operation itself, these steps can be enabled or disabled by top level
 ``do<ACTION>`` configuration options. These top-level configuration
-options are summarised in Figures 1 and 2 (`flowchart source and
-standalone pdf version
+options are summarised in :numref:`flowchart1` and
+:numref:`flowchart2` (`flowchart source and standalone pdf version
 <https://github.com/lsst-dm/diffimTests/tree/master/figure_subtasks>`_
 ). Some of these top level configuration options are also passed on to
 invoked subtasks and influence their functionality. They may not be
 specified for the subtasks directly.
 
+.. _flowchart1:
+
 .. figure:: figures/ImageDifference_flowchart.draw.io-Page-1.svg
     :align: center
     :alt: Subtasks page 1
@@ -47,7 +49,10 @@ specified for the subtasks directly.
     Following the reading of the template and science images, the task
     starts with the preprocessing of the science exposure on the top
     and ends with post processing steps following the subtraction on
-    the bottom.
+    the bottom. Figure source
+    :download:`here <figures/ImageDifference_flowchart.drawio>`.
+
+.. _flowchart2:
 
 .. figure:: figures/ImageDifference_flowchart.draw.io-Page-2.svg
     :align: center

doc/lsst.ip.diffim/figures/AL_kernel_solution_classes.drawio

@@ -0,0 +1 @@
+<mxfile modified="2019-05-10T03:28:48.601Z" host="www.draw.io" agent="Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/9.3.1 Chrome/66.0.3359.181 Electron/3.0.6 Safari/537.36" etag="6kydxjYktaeKT2ujhCPw" version="10.6.7" type="device"><diagram id="69lAjpG5EIWKgLmL0d2c" name="Page-1">7V3ZcqO4Gn4aV3Uu4mK3fZk4mZw+ne5kkvQyV1MEZFunMaIBJ3Y//dGCWCQZ40RkmXFXl4OEFpB+ff8qMbCny/VF6ieLzygE0cAywvXAPhtYlmmMHPyH5Gx4jmGxnHkKwyKvyriFvwEvWOSuYAiyRsEcoSiHSTMzQHEMgryR56cpemwWm6Go2Wviz4GUcRv4kZz7HYb5guWOXaPK/w+A8wXvGb8gu7P0eeEiI1v4IXqsZdnnA3uaIpSzq+V6CiIyenxcWL0/ttwtHywFcd6lAp6Eb1Mz/PlfdDO5ebi+ukjR9bFVzMaDH62KN75N/Bz6Eal9C/Li2fMNH5DsES4jP8ap0xmK89vijonTfgTnMb4O8BOBFGc8gDSHeCxPihs5SnBusIBReOlv0Io8d5b7wU+eOl2gFP7GzfpR0Sa+neYFWVhuo8QtqYmzDZybggyXueaDYQpZn/11o+Cln+VFRoCiyE8yeF++xtJP5zA+RXmOlkUheayL4SdvCNa1rGLsLwBagjzd4CLFXXvEahQLgQ/8Y0VUFieqRY2gXE7LBR3Py4arqcYXxWzvM/O2NPPHOBlBPDKWgWb4p0YJEhngt87p9KToJ5iiCOH5PosRowsYRUIWJ40IzPKthJElfgDj+SUtc+ZUOTfFeJAshOvOIrqQFjAMQUwmFeV+7rMZJNOVIBjndMDcU/wfD+vUGLoDFz/4FKfNKo3/k+JpPkUxfhcf0vkFmDweASERxcy3r6Pd9MAJwO1GAJbXEwHYnkQAn0Aag2jqxyEM/Ry8qbXv/XPWvsmZ2665Nx2jr8kfKVc/HvF8hfs6wdcXV1d0vVjG6Ulx8fXLpy9X378MCE3aJhkcUuk3SBGBjJT8YOLAwIGbgQGjplsUrXKIYqFSVmYbMMM/VzcfL1gn19MT3trN+fT8y51Ys/aM9ysYEbgig2EZMcqrvEYdkhd+OCIlaT3c+wOg4gDpJl+Q65/0ccnzzMq2Aj9YgLDWmNDuGUyJ0IEnKAhAllHhpN7UI8wXMC67yOSxaFweMLYNY0vA6o6xXsd1Zjo9rTPHlNbZJYyBn07R8h7GPqGGTwW1HMBWC9haTbB1HVnScmwVEUz6AltHFrWk2Y4gnWm21rmmYT5poS/xko1AtbLvCDGcHZvS6rfl1W8rVnrk34PoGmWQQpd9lrKyAgJoWuRO9/nevsiV8zvua3ZlFUoUpE8uCQ/yM8rpGHMgVx/+gGsQsuV/JHAWLISROpRrMDE82tBO0g2eLs4gMUEcF1lDygYbbeDB9vE7gpR2FiAwm8EA4vHPjhgbnJH+yb14tbzHxY54u+ouP1hnpNIqDgglZEfb+eKnkpcuiX5NGGn5OphSg1WExUvSs0+ZZ0z6+bDGLGNDHiHhtNZs9Lb2UCFIQBwSUiKqf/VirKtCECC9hZRD0zoBWiYr1u09SSebfEEFEPxE0VDojAJXbQj5BGfNl0H3OSZ1EIrVi2etD5e6In198gtI88WoZID0m4NoMzyICi+DIrbTkUu4feljjivhyCkRW2/xiEeArahvEK8MskIPsoIWWcFxx03NjFtF6hKjqSIDozcykNVyDeyEKDIMxBnaEAWIpRMUwWBzwJbu2OJ2prIWNURFVG5vWshYSVNJiojWWtp7iIJ8mHu9cz9WmHpUc9+bmc+ZSHPfZtc9cJIncpKRYN8fyvbdUsOsT3yb1ln0doMlWJ8IAbXuuGzCxZeR3Bu3MLeRmR/hqYzx2j9FqzjMJGorX/TpBMgNLi0MTbA7Zwdu1B2RJp0pdTs3UhKmrYEbefbF5eiX7X/++090fRemxvcgP7YcJUF8xJRIh5JYS1O0mi8G3G4aHCjjJSjDk3mVmjI0GMqUlOHKKtCU8AJItGcvIi9/j5Ufb06uGlJwDU3EchKdgHAOOF/Dw7xAcxT70XmViycPgyEIi6mrylwiakojXOR/IM83BePyVzkirCtfcrYG1jD/QaoPR26R/KtojVyfreuJDU/EeAx/1BO0zpC6vWm6qkdTvCJl1OxRTKuNnWVolQagbWmycpgpz0HeUs4sGiRD2UpZKYiwrPHQDKbQjyiyiKP0QB1knT4s7LbR0c2iw8CunP/doQxPUJmXfp7CdakxP2AhjNpGMR/yebMR9eUwVQqP+lJoovQDKhx/hRP0gWpdzG0I1iBYEbe74UdzTBv5YnngdyqKbQUB7d5Bry9ByJSdgzV6PMCVrugLz2qqZq5CzvGURj4NPiP1zKsRy69jRuE6UZn/G2iWoGgToyXzNtQKFS6SgOAJrAvONT8IXmuolLITTHtRdoCZdpgxdTgoleQ26g1nZJOyNMv/Qv+zpuktI353zW9vYDJ5TRWHXzNdhWs7ag2nk2LSGq1a10xarU36NBNa9SRN/U2tQEF7VcvXJKPGcwQZ2RNis3cULyJGKqpg/asrO56hqlw+KxuzotKz7HtquUuWYrYr7cZw2EFDL1hOjdqo8LAV98+6iEAlLFE8OcUSzJwuAg5umD3O6L89wGKbiaXcj1A88qAe8q+SUYyhPXaFKFGW2o9OJeI4NhuNlpIHbwHNZhnohy7c14QlYw9YKg0v5tAwnKbxxbTtHdYXmroGKWRxIs9FOu6r24l0PBzyrdhgZCGDRisUcS/cfkvlWGH9T1PAbLvcOY0VbZos1WZVra9JYRGua9dCGfM9gA2Vd5hooyad1uWlBXscp4k9XJl+JvgIFXrEmtGrYE1psB26DdTYhRilddgZeQ2cmoyddqTCCX1Qw52GO8291puQqUxX0OMLqFRKSdooSw6daJVu3gXidIUYvklEA8R4Ix7YxGdPj3jTJIkSFfuHHFt2az9X7H3DkPUMeWaiW57hU49fYSSE7xlaiMoxmwqVJ/gt2DtLCtXeOuJI2c1WHVEoPnFbi5uC1bNZfJtK+QLrhtvpxTBXFpR0iHPtJTrJFFx2YzleSB3nynFav8+uF6fdIc5VmyW0XKhvMs5V/chy/IgyzlWINxOJiLpl2Ta9818rqpt9ODpQTo+U0zVK1ukLjRSbuLZZJRQ2BvV+qapgjVa8Xyty4AdlX8cZZS1kJ5ZpJevq5puQEMsEU02Nyb5hSE83jKmpyusoSPLNP10FSaycuIbb1CPsZ0qS2+WofaXEkSCMuk0xsR+tl29yegu2W6vdeKsv+E2/QfVJE+46TVrEFNAq6EvlvU6SPq/NBTyh9kt5j1z5IIrtarT5PrxHu0myddE939xCVGPRWPY8PCu9SYKd+OXcSYrjahoq4yHKU7PO6E6EHS3cM/V6YZ7qQ2sOcZ7/Lu2it2Ngegv0tMevKk4Nu+9B6CQ9PU/q7ld6KlVEHgwz3mFVFTjavuWtcYsTSq4thqKOhaCwnqUtW94M+u6lrZr7/HiyhwOdL0od3q2RxQPUNAtY3HvSv3zlOK+LUXtslOoEUq2Wqbem4u0FIpLK9sIgojio5N2DSEfM4ItEB2ZY46ZDXFO8Hwci3qo3fDGXOJep3oLVqGPEH54GwxMMncZkh6FTd8SfwrD5Qrsut0gpY0HGcQWNbgvA7O0AnzRh0Bx522FQG5HuE678z0Kv8iQGHQYma2Q2kUaPADQR0OvFsMuSyOJVwKwCpslEACZ78qRQZI0GclcGqtbAwj1CeWxTYFzHVm8uGG00Ix9kcef/pC7D8uRiATsyui+9aSl68/AiUlBN49pH4SqXmBb4sYw+xCfOYjfNpFZaunv4evdrcw9Wx1n0bfzrz/X87Ixv03hZtHm6zNJ5hferSAlyitFQo3aU9ox9/GQCY+KVe9C5lOQhc6frFD3QT6PwM05lJGkcSv7et0ZsQZq2xaRDS6s+WKNVzBGP7+pDR/v09fPi+mzkjm+md+4s/DG+iv7ihPQGcaY1SrBvlccTVBHb7mZT2RewPGGjjeO0m5vF8vYO87T0HmO90SJKmpIjOcvYKc9fEndTfJ8lA9Umrp2hVe8Ko7bZn9tWog6QMhwhfkkPSFnmUGi2D3VMOTiqEwt4dFyhri/p6d9VyBz7jghcJhEgt6pzvgl1sRjgqni2wNMT/p3kaZV3nCmU/S/sZHTSzAI9Dujh5Mz3bhnLFfkOWQQDmEOgqizn3LHvgPhL+gfs8lvjh6Q9k4clySBFjKGzY/EiP8tox8ZH0gj7cgptYpXRz5aIqw3ficHjQDo3jV4HdNHSfmbsxHPWDYxhcZZ56ecOqk8ElZpMF0sJfpOEXJIT11n8xOMC5uA28SkXeEz9pLmWk5oyW1Sq6ben/KR4a8B205GZLUIv7Nptetdt5UTd4yK8SZN3m4rjPZQfctLxkRHlYjFf1+46dMbCLsYJv695F2Pb4V99Swim09xr5IqxyJokBNNpahqON2nl+OJzCeX74fiKk6s6h0sPh6RmTL+8oIybfsMMfh/+bmq1tk48wTKhhcPbhsDhTQ2KCE5WX7VkxauPg9rn/wc=</diagram></mxfile>