Coverage for python/lsst/ip/diffim/dipoleFitTask.py: 10%
473 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-05-29 08:41 +0000
1#
2# LSST Data Management System
3# Copyright 2008-2016 AURA/LSST.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
23import math
24import logging
25import numpy as np
26import warnings
28import lsst.afw.image as afwImage
29import lsst.meas.base as measBase
30import lsst.afw.detection as afwDet
31import lsst.geom as geom
32import lsst.pex.exceptions as pexExcept
33import lsst.pex.config as pexConfig
34from lsst.pipe.base import Struct
35from lsst.utils.timer import timeMethod
37__all__ = ("DipoleFitTask", "DipoleFitPlugin", "DipoleFitTaskConfig", "DipoleFitPluginConfig",
38 "DipoleFitAlgorithm")
41# Create a new measurement task (`DipoleFitTask`) that can handle all other SFM tasks but can
42# pass a separate pos- and neg- exposure/image to the `DipoleFitPlugin`s `run()` method.
45class DipoleFitPluginConfig(measBase.SingleFramePluginConfig):
46 """Configuration for DipoleFitPlugin
47 """
49 fitAllDiaSources = pexConfig.Field(
50 dtype=bool, default=False,
51 doc="""Attempte dipole fit of all diaSources (otherwise just the ones consisting of overlapping
52 positive and negative footprints)""")
54 maxSeparation = pexConfig.Field(
55 dtype=float, default=5.,
56 doc="Assume dipole is not separated by more than maxSeparation * psfSigma")
58 relWeight = pexConfig.Field(
59 dtype=float, default=0.5,
60 doc="""Relative weighting of pre-subtraction images (higher -> greater influence of pre-sub.
61 images on fit)""")
63 tolerance = pexConfig.Field(
64 dtype=float, default=1e-7,
65 doc="Fit tolerance")
67 fitBackground = pexConfig.Field(
68 dtype=int, default=1,
69 doc="Set whether and how to fit for linear gradient in pre-sub. images. Possible values:"
70 "0: do not fit background at all"
71 "1 (default): pre-fit the background using linear least squares and then do not fit it as part"
72 "of the dipole fitting optimization"
73 "2: pre-fit the background using linear least squares (as in 1), and use the parameter"
74 "estimates from that fit as starting parameters for an integrated re-fit of the background")
76 fitSeparateNegParams = pexConfig.Field(
77 dtype=bool, default=False,
78 doc="Include parameters to fit for negative values (flux, gradient) separately from pos.")
80 # Config params for classification of detected diaSources as dipole or not
81 minSn = pexConfig.Field(
82 dtype=float, default=math.sqrt(2) * 5.0,
83 doc="Minimum quadrature sum of positive+negative lobe S/N to be considered a dipole")
85 maxFluxRatio = pexConfig.Field(
86 dtype=float, default=0.65,
87 doc="Maximum flux ratio in either lobe to be considered a dipole")
89 maxChi2DoF = pexConfig.Field(
90 dtype=float, default=0.05,
91 doc="""Maximum Chi2/DoF significance of fit to be considered a dipole.
92 Default value means \"Choose a chi2DoF corresponding to a significance level of at most 0.05\"
93 (note this is actually a significance, not a chi2 value).""")
95 maxFootprintArea = pexConfig.Field(
96 dtype=int, default=1_200,
97 doc=("Maximum area for footprints before they are ignored as large; "
98 "non-positive means no threshold applied"
99 "Threshold chosen for HSC and DECam data, see DM-38741 for details."))
102class DipoleFitTaskConfig(measBase.SingleFrameMeasurementConfig):
104 def setDefaults(self):
105 measBase.SingleFrameMeasurementConfig.setDefaults(self)
107 self.plugins.names = ["base_SdssCentroid",
108 "ip_diffim_DipoleFit",
109 "base_CircularApertureFlux",
110 "base_PixelFlags",
111 "base_SkyCoord",
112 "base_PsfFlux",
113 "base_SdssShape",
114 ]
115 # Only measure the apertures we need to report in the alert stream.
116 self.plugins["base_CircularApertureFlux"].radii = [12.0]
118 self.slots.calibFlux = None
119 self.slots.modelFlux = None
120 self.slots.gaussianFlux = None
121 self.slots.shape = "base_SdssShape"
122 # This will be switched to "ip_diffim_DipoleFit" as this task runs.
123 self.slots.centroid = "base_SdssCentroid"
124 self.doReplaceWithNoise = False
127class DipoleFitTask(measBase.SingleFrameMeasurementTask):
128 """A task that fits a dipole to a difference image, with an optional
129 separate detection image.
131 Because it subclasses SingleFrameMeasurementTask, and calls
132 SingleFrameMeasurementTask.run() from its run() method, it still
133 can be used identically to a standard SingleFrameMeasurementTask.
134 """
136 ConfigClass = DipoleFitTaskConfig
137 _DefaultName = "dipoleFit"
139 def __init__(self, schema, algMetadata=None, **kwargs):
140 super().__init__(schema, algMetadata, **kwargs)
142 # Enforce a specific plugin order, so that DipoleFit can fall back on
143 # SdssCentroid for non-dipoles
144 self.plugins_pre = self.plugins.copy()
145 self.plugins_post = self.plugins.copy()
146 self.plugins_pre.clear()
147 self.plugins_pre["base_SdssCentroid"] = self.plugins["base_SdssCentroid"]
148 self.plugins_post.pop("base_SdssCentroid")
149 self.dipoleFit = self.plugins_post.pop("ip_diffim_DipoleFit")
150 del self.plugins
152 @timeMethod
153 def run(self, sources, exposure, posExp=None, negExp=None, **kwargs):
154 """Run dipole measurement and classification.
156 Run SdssCentroid first, then switch the centroid slot, then DipoleFit
157 then the rest; DipoleFit will fall back on SdssCentroid for sources
158 not containing positive+negative peaks.
160 Parameters
161 ----------
162 sources : `lsst.afw.table.SourceCatalog`
163 ``diaSources`` that will be measured using dipole measurement.
164 exposure : `lsst.afw.image.Exposure`
165 The difference exposure on which the ``sources`` were detected.
166 If neither ``posExp`` nor ``negExp`` are set, then the dipole is also
167 fitted directly to this difference image.
168 posExp : `lsst.afw.image.Exposure`, optional
169 "Positive" exposure, typically a science exposure, or None if unavailable
170 When `posExp` is `None`, will compute `posImage = exposure + negExp`.
171 negExp : `lsst.afw.image.Exposure`, optional
172 "Negative" exposure, typically a template exposure, or None if unavailable
173 When `negExp` is `None`, will compute `negImage = posExp - exposure`.
174 **kwargs
175 Additional keyword arguments for `lsst.meas.base.sfm.SingleFrameMeasurementTask`.
176 """
177 # Run plugins in a very specific order, so DipoleFitPlugin has a
178 # centroid to fall back on.
179 self.plugins = self.plugins_pre
180 super().run(sources, exposure, **kwargs)
182 for source in sources:
183 self.dipoleFit.measureDipoles(source, exposure, posExp, negExp)
184 # Use the new DipoleFit outputs for subsequent measurements, now that
185 # non-dipoles have been filled in with the earlier centroid values.
186 sources.schema.getAliasMap().set("slot_Centroid", "ip_diffim_DipoleFit")
188 self.plugins = self.plugins_post
189 super().run(sources, exposure, **kwargs)
192class DipoleModel:
193 """Lightweight class containing methods for generating a dipole model for fitting
194 to sources in diffims, used by DipoleFitAlgorithm.
196 See also:
197 `DMTN-007: Dipole characterization for image differencing <https://dmtn-007.lsst.io>`_.
198 """
200 def __init__(self):
201 import lsstDebug
202 self.debug = lsstDebug.Info(__name__).debug
203 self.log = logging.getLogger(__name__)
205 def makeBackgroundModel(self, in_x, pars=None):
206 """Generate gradient model (2-d array) with up to 2nd-order polynomial
208 Parameters
209 ----------
210 in_x : `numpy.array`
211 (2, w, h)-dimensional `numpy.array`, containing the
212 input x,y meshgrid providing the coordinates upon which to
213 compute the gradient. This will typically be generated via
214 `_generateXYGrid()`. `w` and `h` correspond to the width and
215 height of the desired grid.
216 pars : `list` of `float`, optional
217 Up to 6 floats for up
218 to 6 2nd-order 2-d polynomial gradient parameters, in the
219 following order: (intercept, x, y, xy, x**2, y**2). If `pars`
220 is emtpy or `None`, do nothing and return `None` (for speed).
222 Returns
223 -------
224 result : `None` or `numpy.array`
225 return None, or 2-d numpy.array of width/height matching
226 input bbox, containing computed gradient values.
227 """
229 # Don't fit for other gradient parameters if the intercept is not included.
230 if (pars is None) or (len(pars) <= 0) or (pars[0] is None):
231 return
233 y, x = in_x[0, :], in_x[1, :]
234 gradient = np.full_like(x, pars[0], dtype='float64')
235 if len(pars) > 1 and pars[1] is not None:
236 gradient += pars[1] * x
237 if len(pars) > 2 and pars[2] is not None:
238 gradient += pars[2] * y
239 if len(pars) > 3 and pars[3] is not None:
240 gradient += pars[3] * (x * y)
241 if len(pars) > 4 and pars[4] is not None:
242 gradient += pars[4] * (x * x)
243 if len(pars) > 5 and pars[5] is not None:
244 gradient += pars[5] * (y * y)
246 return gradient
248 def _generateXYGrid(self, bbox):
249 """Generate a meshgrid covering the x,y coordinates bounded by bbox
251 Parameters
252 ----------
253 bbox : `lsst.geom.Box2I`
254 input Bounding Box defining the coordinate limits
256 Returns
257 -------
258 in_x : `numpy.array`
259 (2, w, h)-dimensional numpy array containing the grid indexing over x- and
260 y- coordinates
261 """
263 x, y = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
264 in_x = np.array([y, x]).astype(np.float64)
265 in_x[0, :] -= np.mean(in_x[0, :])
266 in_x[1, :] -= np.mean(in_x[1, :])
267 return in_x
269 def _getHeavyFootprintSubimage(self, fp, badfill=np.nan, grow=0):
270 """Extract the image from a ``~lsst.afw.detection.HeavyFootprint``
271 as an `lsst.afw.image.ImageF`.
273 Parameters
274 ----------
275 fp : `lsst.afw.detection.HeavyFootprint`
276 HeavyFootprint to use to generate the subimage
277 badfill : `float`, optional
278 Value to fill in pixels in extracted image that are outside the footprint
279 grow : `int`
280 Optionally grow the footprint by this amount before extraction
282 Returns
283 -------
284 subim2 : `lsst.afw.image.ImageF`
285 An `~lsst.afw.image.ImageF` containing the subimage.
286 """
287 bbox = fp.getBBox()
288 if grow > 0:
289 bbox.grow(grow)
291 subim2 = afwImage.ImageF(bbox, badfill)
292 fp.getSpans().unflatten(subim2.array, fp.getImageArray(), bbox.getCorners()[0])
293 return subim2
295 def fitFootprintBackground(self, source, posImage, order=1):
296 """Fit a linear (polynomial) model of given order (max 2) to the background of a footprint.
298 Only fit the pixels OUTSIDE of the footprint, but within its bounding box.
300 Parameters
301 ----------
302 source : `lsst.afw.table.SourceRecord`
303 SourceRecord, the footprint of which is to be fit
304 posImage : `lsst.afw.image.Exposure`
305 The exposure from which to extract the footprint subimage
306 order : `int`
307 Polynomial order of background gradient to fit.
309 Returns
310 -------
311 pars : `tuple` of `float`
312 `tuple` of length (1 if order==0; 3 if order==1; 6 if order == 2),
313 containing the resulting fit parameters
314 """
316 # TODO look into whether to use afwMath background methods -- see
317 # http://lsst-web.ncsa.illinois.edu/doxygen/x_masterDoxyDoc/_background_example.html
318 fp = source.getFootprint()
319 bbox = fp.getBBox()
320 bbox.grow(3)
321 posImg = afwImage.ImageF(posImage.image, bbox, afwImage.PARENT)
323 # This code constructs the footprint image so that we can identify the pixels that are
324 # outside the footprint (but within the bounding box). These are the pixels used for
325 # fitting the background.
326 posHfp = afwDet.HeavyFootprintF(fp, posImage.getMaskedImage())
327 posFpImg = self._getHeavyFootprintSubimage(posHfp, grow=3)
329 isBg = np.isnan(posFpImg.array).ravel()
331 data = posImg.array.ravel()
332 data = data[isBg]
333 B = data
335 x, y = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
336 x = x.astype(np.float64).ravel()
337 x -= np.mean(x)
338 x = x[isBg]
339 y = y.astype(np.float64).ravel()
340 y -= np.mean(y)
341 y = y[isBg]
342 b = np.ones_like(x, dtype=np.float64)
344 M = np.vstack([b]).T # order = 0
345 if order == 1:
346 M = np.vstack([b, x, y]).T
347 elif order == 2:
348 M = np.vstack([b, x, y, x**2., y**2., x*y]).T
350 pars = np.linalg.lstsq(M, B, rcond=-1)[0]
351 return pars
353 def makeStarModel(self, bbox, psf, xcen, ycen, flux):
354 """Generate a 2D image model of a single PDF centered at the given coordinates.
356 Parameters
357 ----------
358 bbox : `lsst.geom.Box`
359 Bounding box marking pixel coordinates for generated model
360 psf : TODO: DM-17458
361 Psf model used to generate the 'star'
362 xcen : `float`
363 Desired x-centroid of the 'star'
364 ycen : `float`
365 Desired y-centroid of the 'star'
366 flux : `float`
367 Desired flux of the 'star'
369 Returns
370 -------
371 p_Im : `lsst.afw.image.Image`
372 2-d stellar image of width/height matching input ``bbox``,
373 containing PSF with given centroid and flux
374 """
376 # Generate the psf image, normalize to flux
377 psf_img = psf.computeImage(geom.Point2D(xcen, ycen)).convertF()
378 psf_img_sum = np.nansum(psf_img.array)
379 psf_img *= (flux/psf_img_sum)
381 # Clip the PSF image bounding box to fall within the footprint bounding box
382 psf_box = psf_img.getBBox()
383 psf_box.clip(bbox)
384 psf_img = afwImage.ImageF(psf_img, psf_box, afwImage.PARENT)
386 # Then actually crop the psf image.
387 # Usually not necessary, but if the dipole is near the edge of the image...
388 # Would be nice if we could compare original pos_box with clipped pos_box and
389 # see if it actually was clipped.
390 p_Im = afwImage.ImageF(bbox)
391 tmpSubim = afwImage.ImageF(p_Im, psf_box, afwImage.PARENT)
392 tmpSubim += psf_img
394 return p_Im
396 def makeModel(self, x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None,
397 b=None, x1=None, y1=None, xy=None, x2=None, y2=None,
398 bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None,
399 **kwargs):
400 """Generate dipole model with given parameters.
402 This is the function whose sum-of-squared difference from data
403 is minimized by `lmfit`.
405 x : TODO: DM-17458
406 Input independent variable. Used here as the grid on
407 which to compute the background gradient model.
408 flux : `float`
409 Desired flux of the positive lobe of the dipole
410 xcenPos, ycenPos : `float`
411 Desired x,y-centroid of the positive lobe of the dipole
412 xcenNeg, ycenNeg : `float`
413 Desired x,y-centroid of the negative lobe of the dipole
414 fluxNeg : `float`, optional
415 Desired flux of the negative lobe of the dipole, set to 'flux' if None
416 b, x1, y1, xy, x2, y2 : `float`
417 Gradient parameters for positive lobe.
418 bNeg, x1Neg, y1Neg, xyNeg, x2Neg, y2Neg : `float`, optional
419 Gradient parameters for negative lobe.
420 They are set to the corresponding positive values if None.
422 **kwargs : `dict` [`str`]
423 Keyword arguments passed through ``lmfit`` and
424 used by this function. These must include:
426 - ``psf`` Psf model used to generate the 'star'
427 - ``rel_weight`` Used to signify least-squares weighting of posImage/negImage
428 relative to diffim. If ``rel_weight == 0`` then posImage/negImage are ignored.
429 - ``bbox`` Bounding box containing region to be modelled
431 Returns
432 -------
433 zout : `numpy.array`
434 Has width and height matching the input bbox, and
435 contains the dipole model with given centroids and flux(es). If
436 ``rel_weight`` = 0, this is a 2-d array with dimensions matching
437 those of bbox; otherwise a stack of three such arrays,
438 representing the dipole (diffim), positive, and negative images
439 respectively.
440 """
442 psf = kwargs.get('psf')
443 rel_weight = kwargs.get('rel_weight') # if > 0, we're including pre-sub. images
444 fp = kwargs.get('footprint')
445 bbox = fp.getBBox()
447 if fluxNeg is None:
448 fluxNeg = flux
450 self.log.debug('flux: %.2f fluxNeg: %.2f x+: %.2f x-: %.2f y+: %.2f y-: %.2f ',
451 flux, fluxNeg, xcenPos, xcenNeg, ycenPos, ycenNeg)
452 if x1 is not None:
453 self.log.debug(' b: %.2f x1: %.2f y1: %.2f', b, x1, y1)
454 if xy is not None:
455 self.log.debug(' xy: %.2f x2: %.2f y2: %.2f', xy, x2, y2)
457 posIm = self.makeStarModel(bbox, psf, xcenPos, ycenPos, flux)
458 negIm = self.makeStarModel(bbox, psf, xcenNeg, ycenNeg, fluxNeg)
460 in_x = x
461 if in_x is None: # use the footprint to generate the input grid
462 y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
463 in_x = np.array([x, y]) * 1.
464 in_x[0, :] -= in_x[0, :].mean() # center it!
465 in_x[1, :] -= in_x[1, :].mean()
467 if b is not None:
468 gradient = self.makeBackgroundModel(in_x, (b, x1, y1, xy, x2, y2))
470 # If bNeg is None, then don't fit the negative background separately
471 if bNeg is not None:
472 gradientNeg = self.makeBackgroundModel(in_x, (bNeg, x1Neg, y1Neg, xyNeg, x2Neg, y2Neg))
473 else:
474 gradientNeg = gradient
476 posIm.array[:, :] += gradient
477 negIm.array[:, :] += gradientNeg
479 # Generate the diffIm model
480 diffIm = afwImage.ImageF(bbox)
481 diffIm += posIm
482 diffIm -= negIm
484 zout = diffIm.array
485 if rel_weight > 0.:
486 zout = np.append([zout], [posIm.array, negIm.array], axis=0)
488 return zout
491class DipoleFitAlgorithm:
492 """Fit a dipole model using an image difference.
494 See also:
495 `DMTN-007: Dipole characterization for image differencing <https://dmtn-007.lsst.io>`_.
496 """
498 # This is just a private version number to sync with the ipython notebooks that I have been
499 # using for algorithm development.
500 _private_version_ = '0.0.5'
502 # Below is a (somewhat incomplete) list of improvements
503 # that would be worth investigating, given the time:
505 # todo 1. evaluate necessity for separate parameters for pos- and neg- images
506 # todo 2. only fit background OUTSIDE footprint (DONE) and dipole params INSIDE footprint (NOT DONE)?
507 # todo 3. correct normalization of least-squares weights based on variance planes
508 # todo 4. account for PSFs that vary across the exposures (should be happening by default?)
509 # todo 5. correctly account for NA/masks (i.e., ignore!)
510 # todo 6. better exception handling in the plugin
511 # todo 7. better classification of dipoles (e.g. by comparing chi2 fit vs. monopole?)
512 # todo 8. (DONE) Initial fast estimate of background gradient(s) params -- perhaps using numpy.lstsq
513 # todo 9. (NOT NEEDED - see (2)) Initial fast test whether a background gradient needs to be fit
514 # todo 10. (DONE) better initial estimate for flux when there's a strong gradient
515 # todo 11. (DONE) requires a new package `lmfit` -- investiate others? (astropy/scipy/iminuit?)
517 def __init__(self, diffim, posImage=None, negImage=None):
518 """Algorithm to run dipole measurement on a diaSource
520 Parameters
521 ----------
522 diffim : `lsst.afw.image.Exposure`
523 Exposure on which the diaSources were detected
524 posImage : `lsst.afw.image.Exposure`
525 "Positive" exposure from which the template was subtracted
526 negImage : `lsst.afw.image.Exposure`
527 "Negative" exposure which was subtracted from the posImage
528 """
530 self.diffim = diffim
531 self.posImage = posImage
532 self.negImage = negImage
533 self.psfSigma = None
534 if diffim is not None:
535 diffimPsf = diffim.getPsf()
536 diffimAvgPos = diffimPsf.getAveragePosition()
537 self.psfSigma = diffimPsf.computeShape(diffimAvgPos).getDeterminantRadius()
539 self.log = logging.getLogger(__name__)
541 import lsstDebug
542 self.debug = lsstDebug.Info(__name__).debug
544 def fitDipoleImpl(self, source, tol=1e-7, rel_weight=0.5,
545 fitBackground=1, bgGradientOrder=1, maxSepInSigma=5.,
546 separateNegParams=True, verbose=False):
547 """Fit a dipole model to an input difference image.
549 Actually, fits the subimage bounded by the input source's
550 footprint) and optionally constrain the fit using the
551 pre-subtraction images posImage and negImage.
553 Parameters
554 ----------
555 source : TODO: DM-17458
556 TODO: DM-17458
557 tol : float, optional
558 TODO: DM-17458
559 rel_weight : `float`, optional
560 TODO: DM-17458
561 fitBackground : `int`, optional
562 TODO: DM-17458
563 bgGradientOrder : `int`, optional
564 TODO: DM-17458
565 maxSepInSigma : `float`, optional
566 TODO: DM-17458
567 separateNegParams : `bool`, optional
568 TODO: DM-17458
569 verbose : `bool`, optional
570 TODO: DM-17458
572 Returns
573 -------
574 result : `lmfit.MinimizerResult`
575 return `lmfit.MinimizerResult` object containing the fit
576 parameters and other information.
577 """
579 # Only import lmfit if someone wants to use the new DipoleFitAlgorithm.
580 import lmfit
582 fp = source.getFootprint()
583 bbox = fp.getBBox()
584 subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)
586 z = diArr = subim.image.array
587 # Make sure we don't overwrite buffers.
588 z = z.copy()
589 weights = 1. / subim.variance.array # get the weights (=1/variance)
591 if rel_weight > 0. and ((self.posImage is not None) or (self.negImage is not None)):
592 if self.negImage is not None:
593 negSubim = afwImage.MaskedImageF(self.negImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
594 if self.posImage is not None:
595 posSubim = afwImage.MaskedImageF(self.posImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
596 if self.posImage is None: # no science image provided; generate it from diffim + negImage
597 posSubim = subim.clone()
598 posSubim += negSubim
599 if self.negImage is None: # no template provided; generate it from the posImage - diffim
600 negSubim = posSubim.clone()
601 negSubim -= subim
603 z = np.append([z], [posSubim.image.array,
604 negSubim.image.array], axis=0)
605 # Weight the pos/neg images by rel_weight relative to the diffim
606 weights = np.append([weights], [1. / posSubim.variance.array * rel_weight,
607 1. / negSubim.variance.array * rel_weight], axis=0)
608 else:
609 rel_weight = 0. # a short-cut for "don't include the pre-subtraction data"
611 # It seems that `lmfit` requires a static functor as its optimized method, which eliminates
612 # the ability to pass a bound method or other class method. Here we write a wrapper which
613 # makes this possible.
614 def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None,
615 b=None, x1=None, y1=None, xy=None, x2=None, y2=None,
616 bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None,
617 **kwargs):
618 """Generate dipole model with given parameters.
620 It simply defers to `modelObj.makeModel()`, where `modelObj` comes
621 out of `kwargs['modelObj']`.
622 """
623 modelObj = kwargs.pop('modelObj')
624 return modelObj.makeModel(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=fluxNeg,
625 b=b, x1=x1, y1=y1, xy=xy, x2=x2, y2=y2,
626 bNeg=bNeg, x1Neg=x1Neg, y1Neg=y1Neg, xyNeg=xyNeg,
627 x2Neg=x2Neg, y2Neg=y2Neg, **kwargs)
629 dipoleModel = DipoleModel()
631 modelFunctor = dipoleModelFunctor # dipoleModel.makeModel does not work for now.
632 # Create the lmfit model (lmfit uses scipy 'leastsq' option by default - Levenberg-Marquardt)
633 # We have to (later) filter out the nans by hand in our input to gmod.fit().
634 # The only independent variable in the model is "x"; lmfit tries to
635 # introspect variables and parameters from the function signature, but
636 # gets it wrong for the model signature above.
637 gmod = lmfit.Model(modelFunctor, independent_vars=["x"], verbose=verbose)
639 # Add the constraints for centroids, fluxes.
640 # starting constraint - near centroid of footprint
641 fpCentroid = np.array([fp.getCentroid().getX(), fp.getCentroid().getY()])
642 cenNeg = cenPos = fpCentroid
644 pks = fp.getPeaks()
646 if len(pks) >= 1:
647 cenPos = pks[0].getF() # if individual (merged) peaks were detected, use those
648 if len(pks) >= 2: # peaks are already sorted by centroid flux so take the most negative one
649 cenNeg = pks[-1].getF()
651 # For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
652 # First assume dipole is not separated by more than 5*psfSigma.
653 maxSep = self.psfSigma * maxSepInSigma
655 # As an initial guess -- assume the dipole is close to the center of the footprint.
656 if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
657 cenPos = fpCentroid
658 if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
659 cenPos = fpCentroid
661 # parameter hints/constraints: https://lmfit.github.io/lmfit-py/model.html#model-param-hints-section
662 # might make sense to not use bounds -- see http://lmfit.github.io/lmfit-py/bounds.html
663 # also see this discussion -- https://github.com/scipy/scipy/issues/3129
664 gmod.set_param_hint('xcenPos', value=cenPos[0],
665 min=cenPos[0]-maxSep, max=cenPos[0]+maxSep)
666 gmod.set_param_hint('ycenPos', value=cenPos[1],
667 min=cenPos[1]-maxSep, max=cenPos[1]+maxSep)
668 gmod.set_param_hint('xcenNeg', value=cenNeg[0],
669 min=cenNeg[0]-maxSep, max=cenNeg[0]+maxSep)
670 gmod.set_param_hint('ycenNeg', value=cenNeg[1],
671 min=cenNeg[1]-maxSep, max=cenNeg[1]+maxSep)
673 # Use the (flux under the dipole)*5 for an estimate.
674 # Lots of testing showed that having startingFlux be too high was better than too low.
675 startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
676 posFlux = negFlux = startingFlux
678 # TBD: set max. flux limit?
679 gmod.set_param_hint('flux', value=posFlux, min=0.1)
681 if separateNegParams:
682 # TBD: set max negative lobe flux limit?
683 gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
685 # Fixed parameters (don't fit for them if there are no pre-sub images or no gradient fit requested):
686 # Right now (fitBackground == 1), we fit a linear model to the background and then subtract
687 # it from the data and then don't fit the background again (this is faster).
688 # A slower alternative (fitBackground == 2) is to use the estimated background parameters as
689 # starting points in the integrated model fit. That is currently not performed by default,
690 # but might be desirable in some cases.
691 bgParsPos = bgParsNeg = (0., 0., 0.)
692 if ((rel_weight > 0.) and (fitBackground != 0) and (bgGradientOrder >= 0)):
693 pbg = 0.
694 bgFitImage = self.posImage if self.posImage is not None else self.negImage
695 # Fit the gradient to the background (linear model)
696 bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
697 order=bgGradientOrder)
699 # Generate the gradient and subtract it from the pre-subtraction image data
700 if fitBackground == 1:
701 in_x = dipoleModel._generateXYGrid(bbox)
702 pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsPos))
703 z[1, :] -= pbg
704 z[1, :] -= np.nanmedian(z[1, :])
705 posFlux = np.nansum(z[1, :])
706 gmod.set_param_hint('flux', value=posFlux*1.5, min=0.1)
708 if separateNegParams and self.negImage is not None:
709 bgParsNeg = dipoleModel.fitFootprintBackground(source, self.negImage,
710 order=bgGradientOrder)
711 pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsNeg))
712 z[2, :] -= pbg
713 z[2, :] -= np.nanmedian(z[2, :])
714 if separateNegParams:
715 negFlux = np.nansum(z[2, :])
716 gmod.set_param_hint('fluxNeg', value=negFlux*1.5, min=0.1)
718 # Do not subtract the background from the images but include the background parameters in the fit
719 if fitBackground == 2:
720 if bgGradientOrder >= 0:
721 gmod.set_param_hint('b', value=bgParsPos[0])
722 if separateNegParams:
723 gmod.set_param_hint('bNeg', value=bgParsNeg[0])
724 if bgGradientOrder >= 1:
725 gmod.set_param_hint('x1', value=bgParsPos[1])
726 gmod.set_param_hint('y1', value=bgParsPos[2])
727 if separateNegParams:
728 gmod.set_param_hint('x1Neg', value=bgParsNeg[1])
729 gmod.set_param_hint('y1Neg', value=bgParsNeg[2])
730 if bgGradientOrder >= 2:
731 gmod.set_param_hint('xy', value=bgParsPos[3])
732 gmod.set_param_hint('x2', value=bgParsPos[4])
733 gmod.set_param_hint('y2', value=bgParsPos[5])
734 if separateNegParams:
735 gmod.set_param_hint('xyNeg', value=bgParsNeg[3])
736 gmod.set_param_hint('x2Neg', value=bgParsNeg[4])
737 gmod.set_param_hint('y2Neg', value=bgParsNeg[5])
739 y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
740 in_x = np.array([x, y]).astype(np.float64)
741 in_x[0, :] -= in_x[0, :].mean() # center it!
742 in_x[1, :] -= in_x[1, :].mean()
744 # Instead of explicitly using a mask to ignore flagged pixels, just set the ignored pixels'
745 # weights to 0 in the fit. TBD: need to inspect mask planes to set this mask.
746 mask = np.ones_like(z, dtype=bool) # TBD: set mask values to False if the pixels are to be ignored
748 # I'm not sure about the variance planes in the diffim (or convolved pre-sub. images
749 # for that matter) so for now, let's just do an un-weighted least-squares fit
750 # (override weights computed above).
751 weights = mask.astype(np.float64)
752 if self.posImage is not None and rel_weight > 0.:
753 weights = np.array([np.ones_like(diArr), np.ones_like(diArr)*rel_weight,
754 np.ones_like(diArr)*rel_weight])
756 # Set the weights to zero if mask is False
757 if np.any(~mask):
758 weights[~mask] = 0.
760 # Filter out any nans, and make the weights 0.
761 nans = (np.isnan(z) | np.isnan(weights))
762 nNans = nans.sum()
763 if nNans > 0:
764 if nNans < len(z):
765 z[nans] = np.nanmedian(z)
766 else:
767 z[nans] = 0
768 weights[nans] = 0
770 # Note that although we can, we're not required to set initial values for params here,
771 # since we set their param_hint's above.
772 # Can add "method" param to not use 'leastsq' (==levenberg-marquardt), e.g. "method='nelder'"
773 with warnings.catch_warnings():
774 # Ignore lmfit unknown argument warnings:
775 # "psf, rel_weight, footprint, modelObj" all become pass-through kwargs for makeModel.
776 warnings.filterwarnings("ignore", "The keyword argument .* does not match", UserWarning)
777 result = gmod.fit(z, weights=weights, x=in_x, max_nfev=250,
778 method="leastsq", # TODO: try using `least_squares` here for speed/robustness
779 verbose=verbose,
780 # see scipy docs for the meaning of these keywords
781 fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
782 # Our model is float32 internally, so we need a larger epsfcn.
783 'epsfcn': 1e-8},
784 psf=self.diffim.getPsf(), # hereon: kwargs that get passed to makeModel()
785 rel_weight=rel_weight,
786 footprint=fp,
787 modelObj=dipoleModel)
789 if verbose: # the ci_report() seems to fail if neg params are constrained -- TBD why.
790 # Never wanted in production - this takes a long time (longer than the fit!)
791 # This is how to get confidence intervals out:
792 # https://lmfit.github.io/lmfit-py/confidence.html and
793 # http://cars9.uchicago.edu/software/python/lmfit/model.html
794 print(result.fit_report(show_correl=False))
795 if separateNegParams:
796 print(result.ci_report())
798 return result
800 def fitDipole(self, source, tol=1e-7, rel_weight=0.1,
801 fitBackground=1, maxSepInSigma=5., separateNegParams=True,
802 bgGradientOrder=1, verbose=False, display=False):
803 """Fit a dipole model to an input ``diaSource`` (wraps `fitDipoleImpl`).
805 Actually, fits the subimage bounded by the input source's
806 footprint) and optionally constrain the fit using the
807 pre-subtraction images self.posImage (science) and
808 self.negImage (template). Wraps the output into a
809 `pipeBase.Struct` named tuple after computing additional
810 statistics such as orientation and SNR.
812 Parameters
813 ----------
814 source : `lsst.afw.table.SourceRecord`
815 Record containing the (merged) dipole source footprint detected on the diffim
816 tol : `float`, optional
817 Tolerance parameter for scipy.leastsq() optimization
818 rel_weight : `float`, optional
819 Weighting of posImage/negImage relative to the diffim in the fit
820 fitBackground : `int`, {0, 1, 2}, optional
821 How to fit linear background gradient in posImage/negImage
823 - 0: do not fit background at all
824 - 1 (default): pre-fit the background using linear least squares and then do not fit it
825 as part of the dipole fitting optimization
826 - 2: pre-fit the background using linear least squares (as in 1), and use the parameter
827 estimates from that fit as starting parameters for an integrated "re-fit" of the
828 background as part of the overall dipole fitting optimization.
829 maxSepInSigma : `float`, optional
830 Allowed window of centroid parameters relative to peak in input source footprint
831 separateNegParams : `bool`, optional
832 Fit separate parameters to the flux and background gradient in
833 bgGradientOrder : `int`, {0, 1, 2}, optional
834 Desired polynomial order of background gradient
835 verbose: `bool`, optional
836 Be verbose
837 display
838 Display input data, best fit model(s) and residuals in a matplotlib window.
840 Returns
841 -------
842 result : `struct`
843 `pipeBase.Struct` object containing the fit parameters and other information.
845 result : `callable`
846 `lmfit.MinimizerResult` object for debugging and error estimation, etc.
848 Notes
849 -----
850 Parameter `fitBackground` has three options, thus it is an integer:
852 """
854 fitResult = self.fitDipoleImpl(
855 source, tol=tol, rel_weight=rel_weight, fitBackground=fitBackground,
856 maxSepInSigma=maxSepInSigma, separateNegParams=separateNegParams,
857 bgGradientOrder=bgGradientOrder, verbose=verbose)
859 # Display images, model fits and residuals (currently uses matplotlib display functions)
860 if display:
861 fp = source.getFootprint()
862 self.displayFitResults(fp, fitResult)
864 # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit.
865 if fitResult.params['flux'].value <= 1.:
866 self.log.debug("Fitted flux too small for id=%d; ModelResult.message='%s'",
867 source["id"], fitResult.message)
868 return None, fitResult
869 if not fitResult.result.errorbars:
870 self.log.debug("Could not estimate error bars for id=%d; ModelResult.message='%s'",
871 source["id"], fitResult.message)
872 return None, fitResult
874 # TODO: We could include covariances, which could be derived from
875 # `fitResult.params[name].correl`, but those are correlations.
876 posCentroid = measBase.CentroidResult(fitResult.params['xcenPos'].value,
877 fitResult.params['ycenPos'].value,
878 fitResult.params['xcenPos'].stderr,
879 fitResult.params['ycenPos'].stderr)
880 negCentroid = measBase.CentroidResult(fitResult.params['xcenNeg'].value,
881 fitResult.params['ycenNeg'].value,
882 fitResult.params['xcenNeg'].stderr,
883 fitResult.params['ycenNeg'].stderr)
884 xposIdx = fitResult.var_names.index("xcenPos")
885 yposIdx = fitResult.var_names.index("ycenPos")
886 xnegIdx = fitResult.var_names.index("xcenNeg")
887 ynegIdx = fitResult.var_names.index("ycenNeg")
888 centroid = measBase.CentroidResult((fitResult.params['xcenPos'] + fitResult.params['xcenNeg']) / 2,
889 (fitResult.params['ycenPos'] + fitResult.params['ycenNeg']) / 2.,
890 math.sqrt(posCentroid.xErr**2 + negCentroid.xErr**2
891 + 2*fitResult.covar[xposIdx, xnegIdx]) / 2.,
892 math.sqrt(posCentroid.yErr**2 + negCentroid.yErr**2
893 + 2*fitResult.covar[yposIdx, ynegIdx]) / 2.)
894 dx = fitResult.params['xcenPos'].value - fitResult.params['xcenNeg'].value
895 dy = fitResult.params['ycenPos'].value - fitResult.params['ycenNeg'].value
896 angle = np.arctan2(dy, dx)
898 # Extract flux value, compute signalToNoise from flux/variance_within_footprint
899 # Also extract the stderr of flux estimate.
900 # TODO: should this instead use the lmfit-computed uncertainty from
901 # `lmfitResult.result.uvars['flux'].std_dev`?
902 def computeSumVariance(exposure, footprint):
903 return math.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array))
905 # NOTE: These will all be the same unless separateNegParams=True!
906 flux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr)
907 posFlux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr)
908 negFlux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr)
909 if self.posImage is not None:
910 fluxVar = computeSumVariance(self.posImage, source.getFootprint())
911 else:
912 fluxVar = computeSumVariance(self.diffim, source.getFootprint())
913 fluxVarNeg = fluxVar
915 if separateNegParams:
916 negFlux.instFlux = fitResult.params['fluxNeg'].value
917 negFlux.instFluxErr = fitResult.params['fluxNeg'].stderr
918 if self.negImage is not None:
919 fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint())
921 try:
922 signalToNoise = math.sqrt((posFlux.instFlux/fluxVar)**2 + (negFlux.instFlux/fluxVarNeg)**2)
923 except ZeroDivisionError: # catch divide by zero - should never happen.
924 signalToNoise = np.nan
926 out = Struct(posCentroid=posCentroid, negCentroid=negCentroid, centroid=centroid,
927 posFlux=posFlux, negFlux=negFlux, flux=flux, orientation=angle,
928 signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi,
929 nData=fitResult.ndata)
931 # fitResult may be returned for debugging
932 return out, fitResult
934 def displayFitResults(self, footprint, result):
935 """Display data, model fits and residuals (currently uses matplotlib display functions).
937 Parameters
938 ----------
939 footprint : `lsst.afw.detection.Footprint`
940 Footprint containing the dipole that was fit
941 result : `lmfit.MinimizerResult`
942 `lmfit.MinimizerResult` object returned by `lmfit` optimizer
943 """
944 try:
945 import matplotlib.pyplot as plt
946 except ImportError as err:
947 self.log.warning('Unable to import matplotlib: %s', err)
948 raise err
950 def display2dArray(ax, arr, x, y, xErr, yErr, title, extent=None):
951 """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range.
952 """
953 fig = ax.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent)
954 ax.set_title(title)
955 ax.errorbar(x["total"], y["total"], xErr["total"], yErr["total"], c="cyan")
956 ax.errorbar(x["Pos"], y["Pos"], xErr["Pos"], yErr["Pos"], c="green")
957 ax.errorbar(x["Neg"], y["Neg"], xErr["Neg"], yErr["Neg"], c="red")
958 return fig
960 z = result.data
961 fit = result.best_fit
962 bbox = footprint.getBBox()
963 extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY())
965 if z.shape[0] == 3:
966 x, y, xErr, yErr = {}, {}, {}, {}
967 for name in ("Pos", "Neg"):
968 x[name] = result.best_values[f"xcen{name}"]
969 y[name] = result.best_values[f"ycen{name}"]
970 xErr[name] = result.params[f"xcen{name}"].stderr
971 yErr[name] = result.params[f"ycen{name}"].stderr
972 x["total"] = (result.best_values["xcenPos"] + result.best_values["xcenNeg"])/2
973 y["total"] = (result.best_values["ycenPos"] + result.best_values["ycenNeg"])/2
974 xErr["total"] = math.sqrt(xErr["Pos"]**2 + xErr["Neg"]**2)
975 yErr["total"] = math.sqrt(yErr["Pos"]**2 + yErr["Neg"]**2)
977 fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True, figsize=(8, 8))
978 for i, label in enumerate(("total", "Pos", "Neg")):
979 display2dArray(axes[i][0], z[i, :], x, y, xErr, yErr,
980 f'Data {label}', extent=extent)
981 display2dArray(axes[i][1], fit[i, :], x, y, xErr, yErr,
982 f'Model {label}', extent=extent)
983 display2dArray(axes[i][2], z[i, :] - fit[i, :], x, y, xErr, yErr,
984 f'Residual {label}', extent=extent)
986 plt.setp(axes[i][1].get_yticklabels(), visible=False)
987 plt.setp(axes[i][2].get_yticklabels(), visible=False)
988 if i != 2: # remove top two row x-axis labels
989 plt.setp(axes[i][0].get_xticklabels(), visible=False)
990 plt.setp(axes[i][1].get_xticklabels(), visible=False)
991 plt.setp(axes[i][2].get_xticklabels(), visible=False)
992 else:
993 fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True, figsize=(10, 2.5))
994 display2dArray(axes[0], z, 'Data', extent=extent)
995 display2dArray(axes[1], 'Model', extent=extent)
996 display2dArray(axes[2], z - fit, 'Residual', extent=extent)
998 fig.tight_layout(pad=0, w_pad=0, h_pad=0)
999 plt.show()
1002@measBase.register("ip_diffim_DipoleFit")
1003class DipoleFitPlugin(measBase.SingleFramePlugin):
1004 """A single frame measurement plugin that fits dipoles to all merged (two-peak) ``diaSources``.
1006 This measurement plugin accepts up to three input images in
1007 its `measure` method. If these are provided, it includes data
1008 from the pre-subtraction posImage (science image) and optionally
1009 negImage (template image) to constrain the fit. The meat of the
1010 fitting routines are in the class `~lsst.module.name.DipoleFitAlgorithm`.
1012 Notes
1013 -----
1014 The motivation behind this plugin and the necessity for including more than
1015 one exposure are documented in DMTN-007 (http://dmtn-007.lsst.io).
1017 This class is named `ip_diffim_DipoleFit` so that it may be used alongside
1018 the existing `ip_diffim_DipoleMeasurement` classes until such a time as those
1019 are deemed to be replaceable by this.
1020 """
1022 ConfigClass = DipoleFitPluginConfig
1023 DipoleFitAlgorithmClass = DipoleFitAlgorithm # Pointer to the class that performs the fit
1025 FAILURE_EDGE = 1 # too close to the edge
1026 FAILURE_FIT = 2 # failure in the fitting
1027 FAILURE_NOT_DIPOLE = 4 # input source is not a putative dipole to begin with
1028 FAILURE_TOO_LARGE = 8 # input source is too large to be fit
1030 @classmethod
1031 def getExecutionOrder(cls):
1032 """This algorithm simultaneously fits the centroid and flux, and does
1033 not require any previous centroid fit.
1034 """
1035 return cls.CENTROID_ORDER
1037 def __init__(self, config, name, schema, metadata, logName=None):
1038 if logName is None:
1039 logName = name
1040 measBase.SingleFramePlugin.__init__(self, config, name, schema, metadata, logName=logName)
1042 self.log = logging.getLogger(logName)
1044 self._setupSchema(config, name, schema, metadata)
1046 def _setupSchema(self, config, name, schema, metadata):
1047 """Add fields for the outputs, and save the keys for fast assignment.
1048 """
1049 self.posFluxKey = measBase.FluxResultKey.addFields(schema,
1050 schema.join(name, "pos"),
1051 "Dipole positive lobe instrumental flux.")
1052 self.negFluxKey = measBase.FluxResultKey.addFields(schema,
1053 schema.join(name, "neg"),
1054 "Dipole negative lobe instrumental flux.")
1055 doc = "Dipole overall instrumental flux (mean of absolute value of positive and negative lobes)."
1056 self.fluxKey = measBase.FluxResultKey.addFields(schema, name, doc)
1058 self.posCentroidKey = measBase.CentroidResultKey.addFields(schema,
1059 schema.join(name, "pos"),
1060 "Dipole positive lobe centroid position.",
1061 measBase.UncertaintyEnum.SIGMA_ONLY)
1062 self.negCentroidKey = measBase.CentroidResultKey.addFields(schema,
1063 schema.join(name, "neg"),
1064 "Dipole negative lobe centroid position.",
1065 measBase.UncertaintyEnum.SIGMA_ONLY)
1066 self.centroidKey = measBase.CentroidResultKey.addFields(schema,
1067 name,
1068 "Dipole centroid position.",
1069 measBase.UncertaintyEnum.SIGMA_ONLY)
1071 self.orientationKey = schema.addField(
1072 schema.join(name, "orientation"), type=float, units="rad",
1073 doc="Dipole orientation. Convention is CCW from +x on image.")
1075 self.separationKey = schema.addField(
1076 schema.join(name, "separation"), type=float, units="pixel",
1077 doc="Pixel separation between positive and negative lobes of dipole")
1079 self.chi2dofKey = schema.addField(
1080 schema.join(name, "chi2dof"), type=float,
1081 doc="Chi2 per degree of freedom (chi2/(nData-nVariables)) of dipole fit")
1083 self.nDataKey = schema.addField(
1084 schema.join(name, "nData"), type=np.int64,
1085 doc="Number of data points in the dipole fit")
1087 self.signalToNoiseKey = schema.addField(
1088 schema.join(name, "signalToNoise"), type=float,
1089 doc="Estimated signal-to-noise of dipole fit")
1091 self.classificationFlagKey = schema.addField(
1092 schema.join(name, "classification"), type="Flag",
1093 doc="Flag indicating diaSource is classified as a dipole")
1095 self.classificationAttemptedFlagKey = schema.addField(
1096 schema.join(name, "classificationAttempted"), type="Flag",
1097 doc="Flag indicating diaSource was attempted to be classified as a dipole")
1099 self.flagKey = schema.addField(
1100 schema.join(name, "flag"), type="Flag",
1101 doc="General failure flag for dipole fit")
1103 self.edgeFlagKey = schema.addField(
1104 schema.join(name, "flag", "edge"), type="Flag",
1105 doc="Flag set when dipole is too close to edge of image")
1107 def measureDipoles(self, measRecord, exposure, posExp=None, negExp=None):
1108 """Perform the non-linear least squares minimization on the putative dipole source.
1110 Parameters
1111 ----------
1112 measRecord : `lsst.afw.table.SourceRecord`
1113 diaSources that will be measured using dipole measurement
1114 exposure : `lsst.afw.image.Exposure`
1115 Difference exposure on which the diaSources were detected; `exposure = posExp-negExp`
1116 If both `posExp` and `negExp` are `None`, will attempt to fit the
1117 dipole to just the `exposure` with no constraint.
1118 posExp : `lsst.afw.image.Exposure`, optional
1119 "Positive" exposure, typically a science exposure, or None if unavailable
1120 When `posExp` is `None`, will compute `posImage = exposure + negExp`.
1121 negExp : `lsst.afw.image.Exposure`, optional
1122 "Negative" exposure, typically a template exposure, or None if unavailable
1123 When `negExp` is `None`, will compute `negImage = posExp - exposure`.
1125 Notes
1126 -----
1127 The main functionality of this routine was placed outside of
1128 this plugin (into `DipoleFitAlgorithm.fitDipole()`) so that
1129 `DipoleFitAlgorithm.fitDipole()` can be called separately for
1130 testing (@see `tests/testDipoleFitter.py`)
1131 """
1132 result = None
1133 pks = measRecord.getFootprint().getPeaks()
1135 # Check if the footprint consists of a putative dipole - else don't fit it.
1136 if (
1137 # One peak in the footprint (not a dipole)
1138 ((nPeaks := len(pks)) <= 1)
1139 # Peaks are the same sign (not a dipole); peaks are ordered
1140 # from highest to lowest.
1141 or (nPeaks > 1 and (np.sign(pks[0].getPeakValue())
1142 == np.sign(pks[-1].getPeakValue())))
1143 ):
1144 if not self.config.fitAllDiaSources:
1145 # Non-dipoles fall back on the centroid slot for positions,
1146 # errors, and the failure flag, if we're not fitting them.
1147 measRecord[self.centroidKey.getX()] = measRecord.getX()
1148 measRecord[self.centroidKey.getY()] = measRecord.getY()
1149 self.centroidKey.getCentroidErr().set(measRecord, measRecord.getCentroidErr())
1150 measRecord[self.flagKey] = measRecord.getCentroidFlag()
1151 return
1153 # Footprint is too large (not a dipole).
1154 if ((area := measRecord.getFootprint().getArea()) > self.config.maxFootprintArea):
1155 self.fail(measRecord, measBase.MeasurementError(f"{area} > {self.config.maxFootprintArea}",
1156 self.FAILURE_TOO_LARGE))
1157 return
1159 try:
1160 alg = self.DipoleFitAlgorithmClass(exposure, posImage=posExp, negImage=negExp)
1161 result, _ = alg.fitDipole(
1162 measRecord, rel_weight=self.config.relWeight,
1163 tol=self.config.tolerance,
1164 maxSepInSigma=self.config.maxSeparation,
1165 fitBackground=self.config.fitBackground,
1166 separateNegParams=self.config.fitSeparateNegParams,
1167 verbose=False, display=False)
1168 except pexExcept.LengthError:
1169 self.fail(measRecord, measBase.MeasurementError('edge failure', self.FAILURE_EDGE))
1170 except Exception as e:
1171 errorMessage = f"Exception in dipole fit. {e.__class__.__name__}: {e}"
1172 self.fail(measRecord, measBase.MeasurementError(errorMessage, self.FAILURE_FIT))
1174 self.log.debug("Dipole fit result: %d %s", measRecord.getId(), str(result))
1176 if result is None:
1177 self.fail(measRecord, measBase.MeasurementError("bad dipole fit", self.FAILURE_FIT))
1178 return
1180 # TODO: add chi2, dipole classification
1181 self.posFluxKey.set(measRecord, result.posFlux)
1182 self.posCentroidKey.set(measRecord, result.posCentroid)
1184 self.negFluxKey.set(measRecord, result.negFlux)
1185 self.negCentroidKey.set(measRecord, result.negCentroid)
1187 self.fluxKey.set(measRecord, result.flux)
1188 self.centroidKey.set(measRecord, result.centroid)
1190 measRecord[self.orientationKey] = result.orientation
1191 measRecord[self.separationKey] = math.sqrt((result.posCentroid.x - result.negCentroid.x)**2
1192 + (result.posCentroid.y - result.negCentroid.y)**2)
1194 measRecord[self.signalToNoiseKey] = result.signalToNoise
1195 measRecord[self.chi2dofKey] = result.redChi2
1197 if result.nData >= 1:
1198 measRecord[self.nDataKey] = result.nData
1199 else:
1200 measRecord[self.nDataKey] = 0
1202 self.doClassify(measRecord, result.chi2)
1204 def doClassify(self, measRecord, chi2val):
1205 """Classify a source as a dipole.
1207 Parameters
1208 ----------
1209 measRecord : TODO: DM-17458
1210 TODO: DM-17458
1211 chi2val : TODO: DM-17458
1212 TODO: DM-17458
1214 Notes
1215 -----
1216 Sources are classified as dipoles, or not, according to three criteria:
1218 1. Does the total signal-to-noise surpass the ``minSn``?
1219 2. Are the pos/neg fluxes greater than 1.0 and no more than 0.65 (``maxFluxRatio``)
1220 of the total flux? By default this will never happen since ``posFlux == negFlux``.
1221 3. Is it a good fit (``chi2dof`` < 1)? (Currently not used.)
1222 """
1224 # First, does the total signal-to-noise surpass the minSn?
1225 passesSn = measRecord[self.signalToNoiseKey] > self.config.minSn
1227 # Second, are the pos/neg fluxes greater than 1.0 and no more than 0.65 (param maxFluxRatio)
1228 # of the total flux? By default this will never happen since posFlux = negFlux.
1229 passesFluxPos = (abs(measRecord[self.posFluxKey.getInstFlux()])
1230 / (measRecord[self.fluxKey.getInstFlux()]*2.)) < self.config.maxFluxRatio
1231 passesFluxPos &= (abs(measRecord[self.posFluxKey.getInstFlux()]) >= 1.0)
1232 passesFluxNeg = (abs(measRecord[self.negFluxKey.getInstFlux()])
1233 / (measRecord[self.fluxKey.getInstFlux()]*2.)) < self.config.maxFluxRatio
1234 passesFluxNeg &= (abs(measRecord[self.negFluxKey.getInstFlux()]) >= 1.0)
1235 allPass = (passesSn and passesFluxPos and passesFluxNeg) # and passesChi2)
1237 # Third, is it a good fit (chi2dof < 1)?
1238 # Use scipy's chi2 cumulative distrib to estimate significance
1239 # This doesn't really work since I don't trust the values in the variance plane (which
1240 # affects the least-sq weights, which affects the resulting chi2).
1241 # But I'm going to keep this here for future use.
1242 if False:
1243 from scipy.stats import chi2
1244 ndof = chi2val / measRecord[self.chi2dofKey]
1245 significance = chi2.cdf(chi2val, ndof)
1246 passesChi2 = significance < self.config.maxChi2DoF
1247 allPass = allPass and passesChi2
1249 measRecord.set(self.classificationAttemptedFlagKey, True)
1251 if allPass: # Note cannot pass `allPass` into the `measRecord.set()` call below...?
1252 measRecord.set(self.classificationFlagKey, True)
1253 else:
1254 measRecord.set(self.classificationFlagKey, False)
1256 def fail(self, measRecord, error=None):
1257 """Catch failures and set the correct flags.
1259 Fallback on the current slot centroid positions, but set the dipole
1260 failure flag, since we attempted to fit the source.
1261 """
1262 measRecord[self.centroidKey.getX()] = measRecord.getX()
1263 measRecord[self.centroidKey.getY()] = measRecord.getY()
1264 self.centroidKey.getCentroidErr().set(measRecord, measRecord.getCentroidErr())
1266 measRecord.set(self.flagKey, True)
1267 if error is not None:
1268 if error.getFlagBit() == self.FAILURE_EDGE:
1269 self.log.debug('DipoleFitPlugin not run on record %d: %s', measRecord.getId(), str(error))
1270 measRecord.set(self.edgeFlagKey, True)
1271 if error.getFlagBit() == self.FAILURE_FIT:
1272 self.log.warning('DipoleFitPlugin failed on record %d: %s', measRecord.getId(), str(error))
1273 if error.getFlagBit() == self.FAILURE_TOO_LARGE:
1274 self.log.debug('DipoleFitPlugin not run on record with too large footprint %d: %s',
1275 measRecord.getId(), str(error))
1276 else:
1277 self.log.warning('DipoleFitPlugin failed on record %d', measRecord.getId())