LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
Public Member Functions | Public Attributes | List of all members
lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm Class Reference
Inheritance diagram for lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm:

Public Member Functions

def __init__ (self, diffim, posImage=None, negImage=None)
 
def fitDipoleImpl (self, source, tol=1e-7, rel_weight=0.5, fitBackground=1, bgGradientOrder=1, maxSepInSigma=5., separateNegParams=True, verbose=False)
 
def fitDipole (self, source, tol=1e-7, rel_weight=0.1, fitBackground=1, maxSepInSigma=5., separateNegParams=True, bgGradientOrder=1, verbose=False, display=False)
 
def displayFitResults (self, footprint, result)
 

Public Attributes

 diffim
 
 posImage
 
 negImage
 
 psfSigma
 
 log
 
 debug
 

Detailed Description

Fit a dipole model using an image difference.

See also:
`DMTN-007: Dipole characterization for image differencing  <https://dmtn-007.lsst.io>`_.

Definition at line 483 of file dipoleFitTask.py.

Constructor & Destructor Documentation

◆ __init__()

def lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.__init__ (   self,
  diffim,
  posImage = None,
  negImage = None 
)
Algorithm to run dipole measurement on a diaSource

Parameters
----------
diffim : `lsst.afw.image.Exposure`
    Exposure on which the diaSources were detected
posImage : `lsst.afw.image.Exposure`
    "Positive" exposure from which the template was subtracted
negImage : `lsst.afw.image.Exposure`
    "Negative" exposure which was subtracted from the posImage

Definition at line 509 of file dipoleFitTask.py.

509  def __init__(self, diffim, posImage=None, negImage=None):
510  """Algorithm to run dipole measurement on a diaSource
511 
512  Parameters
513  ----------
514  diffim : `lsst.afw.image.Exposure`
515  Exposure on which the diaSources were detected
516  posImage : `lsst.afw.image.Exposure`
517  "Positive" exposure from which the template was subtracted
518  negImage : `lsst.afw.image.Exposure`
519  "Negative" exposure which was subtracted from the posImage
520  """
521 
522  self.diffim = diffim
523  self.posImage = posImage
524  self.negImage = negImage
525  self.psfSigma = None
526  if diffim is not None:
527  self.psfSigma = diffim.getPsf().computeShape().getDeterminantRadius()
528 
529  self.log = Log.getLogger(__name__)
530 
531  import lsstDebug
532  self.debug = lsstDebug.Info(__name__).debug
533 

Member Function Documentation

◆ displayFitResults()

def lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.displayFitResults (   self,
  footprint,
  result 
)
Display data, model fits and residuals (currently uses matplotlib display functions).

Parameters
----------
footprint : TODO: DM-17458
    Footprint containing the dipole that was fit
result : `lmfit.MinimizerResult`
    `lmfit.MinimizerResult` object returned by `lmfit` optimizer

Returns
-------
fig : `matplotlib.pyplot.plot`

Definition at line 884 of file dipoleFitTask.py.

884  def displayFitResults(self, footprint, result):
885  """Display data, model fits and residuals (currently uses matplotlib display functions).
886 
887  Parameters
888  ----------
889  footprint : TODO: DM-17458
890  Footprint containing the dipole that was fit
891  result : `lmfit.MinimizerResult`
892  `lmfit.MinimizerResult` object returned by `lmfit` optimizer
893 
894  Returns
895  -------
896  fig : `matplotlib.pyplot.plot`
897  """
898  try:
899  import matplotlib.pyplot as plt
900  except ImportError as err:
901  self.log.warning('Unable to import matplotlib: %s', err)
902  raise err
903 
904  def display2dArray(arr, title='Data', extent=None):
905  """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range.
906  """
907  fig = plt.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent)
908  plt.title(title)
909  plt.colorbar(fig, cmap='gray')
910  return fig
911 
912  z = result.data
913  fit = result.best_fit
914  bbox = footprint.getBBox()
915  extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY())
916  if z.shape[0] == 3:
917  fig = plt.figure(figsize=(8, 8))
918  for i in range(3):
919  plt.subplot(3, 3, i*3+1)
920  display2dArray(z[i, :], 'Data', extent=extent)
921  plt.subplot(3, 3, i*3+2)
922  display2dArray(fit[i, :], 'Model', extent=extent)
923  plt.subplot(3, 3, i*3+3)
924  display2dArray(z[i, :] - fit[i, :], 'Residual', extent=extent)
925  return fig
926  else:
927  fig = plt.figure(figsize=(8, 2.5))
928  plt.subplot(1, 3, 1)
929  display2dArray(z, 'Data', extent=extent)
930  plt.subplot(1, 3, 2)
931  display2dArray(fit, 'Model', extent=extent)
932  plt.subplot(1, 3, 3)
933  display2dArray(z - fit, 'Residual', extent=extent)
934  return fig
935 
936  plt.show()
937 
938 
939 @measBase.register("ip_diffim_DipoleFit")

◆ fitDipole()

def lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.fitDipole (   self,
  source,
  tol = 1e-7,
  rel_weight = 0.1,
  fitBackground = 1,
  maxSepInSigma = 5.,
  separateNegParams = True,
  bgGradientOrder = 1,
  verbose = False,
  display = False 
)
Fit a dipole model to an input ``diaSource`` (wraps `fitDipoleImpl`).

Actually, fits the subimage bounded by the input source's
footprint) and optionally constrain the fit using the
pre-subtraction images self.posImage (science) and
self.negImage (template). Wraps the output into a
`pipeBase.Struct` named tuple after computing additional
statistics such as orientation and SNR.

Parameters
----------
source : `lsst.afw.table.SourceRecord`
    Record containing the (merged) dipole source footprint detected on the diffim
tol : `float`, optional
    Tolerance parameter for scipy.leastsq() optimization
rel_weight : `float`, optional
    Weighting of posImage/negImage relative to the diffim in the fit
fitBackground : `int`, {0, 1, 2}, optional
    How to fit linear background gradient in posImage/negImage

        - 0: do not fit background at all
        - 1 (default): pre-fit the background using linear least squares and then do not fit it
          as part of the dipole fitting optimization
        - 2: pre-fit the background using linear least squares (as in 1), and use the parameter
          estimates from that fit as starting parameters for an integrated "re-fit" of the
          background as part of the overall dipole fitting optimization.
maxSepInSigma : `float`, optional
    Allowed window of centroid parameters relative to peak in input source footprint
separateNegParams : `bool`, optional
    Fit separate parameters to the flux and background gradient in
bgGradientOrder : `int`, {0, 1, 2}, optional
    Desired polynomial order of background gradient
verbose: `bool`, optional
    Be verbose
display
    Display input data, best fit model(s) and residuals in a matplotlib window.

Returns
-------
result : `struct`
    `pipeBase.Struct` object containing the fit parameters and other information.

result : `callable`
    `lmfit.MinimizerResult` object for debugging and error estimation, etc.

Notes
-----
Parameter `fitBackground` has three options, thus it is an integer:

Definition at line 771 of file dipoleFitTask.py.

773  bgGradientOrder=1, verbose=False, display=False):
774  """Fit a dipole model to an input ``diaSource`` (wraps `fitDipoleImpl`).
775 
776  Actually, fits the subimage bounded by the input source's
777  footprint) and optionally constrain the fit using the
778  pre-subtraction images self.posImage (science) and
779  self.negImage (template). Wraps the output into a
780  `pipeBase.Struct` named tuple after computing additional
781  statistics such as orientation and SNR.
782 
783  Parameters
784  ----------
785  source : `lsst.afw.table.SourceRecord`
786  Record containing the (merged) dipole source footprint detected on the diffim
787  tol : `float`, optional
788  Tolerance parameter for scipy.leastsq() optimization
789  rel_weight : `float`, optional
790  Weighting of posImage/negImage relative to the diffim in the fit
791  fitBackground : `int`, {0, 1, 2}, optional
792  How to fit linear background gradient in posImage/negImage
793 
794  - 0: do not fit background at all
795  - 1 (default): pre-fit the background using linear least squares and then do not fit it
796  as part of the dipole fitting optimization
797  - 2: pre-fit the background using linear least squares (as in 1), and use the parameter
798  estimates from that fit as starting parameters for an integrated "re-fit" of the
799  background as part of the overall dipole fitting optimization.
800  maxSepInSigma : `float`, optional
801  Allowed window of centroid parameters relative to peak in input source footprint
802  separateNegParams : `bool`, optional
803  Fit separate parameters to the flux and background gradient in
804  bgGradientOrder : `int`, {0, 1, 2}, optional
805  Desired polynomial order of background gradient
806  verbose: `bool`, optional
807  Be verbose
808  display
809  Display input data, best fit model(s) and residuals in a matplotlib window.
810 
811  Returns
812  -------
813  result : `struct`
814  `pipeBase.Struct` object containing the fit parameters and other information.
815 
816  result : `callable`
817  `lmfit.MinimizerResult` object for debugging and error estimation, etc.
818 
819  Notes
820  -----
821  Parameter `fitBackground` has three options, thus it is an integer:
822 
823  """
824 
825  fitResult = self.fitDipoleImpl(
826  source, tol=tol, rel_weight=rel_weight, fitBackground=fitBackground,
827  maxSepInSigma=maxSepInSigma, separateNegParams=separateNegParams,
828  bgGradientOrder=bgGradientOrder, verbose=verbose)
829 
830  # Display images, model fits and residuals (currently uses matplotlib display functions)
831  if display:
832  fp = source.getFootprint()
833  self.displayFitResults(fp, fitResult)
834 
835  fitParams = fitResult.best_values
836  if fitParams['flux'] <= 1.: # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit.
837  out = Struct(posCentroidX=np.nan, posCentroidY=np.nan,
838  negCentroidX=np.nan, negCentroidY=np.nan,
839  posFlux=np.nan, negFlux=np.nan, posFluxErr=np.nan, negFluxErr=np.nan,
840  centroidX=np.nan, centroidY=np.nan, orientation=np.nan,
841  signalToNoise=np.nan, chi2=np.nan, redChi2=np.nan)
842  return out, fitResult
843 
844  centroid = ((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2.,
845  (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2.)
846  dx, dy = fitParams['xcenPos'] - fitParams['xcenNeg'], fitParams['ycenPos'] - fitParams['ycenNeg']
847  angle = np.arctan2(dy, dx) / np.pi * 180. # convert to degrees (should keep as rad?)
848 
849  # Exctract flux value, compute signalToNoise from flux/variance_within_footprint
850  # Also extract the stderr of flux estimate.
851  def computeSumVariance(exposure, footprint):
852  box = footprint.getBBox()
853  subim = afwImage.MaskedImageF(exposure.getMaskedImage(), box, origin=afwImage.PARENT)
854  return np.sqrt(np.nansum(subim.getArrays()[1][:, :]))
855 
856  fluxVal = fluxVar = fitParams['flux']
857  fluxErr = fluxErrNeg = fitResult.params['flux'].stderr
858  if self.posImage is not None:
859  fluxVar = computeSumVariance(self.posImage, source.getFootprint())
860  else:
861  fluxVar = computeSumVariance(self.diffim, source.getFootprint())
862 
863  fluxValNeg, fluxVarNeg = fluxVal, fluxVar
864  if separateNegParams:
865  fluxValNeg = fitParams['fluxNeg']
866  fluxErrNeg = fitResult.params['fluxNeg'].stderr
867  if self.negImage is not None:
868  fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint())
869 
870  try:
871  signalToNoise = np.sqrt((fluxVal/fluxVar)**2 + (fluxValNeg/fluxVarNeg)**2)
872  except ZeroDivisionError: # catch divide by zero - should never happen.
873  signalToNoise = np.nan
874 
875  out = Struct(posCentroidX=fitParams['xcenPos'], posCentroidY=fitParams['ycenPos'],
876  negCentroidX=fitParams['xcenNeg'], negCentroidY=fitParams['ycenNeg'],
877  posFlux=fluxVal, negFlux=-fluxValNeg, posFluxErr=fluxErr, negFluxErr=fluxErrNeg,
878  centroidX=centroid[0], centroidY=centroid[1], orientation=angle,
879  signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi)
880 
881  # fitResult may be returned for debugging
882  return out, fitResult
883 

◆ fitDipoleImpl()

def lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.fitDipoleImpl (   self,
  source,
  tol = 1e-7,
  rel_weight = 0.5,
  fitBackground = 1,
  bgGradientOrder = 1,
  maxSepInSigma = 5.,
  separateNegParams = True,
  verbose = False 
)
Fit a dipole model to an input difference image.

Actually, fits the subimage bounded by the input source's
footprint) and optionally constrain the fit using the
pre-subtraction images posImage and negImage.

Parameters
----------
source : TODO: DM-17458
    TODO: DM-17458
tol : float, optional
    TODO: DM-17458
rel_weight : `float`, optional
    TODO: DM-17458
fitBackground : `int`, optional
    TODO: DM-17458
bgGradientOrder : `int`, optional
    TODO: DM-17458
maxSepInSigma : `float`, optional
    TODO: DM-17458
separateNegParams : `bool`, optional
    TODO: DM-17458
verbose : `bool`, optional
    TODO: DM-17458

Returns
-------
result : `lmfit.MinimizerResult`
    return `lmfit.MinimizerResult` object containing the fit
    parameters and other information.

Definition at line 534 of file dipoleFitTask.py.

536  separateNegParams=True, verbose=False):
537  """Fit a dipole model to an input difference image.
538 
539  Actually, fits the subimage bounded by the input source's
540  footprint) and optionally constrain the fit using the
541  pre-subtraction images posImage and negImage.
542 
543  Parameters
544  ----------
545  source : TODO: DM-17458
546  TODO: DM-17458
547  tol : float, optional
548  TODO: DM-17458
549  rel_weight : `float`, optional
550  TODO: DM-17458
551  fitBackground : `int`, optional
552  TODO: DM-17458
553  bgGradientOrder : `int`, optional
554  TODO: DM-17458
555  maxSepInSigma : `float`, optional
556  TODO: DM-17458
557  separateNegParams : `bool`, optional
558  TODO: DM-17458
559  verbose : `bool`, optional
560  TODO: DM-17458
561 
562  Returns
563  -------
564  result : `lmfit.MinimizerResult`
565  return `lmfit.MinimizerResult` object containing the fit
566  parameters and other information.
567  """
568 
569  # Only import lmfit if someone wants to use the new DipoleFitAlgorithm.
570  import lmfit
571 
572  fp = source.getFootprint()
573  bbox = fp.getBBox()
574  subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)
575 
576  z = diArr = subim.getArrays()[0]
577  weights = 1. / subim.getArrays()[2] # get the weights (=1/variance)
578 
579  if rel_weight > 0. and ((self.posImage is not None) or (self.negImage is not None)):
580  if self.negImage is not None:
581  negSubim = afwImage.MaskedImageF(self.negImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
582  if self.posImage is not None:
583  posSubim = afwImage.MaskedImageF(self.posImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
584  if self.posImage is None: # no science image provided; generate it from diffim + negImage
585  posSubim = subim.clone()
586  posSubim += negSubim
587  if self.negImage is None: # no template provided; generate it from the posImage - diffim
588  negSubim = posSubim.clone()
589  negSubim -= subim
590 
591  z = np.append([z], [posSubim.getArrays()[0],
592  negSubim.getArrays()[0]], axis=0)
593  # Weight the pos/neg images by rel_weight relative to the diffim
594  weights = np.append([weights], [1. / posSubim.getArrays()[2] * rel_weight,
595  1. / negSubim.getArrays()[2] * rel_weight], axis=0)
596  else:
597  rel_weight = 0. # a short-cut for "don't include the pre-subtraction data"
598 
599  # It seems that `lmfit` requires a static functor as its optimized method, which eliminates
600  # the ability to pass a bound method or other class method. Here we write a wrapper which
601  # makes this possible.
602  def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None,
603  b=None, x1=None, y1=None, xy=None, x2=None, y2=None,
604  bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None,
605  **kwargs):
606  """Generate dipole model with given parameters.
607 
608  It simply defers to `modelObj.makeModel()`, where `modelObj` comes
609  out of `kwargs['modelObj']`.
610  """
611  modelObj = kwargs.pop('modelObj')
612  return modelObj.makeModel(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=fluxNeg,
613  b=b, x1=x1, y1=y1, xy=xy, x2=x2, y2=y2,
614  bNeg=bNeg, x1Neg=x1Neg, y1Neg=y1Neg, xyNeg=xyNeg,
615  x2Neg=x2Neg, y2Neg=y2Neg, **kwargs)
616 
617  dipoleModel = DipoleModel()
618 
619  modelFunctor = dipoleModelFunctor # dipoleModel.makeModel does not work for now.
620  # Create the lmfit model (lmfit uses scipy 'leastsq' option by default - Levenberg-Marquardt)
621  # Note we can also tell it to drop missing values from the data.
622  gmod = lmfit.Model(modelFunctor, verbose=verbose, missing='drop')
623  # independent_vars=independent_vars) #, param_names=param_names)
624 
625  # Add the constraints for centroids, fluxes.
626  # starting constraint - near centroid of footprint
627  fpCentroid = np.array([fp.getCentroid().getX(), fp.getCentroid().getY()])
628  cenNeg = cenPos = fpCentroid
629 
630  pks = fp.getPeaks()
631 
632  if len(pks) >= 1:
633  cenPos = pks[0].getF() # if individual (merged) peaks were detected, use those
634  if len(pks) >= 2: # peaks are already sorted by centroid flux so take the most negative one
635  cenNeg = pks[-1].getF()
636 
637  # For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
638  # First assume dipole is not separated by more than 5*psfSigma.
639  maxSep = self.psfSigma * maxSepInSigma
640 
641  # As an initial guess -- assume the dipole is close to the center of the footprint.
642  if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
643  cenPos = fpCentroid
644  if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
645  cenPos = fpCentroid
646 
647  # parameter hints/constraints: https://lmfit.github.io/lmfit-py/model.html#model-param-hints-section
648  # might make sense to not use bounds -- see http://lmfit.github.io/lmfit-py/bounds.html
649  # also see this discussion -- https://github.com/scipy/scipy/issues/3129
650  gmod.set_param_hint('xcenPos', value=cenPos[0],
651  min=cenPos[0]-maxSep, max=cenPos[0]+maxSep)
652  gmod.set_param_hint('ycenPos', value=cenPos[1],
653  min=cenPos[1]-maxSep, max=cenPos[1]+maxSep)
654  gmod.set_param_hint('xcenNeg', value=cenNeg[0],
655  min=cenNeg[0]-maxSep, max=cenNeg[0]+maxSep)
656  gmod.set_param_hint('ycenNeg', value=cenNeg[1],
657  min=cenNeg[1]-maxSep, max=cenNeg[1]+maxSep)
658 
659  # Use the (flux under the dipole)*5 for an estimate.
660  # Lots of testing showed that having startingFlux be too high was better than too low.
661  startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
662  posFlux = negFlux = startingFlux
663 
664  # TBD: set max. flux limit?
665  gmod.set_param_hint('flux', value=posFlux, min=0.1)
666 
667  if separateNegParams:
668  # TBD: set max negative lobe flux limit?
669  gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
670 
671  # Fixed parameters (don't fit for them if there are no pre-sub images or no gradient fit requested):
672  # Right now (fitBackground == 1), we fit a linear model to the background and then subtract
673  # it from the data and then don't fit the background again (this is faster).
674  # A slower alternative (fitBackground == 2) is to use the estimated background parameters as
675  # starting points in the integrated model fit. That is currently not performed by default,
676  # but might be desirable in some cases.
677  bgParsPos = bgParsNeg = (0., 0., 0.)
678  if ((rel_weight > 0.) and (fitBackground != 0) and (bgGradientOrder >= 0)):
679  pbg = 0.
680  bgFitImage = self.posImage if self.posImage is not None else self.negImage
681  # Fit the gradient to the background (linear model)
682  bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
683  order=bgGradientOrder)
684 
685  # Generate the gradient and subtract it from the pre-subtraction image data
686  if fitBackground == 1:
687  in_x = dipoleModel._generateXYGrid(bbox)
688  pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsPos))
689  z[1, :] -= pbg
690  z[1, :] -= np.nanmedian(z[1, :])
691  posFlux = np.nansum(z[1, :])
692  gmod.set_param_hint('flux', value=posFlux*1.5, min=0.1)
693 
694  if separateNegParams and self.negImage is not None:
695  bgParsNeg = dipoleModel.fitFootprintBackground(source, self.negImage,
696  order=bgGradientOrder)
697  pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsNeg))
698  z[2, :] -= pbg
699  z[2, :] -= np.nanmedian(z[2, :])
700  if separateNegParams:
701  negFlux = np.nansum(z[2, :])
702  gmod.set_param_hint('fluxNeg', value=negFlux*1.5, min=0.1)
703 
704  # Do not subtract the background from the images but include the background parameters in the fit
705  if fitBackground == 2:
706  if bgGradientOrder >= 0:
707  gmod.set_param_hint('b', value=bgParsPos[0])
708  if separateNegParams:
709  gmod.set_param_hint('bNeg', value=bgParsNeg[0])
710  if bgGradientOrder >= 1:
711  gmod.set_param_hint('x1', value=bgParsPos[1])
712  gmod.set_param_hint('y1', value=bgParsPos[2])
713  if separateNegParams:
714  gmod.set_param_hint('x1Neg', value=bgParsNeg[1])
715  gmod.set_param_hint('y1Neg', value=bgParsNeg[2])
716  if bgGradientOrder >= 2:
717  gmod.set_param_hint('xy', value=bgParsPos[3])
718  gmod.set_param_hint('x2', value=bgParsPos[4])
719  gmod.set_param_hint('y2', value=bgParsPos[5])
720  if separateNegParams:
721  gmod.set_param_hint('xyNeg', value=bgParsNeg[3])
722  gmod.set_param_hint('x2Neg', value=bgParsNeg[4])
723  gmod.set_param_hint('y2Neg', value=bgParsNeg[5])
724 
725  y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
726  in_x = np.array([x, y]).astype(np.float)
727  in_x[0, :] -= in_x[0, :].mean() # center it!
728  in_x[1, :] -= in_x[1, :].mean()
729 
730  # Instead of explicitly using a mask to ignore flagged pixels, just set the ignored pixels'
731  # weights to 0 in the fit. TBD: need to inspect mask planes to set this mask.
732  mask = np.ones_like(z, dtype=bool) # TBD: set mask values to False if the pixels are to be ignored
733 
734  # I'm not sure about the variance planes in the diffim (or convolved pre-sub. images
735  # for that matter) so for now, let's just do an un-weighted least-squares fit
736  # (override weights computed above).
737  weights = mask.astype(np.float64)
738  if self.posImage is not None and rel_weight > 0.:
739  weights = np.array([np.ones_like(diArr), np.ones_like(diArr)*rel_weight,
740  np.ones_like(diArr)*rel_weight])
741 
742  # Set the weights to zero if mask is False
743  if np.any(~mask):
744  weights[~mask] = 0.
745 
746  # Note that although we can, we're not required to set initial values for params here,
747  # since we set their param_hint's above.
748  # Can add "method" param to not use 'leastsq' (==levenberg-marquardt), e.g. "method='nelder'"
749  with warnings.catch_warnings():
750  warnings.simplefilter("ignore") # temporarily turn off silly lmfit warnings
751  result = gmod.fit(z, weights=weights, x=in_x,
752  verbose=verbose,
753  fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
754  'maxfev': 250}, # see scipy docs
755  psf=self.diffim.getPsf(), # hereon: kwargs that get passed to genDipoleModel()
756  rel_weight=rel_weight,
757  footprint=fp,
758  modelObj=dipoleModel)
759 
760  if verbose: # the ci_report() seems to fail if neg params are constrained -- TBD why.
761  # Never wanted in production - this takes a long time (longer than the fit!)
762  # This is how to get confidence intervals out:
763  # https://lmfit.github.io/lmfit-py/confidence.html and
764  # http://cars9.uchicago.edu/software/python/lmfit/model.html
765  print(result.fit_report(show_correl=False))
766  if separateNegParams:
767  print(result.ci_report())
768 
769  return result
770 

Member Data Documentation

◆ debug

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.debug

Definition at line 532 of file dipoleFitTask.py.

◆ diffim

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.diffim

Definition at line 522 of file dipoleFitTask.py.

◆ log

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.log

Definition at line 529 of file dipoleFitTask.py.

◆ negImage

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.negImage

Definition at line 524 of file dipoleFitTask.py.

◆ posImage

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.posImage

Definition at line 523 of file dipoleFitTask.py.

◆ psfSigma

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.psfSigma

Definition at line 525 of file dipoleFitTask.py.


The documentation for this class was generated from the following file: