diff --git a/CEP/DP3/DPPP/include/DPPP/MedFlagger.h b/CEP/DP3/DPPP/include/DPPP/MedFlagger.h
index cacd9d6e51947dfa49ffbb2b83cae26078dcd596..cc9853c3f350a861b8fae694e630de3e9ad9f8f7 100644
--- a/CEP/DP3/DPPP/include/DPPP/MedFlagger.h
+++ b/CEP/DP3/DPPP/include/DPPP/MedFlagger.h
@@ -119,6 +119,9 @@ namespace LOFAR {
uint itsNTimes;
uint itsNTimesDone;
vector<uint> itsFlagCorr;
+ bool itsApplyAutoCorr;
+ vector<int> itsAutoCorrIndex; //# baseline index of autocorrelations
+ uint itsNrAutoCorr;
vector<DPBuffer> itsBuf;
FlagCounter itsFlagCounter;
NSTimer itsTimer;
diff --git a/CEP/DP3/DPPP/src/FlagCounter.cc b/CEP/DP3/DPPP/src/FlagCounter.cc
index 0d7e6dd80be9d16d541143236c7bfd02bb1720e2..281fb05637fef7f62854dca184384b2235d5a8e7 100644
--- a/CEP/DP3/DPPP/src/FlagCounter.cc
+++ b/CEP/DP3/DPPP/src/FlagCounter.cc
@@ -47,7 +47,6 @@ namespace LOFAR {
const casa::Vector<int>& ant2,
int64 ntimes) const
{
- os << "newshow"<<endl;
int64 npoints = ntimes * itsChanCounts.size();
os << endl << "Percentage of visibilities flagged per baseline"
" (antenna pair):";
diff --git a/CEP/DP3/DPPP/src/MedFlagger.cc b/CEP/DP3/DPPP/src/MedFlagger.cc
index b9bede098f13c60d6c419c58c2f05b2b9b802b27..d514633c9ac23f48533defd75a229ee1647b266b 100644
--- a/CEP/DP3/DPPP/src/MedFlagger.cc
+++ b/CEP/DP3/DPPP/src/MedFlagger.cc
@@ -52,6 +52,24 @@ namespace LOFAR {
itsBuf.resize (itsTimeWindow);
itsFlagCorr = parset.getUintVector (prefix+"correlations",
vector<uint>());
+ itsApplyAutoCorr = parset.getBool (prefix+"applyautocorr", false);
+ // Determine baseline indices of autocorrelations.
+ const Vector<int>& ant1 = itsInput->getAnt1();
+ const Vector<int>& ant2 = itsInput->getAnt2();
+ int nant = 1 + std::max (max(ant1), max(ant2));
+ itsAutoCorrIndex.resize (nant);
+ std::fill (itsAutoCorrIndex.begin(), itsAutoCorrIndex.end(), -1);
+ itsNrAutoCorr = 0;
+ for (uint i=0; i<ant1.size(); ++i) {
+ if (ant1[i] == ant2[i]) {
+ itsAutoCorrIndex[ant1[i]] = i;
+ itsNrAutoCorr++;
+ }
+ }
+ if (itsApplyAutoCorr) {
+ ASSERTSTR (itsNrAutoCorr > 0, "applyautocorr=True cannot be used if "
+ "the data does not contain autocorrelations");
+ }
}
MedFlagger::~MedFlagger()
@@ -64,6 +82,8 @@ namespace LOFAR {
os << " timewindow: " << itsTimeWindow << std::endl;
os << " threshold: " << itsThreshold << std::endl;
os << " correlations: " << itsFlagCorr << std::endl;
+ os << " applyautocorr: " << itsApplyAutoCorr
+ << " (nautocorr = " << itsNrAutoCorr << ')' << std::endl;
}
void MedFlagger::showCounts (std::ostream& os) const
@@ -222,6 +242,9 @@ namespace LOFAR {
void MedFlagger::flag (uint index, const vector<uint>& timeEntries)
{
///cout << "flag: " <<itsNTimes<<' '<<itsNTimesDone<<' ' <<index << timeEntries << endl;
+ // Get antenna numbers in case applyautocorr is true.
+ const Vector<int>& ant1 = itsInput->getAnt1();
+ const Vector<int>& ant2 = itsInput->getAnt2();
// Result is 'copy' of the entry at the given time index.
DPBuffer buf (itsBuf[index]);
IPosition shp = buf.getData().shape();
@@ -238,36 +261,80 @@ namespace LOFAR {
float MAD = 1.4826; //# constant determined by Pandey
// Now flag each baseline, channel and correlation for this time window.
for (uint ib=0; ib<nrbl; ++ib) {
- for (uint ic=0; ic<nchan; ++ic) {
- bool corrIsFlagged = false;
- // Iterate over given correlations.
- for (vector<uint>::const_iterator iter = itsFlagCorr.begin();
- iter != itsFlagCorr.end(); ++iter) {
- uint ip = *iter;
- // If one correlation is flagged, all of them will be flagged.
- // So no need to check others.
- if (flagPtr[ip]) {
- corrIsFlagged = true;
- break;
+ // Do only autocorrelations if told so.
+ if (!itsApplyAutoCorr || ant1[ib] == ant2[ib]) {
+ for (uint ic=0; ic<nchan; ++ic) {
+ bool corrIsFlagged = false;
+ // Iterate over given correlations.
+ for (vector<uint>::const_iterator iter = itsFlagCorr.begin();
+ iter != itsFlagCorr.end(); ++iter) {
+ uint ip = *iter;
+ // If one correlation is flagged, all of them will be flagged.
+ // So no need to check others.
+ if (flagPtr[ip]) {
+ corrIsFlagged = true;
+ break;
+ }
+ // Calculate values from the median.
+ computeFactors (timeEntries, ib, ic, ip, nchan, ncorr,
+ Z1, Z2, tempBuf.storage());
+ if (dataPtr[ip] > Z1 + itsThreshold * Z2 * MAD) {
+ corrIsFlagged = true;
+ itsFlagCounter.incrBaseline(ib);
+ itsFlagCounter.incrChannel(ic);
+ itsFlagCounter.incrCorrelation(ip);
+ break;
+ }
}
- // Calculate values from the median.
- computeFactors (timeEntries, ib, ic, ip, nchan, ncorr,
- Z1, Z2, tempBuf.storage());
- if (dataPtr[ip] > Z1 + itsThreshold * Z2 * MAD) {
- corrIsFlagged = true;
- itsFlagCounter.incrBaseline(ib);
- itsFlagCounter.incrChannel(ic);
- itsFlagCounter.incrCorrelation(ip);
- break;
+ if (corrIsFlagged) {
+ for (uint ip=0; ip<ncorr; ++ip) {
+ flagPtr[ip] = true;
+ }
}
+ dataPtr += ncorr;
+ flagPtr += ncorr;
}
- if (corrIsFlagged) {
- for (uint ip=0; ip<ncorr; ++ip) {
- flagPtr[ip] = true;
+ } else {
+ dataPtr += nchan*ncorr;
+ flagPtr += nchan*ncorr;
+ }
+ }
+ // Apply autocorrelations flags if needed.
+ if (itsApplyAutoCorr) {
+ flagPtr = buf.getFlags().data();
+ for (uint ib=0; ib<nrbl; ++ib) {
+ // Flag crosscorr if at least one autocorr is present.
+ if (ant1[ib] != ant2[ib]) {
+ int inx1 = itsAutoCorrIndex[ant1[ib]];
+ int inx2 = itsAutoCorrIndex[ant2[ib]];
+ if (inx1 >= 0 || inx2 >= 0) {
+ // Find flags of the autocorrelations of both antennae.
+ // Use other autocorr if one autocorr does not exist.
+ // In this way inx does not need to be tested in the inner loop.
+ if (inx1 < 0) {
+ inx1 = inx2;
+ } else if (inx2 < 0) {
+ inx2 = inx1;
+ }
+ bool* flagAnt1 = buf.getFlags().data() + inx1*nchan*ncorr;
+ bool* flagAnt2 = buf.getFlags().data() + inx2*nchan*ncorr;
+ // Flag if not flagged yet and if one of autocorr is flagged.
+ for (uint ic=0; ic<nchan; ++ic) {
+ if (!*flagPtr && (*flagAnt1 || *flagAnt2)) {
+ for (uint ip=0; ip<ncorr; ++ip) {
+ flagPtr[ip] = true;
+ }
+ itsFlagCounter.incrBaseline(ib);
+ itsFlagCounter.incrChannel(ic);
+ }
+ flagPtr += ncorr;
+ flagAnt1 += ncorr;
+ flagAnt2 += ncorr;
+ }
}
+ } else {
+ flagPtr += nchan*ncorr;
}
- dataPtr += ncorr;
- flagPtr += ncorr;
}
}
// Process the result in the next step.
diff --git a/CEP/DP3/DPPP/src/NDPPP.parset b/CEP/DP3/DPPP/src/NDPPP.parset
index d55c941da5c8e7c7f2b8ca35cc3a388c9e087d6f..c21ab77baae7cb602da8a0fbf921d907bb1662a2 100644
--- a/CEP/DP3/DPPP/src/NDPPP.parset
+++ b/CEP/DP3/DPPP/src/NDPPP.parset
@@ -12,6 +12,8 @@ msout.overwrite=True
#steps = [flag,avg]
steps=[flag]
+preflag.corrtype = auto
+
avg1.type = squash
avg1.freqstep = 4
avg1.timestep = 2
@@ -24,6 +26,7 @@ flag.type=madflagger
flag.threshold=1
flag.freqwindow=31
flag.timewindow=5
+flag.applyautocorr = True
avg.type = squash
avg.freqstep = 1024
diff --git a/CEP/DP3/DPPP/test/tMedFlagger.cc b/CEP/DP3/DPPP/test/tMedFlagger.cc
index 8573bd5e94f9b13654271e0d485dcf4b827fbccd..35242fb49e96c672f255dccfe52e162586a604a5 100644
--- a/CEP/DP3/DPPP/test/tMedFlagger.cc
+++ b/CEP/DP3/DPPP/test/tMedFlagger.cc
@@ -49,17 +49,46 @@ public:
: itsCount(0), itsNTime(ntime), itsNBl(nant*(nant+1)/2), itsNChan(nchan),
itsNCorr(ncorr), itsFlag(flag)
{
+ // Fill the baseline stations; use 4 stations.
+ // So they are called 00 01 02 03 10 11 12 13 20, etc.
itsAnt1.resize (itsNBl);
itsAnt2.resize (itsNBl);
- int i=0;
- for (int st1=0; st1<nant; ++st1) {
- for (int st2=st1; st2<nant; ++st2) {
- itsAnt1[i] = st1;
- itsAnt2[i] = st2;
- i++;
+ int st1 = 0;
+ int st2 = 0;
+ for (int i=0; i<itsNBl; ++i) {
+ itsAnt1[i] = st1;
+ itsAnt2[i] = st2;
+ if (++st2 == 4) {
+ st2 = 0;
+ if (++st1 == 4) {
+ st1 = 0;
+ }
}
}
-}
+ itsAntNames.resize(4);
+ itsAntNames[0] = "rs01.s01";
+ itsAntNames[1] = "rs02.s01";
+ itsAntNames[2] = "cs01.s01";
+ itsAntNames[3] = "cs01.s02";
+ // Define their positions (more or less WSRT RT0-3).
+ itsAntPos.resize (4);
+ Vector<double> vals(3);
+ vals[0] = 3828763; vals[1] = 442449; vals[2] = 5064923;
+ itsAntPos[0] = MPosition(Quantum<Vector<double> >(vals,"m"),
+ MPosition::ITRF);
+ vals[0] = 3828746; vals[1] = 442592; vals[2] = 5064924;
+ itsAntPos[1] = MPosition(Quantum<Vector<double> >(vals,"m"),
+ MPosition::ITRF);
+ vals[0] = 3828729; vals[1] = 442735; vals[2] = 5064925;
+ itsAntPos[2] = MPosition(Quantum<Vector<double> >(vals,"m"),
+ MPosition::ITRF);
+ vals[0] = 3828713; vals[1] = 442878; vals[2] = 5064926;
+ itsAntPos[3] = MPosition(Quantum<Vector<double> >(vals,"m"),
+ MPosition::ITRF);
+ // Define the frequencies.
+ itsChanFreqs.resize (nchan);
+ indgen (itsChanFreqs, 1050000., 100000.);
+ }
private:
virtual bool process (const DPBuffer&)
{
@@ -100,15 +129,16 @@ private:
bool itsFlag;
};
-// Class to check result of flagged, unaveraged TestInput.
+// Class to check result.
class TestOutput: public DPStep
{
public:
TestOutput(int ntime, int nant, int nchan, int ncorr,
- bool flag)
+ bool flag, bool flagged)
: itsCount(0), itsNTime(ntime), itsNBl(nant*(nant+1)/2), itsNChan(nchan),
itsNCorr(ncorr),
- itsFlag(flag)
+ itsFlag(flag),
+ itsFlagged(flagged)
{}
private:
virtual bool process (const DPBuffer& buf)
@@ -121,7 +151,16 @@ private:
// Check the result.
ASSERT (allNear(real(buf.getData()), real(result), 1e-10));
ASSERT (allNear(imag(buf.getData()), imag(result), 1e-10));
- ASSERT (allEQ(buf.getFlags(), itsFlag));
+ Cube<bool> expFlag(itsNCorr,itsNChan,itsNBl);
+ expFlag = itsFlag;
+ if (itsFlagged) {
+ for (int i=0; i<itsNBl; ++i) {
+ for (int j=0; j<itsNCorr; ++j) {
+ expFlag(j,0,i) = itsFlag || itsCount==0;
+ expFlag(j,itsNChan-1,i) = true;
+ }
+ }
+ }
ASSERT (near(buf.getTime(), 2+5.*itsCount));
++itsCount;
return true;
@@ -142,7 +181,7 @@ private:
int itsCount;
int itsNTime, itsNBl, itsNChan, itsNCorr, itsNAvgTime, itsNAvgChan;
- bool itsFlag;
+ bool itsFlag, itsFlagged;
};
@@ -180,7 +219,27 @@ void test1(int ntime, int nant, int nchan, int ncorr, bool flag, int threshold)
parset.add ("timewindow", "1");
parset.add ("threshold", toString(threshold));
DPStep::ShPtr step2(new MedFlagger(in, parset, ""));
- DPStep::ShPtr step3(new TestOutput(ntime, nant, nchan, ncorr, flag));
+ DPStep::ShPtr step3(new TestOutput(ntime, nant, nchan, ncorr, flag, false));
+ step1->setNextStep (step2);
+ step2->setNextStep (step3);
+ execute (step1);
+}
+
+// Test applyautocorr flagging with or without preflagged points.
+void test2(int ntime, int nant, int nchan, int ncorr, bool flag, int threshold)
+{
+ cout << "test2: ntime=" << ntime << " nrant=" << nant << " nchan=" << nchan
+ << " ncorr=" << ncorr << " threshold=" << threshold << endl;
+ // Create the steps.
+ TestInput* in = new TestInput(ntime, nant, nchan, ncorr, flag);
+ DPStep::ShPtr step1(in);
+ ParameterSet parset;
+ parset.add ("freqwindow", "3");
+ parset.add ("timewindow", "1");
+ parset.add ("threshold", toString(threshold));
+ parset.add ("applyautocorr", "True");
+ DPStep::ShPtr step2(new MedFlagger(in, parset, ""));
+ DPStep::ShPtr step3(new TestOutput(ntime, nant, nchan, ncorr, flag, true));
step1->setNextStep (step2);
step2->setNextStep (step3);
execute (step1);
@@ -194,6 +253,9 @@ int main()
test1(10, 2, 32, 4, false, 1);
test1(10, 5, 32, 4, true, 1);
+ test2( 4, 2, 8, 4, false, 100);
+ test2(10, 5, 32, 4, true, 1);
+ test2( 4, 2, 8, 4, false, 100);
} catch (std::exception& x) {
cout << "Unexpected exception: " << x.what() << endl;
return 1;