lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm Class Reference

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

Public Attributes
	diffim = diffim
	posImage = posImage
	negImage = negImage
	psfSigma = None
	log = logging.getLogger(__name__)
	debug = lsstDebug.Info(__name__).debug

Static Protected Attributes
str	_private_version_ = '0.0.5'

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 490 of file dipoleFitTask.py.

Constructor & Destructor Documentation

__init__()

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 516 of file dipoleFitTask.py.

    def __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
        """
 
        self.diffim = diffim
        self.posImage = posImage
        self.negImage = negImage
        self.psfSigma = None
        if diffim is not None:
            diffimPsf = diffim.getPsf()
            diffimAvgPos = diffimPsf.getAveragePosition()
            self.psfSigma = diffimPsf.computeShape(diffimAvgPos).getDeterminantRadius()
 
        self.log = logging.getLogger(__name__)
 
        import lsstDebug
        self.debug = lsstDebug.Info(__name__).debug
 

Member Function Documentation

displayFitResults()

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.displayFitResults	(		self,
			footprint,
			result )

Display data, model fits and residuals (currently uses matplotlib display functions).

Parameters
----------
footprint : `lsst.afw.detection.Footprint`
    Footprint containing the dipole that was fit
result : `lmfit.MinimizerResult`
    `lmfit.MinimizerResult` object returned by `lmfit` optimizer

Definition at line 933 of file dipoleFitTask.py.

    def displayFitResults(self, footprint, result):
        """Display data, model fits and residuals (currently uses matplotlib display functions).
 
        Parameters
        ----------
        footprint : `lsst.afw.detection.Footprint`
            Footprint containing the dipole that was fit
        result : `lmfit.MinimizerResult`
            `lmfit.MinimizerResult` object returned by `lmfit` optimizer
        """
        try:
            import matplotlib.pyplot as plt
        except ImportError as err:
            self.log.warning('Unable to import matplotlib: %s', err)
            raise err
 
        def display2dArray(ax, arr, x, y, xErr, yErr, title, extent=None):
            """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range.
            """
            fig = ax.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent)
            ax.set_title(title)
            ax.errorbar(x["total"], y["total"], xErr["total"], yErr["total"], c="cyan")
            ax.errorbar(x["Pos"], y["Pos"], xErr["Pos"], yErr["Pos"], c="green")
            ax.errorbar(x["Neg"], y["Neg"], xErr["Neg"], yErr["Neg"], c="red")
            return fig
 
        z = result.data
        fit = result.best_fit
        bbox = footprint.getBBox()
        extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY())
 
        if z.shape[0] == 3:
            x, y, xErr, yErr = {}, {}, {}, {}
            for name in ("Pos", "Neg"):
                x[name] = result.best_values[f"xcen{name}"]
                y[name] = result.best_values[f"ycen{name}"]
                xErr[name] = result.params[f"xcen{name}"].stderr
                yErr[name] = result.params[f"ycen{name}"].stderr
            x["total"] = (result.best_values["xcenPos"] + result.best_values["xcenNeg"])/2
            y["total"] = (result.best_values["ycenPos"] + result.best_values["ycenNeg"])/2
            xErr["total"] = math.sqrt(xErr["Pos"]**2 + xErr["Neg"]**2)
            yErr["total"] = math.sqrt(yErr["Pos"]**2 + yErr["Neg"]**2)
 
            fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True, figsize=(8, 8))
            for i, label in enumerate(("total", "Pos", "Neg")):
                display2dArray(axes[i][0], z[i, :], x, y, xErr, yErr,
                               f'Data {label}', extent=extent)
                display2dArray(axes[i][1], fit[i, :], x, y, xErr, yErr,
                               f'Model {label}', extent=extent)
                display2dArray(axes[i][2], z[i, :] - fit[i, :], x, y, xErr, yErr,
                               f'Residual {label}', extent=extent)
 
                plt.setp(axes[i][1].get_yticklabels(), visible=False)
                plt.setp(axes[i][2].get_yticklabels(), visible=False)
                if i != 2:  # remove top two row x-axis labels
                    plt.setp(axes[i][0].get_xticklabels(), visible=False)
                    plt.setp(axes[i][1].get_xticklabels(), visible=False)
                    plt.setp(axes[i][2].get_xticklabels(), visible=False)
        else:
            fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True, figsize=(10, 2.5))
            display2dArray(axes[0], z, 'Data', extent=extent)
            display2dArray(axes[1], 'Model', extent=extent)
            display2dArray(axes[2], z - fit, 'Residual', extent=extent)
 
        fig.tight_layout(pad=0, w_pad=0, h_pad=0)
        plt.show()
 
 
@measBase.register("ip_diffim_DipoleFit")

fitDipole()

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 799 of file dipoleFitTask.py.

                  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:
 
        """
 
        fitResult = self.fitDipoleImpl(
            source, tol=tol, rel_weight=rel_weight, fitBackground=fitBackground,
            maxSepInSigma=maxSepInSigma, separateNegParams=separateNegParams,
            bgGradientOrder=bgGradientOrder, verbose=verbose)
 
        # Display images, model fits and residuals (currently uses matplotlib display functions)
        if display:
            fp = source.getFootprint()
            self.displayFitResults(fp, fitResult)
 
        # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit.
        if fitResult.params['flux'].value <= 1.:
            self.log.debug("Fitted flux too small for id=%d; ModelResult.message='%s'",
                           source["id"], fitResult.message)
            return None, fitResult
        if not fitResult.result.errorbars:
            self.log.debug("Could not estimate error bars for id=%d; ModelResult.message='%s'",
                           source["id"], fitResult.message)
            return None, fitResult
 
        # TODO: We could include covariances, which could be derived from
        # `fitResult.params[name].correl`, but those are correlations.
        posCentroid = measBase.CentroidResult(fitResult.params['xcenPos'].value,
                                              fitResult.params['ycenPos'].value,
                                              fitResult.params['xcenPos'].stderr,
                                              fitResult.params['ycenPos'].stderr)
        negCentroid = measBase.CentroidResult(fitResult.params['xcenNeg'].value,
                                              fitResult.params['ycenNeg'].value,
                                              fitResult.params['xcenNeg'].stderr,
                                              fitResult.params['ycenNeg'].stderr)
        xposIdx = fitResult.var_names.index("xcenPos")
        yposIdx = fitResult.var_names.index("ycenPos")
        xnegIdx = fitResult.var_names.index("xcenNeg")
        ynegIdx = fitResult.var_names.index("ycenNeg")
        centroid = measBase.CentroidResult((fitResult.params['xcenPos'] + fitResult.params['xcenNeg']) / 2,
                                           (fitResult.params['ycenPos'] + fitResult.params['ycenNeg']) / 2.,
                                           math.sqrt(posCentroid.xErr**2 + negCentroid.xErr**2
                                                     + 2*fitResult.covar[xposIdx, xnegIdx]) / 2.,
                                           math.sqrt(posCentroid.yErr**2 + negCentroid.yErr**2
                                                     + 2*fitResult.covar[yposIdx, ynegIdx]) / 2.)
        dx = fitResult.params['xcenPos'].value - fitResult.params['xcenNeg'].value
        dy = fitResult.params['ycenPos'].value - fitResult.params['ycenNeg'].value
        angle = np.arctan2(dy, dx)
 
        # Extract flux value, compute signalToNoise from flux/variance_within_footprint
        # Also extract the stderr of flux estimate.
        # TODO: should this instead use the lmfit-computed uncertainty from
        # `lmfitResult.result.uvars['flux'].std_dev`?
        def computeSumVariance(exposure, footprint):
            return math.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array))
 
        # NOTE: These will all be the same unless separateNegParams=True!
        flux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr)
        posFlux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr)
        negFlux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr)
        if self.posImage is not None:
            fluxVar = computeSumVariance(self.posImage, source.getFootprint())
        else:
            fluxVar = computeSumVariance(self.diffim, source.getFootprint())
        fluxVarNeg = fluxVar
 
        if separateNegParams:
            negFlux.instFlux = fitResult.params['fluxNeg'].value
            negFlux.instFluxErr = fitResult.params['fluxNeg'].stderr
        if self.negImage is not None:
            fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint())
 
        try:
            signalToNoise = math.sqrt((posFlux.instFlux/fluxVar)**2 + (negFlux.instFlux/fluxVarNeg)**2)
        except ZeroDivisionError:  # catch divide by zero - should never happen.
            signalToNoise = np.nan
 
        out = Struct(posCentroid=posCentroid, negCentroid=negCentroid, centroid=centroid,
                     posFlux=posFlux, negFlux=negFlux, flux=flux, orientation=angle,
                     signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi,
                     nData=fitResult.ndata)
 
        # fitResult may be returned for debugging
        return out, fitResult
 

fitDipoleImpl()

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 543 of file dipoleFitTask.py.

                      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.
        """
 
        # Only import lmfit if someone wants to use the new DipoleFitAlgorithm.
        import lmfit
 
        fp = source.getFootprint()
        bbox = fp.getBBox()
        subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)
 
        z = diArr = subim.image.array
        # Make sure we don't overwrite buffers.
        z = z.copy()
        weights = 1. / subim.variance.array  # get the weights (=1/variance)
 
        if rel_weight > 0. and ((self.posImage is not None) or (self.negImage is not None)):
            if self.negImage is not None:
                negSubim = afwImage.MaskedImageF(self.negImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
            if self.posImage is not None:
                posSubim = afwImage.MaskedImageF(self.posImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
            if self.posImage is None:  # no science image provided; generate it from diffim + negImage
                posSubim = subim.clone()
                posSubim += negSubim
            if self.negImage is None:  # no template provided; generate it from the posImage - diffim
                negSubim = posSubim.clone()
                negSubim -= subim
 
            z = np.append([z], [posSubim.image.array,
                                negSubim.image.array], axis=0)
            # Weight the pos/neg images by rel_weight relative to the diffim
            weights = np.append([weights], [1. / posSubim.variance.array * rel_weight,
                                            1. / negSubim.variance.array * rel_weight], axis=0)
        else:
            rel_weight = 0.  # a short-cut for "don't include the pre-subtraction data"
 
        # It seems that `lmfit` requires a static functor as its optimized method, which eliminates
        # the ability to pass a bound method or other class method. Here we write a wrapper which
        # makes this possible.
        def dipoleModelFunctor(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):
            """Generate dipole model with given parameters.
 
            It simply defers to `modelObj.makeModel()`, where `modelObj` comes
            out of `kwargs['modelObj']`.
            """
            modelObj = kwargs.pop('modelObj')
            return modelObj.makeModel(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=fluxNeg,
                                      b=b, x1=x1, y1=y1, xy=xy, x2=x2, y2=y2,
                                      bNeg=bNeg, x1Neg=x1Neg, y1Neg=y1Neg, xyNeg=xyNeg,
                                      x2Neg=x2Neg, y2Neg=y2Neg, **kwargs)
 
        dipoleModel = DipoleModel()
 
        modelFunctor = dipoleModelFunctor  # dipoleModel.makeModel does not work for now.
        # Create the lmfit model (lmfit uses scipy 'leastsq' option by default - Levenberg-Marquardt)
        # We have to (later) filter out the nans by hand in our input to gmod.fit().
        # The only independent variable in the model is "x"; lmfit tries to
        # introspect variables and parameters from the function signature, but
        # gets it wrong for the model signature above.
        gmod = lmfit.Model(modelFunctor, independent_vars=["x"], verbose=verbose)
 
        # Add the constraints for centroids, fluxes.
        # starting constraint - near centroid of footprint
        fpCentroid = np.array([fp.getCentroid().getX(), fp.getCentroid().getY()])
        cenNeg = cenPos = fpCentroid
 
        pks = fp.getPeaks()
 
        if len(pks) >= 1:
            cenPos = pks[0].getF()    # if individual (merged) peaks were detected, use those
        if len(pks) >= 2:    # peaks are already sorted by centroid flux so take the most negative one
            cenNeg = pks[-1].getF()
 
        # For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
        # First assume dipole is not separated by more than 5*psfSigma.
        maxSep = self.psfSigma * maxSepInSigma
 
        # As an initial guess -- assume the dipole is close to the center of the footprint.
        if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
            cenPos = fpCentroid
        if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
            cenPos = fpCentroid
 
        # parameter hints/constraints: https://lmfit.github.io/lmfit-py/model.html#model-param-hints-section
        # might make sense to not use bounds -- see http://lmfit.github.io/lmfit-py/bounds.html
        # also see this discussion -- https://github.com/scipy/scipy/issues/3129
        gmod.set_param_hint('xcenPos', value=cenPos[0],
                            min=cenPos[0]-maxSep, max=cenPos[0]+maxSep)
        gmod.set_param_hint('ycenPos', value=cenPos[1],
                            min=cenPos[1]-maxSep, max=cenPos[1]+maxSep)
        gmod.set_param_hint('xcenNeg', value=cenNeg[0],
                            min=cenNeg[0]-maxSep, max=cenNeg[0]+maxSep)
        gmod.set_param_hint('ycenNeg', value=cenNeg[1],
                            min=cenNeg[1]-maxSep, max=cenNeg[1]+maxSep)
 
        # Use the (flux under the dipole)*5 for an estimate.
        # Lots of testing showed that having startingFlux be too high was better than too low.
        startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
        posFlux = negFlux = startingFlux
 
        # TBD: set max. flux limit?
        gmod.set_param_hint('flux', value=posFlux, min=0.1)
 
        if separateNegParams:
            # TBD: set max negative lobe flux limit?
            gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
 
        # Fixed parameters (don't fit for them if there are no pre-sub images or no gradient fit requested):
        # Right now (fitBackground == 1), we fit a linear model to the background and then subtract
        # it from the data and then don't fit the background again (this is faster).
        # A slower alternative (fitBackground == 2) is to use the estimated background parameters as
        # starting points in the integrated model fit. That is currently not performed by default,
        # but might be desirable in some cases.
        bgParsPos = bgParsNeg = (0., 0., 0.)
        if ((rel_weight > 0.) and (fitBackground != 0) and (bgGradientOrder >= 0)):
            pbg = 0.
            bgFitImage = self.posImage if self.posImage is not None else self.negImage
            # Fit the gradient to the background (linear model)
            bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
                                                                       order=bgGradientOrder)
 
            # Generate the gradient and subtract it from the pre-subtraction image data
            if fitBackground == 1:
                in_x = dipoleModel._generateXYGrid(bbox)
                pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsPos))
                z[1, :] -= pbg
                z[1, :] -= np.nanmedian(z[1, :])
                posFlux = np.nansum(z[1, :])
                gmod.set_param_hint('flux', value=posFlux*1.5, min=0.1)
 
                if separateNegParams and self.negImage is not None:
                    bgParsNeg = dipoleModel.fitFootprintBackground(source, self.negImage,
                                                                   order=bgGradientOrder)
                    pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsNeg))
                z[2, :] -= pbg
                z[2, :] -= np.nanmedian(z[2, :])
                if separateNegParams:
                    negFlux = np.nansum(z[2, :])
                    gmod.set_param_hint('fluxNeg', value=negFlux*1.5, min=0.1)
 
            # Do not subtract the background from the images but include the background parameters in the fit
            if fitBackground == 2:
                if bgGradientOrder >= 0:
                    gmod.set_param_hint('b', value=bgParsPos[0])
                    if separateNegParams:
                        gmod.set_param_hint('bNeg', value=bgParsNeg[0])
                if bgGradientOrder >= 1:
                    gmod.set_param_hint('x1', value=bgParsPos[1])
                    gmod.set_param_hint('y1', value=bgParsPos[2])
                    if separateNegParams:
                        gmod.set_param_hint('x1Neg', value=bgParsNeg[1])
                        gmod.set_param_hint('y1Neg', value=bgParsNeg[2])
                if bgGradientOrder >= 2:
                    gmod.set_param_hint('xy', value=bgParsPos[3])
                    gmod.set_param_hint('x2', value=bgParsPos[4])
                    gmod.set_param_hint('y2', value=bgParsPos[5])
                    if separateNegParams:
                        gmod.set_param_hint('xyNeg', value=bgParsNeg[3])
                        gmod.set_param_hint('x2Neg', value=bgParsNeg[4])
                        gmod.set_param_hint('y2Neg', value=bgParsNeg[5])
 
        y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
        in_x = np.array([x, y]).astype(np.float64)
        in_x[0, :] -= in_x[0, :].mean()  # center it!
        in_x[1, :] -= in_x[1, :].mean()
 
        # Instead of explicitly using a mask to ignore flagged pixels, just set the ignored pixels'
        #  weights to 0 in the fit. TBD: need to inspect mask planes to set this mask.
        mask = np.ones_like(z, dtype=bool)  # TBD: set mask values to False if the pixels are to be ignored
 
        # I'm not sure about the variance planes in the diffim (or convolved pre-sub. images
        # for that matter) so for now, let's just do an un-weighted least-squares fit
        # (override weights computed above).
        weights = mask.astype(np.float64)
        if self.posImage is not None and rel_weight > 0.:
            weights = np.array([np.ones_like(diArr), np.ones_like(diArr)*rel_weight,
                                np.ones_like(diArr)*rel_weight])
 
        # Set the weights to zero if mask is False
        if np.any(~mask):
            weights[~mask] = 0.
 
        # Filter out any nans, and make the weights 0.
        nans = (np.isnan(z) | np.isnan(weights))
        nNans = nans.sum()
        if nNans > 0:
            if nNans < len(z):
                z[nans] = np.nanmedian(z)
            else:
                z[nans] = 0
        weights[nans] = 0
 
        # Note that although we can, we're not required to set initial values for params here,
        # since we set their param_hint's above.
        # Can add "method" param to not use 'leastsq' (==levenberg-marquardt), e.g. "method='nelder'"
        with warnings.catch_warnings():
            # Ignore lmfit unknown argument warnings:
            # "psf, rel_weight, footprint, modelObj" all become pass-through kwargs for makeModel.
            warnings.filterwarnings("ignore", "The keyword argument .* does not match", UserWarning)
            result = gmod.fit(z, weights=weights, x=in_x, max_nfev=250,
                              method="leastsq",  # TODO: try using `least_squares` here for speed/robustness
                              verbose=verbose,
                              # see scipy docs for the meaning of these keywords
                              fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
                                       # Our model is float32 internally, so we need a larger epsfcn.
                                       'epsfcn': 1e-8},
                              psf=self.diffim.getPsf(),  # hereon: kwargs that get passed to makeModel()
                              rel_weight=rel_weight,
                              footprint=fp,
                              modelObj=dipoleModel)
 
        if verbose:  # the ci_report() seems to fail if neg params are constrained -- TBD why.
            # Never wanted in production - this takes a long time (longer than the fit!)
            # This is how to get confidence intervals out:
            #    https://lmfit.github.io/lmfit-py/confidence.html and
            #    http://cars9.uchicago.edu/software/python/lmfit/model.html
            print(result.fit_report(show_correl=False))
            if separateNegParams:
                print(result.ci_report())
 
        return result
 

Member Data Documentation

_private_version_

str lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm._private_version_ = '0.0.5'

staticprotected

Definition at line 499 of file dipoleFitTask.py.

debug

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.debug = lsstDebug.Info(__name__).debug

Definition at line 541 of file dipoleFitTask.py.

diffim

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.diffim = diffim

Definition at line 529 of file dipoleFitTask.py.

log

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.log = logging.getLogger(__name__)

Definition at line 538 of file dipoleFitTask.py.

negImage

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.negImage = negImage

Definition at line 531 of file dipoleFitTask.py.

posImage

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.posImage = posImage

Definition at line 530 of file dipoleFitTask.py.

psfSigma

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.psfSigma = None

Definition at line 532 of file dipoleFitTask.py.

---

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

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-12.0.0/Linux/ip_diffim/g5b326b94bb+db962c32ee/python/lsst/ip/diffim/dipoleFitTask.py