From d0a93a43e5470c7043d16897ae5ed16505e3ddf6 Mon Sep 17 00:00:00 2001
From: Tammo Jan Dijkema <dijkema@astron.nl>
Date: Tue, 3 Apr 2018 10:58:19 +0000
Subject: [PATCH] Task #11536: add stalling detection to DPPP_DDECal
---
.../include/DPPP_DDECal/MultiDirSolver.h | 10 +++-
CEP/DP3/DPPP_DDECal/src/DDECal.cc | 1 +
CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc | 50 +++++++++++++++----
3 files changed, 50 insertions(+), 11 deletions(-)
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h
index a18de7c186e..aa217c3fda5 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h
+++ b/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;
diff --git a/CEP/DP3/DPPP_DDECal/src/DDECal.cc b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
index 522049c290d..eae0aa3352b 100644
--- a/CEP/DP3/DPPP_DDECal/src/DDECal.cc
+++ b/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));
diff --git a/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc b/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc
index 367432b02bb..2e432c30fa3 100644
--- a/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc
+++ b/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,
--
GitLab