544 separateNegParams=True, verbose=False):
545 """Fit a dipole model to an input difference image.
546
547 Actually, fits the subimage bounded by the input source's
548 footprint) and optionally constrain the fit using the
549 pre-subtraction images posImage and negImage.
550
551 Parameters
552 ----------
553 source : TODO: DM-17458
554 TODO: DM-17458
555 tol : float, optional
556 TODO: DM-17458
557 rel_weight : `float`, optional
558 TODO: DM-17458
559 fitBackground : `int`, optional
560 TODO: DM-17458
561 bgGradientOrder : `int`, optional
562 TODO: DM-17458
563 maxSepInSigma : `float`, optional
564 TODO: DM-17458
565 separateNegParams : `bool`, optional
566 TODO: DM-17458
567 verbose : `bool`, optional
568 TODO: DM-17458
569
570 Returns
571 -------
572 result : `lmfit.MinimizerResult`
573 return `lmfit.MinimizerResult` object containing the fit
574 parameters and other information.
575 """
576
577
578 import lmfit
579
580 fp = source.getFootprint()
581 bbox = fp.getBBox()
582 subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)
583
584 z = diArr = subim.image.array
585
586 z = z.copy()
587 weights = 1. / subim.variance.array
588
589 if rel_weight > 0. and ((self.posImage is not None) or (self.negImage is not None)):
590 if self.negImage is not None:
591 negSubim = afwImage.MaskedImageF(self.negImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
592 if self.posImage is not None:
593 posSubim = afwImage.MaskedImageF(self.posImage.getMaskedImage(), bbox, origin=afwImage.PARENT)
594 if self.posImage is None:
595 posSubim = subim.clone()
596 posSubim += negSubim
597 if self.negImage is None:
598 negSubim = posSubim.clone()
599 negSubim -= subim
600
601 z = np.append([z], [posSubim.image.array,
602 negSubim.image.array], axis=0)
603
604 weights = np.append([weights], [1. / posSubim.variance.array * rel_weight,
605 1. / negSubim.variance.array * rel_weight], axis=0)
606 else:
607 rel_weight = 0.
608
609
610
611
612 def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None,
613 b=None, x1=None, y1=None, xy=None, x2=None, y2=None,
614 bNeg=None, x1Neg=None, y1Neg=None, xyNeg=None, x2Neg=None, y2Neg=None,
615 **kwargs):
616 """Generate dipole model with given parameters.
617
618 It simply defers to `modelObj.makeModel()`, where `modelObj` comes
619 out of `kwargs['modelObj']`.
620 """
621 modelObj = kwargs.pop('modelObj')
622 return modelObj.makeModel(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=fluxNeg,
623 b=b, x1=x1, y1=y1, xy=xy, x2=x2, y2=y2,
624 bNeg=bNeg, x1Neg=x1Neg, y1Neg=y1Neg, xyNeg=xyNeg,
625 x2Neg=x2Neg, y2Neg=y2Neg, **kwargs)
626
627 dipoleModel = DipoleModel()
628
629 modelFunctor = dipoleModelFunctor
630
631
632
633
634
635 gmod = lmfit.Model(modelFunctor, independent_vars=["x"], verbose=verbose)
636
637
638
639 fpCentroid = np.array([fp.getCentroid().getX(), fp.getCentroid().getY()])
640 cenNeg = cenPos = fpCentroid
641
642 pks = fp.getPeaks()
643
644 if len(pks) >= 1:
645 cenPos = pks[0].getF()
646 if len(pks) >= 2:
647 cenNeg = pks[-1].getF()
648
649
650
651 maxSep = self.psfSigma * maxSepInSigma
652
653
654 if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
655 cenPos = fpCentroid
656 if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
657 cenPos = fpCentroid
658
659
660
661
662 gmod.set_param_hint('xcenPos', value=cenPos[0],
663 min=cenPos[0]-maxSep, max=cenPos[0]+maxSep)
664 gmod.set_param_hint('ycenPos', value=cenPos[1],
665 min=cenPos[1]-maxSep, max=cenPos[1]+maxSep)
666 gmod.set_param_hint('xcenNeg', value=cenNeg[0],
667 min=cenNeg[0]-maxSep, max=cenNeg[0]+maxSep)
668 gmod.set_param_hint('ycenNeg', value=cenNeg[1],
669 min=cenNeg[1]-maxSep, max=cenNeg[1]+maxSep)
670
671
672
673 startingFlux = np.nansum(np.abs(diArr) - np.nanmedian(np.abs(diArr))) * 5.
674 posFlux = negFlux = startingFlux
675
676
677 gmod.set_param_hint('flux', value=posFlux, min=0.1)
678
679 if separateNegParams:
680
681 gmod.set_param_hint('fluxNeg', value=np.abs(negFlux), min=0.1)
682
683
684
685
686
687
688
689 bgParsPos = bgParsNeg = (0., 0., 0.)
690 if ((rel_weight > 0.) and (fitBackground != 0) and (bgGradientOrder >= 0)):
691 pbg = 0.
692 bgFitImage = self.posImage if self.posImage is not None else self.negImage
693
694 bgParsPos = bgParsNeg = dipoleModel.fitFootprintBackground(source, bgFitImage,
695 order=bgGradientOrder)
696
697
698 if fitBackground == 1:
699 in_x = dipoleModel._generateXYGrid(bbox)
700 pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsPos))
701 z[1, :] -= pbg
702 z[1, :] -= np.nanmedian(z[1, :])
703 posFlux = np.nansum(z[1, :])
704 gmod.set_param_hint('flux', value=posFlux*1.5, min=0.1)
705
706 if separateNegParams and self.negImage is not None:
707 bgParsNeg = dipoleModel.fitFootprintBackground(source, self.negImage,
708 order=bgGradientOrder)
709 pbg = dipoleModel.makeBackgroundModel(in_x, tuple(bgParsNeg))
710 z[2, :] -= pbg
711 z[2, :] -= np.nanmedian(z[2, :])
712 if separateNegParams:
713 negFlux = np.nansum(z[2, :])
714 gmod.set_param_hint('fluxNeg', value=negFlux*1.5, min=0.1)
715
716
717 if fitBackground == 2:
718 if bgGradientOrder >= 0:
719 gmod.set_param_hint('b', value=bgParsPos[0])
720 if separateNegParams:
721 gmod.set_param_hint('bNeg', value=bgParsNeg[0])
722 if bgGradientOrder >= 1:
723 gmod.set_param_hint('x1', value=bgParsPos[1])
724 gmod.set_param_hint('y1', value=bgParsPos[2])
725 if separateNegParams:
726 gmod.set_param_hint('x1Neg', value=bgParsNeg[1])
727 gmod.set_param_hint('y1Neg', value=bgParsNeg[2])
728 if bgGradientOrder >= 2:
729 gmod.set_param_hint('xy', value=bgParsPos[3])
730 gmod.set_param_hint('x2', value=bgParsPos[4])
731 gmod.set_param_hint('y2', value=bgParsPos[5])
732 if separateNegParams:
733 gmod.set_param_hint('xyNeg', value=bgParsNeg[3])
734 gmod.set_param_hint('x2Neg', value=bgParsNeg[4])
735 gmod.set_param_hint('y2Neg', value=bgParsNeg[5])
736
737 y, x = np.mgrid[bbox.getBeginY():bbox.getEndY(), bbox.getBeginX():bbox.getEndX()]
738 in_x = np.array([x, y]).astype(np.float64)
739 in_x[0, :] -= in_x[0, :].mean()
740 in_x[1, :] -= in_x[1, :].mean()
741
742
743
744 mask = np.ones_like(z, dtype=bool)
745
746
747
748
749 weights = mask.astype(np.float64)
750 if self.posImage is not None and rel_weight > 0.:
751 weights = np.array([np.ones_like(diArr), np.ones_like(diArr)*rel_weight,
752 np.ones_like(diArr)*rel_weight])
753
754
755 if np.any(~mask):
756 weights[~mask] = 0.
757
758
759 nans = (np.isnan(z) | np.isnan(weights))
760 nNans = nans.sum()
761 if nNans > 0:
762 if nNans < len(z):
763 z[nans] = np.nanmedian(z)
764 else:
765 z[nans] = 0
766 weights[nans] = 0
767
768
769
770
771 with warnings.catch_warnings():
772
773
774 warnings.filterwarnings("ignore", "The keyword argument .* does not match", UserWarning)
775 result = gmod.fit(z, weights=weights, x=in_x, max_nfev=250,
776 method="leastsq",
777 verbose=verbose,
778
779 fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
780
781 'epsfcn': 1e-8},
782 psf=self.diffim.getPsf(),
783 rel_weight=rel_weight,
784 footprint=fp,
785 modelObj=dipoleModel)
786
787 if verbose:
788
789
790
791
792 print(result.fit_report(show_correl=False))
793 if separateNegParams:
794 print(result.ci_report())
795
796 return result
797