LSST Applications g0265f82a02+d6b5cd48b5,g02d81e74bb+a41d3748ce,g1470d8bcf6+6be6c9203b,g2079a07aa2+14824f138e,g212a7c68fe+a4f2ea4efa,g2305ad1205+72971fe858,g295015adf3+ab2c85acae,g2bbee38e9b+d6b5cd48b5,g337abbeb29+d6b5cd48b5,g3ddfee87b4+31b3a28dff,g487adcacf7+082e807817,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+b2918d57ae,g5a732f18d5+66d966b544,g64a986408d+a41d3748ce,g858d7b2824+a41d3748ce,g8a8a8dda67+a6fc98d2e7,g99cad8db69+7fe4acdf18,g9ddcbc5298+d4bad12328,ga1e77700b3+246acaaf9c,ga8c6da7877+84af8b3ff8,gb0e22166c9+3863383f4c,gb6a65358fc+d6b5cd48b5,gba4ed39666+9664299f35,gbb8dafda3b+d8d527deb2,gc07e1c2157+b2dbe6b631,gc120e1dc64+61440b2abb,gc28159a63d+d6b5cd48b5,gcf0d15dbbd+31b3a28dff,gdaeeff99f8+a38ce5ea23,ge6526c86ff+39927bb362,ge79ae78c31+d6b5cd48b5,gee10cc3b42+a6fc98d2e7,gf1cff7945b+a41d3748ce,v24.1.5.rc1
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | Static Protected Attributes | List of all members
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
 
 posImage
 
 negImage
 
 psfSigma
 
 log
 
 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 491 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 517 of file dipoleFitTask.py.

517 def __init__(self, diffim, posImage=None, negImage=None):
518 """Algorithm to run dipole measurement on a diaSource
519
520 Parameters
521 ----------
522 diffim : `lsst.afw.image.Exposure`
523 Exposure on which the diaSources were detected
524 posImage : `lsst.afw.image.Exposure`
525 "Positive" exposure from which the template was subtracted
526 negImage : `lsst.afw.image.Exposure`
527 "Negative" exposure which was subtracted from the posImage
528 """
529
530 self.diffim = diffim
531 self.posImage = posImage
532 self.negImage = negImage
533 self.psfSigma = None
534 if diffim is not None:
535 diffimPsf = diffim.getPsf()
536 diffimAvgPos = diffimPsf.getAveragePosition()
537 self.psfSigma = diffimPsf.computeShape(diffimAvgPos).getDeterminantRadius()
538
539 self.log = logging.getLogger(__name__)
540
541 import lsstDebug
542 self.debug = lsstDebug.Info(__name__).debug
543

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

904 def displayFitResults(self, footprint, result):
905 """Display data, model fits and residuals (currently uses matplotlib display functions).
906
907 Parameters
908 ----------
909 footprint : `lsst.afw.detection.Footprint`
910 Footprint containing the dipole that was fit
911 result : `lmfit.MinimizerResult`
912 `lmfit.MinimizerResult` object returned by `lmfit` optimizer
913 """
914 try:
915 import matplotlib.pyplot as plt
916 except ImportError as err:
917 self.log.warning('Unable to import matplotlib: %s', err)
918 raise err
919
920 def display2dArray(arr, title='Data', extent=None):
921 """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range.
922 """
923 fig = plt.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent)
924 plt.title(title)
925 plt.colorbar(fig, cmap='gray')
926 return fig
927
928 z = result.data
929 fit = result.best_fit
930 bbox = footprint.getBBox()
931 extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY())
932 if z.shape[0] == 3:
933 plt.figure(figsize=(8, 8))
934 for i in range(3):
935 plt.subplot(3, 3, i*3+1)
936 display2dArray(z[i, :], 'Data', extent=extent)
937 plt.subplot(3, 3, i*3+2)
938 display2dArray(fit[i, :], 'Model', extent=extent)
939 plt.subplot(3, 3, i*3+3)
940 display2dArray(z[i, :] - fit[i, :], 'Residual', extent=extent)
941 else:
942 plt.figure(figsize=(8, 2.5))
943 plt.subplot(1, 3, 1)
944 display2dArray(z, 'Data', extent=extent)
945 plt.subplot(1, 3, 2)
946 display2dArray(fit, 'Model', extent=extent)
947 plt.subplot(1, 3, 3)
948 display2dArray(z - fit, 'Residual', extent=extent)
949
950 plt.show()
951
952
953@measBase.register("ip_diffim_DipoleFit")

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

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

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

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

Member Data Documentation

◆ _private_version_

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

Definition at line 500 of file dipoleFitTask.py.

◆ debug

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.debug

Definition at line 542 of file dipoleFitTask.py.

◆ diffim

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.diffim

Definition at line 530 of file dipoleFitTask.py.

◆ log

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.log

Definition at line 539 of file dipoleFitTask.py.

◆ negImage

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.negImage

Definition at line 532 of file dipoleFitTask.py.

◆ posImage

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.posImage

Definition at line 531 of file dipoleFitTask.py.

◆ psfSigma

lsst.ip.diffim.dipoleFitTask.DipoleFitAlgorithm.psfSigma

Definition at line 533 of file dipoleFitTask.py.


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