From 72a72afc21363f8284c56d3cdcaebedd0ade7ff9 Mon Sep 17 00:00:00 2001
From: Joris van Zwieten <zwieten@astron.nl>
Date: Tue, 28 Aug 2007 09:25:23 +0000
Subject: [PATCH] BugID: 1084 - Improved polarization handling. - Fixed bug
where the next chunk was allocated before the previous was freed. - Updated
VisData interface.
---
.../BBSKernel/include/BBSKernel/Measurement.h | 2 +-
.../include/BBSKernel/MeasurementAIPS.h | 2 +-
.../BBSKernel/include/BBSKernel/Prediffer.h | 17 +-
CEP/BB/BBSKernel/include/BBSKernel/VisData.h | 16 +-
CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc | 5 +
CEP/BB/BBSKernel/src/MeasurementAIPS.cc | 25 +-
CEP/BB/BBSKernel/src/Model.cc | 55 ++-
CEP/BB/BBSKernel/src/Prediffer.cc | 348 ++++++++++--------
CEP/BB/BBSKernel/src/VisData.cc | 34 +-
9 files changed, 274 insertions(+), 230 deletions(-)
diff --git a/CEP/BB/BBSKernel/include/BBSKernel/Measurement.h b/CEP/BB/BBSKernel/include/BBSKernel/Measurement.h
index f1d7ef11a16..271875678df 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/Measurement.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/Measurement.h
@@ -67,7 +67,7 @@ public:
virtual VisGrid grid(const VisSelection &selection) const = 0;
- virtual void read(const VisSelection &selection, VisData::Pointer buffer,
+ virtual VisData::Pointer read(const VisSelection &selection,
const string &column = "DATA", bool readUVW = true) const = 0;
virtual void write(const VisSelection &selection, VisData::Pointer buffer,
diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MeasurementAIPS.h b/CEP/BB/BBSKernel/include/BBSKernel/MeasurementAIPS.h
index 9eaa11bd722..105b23eba61 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/MeasurementAIPS.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/MeasurementAIPS.h
@@ -57,7 +57,7 @@ public:
virtual VisGrid grid(const VisSelection &selection) const;
- virtual void read(const VisSelection &selection, VisData::Pointer buffer,
+ virtual VisData::Pointer read(const VisSelection &selection,
const string &column = "DATA", bool readUVW = true) const;
virtual void write(const VisSelection &selection, VisData::Pointer buffer,
diff --git a/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h b/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h
index 32f038fdc0f..ccdfe1d78fc 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h
@@ -88,8 +88,8 @@ struct ProcessingContext
{}
set<baseline_t> baselines;
- set<pair<size_t, size_t> > polarizations;
-
+ set<size_t> polarizations;
+ string outputColumn;
// Sum of the maximal number (over all solve domains) of partial derivatives
// of each parameter.
size_t derivativeCount;
@@ -196,8 +196,8 @@ private:
// </group>
bool setContext(const Context &context);
-
- void initializeSolveDomains(std::pair<size_t, size_t> size);
+ void initSolveDomains(std::pair<size_t, size_t> size);
+ void initPolarizationMap();
// Make all parameters non-solvable.
void clearSolvableParms();
@@ -216,7 +216,7 @@ private:
BaselineProcessor processor, void *arguments);
// Write the predicted data of a baseline.
- void predictBaseline(int threadnr, void* arguments, VisData::Pointer chunk,
+ void copyBaseline(int threadnr, void* arguments, VisData::Pointer chunk,
pair<size_t, size_t> offset, const MeqRequest& request,
baseline_t baseline, bool showd = false);
@@ -225,11 +225,6 @@ private:
pair<size_t, size_t> offset, const MeqRequest& request,
baseline_t baseline, bool showd = false);
- // Correct the data of a baseline.
- void correctBaseline(int threadnr, void* arguments, VisData::Pointer chunk,
- pair<size_t, size_t> offset, const MeqRequest& request,
- baseline_t baseline, bool showd = false);
-
// Generate equations for a baseline.
void generateBaseline(int threadnr, void* arguments, VisData::Pointer chunk,
pair<size_t, size_t> offset, const MeqRequest& request,
@@ -242,8 +237,8 @@ private:
size_t itsSubband;
string itsInputColumn;
- string itsOutputColumn;
bool itsReadUVW;
+ vector<int> itsPolarizationMap;
scoped_ptr<ParmDB::ParmDB> itsInstrumentDBase;
scoped_ptr<ParmDB::ParmDB> itsSkyDBase;
diff --git a/CEP/BB/BBSKernel/include/BBSKernel/VisData.h b/CEP/BB/BBSKernel/include/BBSKernel/VisData.h
index 571a98ea96e..b68e04c855a 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/VisData.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/VisData.h
@@ -70,11 +70,8 @@ public:
N_TimeslotFlag
};
- VisData()
- {}
- VisData(const VisGrid &grid);
-
- void reset(const VisGrid &in);
+ VisData(const VisGrid &visGrid);
+ ~VisData();
bool hasBaseline(baseline_t baseline) const;
size_t getBaselineIndex(baseline_t baseline) const;
@@ -85,15 +82,16 @@ public:
// Grid on which the visibility data is sampled.
VisGrid grid;
- // Indexes
- map<baseline_t, size_t> baselines;
- map<string, size_t> polarizations;
-
// Data
boost::multi_array<double, 3> uvw;
boost::multi_array<tslot_flag_t, 2> tslot_flag;
boost::multi_array<flag_t, 4> vis_flag;
boost::multi_array<sample_t, 4> vis_data;
+
+private:
+ // Indexes
+ map<baseline_t, size_t> baselines;
+ map<string, size_t> polarizations;
};
} //# namespace BBS
diff --git a/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc b/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc
index d0c88f4084a..6689507c7b7 100644
--- a/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc
+++ b/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc
@@ -105,6 +105,7 @@ MeqResult MeqParmFunklet::getResult (const MeqRequest& request)
}
return res;
}
+
// Get the domain, etc.
const MeqDomain& domain = request.domain();
MeqResult result(0);
@@ -125,6 +126,10 @@ MeqResult MeqParmFunklet::getResult (const MeqRequest& request)
result.setPerturbedValue (itsPertInx+ip, matp);
result.setPerturbedParm (itsPertInx+ip, this);
}
+
+// LOG_DEBUG_STR("Name: " << getName() << " nrpert: " << nrpert << " ["
+// << itsPertInx << " - " << itsPertInx + nrpert - 1 << "]");
+
// Iterate over all funklets.
// Evaluate one if its domain overlaps the request domain.
int sty = 0;
diff --git a/CEP/BB/BBSKernel/src/MeasurementAIPS.cc b/CEP/BB/BBSKernel/src/MeasurementAIPS.cc
index 34919eb626f..0a714989216 100644
--- a/CEP/BB/BBSKernel/src/MeasurementAIPS.cc
+++ b/CEP/BB/BBSKernel/src/MeasurementAIPS.cc
@@ -138,8 +138,8 @@ VisGrid MeasurementAIPS::grid(const VisSelection &selection) const
}
-void MeasurementAIPS::read(const VisSelection &selection,
- VisData::Pointer buffer, const string &column, bool readUVW) const
+VisData::Pointer MeasurementAIPS::read(const VisSelection &selection,
+ const string &column, bool readUVW) const
{
NSTimer readTimer, copyTimer;
@@ -148,7 +148,8 @@ void MeasurementAIPS::read(const VisSelection &selection,
Table tab_selection = itsMS(taqlExpr);
ASSERTSTR(tab_selection.nrow() > 0, "Data selection empty!");
- buffer->reset(grid(tab_selection, slicer));
+ VisGrid visGrid(grid(tab_selection, slicer));
+ VisData::Pointer buffer(new VisData(visGrid));
size_t nChannels = buffer->grid.freq.size();
size_t nPolarizations = buffer->grid.polarization.size();
@@ -219,15 +220,12 @@ void MeasurementAIPS::read(const VisSelection &selection,
for(uInt i = 0; i < nRows; ++i)
{
// Get time sequence for this baseline.
- map<baseline_t, size_t>::iterator it =
- buffer->baselines.find(baseline_t(aips_antenna1[i],
+ size_t baseline =
+ buffer->getBaselineIndex(baseline_t(aips_antenna1[i],
aips_antenna2[i]));
- ASSERTSTR(it != buffer->baselines.end(), "Unknown baseline!");
- size_t baseline = it->second;
-
// Flag timeslot as available.
- buffer->tslot_flag[baseline][tslot] = 0;
+ buffer->tslot_flag[baseline][tslot] &= ~VisData::UNAVAILABLE;
// Copy row (timeslot) flag.
if(aips_flag_row(i))
@@ -259,6 +257,8 @@ void MeasurementAIPS::read(const VisSelection &selection,
LOG_DEBUG_STR("Read time: " << readTimer);
LOG_DEBUG_STR("Copy time: " << copyTimer);
+
+ return buffer;
}
@@ -339,13 +339,10 @@ void MeasurementAIPS::write(const VisSelection &selection,
for(uInt i = 0; i < nRows; ++i)
{
// Get time sequence for this baseline.
- map<baseline_t, size_t>::iterator it =
- buffer->baselines.find(baseline_t(aips_antenna1[i],
+ size_t baseline =
+ buffer->getBaselineIndex(baseline_t(aips_antenna1[i],
aips_antenna2[i]));
- ASSERTSTR(it != buffer->baselines.end(), "Unknown baseline!");
- size_t baseline = it->second;
-
if(writeFlags)
{
// Write row flag.
diff --git a/CEP/BB/BBSKernel/src/Model.cc b/CEP/BB/BBSKernel/src/Model.cc
index bd5fa25e4dd..e1aea1240c5 100644
--- a/CEP/BB/BBSKernel/src/Model.cc
+++ b/CEP/BB/BBSKernel/src/Model.cc
@@ -88,15 +88,26 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
++it)
{
if(*it == "GAIN")
+ {
mask[GAIN] = true;
+ }
else if(*it == "DIRECTIONAL_GAIN")
+ {
mask[DIRECTIONAL_GAIN] = true;
+ }
else if(*it == "DIPOLE_BEAM")
+ {
mask[DIPOLE_BEAM] = true;
+ LOG_DEBUG_STR("Using dipole beam....");
+ }
else if(*it == "BANDPASS")
+ {
mask[BANDPASS] = true;
+ }
else if(*it == "PHASORS")
+ {
mask[PHASORS] = true;
+ }
}
string part1("real:");
@@ -151,7 +162,7 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
for(size_t i = 0; i < nSources; ++i)
{
azel[i] =
- new MeqAzEl(*(itsSourceNodes[i]), *(itsStationNodes[i].get()));
+ new MeqAzEl(*(itsSourceNodes[i]), *(itsStationNodes[0].get()));
dipole_beam[i] = new MeqDipoleBeam(azel[i], 1.706, 1.38,
casa::C::pi / 4.001, -casa::C::pi_4);
@@ -575,7 +586,7 @@ void Model::setStationUVW(const Instrument &instrument, VisData::Pointer buffer)
fill(statDone.begin(), statDone.end(), false);
// Set UVW of first station used to 0 (UVW coordinates are relative!).
- size_t station0 = buffer->baselines.begin()->first.first;
+ size_t station0 = buffer->grid.baseline[0].first;
statUVW[3 * station0] = 0.0;
statUVW[3 * station0 + 1] = 0.0;
statUVW[3 * station0 + 2] = 0.0;
@@ -586,38 +597,38 @@ void Model::setStationUVW(const Instrument &instrument, VisData::Pointer buffer)
do
{
nDone = 0;
- for (map<baseline_t, size_t>::const_iterator it =
- buffer->baselines.begin();
- it != buffer->baselines.end();
- ++it)
+ for(size_t idx = 0; idx < buffer->grid.baseline.size(); ++idx)
{
- size_t blindex = it->second;
- if(buffer->tslot_flag[blindex][tslot])
+ // If the contents of the UVW column is uninitialized, skip
+ // it.
+ if(buffer->tslot_flag[idx][tslot] & VisData::UNAVAILABLE)
+ {
continue;
-
- size_t statA = it->first.first;
- size_t statB = it->first.second;
- if (statDone[statA] && !statDone[statB])
+ }
+
+ size_t statA = buffer->grid.baseline[idx].first;
+ size_t statB = buffer->grid.baseline[idx].second;
+ if(statDone[statA] && !statDone[statB])
{
statUVW[3 * statB] =
- buffer->uvw[blindex][tslot][0] - statUVW[3 * statA];
+ buffer->uvw[idx][tslot][0] - statUVW[3 * statA];
statUVW[3 * statB + 1] =
- buffer->uvw[blindex][tslot][1] - statUVW[3 * statA + 1];
+ buffer->uvw[idx][tslot][1] - statUVW[3 * statA + 1];
statUVW[3 * statB + 2] =
- buffer->uvw[blindex][tslot][2] - statUVW[3 * statA + 2];
+ buffer->uvw[idx][tslot][2] - statUVW[3 * statA + 2];
statDone[statB] = true;
itsUVWNodes[statB]->set(time, statUVW[3 * statB],
statUVW[3 * statB + 1], statUVW[3 * statB + 2]);
++nDone;
}
- else if (statDone[statB] && !statDone[statA])
+ else if(statDone[statB] && !statDone[statA])
{
statUVW[3 * statA] =
- statUVW[3 * statB] - buffer->uvw[blindex][tslot][0];
+ statUVW[3 * statB] - buffer->uvw[idx][tslot][0];
statUVW[3 * statA + 1] =
- statUVW[3 * statB + 1] - buffer->uvw[blindex][tslot][1];
+ statUVW[3 * statB + 1] - buffer->uvw[idx][tslot][1];
statUVW[3 * statA + 2] =
- statUVW[3 * statB + 2] - buffer->uvw[blindex][tslot][2];
+ statUVW[3 * statB + 2] - buffer->uvw[idx][tslot][2];
statDone[statA] = true;
itsUVWNodes[statA]->set(time, statUVW[3 * statA],
statUVW[3 * statA + 1], statUVW[3 * statA + 2]);
@@ -626,6 +637,12 @@ void Model::setStationUVW(const Instrument &instrument, VisData::Pointer buffer)
}
}
while(nDone > 0);
+
+ for(size_t i = 0; i < statDone.size(); ++i)
+ {
+ ASSERTSTR(statDone[i], "UVW of station " << i << " could not be"
+ " determined for timeslot " << tslot);
+ }
}
}
diff --git a/CEP/BB/BBSKernel/src/Prediffer.cc b/CEP/BB/BBSKernel/src/Prediffer.cc
index 19b13f8e7e4..ba7e6a8b390 100644
--- a/CEP/BB/BBSKernel/src/Prediffer.cc
+++ b/CEP/BB/BBSKernel/src/Prediffer.cc
@@ -75,6 +75,7 @@ using LOFAR::operator<<;
using LOFAR::ParmDB::ParmDB;
using LOFAR::ParmDB::ParmDBMeta;
+
Prediffer::Prediffer(const string &measurement, size_t subband,
const string &inputColumn, const string &skyDBase,
const string &instrumentDBase, const string &historyDBase, bool readUVW)
@@ -109,6 +110,9 @@ Prediffer::Prediffer(const string &measurement, size_t subband,
{
LOG_INFO_STR("Input: " << measurement << "::" << itsInputColumn);
itsMeasurement.reset(new MeasurementAIPS(measurement, 0, subband, 0));
+
+ // Initialize the polarization map.
+ initPolarizationMap();
}
catch(casa::AipsError &_ex)
{
@@ -165,11 +169,9 @@ Prediffer::Prediffer(const string &measurement, size_t subband,
LOG_INFO_STR("UVW source: " << (readUVW ? "read" : "computed"));
LOG_INFO_STR("UVW convention: " << (readUVW ? "as recorded in the input"
" measurement" : "ANTENNA2 - ANTENNA1 (AIPS++ convention)"));
-// ASSERTSTR(computeUVW, "Reading of UVW coordinates from the input"
-// " measurement temporarily disabled.");
// Initialize chunk.
- itsChunkData.reset(new VisData());
+// itsChunkData.reset(new VisData());
// Initialize model.
casa::MEpoch startTimeMeas = itsMeasurement->getTimeRange().first;
@@ -203,6 +205,48 @@ Prediffer::~Prediffer()
}
+void Prediffer::initPolarizationMap()
+{
+ // Create a lookup table that maps external (measurement) polarization
+ // indices to internal polarization indices. The internal indices are:
+ // 0 - XX or LL
+ // 1 - XY or LR
+ // 2 - YX or RL
+ // 3 - YY or RR
+ // No assumptions can be made regarding the external polarization indices,
+ // which is why we need this table in the first place.
+
+ const vector<string> &polarizations = itsMeasurement->getPolarizations();
+
+ itsPolarizationMap.resize(polarizations.size());
+ for(size_t i = 0; i < polarizations.size(); ++i)
+ {
+ if(polarizations[i] == "XX" || polarizations[i] == "LL")
+ {
+ itsPolarizationMap[i] = 0;
+ }
+ else if(polarizations[i] == "XY" || polarizations[i] == "LR")
+ {
+ itsPolarizationMap[i] = 1;
+ }
+ else if(polarizations[i] == "YX" || polarizations[i] == "RL")
+ {
+ itsPolarizationMap[i] = 2;
+ }
+ else if(polarizations[i] == "YY" || polarizations[i] == "RR")
+ {
+ itsPolarizationMap[i] = 3;
+ }
+ else
+ {
+ LOG_WARN_STR("Don't know how to process polarization "
+ << polarizations[i] << "; will be skipped.");
+ itsPolarizationMap[i] = -1;
+ }
+ }
+}
+
+
void Prediffer::setSelection(VisSelection selection)
{
itsChunkSelection = selection;
@@ -251,10 +295,15 @@ bool Prediffer::nextChunk()
// Clear equations.
itsModel->clearEquations();
+ // Deallocate chunk.
+ ASSERTSTR(itsChunkData.use_count() == 0 || itsChunkData.use_count() == 1,
+ "itsChunkData shoud be unique (or uninitialized) by now.");
+ itsChunkData.reset();
+
// Read data.
LOG_DEBUG("Reading chunk...");
- itsMeasurement->read(itsChunkSelection, itsChunkData, itsInputColumn,
- itsReadUVW);
+ itsChunkData =
+ itsMeasurement->read(itsChunkSelection, itsInputColumn, itsReadUVW);
if(itsReadUVW)
itsModel->setStationUVW(itsMeasurement->getInstrument(), itsChunkData);
@@ -309,35 +358,36 @@ bool Prediffer::setContext(const Context &context)
itsContext = ProcessingContext();
// Select polarizations.
- vector<string> requested;
if(context.correlation.type.empty())
- requested = itsMeasurement->getPolarizations();
- else
- requested = context.correlation.type;
-
- for(vector<string>::const_iterator it = requested.begin();
- it != requested.end();
- ++it)
{
- size_t idx;
-
- try
+ for(size_t i = 0; i < itsPolarizationMap.size(); ++i)
{
- idx = itsChunkData->getPolarizationIndex(*it);
- if(*it == "XX" || *it == "LL")
- itsContext.polarizations.insert(make_pair(idx, 0));
- else if(*it == "XY" || *it == "LR")
- itsContext.polarizations.insert(make_pair(idx, 1));
- else if(*it == "YX" || *it == "RL")
- itsContext.polarizations.insert(make_pair(idx, 2));
- else if(*it == "YY" || *it == "RR")
- itsContext.polarizations.insert(make_pair(idx, 3));
- else
- LOG_WARN_STR("Don't know how to process polarization " << *it
- << "; skipping");
+ if(itsPolarizationMap[i] >= 0)
+ {
+ itsContext.polarizations.insert(i);
+ }
}
- catch(BBSKernelException &ex)
+ }
+ else
+ {
+ const vector<string> &requested = context.correlation.type;
+ const vector<string> &available = itsMeasurement->getPolarizations();
+
+ for(size_t i = 0; i < available.size(); ++i)
{
+ // Skip all unknown polarizations.
+ if(itsPolarizationMap[i] >= 0)
+ {
+ // Check if this polarization needs to be processed.
+ for(size_t j = 0; j < requested.size(); ++j)
+ {
+ if(available[i] == requested[j])
+ {
+ itsContext.polarizations.insert(i);
+ break;
+ }
+ }
+ }
}
}
@@ -353,12 +403,12 @@ bool Prediffer::setContext(const Context &context)
{
// If no baselines are speficied, use all baselines selected in the
// strategy that match the correlation selection of this context.
- map<baseline_t, size_t>::const_iterator it =
- itsChunkData->baselines.begin();
-
- while(it != itsChunkData->baselines.end())
+ for(vector<baseline_t>::const_iterator it =
+ itsChunkData->grid.baseline.begin();
+ it != itsChunkData->grid.baseline.end();
+ ++it)
{
- const baseline_t &baseline = it->first;
+ const baseline_t &baseline = *it;
if(context.correlation.selection == "ALL"
|| (baseline.first == baseline.second
@@ -366,7 +416,7 @@ bool Prediffer::setContext(const Context &context)
|| (baseline.first != baseline.second
&& context.correlation.selection == "CROSS"))
{
- itsContext.baselines.insert(it->first);
+ itsContext.baselines.insert(baseline);
}
++it;
}
@@ -423,8 +473,10 @@ bool Prediffer::setContext(const Context &context)
{
baseline_t baseline(*it1, *it2);
- if (itsChunkData->hasBaseline(baseline))
+ if(itsChunkData->hasBaseline(baseline))
+ {
itsContext.baselines.insert(baseline);
+ }
}
}
}
@@ -464,7 +516,7 @@ bool Prediffer::setContext(const PredictContext &context)
it->second.fillFunklets(itsParmValues, itsWorkDomain);
}
- itsOutputColumn = context.outputColumn;
+ itsContext.outputColumn = context.outputColumn;
return true;
}
@@ -489,7 +541,7 @@ bool Prediffer::setContext(const SubtractContext &context)
it->second.fillFunklets(itsParmValues, itsWorkDomain);
}
- itsOutputColumn = context.outputColumn;
+ itsContext.outputColumn = context.outputColumn;
return true;
}
@@ -514,8 +566,7 @@ bool Prediffer::setContext(const CorrectContext &context)
it->second.fillFunklets(itsParmValues, itsWorkDomain);
}
-
- itsOutputColumn = context.outputColumn;
+ itsContext.outputColumn = context.outputColumn;
return true;
}
@@ -618,7 +669,7 @@ bool Prediffer::setContext(const GenerateContext &context)
LOG_DEBUG_STR("Nominal solve domain size: " << domainSize.first
<< " channels by " << domainSize.second << " timeslots.");
- initializeSolveDomains(domainSize);
+ initSolveDomains(domainSize);
return (itsContext.derivativeCount > 0);
}
@@ -629,12 +680,14 @@ void Prediffer::predict()
process(false, true, false, make_pair(0, 0),
make_pair(itsChunkData->grid.freq.size() - 1,
itsChunkData->grid.time.size() - 1),
- &Prediffer::predictBaseline, 0);
-
- if(!itsOutputColumn.empty())
- itsMeasurement->write(itsChunkSelection, itsChunkData, itsOutputColumn,
- false);
+ &Prediffer::copyBaseline, 0);
+ if(!itsContext.outputColumn.empty())
+ {
+ itsMeasurement->write(itsChunkSelection, itsChunkData,
+ itsContext.outputColumn, false);
+ }
+
LOG_DEBUG_STR("Predict: " << itsPredTimer);
LOG_DEBUG_STR("Copy: " << itsEqTimer);
itsPredTimer.reset();
@@ -649,9 +702,11 @@ void Prediffer::subtract()
itsChunkData->grid.time.size() - 1),
&Prediffer::subtractBaseline, 0);
- if(!itsOutputColumn.empty())
- itsMeasurement->write(itsChunkSelection, itsChunkData, itsOutputColumn,
- false);
+ if(!itsContext.outputColumn.empty())
+ {
+ itsMeasurement->write(itsChunkSelection, itsChunkData,
+ itsContext.outputColumn, false);
+ }
LOG_DEBUG_STR("Predict: " << itsPredTimer);
LOG_DEBUG_STR("Subtract: " << itsEqTimer);
@@ -665,11 +720,13 @@ void Prediffer::correct()
process(false, true, false, make_pair(0, 0),
make_pair(itsChunkData->grid.freq.size() - 1,
itsChunkData->grid.time.size() - 1),
- &Prediffer::correctBaseline, 0);
+ &Prediffer::copyBaseline, 0);
- if(!itsOutputColumn.empty())
- itsMeasurement->write(itsChunkSelection, itsChunkData, itsOutputColumn,
- false);
+ if(!itsContext.outputColumn.empty())
+ {
+ itsMeasurement->write(itsChunkSelection, itsChunkData,
+ itsContext.outputColumn, false);
+ }
LOG_DEBUG_STR("Correct: " << itsEqTimer);
LOG_DEBUG_STR("Copy: " << itsPredTimer);
@@ -812,7 +869,7 @@ void Prediffer::readParms()
}
-void Prediffer::initializeSolveDomains(pair<size_t, size_t> size)
+void Prediffer::initSolveDomains(pair<size_t, size_t> size)
{
// Need to implement some way to select local solve domains
// from global solve domains, see initSolvableParms().
@@ -901,6 +958,7 @@ void Prediffer::initializeSolveDomains(pair<size_t, size_t> size)
}
itsContext.derivativeCount = nDerivatives;
LOG_DEBUG_STR("nDerivatives: " << nDerivatives);
+ LOG_DEBUG_STR("domainUnknownCount: " << domainUnknownCount);
}
@@ -928,7 +986,7 @@ void Prediffer::process(bool useFlags, bool precalc, bool derivatives,
int nChannels = end.first - start.first + 1;
int nTimeslots = end.second - start.second + 1;
- // NOTE: Temporary vector should be removed after MNS overhaul.
+ // NOTE: Temporary vector; should be removed after MNS overhaul.
vector<double> timeAxis(nTimeslots + 1);
for(size_t i = 0; i < nTimeslots; ++i)
timeAxis[i] = itsChunkData->grid.time.lower(start.second + i);
@@ -951,10 +1009,11 @@ void Prediffer::process(bool useFlags, bool precalc, bool derivatives,
request.setOffset(start);
// Precalculate if required.
- // TODO: Fix this when doing a correct!
if(precalc)
+ {
itsModel->precalculate(request);
-
+ }
+
/************** DEBUG DEBUG DEBUG **************/
/*
cout << "Processing correlations: ";
@@ -975,28 +1034,25 @@ void Prediffer::process(bool useFlags, bool precalc, bool derivatives,
#pragma omp parallel
{
#pragma omp for schedule(dynamic)
- // Process all selected baselines.
- for(set<baseline_t>::const_iterator it = itsContext.baselines.begin();
- it != itsContext.baselines.end();
- ++it)
- {
-// LOG_DEBUG_STR("Processing baseline " << it->first << " - "
-// << it->second);
-
+ // Process all selected baselines.
+ for(set<baseline_t>::const_iterator it = itsContext.baselines.begin();
+ it != itsContext.baselines.end();
+ ++it)
+ {
#if defined _OPENMP
- size_t thread = omp_get_thread_num();
+ size_t thread = omp_get_thread_num();
#else
- size_t thread = 0;
+ size_t thread = 0;
#endif
- // Process baseline.
- (this->*processor)
- (thread, arguments, itsChunkData, start, request, *it, false);
- }
+ // Process baseline.
+ (this->*processor)
+ (thread, arguments, itsChunkData, start, request, *it, false);
+ }
} // omp parallel
}
-void Prediffer::predictBaseline(int, void*, VisData::Pointer chunk,
+void Prediffer::copyBaseline(int, void*, VisData::Pointer chunk,
pair<size_t, size_t> offset, const MeqRequest& request, baseline_t baseline,
bool showd)
{
@@ -1021,62 +1077,9 @@ void Prediffer::predictBaseline(int, void*, VisData::Pointer chunk,
size_t nChannels = request.nx();
size_t nTimeslots = request.ny();
- // Get baseline index.
- size_t basel = chunk->getBaselineIndex(baseline);
-
- // Copy the data to the right location.
- for(size_t tslot = offset.second;
- tslot < offset.second + nTimeslots;
- ++tslot)
- {
- for(set<pair<size_t, size_t> >::const_iterator it =
- itsContext.polarizations.begin();
- it != itsContext.polarizations.end();
- ++it)
- {
- size_t pol = it->first;
-
- const double* re = resultRe[it->second];
- const double* im = resultIm[it->second];
-
- for(size_t chan = 0; chan < nChannels; ++chan)
- {
- chunk->vis_data[basel][tslot][offset.first + chan][pol] =
- makefcomplex(re[chan], im[chan]);
- }
-
- resultRe[pol] += nChannels;
- resultIm[pol] += nChannels;
- }
- }
- itsEqTimer.stop();
-}
-
-
-void Prediffer::subtractBaseline(int, void*, VisData::Pointer chunk,
- pair<size_t, size_t> offset, const MeqRequest& request, baseline_t baseline,
- bool showd)
-{
- // Do the actual correct.
- itsPredTimer.start();
- MeqJonesResult jresult = itsModel->evaluate(baseline, request);
- itsPredTimer.stop();
-
- itsEqTimer.start();
- // Put the results in a single array for easier handling.
- const double *resultRe[4], *resultIm[4];
- jresult.getResult11().getValue().dcomplexStorage(resultRe[0],
- resultIm[0]);
- jresult.getResult12().getValue().dcomplexStorage(resultRe[1],
- resultIm[1]);
- jresult.getResult21().getValue().dcomplexStorage(resultRe[2],
- resultIm[2]);
- jresult.getResult22().getValue().dcomplexStorage(resultRe[3],
- resultIm[3]);
-
- // Get information about request grid.
- size_t nChannels = request.nx();
- size_t nTimeslots = request.ny();
+ // Check preconditions.
+ ASSERT(offset.first + nChannels <= chunk->grid.freq.size());
+ ASSERT(offset.second + nTimeslots <= chunk->grid.time.size());
// Get baseline index.
size_t basel = chunk->getBaselineIndex(baseline);
@@ -1086,31 +1089,32 @@ void Prediffer::subtractBaseline(int, void*, VisData::Pointer chunk,
tslot < offset.second + nTimeslots;
++tslot)
{
- for(set<pair<size_t, size_t> >::const_iterator it =
- itsContext.polarizations.begin();
+ for(set<size_t>::const_iterator it = itsContext.polarizations.begin();
it != itsContext.polarizations.end();
++it)
{
- size_t pol = it->first;
-
- const double* re = resultRe[it->second];
- const double* im = resultIm[it->second];
+ // External (Measurement) polarization index.
+ size_t ext_pol = *it;
+ // Internal (MEQ) polarization index.
+ int int_pol = itsPolarizationMap[ext_pol];
+ const double *re = resultRe[int_pol];
+ const double *im = resultIm[int_pol];
for(size_t chan = 0; chan < nChannels; ++chan)
{
- chunk->vis_data[basel][tslot][offset.first + chan][pol] -=
+ chunk->vis_data[basel][tslot][offset.first + chan][ext_pol] =
makefcomplex(re[chan], im[chan]);
}
- resultRe[pol] += nChannels;
- resultIm[pol] += nChannels;
+ resultRe[int_pol] += nChannels;
+ resultIm[int_pol] += nChannels;
}
}
itsEqTimer.stop();
}
-void Prediffer::correctBaseline(int, void*, VisData::Pointer chunk,
+void Prediffer::subtractBaseline(int, void*, VisData::Pointer chunk,
pair<size_t, size_t> offset, const MeqRequest& request, baseline_t baseline,
bool showd)
{
@@ -1134,40 +1138,45 @@ void Prediffer::correctBaseline(int, void*, VisData::Pointer chunk,
// Get information about request grid.
size_t nChannels = request.nx();
size_t nTimeslots = request.ny();
+
+ // Check preconditions.
+ ASSERT(offset.first + nChannels <= chunk->grid.freq.size());
+ ASSERT(offset.second + nTimeslots <= chunk->grid.time.size());
// Get baseline index.
size_t basel = chunk->getBaselineIndex(baseline);
-
+
// Copy the data to the right location.
for(size_t tslot = offset.second;
tslot < offset.second + nTimeslots;
++tslot)
{
- for(set<pair<size_t, size_t> >::const_iterator it =
- itsContext.polarizations.begin();
+ for(set<size_t>::const_iterator it = itsContext.polarizations.begin();
it != itsContext.polarizations.end();
++it)
{
- size_t pol = it->first;
-
- const double* re = resultRe[it->second];
- const double* im = resultIm[it->second];
+ // External (Measurement) polarization index.
+ size_t ext_pol = *it;
+ // Internal (MEQ) polarization index.
+ int int_pol = itsPolarizationMap[ext_pol];
+ const double *re = resultRe[int_pol];
+ const double *im = resultIm[int_pol];
for(size_t chan = 0; chan < nChannels; ++chan)
{
- chunk->vis_data[basel][tslot][offset.first + chan][pol] =
+ chunk->vis_data[basel][tslot][offset.first + chan][ext_pol] -=
makefcomplex(re[chan], im[chan]);
}
- resultRe[pol] += nChannels;
- resultIm[pol] += nChannels;
+ resultRe[int_pol] += nChannels;
+ resultIm[int_pol] += nChannels;
}
}
itsEqTimer.stop();
}
-void Prediffer::generateBaseline(int threadnr, void* arguments,
+void Prediffer::generateBaseline(int threadnr, void *arguments,
VisData::Pointer chunk, pair<size_t, size_t> offset,
const MeqRequest& request, baseline_t baseline, bool showd)
{
@@ -1205,20 +1214,26 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
size_t nChannels = request.nx();
size_t nTimeslots = request.ny();
+ // Check preconditions.
+ ASSERT(offset.first + nChannels <= chunk->grid.freq.size());
+ ASSERT(offset.second + nTimeslots <= chunk->grid.time.size());
+
// Get baseline index.
size_t basel = chunk->getBaselineIndex(baseline);
// To avoid having to use large temporary arrays, step through the
- // data by timestep and correlation.
- for(set<pair<size_t, size_t> >::const_iterator it =
- itsContext.polarizations.begin();
+ // data by timeslot and correlation.
+ for(set<size_t>::const_iterator it = itsContext.polarizations.begin();
it != itsContext.polarizations.end();
++it)
{
- size_t pol = it->first;
+ // External (Measurement) polarization index.
+ size_t ext_pol = *it;
+ // Internal (MEQ) polarization index.
+ int int_pol = itsPolarizationMap[ext_pol];
// Get the results for this correlation.
- const MeqResult &tcres = *predictResults[it->second];
+ const MeqResult &tcres = *predictResults[int_pol];
// Get pointers to the main data.
const MeqMatrix &val = tcres.getValue();
@@ -1236,6 +1251,7 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
{
if (tcres.isDefined(scinx))
{
+// LOG_DEBUG_STR("scinx: " << scinx << " #params: " << nrParamsFound);
indices[nrParamsFound] = scinx;
const MeqMatrix& valp = tcres.getPerturbedValue(scinx);
valp.dcomplexStorage (pertReal[nrParamsFound],
@@ -1259,6 +1275,20 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
// }
/************** DEBUG DEBUG DEBUG **************/
+ // Skip timeslot if flagged.
+ if(chunk->tslot_flag[basel][tslot])
+ {
+ // Move pointers to the next timeslot.
+ for(int scinx = 0; scinx < nrParamsFound; ++scinx)
+ {
+ pertReal[scinx] += nChannels;
+ pertImag[scinx] += nChannels;
+ }
+ realVals += nChannels;
+ imagVals += nChannels;
+ continue;
+ }
+
// Compute relative solver index (time part).
size_t solverIndexTime =
((tslot - offset.second) / itsContext.domainSize.second)
@@ -1270,7 +1300,7 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
ch < offset.first + nChannels;
++ch)
{
- if (!chunk->vis_flag[basel][tslot][ch][pol])
+ if (!chunk->vis_flag[basel][tslot][ch][ext_pol])
{
// Compute relative solver index.
size_t solverIndex = solverIndexTime +
@@ -1278,10 +1308,10 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
// Compute right hand side of the equation pair.
double diffr =
- real(chunk->vis_data[basel][tslot][ch][pol])
+ real(chunk->vis_data[basel][tslot][ch][ext_pol])
- *realVals;
double diffi =
- imag(chunk->vis_data[basel][tslot][ch][pol])
+ imag(chunk->vis_data[basel][tslot][ch][ext_pol])
- *imagVals;
#ifdef COMPUTE_SQUARED_ERROR
@@ -1308,7 +1338,7 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
// Now add the equations to the correct fitter
// object.
- if (nrParamsFound != itsContext.derivativeCount)
+ if(nrParamsFound != itsContext.derivativeCount)
{
solvers[solverIndex]->makeNorm(nrParamsFound,
indices, resultr, 1.0, diffr);
@@ -1341,7 +1371,7 @@ void Prediffer::generateBaseline(int threadnr, void* arguments,
}
*/
/************** DEBUG DEBUG DEBUG **************/
- } // !chunk->vis_flag[basel][tslot][ch][pol]
+ } // !chunk->vis_flag[basel][tslot][ch][ext_pol]
// Move pointers to the next channel.
for (int scinx = 0; scinx < nrParamsFound; ++scinx)
diff --git a/CEP/BB/BBSKernel/src/VisData.cc b/CEP/BB/BBSKernel/src/VisData.cc
index 1323dafc32f..8ebf8583cdd 100644
--- a/CEP/BB/BBSKernel/src/VisData.cc
+++ b/CEP/BB/BBSKernel/src/VisData.cc
@@ -32,16 +32,9 @@ namespace LOFAR
namespace BBS
{
-VisData::VisData(const VisGrid &grid)
+VisData::VisData(const VisGrid &visGrid)
+ : grid(visGrid)
{
- reset(grid);
-}
-
-
-void VisData::reset(const VisGrid &grid)
-{
- this->grid = grid;
-
size_t nTimeslots = grid.time.size();
size_t nBaselines = grid.baseline.size();
size_t nChannels = grid.freq.size();
@@ -95,6 +88,12 @@ void VisData::reset(const VisGrid &grid)
}
+VisData::~VisData()
+{
+ LOG_DEBUG("VisData destructor called.");
+}
+
+
bool VisData::hasBaseline(baseline_t baseline) const
{
map<baseline_t, size_t>::const_iterator it = baselines.find(baseline);
@@ -106,10 +105,12 @@ size_t VisData::getBaselineIndex(baseline_t baseline) const
{
map<baseline_t, size_t>::const_iterator it = baselines.find(baseline);
if(it != baselines.end())
+ {
return it->second;
- else
- THROW(BBSKernelException, "Request for index of unknown baseline "
- << baseline.first << " - " << baseline.second);
+ }
+
+ THROW(BBSKernelException, "Request for index of unknown baseline "
+ << baseline.first << " - " << baseline.second);
}
@@ -123,12 +124,13 @@ bool VisData::hasPolarization(const string &polarization) const
size_t VisData::getPolarizationIndex(const string &polarization) const
{
map<string, size_t>::const_iterator it = polarizations.find(polarization);
-
if(it != polarizations.end())
+ {
return it->second;
- else
- THROW(BBSKernelException, "Request for index of unknown polarization "
- << polarization);
+ }
+
+ THROW(BBSKernelException, "Request for index of unknown polarization "
+ << polarization);
}
} //# namespace BBS
--
GitLab