Merge pull request #362 from cyriltasse/issue-285

Issue 285

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()

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():
 

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)<np.abs(PreviousFlux):
-                        PreviousFlux=ThisFlux
+                if 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" %

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