From c6b6f73cd20049b23f2ddc43eb93a3bcdb1cd4f4 Mon Sep 17 00:00:00 2001
From: Adriaan Renting <renting@astron.nl>
Date: Tue, 24 Jun 2008 23:22:15 +0000
Subject: [PATCH] Bug #1113: initial working version
---
.gitattributes | 1 +
.../CS1_pp_lib/src/ComplexMedianFlagger.cc | 7 ++--
Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.cc | 9 +++++
Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.h | 1 -
Appl/CEP/CS1/CS1_pp_lib/src/DataSquasher.cc | 22 ++++++++++---
Appl/CEP/CS1/CS1_pp_lib/src/IDPPP.cc | 1 +
Appl/CEP/CS1/CS1_pp_lib/src/MsFile.cc | 6 ++--
Appl/CEP/CS1/CS1_pp_lib/src/Pipeline.cc | 28 ++++++++++------
.../CS1_pp_lib/src/PipelineProcessControl.cc | 20 ++++++-----
.../CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop | 33 +++++++++++++++++++
Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.parset | 6 ++--
11 files changed, 100 insertions(+), 34 deletions(-)
create mode 100644 Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop
diff --git a/.gitattributes b/.gitattributes
index ad1858c8021..26ce8460525 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -123,6 +123,7 @@ Appl/CEP/CS1/CS1_pp_lib/src/PipelineProcessControl.cc -text
Appl/CEP/CS1/CS1_pp_lib/src/RunDetails.cc -text
Appl/CEP/CS1/CS1_pp_lib/src/RunDetails.h -text
Appl/CEP/CS1/CS1_pp_lib/src/versioncs1_idppp.cc -text
+Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop -text
Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.parset -text
Appl/CEP/CS1/CS1_pp_lib/test/Makefile.am -text
Appl/CMakeLists.txt -text
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/ComplexMedianFlagger.cc b/Appl/CEP/CS1/CS1_pp_lib/src/ComplexMedianFlagger.cc
index e5c58cc2c17..3c93d89a1b8 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/ComplexMedianFlagger.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/ComplexMedianFlagger.cc
@@ -100,7 +100,8 @@ void ComplexMedianFlagger::ProcessTimeslot(DataBuffer& data,
FlaggerStatistics& stats)
{
//Data.Position is the last filled timeslot, the middle is 1/2 a window behind it
- int pos = (data.Position - (data.WindowSize-1)/2) % data.WindowSize;
+ int pos = (data.Position + (data.WindowSize+1)/2) % data.WindowSize;
+ cout << pos << endl;
NumChannels = info.NumChannels;
NumPolarizations = info.NumPolarizations;
WindowSize = data.WindowSize;
@@ -115,8 +116,8 @@ void ComplexMedianFlagger::ProcessTimeslot(DataBuffer& data,
* info.BaselineLengths[info.BaselineIndex[baseline_t(j, k)]]
/ info.MaxBaselineLength
) * info.NoiseLevel;
- Matrix<Bool> flags;
- data.Flags[index].xyPlane(pos).reference(flags);
+ Matrix<Bool> flags = data.Flags[index].xyPlane(pos);
+ //data.Flags[index].xyPlane(pos).reference(flags);
// if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
stats(i, j, k) = FlagBaselineBand(flags,
data.Data[index],
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.cc b/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.cc
index 826037a67d7..d56e8dc6b85 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.cc
@@ -134,6 +134,15 @@ void DataBuffer::Init(bool Columns)
void DataBuffer::PrintInfo(void)
{
+ std::cout << "Baselines: " << Data.size() << std::endl;
+ std::cout << "Data: " << Data[0].shape() << std::endl;
+ std::cout << "Flags: " << Flags[0].shape() << std::endl;
+ std::cout << "Weights: " << Weights[0].shape() << std::endl;
+ if (ModelData.size())
+ {
+ std::cout << "ModelData: " << ModelData[0].shape() << std::endl;
+ std::cout << "CorrectedData: " << CorrectedData[0].shape() << std::endl;
+ }
}
//===============>>> DataBuffer <<<===============
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.h b/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.h
index 84ebfa24f40..cc42c443335 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.h
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/DataBuffer.h
@@ -45,7 +45,6 @@ namespace LOFAR
std::vector< casa::Cube<casa::Float> > Weights;
void DeterminePolarizationsToCheck(bool UseOnlyXpolarizations);
void PrintInfo(void);
- void UpdateData();
private:
MsInfo* myInfo;
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/DataSquasher.cc b/Appl/CEP/CS1/CS1_pp_lib/src/DataSquasher.cc
index df5d2d9f82d..85e8aeb8b3b 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/DataSquasher.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/DataSquasher.cc
@@ -64,12 +64,11 @@ void DataSquasher::Squash(Matrix<Complex>& oldData, Matrix<Complex>& newData,
incounter++;
if ((incounter) % Step == 0)
{
- if (flagnew)
+/* if (flagnew)
{
newWeights.column(outcounter) = 1.0;
for (int i = 0; i < itsNumPolarizations; i++)
newData(i, outcounter) = values(i)/Step;
- newFlags.column(outcounter) = true;
}
else
{
@@ -79,7 +78,13 @@ void DataSquasher::Squash(Matrix<Complex>& oldData, Matrix<Complex>& newData,
}
newData.column(outcounter) = values;
newFlags.column(outcounter) = flagnew;
+ }*/
+ for (int i = 0; i < itsNumPolarizations; i++)
+ { values(i) = values(i) / weights(i);
+ newWeights(i, outcounter) = abs(weights(i)) / Step;
}
+ newData.column(outcounter) = values;
+ newFlags.column(outcounter) = flagnew;
values = 0;
weights = 0;
outcounter++;
@@ -94,12 +99,14 @@ void DataSquasher::ProcessTimeslot(DataBuffer& InData, DataBuffer& OutData,
MsInfo& Info, RunDetails& Details)
{
//Data.Position is the last filled timeslot, the middle is 1/2 a window behind it
- int pos = (InData.Position - (InData.WindowSize-1)/2) % InData.WindowSize;
+ int pos = (InData.Position + (InData.WindowSize+1)/2) % InData.WindowSize;
Matrix<Complex> myOldData;
Matrix<Complex> myNewData;
Matrix<Bool> myOldFlags;
Matrix<Bool> myNewFlags;
Matrix<Float> NewWeights;
+ //InData.PrintInfo();
+ //OutData.PrintInfo();
for (int i = 0; i < Info.NumBands; i++)
{
for(int j = 0; j < Info.NumAntennae; j++)
@@ -108,11 +115,16 @@ void DataSquasher::ProcessTimeslot(DataBuffer& InData, DataBuffer& OutData,
{
int index = i * Info.NumPairs + Info.BaselineIndex[baseline_t(j, k)];
- InData.Data[index].xyPlane(pos).reference(myOldData);
+/* InData.Data[index].xyPlane(pos).reference(myOldData);
OutData.Data[index].xyPlane(pos).reference(myNewData);
InData.Flags[index].xyPlane(pos).reference(myOldFlags);
OutData.Flags[index].xyPlane(pos).reference(myNewFlags);
- OutData.Weights[index].xyPlane(pos).reference(NewWeights);
+ OutData.Weights[index].xyPlane(pos).reference(NewWeights);*/
+ myOldData = InData.Data[index].xyPlane(pos);
+ myNewData = OutData.Data[index].xyPlane(OutData.Position);
+ myOldFlags = InData.Flags[index].xyPlane(pos);
+ myNewFlags = OutData.Flags[index].xyPlane(OutData.Position);
+ NewWeights = OutData.Weights[index].xyPlane(OutData.Position);
Squash(myOldData, myNewData, myOldFlags, myNewFlags, NewWeights,
Info.NumPolarizations, Details.Start, Details.Step, Details.NChan);
if (Details.Columns)
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/IDPPP.cc b/Appl/CEP/CS1/CS1_pp_lib/src/IDPPP.cc
index a1280bdbb42..7575f44b161 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/IDPPP.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/IDPPP.cc
@@ -26,6 +26,7 @@
#include <PLC/ACCmain.h>
#include <casa/Exceptions.h>
#include <CS1_IDPPP/PipelineProcessControl.h>
+#include <Common/LofarLogger.h>
int main(int argc, char *argv[])
{
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/MsFile.cc b/Appl/CEP/CS1/CS1_pp_lib/src/MsFile.cc
index 6cc6f973996..60a4d5cfece 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/MsFile.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/MsFile.cc
@@ -87,7 +87,7 @@ void MsFile::Init(MsInfo& Info, RunDetails& Details)
std::cout << "New shape: " << temp_pos(0) << ":" << temp_pos(1) << std::endl;
IPosition data_ipos(temp_pos);
- tdesc.removeColumn("WEIGHT_SPECTRUM");
+ //tdesc.removeColumn("WEIGHT_SPECTRUM");
tdesc.addColumn(ArrayColumnDesc<Float>("WEIGHT_SPECTRUM", "Added by datasquasher",
data_ipos, ColumnDesc::FixedShape));
@@ -258,7 +258,8 @@ void MsFile::UpdateTimeslotData(casa::TableIterator Data_iter,
data.getColumn(tempData); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
flags.getColumn(tempFlags);
- Buffer.Position = Buffer.Position++ % Buffer.WindowSize;
+ Buffer.Position = ++(Buffer.Position) % Buffer.WindowSize;
+ cout << Buffer.Position << endl;
for (int i = 0; i < rowcount; i++)
{
int bi = Info.BaselineIndex[baseline_t(antenna1(i), antenna2(i))];
@@ -289,7 +290,6 @@ void MsFile::WriteData(casa::TableIterator Data_iter,
int bi = Info.BaselineIndex[baseline_t(antenna1(i), antenna2(i))];
int band = bandnr(i);
int index = (band % Info.NumBands) * Info.NumPairs + bi;
-
data.put(i, Buffer.Data[index].xyPlane(Buffer.Position));
flags.put(i, Buffer.Flags[index].xyPlane(Buffer.Position));
}
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/Pipeline.cc b/Appl/CEP/CS1/CS1_pp_lib/src/Pipeline.cc
index f0072a3eb64..1603371da59 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/Pipeline.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/Pipeline.cc
@@ -72,6 +72,7 @@ void Pipeline::Run(MsInfo* SquashedInfo, bool Columns)
// else
// { FlaggerData = BandpassData;
// }
+ FlaggerData = BandpassData;
if (mySquasher)
{ SquasherData = new DataBuffer(SquashedInfo, myDetails->TimeWindow, Columns);
}
@@ -82,33 +83,40 @@ void Pipeline::Run(MsInfo* SquashedInfo, bool Columns)
TableIterator read_iter = (*myFile).ReadIterator();
TableIterator write_iter = (*myFile).WriteIterator();
int TimeCounter = 0;
- int step = myInfo->NumTimeslots / 10 + 1; //not exact but it'll do
- int row = 0;
+ int step = myInfo->NumTimeslots / 10 + 1; //not exact but it'll do
+ int row = 0;
while (!read_iter.pastEnd())
- {
+ { cout << "reading data" << endl;
myFile->UpdateTimeslotData(read_iter, *myInfo, *BandpassData);
read_iter++;
if (myBandpass)
- { myBandpass->ProcessTimeslot(*BandpassData, *myInfo, *myDetails);
+ { cout << "Running Bandpass correction" << endl;
+ myBandpass->ProcessTimeslot(*BandpassData, *myInfo, *myDetails);
}
- if (TimeCounter > (BandpassData->WindowSize - 1)/2)
+ if (TimeCounter >= (BandpassData->WindowSize - 1)/2)
{
if (myFlagger)
- { myFlagger->ProcessTimeslot(*BandpassData, *myInfo, *myDetails, *myStatistics);
+ { cout << "Running flagger" << endl;
+ myFlagger->ProcessTimeslot(*BandpassData, *myInfo, *myDetails, *myStatistics);
}
if (mySquasher)
- { mySquasher->ProcessTimeslot(*BandpassData, *SquasherData, *myInfo, *myDetails);
+ { cout << "Running Data reduction" << endl;
+ SquasherData->Position = (SquasherData->Position + 1) % SquasherData->WindowSize;
+ cout << SquasherData->Position << endl;
+ mySquasher->ProcessTimeslot(*BandpassData, *SquasherData, *myInfo, *myDetails);
}
+ cout << "writing data" << endl;
myFile->WriteData(write_iter, *myInfo, *SquasherData);
write_iter++;
- if (row++ % step == 0) // to tell the user how much % we have processed,
- { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
- }
}
else
{
initBuffer(*BandpassData, *myInfo);
}
+ TimeCounter++;
+ if (row++ % step == 0) // to tell the user how much % we have processed,
+ { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
+ }
}
myStatistics->PrintStatistics(std::cout);
}
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/PipelineProcessControl.cc b/Appl/CEP/CS1/CS1_pp_lib/src/PipelineProcessControl.cc
index 293d92fc34b..dec1d87bfc9 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/src/PipelineProcessControl.cc
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/PipelineProcessControl.cc
@@ -45,10 +45,10 @@ namespace LOFAR
: ProcessControl()
{
myPipeline = NULL;
- myFile = NULL;
- myInfo = NULL;
+ myFile = NULL;
+ myInfo = NULL;
myBandpass = NULL;
- myFlagger = NULL;
+ myFlagger = NULL;
mySquasher = NULL;
}
@@ -94,6 +94,8 @@ namespace LOFAR
std::cout << "Runnning pipeline please wait..." << std::endl;
myFile->Init(*myInfo, *myDetails);
MsInfo* outInfo = new MsInfo(itsOutMS);
+ outInfo->Update();
+ outInfo->PrintInfo();
myPipeline->Run(outInfo, myDetails->Columns);
delete outInfo;
}
@@ -123,18 +125,18 @@ namespace LOFAR
myInfo->PrintInfo();
switch (itsBandpass)
{
- case 1: myBandpass = new BandpassCorrector();
+ case 1: myBandpass = new BandpassCorrector(); break;
}
switch (itsFlagger)
{
- case 1: myFlagger = new ComplexMedianFlagger();
- // case 2: myFlagger = new ComplexMedian2Flagger();
- // case 3: myFlagger = new FrequencyFlagger();
- case 4: myFlagger = new MADFlagger();
+ case 1: myFlagger = new ComplexMedianFlagger(); break;
+ // case 2: myFlagger = new ComplexMedian2Flagger(); break;
+ // case 3: myFlagger = new FrequencyFlagger(); break;
+ case 4: myFlagger = new MADFlagger(); break;
}
switch (itsSquasher)
{
- case 1: mySquasher = new DataSquasher();
+ case 1: mySquasher = new DataSquasher(); break;
}
myPipeline = new Pipeline(myInfo, myFile, myDetails,
myBandpass, myFlagger, mySquasher);
diff --git a/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop b/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop
new file mode 100644
index 00000000000..ff008b5b7d7
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop
@@ -0,0 +1,33 @@
+#
+# setup the right levels for logging and tracing
+#
+log4cplus.rootLogger=DEBUG, STDOUT
+log4cplus.logger.TRC=DEBUG, STDOUT
+log4cplus.additivity.TRC=FALSE
+
+#
+# define the output channels
+#
+log4cplus.appender.STDOUT=log4cplus::ConsoleAppender
+log4cplus.appender.STDOUT.layout=log4cplus::PatternLayout
+log4cplus.appender.STDOUT.layout.ConversionPattern=%D{%d-%m-%y %H:%M:%S} %-5p %c{3} - %m [%.25l]%n
+
+#log4cplus.appender.MACSTDERR=log4cplus::ConsoleAppender
+#log4cplus.appender.MACSTDERR.layout=log4cplus::PatternLayout
+#log4cplus.appender.MACSTDERR.layout.ConversionPattern=%x %D{%d-%m %H:%M:%S.%q} %-5p %c{4} - %m [%.25l]%n
+#log4cplus.appender.MACSTDERR.logToStdErr=true
+
+#log4cplus.appender.MACCLP=log4cplus::SocketAppender
+#log4cplus.appender.MACCLP.port=23999
+#log4cplus.appender.MACCLP.host=mcu001t
+#log4cplus.appender.MACCLP.Threshold=INFO
+
+#log4cplus.appender.FILE=log4cplus::RollingFileAppender
+#log4cplus.appender.FILE.File=./test.log
+#log4cplus.appender.FILE.MaxFileSize=10MB
+#log4cplus.appender.FILE.MaxBackupIndex=9
+#log4cplus.appender.FILE.layout=log4cplus::PatternLayout
+#log4cplus.appender.FILE.layout.ConversionPattern=%x %D{%d-%m-%y %H:%M:%S} %-5p %c{3} - %m [%.25l]%n
+
+#log4cplus.appender.DUMP=log4cplus::NullAppender
+
diff --git a/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.parset b/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.parset
index 2c2de0a302f..0cbd3a1a5d1 100644
--- a/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.parset
+++ b/Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.parset
@@ -4,8 +4,8 @@ timewindow=7
treshold=1.0
min=10
max=12
-existing=True
-nchan=16
+existing=False
+nchan=128
start=32
step=16
skipflags=True
@@ -14,4 +14,4 @@ msin=test.MS
msout=test_out.MS
bandpass=1
flagger=1
-squasher=1
\ No newline at end of file
+squasher=1
--
GitLab