diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
index 2bbb25cf5fcdbd8f436a2fb6e80394750e54055d..5d87ba0731ed81dc7f327760b7806b5489db4b3d 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
@@ -157,6 +157,8 @@ namespace LOFAR {
MultiDirSolver itsMultiDirSolver;
bool itsFullMatrixMinimalization;
bool itsApproximateTEC;
+ std::string itsStatFilename;
+ std::unique_ptr<std::ofstream> itsStatStream;
};
} //# end namespace
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h
index 517b78c276df1ca9db2b5789be6ab9172da9503d..a18de7c186ebbd7b66824c2022f1f45bcd33eb49 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/MultiDirSolver.h
@@ -55,8 +55,10 @@ public:
// data[i] is een pointer naar de data voor tijdstap i, vanaf die pointer staat het in volgorde als in MS (bl, chan, pol)
// mdata[i] is een pointer voor tijdstap i naar arrays van ndir model data pointers (elk van die data pointers staat in zelfde volgorde als data)
// solutions[ch] is een pointer voor channelblock ch naar antenna x directions oplossingen.
- SolveResult processScalar(std::vector<Complex*>& data, std::vector<std::vector<Complex* > >& modelData,
- std::vector<std::vector<DComplex> >& solutions, double time);
+ SolveResult processScalar(std::vector<Complex*>& data,
+ std::vector<std::vector<Complex* > >& modelData,
+ std::vector<std::vector<DComplex> >& solutions, double time,
+ std::ostream* statStream);
/**
* Same as @ref processScalar(), but solves full Jones matrices.
@@ -67,7 +69,8 @@ public:
*/
SolveResult processFullMatrix(std::vector<Complex *>& data,
std::vector<std::vector<Complex *> >& modelData,
- std::vector<std::vector<DComplex> >& solutions, double time);
+ std::vector<std::vector<DComplex> >& solutions, double time,
+ std::ostream* statStream);
void set_phase_only(bool phaseOnly) { _phaseOnly = phaseOnly; }
@@ -114,8 +117,12 @@ private:
* Assign the solutions in nextSolutions to the solutions.
* @returns whether the solutions have been converged.
*/
- bool assignSolutions(std::vector<std::vector<DComplex> >& solutions,
- std::vector<std::vector<DComplex> >& nextSolutions, bool useConstraintAccuracy) const;
+ bool assignSolutions(
+ std::vector<std::vector<DComplex> >& solutions,
+ std::vector<std::vector<DComplex> >& nextSolutions,
+ bool useConstraintAccuracy,
+ double& sum, double& normSum
+ ) const;
size_t _nAntennas, _nDirections, _nChannels, _nChannelBlocks;
std::vector<int> _ant1, _ant2;
diff --git a/CEP/DP3/DPPP_DDECal/src/DDECal.cc b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
index 9678ccab75844a00cc77b844be6264248ae4dca3..5213942dc7191cd11263b0d5d1ca77b92f7ba53d 100644
--- a/CEP/DP3/DPPP_DDECal/src/DDECal.cc
+++ b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
@@ -92,7 +92,8 @@ namespace LOFAR {
itsCoreConstraint(parset.getDouble (prefix + "coreconstraint", 0.0)),
itsScreenCoreConstraint(parset.getDouble (prefix + "tecscreen.coreconstraint", 0.0)),
itsFullMatrixMinimalization(false),
- itsApproximateTEC(false)
+ itsApproximateTEC(false),
+ itsStatFilename(parset.getString(prefix + "statfilename", ""))
{
stringstream ss;
ss << parset;
@@ -108,6 +109,9 @@ namespace LOFAR {
itsMultiDirSolver.set_constraint_accuracy(parset.getDouble(prefix + "approxtolerance", tolerance*10.0));
itsMultiDirSolver.set_step_size(parset.getDouble(prefix + "stepsize", 0.2));
+ if(!itsStatFilename.empty())
+ itsStatStream.reset(new std::ofstream(itsStatFilename));
+
// Default directions are all patches
if (strDirections.empty()) {
string sourceDBName = parset.getString(prefix+"sourcedb");
@@ -509,12 +513,12 @@ namespace LOFAR {
{
solveResult = itsMultiDirSolver.processFullMatrix(itsDataPtrs, itsModelDataPtrs,
itsSols[itsTimeStep/itsSolInt],
- itsAvgTime / itsSolInt);
+ itsAvgTime / itsSolInt, itsStatStream.get());
}
else {
solveResult = itsMultiDirSolver.processScalar(itsDataPtrs, itsModelDataPtrs,
itsSols[itsTimeStep/itsSolInt],
- itsAvgTime / itsSolInt);
+ itsAvgTime / itsSolInt, itsStatStream.get());
}
itsTimerSolve.stop();
diff --git a/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc b/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc
index c73b0c4002ce5b98e12993a4d341acc5cd2d1e55..a088e164a62ce5a4bd9790f9dcaf7a20459bda65 100644
--- a/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc
+++ b/CEP/DP3/DPPP_DDECal/src/MultiDirSolver.cc
@@ -86,9 +86,11 @@ 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) const
+ std::vector<std::vector<DComplex> >& nextSolutions, bool useConstraintAccuracy,
+ double& sum, double& normSum) const
{
- double normSum = 0.0, sum = 0.0;
+ sum = 0.0;
+ normSum = 0.0;
// Calculate the norm of the difference between the old and new solutions
size_t n = 0;
for(size_t chBlock=0; chBlock<_nChannelBlocks; ++chBlock)
@@ -103,8 +105,8 @@ bool MultiDirSolver::assignSolutions(std::vector<std::vector<DComplex> >& soluti
solutions[chBlock][i] = nextSolutions[chBlock][i];
}
}
- normSum /= n;
sum /= n;
+ normSum /= n;
if(useConstraintAccuracy)
return normSum*_stepSize/sum <= _constraintAccuracy;
else
@@ -113,7 +115,8 @@ bool MultiDirSolver::assignSolutions(std::vector<std::vector<DComplex> >& soluti
MultiDirSolver::SolveResult MultiDirSolver::processScalar(std::vector<Complex *>& data,
std::vector<std::vector<Complex *> >& modelData,
- std::vector<std::vector<DComplex> >& solutions, double time)
+ std::vector<std::vector<DComplex> >& solutions, double time,
+ std::ostream* statStream)
{
const size_t nTimes = data.size();
SolveResult result;
@@ -200,7 +203,12 @@ MultiDirSolver::SolveResult MultiDirSolver::processScalar(std::vector<Complex *>
if(!constraintsSatisfied)
constrainedIterations = iteration+1;
- hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied);
+ double sum, normSum;
+ hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied, sum, normSum);
+ if(statStream != nullptr)
+ {
+ (*statStream) << iteration << '\t' << normSum*_stepSize/sum << '\t' << normSum << '\n';
+ }
iteration++;
hasPreviouslyConverged = hasConverged || hasPreviouslyConverged;
@@ -304,7 +312,8 @@ void MultiDirSolver::performScalarIteration(size_t channelBlockIndex,
MultiDirSolver::SolveResult MultiDirSolver::processFullMatrix(std::vector<Complex *>& data,
std::vector<std::vector<Complex *> >& modelData,
- std::vector<std::vector<DComplex> >& solutions, double time)
+ std::vector<std::vector<DComplex> >& solutions, double time,
+ std::ostream* statStream)
{
// This algorithm is basically the same as the scalar algorithm,
// but visibility values are extended to 2x2 matrices and concatenated
@@ -412,7 +421,12 @@ MultiDirSolver::SolveResult MultiDirSolver::processFullMatrix(std::vector<Comple
if(!constraintsSatisfied)
constrainedIterations = iteration+1;
- hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied);
+ double sum, normSum;
+ hasConverged = assignSolutions(solutions, nextSolutions, !constraintsSatisfied, sum, normSum);
+ if(statStream != nullptr)
+ {
+ (*statStream) << iteration << '\t' << normSum*_stepSize/sum << '\t' << normSum << '\n';
+ }
iteration++;
hasPreviouslyConverged = hasConverged || hasPreviouslyConverged;