diff --git a/.gitattributes b/.gitattributes
index eb7e5409fb7a83c5efb457e379dfbdd4e0771bd6..1248179f7873e9cfa1846b6416dca9e46d2bbfc0 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -226,6 +226,7 @@ CEP/DP3/AOFlagger/include/AOFlagger/msio/system.h -text
CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencydata.h -text
CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencyimager.h -text
CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencymetadata.h -text
+CEP/DP3/AOFlagger/include/AOFlagger/msio/timestepaccessor.h -text
CEP/DP3/AOFlagger/include/AOFlagger/msio/types.h -text
CEP/DP3/AOFlagger/include/AOFlagger/ref/concatenatescript.h -text
CEP/DP3/AOFlagger/include/AOFlagger/ref/copyallscript.h -text
@@ -401,6 +402,7 @@ CEP/DP3/AOFlagger/src/msio/segmentedimage.cpp -text
CEP/DP3/AOFlagger/src/msio/stokesimager.cpp -text
CEP/DP3/AOFlagger/src/msio/timefrequencydata.cpp -text
CEP/DP3/AOFlagger/src/msio/timefrequencyimager.cpp -text
+CEP/DP3/AOFlagger/src/msio/timestepaccessor.cpp -text
CEP/DP3/AOFlagger/src/ns2bbs.cpp -text
CEP/DP3/AOFlagger/src/rfi/antennaflagcountplot.cpp -text
CEP/DP3/AOFlagger/src/rfi/eigenvalue.cpp -text
diff --git a/CEP/DP3/AOFlagger/detectrfi.kdevelop b/CEP/DP3/AOFlagger/detectrfi.kdevelop
index a192725b1bf01e2cdd9af6609ff400fe02aa7b84..768f36ed3d0416d27ee98a8fb8240b772d4b034b 100644
--- a/CEP/DP3/AOFlagger/detectrfi.kdevelop
+++ b/CEP/DP3/AOFlagger/detectrfi.kdevelop
@@ -47,7 +47,7 @@
<default/>
</environments>
<prio>0</prio>
- <defaulttarget>rfigui</defaulttarget>
+ <defaulttarget>aofrequencyfilter</defaulttarget>
<makeoptions>-s</makeoptions>
</make>
<blacklist/>
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/msio/timestepaccessor.h b/CEP/DP3/AOFlagger/include/AOFlagger/msio/timestepaccessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..f81c01a9a0f76cdea1db1231dbf258bd1f0b2ae0
--- /dev/null
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/msio/timestepaccessor.h
@@ -0,0 +1,203 @@
+/***************************************************************************
+ * Copyright (C) 2011 by A.R. Offringa *
+ * offringa@astro.rug.nl *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#ifndef MSIO_TIMESTEP_ACCESSOR_H
+#define MSIO_TIMESTEP_ACCESSOR_H
+
+#include <stdexcept>
+#include <vector>
+
+#include <ms/MeasurementSets/MSTable.h>
+#include <tables/Tables/TableIter.h>
+
+#include <AOFlagger/msio/types.h>
+
+class TimestepAccessorException : public std::runtime_error
+{
+ public: TimestepAccessorException(const std::string &str) : std::runtime_error(str) { }
+};
+
+class TimestepAccessor
+{
+ public:
+ class TimestepIndex
+ {
+ public:
+ TimestepIndex(unsigned tableCount) : tables(new casa::Table*[tableCount]), _tableCount(tableCount)
+ {
+ for(unsigned i=0;i<tableCount;++i)
+ tables[i] = 0;
+ }
+ ~TimestepIndex()
+ {
+ FreeTables();
+ delete[] tables;
+ }
+ void FreeTables()
+ {
+ for(unsigned i=0;i<_tableCount;++i)
+ {
+ if(tables[i] != 0)
+ {
+ delete tables[i];
+ tables[i] = 0;
+ }
+ }
+ }
+
+ casa::Table **tables;
+ private:
+ unsigned _tableCount;
+ TimestepIndex(TimestepIndex &) { }
+ void operator=(TimestepIndex &) { }
+ };
+ struct TimestepData
+ {
+ num_t **realData, **imagData;
+ unsigned antenna1, antenna2;
+ double u,v;
+
+ void Allocate(unsigned polarizationCount, unsigned channelCount)
+ {
+ realData = new num_t*[polarizationCount];
+ imagData = new num_t*[polarizationCount];
+ for(unsigned p=0;p<polarizationCount;++p)
+ {
+ realData[p] = new num_t[channelCount];
+ imagData[p] = new num_t[channelCount];
+ }
+ }
+ void Free(unsigned polarizationCount)
+ {
+ for(unsigned p=0;p<polarizationCount;++p)
+ {
+ delete[] realData[p];
+ delete[] imagData[p];
+ }
+ delete[] realData;
+ delete[] imagData;
+ }
+ };
+
+ TimestepAccessor() : _isOpen(false), _polarizationCount(0), _totalChannelCount(0)
+ {
+ }
+
+ ~TimestepAccessor()
+ {
+ if(_isOpen)
+ Close();
+ }
+
+ void Open();
+
+ void Close();
+
+ void AddMS(const std::string &path)
+ {
+ assertNotOpen();
+ SetInfo newInfo;
+ newInfo.path = path;
+ _sets.push_back(newInfo);
+ }
+
+ unsigned TotalChannelCount() const
+ {
+ assertOpen();
+ return _totalChannelCount;
+ }
+
+ unsigned PolarizationCount() const
+ {
+ assertOpen();
+ return _polarizationCount;
+ }
+
+ double LowestFrequency() const
+ {
+ assertOpen();
+ return _lowestFrequency;
+ }
+
+ double HighestFrequency() const
+ {
+ assertOpen();
+ return _highestFrequency;
+ }
+
+ unsigned TableCount() const
+ {
+ return _sets.size();
+ }
+
+ bool ReadNext(TimestepIndex &index, TimestepData &data);
+
+ void Write(TimestepIndex &index, const TimestepData &data);
+
+ private:
+ struct SetInfo
+ {
+ SetInfo() : tableIter(0)
+ {
+ }
+ SetInfo(const SetInfo &source) :
+ path(source.path),
+ tableIter(source.tableIter),
+ bandCount(source.bandCount),
+ channelsPerBand(source.channelsPerBand),
+ highestFrequency(source.highestFrequency),
+ lowestFrequency(source.lowestFrequency)
+ {
+ }
+ void operator=(const SetInfo &source)
+ {
+ path = source.path;
+ tableIter = source.tableIter;
+ bandCount = source.bandCount;
+ channelsPerBand = source.channelsPerBand;
+ highestFrequency = source.highestFrequency;
+ lowestFrequency = source.lowestFrequency;
+ }
+ std::string path;
+ casa::TableIterator *tableIter;
+ unsigned bandCount, channelsPerBand;
+ double highestFrequency, lowestFrequency;
+ };
+
+ typedef std::vector<SetInfo> SetInfoVector;
+
+ bool _isOpen;
+ unsigned _polarizationCount, _totalChannelCount;
+ SetInfoVector _sets;
+ double _highestFrequency, _lowestFrequency;
+
+ void assertOpen() const
+ {
+ if(!_isOpen)
+ throw TimestepAccessorException("Timestep accessor has not been opened yet");
+ }
+ void assertNotOpen() const
+ {
+ if(_isOpen)
+ throw TimestepAccessorException("Timestep accessor has already been opened");
+ }
+};
+
+#endif
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/ref/concatenatescript.h b/CEP/DP3/AOFlagger/include/AOFlagger/ref/concatenatescript.h
index 441ad196a28576718f3452f9b697941f25260791..41b7846023c671db2d5d1ed04fe238c2bb63eab6 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/ref/concatenatescript.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/ref/concatenatescript.h
@@ -38,12 +38,12 @@ namespace AOTools
stream <<
"#!/usr/bin/env python\n\n"
"import pyrap.tables as pt\n\n"
- "pt.concatms([";
+ "pt.msconcat([";
RefFile::const_iterator i=file.begin();
if(i!=file.end())
{
- stream << '\"' << i->Path() << '\"';
+ stream << "\"/net/" << i->Node() << i->Path() << '\"';
++i;
}
diff --git a/CEP/DP3/AOFlagger/src/CMakeLists.txt b/CEP/DP3/AOFlagger/src/CMakeLists.txt
index 73104a9edd48f1024e6cb8164ffc8cc8759c144d..f95292b391116c9829bb77cfd3d6e796174f759b 100644
--- a/CEP/DP3/AOFlagger/src/CMakeLists.txt
+++ b/CEP/DP3/AOFlagger/src/CMakeLists.txt
@@ -4,7 +4,7 @@ include(LofarPackageVersion)
lofar_add_bin_program(rfihistory rfihistory.cpp)
-lofar_add_bin_program(aofrequencyfilter aofrequencyfilter.cpp rfi/thresholdtools.cpp util/rng.cpp msio/image2d.cpp msio/fitsfile.cpp util/aologger.cpp)
+lofar_add_bin_program(aofrequencyfilter aofrequencyfilter.cpp rfi/thresholdtools.cpp util/rng.cpp msio/image2d.cpp msio/fitsfile.cpp msio/timestepaccessor.cpp util/aologger.cpp)
lofar_add_bin_program(aorefscript aorefscript.cpp)
diff --git a/CEP/DP3/AOFlagger/src/aofrequencyfilter.cpp b/CEP/DP3/AOFlagger/src/aofrequencyfilter.cpp
index 253b15919b4cc820a22b2fa30a7ef97e65e0afd8..5b370e4998132771be1d467d2855ac059699fd54 100644
--- a/CEP/DP3/AOFlagger/src/aofrequencyfilter.cpp
+++ b/CEP/DP3/AOFlagger/src/aofrequencyfilter.cpp
@@ -14,7 +14,9 @@
#include <AOFlagger/rfi/thresholdtools.h>
#include <AOFlagger/imaging/uvimager.h>
+
#include <AOFlagger/msio/system.h>
+#include <AOFlagger/msio/timestepaccessor.h>
using namespace std;
@@ -30,90 +32,50 @@ double uvDist(double u, double v, double firstFrequency, double lastFrequency)
return sqrt(ud * ud + vd * vd) / UVImager::SpeedOfLight();
}
-// This represents data that is constant for all threads
-struct SetInfo
-{
- unsigned bandCount;
- unsigned frequencyCount;
- unsigned polarizationCount;
-};
-
// This represents the data that is specific for a single thread
struct TaskInfo
{
- TaskInfo() : length(0), convolutionSize(0.0), table(0), dataColumn(0), realData(0), imagData(0) { }
+ TaskInfo() : length(0), convolutionSize(0.0), index(0), data() { }
TaskInfo(const TaskInfo &source) :
length(source.length),
convolutionSize(source.convolutionSize),
- table(source.table),
- dataColumn(source.dataColumn),
- realData(source.realData),
- imagData(source.imagData)
+ index(source.index),
+ data(source.data)
{
}
void operator=(const TaskInfo &source)
{
length = source.length;
convolutionSize = source.convolutionSize;
- table = source.table;
- dataColumn = source.dataColumn;
- realData = source.realData;
- imagData = source.imagData;
+ index = source.index;
+ data = source.data;
}
unsigned length;
double convolutionSize;
- casa::Table *table;
- casa::ArrayColumn<casa::Complex> *dataColumn;
- float
- **realData,
- **imagData;
+ TimestepAccessor::TimestepIndex *index;
+ TimestepAccessor::TimestepData data;
};
-void performAndWriteConvolution(const SetInfo &set, TaskInfo task, boost::mutex &mutex)
+void performAndWriteConvolution(TimestepAccessor &accessor, TaskInfo task, boost::mutex &mutex)
{
+ unsigned polarizationCount = accessor.PolarizationCount();
if(task.convolutionSize > 1.0)
{
// Convolve the data
- for(unsigned p=0;p<set.polarizationCount;++p)
+ for(unsigned p=0;p<polarizationCount;++p)
{
- ThresholdTools::OneDimensionalSincConvolution(task.realData[p], task.length, task.convolutionSize);
- ThresholdTools::OneDimensionalSincConvolution(task.imagData[p], task.length, task.convolutionSize);
+ ThresholdTools::OneDimensionalSincConvolution(task.data.realData[p], task.length, task.convolutionSize);
+ ThresholdTools::OneDimensionalSincConvolution(task.data.imagData[p], task.length, task.convolutionSize);
}
- boost::mutex::scoped_lock lock(mutex);
// Copy data back to tables
- for(unsigned i=0;i<set.bandCount;++i)
- {
- casa::Array<casa::Complex> dataArray = (*task.dataColumn)(i);
- casa::Array<casa::Complex>::iterator dataIterator = dataArray.begin();
- unsigned index = i * set.frequencyCount;
- for(unsigned f=0;f<set.frequencyCount;++f)
- {
- for(unsigned p=0;p<set.polarizationCount;++p)
- {
- *dataIterator = casa::Complex(task.realData[p][index], task.imagData[p][index]);
- ++dataIterator;
- }
- ++index;
- }
- task.dataColumn->basePut(i, dataArray);
- }
+ boost::mutex::scoped_lock lock(mutex);
+ accessor.Write(*task.index, task.data);
lock.unlock();
}
- // Free memory
- boost::mutex::scoped_lock lock(mutex);
- for(unsigned p=0;p<set.polarizationCount;++p)
- {
- delete[] task.realData[p];
- delete[] task.imagData[p];
- }
- delete[] task.realData;
- delete[] task.imagData;
-
- delete task.dataColumn;
- delete task.table;
+ task.data.Free(polarizationCount);
}
struct ThreadFunction
@@ -126,8 +88,8 @@ struct ThreadFunction
class ThreadControl
{
public:
- ThreadControl(unsigned threadCount, const struct SetInfo &setInfo)
- : _setInfo(setInfo), _threadCount(threadCount), _isFinishing(false)
+ ThreadControl(unsigned threadCount, TimestepAccessor &accessor)
+ : _accessor(accessor), _threadCount(threadCount), _isFinishing(false)
{
for(unsigned i=0;i<threadCount;++i)
{
@@ -173,10 +135,10 @@ class ThreadControl
_dataAvailableCondition.notify_all();
_threadGroup.join_all();
}
- const struct SetInfo &SetInfo() const { return _setInfo; }
boost::mutex &WriteMutex() { return _writeMutex; }
+ TimestepAccessor &Accessor() { return _accessor; }
private:
- const struct SetInfo _setInfo;
+ TimestepAccessor &_accessor;
boost::thread_group _threadGroup;
unsigned _threadCount;
bool _isFinishing;
@@ -194,7 +156,7 @@ void ThreadFunction::operator()()
bool hasTask = threadControl->WaitForTask(task);
while(hasTask)
{
- performAndWriteConvolution(threadControl->SetInfo(), task, threadControl->WriteMutex());
+ performAndWriteConvolution(threadControl->Accessor(), task, threadControl->WriteMutex());
hasTask = threadControl->WaitForTask(task);
}
cout << "Thread " << number << " finished\n";
@@ -202,169 +164,81 @@ void ThreadFunction::operator()()
int main(int argc, char *argv[])
{
- if(argc != 3)
+ if(argc < 3)
{
- cerr << "Syntax: " << argv[0] << " <fringe size> <MS>\n";
+ cerr << "Syntax: " << argv[0] << " <fringe size> <MS1> [<MS2> [..]]\n";
cerr << " fringe size is a double and should be given in units of wavelength / fringe.\n";
} else {
const double fringeSize = atof(argv[1]);
- const string msFilename = argv[2];
cout << "Fringe size: " << fringeSize << '\n';
- casa::MeasurementSet ms(msFilename);
- casa::Table mainTable(msFilename, casa::Table::Update);
- class SetInfo setInfo;
-
- //count number of polarizations
- casa::Table polTable = ms.polarization();
- casa::ROArrayColumn<int> corTypeColumn(polTable, "CORR_TYPE");
- setInfo.polarizationCount = corTypeColumn(0).shape()[0];
- cout << "Number of polarizations: " << setInfo.polarizationCount << '\n';
+ TimestepAccessor accessor;
+ for(int i=2;i<argc;++i)
+ accessor.AddMS(argv[i]);
+ accessor.Open();
- // Find lowest and highest frequency and check order
- double lowestFrequency = 0.0, highestFrequency = 0.0;
- casa::Table spectralWindowTable = ms.spectralWindow();
- casa::ROArrayColumn<double> frequencyCol(spectralWindowTable, "CHAN_FREQ");
- for(unsigned b=0;b<spectralWindowTable.nrow();++b)
- {
- casa::Array<double> frequencyArray = frequencyCol(b);
- casa::Array<double>::const_iterator frequencyIterator = frequencyArray.begin();
- while(frequencyIterator != frequencyArray.end())
- {
- double frequency = *frequencyIterator;
- if(lowestFrequency == 0.0) lowestFrequency = frequency;
- if(frequency < lowestFrequency || frequency <= highestFrequency)
- {
- cerr << "ERROR: Channels are not ordered in increasing frequency!\n";
- abort();
- }
- highestFrequency = frequency;
- ++frequencyIterator;
- }
- }
- setInfo.bandCount = spectralWindowTable.nrow();
+ cout << "Number of polarizations: " << accessor.PolarizationCount() << '\n';
cout
- << "Number of bands: " << setInfo.bandCount
- << " (" << round(lowestFrequency/1e6) << " MHz - "
- << round(highestFrequency/1e6) << " MHz)\n";
+ << "Number of channels: " << accessor.TotalChannelCount()
+ << " (" << round(accessor.LowestFrequency()/1e6) << " MHz - "
+ << round(accessor.HighestFrequency()/1e6) << " MHz)\n";
- setInfo.frequencyCount =
- casa::ROArrayColumn<casa::Complex>(mainTable, "DATA")(0).shape()[1];
- cout << "Channels per band: " << setInfo.frequencyCount << '\n';
-
- const unsigned long totalIterations =
- mainTable.nrow() / spectralWindowTable.nrow();
+ const unsigned long totalIterations = 1; // TODO
cout << "Total iterations: " << totalIterations << '\n';
- unsigned processorCount = System::ProcessorCount();
+ const unsigned processorCount = System::ProcessorCount();
cout << "CPUs: " << processorCount << '\n';
- ThreadControl threads(processorCount, setInfo);
+ ThreadControl threads(processorCount, accessor);
- // Create the sorted table and iterate over it
- casa::Block<casa::String> names(4);
- names[0] = "TIME";
- names[1] = "ANTENNA1";
- names[2] = "ANTENNA2";
- names[3] = "DATA_DESC_ID";
- cout << "Sorting...\n";
- casa::Table sortab = mainTable.sort(names);
- cout << "Iterating...\n";
- unsigned long iterSteps = 0;
- names.resize(3, true, true);
- casa::TableIterator iter (sortab, names, casa::TableIterator::Ascending, casa::TableIterator::NoSort);
double maxFringeChannels = 0.0, minFringeChannels = 1e100;
- while (! iter.pastEnd()) {
+ const unsigned totalChannels = accessor.TotalChannelCount();
+
+ TimestepAccessor::TimestepIndex *index = new TimestepAccessor::TimestepIndex(accessor.TableCount());
+ TimestepAccessor::TimestepData data;
+
+ data.Allocate(accessor.PolarizationCount(), totalChannels);
+
+ unsigned iterSteps = 0;
+
+ boost::mutex::scoped_lock lock(threads.WriteMutex());
+ while(accessor.ReadNext(*index, data)) {
+ lock.unlock();
+
TaskInfo task;
- task.table = new casa::Table(iter.table());
-
- int antenna1, antenna2;
- // we start a new block to let the columns go out of scope (table might be freed
- // before end of parent block)
- {
- casa::ROScalarColumn<int> antenna1Column =
- casa::ROScalarColumn<int>(*task.table, "ANTENNA1");
- casa::ROScalarColumn<int> antenna2Column =
- casa::ROScalarColumn<int>(*task.table, "ANTENNA2");
- antenna1 = antenna1Column(0);
- antenna2 = antenna2Column(0);
- }
+ task.index = index;
+ task.data = data;
// Skip autocorrelations
- if(antenna1 == antenna2)
+ if(data.antenna1 == data.antenna2)
{
- delete task.table;
+ data.Free(accessor.PolarizationCount());
+ delete index;
} else
{
- task.dataColumn = new casa::ArrayColumn<casa::Complex>(*task.table, "DATA");
-
- // Check number of channels & bands
- const casa::IPosition &dataShape = task.dataColumn->shape(0);
- if(dataShape[1] != setInfo.frequencyCount) {
- std::cerr << "ERROR: bands do not have equal number of channels!\n";
- abort();
- }
- if(task.table->nrow() != setInfo.bandCount) {
- std::cerr << "ERROR: inconsistent band information in specific correlation/time step\n"
- " (rows in table iterator's table: " << task.table->nrow() <<
- ", in set: " << setInfo.bandCount << ")\n";
- abort();
- }
-
- // Retrieve uv info and calculate the convolution size
- double u, v;
- {
- casa::ROArrayColumn<double> uvwColumn =
- casa::ROArrayColumn<double>(*task.table, "UVW");
- casa::Array<double> uvwArray = uvwColumn(0);
- casa::Array<double>::const_iterator uvwIterator = uvwArray.begin();
- u = *uvwIterator;
- ++uvwIterator;
- v = *uvwIterator;
- }
- task.length = setInfo.bandCount * setInfo.frequencyCount;
- task.convolutionSize = fringeSize * (double) task.length / uvDist(u, v, lowestFrequency, highestFrequency);
+ // Calculate the convolution size
+ double u = data.u, v = data.v;
+ task.length = totalChannels;
+ task.convolutionSize = fringeSize * (double) totalChannels / uvDist(u, v, accessor.LowestFrequency(), accessor.HighestFrequency());
if(task.convolutionSize > maxFringeChannels) maxFringeChannels = task.convolutionSize;
if(task.convolutionSize < minFringeChannels) minFringeChannels = task.convolutionSize;
+ task.data = data;
- // Allocate memory for putting the data of all channels in an array, for each polarization
- task.realData = new float*[setInfo.polarizationCount];
- task.imagData = new float*[setInfo.polarizationCount];
- for(unsigned p=0;p<setInfo.polarizationCount;++p)
- {
- task.realData[p] = new float[task.length];
- task.imagData[p] = new float[task.length];
- }
-
- boost::mutex::scoped_lock lock(threads.WriteMutex());
- // Copy data from tables in arrays
- for(unsigned i=0;i<setInfo.bandCount;++i)
- {
- casa::Array<casa::Complex> dataArray = (*task.dataColumn)(i);
- casa::Array<casa::Complex>::const_iterator dataIterator = dataArray.begin();
- unsigned index = i * setInfo.frequencyCount;
- for(unsigned f=0;f<setInfo.frequencyCount;++f)
- {
- for(unsigned p=0;p<setInfo.polarizationCount;++p)
- {
- task.realData[p][index] = (*dataIterator).real();
- task.imagData[p][index] = (*dataIterator).imag();
- ++dataIterator;
- }
- ++index;
- }
- }
- lock.unlock();
-
+ // Add task
threads.PushTask(task);
}
- iter.next();
+ index = new TimestepAccessor::TimestepIndex(accessor.TableCount());
+ data.Allocate(accessor.PolarizationCount(), totalChannels);
+
+ lock.lock();
++iterSteps;
- if((iterSteps * 100UL) % totalIterations < 100UL)
- cout << '.' << flush;
- if((iterSteps * 10UL) % totalIterations < 10UL)
- cout << (iterSteps*100/totalIterations) << '%' << flush;
+
+ if(iterSteps%100==0) cout << '.' << flush;
}
+ lock.unlock();
+
+ data.Free(accessor.PolarizationCount());
+
cout << '\n';
threads.Finish();
cout
diff --git a/CEP/DP3/AOFlagger/src/msio/timestepaccessor.cpp b/CEP/DP3/AOFlagger/src/msio/timestepaccessor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bceb4f370eba9fd880205a8924f1bc968e5898ba
--- /dev/null
+++ b/CEP/DP3/AOFlagger/src/msio/timestepaccessor.cpp
@@ -0,0 +1,177 @@
+#include <AOFlagger/msio/timestepaccessor.h>
+
+#include <ms/MeasurementSets/MeasurementSet.h>
+#include <ms/MeasurementSets/MSColumns.h>
+
+void TimestepAccessor::Open()
+{
+ assertNotOpen();
+
+ _totalChannelCount = 0;
+ _highestFrequency = 0.0;
+ _lowestFrequency = 0.0;
+
+ for(SetInfoVector::iterator i=_sets.begin(); i!=_sets.end(); ++i)
+ {
+ SetInfo &set = *i;
+
+ casa::MeasurementSet ms(set.path);
+ casa::Table mainTable(set.path, casa::Table::Update);
+
+ // Create the sorted table and iterator
+ casa::Block<casa::String> names(4);
+ names[0] = "TIME";
+ names[1] = "ANTENNA1";
+ names[2] = "ANTENNA2";
+ names[3] = "DATA_DESC_ID";
+ casa::Table sortab = mainTable.sort(names);
+ names.resize(3, true, true);
+ set.tableIter = new casa::TableIterator(sortab, names, casa::TableIterator::Ascending, casa::TableIterator::NoSort);
+
+ // Check number of polarizations
+ casa::Table polTable = ms.polarization();
+ casa::ROArrayColumn<int> corTypeColumn(polTable, "CORR_TYPE");
+ if(_polarizationCount==0 && i==_sets.begin())
+ _polarizationCount = corTypeColumn(0).shape()[0];
+ else if(_polarizationCount != corTypeColumn(0).shape()[0])
+ throw TimestepAccessorException("Number of polarizations don't match!");
+
+ // Find lowest and highest frequency and check order
+ set.lowestFrequency = 0.0;
+ set.highestFrequency = 0.0;
+ casa::Table spectralWindowTable = ms.spectralWindow();
+ casa::ROArrayColumn<double> frequencyCol(spectralWindowTable, "CHAN_FREQ");
+ for(unsigned b=0;b<spectralWindowTable.nrow();++b)
+ {
+ casa::Array<double> frequencyArray = frequencyCol(b);
+ casa::Array<double>::const_iterator frequencyIterator = frequencyArray.begin();
+ while(frequencyIterator != frequencyArray.end())
+ {
+ double frequency = *frequencyIterator;
+ if(set.lowestFrequency == 0.0) set.lowestFrequency = frequency;
+ if(frequency < set.lowestFrequency || frequency <= set.highestFrequency)
+ throw TimestepAccessorException("Channels or bands are not ordered in increasing frequency");
+ set.highestFrequency = frequency;
+ ++frequencyIterator;
+ }
+ }
+ if(set.lowestFrequency < _highestFrequency)
+ throw TimestepAccessorException("Sub-bands are not given in order of increasing frequency");
+ if(_lowestFrequency == 0.0) _lowestFrequency = set.lowestFrequency;
+ _highestFrequency = set.highestFrequency;
+
+ set.bandCount = spectralWindowTable.nrow();
+ set.channelsPerBand = casa::ROArrayColumn<casa::Complex>(mainTable, "DATA")(0).shape()[1];
+ _totalChannelCount += set.bandCount * set.channelsPerBand;
+ }
+ _isOpen = true;
+}
+
+void TimestepAccessor::Close()
+{
+ assertOpen();
+
+ for(SetInfoVector::iterator i=_sets.begin(); i!=_sets.end(); ++i)
+ {
+ delete i->tableIter;
+ }
+ _isOpen = false;
+}
+
+bool TimestepAccessor::ReadNext(TimestepAccessor::TimestepIndex &index, TimestepAccessor::TimestepData &data)
+{
+ assertOpen();
+
+ index.FreeTables();
+
+ double timeStep = 0.0;
+ unsigned valIndex = 0;
+ casa::Table **tablePtr = index.tables;
+
+ for(SetInfoVector::iterator i=_sets.begin(); i!=_sets.end(); ++i)
+ {
+ SetInfo &set = *i;
+
+ if(set.tableIter->pastEnd())
+ return false;
+ casa::Table table(set.tableIter->table());
+ *tablePtr = new casa::Table(table);
+
+ // Check timestep & read u,v coordinates & antenna's
+ casa::ROScalarColumn<double> timeColumn = casa::ROScalarColumn<double>(table, "TIME");
+ casa::ROArrayColumn<double> uvwColumn = casa::ROArrayColumn<double>(table, "UVW");
+ if(timeStep == 0.0) {
+ timeStep = timeColumn(0);
+ casa::Array<double> uvwArray = uvwColumn(0);
+ casa::Array<double>::const_iterator uvwIterator = uvwArray.begin();
+ data.u = *uvwIterator;
+ ++uvwIterator;
+ data.v = *uvwIterator;
+ casa::ROScalarColumn<int>
+ antenna1Column = casa::ROScalarColumn<int>(table, "ANTENNA1"),
+ antenna2Column = casa::ROScalarColumn<int>(table, "ANTENNA2");
+ data.antenna1 = antenna1Column(0);
+ data.antenna2 = antenna2Column(0);
+ }
+ else if(timeStep != timeColumn(0))
+ throw TimestepAccessorException("Sets do not have same time steps");
+
+ // Copy data from tables in arrays
+ casa::ROArrayColumn<casa::Complex> dataColumn(table, "DATA");
+ for(unsigned band=0;band<set.bandCount;++band)
+ {
+ casa::Array<casa::Complex> dataArray = dataColumn(band);
+ casa::Array<casa::Complex>::const_iterator dataIterator = dataArray.begin();
+ for(unsigned f=0;f<set.channelsPerBand;++f)
+ {
+ for(unsigned p=0;p<_polarizationCount;++p)
+ {
+ data.realData[p][valIndex] = (*dataIterator).real();
+ data.imagData[p][valIndex] = (*dataIterator).imag();
+ ++dataIterator;
+ }
+ ++valIndex;
+ }
+ }
+
+ set.tableIter->next();
+ ++tablePtr;
+ }
+ return true;
+}
+
+void TimestepAccessor::Write(TimestepAccessor::TimestepIndex &index, const TimestepAccessor::TimestepData &data)
+{
+ assertOpen();
+
+ casa::Table **tablePtr = index.tables;
+ unsigned valIndex = 0;
+
+ for(SetInfoVector::iterator i=_sets.begin(); i!=_sets.end(); ++i)
+ {
+ const SetInfo &set = *i;
+
+ casa::Table &table = **tablePtr;
+
+ // Copy data from arrays in tables
+ casa::ArrayColumn<casa::Complex> dataColumn(table, "DATA");
+ for(unsigned band=0;band<set.bandCount;++band)
+ {
+ casa::Array<casa::Complex> dataArray = dataColumn(band);
+ casa::Array<casa::Complex>::iterator dataIterator = dataArray.begin();
+ for(unsigned f=0;f<set.channelsPerBand;++f)
+ {
+ for(unsigned p=0;p<_polarizationCount;++p)
+ {
+ (*dataIterator).real() = data.realData[p][valIndex];
+ (*dataIterator).imag() = data.imagData[p][valIndex];
+ ++dataIterator;
+ }
+ ++valIndex;
+ }
+ dataColumn.basePut(band, dataArray);
+ }
+ set.tableIter->next();
+ }
+ index.FreeTables();
+}