diff --git a/DDFacet/Imager/ClassDeconvMachine.py b/DDFacet/Imager/ClassDeconvMachine.py index aa8a48c0..5884118b 100644 --- a/DDFacet/Imager/ClassDeconvMachine.py +++ b/DDFacet/Imager/ClassDeconvMachine.py @@ -164,7 +164,7 @@ def Init(self): del(self.GD["Deconv"]["MaxMajorIter"]) # If we do the deconvolution construct a model according to what is in MinorCycleConfig MinorCycleConfig=dict(self.GD["Deconv"]) - MinorCycleConfig["NCPU"]=self.GD["Parallel"]["NCPU"] + MinorCycleConfig["NCPU"] = self.GD["Parallel"]["NCPU"] MinorCycleConfig["NBand"]=MinorCycleConfig["NFreqBands"]=self.VS.NFreqBands MinorCycleConfig["GD"] = self.GD MinorCycleConfig["ImagePolDescriptor"] = self.VS.StokesConverter.RequiredStokesProducts() diff --git a/DDFacet/Imager/ClassFacetMachine.py b/DDFacet/Imager/ClassFacetMachine.py index b1fd26b3..2610829c 100644 --- a/DDFacet/Imager/ClassFacetMachine.py +++ b/DDFacet/Imager/ClassFacetMachine.py @@ -1190,7 +1190,7 @@ def FacetsToIm_Channel(self, kind="Dirty",ChanSel=None): NFacets=len(self.DicoImager.keys()) pBAR.render(0, NFacets) - numexpr.set_num_threads(self.GD["Parallel"]["NCPU"]) + # numexpr.set_num_threads(self.GD["Parallel"]["NCPU"]) # done in DDF.py for iFacet in self.DicoImager.keys(): diff --git a/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py b/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py index f26f343d..40ba393d 100644 --- a/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py +++ b/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py @@ -88,7 +88,7 @@ def __init__(self, Gain=0.3, if self.ModelMachine.DicoSMStacked["Type"] not in ("MSMF", "HMP"): raise ValueError("ModelMachine Type should be HMP") self.PSFHasChanged=False - # reset overall iteration counter + self._previous_initial_peak = None self.maincache = MainCache # reset overall iteration counter self._niter = 0 @@ -551,7 +551,7 @@ def updateRMS(self): def setMask(self,Mask): self.CurrentNegMask=Mask - def Deconvolve(self, ch=0,UpdateRMS=True): + def Deconvolve(self, ch=0, UpdateRMS=True): """ Runs minor cycle over image channel 'ch'. initMinor is number of minor iteration (keeps continuous count through major iterations) @@ -560,7 +560,7 @@ def Deconvolve(self, ch=0,UpdateRMS=True): Returns tuple of: return_code,continue,updated where return_code is a status string; continue is True if another cycle should be executed; - update is True if model has been updated (note that1-sparse1.model01.fitst update=False implies continue=False) + update is True if model has been updated (note update=False implies continue=False) """ if self._niter >= self.MaxMinorIter: return "MaxIter", False, False @@ -589,16 +589,16 @@ def Deconvolve(self, ch=0,UpdateRMS=True): if self.CurrentNegMask is not None: - print>>log,"Using externally defined Mask (self.CurrentNegMask)" + print>>log," using externally defined Mask (self.CurrentNegMask)" CurrentNegMask=self.CurrentNegMask elif self.MaskMachine: - print>>log,"Using MaskMachine Mask" + print>>log," using MaskMachine Mask" CurrentNegMask=self.MaskMachine.CurrentNegMask elif self._MaskArray is not None: - print>>log,"Using externally defined Mask (self._MaskArray)" + print>>log," using externally defined Mask (self._MaskArray)" CurrentNegMask=self._MaskArray else: - print>>log,"Not using a mask" + print>>log," not using a mask" CurrentNegMask=None x,y,MaxDirty = NpParallel.A_whereMax(self._PeakSearchImage,NCPU=self.NCPU,DoAbs=DoAbs,Mask=CurrentNegMask) @@ -609,7 +609,16 @@ def Deconvolve(self, ch=0,UpdateRMS=True): ThisFlux = abs(ThisFlux) # in weighted or noisemap mode, look up the true max as well trueMaxDirty = MaxDirty if self._peakMode is "normal" else ThisFlux + # return condition indicating cleaning is to be continued + cont = True + + if self._previous_initial_peak is not None and abs(ThisFlux) > self.GD["HMP"]["MajorStallThreshold"]*self._previous_initial_peak: + print>>log,ModColor.Str("STALL! dirty image peak %10.6g Jy, was %10.6g at previous major cycle." + % (ThisFlux, self._previous_initial_peak), col="red") + print>>log,ModColor.Str("This will be the last major cycle") + cont = False + self._previous_initial_peak = abs(ThisFlux) #x,y,MaxDirty=NpParallel.A_whereMax(self._MeanDirty.copy(),NCPU=1,DoAbs=DoAbs,Mask=self._MaskArray.copy()) #A=self._MeanDirty.copy() #A.flat[:]=np.arange(A.size)[:] @@ -697,6 +706,8 @@ def Deconvolve(self, ch=0,UpdateRMS=True): # pBAR.render(0,"g=%3.3f"%self.GainMachine.GiveGain()) PreviousFlux=ThisFlux + divergence_factor = 1 + max(self.GD["HMP"]["AllowResidIncrease"],0) + def GivePercentDone(ThisMaxFlux): fracDone = 1.-(ThisMaxFlux-StopFlux)/(MaxDirty-StopFlux) return max(int(round(100*fracDone)), 100) @@ -725,14 +736,13 @@ def GivePercentDone(ThisMaxFlux): # stop T.timeit("max0") - if not self.GD["HMP"]["AllowResidIncrease"]: - if np.abs(ThisFlux)>1.1*np.abs(PreviousFlux): - print>>log, ModColor.Str( - " [iter=%i] peak of %.3g Jy higher than previous one of %.3g Jy " % - (i, ThisFlux, PreviousFlux), col="red") - return "Diverging", True, True - elif np.abs(ThisFlux) divergence_factor*np.abs(PreviousFlux): + print>>log, ModColor.Str( + " [iter=%i] peak of %.3g Jy diverging w.r.t. floor of %.3g Jy " % + (i, ThisFlux, PreviousFlux), col="red") + return "Diverging", False, True + if np.abs(ThisFlux) < np.abs(PreviousFlux): + PreviousFlux = ThisFlux ThisPNR=ThisFlux/rms @@ -749,7 +759,7 @@ def GivePercentDone(ThisMaxFlux): (i, ThisPNR, ThisFlux), col="green") - cont = ThisFlux > self.FluxThreshold + cont = cont and ThisFlux > self.FluxThreshold if not cont: print>>log, ModColor.Str( " [iter=%i] absolute flux threshold of %.3g Jy has been reached, PNR %.3g" % diff --git a/DDFacet/Parset/DefaultParset.cfg b/DDFacet/Parset/DefaultParset.cfg index 9e040139..f165ced6 100755 --- a/DDFacet/Parset/DefaultParset.cfg +++ b/DDFacet/Parset/DefaultParset.cfg @@ -234,7 +234,11 @@ Scales = [0] # List of scales to use. #metavar:LIST Ratios = [] # @cyriltasse please document NTheta = 6 # Number of PA steps to use. #metavar:N SolverMode = PI # Solver mode: pseudoinverse, or non-negative least squares. #options:PI|NNLS -AllowResidIncrease = True # Allow the maximum residual to increase within the minor cycle +AllowResidIncrease = 0.1 # Allow the maximum residual to increase by at most this much relative to + the lowest residual, before bailing out due to divergence.#metavar:FACTOR +MajorStallThreshold = 0.8 # Major cycle stall threshold. If the residual at the beginning of a major cycle + is above X*residual at the beginning of the previous major cycle, then we consider the deconvolution stalled + and bail out. #metavar:X #type:float Taper = 0 # Weighting taper size for HMP fit. If 0, determined automatically. #type:int Support = 0 # Basis function support size. If 0, determined automatically. #type:int PeakWeightImage = None # weigh the peak finding by given image