LSSTApplications  18.1.0
LSSTDataManagementBasePackage
dipoleFitTask.py
Go to the documentation of this file.
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 #
22 
23 import numpy as np
24 import warnings
25 
26 import lsst.afw.geom as afwGeom
27 import lsst.afw.image as afwImage
28 import lsst.meas.base as measBase
29 import lsst.afw.table as afwTable
30 import lsst.afw.detection as afwDet
31 from lsst.log import Log
32 import lsst.pex.exceptions as pexExcept
33 import lsst.pex.config as pexConfig
34 from lsst.pipe.base import Struct, 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 
44 class DipoleFitPluginConfig(measBase.SingleFramePluginConfig):
45  """Configuration for DipoleFitPlugin
46  """
47 
48  fitAllDiaSources = pexConfig.Field(
49  dtype=float, 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 
95 class DipoleFitTaskConfig(measBase.SingleFrameMeasurementConfig):
96  """Measurement of detected diaSources as dipoles
97 
98  Currently we keep the "old" DipoleMeasurement algorithms turned on.
99  """
100 
101  def setDefaults(self):
102  measBase.SingleFrameMeasurementConfig.setDefaults(self)
103 
104  self.plugins.names = ["base_CircularApertureFlux",
105  "base_PixelFlags",
106  "base_SkyCoord",
107  "base_PsfFlux",
108  "base_SdssCentroid",
109  "base_SdssShape",
110  "base_SdssCentroid",
111  "base_GaussianFlux",
112  "base_PeakLikelihoodFlux",
113  "base_PeakCentroid",
114  "base_NaiveCentroid",
115  "ip_diffim_NaiveDipoleCentroid",
116  "ip_diffim_NaiveDipoleFlux",
117  "ip_diffim_PsfDipoleFlux",
118  "ip_diffim_ClassificationDipole",
119  ]
120 
121  self.slots.calibFlux = None
122  self.slots.modelFlux = None
123  self.slots.gaussianFlux = None
124  self.slots.shape = "base_SdssShape"
125  self.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
126  self.doReplaceWithNoise = False
127 
128 
129 class DipoleFitTask(measBase.SingleFrameMeasurementTask):
130  """A task that fits a dipole to a difference image, with an optional separate detection image.
131 
132  Because it subclasses SingleFrameMeasurementTask, and calls
133  SingleFrameMeasurementTask.run() from its run() method, it still
134  can be used identically to a standard SingleFrameMeasurementTask.
135  """
136 
137  ConfigClass = DipoleFitTaskConfig
138  _DefaultName = "ip_diffim_DipoleFit"
139 
140  def __init__(self, schema, algMetadata=None, **kwargs):
141 
142  measBase.SingleFrameMeasurementTask.__init__(self, schema, algMetadata, **kwargs)
143 
144  dpFitPluginConfig = self.config.plugins['ip_diffim_DipoleFit']
145 
146  self.dipoleFitter = DipoleFitPlugin(dpFitPluginConfig, name=self._DefaultName,
147  schema=schema, metadata=algMetadata)
148 
149  @timeMethod
150  def run(self, sources, exposure, posExp=None, negExp=None, **kwargs):
151  """Run dipole measurement and classification
152 
153  Parameters
154  ----------
155  sources : `lsst.afw.table.SourceCatalog`
156  ``diaSources`` that will be measured using dipole measurement
157  exposure : `lsst.afw.image.Exposure`
158  The difference exposure on which the ``diaSources`` of the ``sources`` parameter
159  were detected. If neither ``posExp`` nor ``negExp`` are set, then the dipole is also
160  fitted directly to this difference image.
161  posExp : `lsst.afw.image.Exposure`, optional
162  "Positive" exposure, typically a science exposure, or None if unavailable
163  When `posExp` is `None`, will compute `posImage = exposure + negExp`.
164  negExp : `lsst.afw.image.Exposure`, optional
165  "Negative" exposure, typically a template exposure, or None if unavailable
166  When `negExp` is `None`, will compute `negImage = posExp - exposure`.
167  **kwargs
168  Additional keyword arguments for `lsst.meas.base.sfm.SingleFrameMeasurementTask`.
169  """
170 
171  measBase.SingleFrameMeasurementTask.run(self, sources, exposure, **kwargs)
172 
173  if not sources:
174  return
175 
176  for source in sources:
177  self.dipoleFitter.measure(source, exposure, posExp, negExp)
178 
179 
181  """Lightweight class containing methods for generating a dipole model for fitting
182  to sources in diffims, used by DipoleFitAlgorithm.
183 
184  See also:
185  `DMTN-007: Dipole characterization for image differencing <https://dmtn-007.lsst.io>`_.
186  """
187 
188  def __init__(self):
189  import lsstDebug
190  self.debug = lsstDebug.Info(__name__).debug
191  self.log = Log.getLogger(__name__)
192 
193  def makeBackgroundModel(self, in_x, pars=None):
194  """Generate gradient model (2-d array) with up to 2nd-order polynomial
195 
196  Parameters
197  ----------
198  in_x : `numpy.array`
199  (2, w, h)-dimensional `numpy.array`, containing the
200  input x,y meshgrid providing the coordinates upon which to
201  compute the gradient. This will typically be generated via
202  `_generateXYGrid()`. `w` and `h` correspond to the width and
203  height of the desired grid.
204  pars : `list` of `float`, optional
205  Up to 6 floats for up
206  to 6 2nd-order 2-d polynomial gradient parameters, in the
207  following order: (intercept, x, y, xy, x**2, y**2). If `pars`
208  is emtpy or `None`, do nothing and return `None` (for speed).
209 
210  Returns
211  -------
212  result : `None` or `numpy.array`
213  return None, or 2-d numpy.array of width/height matching
214  input bbox, containing computed gradient values.
215  """
216 
217  # Don't fit for other gradient parameters if the intercept is not included.
218  if (pars is None) or (len(pars) <= 0) or (pars[0] is None):
219  return
220 
221  y, x = in_x[0, :], in_x[1, :]
222  gradient = np.full_like(x, pars[0], dtype='float64')
223  if len(pars) > 1 and pars[1] is not None:
224  gradient += pars[1] * x
225  if len(pars) > 2 and pars[2] is not None:
226  gradient += pars[2] * y
227  if len(pars) > 3 and pars[3] is not None:
228  gradient += pars[3] * (x * y)
229  if len(pars) > 4 and pars[4] is not None:
230  gradient += pars[4] * (x * x)
231  if len(pars) > 5 and pars[5] is not None:
232  gradient += pars[5] * (y * y)
233 
234  return gradient
235 
236  def _generateXYGrid(self, bbox):
237  """Generate a meshgrid covering the x,y coordinates bounded by bbox
238 
239  Parameters
240  ----------
241  bbox : `lsst.geom.Box`
242  input BoundingBox defining the coordinate limits
243 
244  Returns
245  -------
246  in_x : `numpy.array`
247  (2, w, h)-dimensional numpy array containing the grid indexing over x- and
248  y- coordinates
249  """
250 
251  x, y = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
252  in_x = np.array([y, x]).astype(np.float64)
253  in_x[0, :] -= np.mean(in_x[0, :])
254  in_x[1, :] -= np.mean(in_x[1, :])
255  return in_x
256 
257  def _getHeavyFootprintSubimage(self, fp, badfill=np.nan, grow=0):
258  """Extract the image from a ``~lsst.afw.detection.HeavyFootprint``
259  as an `lsst.afw.image.ImageF`.
260 
261  Parameters
262  ----------
263  fp : `lsst.afw.detection.HeavyFootprint`
264  HeavyFootprint to use to generate the subimage
265  badfill : `float`, optional
266  Value to fill in pixels in extracted image that are outside the footprint
267  grow : `int`
268  Optionally grow the footprint by this amount before extraction
269 
270  Returns
271  -------
272  subim2 : `lsst.afw.image.ImageF`
273  An `~lsst.afw.image.ImageF` containing the subimage.
274  """
275  bbox = fp.getBBox()
276  if grow > 0:
277  bbox.grow(grow)
278 
279  subim2 = afwImage.ImageF(bbox, badfill)
280  fp.getSpans().unflatten(subim2.getArray(), fp.getImageArray(), bbox.getCorners()[0])
281  return subim2
282 
283  def fitFootprintBackground(self, source, posImage, order=1):
284  """Fit a linear (polynomial) model of given order (max 2) to the background of a footprint.
285 
286  Only fit the pixels OUTSIDE of the footprint, but within its bounding box.
287 
288  Parameters
289  ----------
290  source : `lsst.afw.table.SourceRecord`
291  SourceRecord, the footprint of which is to be fit
292  posImage : `lsst.afw.image.Exposure`
293  The exposure from which to extract the footprint subimage
294  order : `int`
295  Polynomial order of background gradient to fit.
296 
297  Returns
298  -------
299  pars : `tuple` of `float`
300  `tuple` of length (1 if order==0; 3 if order==1; 6 if order == 2),
301  containing the resulting fit parameters
302  """
303 
304  # TODO look into whether to use afwMath background methods -- see
305  # http://lsst-web.ncsa.illinois.edu/doxygen/x_masterDoxyDoc/_background_example.html
306  fp = source.getFootprint()
307  bbox = fp.getBBox()
308  bbox.grow(3)
309  posImg = afwImage.ImageF(posImage.getMaskedImage().getImage(), bbox, afwImage.PARENT)
310 
311  # This code constructs the footprint image so that we can identify the pixels that are
312  # outside the footprint (but within the bounding box). These are the pixels used for
313  # fitting the background.
314  posHfp = afwDet.HeavyFootprintF(fp, posImage.getMaskedImage())
315  posFpImg = self._getHeavyFootprintSubimage(posHfp, grow=3)
316 
317  isBg = np.isnan(posFpImg.getArray()).ravel()
318 
319  data = posImg.getArray().ravel()
320  data = data[isBg]
321  B = data
322 
323  x, y = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
324  x = x.astype(np.float64).ravel()
325  x -= np.mean(x)
326  x = x[isBg]
327  y = y.astype(np.float64).ravel()
328  y -= np.mean(y)
329  y = y[isBg]
330  b = np.ones_like(x, dtype=np.float64)
331 
332  M = np.vstack([b]).T # order = 0
333  if order == 1:
334  M = np.vstack([b, x, y]).T
335  elif order == 2:
336  M = np.vstack([b, x, y, x**2., y**2., x*y]).T
337 
338  pars = np.linalg.lstsq(M, B, rcond=-1)[0]
339  return pars
340 
341  def makeStarModel(self, bbox, psf, xcen, ycen, flux):
342  """Generate a 2D image model of a single PDF centered at the given coordinates.
343 
344  Parameters
345  ----------
346  bbox : `lsst.geom.Box`
347  Bounding box marking pixel coordinates for generated model
348  psf : TODO: DM-17458
349  Psf model used to generate the 'star'
350  xcen : `float`
351  Desired x-centroid of the 'star'
352  ycen : `float`
353  Desired y-centroid of the 'star'
354  flux : `float`
355  Desired flux of the 'star'
356 
357  Returns
358  -------
359  p_Im : `lsst.afw.image.Image`
360  2-d stellar image of width/height matching input ``bbox``,
361  containing PSF with given centroid and flux
362  """
363 
364  # Generate the psf image, normalize to flux
365  psf_img = psf.computeImage(afwGeom.Point2D(xcen, ycen)).convertF()
366  psf_img_sum = np.nansum(psf_img.getArray())
367  psf_img *= (flux/psf_img_sum)
368 
369  # Clip the PSF image bounding box to fall within the footprint bounding box
370  psf_box = psf_img.getBBox()
371  psf_box.clip(bbox)
372  psf_img = afwImage.ImageF(psf_img, psf_box, afwImage.PARENT)
373 
374  # Then actually crop the psf image.
375  # Usually not necessary, but if the dipole is near the edge of the image...
376  # Would be nice if we could compare original pos_box with clipped pos_box and
377  # see if it actually was clipped.
378  p_Im = afwImage.ImageF(bbox)
379  tmpSubim = afwImage.ImageF(p_Im, psf_box, afwImage.PARENT)
380  tmpSubim += psf_img
381 
382  return p_Im
383 
384  def makeModel(self, x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None,
385  b=None, x1=None, y1=None, xy=None, x2=None, y2=None,
386  bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None,
387  **kwargs):
388  """Generate dipole model with given parameters.
389 
390  This is the function whose sum-of-squared difference from data
391  is minimized by `lmfit`.
392 
393  x : TODO: DM-17458
394  Input independent variable. Used here as the grid on
395  which to compute the background gradient model.
396  flux : `float`
397  Desired flux of the positive lobe of the dipole
398  xcenPos : `float`
399  Desired x-centroid of the positive lobe of the dipole
400  ycenPos : `float`
401  Desired y-centroid of the positive lobe of the dipole
402  xcenNeg : `float`
403  Desired x-centroid of the negative lobe of the dipole
404  ycenNeg : `float`
405  Desired y-centroid of the negative lobe of the dipole
406  fluxNeg : `float`, optional
407  Desired flux of the negative lobe of the dipole, set to 'flux' if None
408  b, x1, y1, xy, x2, y2 : `float`
409  Gradient parameters for positive lobe.
410  bNeg, x1Neg, y1Neg, xyNeg, x2Neg, y2Neg : `float`, optional
411  Gradient parameters for negative lobe.
412  They are set to the corresponding positive values if None.
413 
414  **kwargs
415  Keyword arguments passed through ``lmfit`` and
416  used by this function. These must include:
417 
418  - ``psf`` Psf model used to generate the 'star'
419  - ``rel_weight`` Used to signify least-squares weighting of posImage/negImage
420  relative to diffim. If ``rel_weight == 0`` then posImage/negImage are ignored.
421  - ``bbox`` Bounding box containing region to be modelled
422 
423  Returns
424  -------
425  zout : `numpy.array`
426  Has width and height matching the input bbox, and
427  contains the dipole model with given centroids and flux(es). If
428  ``rel_weight`` = 0, this is a 2-d array with dimensions matching
429  those of bbox; otherwise a stack of three such arrays,
430  representing the dipole (diffim), positive and negative images
431  respectively.
432  """
433 
434  psf = kwargs.get('psf')
435  rel_weight = kwargs.get('rel_weight') # if > 0, we're including pre-sub. images
436  fp = kwargs.get('footprint')
437  bbox = fp.getBBox()
438 
439  if fluxNeg is None:
440  fluxNeg = flux
441 
442  if self.debug:
443  self.log.debug('%.2f %.2f %.2f %.2f %.2f %.2f',
444  flux, fluxNeg, xcenPos, ycenPos, xcenNeg, ycenNeg)
445  if x1 is not None:
446  self.log.debug(' %.2f %.2f %.2f', b, x1, y1)
447  if xy is not None:
448  self.log.debug(' %.2f %.2f %.2f', xy, x2, y2)
449 
450  posIm = self.makeStarModel(bbox, psf, xcenPos, ycenPos, flux)
451  negIm = self.makeStarModel(bbox, psf, xcenNeg, ycenNeg, fluxNeg)
452 
453  in_x = x
454  if in_x is None: # use the footprint to generate the input grid
455  y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
456  in_x = np.array([x, y]) * 1.
457  in_x[0, :] -= in_x[0, :].mean() # center it!
458  in_x[1, :] -= in_x[1, :].mean()
459 
460  if b is not None:
461  gradient = self.makeBackgroundModel(in_x, (b, x1, y1, xy, x2, y2))
462 
463  # If bNeg is None, then don't fit the negative background separately
464  if bNeg is not None:
465  gradientNeg = self.makeBackgroundModel(in_x, (bNeg, x1Neg, y1Neg, xyNeg, x2Neg, y2Neg))
466  else:
467  gradientNeg = gradient
468 
469  posIm.getArray()[:, :] += gradient
470  negIm.getArray()[:, :] += gradientNeg
471 
472  # Generate the diffIm model
473  diffIm = afwImage.ImageF(bbox)
474  diffIm += posIm
475  diffIm -= negIm
476 
477  zout = diffIm.getArray()
478  if rel_weight > 0.:
479  zout = np.append([zout], [posIm.getArray(), negIm.getArray()], axis=0)
480 
481  return zout
482 
483 
485  """Fit a dipole model using an image difference.
486 
487  See also:
488  `DMTN-007: Dipole characterization for image differencing <https://dmtn-007.lsst.io>`_.
489  """
490 
491  # This is just a private version number to sync with the ipython notebooks that I have been
492  # using for algorithm development.
493  _private_version_ = '0.0.5'
494 
495  # Below is a (somewhat incomplete) list of improvements
496  # that would be worth investigating, given the time:
497 
498  # todo 1. evaluate necessity for separate parameters for pos- and neg- images
499  # todo 2. only fit background OUTSIDE footprint (DONE) and dipole params INSIDE footprint (NOT DONE)?
500  # todo 3. correct normalization of least-squares weights based on variance planes
501  # todo 4. account for PSFs that vary across the exposures (should be happening by default?)
502  # todo 5. correctly account for NA/masks (i.e., ignore!)
503  # todo 6. better exception handling in the plugin
504  # todo 7. better classification of dipoles (e.g. by comparing chi2 fit vs. monopole?)
505  # todo 8. (DONE) Initial fast estimate of background gradient(s) params -- perhaps using numpy.lstsq
506  # todo 9. (NOT NEEDED - see (2)) Initial fast test whether a background gradient needs to be fit
507  # todo 10. (DONE) better initial estimate for flux when there's a strong gradient
508  # todo 11. (DONE) requires a new package `lmfit` -- investiate others? (astropy/scipy/iminuit?)
509 
510  def __init__(self, diffim, posImage=None, negImage=None):
511  """Algorithm to run dipole measurement on a diaSource
512 
513  Parameters
514  ----------
515  diffim : `lsst.afw.image.Exposure`
516  Exposure on which the diaSources were detected
517  posImage : `lsst.afw.image.Exposure`
518  "Positive" exposure from which the template was subtracted
519  negImage : `lsst.afw.image.Exposure`
520  "Negative" exposure which was subtracted from the posImage
521  """
522 
523  self.diffim = diffim
524  self.posImage = posImage
525  self.negImage = negImage
526  self.psfSigma = None
527  if diffim is not None:
528  self.psfSigma = diffim.getPsf().computeShape().getDeterminantRadius()
529 
530  self.log = Log.getLogger(__name__)
531 
532  import lsstDebug
533  self.debug = lsstDebug.Info(__name__).debug
534 
535  def fitDipoleImpl(self, source, tol=1e-7, rel_weight=0.5,
536  fitBackground=1, bgGradientOrder=1, maxSepInSigma=5.,
537  separateNegParams=True, verbose=False):
538  """Fit a dipole model to an input difference image.
539 
540  Actually, fits the subimage bounded by the input source's
541  footprint) and optionally constrain the fit using the
542  pre-subtraction images posImage and negImage.
543 
544  Parameters
545  ----------
546  source : TODO: DM-17458
547  TODO: DM-17458
548  tol : float, optional
549  TODO: DM-17458
550  rel_weight : `float`, optional
551  TODO: DM-17458
552  fitBackground : `int`, optional
553  TODO: DM-17458
554  bgGradientOrder : `int`, optional
555  TODO: DM-17458
556  maxSepInSigma : `float`, optional
557  TODO: DM-17458
558  separateNegParams : `bool`, optional
559  TODO: DM-17458
560  verbose : `bool`, optional
561  TODO: DM-17458
562 
563  Returns
564  -------
565  result : `lmfit.MinimizerResult`
566  return `lmfit.MinimizerResult` object containing the fit
567  parameters and other information.
568  """
569 
570  # Only import lmfit if someone wants to use the new DipoleFitAlgorithm.
571  import lmfit
572 
573  fp = source.getFootprint()
574  bbox = fp.getBBox()
575  subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)
576 
577  z = diArr = subim.getArrays()[0]
578  weights = 1. / subim.getArrays()[2] # get the weights (=1/variance)
579 
580  if rel_weight > 0. and ((self.posImage is not None) or (self.negImage is not None)):
581  if self.negImage is not None:
582  negSubim = afwImage.MaskedImageF(self.negImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
583  if self.posImage is not None:
584  posSubim = afwImage.MaskedImageF(self.posImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
585  if self.posImage is None: # no science image provided; generate it from diffim + negImage
586  posSubim = subim.clone()
587  posSubim += negSubim
588  if self.negImage is None: # no template provided; generate it from the posImage - diffim
589  negSubim = posSubim.clone()
590  negSubim -= subim
591 
592  z = np.append([z], [posSubim.getArrays()[0],
593  negSubim.getArrays()[0]], axis=0)
594  # Weight the pos/neg images by rel_weight relative to the diffim
595  weights = np.append([weights], [1. / posSubim.getArrays()[2] * rel_weight,
596  1. / negSubim.getArrays()[2] * rel_weight], axis=0)
597  else:
598  rel_weight = 0. # a short-cut for "don't include the pre-subtraction data"
599 
600  # It seems that `lmfit` requires a static functor as its optimized method, which eliminates
601  # the ability to pass a bound method or other class method. Here we write a wrapper which
602  # makes this possible.
603  def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None,
604  b=None, x1=None, y1=None, xy=None, x2=None, y2=None,
605  bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None,
606  **kwargs):
607  """Generate dipole model with given parameters.
608 
609  It simply defers to `modelObj.makeModel()`, where `modelObj` comes
610  out of `kwargs['modelObj']`.
611  """
612  modelObj = kwargs.pop('modelObj')
613  return modelObj.makeModel(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=fluxNeg,
614  b=b, x1=x1, y1=y1, xy=xy, x2=x2, y2=y2,
615  bNeg=bNeg, x1Neg=x1Neg, y1Neg=y1Neg, xyNeg=xyNeg,
616  x2Neg=x2Neg, y2Neg=y2Neg, **kwargs)
617 
618  dipoleModel = DipoleModel()
619 
620  modelFunctor = dipoleModelFunctor # dipoleModel.makeModel does not work for now.
621  # Create the lmfit model (lmfit uses scipy 'leastsq' option by default - Levenberg-Marquardt)
622  # Note we can also tell it to drop missing values from the data.
623  gmod = lmfit.Model(modelFunctor, verbose=verbose, missing='drop')
624  # independent_vars=independent_vars) #, param_names=param_names)
625 
626  # Add the constraints for centroids, fluxes.
627  # starting constraint - near centroid of footprint
628  fpCentroid = np.array([fp.getCentroid().getX(), fp.getCentroid().getY()])
629  cenNeg = cenPos = fpCentroid
630 
631  pks = fp.getPeaks()
632 
633  if len(pks) >= 1:
634  cenPos = pks[0].getF() # if individual (merged) peaks were detected, use those
635  if len(pks) >= 2: # peaks are already sorted by centroid flux so take the most negative one
636  cenNeg = pks[-1].getF()
637 
638  # For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
639  # First assume dipole is not separated by more than 5*psfSigma.
640  maxSep = self.psfSigma * maxSepInSigma
641 
642  # As an initial guess -- assume the dipole is close to the center of the footprint.
643  if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
644  cenPos = fpCentroid
645  if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
646  cenPos = fpCentroid
647 
648  # parameter hints/constraints: https://lmfit.github.io/lmfit-py/model.html#model-param-hints-section
649  # might make sense to not use bounds -- see http://lmfit.github.io/lmfit-py/bounds.html
650  # also see this discussion -- https://github.com/scipy/scipy/issues/3129
651  gmod.set_param_hint('xcenPos', value=cenPos[0],
652  min=cenPos[0]-maxSep, max=cenPos[0]+maxSep)
653  gmod.set_param_hint('ycenPos', value=cenPos[1],
654  min=cenPos[1]-maxSep, max=cenPos[1]+maxSep)
655  gmod.set_param_hint('xcenNeg', value=cenNeg[0],
656  min=cenNeg[0]-maxSep, max=cenNeg[0]+maxSep)
657  gmod.set_param_hint('ycenNeg', value=cenNeg[1],
658  min=cenNeg[1]-maxSep, max=cenNeg[1]+maxSep)
659 
660  # Use the (flux under the dipole)*5 for an estimate.
661  # Lots of testing showed that having startingFlux be too high was better than too low.
662  startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
663  posFlux = negFlux = startingFlux
664 
665  # TBD: set max. flux limit?
666  gmod.set_param_hint('flux', value=posFlux, min=0.1)
667 
668  if separateNegParams:
669  # TBD: set max negative lobe flux limit?
670  gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
671 
672  # Fixed parameters (don't fit for them if there are no pre-sub images or no gradient fit requested):
673  # Right now (fitBackground == 1), we fit a linear model to the background and then subtract
674  # it from the data and then don't fit the background again (this is faster).
675  # A slower alternative (fitBackground == 2) is to use the estimated background parameters as
676  # starting points in the integrated model fit. That is currently not performed by default,
677  # but might be desirable in some cases.
678  bgParsPos = bgParsNeg = (0., 0., 0.)
679  if ((rel_weight > 0.) and (fitBackground != 0) and (bgGradientOrder >= 0)):
680  pbg = 0.
681  bgFitImage = self.posImage if self.posImage is not None else self.negImage
682  # Fit the gradient to the background (linear model)
683  bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
684  order=bgGradientOrder)
685 
686  # Generate the gradient and subtract it from the pre-subtraction image data
687  if fitBackground == 1:
688  in_x = dipoleModel._generateXYGrid(bbox)
689  pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsPos))
690  z[1, :] -= pbg
691  z[1, :] -= np.nanmedian(z[1, :])
692  posFlux = np.nansum(z[1, :])
693  gmod.set_param_hint('flux', value=posFlux*1.5, min=0.1)
694 
695  if separateNegParams and self.negImage is not None:
696  bgParsNeg = dipoleModel.fitFootprintBackground(source, self.negImage,
697  order=bgGradientOrder)
698  pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsNeg))
699  z[2, :] -= pbg
700  z[2, :] -= np.nanmedian(z[2, :])
701  if separateNegParams:
702  negFlux = np.nansum(z[2, :])
703  gmod.set_param_hint('fluxNeg', value=negFlux*1.5, min=0.1)
704 
705  # Do not subtract the background from the images but include the background parameters in the fit
706  if fitBackground == 2:
707  if bgGradientOrder >= 0:
708  gmod.set_param_hint('b', value=bgParsPos[0])
709  if separateNegParams:
710  gmod.set_param_hint('bNeg', value=bgParsNeg[0])
711  if bgGradientOrder >= 1:
712  gmod.set_param_hint('x1', value=bgParsPos[1])
713  gmod.set_param_hint('y1', value=bgParsPos[2])
714  if separateNegParams:
715  gmod.set_param_hint('x1Neg', value=bgParsNeg[1])
716  gmod.set_param_hint('y1Neg', value=bgParsNeg[2])
717  if bgGradientOrder >= 2:
718  gmod.set_param_hint('xy', value=bgParsPos[3])
719  gmod.set_param_hint('x2', value=bgParsPos[4])
720  gmod.set_param_hint('y2', value=bgParsPos[5])
721  if separateNegParams:
722  gmod.set_param_hint('xyNeg', value=bgParsNeg[3])
723  gmod.set_param_hint('x2Neg', value=bgParsNeg[4])
724  gmod.set_param_hint('y2Neg', value=bgParsNeg[5])
725 
726  y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
727  in_x = np.array([x, y]).astype(np.float)
728  in_x[0, :] -= in_x[0, :].mean() # center it!
729  in_x[1, :] -= in_x[1, :].mean()
730 
731  # Instead of explicitly using a mask to ignore flagged pixels, just set the ignored pixels'
732  # weights to 0 in the fit. TBD: need to inspect mask planes to set this mask.
733  mask = np.ones_like(z, dtype=bool) # TBD: set mask values to False if the pixels are to be ignored
734 
735  # I'm not sure about the variance planes in the diffim (or convolved pre-sub. images
736  # for that matter) so for now, let's just do an un-weighted least-squares fit
737  # (override weights computed above).
738  weights = mask.astype(np.float64)
739  if self.posImage is not None and rel_weight > 0.:
740  weights = np.array([np.ones_like(diArr), np.ones_like(diArr)*rel_weight,
741  np.ones_like(diArr)*rel_weight])
742 
743  # Set the weights to zero if mask is False
744  if np.any(~mask):
745  weights[~mask] = 0.
746 
747  # Note that although we can, we're not required to set initial values for params here,
748  # since we set their param_hint's above.
749  # Can add "method" param to not use 'leastsq' (==levenberg-marquardt), e.g. "method='nelder'"
750  with warnings.catch_warnings():
751  warnings.simplefilter("ignore") # temporarily turn off silly lmfit warnings
752  result = gmod.fit(z, weights=weights, x=in_x,
753  verbose=verbose,
754  fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
755  'maxfev': 250}, # see scipy docs
756  psf=self.diffim.getPsf(), # hereon: kwargs that get passed to genDipoleModel()
757  rel_weight=rel_weight,
758  footprint=fp,
759  modelObj=dipoleModel)
760 
761  if verbose: # the ci_report() seems to fail if neg params are constrained -- TBD why.
762  # Never wanted in production - this takes a long time (longer than the fit!)
763  # This is how to get confidence intervals out:
764  # https://lmfit.github.io/lmfit-py/confidence.html and
765  # http://cars9.uchicago.edu/software/python/lmfit/model.html
766  print(result.fit_report(show_correl=False))
767  if separateNegParams:
768  print(result.ci_report())
769 
770  return result
771 
772  def fitDipole(self, source, tol=1e-7, rel_weight=0.1,
773  fitBackground=1, maxSepInSigma=5., separateNegParams=True,
774  bgGradientOrder=1, verbose=False, display=False):
775  """Fit a dipole model to an input ``diaSource`` (wraps `fitDipoleImpl`).
776 
777  Actually, fits the subimage bounded by the input source's
778  footprint) and optionally constrain the fit using the
779  pre-subtraction images self.posImage (science) and
780  self.negImage (template). Wraps the output into a
781  `pipeBase.Struct` named tuple after computing additional
782  statistics such as orientation and SNR.
783 
784  Parameters
785  ----------
786  source : `lsst.afw.table.SourceRecord`
787  Record containing the (merged) dipole source footprint detected on the diffim
788  tol : `float`, optional
789  Tolerance parameter for scipy.leastsq() optimization
790  rel_weight : `float`, optional
791  Weighting of posImage/negImage relative to the diffim in the fit
792  fitBackground : `int`, {0, 1, 2}, optional
793  How to fit linear background gradient in posImage/negImage
794 
795  - 0: do not fit background at all
796  - 1 (default): pre-fit the background using linear least squares and then do not fit it
797  as part of the dipole fitting optimization
798  - 2: pre-fit the background using linear least squares (as in 1), and use the parameter
799  estimates from that fit as starting parameters for an integrated "re-fit" of the
800  background as part of the overall dipole fitting optimization.
801  maxSepInSigma : `float`, optional
802  Allowed window of centroid parameters relative to peak in input source footprint
803  separateNegParams : `bool`, optional
804  Fit separate parameters to the flux and background gradient in
805  bgGradientOrder : `int`, {0, 1, 2}, optional
806  Desired polynomial order of background gradient
807  verbose: `bool`, optional
808  Be verbose
809  display
810  Display input data, best fit model(s) and residuals in a matplotlib window.
811 
812  Returns
813  -------
814  result : `struct`
815  `pipeBase.Struct` object containing the fit parameters and other information.
816 
817  result : `callable`
818  `lmfit.MinimizerResult` object for debugging and error estimation, etc.
819 
820  Notes
821  -----
822  Parameter `fitBackground` has three options, thus it is an integer:
823 
824  """
825 
826  fitResult = self.fitDipoleImpl(
827  source, tol=tol, rel_weight=rel_weight, fitBackground=fitBackground,
828  maxSepInSigma=maxSepInSigma, separateNegParams=separateNegParams,
829  bgGradientOrder=bgGradientOrder, verbose=verbose)
830 
831  # Display images, model fits and residuals (currently uses matplotlib display functions)
832  if display:
833  fp = source.getFootprint()
834  self.displayFitResults(fp, fitResult)
835 
836  fitParams = fitResult.best_values
837  if fitParams['flux'] <= 1.: # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit.
838  out = Struct(posCentroidX=np.nan, posCentroidY=np.nan,
839  negCentroidX=np.nan, negCentroidY=np.nan,
840  posFlux=np.nan, negFlux=np.nan, posFluxErr=np.nan, negFluxErr=np.nan,
841  centroidX=np.nan, centroidY=np.nan, orientation=np.nan,
842  signalToNoise=np.nan, chi2=np.nan, redChi2=np.nan)
843  return out, fitResult
844 
845  centroid = ((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2.,
846  (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2.)
847  dx, dy = fitParams['xcenPos'] - fitParams['xcenNeg'], fitParams['ycenPos'] - fitParams['ycenNeg']
848  angle = np.arctan2(dy, dx) / np.pi * 180. # convert to degrees (should keep as rad?)
849 
850  # Exctract flux value, compute signalToNoise from flux/variance_within_footprint
851  # Also extract the stderr of flux estimate.
852  def computeSumVariance(exposure, footprint):
853  box = footprint.getBBox()
854  subim = afwImage.MaskedImageF(exposure.getMaskedImage(), box, origin=afwImage.PARENT)
855  return np.sqrt(np.nansum(subim.getArrays()[1][:, :]))
856 
857  fluxVal = fluxVar = fitParams['flux']
858  fluxErr = fluxErrNeg = fitResult.params['flux'].stderr
859  if self.posImage is not None:
860  fluxVar = computeSumVariance(self.posImage, source.getFootprint())
861  else:
862  fluxVar = computeSumVariance(self.diffim, source.getFootprint())
863 
864  fluxValNeg, fluxVarNeg = fluxVal, fluxVar
865  if separateNegParams:
866  fluxValNeg = fitParams['fluxNeg']
867  fluxErrNeg = fitResult.params['fluxNeg'].stderr
868  if self.negImage is not None:
869  fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint())
870 
871  try:
872  signalToNoise = np.sqrt((fluxVal/fluxVar)**2 + (fluxValNeg/fluxVarNeg)**2)
873  except ZeroDivisionError: # catch divide by zero - should never happen.
874  signalToNoise = np.nan
875 
876  out = Struct(posCentroidX=fitParams['xcenPos'], posCentroidY=fitParams['ycenPos'],
877  negCentroidX=fitParams['xcenNeg'], negCentroidY=fitParams['ycenNeg'],
878  posFlux=fluxVal, negFlux=-fluxValNeg, posFluxErr=fluxErr, negFluxErr=fluxErrNeg,
879  centroidX=centroid[0], centroidY=centroid[1], orientation=angle,
880  signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi)
881 
882  # fitResult may be returned for debugging
883  return out, fitResult
884 
885  def displayFitResults(self, footprint, result):
886  """Display data, model fits and residuals (currently uses matplotlib display functions).
887 
888  Parameters
889  ----------
890  footprint : TODO: DM-17458
891  Footprint containing the dipole that was fit
892  result : `lmfit.MinimizerResult`
893  `lmfit.MinimizerResult` object returned by `lmfit` optimizer
894 
895  Returns
896  -------
897  fig : `matplotlib.pyplot.plot`
898  """
899  try:
900  import matplotlib.pyplot as plt
901  except ImportError as err:
902  self.log.warn('Unable to import matplotlib: %s', err)
903  raise err
904 
905  def display2dArray(arr, title='Data', extent=None):
906  """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range.
907  """
908  fig = plt.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent)
909  plt.title(title)
910  plt.colorbar(fig, cmap='gray')
911  return fig
912 
913  z = result.data
914  fit = result.best_fit
915  bbox = footprint.getBBox()
916  extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY())
917  if z.shape[0] == 3:
918  fig = plt.figure(figsize=(8, 8))
919  for i in range(3):
920  plt.subplot(3, 3, i*3+1)
921  display2dArray(z[i, :], 'Data', extent=extent)
922  plt.subplot(3, 3, i*3+2)
923  display2dArray(fit[i, :], 'Model', extent=extent)
924  plt.subplot(3, 3, i*3+3)
925  display2dArray(z[i, :] - fit[i, :], 'Residual', extent=extent)
926  return fig
927  else:
928  fig = plt.figure(figsize=(8, 2.5))
929  plt.subplot(1, 3, 1)
930  display2dArray(z, 'Data', extent=extent)
931  plt.subplot(1, 3, 2)
932  display2dArray(fit, 'Model', extent=extent)
933  plt.subplot(1, 3, 3)
934  display2dArray(z - fit, 'Residual', extent=extent)
935  return fig
936 
937  plt.show()
938 
939 
940 @measBase.register("ip_diffim_DipoleFit")
941 class DipoleFitPlugin(measBase.SingleFramePlugin):
942  """A single frame measurement plugin that fits dipoles to all merged (two-peak) ``diaSources``.
943 
944  This measurement plugin accepts up to three input images in
945  its `measure` method. If these are provided, it includes data
946  from the pre-subtraction posImage (science image) and optionally
947  negImage (template image) to constrain the fit. The meat of the
948  fitting routines are in the class `~lsst.module.name.DipoleFitAlgorithm`.
949 
950  Notes
951  -----
952  The motivation behind this plugin and the necessity for including more than
953  one exposure are documented in DMTN-007 (http://dmtn-007.lsst.io).
954 
955  This class is named `ip_diffim_DipoleFit` so that it may be used alongside
956  the existing `ip_diffim_DipoleMeasurement` classes until such a time as those
957  are deemed to be replaceable by this.
958  """
959 
960  ConfigClass = DipoleFitPluginConfig
961  DipoleFitAlgorithmClass = DipoleFitAlgorithm # Pointer to the class that performs the fit
962 
963  FAILURE_EDGE = 1 # too close to the edge
964  FAILURE_FIT = 2 # failure in the fitting
965  FAILURE_NOT_DIPOLE = 4 # input source is not a putative dipole to begin with
966 
967  @classmethod
969  """Set execution order to `FLUX_ORDER`.
970 
971  This includes algorithms that require both `getShape()` and `getCentroid()`,
972  in addition to a Footprint and its Peaks.
973  """
974  return cls.FLUX_ORDER
975 
976  def __init__(self, config, name, schema, metadata):
977  measBase.SingleFramePlugin.__init__(self, config, name, schema, metadata)
978 
979  self.log = Log.getLogger(name)
980 
981  self._setupSchema(config, name, schema, metadata)
982 
983  def _setupSchema(self, config, name, schema, metadata):
984  # Get a FunctorKey that can quickly look up the "blessed" centroid value.
985  self.centroidKey = afwTable.Point2DKey(schema["slot_Centroid"])
986 
987  # Add some fields for our outputs, and save their Keys.
988  # Use setattr() to programmatically set the pos/neg named attributes to values, e.g.
989  # self.posCentroidKeyX = 'ip_diffim_DipoleFit_pos_centroid_x'
990 
991  for pos_neg in ['pos', 'neg']:
992 
993  key = schema.addField(
994  schema.join(name, pos_neg, "instFlux"), type=float, units="count",
995  doc="Dipole {0} lobe flux".format(pos_neg))
996  setattr(self, ''.join((pos_neg, 'FluxKey')), key)
997 
998  key = schema.addField(
999  schema.join(name, pos_neg, "instFluxErr"), type=float, units="pixel",
1000  doc="1-sigma uncertainty for {0} dipole flux".format(pos_neg))
1001  setattr(self, ''.join((pos_neg, 'FluxErrKey')), key)
1002 
1003  for x_y in ['x', 'y']:
1004  key = schema.addField(
1005  schema.join(name, pos_neg, "centroid", x_y), type=float, units="pixel",
1006  doc="Dipole {0} lobe centroid".format(pos_neg))
1007  setattr(self, ''.join((pos_neg, 'CentroidKey', x_y.upper())), key)
1008 
1009  for x_y in ['x', 'y']:
1010  key = schema.addField(
1011  schema.join(name, "centroid", x_y), type=float, units="pixel",
1012  doc="Dipole centroid")
1013  setattr(self, ''.join(('centroidKey', x_y.upper())), key)
1014 
1015  self.fluxKey = schema.addField(
1016  schema.join(name, "instFlux"), type=float, units="count",
1017  doc="Dipole overall flux")
1018 
1019  self.orientationKey = schema.addField(
1020  schema.join(name, "orientation"), type=float, units="deg",
1021  doc="Dipole orientation")
1022 
1023  self.separationKey = schema.addField(
1024  schema.join(name, "separation"), type=float, units="pixel",
1025  doc="Pixel separation between positive and negative lobes of dipole")
1026 
1027  self.chi2dofKey = schema.addField(
1028  schema.join(name, "chi2dof"), type=float,
1029  doc="Chi2 per degree of freedom of dipole fit")
1030 
1031  self.signalToNoiseKey = schema.addField(
1032  schema.join(name, "signalToNoise"), type=float,
1033  doc="Estimated signal-to-noise of dipole fit")
1034 
1035  self.classificationFlagKey = schema.addField(
1036  schema.join(name, "flag", "classification"), type="Flag",
1037  doc="Flag indicating diaSource is classified as a dipole")
1038 
1039  self.classificationAttemptedFlagKey = schema.addField(
1040  schema.join(name, "flag", "classificationAttempted"), type="Flag",
1041  doc="Flag indicating diaSource was attempted to be classified as a dipole")
1042 
1043  self.flagKey = schema.addField(
1044  schema.join(name, "flag"), type="Flag",
1045  doc="General failure flag for dipole fit")
1046 
1047  self.edgeFlagKey = schema.addField(
1048  schema.join(name, "flag", "edge"), type="Flag",
1049  doc="Flag set when dipole is too close to edge of image")
1050 
1051  def measure(self, measRecord, exposure, posExp=None, negExp=None):
1052  """Perform the non-linear least squares minimization on the putative dipole source.
1053 
1054  Parameters
1055  ----------
1056  measRecord : `lsst.afw.table.SourceRecord`
1057  diaSources that will be measured using dipole measurement
1058  exposure : `lsst.afw.image.Exposure`
1059  Difference exposure on which the diaSources were detected; `exposure = posExp-negExp`
1060  If both `posExp` and `negExp` are `None`, will attempt to fit the
1061  dipole to just the `exposure` with no constraint.
1062  posExp : `lsst.afw.image.Exposure`, optional
1063  "Positive" exposure, typically a science exposure, or None if unavailable
1064  When `posExp` is `None`, will compute `posImage = exposure + negExp`.
1065  negExp : `lsst.afw.image.Exposure`, optional
1066  "Negative" exposure, typically a template exposure, or None if unavailable
1067  When `negExp` is `None`, will compute `negImage = posExp - exposure`.
1068 
1069  Notes
1070  -----
1071  The main functionality of this routine was placed outside of
1072  this plugin (into `DipoleFitAlgorithm.fitDipole()`) so that
1073  `DipoleFitAlgorithm.fitDipole()` can be called separately for
1074  testing (@see `tests/testDipoleFitter.py`)
1075 
1076  Returns
1077  -------
1078  result : TODO: DM-17458
1079  TODO: DM-17458
1080  """
1081 
1082  result = None
1083  pks = measRecord.getFootprint().getPeaks()
1084 
1085  # Check if the footprint consists of a putative dipole - else don't fit it.
1086  if (
1087  (len(pks) <= 1) or # one peak in the footprint - not a dipole
1088  (len(pks) > 1 and (np.sign(pks[0].getPeakValue()) ==
1089  np.sign(pks[-1].getPeakValue()))) # peaks are same sign - not a dipole
1090  ):
1091  measRecord.set(self.classificationFlagKey, False)
1092  measRecord.set(self.classificationAttemptedFlagKey, False)
1093  self.fail(measRecord, measBase.MeasurementError('not a dipole', self.FAILURE_NOT_DIPOLE))
1094  if not self.config.fitAllDiaSources:
1095  return result
1096 
1097  try:
1098  alg = self.DipoleFitAlgorithmClass(exposure, posImage=posExp, negImage=negExp)
1099  result, _ = alg.fitDipole(
1100  measRecord, rel_weight=self.config.relWeight,
1101  tol=self.config.tolerance,
1102  maxSepInSigma=self.config.maxSeparation,
1103  fitBackground=self.config.fitBackground,
1104  separateNegParams=self.config.fitSeparateNegParams,
1105  verbose=False, display=False)
1106  except pexExcept.LengthError:
1107  self.fail(measRecord, measBase.MeasurementError('edge failure', self.FAILURE_EDGE))
1108  except Exception:
1109  self.fail(measRecord, measBase.MeasurementError('dipole fit failure', self.FAILURE_FIT))
1110 
1111  if result is None:
1112  measRecord.set(self.classificationFlagKey, False)
1113  measRecord.set(self.classificationAttemptedFlagKey, False)
1114  return result
1115 
1116  self.log.debug("Dipole fit result: %d %s", measRecord.getId(), str(result))
1117 
1118  if result.posFlux <= 1.: # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit.
1119  self.fail(measRecord, measBase.MeasurementError('dipole fit failure', self.FAILURE_FIT))
1120 
1121  # add chi2, coord/flux uncertainties (TBD), dipole classification
1122  # Add the relevant values to the measRecord
1123  measRecord[self.posFluxKey] = result.posFlux
1124  measRecord[self.posFluxErrKey] = result.signalToNoise # to be changed to actual sigma!
1125  measRecord[self.posCentroidKeyX] = result.posCentroidX
1126  measRecord[self.posCentroidKeyY] = result.posCentroidY
1127 
1128  measRecord[self.negFluxKey] = result.negFlux
1129  measRecord[self.negFluxErrKey] = result.signalToNoise # to be changed to actual sigma!
1130  measRecord[self.negCentroidKeyX] = result.negCentroidX
1131  measRecord[self.negCentroidKeyY] = result.negCentroidY
1132 
1133  # Dia source flux: average of pos+neg
1134  measRecord[self.fluxKey] = (abs(result.posFlux) + abs(result.negFlux))/2.
1135  measRecord[self.orientationKey] = result.orientation
1136  measRecord[self.separationKey] = np.sqrt((result.posCentroidX - result.negCentroidX)**2. +
1137  (result.posCentroidY - result.negCentroidY)**2.)
1138  measRecord[self.centroidKeyX] = result.centroidX
1139  measRecord[self.centroidKeyY] = result.centroidY
1140 
1141  measRecord[self.signalToNoiseKey] = result.signalToNoise
1142  measRecord[self.chi2dofKey] = result.redChi2
1143 
1144  self.doClassify(measRecord, result.chi2)
1145 
1146  def doClassify(self, measRecord, chi2val):
1147  """Classify a source as a dipole.
1148 
1149  Parameters
1150  ----------
1151  measRecord : TODO: DM-17458
1152  TODO: DM-17458
1153  chi2val : TODO: DM-17458
1154  TODO: DM-17458
1155 
1156  Notes
1157  -----
1158  Sources are classified as dipoles, or not, according to three criteria:
1159 
1160  1. Does the total signal-to-noise surpass the ``minSn``?
1161  2. Are the pos/neg fluxes greater than 1.0 and no more than 0.65 (``maxFluxRatio``)
1162  of the total flux? By default this will never happen since ``posFlux == negFlux``.
1163  3. Is it a good fit (``chi2dof`` < 1)? (Currently not used.)
1164  """
1165 
1166  # First, does the total signal-to-noise surpass the minSn?
1167  passesSn = measRecord[self.signalToNoiseKey] > self.config.minSn
1168 
1169  # Second, are the pos/neg fluxes greater than 1.0 and no more than 0.65 (param maxFluxRatio)
1170  # of the total flux? By default this will never happen since posFlux = negFlux.
1171  passesFluxPos = (abs(measRecord[self.posFluxKey]) /
1172  (measRecord[self.fluxKey]*2.)) < self.config.maxFluxRatio
1173  passesFluxPos &= (abs(measRecord[self.posFluxKey]) >= 1.0)
1174  passesFluxNeg = (abs(measRecord[self.negFluxKey]) /
1175  (measRecord[self.fluxKey]*2.)) < self.config.maxFluxRatio
1176  passesFluxNeg &= (abs(measRecord[self.negFluxKey]) >= 1.0)
1177  allPass = (passesSn and passesFluxPos and passesFluxNeg) # and passesChi2)
1178 
1179  # Third, is it a good fit (chi2dof < 1)?
1180  # Use scipy's chi2 cumulative distrib to estimate significance
1181  # This doesn't really work since I don't trust the values in the variance plane (which
1182  # affects the least-sq weights, which affects the resulting chi2).
1183  # But I'm going to keep this here for future use.
1184  if False:
1185  from scipy.stats import chi2
1186  ndof = chi2val / measRecord[self.chi2dofKey]
1187  significance = chi2.cdf(chi2val, ndof)
1188  passesChi2 = significance < self.config.maxChi2DoF
1189  allPass = allPass and passesChi2
1190 
1191  measRecord.set(self.classificationAttemptedFlagKey, True)
1192 
1193  if allPass: # Note cannot pass `allPass` into the `measRecord.set()` call below...?
1194  measRecord.set(self.classificationFlagKey, True)
1195  else:
1196  measRecord.set(self.classificationFlagKey, False)
1197 
1198  def fail(self, measRecord, error=None):
1199  """Catch failures and set the correct flags.
1200  """
1201 
1202  measRecord.set(self.flagKey, True)
1203  if error is not None:
1204  if error.getFlagBit() == self.FAILURE_EDGE:
1205  self.log.warn('DipoleFitPlugin not run on record %d: %s', measRecord.getId(), str(error))
1206  measRecord.set(self.edgeFlagKey, True)
1207  if error.getFlagBit() == self.FAILURE_FIT:
1208  self.log.warn('DipoleFitPlugin failed on record %d: %s', measRecord.getId(), str(error))
1209  measRecord.set(self.flagKey, True)
1210  if error.getFlagBit() == self.FAILURE_NOT_DIPOLE:
1211  self.log.debug('DipoleFitPlugin not run on record %d: %s',
1212  measRecord.getId(), str(error))
1213  measRecord.set(self.classificationAttemptedFlagKey, False)
1214  measRecord.set(self.flagKey, True)
1215  else:
1216  self.log.warn('DipoleFitPlugin failed on record %d', measRecord.getId())
def _getHeavyFootprintSubimage(self, fp, badfill=np.nan, grow=0)
Angle abs(Angle const &a)
Definition: Angle.h:106
def 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)
def makeStarModel(self, bbox, psf, xcen, ycen, flux)
def fitDipole(self, source, tol=1e-7, rel_weight=0.1, fitBackground=1, maxSepInSigma=5., separateNegParams=True, bgGradientOrder=1, verbose=False, display=False)
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
def __init__(self, diffim, posImage=None, negImage=None)
def displayFitResults(self, footprint, result)
Definition: Log.h:691
def run(self, sources, exposure, posExp=None, negExp=None, kwargs)
def fitFootprintBackground(self, source, posImage, order=1)
def __init__(self, config, name, schema, metadata)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
def makeBackgroundModel(self, in_x, pars=None)
def __init__(self, schema, algMetadata=None, kwargs)
def _setupSchema(self, config, name, schema, metadata)
def measure(mi, x, y, size, statistic, stats)
Definition: fringe.py:357
def fitDipoleImpl(self, source, tol=1e-7, rel_weight=0.5, fitBackground=1, bgGradientOrder=1, maxSepInSigma=5., separateNegParams=True, verbose=False)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
def doClassify(self, measRecord, chi2val)
def fail(self, measRecord, error=None)
def measure(self, measRecord, exposure, posExp=None, negExp=None)