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
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
588 z = z.copy()
589 weights = 1. / subim.variance.array
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:
597 posSubim = subim.clone()
598 posSubim += negSubim
599 if self.negImage is None:
600 negSubim = posSubim.clone()
601 negSubim -= subim
602
603 z = np.append([z], [posSubim.image.array,
604 negSubim.image.array], axis=0)
605
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.
610
611
612
613
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
632
633
634 gmod = lmfit.Model(modelFunctor, verbose=verbose)
635
636
637
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()
645 if len(pks) >= 2:
646 cenNeg = pks[-1].getF()
647
648
649
650 maxSep = self.psfSigma * maxSepInSigma
651
652
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
659
660
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
671
672 startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
673 posFlux = negFlux = startingFlux
674
675
676 gmod.set_param_hint('flux', value=posFlux, min=0.1)
677
678 if separateNegParams:
679
680 gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
681
682
683
684
685
686
687
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
693 bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
694 order=bgGradientOrder)
695
696
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
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()
739 in_x[1, :] -= in_x[1, :].mean()
740
741
742
743 mask = np.ones_like(z, dtype=bool)
744
745
746
747
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
754 if np.any(~mask):
755 weights[~mask] = 0.
756
757
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
768
769
770 with warnings.catch_warnings():
771
772
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",
776 verbose=verbose,
777
778 fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
779
780 'epsfcn': 1e-8},
781 psf=self.diffim.getPsf(),
782 rel_weight=rel_weight,
783 footprint=fp,
784 modelObj=dipoleModel)
785
786 if verbose:
787
788
789
790
791 print(result.fit_report(show_correl=False))
792 if separateNegParams:
793 print(result.ci_report())
794
795 return result
796