Task #11536: add stalling detection to DPPP_DDECal

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h

@@ -86,6 +86,8 @@ public:
     _constraintAccuracy = constraintAccuracy;
   }
   void set_step_size(double stepSize) { _stepSize = stepSize; }
+
+  void set_detect_stalling(bool detectStalling) { _detectStalling = detectStalling; }
   
   void add_constraint(Constraint* constraint) { _constraints.push_back(constraint); }
   
@@ -110,18 +112,21 @@ private:
 
   void makeStep(const std::vector<std::vector<DComplex> >& solutions,
     std::vector<std::vector<DComplex> >& nextSolutions) const;
+
+  bool detectStall(size_t iteration, const std::vector<double>& step_magnitudes) const;
                 
   void makeSolutionsFinite(std::vector<std::vector<DComplex> >& solutions, size_t perPol) const;
                 
   /**
    * Assign the solutions in nextSolutions to the solutions.
-   * @returns whether the solutions have been converged.
+   * @returns whether the solutions have converged. Appends the current step magnitude to step_magnitudes
    */
   bool assignSolutions(
     std::vector<std::vector<DComplex> >& solutions,
     std::vector<std::vector<DComplex> >& nextSolutions,
     bool useConstraintAccuracy,
-    double& sum, double& normSum
+    double& sum, double& normSum,
+    std::vector<double>& step_magnitudes
   ) const;
                              
   size_t _nAntennas, _nDirections, _nChannels, _nChannelBlocks;
@@ -131,6 +136,7 @@ private:
   size_t _maxIterations;
   double _accuracy, _constraintAccuracy;
   double _stepSize;
+  bool _detectStalling;
   bool _phaseOnly;
   std::vector<Constraint*> _constraints;
 

CEP/DP3/DPPP_DDECal/src/DDECal.cc

@@ -111,6 +111,7 @@ namespace LOFAR {
       itsMultiDirSolver.set_accuracy(tolerance);
       itsMultiDirSolver.set_constraint_accuracy(parset.getDouble(prefix + "approxtolerance", tolerance*10.0));
       itsMultiDirSolver.set_step_size(parset.getDouble(prefix + "stepsize", 0.2));
+      itsMultiDirSolver.set_detect_stalling(parset.getBool(prefix + "detectstalling", true));
 
       if(!itsStatFilename.empty())
 	itsStatStream.reset(new std::ofstream(itsStatFilename));

CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc

@@ -87,7 +87,7 @@ void MultiDirSolver::makeSolutionsFinite(std::vector<std::vector<DComplex> >& so
 
 bool MultiDirSolver::assignSolutions(std::vector<std::vector<DComplex> >& solutions,
   std::vector<std::vector<DComplex> >& nextSolutions, bool useConstraintAccuracy,
-  double& sum, double& normSum) const
+  double& sum, double& normSum, std::vector<double>& step_magnitudes) const
 {
   sum = 0.0;
   normSum = 0.0;
@@ -107,6 +107,9 @@ bool MultiDirSolver::assignSolutions(std::vector<std::vector<DComplex> >& soluti
   }
   sum /= n;
   normSum /= n;
+
+  step_magnitudes.push_back(sqrt(normSum/sum)*_stepSize);
+
   if(useConstraintAccuracy)
     return sqrt(normSum/sum)*_stepSize <= _constraintAccuracy;
   else
@@ -173,7 +176,12 @@ MultiDirSolver::SolveResult MultiDirSolver::processScalar(std::vector<Complex *>
   bool
     hasConverged = false,
     hasPreviouslyConverged = false,
-    constraintsSatisfied = false;
+    constraintsSatisfied = false,
+    hasStalled = false;
+
+  std::vector<double> step_magnitudes;
+  step_magnitudes.reserve(_maxIterations);
+
   do {
     makeSolutionsFinite(solutions, 1);
     
@@ -204,7 +212,7 @@ MultiDirSolver::SolveResult MultiDirSolver::processScalar(std::vector<Complex *>
       constrainedIterations = iteration+1;
     
     double sum, normSum;
-    hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied, sum, normSum);
+    hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied, sum, normSum, step_magnitudes);
     if(statStream != nullptr)
     {
       (*statStream) << iteration << '\t' << normSum*_stepSize/sum << '\t' << normSum << '\n';
@@ -212,8 +220,11 @@ MultiDirSolver::SolveResult MultiDirSolver::processScalar(std::vector<Complex *>
     iteration++;
     
     hasPreviouslyConverged = hasConverged || hasPreviouslyConverged;
-    
-  } while(iteration < _maxIterations && (!hasConverged || !constraintsSatisfied));
+
+    if (_detectStalling)
+      hasStalled = detectStall(iteration, step_magnitudes);
+
+  } while(iteration < _maxIterations && (!hasConverged || !constraintsSatisfied) && !hasStalled);
   
   if(hasConverged)
     result.iterations = iteration;
@@ -393,7 +404,14 @@ MultiDirSolver::SolveResult MultiDirSolver::processFullMatrix(std::vector<Comple
   /// Start iterating
   ///
   size_t iteration = 0, constrainedIterations = 0;
-  bool hasConverged = false, hasPreviouslyConverged = false, constraintsSatisfied = false;
+  bool hasConverged = false,
+       hasPreviouslyConverged = false,
+       constraintsSatisfied = false,
+       hasStalled = false;
+
+  std::vector<double> step_magnitudes;
+  step_magnitudes.reserve(_maxIterations);
+
   do {
     makeSolutionsFinite(solutions, 4);
     
@@ -419,7 +437,7 @@ MultiDirSolver::SolveResult MultiDirSolver::processFullMatrix(std::vector<Comple
       constrainedIterations = iteration+1;
     
     double sum, normSum;
-    hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied, sum, normSum);
+    hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied, sum, normSum, step_magnitudes);
     if(statStream != nullptr)
     {
       (*statStream) << iteration << '\t' << normSum*_stepSize/sum << '\t' << normSum << '\n';
@@ -427,8 +445,12 @@ MultiDirSolver::SolveResult MultiDirSolver::processFullMatrix(std::vector<Comple
     iteration++;
     
     hasPreviouslyConverged = hasConverged || hasPreviouslyConverged;
-  } while(iteration < _maxIterations && (!hasConverged || !constraintsSatisfied));
-  
+
+    if (_detectStalling)
+      hasStalled = detectStall(iteration, step_magnitudes);
+
+  } while(iteration < _maxIterations && (!hasConverged || !constraintsSatisfied) && !hasStalled);
+ 
   if(hasConverged)
     result.iterations = iteration;
   else
@@ -437,6 +459,16 @@ MultiDirSolver::SolveResult MultiDirSolver::processFullMatrix(std::vector<Comple
   return result;
 }
 
+bool MultiDirSolver::detectStall(size_t iteration, const std::vector<double>& step_magnitudes) const
+{
+  if (iteration<30) {
+    return false;
+  } else {
+    return std::abs(step_magnitudes[iteration-1]/step_magnitudes[iteration-2]-1) < 1.e-4 &&
+           std::abs(step_magnitudes[iteration-2]/step_magnitudes[iteration-3]-1) < 1.e-4;
+  }
+}
+
 void MultiDirSolver::performFullMatrixIteration(size_t channelBlockIndex,
                              std::vector<Matrix>& gTimesCs,
                              std::vector<Matrix>& vs,