LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
dipoleFitTask.py
Go to the documentation of this file.
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#
22
23import logging
24import numpy as np
25import warnings
26
27import lsst.afw.image as afwImage
28import lsst.meas.base as measBase
29import lsst.afw.detection as afwDet
30import lsst.geom as geom
31import lsst.pex.exceptions as pexExcept
32import lsst.pex.config as pexConfig
33from lsst.pipe.base import Struct
34from lsst.utils.timer import timeMethod
35
36__all__ = ("DipoleFitTask", "DipoleFitPlugin", "DipoleFitTaskConfig", "DipoleFitPluginConfig",
37 "DipoleFitAlgorithm")
38
39
40# Create a new measurement task (`DipoleFitTask`) that can handle all other SFM tasks but can
41# pass a separate pos- and neg- exposure/image to the `DipoleFitPlugin`s `run()` method.
42
43
44class DipoleFitPluginConfig(measBase.SingleFramePluginConfig):
45 """Configuration for DipoleFitPlugin
46 """
47
48 fitAllDiaSources = pexConfig.Field(
49 dtype=bool, default=False,
50 doc="""Attempte dipole fit of all diaSources (otherwise just the ones consisting of overlapping
51 positive and negative footprints)""")
52
53 maxSeparation = pexConfig.Field(
54 dtype=float, default=5.,
55 doc="Assume dipole is not separated by more than maxSeparation * psfSigma")
56
57 relWeight = pexConfig.Field(
58 dtype=float, default=0.5,
59 doc="""Relative weighting of pre-subtraction images (higher -> greater influence of pre-sub.
60 images on fit)""")
61
62 tolerance = pexConfig.Field(
63 dtype=float, default=1e-7,
64 doc="Fit tolerance")
65
66 fitBackground = pexConfig.Field(
67 dtype=int, default=1,
68 doc="Set whether and how to fit for linear gradient in pre-sub. images. Possible values:"
69 "0: do not fit background at all"
70 "1 (default): pre-fit the background using linear least squares and then do not fit it as part"
71 "of the dipole fitting optimization"
72 "2: pre-fit the background using linear least squares (as in 1), and use the parameter"
73 "estimates from that fit as starting parameters for an integrated re-fit of the background")
74
75 fitSeparateNegParams = pexConfig.Field(
76 dtype=bool, default=False,
77 doc="Include parameters to fit for negative values (flux, gradient) separately from pos.")
78
79 # Config params for classification of detected diaSources as dipole or not
80 minSn = pexConfig.Field(
81 dtype=float, default=np.sqrt(2) * 5.0,
82 doc="Minimum quadrature sum of positive+negative lobe S/N to be considered a dipole")
83
84 maxFluxRatio = pexConfig.Field(
85 dtype=float, default=0.65,
86 doc="Maximum flux ratio in either lobe to be considered a dipole")
87
88 maxChi2DoF = pexConfig.Field(
89 dtype=float, default=0.05,
90 doc="""Maximum Chi2/DoF significance of fit to be considered a dipole.
91 Default value means \"Choose a chi2DoF corresponding to a significance level of at most 0.05\"
92 (note this is actually a significance, not a chi2 value).""")
93
94 maxFootprintArea = pexConfig.Field(
95 dtype=int, default=1_200,
96 doc=("Maximum area for footprints before they are ignored as large; "
97 "non-positive means no threshold applied"
98 "Threshold chosen for HSC and DECam data, see DM-38741 for details."))
99
100
101class DipoleFitTaskConfig(measBase.SingleFrameMeasurementConfig):
102
103 def setDefaults(self):
104 measBase.SingleFrameMeasurementConfig.setDefaults(self)
105
106 self.plugins.names = ["base_SdssCentroid",
107 "ip_diffim_DipoleFit",
108 "base_CircularApertureFlux",
109 "base_PixelFlags",
110 "base_SkyCoord",
111 "base_PsfFlux",
112 "base_SdssShape",
113 "base_GaussianFlux",
114 ]
115 # Only measure the apertures we need to report in the alert stream.
116 self.plugins["base_CircularApertureFlux"].radii = [12.0]
117
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"
125
126
127class DipoleFitTask(measBase.SingleFrameMeasurementTask):
128 """A task that fits a dipole to a difference image, with an optional
129 separate detection image.
130
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 """
135
136 ConfigClass = DipoleFitTaskConfig
137 _DefaultName = "dipoleFit"
138
139 def __init__(self, schema, algMetadata=None, **kwargs):
140 super().__init__(schema, algMetadata, **kwargs)
141
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
151
152 @timeMethod
153 def run(self, sources, exposure, posExp=None, negExp=None, **kwargs):
154 """Run dipole measurement and classification.
155
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.
159
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.
180 super().run(sources, exposure, **kwargs)
181
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")
187
188 self.plugins = self.plugins_post
189 super().run(sources, exposure, **kwargs)
190
191
193 """Lightweight class containing methods for generating a dipole model for fitting
194 to sources in diffims, used by DipoleFitAlgorithm.
195
196 See also:
197 `DMTN-007: Dipole characterization for image differencing <https://dmtn-007.lsst.io>`_.
198 """
199
200 def __init__(self):
201 import lsstDebug
202 self.debug = lsstDebug.Info(__name__).debug
203 self.log = logging.getLogger(__name__)
204
205 def makeBackgroundModel(self, in_x, pars=None):
206 """Generate gradient model (2-d array) with up to 2nd-order polynomial
207
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).
221
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 """
228
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
232
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)
245
246 return gradient
247
248 def _generateXYGrid(self, bbox):
249 """Generate a meshgrid covering the x,y coordinates bounded by bbox
250
251 Parameters
252 ----------
253 bbox : `lsst.geom.Box2I`
254 input Bounding Box defining the coordinate limits
255
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 """
262
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
268
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`.
272
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
281
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)
290
291 subim2 = afwImage.ImageF(bbox, badfill)
292 fp.getSpans().unflatten(subim2.array, fp.getImageArray(), bbox.getCorners()[0])
293 return subim2
294
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.
297
298 Only fit the pixels OUTSIDE of the footprint, but within its bounding box.
299
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.
308
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 """
315
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)
322
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)
328
329 isBg = np.isnan(posFpImg.array).ravel()
330
331 data = posImg.array.ravel()
332 data = data[isBg]
333 B = data
334
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)
343
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
349
350 pars = np.linalg.lstsq(M, B, rcond=-1)[0]
351 return pars
352
353 def makeStarModel(self, bbox, psf, xcen, ycen, flux):
354 """Generate a 2D image model of a single PDF centered at the given coordinates.
355
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'
368
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 """
375
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)
380
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)
385
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
393
394 return p_Im
395
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.
401
402 This is the function whose sum-of-squared difference from data
403 is minimized by `lmfit`.
404
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.
421
422 **kwargs : `dict` [`str`]
423 Keyword arguments passed through ``lmfit`` and
424 used by this function. These must include:
425
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
430
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 """
441
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()
446
447 if fluxNeg is None:
448 fluxNeg = flux
449
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)
456
457 posIm = self.makeStarModel(bbox, psf, xcenPos, ycenPos, flux)
458 negIm = self.makeStarModel(bbox, psf, xcenNeg, ycenNeg, fluxNeg)
459
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()
466
467 if b is not None:
468 gradient = self.makeBackgroundModel(in_x, (b, x1, y1, xy, x2, y2))
469
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
475
476 posIm.array[:, :] += gradient
477 negIm.array[:, :] += gradientNeg
478
479 # Generate the diffIm model
480 diffIm = afwImage.ImageF(bbox)
481 diffIm += posIm
482 diffIm -= negIm
483
484 zout = diffIm.array
485 if rel_weight > 0.:
486 zout = np.append([zout], [posIm.array, negIm.array], axis=0)
487
488 return zout
489
490
492 """Fit a dipole model using an image difference.
493
494 See also:
495 `DMTN-007: Dipole characterization for image differencing <https://dmtn-007.lsst.io>`_.
496 """
497
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'
501
502 # Below is a (somewhat incomplete) list of improvements
503 # that would be worth investigating, given the time:
504
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?)
516
517 def __init__(self, diffim, posImage=None, negImage=None):
518 """Algorithm to run dipole measurement on a diaSource
519
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 """
529
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()
538
539 self.log = logging.getLogger(__name__)
540
541 import lsstDebug
542 self.debug = lsstDebug.Info(__name__).debug
543
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.
548
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.
552
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
571
572 Returns
573 -------
574 result : `lmfit.MinimizerResult`
575 return `lmfit.MinimizerResult` object containing the fit
576 parameters and other information.
577 """
578
579 # Only import lmfit if someone wants to use the new DipoleFitAlgorithm.
580 import lmfit
581
582 fp = source.getFootprint()
583 bbox = fp.getBBox()
584 subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)
585
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)
590
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
602
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"
610
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.
619
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)
628
629 dipoleModel = DipoleModel()
630
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 gmod = lmfit.Model(modelFunctor, verbose=verbose)
635
636 # Add the constraints for centroids, fluxes.
637 # starting constraint - near centroid of footprint
638 fpCentroid = np.array([fp.getCentroid().getX(), fp.getCentroid().getY()])
639 cenNeg = cenPos = fpCentroid
640
641 pks = fp.getPeaks()
642
643 if len(pks) >= 1:
644 cenPos = pks[0].getF() # if individual (merged) peaks were detected, use those
645 if len(pks) >= 2: # peaks are already sorted by centroid flux so take the most negative one
646 cenNeg = pks[-1].getF()
647
648 # For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
649 # First assume dipole is not separated by more than 5*psfSigma.
650 maxSep = self.psfSigma * maxSepInSigma
651
652 # As an initial guess -- assume the dipole is close to the center of the footprint.
653 if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
654 cenPos = fpCentroid
655 if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
656 cenPos = fpCentroid
657
658 # parameter hints/constraints: https://lmfit.github.io/lmfit-py/model.html#model-param-hints-section
659 # might make sense to not use bounds -- see http://lmfit.github.io/lmfit-py/bounds.html
660 # also see this discussion -- https://github.com/scipy/scipy/issues/3129
661 gmod.set_param_hint('xcenPos', value=cenPos[0],
662 min=cenPos[0]-maxSep, max=cenPos[0]+maxSep)
663 gmod.set_param_hint('ycenPos', value=cenPos[1],
664 min=cenPos[1]-maxSep, max=cenPos[1]+maxSep)
665 gmod.set_param_hint('xcenNeg', value=cenNeg[0],
666 min=cenNeg[0]-maxSep, max=cenNeg[0]+maxSep)
667 gmod.set_param_hint('ycenNeg', value=cenNeg[1],
668 min=cenNeg[1]-maxSep, max=cenNeg[1]+maxSep)
669
670 # Use the (flux under the dipole)*5 for an estimate.
671 # Lots of testing showed that having startingFlux be too high was better than too low.
672 startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
673 posFlux = negFlux = startingFlux
674
675 # TBD: set max. flux limit?
676 gmod.set_param_hint('flux', value=posFlux, min=0.1)
677
678 if separateNegParams:
679 # TBD: set max negative lobe flux limit?
680 gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
681
682 # Fixed parameters (don't fit for them if there are no pre-sub images or no gradient fit requested):
683 # Right now (fitBackground == 1), we fit a linear model to the background and then subtract
684 # it from the data and then don't fit the background again (this is faster).
685 # A slower alternative (fitBackground == 2) is to use the estimated background parameters as
686 # starting points in the integrated model fit. That is currently not performed by default,
687 # but might be desirable in some cases.
688 bgParsPos = bgParsNeg = (0., 0., 0.)
689 if ((rel_weight > 0.) and (fitBackground != 0) and (bgGradientOrder >= 0)):
690 pbg = 0.
691 bgFitImage = self.posImage if self.posImage is not None else self.negImage
692 # Fit the gradient to the background (linear model)
693 bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
694 order=bgGradientOrder)
695
696 # Generate the gradient and subtract it from the pre-subtraction image data
697 if fitBackground == 1:
698 in_x = dipoleModel._generateXYGrid(bbox)
699 pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsPos))
700 z[1, :] -= pbg
701 z[1, :] -= np.nanmedian(z[1, :])
702 posFlux = np.nansum(z[1, :])
703 gmod.set_param_hint('flux', value=posFlux*1.5, min=0.1)
704
705 if separateNegParams and self.negImage is not None:
706 bgParsNeg = dipoleModel.fitFootprintBackground(source, self.negImage,
707 order=bgGradientOrder)
708 pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsNeg))
709 z[2, :] -= pbg
710 z[2, :] -= np.nanmedian(z[2, :])
711 if separateNegParams:
712 negFlux = np.nansum(z[2, :])
713 gmod.set_param_hint('fluxNeg', value=negFlux*1.5, min=0.1)
714
715 # Do not subtract the background from the images but include the background parameters in the fit
716 if fitBackground == 2:
717 if bgGradientOrder >= 0:
718 gmod.set_param_hint('b', value=bgParsPos[0])
719 if separateNegParams:
720 gmod.set_param_hint('bNeg', value=bgParsNeg[0])
721 if bgGradientOrder >= 1:
722 gmod.set_param_hint('x1', value=bgParsPos[1])
723 gmod.set_param_hint('y1', value=bgParsPos[2])
724 if separateNegParams:
725 gmod.set_param_hint('x1Neg', value=bgParsNeg[1])
726 gmod.set_param_hint('y1Neg', value=bgParsNeg[2])
727 if bgGradientOrder >= 2:
728 gmod.set_param_hint('xy', value=bgParsPos[3])
729 gmod.set_param_hint('x2', value=bgParsPos[4])
730 gmod.set_param_hint('y2', value=bgParsPos[5])
731 if separateNegParams:
732 gmod.set_param_hint('xyNeg', value=bgParsNeg[3])
733 gmod.set_param_hint('x2Neg', value=bgParsNeg[4])
734 gmod.set_param_hint('y2Neg', value=bgParsNeg[5])
735
736 y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
737 in_x = np.array([x, y]).astype(np.float64)
738 in_x[0, :] -= in_x[0, :].mean() # center it!
739 in_x[1, :] -= in_x[1, :].mean()
740
741 # Instead of explicitly using a mask to ignore flagged pixels, just set the ignored pixels'
742 # weights to 0 in the fit. TBD: need to inspect mask planes to set this mask.
743 mask = np.ones_like(z, dtype=bool) # TBD: set mask values to False if the pixels are to be ignored
744
745 # I'm not sure about the variance planes in the diffim (or convolved pre-sub. images
746 # for that matter) so for now, let's just do an un-weighted least-squares fit
747 # (override weights computed above).
748 weights = mask.astype(np.float64)
749 if self.posImage is not None and rel_weight > 0.:
750 weights = np.array([np.ones_like(diArr), np.ones_like(diArr)*rel_weight,
751 np.ones_like(diArr)*rel_weight])
752
753 # Set the weights to zero if mask is False
754 if np.any(~mask):
755 weights[~mask] = 0.
756
757 # Filter out any nans, and make the weights 0.
758 nans = (np.isnan(z) | np.isnan(weights))
759 nNans = nans.sum()
760 if nNans > 0:
761 if nNans < len(z):
762 z[nans] = np.nanmedian(z)
763 else:
764 z[nans] = 0
765 weights[nans] = 0
766
767 # Note that although we can, we're not required to set initial values for params here,
768 # since we set their param_hint's above.
769 # Can add "method" param to not use 'leastsq' (==levenberg-marquardt), e.g. "method='nelder'"
770 with warnings.catch_warnings():
771 # Ignore lmfit unknown argument warnings:
772 # "psf, rel_weight, footprint, modelObj" all become pass-through kwargs for makeModel.
773 warnings.filterwarnings("ignore", "The keyword argument .* does not match", UserWarning)
774 result = gmod.fit(z, weights=weights, x=in_x, max_nfev=250,
775 method="leastsq", # TODO: try using `least_squares` here for speed/robustness
776 verbose=verbose,
777 # see scipy docs for the meaning of these keywords
778 fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
779 # Our model is float32 internally, so we need a larger epsfcn.
780 'epsfcn': 1e-8},
781 psf=self.diffim.getPsf(), # hereon: kwargs that get passed to makeModel()
782 rel_weight=rel_weight,
783 footprint=fp,
784 modelObj=dipoleModel)
785
786 if verbose: # the ci_report() seems to fail if neg params are constrained -- TBD why.
787 # Never wanted in production - this takes a long time (longer than the fit!)
788 # This is how to get confidence intervals out:
789 # https://lmfit.github.io/lmfit-py/confidence.html and
790 # http://cars9.uchicago.edu/software/python/lmfit/model.html
791 print(result.fit_report(show_correl=False))
792 if separateNegParams:
793 print(result.ci_report())
794
795 return result
796
797 def fitDipole(self, source, tol=1e-7, rel_weight=0.1,
798 fitBackground=1, maxSepInSigma=5., separateNegParams=True,
799 bgGradientOrder=1, verbose=False, display=False):
800 """Fit a dipole model to an input ``diaSource`` (wraps `fitDipoleImpl`).
801
802 Actually, fits the subimage bounded by the input source's
803 footprint) and optionally constrain the fit using the
804 pre-subtraction images self.posImage (science) and
805 self.negImage (template). Wraps the output into a
806 `pipeBase.Struct` named tuple after computing additional
807 statistics such as orientation and SNR.
808
809 Parameters
810 ----------
811 source : `lsst.afw.table.SourceRecord`
812 Record containing the (merged) dipole source footprint detected on the diffim
813 tol : `float`, optional
814 Tolerance parameter for scipy.leastsq() optimization
815 rel_weight : `float`, optional
816 Weighting of posImage/negImage relative to the diffim in the fit
817 fitBackground : `int`, {0, 1, 2}, optional
818 How to fit linear background gradient in posImage/negImage
819
820 - 0: do not fit background at all
821 - 1 (default): pre-fit the background using linear least squares and then do not fit it
822 as part of the dipole fitting optimization
823 - 2: pre-fit the background using linear least squares (as in 1), and use the parameter
824 estimates from that fit as starting parameters for an integrated "re-fit" of the
825 background as part of the overall dipole fitting optimization.
826 maxSepInSigma : `float`, optional
827 Allowed window of centroid parameters relative to peak in input source footprint
828 separateNegParams : `bool`, optional
829 Fit separate parameters to the flux and background gradient in
830 bgGradientOrder : `int`, {0, 1, 2}, optional
831 Desired polynomial order of background gradient
832 verbose: `bool`, optional
833 Be verbose
834 display
835 Display input data, best fit model(s) and residuals in a matplotlib window.
836
837 Returns
838 -------
839 result : `struct`
840 `pipeBase.Struct` object containing the fit parameters and other information.
841
842 result : `callable`
843 `lmfit.MinimizerResult` object for debugging and error estimation, etc.
844
845 Notes
846 -----
847 Parameter `fitBackground` has three options, thus it is an integer:
848
849 """
850
851 fitResult = self.fitDipoleImpl(
852 source, tol=tol, rel_weight=rel_weight, fitBackground=fitBackground,
853 maxSepInSigma=maxSepInSigma, separateNegParams=separateNegParams,
854 bgGradientOrder=bgGradientOrder, verbose=verbose)
855
856 # Display images, model fits and residuals (currently uses matplotlib display functions)
857 if display:
858 fp = source.getFootprint()
859 self.displayFitResults(fp, fitResult)
860
861 fitParams = fitResult.best_values
862 if fitParams['flux'] <= 1.: # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit.
863 return None, fitResult
864
865 centroid = ((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2.,
866 (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2.)
867 dx, dy = fitParams['xcenPos'] - fitParams['xcenNeg'], fitParams['ycenPos'] - fitParams['ycenNeg']
868 angle = np.arctan2(dy, dx) / np.pi * 180. # convert to degrees (should keep as rad?)
869
870 # Exctract flux value, compute signalToNoise from flux/variance_within_footprint
871 # Also extract the stderr of flux estimate.
872 def computeSumVariance(exposure, footprint):
873 return np.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array))
874
875 fluxVal = fluxVar = fitParams['flux']
876 fluxErr = fluxErrNeg = fitResult.params['flux'].stderr
877 if self.posImage is not None:
878 fluxVar = computeSumVariance(self.posImage, source.getFootprint())
879 else:
880 fluxVar = computeSumVariance(self.diffim, source.getFootprint())
881
882 fluxValNeg, fluxVarNeg = fluxVal, fluxVar
883 if separateNegParams:
884 fluxValNeg = fitParams['fluxNeg']
885 fluxErrNeg = fitResult.params['fluxNeg'].stderr
886 if self.negImage is not None:
887 fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint())
888
889 try:
890 signalToNoise = np.sqrt((fluxVal/fluxVar)**2 + (fluxValNeg/fluxVarNeg)**2)
891 except ZeroDivisionError: # catch divide by zero - should never happen.
892 signalToNoise = np.nan
893
894 out = Struct(posCentroidX=fitParams['xcenPos'], posCentroidY=fitParams['ycenPos'],
895 negCentroidX=fitParams['xcenNeg'], negCentroidY=fitParams['ycenNeg'],
896 posFlux=fluxVal, negFlux=-fluxValNeg, posFluxErr=fluxErr, negFluxErr=fluxErrNeg,
897 centroidX=centroid[0], centroidY=centroid[1], orientation=angle,
898 signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi,
899 nData=fitResult.ndata)
900
901 # fitResult may be returned for debugging
902 return out, fitResult
903
904 def displayFitResults(self, footprint, result):
905 """Display data, model fits and residuals (currently uses matplotlib display functions).
906
907 Parameters
908 ----------
909 footprint : `lsst.afw.detection.Footprint`
910 Footprint containing the dipole that was fit
911 result : `lmfit.MinimizerResult`
912 `lmfit.MinimizerResult` object returned by `lmfit` optimizer
913 """
914 try:
915 import matplotlib.pyplot as plt
916 except ImportError as err:
917 self.log.warning('Unable to import matplotlib: %s', err)
918 raise err
919
920 def display2dArray(arr, title='Data', extent=None):
921 """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range.
922 """
923 fig = plt.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent)
924 plt.title(title)
925 plt.colorbar(fig, cmap='gray')
926 return fig
927
928 z = result.data
929 fit = result.best_fit
930 bbox = footprint.getBBox()
931 extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY())
932 if z.shape[0] == 3:
933 plt.figure(figsize=(8, 8))
934 for i in range(3):
935 plt.subplot(3, 3, i*3+1)
936 display2dArray(z[i, :], 'Data', extent=extent)
937 plt.subplot(3, 3, i*3+2)
938 display2dArray(fit[i, :], 'Model', extent=extent)
939 plt.subplot(3, 3, i*3+3)
940 display2dArray(z[i, :] - fit[i, :], 'Residual', extent=extent)
941 else:
942 plt.figure(figsize=(8, 2.5))
943 plt.subplot(1, 3, 1)
944 display2dArray(z, 'Data', extent=extent)
945 plt.subplot(1, 3, 2)
946 display2dArray(fit, 'Model', extent=extent)
947 plt.subplot(1, 3, 3)
948 display2dArray(z - fit, 'Residual', extent=extent)
949
950 plt.show()
951
952
953@measBase.register("ip_diffim_DipoleFit")
954class DipoleFitPlugin(measBase.SingleFramePlugin):
955 """A single frame measurement plugin that fits dipoles to all merged (two-peak) ``diaSources``.
956
957 This measurement plugin accepts up to three input images in
958 its `measure` method. If these are provided, it includes data
959 from the pre-subtraction posImage (science image) and optionally
960 negImage (template image) to constrain the fit. The meat of the
961 fitting routines are in the class `~lsst.module.name.DipoleFitAlgorithm`.
962
963 Notes
964 -----
965 The motivation behind this plugin and the necessity for including more than
966 one exposure are documented in DMTN-007 (http://dmtn-007.lsst.io).
967
968 This class is named `ip_diffim_DipoleFit` so that it may be used alongside
969 the existing `ip_diffim_DipoleMeasurement` classes until such a time as those
970 are deemed to be replaceable by this.
971 """
972
973 ConfigClass = DipoleFitPluginConfig
974 DipoleFitAlgorithmClass = DipoleFitAlgorithm # Pointer to the class that performs the fit
975
976 FAILURE_EDGE = 1 # too close to the edge
977 FAILURE_FIT = 2 # failure in the fitting
978 FAILURE_NOT_DIPOLE = 4 # input source is not a putative dipole to begin with
979 FAILURE_TOO_LARGE = 8 # input source is too large to be fit
980
981 @classmethod
983 """This algorithm simultaneously fits the centroid and flux, and does
984 not require any previous centroid fit.
985 """
986 return cls.CENTROID_ORDER
987
988 def __init__(self, config, name, schema, metadata, logName=None):
989 if logName is None:
990 logName = name
991 measBase.SingleFramePlugin.__init__(self, config, name, schema, metadata, logName=logName)
992
993 self.log = logging.getLogger(logName)
994
995 self._setupSchema(config, name, schema, metadata)
996
997 def _setupSchema(self, config, name, schema, metadata):
998 """Add fields for the outputs, and save the keys for fast assignment.
999 """
1000 self.posFluxKey = measBase.FluxResultKey.addFields(schema,
1001 schema.join(name, "pos"),
1002 "Dipole positive lobe instrumental flux.")
1003 self.negFluxKey = measBase.FluxResultKey.addFields(schema,
1004 schema.join(name, "neg"),
1005 "Dipole negative lobe instrumental flux.")
1006 doc = "Dipole overall instrumental flux (mean of absolute value of positive and negative lobes)."
1007 self.fluxKey = measBase.FluxResultKey.addFields(schema, name, doc)
1008
1009 self.posCentroidKey = measBase.CentroidResultKey.addFields(schema,
1010 schema.join(name, "pos"),
1011 "Dipole positive lobe centroid position.",
1012 measBase.UncertaintyEnum.NO_UNCERTAINTY)
1013 self.negCentroidKey = measBase.CentroidResultKey.addFields(schema,
1014 schema.join(name, "neg"),
1015 "Dipole negative lobe centroid position.",
1016 measBase.UncertaintyEnum.NO_UNCERTAINTY)
1017 self.centroidKey = measBase.CentroidResultKey.addFields(schema,
1018 name,
1019 "Dipole centroid position.",
1020 measBase.UncertaintyEnum.SIGMA_ONLY)
1021
1022 self.orientationKey = schema.addField(
1023 schema.join(name, "orientation"), type=float, units="deg",
1024 doc="Dipole orientation")
1025
1026 self.separationKey = schema.addField(
1027 schema.join(name, "separation"), type=float, units="pixel",
1028 doc="Pixel separation between positive and negative lobes of dipole")
1029
1030 self.chi2dofKey = schema.addField(
1031 schema.join(name, "chi2dof"), type=float,
1032 doc="Chi2 per degree of freedom (chi2/(nData-nVariables)) of dipole fit")
1033
1034 self.nDataKey = schema.addField(
1035 schema.join(name, "nData"), type=np.int64,
1036 doc="Number of data points in the dipole fit")
1037
1038 self.signalToNoiseKey = schema.addField(
1039 schema.join(name, "signalToNoise"), type=float,
1040 doc="Estimated signal-to-noise of dipole fit")
1041
1042 self.classificationFlagKey = schema.addField(
1043 schema.join(name, "flag", "classification"), type="Flag",
1044 doc="Flag indicating diaSource is classified as a dipole")
1045
1046 self.classificationAttemptedFlagKey = schema.addField(
1047 schema.join(name, "flag", "classificationAttempted"), type="Flag",
1048 doc="Flag indicating diaSource was attempted to be classified as a dipole")
1049
1050 self.flagKey = schema.addField(
1051 schema.join(name, "flag"), type="Flag",
1052 doc="General failure flag for dipole fit")
1053
1054 self.edgeFlagKey = schema.addField(
1055 schema.join(name, "flag", "edge"), type="Flag",
1056 doc="Flag set when dipole is too close to edge of image")
1057
1058 def measureDipoles(self, measRecord, exposure, posExp=None, negExp=None):
1059 """Perform the non-linear least squares minimization on the putative dipole source.
1060
1061 Parameters
1062 ----------
1063 measRecord : `lsst.afw.table.SourceRecord`
1064 diaSources that will be measured using dipole measurement
1065 exposure : `lsst.afw.image.Exposure`
1066 Difference exposure on which the diaSources were detected; `exposure = posExp-negExp`
1067 If both `posExp` and `negExp` are `None`, will attempt to fit the
1068 dipole to just the `exposure` with no constraint.
1069 posExp : `lsst.afw.image.Exposure`, optional
1070 "Positive" exposure, typically a science exposure, or None if unavailable
1071 When `posExp` is `None`, will compute `posImage = exposure + negExp`.
1072 negExp : `lsst.afw.image.Exposure`, optional
1073 "Negative" exposure, typically a template exposure, or None if unavailable
1074 When `negExp` is `None`, will compute `negImage = posExp - exposure`.
1075
1076 Notes
1077 -----
1078 The main functionality of this routine was placed outside of
1079 this plugin (into `DipoleFitAlgorithm.fitDipole()`) so that
1080 `DipoleFitAlgorithm.fitDipole()` can be called separately for
1081 testing (@see `tests/testDipoleFitter.py`)
1082 """
1083 result = None
1084 pks = measRecord.getFootprint().getPeaks()
1085
1086 # Check if the footprint consists of a putative dipole - else don't fit it.
1087 if (
1088 # One peak in the footprint (not a dipole)
1089 ((nPeaks := len(pks)) <= 1)
1090 # Peaks are the same sign (not a dipole); peaks are ordered
1091 # from highest to lowest.
1092 or (nPeaks > 1 and (np.sign(pks[0].getPeakValue())
1093 == np.sign(pks[-1].getPeakValue())))
1094 ):
1095 if not self.config.fitAllDiaSources:
1096 # Non-dipoles fall back on the centroid slot for positions,
1097 # errors, and the failure flag, if we're not fitting them.
1098 measRecord[self.centroidKey.getX()] = measRecord.getX()
1099 measRecord[self.centroidKey.getY()] = measRecord.getY()
1100 self.centroidKey.getCentroidErr().set(measRecord, measRecord.getCentroidErr())
1101 measRecord[self.flagKey] = measRecord.getCentroidFlag()
1102 return
1103
1104 # Footprint is too large (not a dipole).
1105 if ((area := measRecord.getFootprint().getArea()) > self.config.maxFootprintArea):
1106 self.fail(measRecord, measBase.MeasurementError(f"{area} > {self.config.maxFootprintArea}",
1108
1109 try:
1110 alg = self.DipoleFitAlgorithmClass(exposure, posImage=posExp, negImage=negExp)
1111 result, _ = alg.fitDipole(
1112 measRecord, rel_weight=self.config.relWeight,
1113 tol=self.config.tolerance,
1114 maxSepInSigma=self.config.maxSeparation,
1115 fitBackground=self.config.fitBackground,
1116 separateNegParams=self.config.fitSeparateNegParams,
1117 verbose=False, display=False)
1118 except pexExcept.LengthError:
1119 self.fail(measRecord, measBase.MeasurementError('edge failure', self.FAILURE_EDGEFAILURE_EDGE))
1120 except Exception as e:
1121 errorMessage = f"Exception in dipole fit. {e.__class__.__name__}: {e}"
1122 self.fail(measRecord, measBase.MeasurementError(errorMessage, self.FAILURE_FITFAILURE_FIT))
1123
1124 self.log.debug("Dipole fit result: %d %s", measRecord.getId(), str(result))
1125
1126 if result is None:
1127 self.fail(measRecord, measBase.MeasurementError("bad dipole fit", self.FAILURE_FITFAILURE_FIT))
1128 return
1129
1130 # add chi2, coord/flux uncertainties (TBD), dipole classification
1131 # Add the relevant values to the measRecord
1132 measRecord[self.posFluxKey.getInstFlux()] = result.posFlux
1133 measRecord[self.posFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma!
1134 measRecord[self.posCentroidKey.getX()] = result.posCentroidX
1135 measRecord[self.posCentroidKey.getY()] = result.posCentroidY
1136
1137 measRecord[self.negFluxKey.getInstFlux()] = result.negFlux
1138 measRecord[self.negFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma!
1139 measRecord[self.negCentroidKey.getX()] = result.negCentroidX
1140 measRecord[self.negCentroidKey.getY()] = result.negCentroidY
1141
1142 # Dia source flux: average of pos+neg
1143 measRecord[self.fluxKey.getInstFlux()] = (abs(result.posFlux) + abs(result.negFlux))/2.
1144 measRecord[self.orientationKey] = result.orientation
1145 measRecord[self.separationKey] = np.sqrt((result.posCentroidX - result.negCentroidX)**2.
1146 + (result.posCentroidY - result.negCentroidY)**2.)
1147 measRecord[self.centroidKey.getX()] = result.centroidX
1148 measRecord[self.centroidKey.getY()] = result.centroidY
1149
1150 measRecord[self.signalToNoiseKey] = result.signalToNoise
1151 measRecord[self.chi2dofKey] = result.redChi2
1152
1153 if result.nData >= 1:
1154 measRecord[self.nDataKey] = result.nData
1155 else:
1156 measRecord[self.nDataKey] = 0
1157
1158 self.doClassify(measRecord, result.chi2)
1159
1160 def doClassify(self, measRecord, chi2val):
1161 """Classify a source as a dipole.
1162
1163 Parameters
1164 ----------
1165 measRecord : TODO: DM-17458
1166 TODO: DM-17458
1167 chi2val : TODO: DM-17458
1168 TODO: DM-17458
1169
1170 Notes
1171 -----
1172 Sources are classified as dipoles, or not, according to three criteria:
1173
1174 1. Does the total signal-to-noise surpass the ``minSn``?
1175 2. Are the pos/neg fluxes greater than 1.0 and no more than 0.65 (``maxFluxRatio``)
1176 of the total flux? By default this will never happen since ``posFlux == negFlux``.
1177 3. Is it a good fit (``chi2dof`` < 1)? (Currently not used.)
1178 """
1179
1180 # First, does the total signal-to-noise surpass the minSn?
1181 passesSn = measRecord[self.signalToNoiseKey] > self.config.minSn
1182
1183 # Second, are the pos/neg fluxes greater than 1.0 and no more than 0.65 (param maxFluxRatio)
1184 # of the total flux? By default this will never happen since posFlux = negFlux.
1185 passesFluxPos = (abs(measRecord[self.posFluxKey.getInstFlux()])
1186 / (measRecord[self.fluxKey.getInstFlux()]*2.)) < self.config.maxFluxRatio
1187 passesFluxPos &= (abs(measRecord[self.posFluxKey.getInstFlux()]) >= 1.0)
1188 passesFluxNeg = (abs(measRecord[self.negFluxKey.getInstFlux()])
1189 / (measRecord[self.fluxKey.getInstFlux()]*2.)) < self.config.maxFluxRatio
1190 passesFluxNeg &= (abs(measRecord[self.negFluxKey.getInstFlux()]) >= 1.0)
1191 allPass = (passesSn and passesFluxPos and passesFluxNeg) # and passesChi2)
1192
1193 # Third, is it a good fit (chi2dof < 1)?
1194 # Use scipy's chi2 cumulative distrib to estimate significance
1195 # This doesn't really work since I don't trust the values in the variance plane (which
1196 # affects the least-sq weights, which affects the resulting chi2).
1197 # But I'm going to keep this here for future use.
1198 if False:
1199 from scipy.stats import chi2
1200 ndof = chi2val / measRecord[self.chi2dofKey]
1201 significance = chi2.cdf(chi2val, ndof)
1202 passesChi2 = significance < self.config.maxChi2DoF
1203 allPass = allPass and passesChi2
1204
1205 measRecord.set(self.classificationAttemptedFlagKey, True)
1206
1207 if allPass: # Note cannot pass `allPass` into the `measRecord.set()` call below...?
1208 measRecord.set(self.classificationFlagKey, True)
1209 else:
1210 measRecord.set(self.classificationFlagKey, False)
1211
1212 def fail(self, measRecord, error=None):
1213 """Catch failures and set the correct flags.
1214
1215 Fallback on the current slot centroid positions, but set the dipole
1216 failure flag, since we attempted to fit the source.
1217 """
1218 measRecord[self.centroidKey.getX()] = measRecord.getX()
1219 measRecord[self.centroidKey.getY()] = measRecord.getY()
1220 self.centroidKey.getCentroidErr().set(measRecord, measRecord.getCentroidErr())
1221
1222 measRecord.set(self.flagKey, True)
1223 if error is not None:
1224 if error.getFlagBit() == self.FAILURE_EDGEFAILURE_EDGE:
1225 self.log.debug('DipoleFitPlugin not run on record %d: %s', measRecord.getId(), str(error))
1226 measRecord.set(self.edgeFlagKey, True)
1227 if error.getFlagBit() == self.FAILURE_FITFAILURE_FIT:
1228 self.log.warning('DipoleFitPlugin failed on record %d: %s', measRecord.getId(), str(error))
1229 if error.getFlagBit() == self.FAILURE_TOO_LARGEFAILURE_TOO_LARGE:
1230 self.log.debug('DipoleFitPlugin not run on record with too large footprint %d: %s',
1231 measRecord.getId(), str(error))
1232 else:
1233 self.log.warning('DipoleFitPlugin failed on record %d', measRecord.getId())
fitDipole(self, source, tol=1e-7, rel_weight=0.1, fitBackground=1, maxSepInSigma=5., separateNegParams=True, bgGradientOrder=1, verbose=False, display=False)
fitDipoleImpl(self, source, tol=1e-7, rel_weight=0.5, fitBackground=1, bgGradientOrder=1, maxSepInSigma=5., separateNegParams=True, verbose=False)
__init__(self, diffim, posImage=None, negImage=None)
measureDipoles(self, measRecord, exposure, posExp=None, negExp=None)
_setupSchema(self, config, name, schema, metadata)
__init__(self, config, name, schema, metadata, logName=None)
__init__(self, schema, algMetadata=None, **kwargs)
run(self, sources, exposure, posExp=None, negExp=None, **kwargs)
makeStarModel(self, bbox, psf, xcen, ycen, flux)
_getHeavyFootprintSubimage(self, fp, badfill=np.nan, grow=0)
makeModel(self, x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None, b=None, x1=None, y1=None, xy=None, x2=None, y2=None, bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None, **kwargs)
fitFootprintBackground(self, source, posImage, order=1)
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
daf::base::PropertySet * set
Definition fits.cc:931