Bug #1113: initial working version

.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

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],

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  <<<===============

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;

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)

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[])
 {

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));
   }

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);
 }

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);

Appl/CEP/CS1/CS1_pp_lib/test/CS1_IDPPP.log_prop

+#
+# 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
+

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
+squasher=1
\ No newline at end of file