diff --git a/.gitattributes b/.gitattributes
index 330172a49cf8a441de3aab2ebff0777b476b7756..cb42d5500f666ff55664a04f7f732c72ea977eb8 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -300,7 +300,7 @@ CEP/DP3/DPPP/etc/DPPP.log_prop -text
CEP/DP3/DPPP/include/DPPP/Apply.h -text
CEP/DP3/DPPP/include/DPPP/ApplyBeam.h -text
CEP/DP3/DPPP/include/DPPP/ApplyBeam.tcc -text
-CEP/DP3/DPPP/include/DPPP/ApplyCal.h -text svneol=unset#text/plain
+CEP/DP3/DPPP/include/DPPP/ApplyCal.h -text
CEP/DP3/DPPP/include/DPPP/Baseline.h -text
CEP/DP3/DPPP/include/DPPP/Cursor.h -text
CEP/DP3/DPPP/include/DPPP/CursorUtilCasa.h -text
@@ -317,6 +317,7 @@ CEP/DP3/DPPP/include/DPPP/H5ParmPredict.h -text
CEP/DP3/DPPP/include/DPPP/ModelComponent.h -text
CEP/DP3/DPPP/include/DPPP/ModelComponentVisitor.h -text
CEP/DP3/DPPP/include/DPPP/MultiMSReader.h -text
+CEP/DP3/DPPP/include/DPPP/OneApplyCal.h -text svneol=unset#text/plain
CEP/DP3/DPPP/include/DPPP/Patch.h -text
CEP/DP3/DPPP/include/DPPP/PhaseFitter.h -text
CEP/DP3/DPPP/include/DPPP/PointSource.h -text
@@ -336,7 +337,7 @@ CEP/DP3/DPPP/share/HBAdefault -text
CEP/DP3/DPPP/share/LBAdefault -text
CEP/DP3/DPPP/src/Apply.cc -text
CEP/DP3/DPPP/src/ApplyBeam.cc -text
-CEP/DP3/DPPP/src/ApplyCal.cc -text svneol=unset#text/plain
+CEP/DP3/DPPP/src/ApplyCal.cc -text
CEP/DP3/DPPP/src/DemixInfo.cc -text
CEP/DP3/DPPP/src/DemixWorker.cc -text
CEP/DP3/DPPP/src/DemixerNew.cc -text
@@ -350,6 +351,7 @@ CEP/DP3/DPPP/src/H5ParmPredict.cc -text
CEP/DP3/DPPP/src/ModelComponent.cc -text
CEP/DP3/DPPP/src/ModelComponentVisitor.cc -text
CEP/DP3/DPPP/src/MultiMSReader.cc -text
+CEP/DP3/DPPP/src/OneApplyCal.cc -text svneol=unset#text/plain
CEP/DP3/DPPP/src/Patch.cc -text
CEP/DP3/DPPP/src/PhaseFitter.cc -text
CEP/DP3/DPPP/src/PointSource.cc -text
@@ -388,6 +390,8 @@ CEP/DP3/DPPP/test/tGainCal.tab.tgz -text
CEP/DP3/DPPP/test/tGainCal_ref -text
CEP/DP3/DPPP/test/tH5Parm.cc -text
CEP/DP3/DPPP/test/tH5Parm.sh -text
+CEP/DP3/DPPP/test/tMultiApplyCal.run -text
+CEP/DP3/DPPP/test/tMultiApplyCal.sh -text
CEP/DP3/DPPP/test/tNDPPP-generic.in_MS.tgz -text svneol=unset#application/x-gzip
CEP/DP3/DPPP/test/tNDPPP.in_MS.tgz -text svneol=unset#application/x-compressed-tar
CEP/DP3/DPPP/test/tPredict.run -text
diff --git a/CEP/DP3/DPPP/include/DPPP/ApplyCal.h b/CEP/DP3/DPPP/include/DPPP/ApplyCal.h
index e047a263b2378713310a719f268a75f6482c0f5f..cf458c442a8b585056348c7a58e617b5b733a981 100644
--- a/CEP/DP3/DPPP/include/DPPP/ApplyCal.h
+++ b/CEP/DP3/DPPP/include/DPPP/ApplyCal.h
@@ -1,4 +1,4 @@
-//# ApplyCal.h: DPPP step class to apply a calibration correction to the data
+//# ApplyCal.h: DPPP step class to ApplyCal visibilities from a source model
//# Copyright (C) 2013
//# ASTRON (Netherlands Institute for Radio Astronomy)
//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
@@ -17,44 +17,43 @@
//# You should have received a copy of the GNU General Public License along
//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
//#
-//# $Id: ApplyCal.h 21598 2012-07-16 08:07:34Z diepen $
+//# $Id:
//#
//# @author Tammo Jan Dijkema
-#ifndef DPPP_APPLYCAL_H
-#define DPPP_APPLYCAL_H
+#ifndef DPPP_ApplyCal_H
+#define DPPP_ApplyCal_H
// @file
-// @brief DPPP step class to apply a calibration correction to the data
+// @brief DPPP step class to apply multiple calibration solutions
#include <DPPP/DPInput.h>
#include <DPPP/DPBuffer.h>
-#include <DPPP/H5Parm.h>
-#include <ParmDB/ParmFacade.h>
-#include <ParmDB/ParmSet.h>
-#include <ParmDB/Parm.h>
-#include <casacore/casa/Arrays/Cube.h>
-#include <casacore/casa/Arrays/ArrayMath.h>
-#include <DPPP/FlagCounter.h>
+
+#include <DPPP/OneApplyCal.h>
+
+#include <utility>
namespace LOFAR {
+
+ class ParameterSet;
+
namespace DPPP {
// @ingroup NDPPP
- // This class is a DPStep class applying calibration parameters to the data.
+ // This class is a DPStep class to apply multiple ParmDB or H5Parm
+ // solutions to data.
+
class ApplyCal: public DPStep
{
public:
-
- enum CorrectType {GAIN, FULLJONES, TEC, CLOCK, ROTATIONANGLE, SCALARPHASE, PHASE,
- ROTATIONMEASURE, SCALARAMPLITUDE, AMPLITUDE};
-
// Construct the object.
// Parameters are obtained from the parset using the given prefix.
ApplyCal (DPInput*, const ParameterSet&, const string& prefix,
bool substep=false, std::string predictDirection="");
- ApplyCal();
+ // Empty constructor
+ ApplyCal ();
virtual ~ApplyCal();
@@ -66,19 +65,20 @@ namespace LOFAR {
// Finish the processing of this step and subsequent steps.
virtual void finish();
- // Update the general info.
- virtual void updateInfo (const DPInfo&);
+ // Set the next step. It squeezes in the actual OneApplyCal steps
+ // between this ApplyCal step and the next step.
+ virtual void setNextStep (DPStep::ShPtr nextStep);
- // Show the step parameters.
- virtual void show (std::ostream&) const;
+ // Show the step. When ApplyCal is a step in the main chain, this does
+ // nothing; the nextStep mechanism in DPRun will call show on the actual
+ // OneApplyCals.
+ virtual void show(std::ostream&) const;
- // Show the timings.
+ // Show the timings. When ApplyCal is a step in the main chain, this does
+ // nothing; the nextStep mechanism in DPRun will call show on the actual
+ // OneApplyCals.
virtual void showTimings (std::ostream&, double duration) const;
- bool invert() {
- return itsInvert;
- }
-
// Invert a 2x2 matrix in place
static void invert (casacore::DComplex* v, double sigmaMMSE=0);
@@ -117,50 +117,11 @@ namespace LOFAR {
float* weight);
private:
- // Read parameters from the associated parmdb and store them in itsParms
- void updateParms (const double bufStartTime);
-
- // If needed, show the flag counts.
- virtual void showCounts (std::ostream&) const;
-
- void initDataArrays();
-
- // Check the number of polarizations in the parmdb or h5parm
- uint nPol(const std::string& parmName);
-
- static std::string correctTypeToString(CorrectType);
- static CorrectType stringToCorrectType(const string&);
-
//# Data members.
- DPInput* itsInput;
- DPBuffer itsBuffer;
- string itsName;
- string itsParmDBName;
- bool itsUseH5Parm;
- boost::shared_ptr<BBS::ParmFacade> itsParmDB;
- H5Parm itsH5Parm;
- H5Parm::SolTab itsSolTab;
- CorrectType itsCorrectType;
- bool itsInvert;
- uint itsTimeSlotsPerParmUpdate;
- double itsSigmaMMSE;
- bool itsUpdateWeights;
-
- uint itsCount; // number of steps
-
- // Expressions to search for in itsParmDB
- vector<casacore::String> itsParmExprs;
-
- // parameters, numparms, antennas, time x frequency
- casacore::Cube<casacore::DComplex> itsParms;
- uint itsTimeStep; // time step within current chunk
- uint itsNCorr;
- double itsTimeInterval;
- double itsLastTime; // last time of current chunk
- FlagCounter itsFlagCounter;
- bool itsUseAP; //# use ampl/phase or real/imag
- hsize_t itsDirection;
- NSTimer itsTimer;
+ bool itsIsSubstep;
+ string itsName;
+
+ std::vector<OneApplyCal::ShPtr> itsApplyCals;
};
} //# end namespace
diff --git a/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt b/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
index 12ddf1a6b3e06286bcff33209af6b3987cd08b1a..ce149ffe8b6b2b4fa133a2a55ae2164c466d0cc1 100644
--- a/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
+++ b/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
@@ -15,7 +15,7 @@ set(inst_HEADERS
ModelComponentVisitor.h
DemixerNew.h DemixInfo.h DemixWorker.h
ApplyBeam.h ApplyBeam.tcc
- Predict.h
+ Predict.h OneApplyCal.h
GainCal.h StefCal.h PhaseFitter.h
StManParsetKeys.h H5Parm.h DummyStep.h H5ParmPredict.h
)
diff --git a/CEP/DP3/DPPP/include/DPPP/DPStep.h b/CEP/DP3/DPPP/include/DPPP/DPStep.h
index 5bb207d2bb0ee4d7eab620fa18abfe121a58b44e..981197b47960df5a71adec51b1d890491ca0991b 100644
--- a/CEP/DP3/DPPP/include/DPPP/DPStep.h
+++ b/CEP/DP3/DPPP/include/DPPP/DPStep.h
@@ -116,7 +116,7 @@ namespace LOFAR {
{ return itsPrevStep; }
// Set the next step.
- void setNextStep (DPStep::ShPtr nextStep)
+ virtual void setNextStep (DPStep::ShPtr nextStep)
{ itsNextStep = nextStep;
nextStep->setPrevStep(this);
}
@@ -129,11 +129,11 @@ namespace LOFAR {
DPInfo& info()
{ return itsInfo; }
- private:
// Update the general info (called by setInfo).
// The default implementation copies the info.
virtual void updateInfo (const DPInfo&);
+ private:
//# Data members.
DPStep::ShPtr itsNextStep;
DPStep* itsPrevStep; // Normal pointer for back links, prevent
diff --git a/CEP/DP3/DPPP/include/DPPP/DummyStep.h b/CEP/DP3/DPPP/include/DPPP/DummyStep.h
index a7bd04060868104cc12646345ac0fce60330e5e1..9611cc0b86ea7153a9b4f2b43f4f360e64f18124 100644
--- a/CEP/DP3/DPPP/include/DPPP/DummyStep.h
+++ b/CEP/DP3/DPPP/include/DPPP/DummyStep.h
@@ -29,16 +29,7 @@
#include <DPPP/DPInput.h>
#include <DPPP/DPBuffer.h>
-#include <DPPP/Patch.h>
-#include <DPPP/SourceDBUtil.h>
-#include <DPPP/ApplyBeam.h>
-#include <DPPP/ModelComponent.h>
-#include <StationResponse/Station.h>
-#include <StationResponse/Types.h>
-#include <casacore/casa/Arrays/Cube.h>
-#include <casacore/casa/Quanta/MVEpoch.h>
-#include <casacore/measures/Measures/MEpoch.h>
-#include <casacore/casa/Arrays/ArrayMath.h>
+
#include <utility>
namespace LOFAR {
@@ -48,10 +39,7 @@ namespace LOFAR {
namespace DPPP {
// @ingroup NDPPP
- // This class is a DPStep class to DummyStep visibilities with optionally beam
-
- typedef std::pair<size_t, size_t> Baseline;
- typedef std::pair<ModelComponent::ConstPtr, Patch::ConstPtr> Source;
+ // This class is an empty DPStep subclass to use as implementation template
class DummyStep: public DPStep
{
diff --git a/CEP/DP3/DPPP/include/DPPP/H5ParmPredict.h b/CEP/DP3/DPPP/include/DPPP/H5ParmPredict.h
index a82968e6098fb4009de7882ddb94aecf7678b5dc..ad1a464583613167d05c7aac11033c1b7213d17d 100644
--- a/CEP/DP3/DPPP/include/DPPP/H5ParmPredict.h
+++ b/CEP/DP3/DPPP/include/DPPP/H5ParmPredict.h
@@ -30,6 +30,9 @@
#include <DPPP/DPInput.h>
#include <DPPP/DPBuffer.h>
#include <DPPP/Predict.h>
+
+#include <DPPP/H5Parm.h>
+
#include <utility>
namespace LOFAR {
diff --git a/CEP/DP3/DPPP/include/DPPP/OneApplyCal.h b/CEP/DP3/DPPP/include/DPPP/OneApplyCal.h
new file mode 100644
index 0000000000000000000000000000000000000000..f7e87b7d0dd7af39181919e3574debc66e8d3d3a
--- /dev/null
+++ b/CEP/DP3/DPPP/include/DPPP/OneApplyCal.h
@@ -0,0 +1,135 @@
+//# OneApplyCal.h: DPPP step class to apply a calibration correction to the data
+//# Copyright (C) 2013
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id: OneApplyCal.h 21598 2012-07-16 08:07:34Z diepen $
+//#
+//# @author Tammo Jan Dijkema
+
+#ifndef DPPP_ONEAPPLYCAL_H
+#define DPPP_ONEAPPLYCAL_H
+
+// @file
+// @brief DPPP step class to apply a calibration correction to the data
+
+#include <DPPP/DPInput.h>
+#include <DPPP/DPBuffer.h>
+#include <DPPP/H5Parm.h>
+#include <ParmDB/ParmFacade.h>
+#include <ParmDB/ParmSet.h>
+#include <ParmDB/Parm.h>
+#include <casacore/casa/Arrays/Cube.h>
+#include <casacore/casa/Arrays/ArrayMath.h>
+#include <DPPP/FlagCounter.h>
+
+namespace LOFAR {
+ namespace DPPP {
+ // @ingroup NDPPP
+
+ // This class is a DPStep class applying calibration parameters to the data.
+ // It only applies one correction.
+
+ class OneApplyCal: public DPStep
+ {
+ public:
+ // Define the shared pointer for this type.
+ typedef shared_ptr<OneApplyCal> ShPtr;
+
+ enum CorrectType {GAIN, FULLJONES, TEC, CLOCK, ROTATIONANGLE, SCALARPHASE, PHASE,
+ ROTATIONMEASURE, SCALARAMPLITUDE, AMPLITUDE};
+
+ // Construct the object.
+ // Parameters are obtained from the parset using the given prefix.
+ OneApplyCal (DPInput*, const ParameterSet&, const std::string& prefix,
+ const std::string& defaultPrefix,
+ bool substep=false, std::string predictDirection="");
+
+ virtual ~OneApplyCal();
+
+ // Process the data.
+ // It keeps the data.
+ // When processed, it invokes the process function of the next step.
+ virtual bool process (const DPBuffer&);
+
+ // Finish the processing of this step and subsequent steps.
+ virtual void finish();
+
+ // Update the general info.
+ virtual void updateInfo (const DPInfo&);
+
+ // Show the step parameters.
+ virtual void show (std::ostream&) const;
+
+ // Show the timings.
+ virtual void showTimings (std::ostream&, double duration) const;
+
+ bool invert() {
+ return itsInvert;
+ }
+
+ private:
+ // Read parameters from the associated parmdb and store them in itsParms
+ void updateParms (const double bufStartTime);
+
+ // If needed, show the flag counts.
+ virtual void showCounts (std::ostream&) const;
+
+ void initDataArrays();
+
+ // Check the number of polarizations in the parmdb or h5parm
+ uint nPol(const std::string& parmName);
+
+ static std::string correctTypeToString(CorrectType);
+ static CorrectType stringToCorrectType(const string&);
+
+ //# Data members.
+ DPInput* itsInput;
+ DPBuffer itsBuffer;
+ string itsName;
+ string itsParmDBName;
+ bool itsUseH5Parm;
+ boost::shared_ptr<BBS::ParmFacade> itsParmDB;
+ H5Parm itsH5Parm;
+ H5Parm::SolTab itsSolTab;
+ CorrectType itsCorrectType;
+ bool itsInvert;
+ uint itsTimeSlotsPerParmUpdate;
+ double itsSigmaMMSE;
+ bool itsUpdateWeights;
+
+ uint itsCount; // number of steps
+
+ // Expressions to search for in itsParmDB
+ vector<casacore::String> itsParmExprs;
+
+ // parameters, numparms, antennas, time x frequency
+ casacore::Cube<casacore::DComplex> itsParms;
+ uint itsTimeStep; // time step within current chunk
+ uint itsNCorr;
+ double itsTimeInterval;
+ double itsLastTime; // last time of current chunk
+ FlagCounter itsFlagCounter;
+ bool itsUseAP; //# use ampl/phase or real/imag
+ hsize_t itsDirection;
+ NSTimer itsTimer;
+ };
+
+ } //# end namespace
+}
+
+#endif
diff --git a/CEP/DP3/DPPP/src/ApplyCal.cc b/CEP/DP3/DPPP/src/ApplyCal.cc
index 00b2f0fc5057463c07c82cd9b5a0975cae2e7d2e..2a41912f49a226abca37b94437f2864cbc717478 100644
--- a/CEP/DP3/DPPP/src/ApplyCal.cc
+++ b/CEP/DP3/DPPP/src/ApplyCal.cc
@@ -1,4 +1,4 @@
-//# ApplyCal.cc: DPPP step class to apply a calibration correction to the data
+//# GainCal.cc: DPPP step class to ApplyCal visibilities
//# Copyright (C) 2013
//# ASTRON (Netherlands Institute for Radio Astronomy)
//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
@@ -17,671 +17,122 @@
//# You should have received a copy of the GNU General Public License along
//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
//#
-//# $Id: ApplyCal.cc 21598 2012-07-16 08:07:34Z diepen $
+//# $Id: GainCal.cc 21598 2012-07-16 08:07:34Z diepen $
//#
//# @author Tammo Jan Dijkema
#include <lofar_config.h>
#include <DPPP/ApplyCal.h>
-#include <DPPP/DPBuffer.h>
-#include <DPPP/DPInfo.h>
-#include <DPPP/MSReader.h>
-#include <Common/ParameterSet.h>
-#include <Common/StringUtil.h>
-#include <Common/LofarLogger.h>
-#include <casacore/casa/Arrays/ArrayMath.h>
-#include <casacore/casa/OS/File.h>
+
#include <iostream>
-#include <iomanip>
+#include <Common/ParameterSet.h>
+#include <Common/ParameterValue.h>
+#include <Common/Timer.h>
-using namespace casacore;
-using namespace LOFAR::BBS;
+#include <stddef.h>
+#include <string>
+#include <sstream>
+#include <utility>
+#include <vector>
-/// Look at BBSKernel MeasurementExprLOFARUtil.cc and Apply.cc
+using namespace casacore;
namespace LOFAR {
namespace DPPP {
ApplyCal::ApplyCal (DPInput* input,
- const ParameterSet& parset,
- const string& prefix,
- bool substep,
- string predictDirection
- )
- : itsInput (input),
- itsName (prefix),
- itsParmDBName (parset.getString (prefix + "parmdb")),
- itsUseH5Parm (itsParmDBName.find(".h5")!=string::npos),
- itsTimeSlotsPerParmUpdate (parset.getInt (prefix +
- "timeslotsperparmupdate", 500)),
- itsSigmaMMSE (parset.getDouble (prefix + "MMSE.Sigma", 0)),
- itsUpdateWeights (parset.getBool (prefix + "updateweights", false)),
- itsCount (0),
- itsTimeStep (0),
- itsNCorr (0),
- itsTimeInterval (-1),
- itsLastTime (-1),
- itsUseAP (false)
+ const ParameterSet& parset,
+ const string& prefix,
+ bool substep,
+ string predictDirection
+ )
+ : itsIsSubstep(substep)
{
- ASSERT (!itsParmDBName.empty());
+ vector<string> subStepNames;
+ ParameterValue namesPar (parset.getString(prefix + "steps", ""));
- if (substep) {
- itsInvert=false;
+ if (namesPar.isVector()) {
+ subStepNames = namesPar.getStringVector();
} else {
- itsInvert = parset.getBool (prefix + "invert", true);
- }
-
- if (itsUseH5Parm) {
- string directionStr;
- if (substep) {
- directionStr = predictDirection;
- } else {
- directionStr = parset.getString (prefix + "direction", "");
- }
- itsH5Parm = H5Parm(itsParmDBName);
- string solTabName = toLower(parset.getString (prefix + "correction"));
- itsSolTab = itsH5Parm.getSolTab(solTabName);
- itsCorrectType = stringToCorrectType(itsSolTab.getType());
- if (itsCorrectType==PHASE && nPol("")==1) {
- itsCorrectType = SCALARPHASE;
- }
- if (itsCorrectType==AMPLITUDE && nPol("")==1) {
- itsCorrectType = SCALARAMPLITUDE;
- }
- if (directionStr=="") {
- ASSERT(!itsSolTab.hasAxis("dir") || itsSolTab.getAxis("dir").size==1);
+ subStepNames.push_back(namesPar.getString());
+ }
+
+ vector<string>::const_iterator subStepNameIter;
+ for (subStepNameIter = subStepNames.begin();
+ subStepNameIter != subStepNames.end();
+ ++subStepNameIter) {
+ string subStepName = (*subStepNameIter);
+ string subStepPrefix;
+ if (subStepName.empty()) {
+ // No substeps given, use parameters of this step
+ subStepPrefix = prefix;
} else {
- itsDirection = itsSolTab.getDirIndex(directionStr);
+ // Substeps given, use named parameters like applycal.applySol.parmdb
+ subStepPrefix = prefix + subStepName + ".";
}
- } else {
- string correctTypeStr = toLower(parset.getString (prefix + "correction", "gain"));
- itsCorrectType = stringToCorrectType(correctTypeStr);
+ itsApplyCals.push_back(OneApplyCal::ShPtr(new OneApplyCal(
+ input, parset, subStepPrefix, prefix, substep,
+ predictDirection)));
}
- if (itsCorrectType==FULLJONES && itsUpdateWeights) {
- ASSERTSTR (itsInvert, "Updating weights has not been implemented for invert=false and fulljones");
+ uint numSteps = itsApplyCals.size();
+ for (uint step=0; step<numSteps-1; ++step) {
+ itsApplyCals[step]->setNextStep(itsApplyCals[step+1]);
}
}
ApplyCal::ApplyCal()
{}
- string ApplyCal::correctTypeToString(CorrectType ct) {
- if (ct==GAIN) return "gain";
- if (ct==FULLJONES) return "fulljones";
- if (ct==TEC) return "tec";
- if (ct==CLOCK) return "clock";
- if (ct==SCALARPHASE) return "scalarphase";
- if (ct==SCALARAMPLITUDE) return "scalaramplitude";
- if (ct==ROTATIONANGLE) return "rotationangle";
- if (ct==ROTATIONMEASURE) return "rotationmeasure";
- if (ct==PHASE) return "phase";
- if (ct==AMPLITUDE) return "amplitude";
- THROW(Exception, "Unknown correction type: "<<ct);
- return "";
- }
-
- ApplyCal::CorrectType ApplyCal::stringToCorrectType(const string& ctStr) {
- if (ctStr=="gain") return GAIN;
- if (ctStr=="fulljones") return FULLJONES;
- if (ctStr=="tec") return TEC;
- if (ctStr=="clock") return CLOCK;
- if (ctStr=="scalarphase" || ctStr=="commonscalarphase") return SCALARPHASE;
- if (ctStr=="scalaramplitude" || ctStr=="commonscalaramplitude") return SCALARAMPLITUDE;
- if (ctStr=="phase") return PHASE;
- if (ctStr=="amplitude") return AMPLITUDE;
- if (ctStr=="rotationangle" || ctStr=="commonrotationangle") return ROTATIONANGLE;
- if (ctStr=="rotationmeasure") return ROTATIONMEASURE;
- THROW(Exception, "Unknown correction type: "<<ctStr);
- return GAIN;
- }
-
ApplyCal::~ApplyCal()
{}
- void ApplyCal::updateInfo (const DPInfo& infoIn)
+ void ApplyCal::setNextStep (DPStep::ShPtr nextStep)
{
- info() = infoIn;
- info().setNeedVisData();
- info().setWriteData();
- info().setWriteFlags();
- if (itsUpdateWeights) {
- info().setWriteWeights();
- }
- itsTimeInterval = infoIn.timeInterval();
- itsNCorr = infoIn.ncorr();
-
- ASSERT(itsNCorr==4);
-
- if (!itsUseH5Parm) { // Use ParmDB
- itsParmDB.reset(new BBS::ParmFacade(itsParmDBName));
- }
-
- // Detect if full jones solutions are present
- if (!itsUseH5Parm &&
- (itsCorrectType == GAIN || itsCorrectType==FULLJONES) &&
- (itsParmDB->getNames("Gain:0:1:*").size() +
- itsParmDB->getDefNames("Gain:0:1:*").size() >0 )) {
- itsCorrectType=FULLJONES;
- }
-
- // Detect if solutions are saved as Real/Imag or Ampl/Phase
- if (itsCorrectType == GAIN || itsCorrectType == FULLJONES ){
- if (itsUseH5Parm) {
- // H5Parm uses amplitude / phase by definition
- itsUseAP = true;
- } else {
- // Determine from values present in parmdb what to use
- if (!itsParmDB->getNames("Gain:0:0:Real*").empty()) {
- // Values with :Real present
- itsUseAP = false;
- } else if (!itsParmDB->getNames("Gain:0:0:Ampl*").empty() ||
- !itsParmDB->getNames("Phase:0:0:Ampl*").empty()) {
- // Values with :Ampl present
- itsUseAP = true;
- } else if (!itsParmDB->getDefNames("Gain:0:0:Real*").empty()) {
- // Defvalues with :Real present
- itsUseAP = false;
- } else if (!itsParmDB->getDefNames("Gain:0:0:Ampl*").empty() ||
- !itsParmDB->getDefNames("Gain:0:0:Phase*").empty()) {
- // Defvalues with :Ampl present
- itsUseAP = true;
- } else {
- THROW (Exception, "No gains found in parmdb "+itsParmDBName);
- }
- }
- }
-
- if (itsCorrectType == GAIN) {
- if (itsUseAP) {
- itsParmExprs.push_back("Gain:0:0:Ampl");
- itsParmExprs.push_back("Gain:0:0:Phase");
- itsParmExprs.push_back("Gain:1:1:Ampl");
- itsParmExprs.push_back("Gain:1:1:Phase");
- } else {
- itsParmExprs.push_back("Gain:0:0:Real");
- itsParmExprs.push_back("Gain:0:0:Imag");
- itsParmExprs.push_back("Gain:1:1:Real");
- itsParmExprs.push_back("Gain:1:1:Imag");
- }
- } else if (itsCorrectType == FULLJONES) {
- if (itsUseAP) {
- itsParmExprs.push_back("Gain:0:0:Ampl");
- itsParmExprs.push_back("Gain:0:0:Phase");
- itsParmExprs.push_back("Gain:0:1:Ampl");
- itsParmExprs.push_back("Gain:0:1:Phase");
- itsParmExprs.push_back("Gain:1:0:Ampl");
- itsParmExprs.push_back("Gain:1:0:Phase");
- itsParmExprs.push_back("Gain:1:1:Ampl");
- itsParmExprs.push_back("Gain:1:1:Phase");
- } else {
- itsParmExprs.push_back("Gain:0:0:Real");
- itsParmExprs.push_back("Gain:0:0:Imag");
- itsParmExprs.push_back("Gain:0:1:Real");
- itsParmExprs.push_back("Gain:0:1:Imag");
- itsParmExprs.push_back("Gain:1:0:Real");
- itsParmExprs.push_back("Gain:1:0:Imag");
- itsParmExprs.push_back("Gain:1:1:Real");
- itsParmExprs.push_back("Gain:1:1:Imag");
- }
- } else if (itsCorrectType == TEC) {
- if (nPol("TEC")==1) {
- itsParmExprs.push_back("TEC");
- }
- else {
- itsParmExprs.push_back("TEC:0");
- itsParmExprs.push_back("TEC:1");
- }
- } else if (itsCorrectType == CLOCK) {
- if (nPol("Clock")==1) {
- itsParmExprs.push_back("Clock");
- }
- else {
- itsParmExprs.push_back("Clock:0");
- itsParmExprs.push_back("Clock:1");
- }
- } else if (itsCorrectType == ROTATIONANGLE) {
- itsParmExprs.push_back("{Common,}RotationAngle");
- } else if (itsCorrectType == SCALARPHASE) {
- itsParmExprs.push_back("{Common,}ScalarPhase");
- } else if (itsCorrectType == ROTATIONMEASURE) {
- itsParmExprs.push_back("RotationMeasure");
- } else if (itsCorrectType == SCALARAMPLITUDE) {
- itsParmExprs.push_back("{Common,}ScalarAmplitude");
- } else if (itsCorrectType == PHASE) {
- ASSERT(itsUseH5Parm);
- itsParmExprs.push_back("Phase:0");
- itsParmExprs.push_back("Phase:1");
- } else if (itsCorrectType == AMPLITUDE) {
- ASSERT(itsUseH5Parm);
- itsParmExprs.push_back("Amplitude:0");
- itsParmExprs.push_back("Amplitude:1");
- } else {
- THROW (Exception, "Correction type "<<
- correctTypeToString(itsCorrectType)<<" is unknown");
- }
-
- initDataArrays();
- itsFlagCounter.init(getInfo());
-
- // Check that channels are evenly spaced
- if (info().nchan()>1) {
- Vector<Double> upFreq = info().chanFreqs()(
- Slicer(IPosition(1,1),
- IPosition(1,info().nchan()-1)));
- Vector<Double> lowFreq = info().chanFreqs()(
- Slicer(IPosition(1,0),
- IPosition(1,info().nchan()-1)));
- Double freqstep0=upFreq(0)-lowFreq(0);
- // Compare up to 1kHz accuracy
- bool regularChannels=allNearAbs(upFreq-lowFreq, freqstep0, 1.e3) &&
- allNearAbs(info().chanWidths(),
- info().chanWidths()(0), 1.e3);
- ASSERTSTR(regularChannels,
- "ApplyCal requires evenly spaced channels.");
- }
+ DPStep::setNextStep(itsApplyCals[0]);
+ itsApplyCals[itsApplyCals.size()-1]->setNextStep(nextStep);
}
- void ApplyCal::show (std::ostream& os) const
+ void ApplyCal::show(std::ostream& os) const
{
- os << "ApplyCal " << itsName << std::endl;
- if (itsUseH5Parm) {
- os << " H5Parm: " << itsParmDBName <<endl;
- } else {
- os << " parmdb: " << itsParmDBName << endl;
- }
- os << " correction: " << correctTypeToString(itsCorrectType) << endl;
- if (itsCorrectType==GAIN || itsCorrectType==FULLJONES) {
- os << " Ampl/Phase: " << boolalpha << itsUseAP << endl;
+ // If not a substep, show will be called by DPRun,
+ // through the nextStep() mechanism
+ if (itsIsSubstep) {
+ vector<OneApplyCal::ShPtr>::const_iterator applycalIter;
+
+ for (applycalIter = itsApplyCals.begin();
+ applycalIter != itsApplyCals.end();
+ applycalIter++) {
+ (*applycalIter)->show(os);
+ }
}
- os << " update weights: " << boolalpha << itsUpdateWeights << endl;
- os << " sigmaMMSE: " << itsSigmaMMSE << endl;
- os << " invert: " << boolalpha << itsInvert <<endl;
- os << " timeSlotsPerParmUpdate: " << itsTimeSlotsPerParmUpdate <<endl;
}
void ApplyCal::showTimings (std::ostream& os, double duration) const
{
- os << " ";
- FlagCounter::showPerc1 (os, itsTimer.getElapsed(), duration);
- os << " ApplyCal " << itsName << endl;
+ if (itsIsSubstep) {
+ vector<OneApplyCal::ShPtr>::const_iterator iter;
+ for (iter = itsApplyCals.begin();
+ iter != itsApplyCals.end();
+ iter++) {
+ (*iter)->showTimings(os, duration);
+ }
+ }
}
bool ApplyCal::process (const DPBuffer& bufin)
{
- itsTimer.start();
- itsBuffer.copy (bufin);
-
- if (bufin.getTime() > itsLastTime) {
- updateParms(bufin.getTime());
- itsTimeStep=0;
- }
- else {
- itsTimeStep++;
- }
-
- // Loop through all baselines in the buffer.
- size_t nbl = itsBuffer.getData().shape()[2];
-
- Complex* data = itsBuffer.getData().data();
-
- itsInput->fetchWeights (bufin, itsBuffer, itsTimer);
- float* weight = itsBuffer.getWeights().data();
-
- bool* flag = itsBuffer.getFlags().data();
-
- size_t nchan = itsBuffer.getData().shape()[1];
-
-#pragma omp parallel for
- for (size_t bl=0; bl<nbl; ++bl) {
- for (size_t chan=0;chan<nchan;chan++) {
- uint timeFreqOffset=(itsTimeStep*info().nchan())+chan;
- uint antA = info().getAnt1()[bl];
- uint antB = info().getAnt2()[bl];
- if (itsParms.shape()[0]>2) {
- applyFull( &itsParms(0, antA, timeFreqOffset),
- &itsParms(0, antB, timeFreqOffset),
- &data[bl * itsNCorr * nchan + chan * itsNCorr ],
- &weight[bl * itsNCorr * nchan + chan * itsNCorr ],
- &flag[ bl * itsNCorr * nchan + chan * itsNCorr ],
- bl, chan, itsUpdateWeights, itsFlagCounter);
- }
- else {
- applyDiag( &itsParms(0, antA, timeFreqOffset),
- &itsParms(0, antB, timeFreqOffset),
- &data[bl * itsNCorr * nchan + chan * itsNCorr ],
- &weight[bl * itsNCorr * nchan + chan * itsNCorr ],
- &flag[ bl * itsNCorr * nchan + chan * itsNCorr ],
- bl, chan, itsUpdateWeights, itsFlagCounter);
- }
- }
- }
-
- itsTimer.stop();
- getNextStep()->process(itsBuffer);
-
- itsCount++;
+ getNextStep()->process(bufin);
return true;
}
+
void ApplyCal::finish()
{
// Let the next steps finish.
getNextStep()->finish();
}
-
- void ApplyCal::updateParms (const double bufStartTime)
- {
- uint numAnts = info().antennaNames().size();
-
- // itsParms contains the parameters to a grid, first for all parameters
- // (e.g. Gain:0:0 and Gain:1:1), next all antennas, next over freq * time
- // as returned by ParmDB
- vector<vector<vector<double> > > parmvalues;
- parmvalues.resize(itsParmExprs.size());
- for (size_t i=0;i<parmvalues.size();++i) {
- parmvalues[i].resize(numAnts);
- }
-
- uint numFreqs (info().chanFreqs().size());
- double freqInterval (info().chanWidths()[0]);
- if (numFreqs>1) { // Handle data with evenly spaced gaps between channels
- freqInterval = info().chanFreqs()[1]-info().chanFreqs()[0];
- }
- double minFreq (info().chanFreqs()[0]-0.5*freqInterval);
- double maxFreq (info().chanFreqs()[numFreqs-1]+0.5*freqInterval);
- itsLastTime = bufStartTime - 0.5*itsTimeInterval +
- itsTimeSlotsPerParmUpdate * itsTimeInterval;
- uint numTimes = itsTimeSlotsPerParmUpdate;
-
- double lastMSTime = info().startTime() + info().ntime() * itsTimeInterval;
- if (itsLastTime > lastMSTime && !nearAbs(itsLastTime, lastMSTime, 1.e-3)) {
- itsLastTime = lastMSTime;
- numTimes = info().ntime() % itsTimeSlotsPerParmUpdate;
- }
-
- map<string, vector<double> > parmMap;
- map<string, vector<double> >::iterator parmIt;
-
- uint tfDomainSize=numTimes*numFreqs;
-
- // Fill parmvalues here, get raw data from H5Parm or ParmDB
- if (itsUseH5Parm) {
-#pragma omp critical(updateH5ParmValues)
-{
- // TODO: understand polarization etc.
- ASSERT(itsParmExprs.size()==1 || itsParmExprs.size()==2);
- hsize_t startTime = itsSolTab.getTimeIndex(bufStartTime);
- hsize_t startFreq = itsSolTab.getFreqIndex(info().chanFreqs()[0]);
- for (uint ant = 0; ant < numAnts; ++ant) {
-
- uint freqUpsampleFactor = numFreqs;
- if (itsSolTab.hasAxis("freq") && itsSolTab.getAxis("freq").size > 1) {
- double h5freqinterval = itsSolTab.getFreqInterval();
- freqUpsampleFactor = h5freqinterval/info().chanWidths()[0] + 0.5; // Round;
- ASSERT(near(h5freqinterval, freqUpsampleFactor*info().chanWidths()[0],1.e-5));
- }
-
- uint timeUpsampleFactor = numTimes;
- if (itsSolTab.getAxis("time").size > 1) {
- double h5timeInterval = itsSolTab.getTimeInterval();
- timeUpsampleFactor = h5timeInterval/itsTimeInterval+0.5; // Round
- ASSERT(near(h5timeInterval,timeUpsampleFactor*itsTimeInterval,1.e-5));
- }
-
- // Figure out whether time or frequency is first axis
- bool freqvariesfastest = true;
- if (itsSolTab.hasAxis("freq") &&
- itsSolTab.getAxisIndex("freq") < itsSolTab.getAxisIndex("time")) {
- freqvariesfastest = false;
- }
- ASSERT(freqvariesfastest);
-
- // Take the ceiling of numTimes/timeUpsampleFactor, same for freq
- uint numTimesInH5Parm = (numTimes+timeUpsampleFactor-1)/timeUpsampleFactor;
- uint numFreqsInH5Parm = (numFreqs+freqUpsampleFactor-1)/freqUpsampleFactor;
-
- for (uint pol=0; pol<itsParmExprs.size(); ++pol) {
- vector<double> rawsols;
- rawsols = itsSolTab.getValues(info().antennaNames()[ant],
- startTime, numTimesInH5Parm, 1,
- startFreq, numFreqsInH5Parm, 1, pol, itsDirection);
-
- parmvalues[pol][ant].resize(tfDomainSize);
-
- size_t tf=0;
- for (uint t=0; t<numTimesInH5Parm; ++t) {
- for (uint ti=0; ti<timeUpsampleFactor; ++ti) {
- for (uint f=0; f<numFreqsInH5Parm; ++f) {
- for (uint fi=0; fi<freqUpsampleFactor; ++fi) {
- if (tf<tfDomainSize) {
- parmvalues[pol][ant][tf++] = rawsols[t*numFreqsInH5Parm + f];
- }
- }
- }
- }
- }
- ASSERT(tf==tfDomainSize);
- }
- }
-} // End pragma omp critical
- } else { // Use ParmDB
- for (uint parmExprNum = 0; parmExprNum<itsParmExprs.size();++parmExprNum) {
- // parmMap contains parameter values for all antennas
- parmMap = itsParmDB->getValuesMap( itsParmExprs[parmExprNum] + "*",
- minFreq, maxFreq, freqInterval,
- bufStartTime-0.5*itsTimeInterval, itsLastTime,
- itsTimeInterval, true);
-
- // Resolve {Common,}Bla to CommonBla or Bla
- if (!parmMap.empty() &&
- itsParmExprs[parmExprNum].find("{") != std::string::npos) {
- uint colonPos = (parmMap.begin()->first).find(":");
- itsParmExprs[parmExprNum] = (parmMap.begin()->first).substr(0, colonPos);
- }
-
- for (uint ant = 0; ant < numAnts; ++ant) {
- parmIt = parmMap.find(
- itsParmExprs[parmExprNum] + ":" + info().antennaNames()[ant]);
-
- if (parmIt != parmMap.end()) {
- parmvalues[parmExprNum][ant].swap(parmIt->second);
- ASSERT(parmvalues[parmExprNum][ant].size()==tfDomainSize);
- } else {// No value found, try default
- Array<double> defValues;
- double defValue;
-
- if (itsParmDB->getDefValues(itsParmExprs[parmExprNum] + ":" +
- info().antennaNames()[ant]).size()==1) { // Default for antenna
- itsParmDB->getDefValues(itsParmExprs[parmExprNum] + ":" +
- info().antennaNames()[ant]).get(0,defValues);
- ASSERT(defValues.size()==1);
- defValue=defValues.data()[0];
- }
- else if (itsParmDB->getDefValues(itsParmExprs[parmExprNum]).size()
- == 1) { //Default value
- itsParmDB->getDefValues(itsParmExprs[parmExprNum]).get(0,defValues);
- ASSERT(defValues.size()==1);
- defValue=defValues.data()[0];
- } else if (itsParmExprs[parmExprNum].substr(0,5)=="Gain:") {
- defValue=0.;
- }
- else {
- THROW (Exception, "No parameter value found for "+
- itsParmExprs[parmExprNum]+":"+info().antennaNames()[ant]);
- }
-
- parmvalues[parmExprNum][ant].resize(tfDomainSize);
- for (uint tf=0; tf<tfDomainSize;++tf) {
- parmvalues[parmExprNum][ant][tf]=defValue;
- }
- }
- }
- }
- }
-
- ASSERT(parmvalues[0][0].size() <= tfDomainSize); // Catches multiple matches
-
- double freq;
-
- // Make parameters complex
- for (uint tf=0;tf<tfDomainSize;++tf) {
- for (uint ant=0;ant<numAnts;++ant) {
-
- freq=info().chanFreqs()[tf % numFreqs];
-
- if (itsCorrectType==GAIN) {
- if (itsUseAP) { // Data as Amplitude / Phase
- itsParms(0, ant, tf) = polar(parmvalues[0][ant][tf],
- parmvalues[1][ant][tf]);
- itsParms(1, ant, tf) = polar(parmvalues[2][ant][tf],
- parmvalues[3][ant][tf]);
- } else { // Data as Real / Imaginary
- itsParms(0, ant, tf) = DComplex(parmvalues[0][ant][tf],
- parmvalues[1][ant][tf]);
- itsParms(1, ant, tf) = DComplex(parmvalues[2][ant][tf],
- parmvalues[3][ant][tf]);
- }
- }
- else if (itsCorrectType==FULLJONES) {
- if (itsUseAP) { // Data as Amplitude / Phase
- itsParms(0, ant, tf) = polar(parmvalues[0][ant][tf],
- parmvalues[1][ant][tf]);
- itsParms(1, ant, tf) = polar(parmvalues[2][ant][tf],
- parmvalues[3][ant][tf]);
- itsParms(2, ant, tf) = polar(parmvalues[4][ant][tf],
- parmvalues[5][ant][tf]);
- itsParms(3, ant, tf) = polar(parmvalues[6][ant][tf],
- parmvalues[7][ant][tf]);
- } else { // Data as Real / Imaginary
- itsParms(0, ant, tf) = DComplex(parmvalues[0][ant][tf],
- parmvalues[1][ant][tf]);
- itsParms(1, ant, tf) = DComplex(parmvalues[2][ant][tf],
- parmvalues[3][ant][tf]);
- itsParms(2, ant, tf) = DComplex(parmvalues[4][ant][tf],
- parmvalues[5][ant][tf]);
- itsParms(3, ant, tf) = DComplex(parmvalues[6][ant][tf],
- parmvalues[7][ant][tf]);
- }
- }
- else if (itsCorrectType==TEC) {
- itsParms(0, ant, tf)=polar(1.,
- parmvalues[0][ant][tf] * -8.44797245e9 / freq);
- if (itsParmExprs.size() == 1) { // No TEC:0, only TEC:
- itsParms(1, ant, tf)=polar(1.,
- parmvalues[0][ant][tf] * -8.44797245e9 / freq);
- }
- else { // TEC:0 and TEC:1
- itsParms(1, ant, tf)=polar(1.,
- parmvalues[1][ant][tf] * -8.44797245e9 / freq);
- }
- }
- else if (itsCorrectType==CLOCK) {
- itsParms(0, ant, tf)=polar(1.,
- parmvalues[0][ant][tf] * freq * casacore::C::_2pi);
- if (itsParmExprs.size() == 1) { // No Clock:0, only Clock:
- itsParms(1, ant, tf)=polar(1.,
- parmvalues[0][ant][tf] * freq * casacore::C::_2pi);
- }
- else { // Clock:0 and Clock:1
- itsParms(1, ant, tf)=polar(1.,
- parmvalues[1][ant][tf] * freq * casacore::C::_2pi);
- }
- }
- else if (itsCorrectType==ROTATIONANGLE) {
- double phi=parmvalues[0][ant][tf];
- if (itsInvert) {
- phi = -phi;
- }
- double sinv=sin(phi);
- double cosv=cos(phi);
- itsParms(0, ant, tf) = cosv;
- itsParms(1, ant, tf) = -sinv;
- itsParms(2, ant, tf) = sinv;
- itsParms(3, ant, tf) = cosv;
- }
- else if (itsCorrectType==ROTATIONMEASURE) {
- double lambda2 = casacore::C::c / freq;
- lambda2 *= lambda2;
- double chi = parmvalues[0][ant][tf] * lambda2;
- if (itsInvert) {
- chi = -chi;
- }
- double sinv = sin(chi);
- double cosv = cos(chi);
- itsParms(0, ant, tf) = cosv;
- itsParms(1, ant, tf) = -sinv;
- itsParms(2, ant, tf) = sinv;
- itsParms(3, ant, tf) = cosv;
- }
- else if (itsCorrectType==PHASE || itsCorrectType==SCALARPHASE) {
- itsParms(0, ant, tf) = polar(1., parmvalues[0][ant][tf]);
- if (itsCorrectType==SCALARPHASE) { // Same value for x and y
- itsParms(1, ant, tf) = polar(1., parmvalues[0][ant][tf]);
- } else { // Different value for x and y
- itsParms(1, ant, tf) = polar(1., parmvalues[1][ant][tf]);
- }
- }
- else if (itsCorrectType==AMPLITUDE || itsCorrectType==SCALARAMPLITUDE) {
- itsParms(0, ant, tf) = parmvalues[0][ant][tf];
- if (itsCorrectType==SCALARAMPLITUDE) { // Same value for x and y
- itsParms(1, ant, tf) = parmvalues[0][ant][tf];
- } else { // Different value for x and y
- itsParms(1, ant, tf) = parmvalues[1][ant][tf];
- }
- }
-
- // Invert
- if (itsInvert) {
- if (itsParms.shape()[0]==2) {
- itsParms(0, ant, tf) = 1./itsParms(0, ant, tf);
- itsParms(1, ant, tf) = 1./itsParms(1, ant, tf);
- } else if (itsCorrectType==FULLJONES) {
- invert(&itsParms(0, ant, tf),itsSigmaMMSE);
- } else {
- ASSERT (itsCorrectType==ROTATIONMEASURE || itsCorrectType==ROTATIONANGLE);
- // rotationmeasure and commonrotationangle are already inverted above
- }
- }
- }
- }
- }
-
- uint ApplyCal::nPol(const string& parmName) {
- if (itsUseH5Parm) {
- if (!itsSolTab.hasAxis("pol")) {
- return 1;
- } else {
- return itsSolTab.getAxis("pol").size;
- }
- } else { // Use ParmDB
- if (itsParmDB->getNames(parmName+":0:*").empty() &&
- itsParmDB->getDefNames(parmName+":0:*").empty() ) {
- return 1;
- } else {
- return 2;
- }
- }
- }
-
- void ApplyCal::initDataArrays() {
- uint numAnts=info().antennaNames().size();
- uint tfDomainSize=itsTimeSlotsPerParmUpdate*info().chanFreqs().size();
-
- uint numParms;
- if (itsCorrectType==FULLJONES ||
- itsCorrectType==ROTATIONANGLE ||
- itsCorrectType==ROTATIONMEASURE) {
- numParms = 4;
- }
- else {
- numParms = 2;
- }
-
- itsParms.resize(numParms, numAnts, tfDomainSize);
- }
-
void ApplyCal::applyDiag (const DComplex* gainA, const DComplex* gainB,
Complex* vis, float* weight, bool* flag,
uint bl, uint chan, bool updateWeights,
@@ -844,13 +295,5 @@ namespace LOFAR {
weight[3]=1./weight[3];
}
- void ApplyCal::showCounts (std::ostream& os) const
- {
- os << endl << "Flags set by ApplyCal " << itsName;
- os << endl << "=======================" << endl;
- itsFlagCounter.showBaseline (os, itsCount);
- itsFlagCounter.showChannel (os, itsCount);
- }
-
} //# end namespace
}
diff --git a/CEP/DP3/DPPP/src/CMakeLists.txt b/CEP/DP3/DPPP/src/CMakeLists.txt
index b2d33f4dcf6c7781c849c34e23265db818acf5e7..0d6115d9b0dc0340ca3a0c5ffced51d7f279e153 100644
--- a/CEP/DP3/DPPP/src/CMakeLists.txt
+++ b/CEP/DP3/DPPP/src/CMakeLists.txt
@@ -18,7 +18,7 @@ lofar_add_library(dppp
ModelComponent.cc PointSource.cc GaussianSource.cc Patch.cc
ModelComponentVisitor.cc GainCal.cc StefCal.cc
DemixerNew.cc DemixInfo.cc DemixWorker.cc
- Predict.cc
+ Predict.cc OneApplyCal.cc
ApplyBeam.cc
PhaseFitter.cc H5Parm.cc SolTab.cc DummyStep.cc H5ParmPredict.cc
)
diff --git a/CEP/DP3/DPPP/src/DummyStep.cc b/CEP/DP3/DPPP/src/DummyStep.cc
index 68bab17a8685e76a142d4f8756ee4edf50017c3b..4fb62624db78903e2161ae0d27f0a0446404d6c9 100644
--- a/CEP/DP3/DPPP/src/DummyStep.cc
+++ b/CEP/DP3/DPPP/src/DummyStep.cc
@@ -35,7 +35,6 @@
#include <vector>
using namespace casacore;
-using namespace LOFAR::BBS;
namespace LOFAR {
namespace DPPP {
@@ -43,6 +42,7 @@ namespace LOFAR {
DummyStep::DummyStep (DPInput* input,
const ParameterSet& parset,
const string& prefix)
+ : itsInput(input)
{
}
diff --git a/CEP/DP3/DPPP/src/H5ParmPredict.cc b/CEP/DP3/DPPP/src/H5ParmPredict.cc
index 8714455202a60bdbd2938acfe262fa0d0978d867..33b0a2dc1b2666d3f0afeb4bc112f2975276781e 100644
--- a/CEP/DP3/DPPP/src/H5ParmPredict.cc
+++ b/CEP/DP3/DPPP/src/H5ParmPredict.cc
@@ -137,6 +137,7 @@ namespace LOFAR {
void H5ParmPredict::finish()
{
// Let the next steps finish.
+ itsPredictSteps[0]->finish();
getNextStep()->finish();
}
} //# end namespace
diff --git a/CEP/DP3/DPPP/src/OneApplyCal.cc b/CEP/DP3/DPPP/src/OneApplyCal.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c602eece041bfaf50d848bf64c48982f0db6dc86
--- /dev/null
+++ b/CEP/DP3/DPPP/src/OneApplyCal.cc
@@ -0,0 +1,713 @@
+//# OneApplyCal.cc: DPPP step class to apply a calibration correction to the data
+//# Copyright (C) 2013
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id: OneApplyCal.cc 21598 2012-07-16 08:07:34Z diepen $
+//#
+//# @author Tammo Jan Dijkema
+
+#include <lofar_config.h>
+#include <DPPP/OneApplyCal.h>
+#include <DPPP/ApplyCal.h>
+#include <DPPP/DPBuffer.h>
+#include <DPPP/DPInfo.h>
+#include <DPPP/MSReader.h>
+#include <Common/ParameterSet.h>
+#include <Common/StringUtil.h>
+#include <Common/LofarLogger.h>
+#include <casacore/casa/Arrays/ArrayMath.h>
+#include <casacore/casa/OS/File.h>
+#include <iostream>
+#include <iomanip>
+
+using namespace casacore;
+using namespace LOFAR::BBS;
+
+/// Look at BBSKernel MeasurementExprLOFARUtil.cc and Apply.cc
+
+namespace LOFAR {
+ namespace DPPP {
+
+ OneApplyCal::OneApplyCal (DPInput* input,
+ const ParameterSet& parset,
+ const string& prefix,
+ const string& defaultPrefix,
+ bool substep,
+ string predictDirection
+ )
+ : itsInput (input),
+ itsName (prefix),
+ itsParmDBName (
+ parset.isDefined(prefix+"parmdb") ?
+ parset.getString(prefix + "parmdb") :
+ parset.getString(defaultPrefix + "parmdb")),
+ itsUseH5Parm (itsParmDBName.find(".h5") != string::npos),
+ itsTimeSlotsPerParmUpdate (
+ parset.isDefined(prefix + "timeslotsperparmupdate") ?
+ parset.getInt (prefix + "timeslotsperparmupdate") :
+ parset.getInt (defaultPrefix + "timeslotsperparmupdate", 500)),
+ itsSigmaMMSE (
+ parset.isDefined(prefix + "MMSE.Sigma") ?
+ parset.getDouble(prefix + "MMSE.Sigma") :
+ parset.getDouble(defaultPrefix + "MMSE.Sigma", 0.)),
+ itsUpdateWeights (
+ parset.isDefined(prefix + "updateweights") ?
+ parset.getBool (prefix + "updateweights") :
+ parset.getBool (defaultPrefix + "updateweights", false)),
+ itsCount (0),
+ itsTimeStep (0),
+ itsNCorr (0),
+ itsTimeInterval (-1),
+ itsLastTime (-1),
+ itsUseAP (false)
+ {
+
+ ASSERT (!itsParmDBName.empty());
+
+ if (substep) {
+ itsInvert=false;
+ } else {
+ itsInvert = (parset.isDefined(prefix + "invert") ?
+ parset.getBool (prefix + "invert") :
+ parset.getBool (defaultPrefix + "invert", true));
+ }
+
+ if (itsUseH5Parm) {
+ string directionStr;
+ directionStr = (parset.isDefined(prefix + "direction") ?
+ parset.getString(prefix + "direction") :
+ parset.getString(defaultPrefix + "direction",
+ predictDirection));
+ itsH5Parm = H5Parm(itsParmDBName);
+
+ string solTabName = (parset.isDefined(prefix + "correction") ?
+ parset.getString(prefix + "correction") :
+ parset.getString(defaultPrefix + "correction"));
+
+ itsSolTab = itsH5Parm.getSolTab(solTabName);
+ itsCorrectType = stringToCorrectType(itsSolTab.getType());
+ if (itsCorrectType==PHASE && nPol("")==1) {
+ itsCorrectType = SCALARPHASE;
+ }
+ if (itsCorrectType==AMPLITUDE && nPol("")==1) {
+ itsCorrectType = SCALARAMPLITUDE;
+ }
+ if (directionStr=="") {
+ ASSERT(!itsSolTab.hasAxis("dir") || itsSolTab.getAxis("dir").size==1);
+ } else {
+ itsDirection = itsSolTab.getDirIndex(directionStr);
+ }
+ } else {
+ string correctTypeStr = toLower(
+ parset.isDefined(prefix + "correction") ?
+ parset.getString(prefix + "correction") :
+ parset.getString (defaultPrefix + "correction", "gain"));
+ itsCorrectType = stringToCorrectType(correctTypeStr);
+ }
+
+ if (itsCorrectType==FULLJONES && itsUpdateWeights) {
+ ASSERTSTR (itsInvert, "Updating weights has not been implemented for invert=false and fulljones");
+ }
+ }
+
+ string OneApplyCal::correctTypeToString(CorrectType ct) {
+ if (ct==GAIN) return "gain";
+ if (ct==FULLJONES) return "fulljones";
+ if (ct==TEC) return "tec";
+ if (ct==CLOCK) return "clock";
+ if (ct==SCALARPHASE) return "scalarphase";
+ if (ct==SCALARAMPLITUDE) return "scalaramplitude";
+ if (ct==ROTATIONANGLE) return "rotationangle";
+ if (ct==ROTATIONMEASURE) return "rotationmeasure";
+ if (ct==PHASE) return "phase";
+ if (ct==AMPLITUDE) return "amplitude";
+ THROW(Exception, "Unknown correction type: "<<ct);
+ return "";
+ }
+
+ OneApplyCal::CorrectType OneApplyCal::stringToCorrectType(const string& ctStr) {
+ if (ctStr=="gain") return GAIN;
+ if (ctStr=="fulljones") return FULLJONES;
+ if (ctStr=="tec") return TEC;
+ if (ctStr=="clock") return CLOCK;
+ if (ctStr=="scalarphase" || ctStr=="commonscalarphase") return SCALARPHASE;
+ if (ctStr=="scalaramplitude" || ctStr=="commonscalaramplitude") return SCALARAMPLITUDE;
+ if (ctStr=="phase") return PHASE;
+ if (ctStr=="amplitude") return AMPLITUDE;
+ if (ctStr=="rotationangle" || ctStr=="commonrotationangle") return ROTATIONANGLE;
+ if (ctStr=="rotationmeasure") return ROTATIONMEASURE;
+ THROW(Exception, "Unknown correction type: "<<ctStr);
+ return GAIN;
+ }
+
+ OneApplyCal::~OneApplyCal()
+ {}
+
+ void OneApplyCal::updateInfo (const DPInfo& infoIn)
+ {
+ info() = infoIn;
+ info().setNeedVisData();
+ info().setWriteData();
+ info().setWriteFlags();
+ if (itsUpdateWeights) {
+ info().setWriteWeights();
+ }
+ itsTimeInterval = infoIn.timeInterval();
+ itsNCorr = infoIn.ncorr();
+
+ ASSERT(itsNCorr==4);
+
+ if (!itsUseH5Parm) { // Use ParmDB
+ itsParmDB.reset(new BBS::ParmFacade(itsParmDBName));
+ }
+
+ // Detect if full jones solutions are present
+ if (!itsUseH5Parm &&
+ (itsCorrectType == GAIN || itsCorrectType==FULLJONES) &&
+ (itsParmDB->getNames("Gain:0:1:*").size() +
+ itsParmDB->getDefNames("Gain:0:1:*").size() >0 )) {
+ itsCorrectType=FULLJONES;
+ }
+
+ // Detect if solutions are saved as Real/Imag or Ampl/Phase
+ if (itsCorrectType == GAIN || itsCorrectType == FULLJONES ){
+ if (itsUseH5Parm) {
+ // H5Parm uses amplitude / phase by definition
+ itsUseAP = true;
+ } else {
+ // Determine from values present in parmdb what to use
+ if (!itsParmDB->getNames("Gain:0:0:Real*").empty()) {
+ // Values with :Real present
+ itsUseAP = false;
+ } else if (!itsParmDB->getNames("Gain:0:0:Ampl*").empty() ||
+ !itsParmDB->getNames("Phase:0:0:Ampl*").empty()) {
+ // Values with :Ampl present
+ itsUseAP = true;
+ } else if (!itsParmDB->getDefNames("Gain:0:0:Real*").empty()) {
+ // Defvalues with :Real present
+ itsUseAP = false;
+ } else if (!itsParmDB->getDefNames("Gain:0:0:Ampl*").empty() ||
+ !itsParmDB->getDefNames("Gain:0:0:Phase*").empty()) {
+ // Defvalues with :Ampl present
+ itsUseAP = true;
+ } else {
+ THROW (Exception, "No gains found in parmdb "+itsParmDBName);
+ }
+ }
+ }
+
+ if (itsCorrectType == GAIN) {
+ if (itsUseAP) {
+ itsParmExprs.push_back("Gain:0:0:Ampl");
+ itsParmExprs.push_back("Gain:0:0:Phase");
+ itsParmExprs.push_back("Gain:1:1:Ampl");
+ itsParmExprs.push_back("Gain:1:1:Phase");
+ } else {
+ itsParmExprs.push_back("Gain:0:0:Real");
+ itsParmExprs.push_back("Gain:0:0:Imag");
+ itsParmExprs.push_back("Gain:1:1:Real");
+ itsParmExprs.push_back("Gain:1:1:Imag");
+ }
+ } else if (itsCorrectType == FULLJONES) {
+ if (itsUseAP) {
+ itsParmExprs.push_back("Gain:0:0:Ampl");
+ itsParmExprs.push_back("Gain:0:0:Phase");
+ itsParmExprs.push_back("Gain:0:1:Ampl");
+ itsParmExprs.push_back("Gain:0:1:Phase");
+ itsParmExprs.push_back("Gain:1:0:Ampl");
+ itsParmExprs.push_back("Gain:1:0:Phase");
+ itsParmExprs.push_back("Gain:1:1:Ampl");
+ itsParmExprs.push_back("Gain:1:1:Phase");
+ } else {
+ itsParmExprs.push_back("Gain:0:0:Real");
+ itsParmExprs.push_back("Gain:0:0:Imag");
+ itsParmExprs.push_back("Gain:0:1:Real");
+ itsParmExprs.push_back("Gain:0:1:Imag");
+ itsParmExprs.push_back("Gain:1:0:Real");
+ itsParmExprs.push_back("Gain:1:0:Imag");
+ itsParmExprs.push_back("Gain:1:1:Real");
+ itsParmExprs.push_back("Gain:1:1:Imag");
+ }
+ } else if (itsCorrectType == TEC) {
+ if (nPol("TEC")==1) {
+ itsParmExprs.push_back("TEC");
+ }
+ else {
+ itsParmExprs.push_back("TEC:0");
+ itsParmExprs.push_back("TEC:1");
+ }
+ } else if (itsCorrectType == CLOCK) {
+ if (nPol("Clock")==1) {
+ itsParmExprs.push_back("Clock");
+ }
+ else {
+ itsParmExprs.push_back("Clock:0");
+ itsParmExprs.push_back("Clock:1");
+ }
+ } else if (itsCorrectType == ROTATIONANGLE) {
+ itsParmExprs.push_back("{Common,}RotationAngle");
+ } else if (itsCorrectType == SCALARPHASE) {
+ itsParmExprs.push_back("{Common,}ScalarPhase");
+ } else if (itsCorrectType == ROTATIONMEASURE) {
+ itsParmExprs.push_back("RotationMeasure");
+ } else if (itsCorrectType == SCALARAMPLITUDE) {
+ itsParmExprs.push_back("{Common,}ScalarAmplitude");
+ } else if (itsCorrectType == PHASE) {
+ ASSERT(itsUseH5Parm);
+ itsParmExprs.push_back("Phase:0");
+ itsParmExprs.push_back("Phase:1");
+ } else if (itsCorrectType == AMPLITUDE) {
+ ASSERT(itsUseH5Parm);
+ itsParmExprs.push_back("Amplitude:0");
+ itsParmExprs.push_back("Amplitude:1");
+ } else {
+ THROW (Exception, "Correction type "<<
+ correctTypeToString(itsCorrectType)<<" is unknown");
+ }
+
+ initDataArrays();
+ itsFlagCounter.init(getInfo());
+
+ // Check that channels are evenly spaced
+ if (info().nchan()>1) {
+ Vector<Double> upFreq = info().chanFreqs()(
+ Slicer(IPosition(1,1),
+ IPosition(1,info().nchan()-1)));
+ Vector<Double> lowFreq = info().chanFreqs()(
+ Slicer(IPosition(1,0),
+ IPosition(1,info().nchan()-1)));
+ Double freqstep0=upFreq(0)-lowFreq(0);
+ // Compare up to 1kHz accuracy
+ bool regularChannels=allNearAbs(upFreq-lowFreq, freqstep0, 1.e3) &&
+ allNearAbs(info().chanWidths(),
+ info().chanWidths()(0), 1.e3);
+ ASSERTSTR(regularChannels,
+ "ApplyCal requires evenly spaced channels.");
+ }
+ }
+
+ void OneApplyCal::show (std::ostream& os) const
+ {
+ os << "ApplyCal " << itsName << std::endl;
+ if (itsUseH5Parm) {
+ os << " H5Parm: " << itsParmDBName <<endl;
+ } else {
+ os << " parmdb: " << itsParmDBName << endl;
+ }
+ os << " correction: " << correctTypeToString(itsCorrectType) << endl;
+ if (itsCorrectType==GAIN || itsCorrectType==FULLJONES) {
+ os << " Ampl/Phase: " << boolalpha << itsUseAP << endl;
+ }
+ os << " update weights: " << boolalpha << itsUpdateWeights << endl;
+ os << " sigmaMMSE: " << itsSigmaMMSE << endl;
+ os << " invert: " << boolalpha << itsInvert <<endl;
+ os << " timeSlotsPerParmUpdate: " << itsTimeSlotsPerParmUpdate <<endl;
+ }
+
+ void OneApplyCal::showTimings (std::ostream& os, double duration) const
+ {
+ os << " ";
+ FlagCounter::showPerc1 (os, itsTimer.getElapsed(), duration);
+ os << " OneApplyCal " << itsName << endl;
+ }
+
+ bool OneApplyCal::process (const DPBuffer& bufin)
+ {
+ itsTimer.start();
+ itsBuffer.copy (bufin);
+
+ if (bufin.getTime() > itsLastTime) {
+ updateParms(bufin.getTime());
+ itsTimeStep=0;
+ }
+ else {
+ itsTimeStep++;
+ }
+
+ // Loop through all baselines in the buffer.
+ size_t nbl = itsBuffer.getData().shape()[2];
+
+ Complex* data = itsBuffer.getData().data();
+
+ itsInput->fetchWeights (bufin, itsBuffer, itsTimer);
+ float* weight = itsBuffer.getWeights().data();
+
+ bool* flag = itsBuffer.getFlags().data();
+
+ size_t nchan = itsBuffer.getData().shape()[1];
+
+#pragma omp parallel for
+ for (size_t bl=0; bl<nbl; ++bl) {
+ for (size_t chan=0;chan<nchan;chan++) {
+ uint timeFreqOffset=(itsTimeStep*info().nchan())+chan;
+ uint antA = info().getAnt1()[bl];
+ uint antB = info().getAnt2()[bl];
+ if (itsParms.shape()[0]>2) {
+ ApplyCal::applyFull( &itsParms(0, antA, timeFreqOffset),
+ &itsParms(0, antB, timeFreqOffset),
+ &data[bl * itsNCorr * nchan + chan * itsNCorr ],
+ &weight[bl * itsNCorr * nchan + chan * itsNCorr ],
+ &flag[ bl * itsNCorr * nchan + chan * itsNCorr ],
+ bl, chan, itsUpdateWeights, itsFlagCounter);
+ }
+ else {
+ ApplyCal::applyDiag( &itsParms(0, antA, timeFreqOffset),
+ &itsParms(0, antB, timeFreqOffset),
+ &data[bl * itsNCorr * nchan + chan * itsNCorr ],
+ &weight[bl * itsNCorr * nchan + chan * itsNCorr ],
+ &flag[ bl * itsNCorr * nchan + chan * itsNCorr ],
+ bl, chan, itsUpdateWeights, itsFlagCounter);
+ }
+ }
+ }
+
+ itsTimer.stop();
+ getNextStep()->process(itsBuffer);
+
+ itsCount++;
+ return true;
+ }
+
+ void OneApplyCal::finish()
+ {
+ // Let the next steps finish.
+ getNextStep()->finish();
+ }
+
+
+ void OneApplyCal::updateParms (const double bufStartTime)
+ {
+ uint numAnts = info().antennaNames().size();
+
+ // itsParms contains the parameters to a grid, first for all parameters
+ // (e.g. Gain:0:0 and Gain:1:1), next all antennas, next over freq * time
+ // as returned by ParmDB
+ vector<vector<vector<double> > > parmvalues;
+ parmvalues.resize(itsParmExprs.size());
+ for (size_t i=0;i<parmvalues.size();++i) {
+ parmvalues[i].resize(numAnts);
+ }
+
+ uint numFreqs (info().chanFreqs().size());
+ double freqInterval (info().chanWidths()[0]);
+ if (numFreqs>1) { // Handle data with evenly spaced gaps between channels
+ freqInterval = info().chanFreqs()[1]-info().chanFreqs()[0];
+ }
+ double minFreq (info().chanFreqs()[0]-0.5*freqInterval);
+ double maxFreq (info().chanFreqs()[numFreqs-1]+0.5*freqInterval);
+ itsLastTime = bufStartTime - 0.5*itsTimeInterval +
+ itsTimeSlotsPerParmUpdate * itsTimeInterval;
+ uint numTimes = itsTimeSlotsPerParmUpdate;
+
+ double lastMSTime = info().startTime() + info().ntime() * itsTimeInterval;
+ if (itsLastTime > lastMSTime && !nearAbs(itsLastTime, lastMSTime, 1.e-3)) {
+ itsLastTime = lastMSTime;
+ numTimes = info().ntime() % itsTimeSlotsPerParmUpdate;
+ }
+
+ map<string, vector<double> > parmMap;
+ map<string, vector<double> >::iterator parmIt;
+
+ uint tfDomainSize=numTimes*numFreqs;
+
+ // Fill parmvalues here, get raw data from H5Parm or ParmDB
+ if (itsUseH5Parm) {
+#pragma omp critical(updateH5ParmValues)
+{
+ // TODO: understand polarization etc.
+ ASSERT(itsParmExprs.size()==1 || itsParmExprs.size()==2);
+ hsize_t startTime = itsSolTab.getTimeIndex(bufStartTime);
+ hsize_t startFreq = itsSolTab.getFreqIndex(info().chanFreqs()[0]);
+ for (uint ant = 0; ant < numAnts; ++ant) {
+
+ uint freqUpsampleFactor = numFreqs;
+ if (itsSolTab.hasAxis("freq") && itsSolTab.getAxis("freq").size > 1) {
+ double h5freqinterval = itsSolTab.getFreqInterval();
+ freqUpsampleFactor = h5freqinterval/info().chanWidths()[0] + 0.5; // Round;
+ ASSERT(near(h5freqinterval, freqUpsampleFactor*info().chanWidths()[0],1.e-5));
+ }
+
+ uint timeUpsampleFactor = numTimes;
+ if (itsSolTab.getAxis("time").size > 1) {
+ double h5timeInterval = itsSolTab.getTimeInterval();
+ timeUpsampleFactor = h5timeInterval/itsTimeInterval+0.5; // Round
+ ASSERT(near(h5timeInterval,timeUpsampleFactor*itsTimeInterval,1.e-5));
+ }
+
+ // Figure out whether time or frequency is first axis
+ bool freqvariesfastest = true;
+ if (itsSolTab.hasAxis("freq") &&
+ itsSolTab.getAxisIndex("freq") < itsSolTab.getAxisIndex("time")) {
+ freqvariesfastest = false;
+ }
+ ASSERT(freqvariesfastest);
+
+ // Take the ceiling of numTimes/timeUpsampleFactor, same for freq
+ uint numTimesInH5Parm = (numTimes+timeUpsampleFactor-1)/timeUpsampleFactor;
+ uint numFreqsInH5Parm = (numFreqs+freqUpsampleFactor-1)/freqUpsampleFactor;
+
+ for (uint pol=0; pol<itsParmExprs.size(); ++pol) {
+ vector<double> rawsols;
+ rawsols = itsSolTab.getValues(info().antennaNames()[ant],
+ startTime, numTimesInH5Parm, 1,
+ startFreq, numFreqsInH5Parm, 1, pol, itsDirection);
+
+ parmvalues[pol][ant].resize(tfDomainSize);
+
+ size_t tf=0;
+ for (uint t=0; t<numTimesInH5Parm; ++t) {
+ for (uint ti=0; ti<timeUpsampleFactor; ++ti) {
+ for (uint f=0; f<numFreqsInH5Parm; ++f) {
+ for (uint fi=0; fi<freqUpsampleFactor; ++fi) {
+ if (tf<tfDomainSize) {
+ parmvalues[pol][ant][tf++] = rawsols[t*numFreqsInH5Parm + f];
+ }
+ }
+ }
+ }
+ }
+ ASSERT(tf==tfDomainSize);
+ }
+ }
+} // End pragma omp critical
+ } else { // Use ParmDB
+ for (uint parmExprNum = 0; parmExprNum<itsParmExprs.size();++parmExprNum) {
+ // parmMap contains parameter values for all antennas
+ parmMap = itsParmDB->getValuesMap( itsParmExprs[parmExprNum] + "*",
+ minFreq, maxFreq, freqInterval,
+ bufStartTime-0.5*itsTimeInterval, itsLastTime,
+ itsTimeInterval, true);
+
+ // Resolve {Common,}Bla to CommonBla or Bla
+ if (!parmMap.empty() &&
+ itsParmExprs[parmExprNum].find("{") != std::string::npos) {
+ uint colonPos = (parmMap.begin()->first).find(":");
+ itsParmExprs[parmExprNum] = (parmMap.begin()->first).substr(0, colonPos);
+ }
+
+ for (uint ant = 0; ant < numAnts; ++ant) {
+ parmIt = parmMap.find(
+ itsParmExprs[parmExprNum] + ":" + info().antennaNames()[ant]);
+
+ if (parmIt != parmMap.end()) {
+ parmvalues[parmExprNum][ant].swap(parmIt->second);
+ ASSERT(parmvalues[parmExprNum][ant].size()==tfDomainSize);
+ } else {// No value found, try default
+ Array<double> defValues;
+ double defValue;
+
+ if (itsParmDB->getDefValues(itsParmExprs[parmExprNum] + ":" +
+ info().antennaNames()[ant]).size()==1) { // Default for antenna
+ itsParmDB->getDefValues(itsParmExprs[parmExprNum] + ":" +
+ info().antennaNames()[ant]).get(0,defValues);
+ ASSERT(defValues.size()==1);
+ defValue=defValues.data()[0];
+ }
+ else if (itsParmDB->getDefValues(itsParmExprs[parmExprNum]).size()
+ == 1) { //Default value
+ itsParmDB->getDefValues(itsParmExprs[parmExprNum]).get(0,defValues);
+ ASSERT(defValues.size()==1);
+ defValue=defValues.data()[0];
+ } else if (itsParmExprs[parmExprNum].substr(0,5)=="Gain:") {
+ defValue=0.;
+ }
+ else {
+ THROW (Exception, "No parameter value found for "+
+ itsParmExprs[parmExprNum]+":"+info().antennaNames()[ant]);
+ }
+
+ parmvalues[parmExprNum][ant].resize(tfDomainSize);
+ for (uint tf=0; tf<tfDomainSize;++tf) {
+ parmvalues[parmExprNum][ant][tf]=defValue;
+ }
+ }
+ }
+ }
+ }
+
+ ASSERT(parmvalues[0][0].size() <= tfDomainSize); // Catches multiple matches
+
+ double freq;
+
+ // Make parameters complex
+ for (uint tf=0;tf<tfDomainSize;++tf) {
+ for (uint ant=0;ant<numAnts;++ant) {
+
+ freq=info().chanFreqs()[tf % numFreqs];
+
+ if (itsCorrectType==GAIN) {
+ if (itsUseAP) { // Data as Amplitude / Phase
+ itsParms(0, ant, tf) = polar(parmvalues[0][ant][tf],
+ parmvalues[1][ant][tf]);
+ itsParms(1, ant, tf) = polar(parmvalues[2][ant][tf],
+ parmvalues[3][ant][tf]);
+ } else { // Data as Real / Imaginary
+ itsParms(0, ant, tf) = DComplex(parmvalues[0][ant][tf],
+ parmvalues[1][ant][tf]);
+ itsParms(1, ant, tf) = DComplex(parmvalues[2][ant][tf],
+ parmvalues[3][ant][tf]);
+ }
+ }
+ else if (itsCorrectType==FULLJONES) {
+ if (itsUseAP) { // Data as Amplitude / Phase
+ itsParms(0, ant, tf) = polar(parmvalues[0][ant][tf],
+ parmvalues[1][ant][tf]);
+ itsParms(1, ant, tf) = polar(parmvalues[2][ant][tf],
+ parmvalues[3][ant][tf]);
+ itsParms(2, ant, tf) = polar(parmvalues[4][ant][tf],
+ parmvalues[5][ant][tf]);
+ itsParms(3, ant, tf) = polar(parmvalues[6][ant][tf],
+ parmvalues[7][ant][tf]);
+ } else { // Data as Real / Imaginary
+ itsParms(0, ant, tf) = DComplex(parmvalues[0][ant][tf],
+ parmvalues[1][ant][tf]);
+ itsParms(1, ant, tf) = DComplex(parmvalues[2][ant][tf],
+ parmvalues[3][ant][tf]);
+ itsParms(2, ant, tf) = DComplex(parmvalues[4][ant][tf],
+ parmvalues[5][ant][tf]);
+ itsParms(3, ant, tf) = DComplex(parmvalues[6][ant][tf],
+ parmvalues[7][ant][tf]);
+ }
+ }
+ else if (itsCorrectType==TEC) {
+ itsParms(0, ant, tf)=polar(1.,
+ parmvalues[0][ant][tf] * -8.44797245e9 / freq);
+ if (itsParmExprs.size() == 1) { // No TEC:0, only TEC:
+ itsParms(1, ant, tf)=polar(1.,
+ parmvalues[0][ant][tf] * -8.44797245e9 / freq);
+ }
+ else { // TEC:0 and TEC:1
+ itsParms(1, ant, tf)=polar(1.,
+ parmvalues[1][ant][tf] * -8.44797245e9 / freq);
+ }
+ }
+ else if (itsCorrectType==CLOCK) {
+ itsParms(0, ant, tf)=polar(1.,
+ parmvalues[0][ant][tf] * freq * casacore::C::_2pi);
+ if (itsParmExprs.size() == 1) { // No Clock:0, only Clock:
+ itsParms(1, ant, tf)=polar(1.,
+ parmvalues[0][ant][tf] * freq * casacore::C::_2pi);
+ }
+ else { // Clock:0 and Clock:1
+ itsParms(1, ant, tf)=polar(1.,
+ parmvalues[1][ant][tf] * freq * casacore::C::_2pi);
+ }
+ }
+ else if (itsCorrectType==ROTATIONANGLE) {
+ double phi=parmvalues[0][ant][tf];
+ if (itsInvert) {
+ phi = -phi;
+ }
+ double sinv=sin(phi);
+ double cosv=cos(phi);
+ itsParms(0, ant, tf) = cosv;
+ itsParms(1, ant, tf) = -sinv;
+ itsParms(2, ant, tf) = sinv;
+ itsParms(3, ant, tf) = cosv;
+ }
+ else if (itsCorrectType==ROTATIONMEASURE) {
+ double lambda2 = casacore::C::c / freq;
+ lambda2 *= lambda2;
+ double chi = parmvalues[0][ant][tf] * lambda2;
+ if (itsInvert) {
+ chi = -chi;
+ }
+ double sinv = sin(chi);
+ double cosv = cos(chi);
+ itsParms(0, ant, tf) = cosv;
+ itsParms(1, ant, tf) = -sinv;
+ itsParms(2, ant, tf) = sinv;
+ itsParms(3, ant, tf) = cosv;
+ }
+ else if (itsCorrectType==PHASE || itsCorrectType==SCALARPHASE) {
+ itsParms(0, ant, tf) = polar(1., parmvalues[0][ant][tf]);
+ if (itsCorrectType==SCALARPHASE) { // Same value for x and y
+ itsParms(1, ant, tf) = polar(1., parmvalues[0][ant][tf]);
+ } else { // Different value for x and y
+ itsParms(1, ant, tf) = polar(1., parmvalues[1][ant][tf]);
+ }
+ }
+ else if (itsCorrectType==AMPLITUDE || itsCorrectType==SCALARAMPLITUDE) {
+ itsParms(0, ant, tf) = parmvalues[0][ant][tf];
+ if (itsCorrectType==SCALARAMPLITUDE) { // Same value for x and y
+ itsParms(1, ant, tf) = parmvalues[0][ant][tf];
+ } else { // Different value for x and y
+ itsParms(1, ant, tf) = parmvalues[1][ant][tf];
+ }
+ }
+
+ // Invert
+ if (itsInvert) {
+ if (itsParms.shape()[0]==2) {
+ itsParms(0, ant, tf) = 1./itsParms(0, ant, tf);
+ itsParms(1, ant, tf) = 1./itsParms(1, ant, tf);
+ } else if (itsCorrectType==FULLJONES) {
+ ApplyCal::invert(&itsParms(0, ant, tf),itsSigmaMMSE);
+ } else {
+ ASSERT (itsCorrectType==ROTATIONMEASURE || itsCorrectType==ROTATIONANGLE);
+ // rotationmeasure and commonrotationangle are already inverted above
+ }
+ }
+ }
+ }
+ }
+
+ uint OneApplyCal::nPol(const string& parmName) {
+ if (itsUseH5Parm) {
+ if (!itsSolTab.hasAxis("pol")) {
+ return 1;
+ } else {
+ return itsSolTab.getAxis("pol").size;
+ }
+ } else { // Use ParmDB
+ if (itsParmDB->getNames(parmName+":0:*").empty() &&
+ itsParmDB->getDefNames(parmName+":0:*").empty() ) {
+ return 1;
+ } else {
+ return 2;
+ }
+ }
+ }
+
+ void OneApplyCal::initDataArrays() {
+ uint numAnts=info().antennaNames().size();
+ uint tfDomainSize=itsTimeSlotsPerParmUpdate*info().chanFreqs().size();
+
+ uint numParms;
+ if (itsCorrectType==FULLJONES ||
+ itsCorrectType==ROTATIONANGLE ||
+ itsCorrectType==ROTATIONMEASURE) {
+ numParms = 4;
+ }
+ else {
+ numParms = 2;
+ }
+
+ itsParms.resize(numParms, numAnts, tfDomainSize);
+ }
+
+ void OneApplyCal::showCounts (std::ostream& os) const
+ {
+ os << endl << "Flags set by OneApplyCal " << itsName;
+ os << endl << "=======================" << endl;
+ itsFlagCounter.showBaseline (os, itsCount);
+ itsFlagCounter.showChannel (os, itsCount);
+ }
+
+ } //# end namespace
+}
diff --git a/CEP/DP3/DPPP/test/CMakeLists.txt b/CEP/DP3/DPPP/test/CMakeLists.txt
index 5a8912b708dc2651e500f7e8bd3ebd5409529d60..4ec9ebecd794b3d76977d3f3ed69aaec27b4ab96 100644
--- a/CEP/DP3/DPPP/test/CMakeLists.txt
+++ b/CEP/DP3/DPPP/test/CMakeLists.txt
@@ -18,6 +18,7 @@ lofar_add_test(tStationAdder tStationAdder.cc)
lofar_add_test(tScaleData tScaleData.cc)
lofar_add_test(tApplyCal tApplyCal.cc)
lofar_add_test(tApplyCal2)
+lofar_add_test(tMultiApplyCal)
lofar_add_test(tFilter tFilter.cc)
#lofar_add_test(tDemixer tDemixer.cc)
lofar_add_test(tNDPPP tNDPPP.cc)
diff --git a/CEP/DP3/DPPP/test/tApplyCal2.run b/CEP/DP3/DPPP/test/tApplyCal2.run
index f54aeb086f3288bff9b1de8aa04fb9bdd48cef4a..72094bc40ed1cadd54b432214d8dea29cf34a7e2 100755
--- a/CEP/DP3/DPPP/test/tApplyCal2.run
+++ b/CEP/DP3/DPPP/test/tApplyCal2.run
@@ -32,14 +32,18 @@ adddef Gain:1:1:Real values=3.
EOL
echo; echo "Testing without updateweights"
-NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 msout.weightcolumn=WEIGHTS_NEW steps=[applycal] applycal.parmdb=tApplyCal.parmdb showcounts=false
+cmd='NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 msout.weightcolumn=WEIGHTS_NEW steps=[applycal] applycal.parmdb=tApplyCal.parmdb showcounts=false'
+echo $cmd
+eval $cmd
$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=9*DATA3))' > taql.out
diff taql.out taql.ref || exit 1
$taqlexe 'select from tNDPPP-generic.MS where not(all(WEIGHTS_NEW~=WEIGHT_SPECTRUM))' > taql.out
diff taql.out taql.ref || exit 1
echo; echo "Testing with updateweights"
-NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 msout.weightcolumn=WEIGHTS_NEW steps=[applycal] applycal.parmdb=tApplyCal.parmdb showcounts=false applycal.updateweights=true
+cmd='NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 msout.weightcolumn=WEIGHTS_NEW steps=[applycal] applycal.parmdb=tApplyCal.parmdb showcounts=false applycal.updateweights=true'
+echo $cmd
+eval $cmd
$taqlexe 'select from tNDPPP-generic.MS where not(all(WEIGHTS_NEW~=81*WEIGHT_SPECTRUM))' > taql.out
diff taql.out taql.ref || exit 1
@@ -49,7 +53,9 @@ parmdbm <<EOL
open table="tApplyCal.parmdb"
adddef CommonScalarPhase values=0
EOL
-NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=commonscalarphase showcounts=false
+cmd='NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=commonscalarphase showcounts=false'
+echo $cmd
+eval $cmd
$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=DATA3))' > taql.out
diff taql.out taql.ref || exit 1
@@ -59,7 +65,9 @@ parmdbm <<EOL
open table="tApplyCal.parmdb"
adddef ScalarPhase values=0
EOL
-NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=commonscalarphase showcounts=false
+cmd='NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=commonscalarphase showcounts=false'
+echo $cmd
+eval $cmd
$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=DATA3))' > taql.out
diff taql.out taql.ref || exit 1
@@ -76,7 +84,9 @@ add ScalarPhase:RS208HBA values=[0.,0.,0.,0.], domain=[10e6, 200e6, 4472025735,
add ScalarPhase:RS305HBA values=[0.,0.,0.,0.], domain=[10e6, 200e6, 4472025735, 4972025795], shape=[2,2], shape=[2,2]
add ScalarPhase:RS307HBA values=[0.,0.,0.,0.], domain=[10e6, 200e6, 4472025735, 4972025795], shape=[2,2], shape=[2,2]
EOL
-NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=commonscalarphase showcounts=false
+cmd='NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=commonscalarphase showcounts=false'
+echo $cmd
+eval $cmd
$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=DATA3))' > taql.out
diff taql.out taql.ref || exit 1
@@ -93,7 +103,9 @@ add ScalarAmplitude:RS208HBA values=[3.,3.,3.,3.], domain=[10e6, 200e6, 44720257
add ScalarAmplitude:RS305HBA values=[3.,3.,3.,3.], domain=[10e6, 200e6, 4472025735, 4972025795], shape=[2,2], shape=[2,2]
add ScalarAmplitude:RS307HBA values=[3.,3.,3.,3.], domain=[10e6, 200e6, 4472025735, 4972025795], shape=[2,2], shape=[2,2]
EOL
-NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=scalaramplitude showcounts=false
+cmd='NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.parmdb=tApplyCal.parmdb applycal.correction=scalaramplitude showcounts=false'
+echo $cmd
+eval $cmd
$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=9*DATA3))' > taql.out
diff taql.out taql.ref || exit 1
diff --git a/CEP/DP3/DPPP/test/tMultiApplyCal.run b/CEP/DP3/DPPP/test/tMultiApplyCal.run
new file mode 100755
index 0000000000000000000000000000000000000000..31762901986ee3d5555093f4f2774538b5cd84e9
--- /dev/null
+++ b/CEP/DP3/DPPP/test/tMultiApplyCal.run
@@ -0,0 +1,46 @@
+#!/bin/bash
+
+# Get the taql executable and srcdir (script created by cmake's CONFIGURE_FILE).
+source findenv.run_script
+echo "srcdirx=$rt_srcdir"
+
+# Set srcdir if not defined (in case run by hand).
+if test "$srcdir" = ""; then
+ srcdir="$rt_srcdir"
+fi
+
+if test ! -f ${srcdir}/tNDPPP-generic.in_MS.tgz; then
+ exit 3 # untested
+fi
+
+set -e # Stop on any error
+
+rm -rf tMultiApplyCal_tmp
+mkdir -p tMultiApplyCal_tmp
+# Unpack the MS and other files and do the DPPP run.
+cd tMultiApplyCal_tmp
+tar zxf ${srcdir}/tNDPPP-generic.in_MS.tgz
+
+# Create expected taql output.
+echo " select result of 0 rows" > taql.ref
+
+echo; echo "Creating parmdb with defvalues 3"
+parmdbm <<EOL
+open table="tApplyCal.parmdb"
+adddef Gain:0:0:Real values=3.
+adddef Gain:1:1:Real values=3.
+adddef CommonScalarPhase values=0
+EOL
+
+cmd='NDPPP msin=tNDPPP-generic.MS checkparset=1 msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.steps="[gain,csp]" applycal.gain.parmdb=tApplyCal.parmdb applycal.gain.correction=gain applycal.csp.parmdb=tApplyCal.parmdb applycal.csp.correction=commonscalarphase showcounts=false'
+echo $cmd
+eval $cmd
+$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=9*DATA3))' > taql.out
+diff taql.out taql.ref || exit 1
+
+cmd='NDPPP msin=tNDPPP-generic.MS checkparset=1 msout=. msout.datacolumn=DATA3 steps=[applycal] applycal.steps="[gain,csp]" applycal.parmdb=tApplyCal.parmdb applycal.gain.correction=gain applycal.csp.correction=commonscalarphase showcounts=false'
+echo $cmd
+eval $cmd
+$taqlexe 'select from tNDPPP-generic.MS where not(all(DATA~=9*DATA3))' > taql.out
+diff taql.out taql.ref || exit 1
+
diff --git a/CEP/DP3/DPPP/test/tMultiApplyCal.sh b/CEP/DP3/DPPP/test/tMultiApplyCal.sh
new file mode 100755
index 0000000000000000000000000000000000000000..12b39f281e6f9ca66c4a064ebd0f89b5422d5b63
--- /dev/null
+++ b/CEP/DP3/DPPP/test/tMultiApplyCal.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+./runctest.sh tMultiApplyCal
diff --git a/CEP/DP3/DPPP_DDECal/test/tDDECal.run b/CEP/DP3/DPPP_DDECal/test/tDDECal.run
index 62b0dc972cabf60394666cce7ba7d4f4417a7e48..ae6ac580fb043ba43d1a8bf022ba098a3f296875 100755
--- a/CEP/DP3/DPPP_DDECal/test/tDDECal.run
+++ b/CEP/DP3/DPPP_DDECal/test/tDDECal.run
@@ -40,7 +40,7 @@ for caltype in complexgain scalarcomplexgain amplitudeonly scalaramplitude
do
for solint in 0 1 2 4
do
- for nchan in 0 1 2 5
+ for nchan in 1
do
rm -rf instrument.h5 # Remove h5parm if it exists
echo "Calibrate on the original sources, caltype=$caltype"
@@ -106,7 +106,8 @@ cmd="NDPPP checkparset=1 msin=tDDECal.MS msout=. numthreads=4\
steps=[ddecal]\
ddecal.sourcedb=tDDECal.MS/sky\
ddecal.directions=[[center,dec_off],[ra_off],[radec_off]]\
- ddecal.applycal.parmdb=instrument.h5 ddecal.applycal.correction=amplitude000\
+ ddecal.applycal.parmdb=instrument.h5 ddecal.applycal.steps=applyampl\
+ ddecal.applycal.applyampl.correction=amplitude000\
ddecal.h5parm=instrument2.h5 ddecal.mode=scalarcomplexgain"
echo $cmd
$cmd