From d980ffab14a892b62f54c2d3b07ef89e8eceded9 Mon Sep 17 00:00:00 2001
From: Tammo Jan Dijkema <dijkema@astron.nl>
Date: Wed, 20 Jul 2016 14:24:45 +0000
Subject: [PATCH] Task #9697: optionally update weights in DPPP applybeam
---
.gitattributes | 2 +
CEP/DP3/DPPP/include/DPPP/ApplyBeam.h | 6 +-
CEP/DP3/DPPP/include/DPPP/ApplyBeam.tcc | 8 ++-
CEP/DP3/DPPP/include/DPPP/ApplyCal.h | 12 ++++
CEP/DP3/DPPP/src/ApplyBeam.cc | 16 +++++-
CEP/DP3/DPPP/src/ApplyCal.cc | 76 ++++++++++++-------------
CEP/DP3/DPPP/src/Predict.cc | 7 ++-
CEP/DP3/DPPP/test/CMakeLists.txt | 1 +
CEP/DP3/DPPP/test/tApplyBeam.run | 6 ++
CEP/DP3/DPPP/test/tApplyCal2.run | 43 ++++++++++++++
CEP/DP3/DPPP/test/tApplyCal2.sh | 2 +
11 files changed, 131 insertions(+), 48 deletions(-)
create mode 100755 CEP/DP3/DPPP/test/tApplyCal2.run
create mode 100755 CEP/DP3/DPPP/test/tApplyCal2.sh
diff --git a/.gitattributes b/.gitattributes
index 4a0935e794d..025ea37510b 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -895,6 +895,8 @@ CEP/DP3/DPPP/test/tApplyCal.cc -text
CEP/DP3/DPPP/test/tApplyCal.in_parmdb.tgz -text svneol=unset#application/x-gzip
CEP/DP3/DPPP/test/tApplyCal.run -text
CEP/DP3/DPPP/test/tApplyCal.sh -text
+CEP/DP3/DPPP/test/tApplyCal2.run -text
+CEP/DP3/DPPP/test/tApplyCal2.sh -text
CEP/DP3/DPPP/test/tApplyCal_parmdbscript -text
CEP/DP3/DPPP/test/tDemix.in_MS.tgz -text
CEP/DP3/DPPP/test/tDemix.run -text
diff --git a/CEP/DP3/DPPP/include/DPPP/ApplyBeam.h b/CEP/DP3/DPPP/include/DPPP/ApplyBeam.h
index 2f80b2f860f..c27c34e2af0 100644
--- a/CEP/DP3/DPPP/include/DPPP/ApplyBeam.h
+++ b/CEP/DP3/DPPP/include/DPPP/ApplyBeam.h
@@ -82,13 +82,14 @@ namespace LOFAR {
template<typename T>
static void applyBeam(
- const DPInfo& info, double time, T* data0,
+ const DPInfo& info, double time, T* data0, float* weight0,
const StationResponse::vector3r_t& srcdir,
const StationResponse::vector3r_t& refdir,
const StationResponse::vector3r_t& tiledir,
const vector<StationResponse::Station::Ptr>& antBeamInfo,
vector<StationResponse::matrix22c_t>& beamValues,
- bool useChannelFreq, bool invert, int mode=DEFAULT);
+ bool useChannelFreq, bool invert, int mode,
+ bool doUpdateWeights=false);
private:
StationResponse::vector3r_t dir2Itrf(
@@ -100,6 +101,7 @@ namespace LOFAR {
string itsName;
DPBuffer itsBuffer;
bool itsInvert;
+ bool itsUpdateWeights;
bool itsUseChannelFreq;
Position itsPhaseRef;
BeamMode itsMode;
diff --git a/CEP/DP3/DPPP/include/DPPP/ApplyBeam.tcc b/CEP/DP3/DPPP/include/DPPP/ApplyBeam.tcc
index ecaeef750c8..d241f43c03c 100644
--- a/CEP/DP3/DPPP/include/DPPP/ApplyBeam.tcc
+++ b/CEP/DP3/DPPP/include/DPPP/ApplyBeam.tcc
@@ -66,13 +66,13 @@ namespace LOFAR {
// applyBeam is templated on the type of the data, could be double or float
template<typename T>
void ApplyBeam::applyBeam(
- const DPInfo& info, double time, T* data0,
+ const DPInfo& info, double time, T* data0, float* weight0,
const StationResponse::vector3r_t& srcdir,
const StationResponse::vector3r_t& refdir,
const StationResponse::vector3r_t& tiledir,
const vector<StationResponse::Station::Ptr>& antBeamInfo,
vector<StationResponse::matrix22c_t>& beamValues, bool useChannelFreq,
- bool invert, int mode)
+ bool invert, int mode, bool doUpdateWeights)
{
// Get the beam values for each station.
uint nCh = info.chanFreqs().size();
@@ -162,6 +162,10 @@ void ApplyBeam::applyBeam(
data[1] = tmp[0] * r[1] + tmp[1] * r[3];
data[2] = tmp[2] * r[0] + tmp[3] * r[2];
data[3] = tmp[2] * r[1] + tmp[3] * r[3];
+
+ if (doUpdateWeights) {
+ ApplyCal::applyWeights(l, r, weight0 + bl * 4 * nCh + ch * 4);
+ }
}
}
}
diff --git a/CEP/DP3/DPPP/include/DPPP/ApplyCal.h b/CEP/DP3/DPPP/include/DPPP/ApplyCal.h
index 113b39c2a5e..b827eaa81d4 100644
--- a/CEP/DP3/DPPP/include/DPPP/ApplyCal.h
+++ b/CEP/DP3/DPPP/include/DPPP/ApplyCal.h
@@ -92,6 +92,18 @@ namespace LOFAR {
uint bl, uint chan, bool updateWeights,
FlagCounter& flagCounter);
+ // Do the same as the combination of BBS + python script
+ // covariance2weight.py (cookbook), except it stores weights per freq.
+ // The diagonal of covariance matrix is transferred to the weights.
+ // Note that the real covariance (mixing of noise terms after which they
+ // are not independent anymore) is not stored.
+ // The input covariance matrix C is assumed to be diagonal with elements
+ // w_i (the weights), the result the diagonal of
+ // (gainA kronecker gainB^H).C.(gainA kronecker gainB^H)^H
+ static void applyWeights (const casa::DComplex* gainA,
+ const casa::DComplex* gainB,
+ float* weight);
+
private:
// Read parameters from the associated parmdb and store them in itsParms
void updateParms (const double bufStartTime);
diff --git a/CEP/DP3/DPPP/src/ApplyBeam.cc b/CEP/DP3/DPPP/src/ApplyBeam.cc
index 7f9b139c69c..8d1484bb1a0 100644
--- a/CEP/DP3/DPPP/src/ApplyBeam.cc
+++ b/CEP/DP3/DPPP/src/ApplyBeam.cc
@@ -54,10 +54,11 @@ namespace LOFAR {
namespace DPPP {
ApplyBeam::ApplyBeam(DPInput* input, const ParameterSet& parset,
- const string& prefix, bool substep)
+ const string& prefix, bool substep)
:
itsInput(input),
itsName(prefix),
+ itsUpdateWeights(parset.getBool(prefix + "updateweights", false)),
itsUseChannelFreq(parset.getBool(prefix + "usechannelfreq", true)),
itsDebugLevel(parset.getInt(prefix + "debuglevel", 0))
{
@@ -94,6 +95,9 @@ namespace LOFAR {
info() = infoIn;
info().setNeedVisData();
info().setWriteData();
+ if (itsUpdateWeights) {
+ info().setWriteWeights();
+ }
MDirection dirJ2000(
MDirection::Convert(infoIn.phaseCenter(), MDirection::J2000)());
@@ -136,6 +140,7 @@ namespace LOFAR {
os << endl;
os << " use channelfreq: " << boolalpha << itsUseChannelFreq << endl;
os << " invert: " << boolalpha << itsInvert << endl;
+ os << " update weights: " << boolalpha << itsUpdateWeights << endl;
}
void ApplyBeam::showTimings(std::ostream& os, double duration) const
@@ -151,6 +156,11 @@ namespace LOFAR {
itsBuffer.copy (bufin);
Complex* data=itsBuffer.getData().data();
+ if (itsUpdateWeights) {
+ itsInput->fetchWeights (bufin, itsBuffer, itsTimer);
+ }
+ float* weight = itsBuffer.getWeights().data();
+
double time = itsBuffer.getTime();
//Set up directions for beam evaluation
@@ -168,9 +178,9 @@ namespace LOFAR {
uint thread = OpenMP::threadNum();
StationResponse::vector3r_t srcdir = refdir;
- applyBeam(info(), time, data, srcdir, refdir, tiledir,
+ applyBeam(info(), time, data, weight, srcdir, refdir, tiledir,
itsAntBeamInfo[thread], itsBeamValues[thread],
- itsUseChannelFreq, itsInvert, itsMode);
+ itsUseChannelFreq, itsInvert, itsMode, itsUpdateWeights);
itsTimer.stop();
getNextStep()->process(itsBuffer);
diff --git a/CEP/DP3/DPPP/src/ApplyCal.cc b/CEP/DP3/DPPP/src/ApplyCal.cc
index 0e4f28cdca3..08f0dc16ea9 100644
--- a/CEP/DP3/DPPP/src/ApplyCal.cc
+++ b/CEP/DP3/DPPP/src/ApplyCal.cc
@@ -554,7 +554,7 @@ namespace LOFAR {
void ApplyCal::applyFull (const DComplex* gainA, const DComplex* gainB,
Complex* vis, float* weight, bool* flag,
- uint bl, uint chan, bool updateWeights,
+ uint bl, uint chan, bool doUpdateWeights,
FlagCounter& flagCounter) {
DComplex gainAxvis[4];
@@ -595,46 +595,44 @@ namespace LOFAR {
}
}
- // The code below does the same as the combination of BBS + python script
- // covariance2weight.py (cookbook), except it stores weights per freq.
- // The diagonal of covariance matrix is transferred to the weights.
- // Note that the real covariance (mixing of noise terms after which they
- // are not independent anymore) is not stored.
- // The input covariance matrix C is assumed to be diagonal with elements
- // w_i (the weights), the result the diagonal of
- // (gainA kronecker gainB^H).C.(gainA kronecker gainB^H)^H
- if (updateWeights) {
- float cov[4], normGainA[4], normGainB[4];
- for (uint i=0;i<4;++i) {
- cov[i]=1./weight[i];
- normGainA[i]=norm(gainA[i]);
- normGainB[i]=norm(gainB[i]);
- }
+ if (doUpdateWeights) {
+ applyWeights(gainA, gainB, weight);
+ }
+ }
- weight[0]=cov[0]*(normGainA[0]*normGainB[0])
- +cov[1]*(normGainA[0]*normGainB[1])
- +cov[2]*(normGainA[1]*normGainB[0])
- +cov[3]*(normGainA[1]*normGainB[1]);
- weight[0]=1./weight[0];
-
- weight[1]=cov[0]*(normGainA[0]*normGainB[2])
- +cov[1]*(normGainA[0]*normGainB[3])
- +cov[2]*(normGainA[1]*normGainB[2])
- +cov[3]*(normGainA[1]*normGainB[3]);
- weight[1]=1./weight[1];
-
- weight[2]=cov[0]*(normGainA[2]*normGainB[0])
- +cov[1]*(normGainA[2]*normGainB[1])
- +cov[2]*(normGainA[3]*normGainB[0])
- +cov[3]*(normGainA[3]*normGainB[1]);
- weight[2]=1./weight[2];
-
- weight[3]=cov[0]*(normGainA[2]*normGainB[2])
- +cov[1]*(normGainA[2]*normGainB[3])
- +cov[2]*(normGainA[3]*normGainB[2])
- +cov[3]*(normGainA[3]*normGainB[3]);
- weight[3]=1./weight[3];
+ void ApplyCal::applyWeights(const DComplex* gainA,
+ const DComplex* gainB,
+ float* weight) {
+ float cov[4], normGainA[4], normGainB[4];
+ for (uint i=0;i<4;++i) {
+ cov[i]=1./weight[i];
+ normGainA[i]=norm(gainA[i]);
+ normGainB[i]=norm(gainB[i]);
}
+
+ weight[0]=cov[0]*(normGainA[0]*normGainB[0])
+ +cov[1]*(normGainA[0]*normGainB[1])
+ +cov[2]*(normGainA[1]*normGainB[0])
+ +cov[3]*(normGainA[1]*normGainB[1]);
+ weight[0]=1./weight[0];
+
+ weight[1]=cov[0]*(normGainA[0]*normGainB[2])
+ +cov[1]*(normGainA[0]*normGainB[3])
+ +cov[2]*(normGainA[1]*normGainB[2])
+ +cov[3]*(normGainA[1]*normGainB[3]);
+ weight[1]=1./weight[1];
+
+ weight[2]=cov[0]*(normGainA[2]*normGainB[0])
+ +cov[1]*(normGainA[2]*normGainB[1])
+ +cov[2]*(normGainA[3]*normGainB[0])
+ +cov[3]*(normGainA[3]*normGainB[1]);
+ weight[2]=1./weight[2];
+
+ weight[3]=cov[0]*(normGainA[2]*normGainB[2])
+ +cov[1]*(normGainA[2]*normGainB[3])
+ +cov[2]*(normGainA[3]*normGainB[2])
+ +cov[3]*(normGainA[3]*normGainB[3]);
+ weight[3]=1./weight[3];
}
void ApplyCal::showCounts (std::ostream& os) const
diff --git a/CEP/DP3/DPPP/src/Predict.cc b/CEP/DP3/DPPP/src/Predict.cc
index 64db694225c..e8b318a82fc 100644
--- a/CEP/DP3/DPPP/src/Predict.cc
+++ b/CEP/DP3/DPPP/src/Predict.cc
@@ -337,9 +337,12 @@ namespace LOFAR {
MDirection::J2000);
StationResponse::vector3r_t srcdir = dir2Itrf(dir,itsMeasConverters[thread]);
- ApplyBeam::applyBeam(info(), time, data0, srcdir, refdir,
+ float* dummyweight = 0;
+
+ ApplyBeam::applyBeam(info(), time, data0, dummyweight, srcdir, refdir,
tiledir, itsAntBeamInfo[thread],
- itsBeamValues[thread], itsUseChannelFreq, false, itsBeamMode);
+ itsBeamValues[thread], itsUseChannelFreq, false,
+ itsBeamMode, false);
//Add temporary buffer to itsModelVis
std::transform(itsModelVis[thread].data(),
diff --git a/CEP/DP3/DPPP/test/CMakeLists.txt b/CEP/DP3/DPPP/test/CMakeLists.txt
index cbe86586311..f5e4caa66ae 100644
--- a/CEP/DP3/DPPP/test/CMakeLists.txt
+++ b/CEP/DP3/DPPP/test/CMakeLists.txt
@@ -17,6 +17,7 @@ lofar_add_test(tPhaseShift tPhaseShift.cc)
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(tFilter tFilter.cc)
#lofar_add_test(tDemixer tDemixer.cc)
lofar_add_test(tNDPPP tNDPPP.cc)
diff --git a/CEP/DP3/DPPP/test/tApplyBeam.run b/CEP/DP3/DPPP/test/tApplyBeam.run
index 8878b0ad9ec..19daa40a2d1 100755
--- a/CEP/DP3/DPPP/test/tApplyBeam.run
+++ b/CEP/DP3/DPPP/test/tApplyBeam.run
@@ -58,3 +58,9 @@ echo; echo "Test with beammode=ELEMENT"; echo
# Compare the DATA column of the output MS with the BBS reference output.
$taqlexe 'select from outinv.ms t1, tApplyBeam.tab t2 where not all(near(t1.DATA,t2.DATA_ELEMENT,1e-5) || (isnan(t1.DATA) && isnan(t2.DATA_ELEMENT)))' > taql.out
diff taql.out taql.ref || exit 1
+
+echo; echo "Test with updateweights=true"; echo
+../../src/NDPPP msin=tNDPPP-generic.MS msout=. steps=[applybeam] applybeam.updateweights=truue msout.weightcolumn=NEW_WEIGHT_SPECTRUM
+# Check that the weights have changed
+$taqlexe 'select from tNDPPP-generic.MS where all(near(WEIGHT_SPECTRUM, NEW_WEIGHT_SPECTRUM))' > taql.out
+diff taql.out taql.ref || exit 1
diff --git a/CEP/DP3/DPPP/test/tApplyCal2.run b/CEP/DP3/DPPP/test/tApplyCal2.run
new file mode 100755
index 00000000000..3e460f1ffa5
--- /dev/null
+++ b/CEP/DP3/DPPP/test/tApplyCal2.run
@@ -0,0 +1,43 @@
+#!/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
+
+rm -rf tApplyCal2_tmp
+mkdir -p tApplyCal2_tmp
+# Unpack the MS and other files and do the DPPP run.
+cd tApplyCal2_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"
+../../../../ParmDB/src/parmdbm <<EOL
+open table="tApplyCal.parmdb"
+adddef Gain:0:0:Real values=3.
+adddef Gain:1:1:Real values=3.
+EOL
+
+echo; echo "Testing without updateweights"
+../../src/NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 msout.weightcolumn=WEIGHTS_NEW steps=[applycal] applycal.parmdb=tApplyCal.parmdb showcounts=false
+$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"
+../../src/NDPPP msin=tNDPPP-generic.MS msout=. msout.datacolumn=DATA3 msout.weightcolumn=WEIGHTS_NEW steps=[applycal] applycal.parmdb=tApplyCal.parmdb showcounts=false applycal.updateweights=true
+$taqlexe 'select from tNDPPP-generic.MS where not(all(WEIGHTS_NEW~=81*WEIGHT_SPECTRUM))' > taql.out
+diff taql.out taql.ref || exit 1
+
diff --git a/CEP/DP3/DPPP/test/tApplyCal2.sh b/CEP/DP3/DPPP/test/tApplyCal2.sh
new file mode 100755
index 00000000000..2c166e188f8
--- /dev/null
+++ b/CEP/DP3/DPPP/test/tApplyCal2.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+./runctest.sh tApplyCal2
--
GitLab